← Back to branch summary

~ubuntu-branches/ubuntu/wily/openms/wily

« back to all changes in this revision

Viewing changes to source/MATH/STATISTICS/GaussFitter.C

  • Committer: Package Import Robot
  • Author(s): Filippo Rusconi
  • Date: 2013-12-20 11:30:16 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: package-import@ubuntu.com-20131220113016-wre5g9bteeheq6he
Tags: 1.11.1-3
* remove version number from libbost development package names;
* ensure that AUTHORS is correctly shipped in all packages.

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/MATH/STATISTICS/GaussFitter.C
1
 
// -*- mode: C++; tab-width: 2; -*-
2
 
// vi: set ts=2:
3
 
//
4
 
// --------------------------------------------------------------------------
5
 
//                   OpenMS Mass Spectrometry Framework
6
 
// --------------------------------------------------------------------------
7
 
//  Copyright (C) 2003-2011 -- Oliver Kohlbacher, Knut Reinert
8
 
//
9
 
//  This library is free software; you can redistribute it and/or
10
 
//  modify it under the terms of the GNU Lesser General Public
11
 
//  License as published by the Free Software Foundation; either
12
 
//  version 2.1 of the License, or (at your option) any later version.
13
 
//
14
 
//  This library is distributed in the hope that it will be useful,
15
 
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
16
 
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17
 
//  Lesser General Public License for more details.
18
 
//
19
 
//  You should have received a copy of the GNU Lesser General Public
20
 
//  License along with this library; if not, write to the Free Software
21
 
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
1
// --------------------------------------------------------------------------
 
2
//                   OpenMS -- Open-Source Mass Spectrometry
 
3
// --------------------------------------------------------------------------
 
4
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
 
5
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
 
6
//
 
7
// This software is released under a three-clause BSD license:
 
8
//  * Redistributions of source code must retain the above copyright
 
9
//    notice, this list of conditions and the following disclaimer.
 
10
//  * Redistributions in binary form must reproduce the above copyright
 
11
//    notice, this list of conditions and the following disclaimer in the
 
12
//    documentation and/or other materials provided with the distribution.
 
13
//  * Neither the name of any author or any participating institution
 
14
//    may be used to endorse or promote products derived from this software
 
15
//    without specific prior written permission.
 
16
// For a full list of authors, refer to the file AUTHORS.
 
17
// --------------------------------------------------------------------------
 
18
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 
19
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 
20
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 
21
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
 
22
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 
23
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 
24
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
 
25
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 
26
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
 
27
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 
28
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
22
29
//
23
30
// --------------------------------------------------------------------------
24
31
// $Maintainer: Andreas Bertsch $
38
45
 
39
46
namespace OpenMS
40
47
{
41
 
        namespace Math
42
 
        {
43
 
                GaussFitter::GaussFitter()
44
 
                {
45
 
                        init_param_.A = 0.06;
46
 
                        init_param_.x0 = 3.0;
47
 
                        init_param_.sigma = 0.5;
48
 
                }
49
 
        
50
 
                GaussFitter::~GaussFitter()
51
 
                {
52
 
                }
53
 
        
54
 
                void GaussFitter::setInitialParameters(const GaussFitResult& param)
55
 
                {
56
 
                        init_param_.A = param.A;
57
 
                        init_param_.x0 = param.x0;
58
 
                        init_param_.sigma = param.sigma;
59
 
                }
60
 
        
61
 
                const String& GaussFitter::getGnuplotFormula() const
62
 
                {
63
 
                        return gnuplot_formula_;
64
 
                }
65
 
        
66
 
                int GaussFitter::gaussFitterf_(const gsl_vector* x, void* params, gsl_vector* f)
67
 
                {
68
 
                        vector<DPosition<2> >* data = static_cast<vector<DPosition<2> >*>(params);
69
 
                        
70
 
                        double A = gsl_vector_get (x, 0);
71
 
                        double x0 = gsl_vector_get (x, 1);
72
 
                        double sig = gsl_vector_get (x, 2);
73
 
                
74
 
                        UInt i = 0;
75
 
                        for (vector<DPosition<2> >::iterator it = data->begin(); it != data->end(); ++it)
76
 
                        {
77
 
                                gsl_vector_set(f, i++,  A * exp(-1.0 * pow(it->getX() - x0, 2) / (2 * pow(sig, 2))) - it->getY());
78
 
                        }
79
 
                
80
 
                        return GSL_SUCCESS;
81
 
                }
82
 
        
83
 
                int GaussFitter::gaussFitterdf_(const gsl_vector* x, void* params, gsl_matrix* J)
84
 
                {
85
 
                        vector<DPosition<2> >* data = static_cast<vector<DPosition<2> >*>(params);
86
 
        
87
 
                        double A = gsl_vector_get (x, 0);
88
 
                        double x0 = gsl_vector_get (x, 1);
89
 
                        double sig = gsl_vector_get (x, 2);
90
 
                
91
 
                        UInt i = 0;
92
 
                        for (vector<DPosition<2> >::iterator it = data->begin(); it != data->end(); ++it)
93
 
                        {
94
 
                                gsl_matrix_set (J, i, 0, exp(-1.0 * pow(it->getX() - x0, 2) / (2 * pow(sig, 2))));
95
 
                                gsl_matrix_set (J, i, 1, (A * exp(-1.0 * pow(it->getX() - x0, 2) / (2 * pow(sig, 2))) * (-1 * (-2 * it->getX() + 2.0 * x0) / (2 * pow(sig, 2)))));
96
 
                                gsl_matrix_set (J, i++, 2, (A * exp(-1.0 * pow(it->getX() - x0, 2) / (2 * pow(sig, 2))) * (pow(it->getX() - x0, 2) / (4 * pow(sig, 3)))));                              
97
 
                        }
98
 
                  return GSL_SUCCESS;
99
 
                }
100
 
        
101
 
                int GaussFitter::gaussFitterfdf_(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* J)
102
 
                {
103
 
                  gaussFitterf_(x, params, f);
104
 
                  gaussFitterdf_(x, params, J);
105
 
                  return GSL_SUCCESS;
106
 
                }
107
 
        
108
 
#ifdef GAUSS_FITTER_VERBOSE
109
 
                void GaussFitter::printState_(size_t iter, gsl_multifit_fdfsolver * s)
110
 
                {
111
 
                  printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
112
 
                          "|f(x)| = %g\n",
113
 
                          iter,
114
 
                          gsl_vector_get(s->x, 0), 
115
 
                          gsl_vector_get(s->x, 1),
116
 
                          gsl_vector_get(s->x, 2), 
117
 
                          /*gsl_blas_dnrm2(s->f)*/0);
118
 
                }
119
 
#endif
120
 
        
121
 
                GaussFitter::GaussFitResult GaussFitter::fit(vector<DPosition<2> >& input)
122
 
                {
123
 
                  const gsl_multifit_fdfsolver_type *T = NULL;
124
 
                  gsl_multifit_fdfsolver *s = NULL;
125
 
                
126
 
                  int status(0);
127
 
                  size_t iter = 0;
128
 
                
129
 
                  const size_t p = 3;
130
 
                
131
 
                  //gsl_matrix *covar = gsl_matrix_alloc (p, p);
132
 
                  gsl_multifit_function_fdf f;
133
 
                  double x_init[3] = { init_param_.A, init_param_.x0, init_param_.sigma };
134
 
                  gsl_vector_view x = gsl_vector_view_array (x_init, p);
135
 
                  const gsl_rng_type * type = NULL;
136
 
                  gsl_rng * r = NULL;
137
 
                
138
 
                  gsl_rng_env_setup();
139
 
                
140
 
                  type = gsl_rng_default;
141
 
                  r = gsl_rng_alloc (type);
142
 
                
143
 
                  f.f = &gaussFitterf_;
144
 
                  f.df = &gaussFitterdf_;
145
 
                  f.fdf = &gaussFitterfdf_;
146
 
                  f.n = input.size();
147
 
                  f.p = p;
148
 
                  f.params = &input;
149
 
                
150
 
                  T = gsl_multifit_fdfsolver_lmsder;
151
 
                  s = gsl_multifit_fdfsolver_alloc (T, input.size(), p);
152
 
                  gsl_multifit_fdfsolver_set (s, &f, &x.vector);
153
 
        
154
 
                        #ifdef GAUSS_FITTER_VERBOSE
155
 
                  printState_(iter, s);
156
 
                        #endif
157
 
                
158
 
                  do
159
 
                  {
160
 
                    iter++;
161
 
                    status = gsl_multifit_fdfsolver_iterate (s);
162
 
                
163
 
                                #ifdef GAUSS_FITTER_VERBOSE
164
 
                    printf ("status = %s\n", gsl_strerror (status));
165
 
                    printState_(iter, s);
166
 
        
167
 
                                cerr << "f(x)=" << gsl_vector_get(s->x, 0) << " * exp(-(x - " << gsl_vector_get(s->x, 1) << ") ** 2 / 2 / (" << gsl_vector_get(s->x, 2) << ") ** 2)";
168
 
                                #endif
169
 
        
170
 
                    if (status)
171
 
                                {
172
 
                            break;
173
 
                                }
174
 
        
175
 
                    status = gsl_multifit_test_delta (s->dx, s->x,
176
 
                                                        1e-4, 1e-4);
177
 
                  }
178
 
                  while (status == GSL_CONTINUE && iter < 500);
179
 
        
180
 
                        if (status!=GSL_SUCCESS)
181
 
                        {
182
 
                                gsl_rng_free(r);
183
 
                                gsl_multifit_fdfsolver_free(s);
184
 
        
185
 
                                throw Exception::UnableToFit(__FILE__,__LINE__,__PRETTY_FUNCTION__,"UnableToFit-GaussFitter","Could not fit the gaussian to the data");
186
 
                        }
187
 
                  
188
 
                        // write the result in a GaussFitResult struct
189
 
                        GaussFitResult result;
190
 
                        result.A = gsl_vector_get(s->x, 0);
191
 
                        result.x0 = gsl_vector_get(s->x, 1);
192
 
                        result.sigma = gsl_vector_get(s->x, 2);
193
 
        
194
 
                        // build a formula with the fitted parameters for gnuplot
195
 
                        stringstream formula;
196
 
                        formula << "f(x)=" << result.A << " * exp(-(x - " << result.x0 << ") ** 2 / 2 / (" << result.sigma << ") ** 2)";
197
 
                        gnuplot_formula_ = formula.str();
198
 
                        #ifdef GAUSS_FITTER_VERBOSE
199
 
                        cout << gnuplot_formula_ << endl;
200
 
                        #endif
201
 
                        
202
 
                        gsl_rng_free(r);
203
 
                        gsl_multifit_fdfsolver_free (s);
204
 
        
205
 
                        return result;
206
 
                }
207
 
                
208
 
        } //namespace Math
 
48
  namespace Math
 
49
  {
 
50
    GaussFitter::GaussFitter()
 
51
    {
 
52
      init_param_.A = 0.06;
 
53
      init_param_.x0 = 3.0;
 
54
      init_param_.sigma = 0.5;
 
55
    }
 
56
 
 
57
    GaussFitter::~GaussFitter()
 
58
    {
 
59
    }
 
60
 
 
61
    void GaussFitter::setInitialParameters(const GaussFitResult & param)
 
62
    {
 
63
      init_param_.A = param.A;
 
64
      init_param_.x0 = param.x0;
 
65
      init_param_.sigma = param.sigma;
 
66
    }
 
67
 
 
68
    const String & GaussFitter::getGnuplotFormula() const
 
69
    {
 
70
      return gnuplot_formula_;
 
71
    }
 
72
 
 
73
    int GaussFitter::gaussFitterf_(const gsl_vector * x, void * params, gsl_vector * f)
 
74
    {
 
75
      vector<DPosition<2> > * data = static_cast<vector<DPosition<2> > *>(params);
 
76
 
 
77
      double A = gsl_vector_get(x, 0);
 
78
      double x0 = gsl_vector_get(x, 1);
 
79
      double sig = gsl_vector_get(x, 2);
 
80
 
 
81
      UInt i = 0;
 
82
      for (vector<DPosition<2> >::iterator it = data->begin(); it != data->end(); ++it)
 
83
      {
 
84
        gsl_vector_set(f, i++, A * exp(-1.0 * pow(it->getX() - x0, 2) / (2 * pow(sig, 2))) - it->getY());
 
85
      }
 
86
 
 
87
      return GSL_SUCCESS;
 
88
    }
 
89
 
 
90
    int GaussFitter::gaussFitterdf_(const gsl_vector * x, void * params, gsl_matrix * J)
 
91
    {
 
92
      vector<DPosition<2> > * data = static_cast<vector<DPosition<2> > *>(params);
 
93
 
 
94
      double A = gsl_vector_get(x, 0);
 
95
      double x0 = gsl_vector_get(x, 1);
 
96
      double sig = gsl_vector_get(x, 2);
 
97
 
 
98
      UInt i = 0;
 
99
      for (vector<DPosition<2> >::iterator it = data->begin(); it != data->end(); ++it)
 
100
      {
 
101
        gsl_matrix_set(J, i, 0, exp(-1.0 * pow(it->getX() - x0, 2) / (2 * pow(sig, 2))));
 
102
        gsl_matrix_set(J, i, 1, (A * exp(-1.0 * pow(it->getX() - x0, 2) / (2 * pow(sig, 2))) * (-1 * (-2 * it->getX() + 2.0 * x0) / (2 * pow(sig, 2)))));
 
103
        gsl_matrix_set(J, i++, 2, (A * exp(-1.0 * pow(it->getX() - x0, 2) / (2 * pow(sig, 2))) * (pow(it->getX() - x0, 2) / (4 * pow(sig, 3)))));
 
104
      }
 
105
      return GSL_SUCCESS;
 
106
    }
 
107
 
 
108
    int GaussFitter::gaussFitterfdf_(const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix * J)
 
109
    {
 
110
      gaussFitterf_(x, params, f);
 
111
      gaussFitterdf_(x, params, J);
 
112
      return GSL_SUCCESS;
 
113
    }
 
114
 
 
115
#ifdef GAUSS_FITTER_VERBOSE
 
116
    void GaussFitter::printState_(size_t iter, gsl_multifit_fdfsolver * s)
 
117
    {
 
118
      printf("iter: %3u x = % 15.8f % 15.8f % 15.8f "
 
119
             "|f(x)| = %g\n",
 
120
             iter,
 
121
             gsl_vector_get(s->x, 0),
 
122
             gsl_vector_get(s->x, 1),
 
123
             gsl_vector_get(s->x, 2),
 
124
             /*gsl_blas_dnrm2(s->f)*/ 0);
 
125
    }
 
126
 
 
127
#endif
 
128
 
 
129
    GaussFitter::GaussFitResult GaussFitter::fit(vector<DPosition<2> > & input)
 
130
    {
 
131
      const gsl_multifit_fdfsolver_type * T = NULL;
 
132
      gsl_multifit_fdfsolver * s = NULL;
 
133
 
 
134
      int status(0);
 
135
      size_t iter = 0;
 
136
 
 
137
      const size_t p = 3;
 
138
 
 
139
      //gsl_matrix *covar = gsl_matrix_alloc (p, p);
 
140
      gsl_multifit_function_fdf f;
 
141
      double x_init[3] = { init_param_.A, init_param_.x0, init_param_.sigma };
 
142
      gsl_vector_view x = gsl_vector_view_array(x_init, p);
 
143
      const gsl_rng_type * type = NULL;
 
144
      gsl_rng * r = NULL;
 
145
 
 
146
      gsl_rng_env_setup();
 
147
 
 
148
      type = gsl_rng_default;
 
149
      r = gsl_rng_alloc(type);
 
150
 
 
151
      f.f = &gaussFitterf_;
 
152
      f.df = &gaussFitterdf_;
 
153
      f.fdf = &gaussFitterfdf_;
 
154
      f.n = input.size();
 
155
      f.p = p;
 
156
      f.params = &input;
 
157
 
 
158
      T = gsl_multifit_fdfsolver_lmsder;
 
159
      s = gsl_multifit_fdfsolver_alloc(T, input.size(), p);
 
160
      gsl_multifit_fdfsolver_set(s, &f, &x.vector);
 
161
 
 
162
#ifdef GAUSS_FITTER_VERBOSE
 
163
      printState_(iter, s);
 
164
#endif
 
165
 
 
166
      do
 
167
      {
 
168
        iter++;
 
169
        status = gsl_multifit_fdfsolver_iterate(s);
 
170
 
 
171
#ifdef GAUSS_FITTER_VERBOSE
 
172
        printf("status = %s\n", gsl_strerror(status));
 
173
        printState_(iter, s);
 
174
 
 
175
        cerr << "f(x)=" << gsl_vector_get(s->x, 0) << " * exp(-(x - " << gsl_vector_get(s->x, 1) << ") ** 2 / 2 / (" << gsl_vector_get(s->x, 2) << ") ** 2)";
 
176
#endif
 
177
 
 
178
        if (status)
 
179
        {
 
180
          break;
 
181
        }
 
182
 
 
183
        status = gsl_multifit_test_delta(s->dx, s->x,
 
184
                                         1e-4, 1e-4);
 
185
      }
 
186
      while (status == GSL_CONTINUE && iter < 500);
 
187
 
 
188
      if (status != GSL_SUCCESS)
 
189
      {
 
190
        gsl_rng_free(r);
 
191
        gsl_multifit_fdfsolver_free(s);
 
192
 
 
193
        throw Exception::UnableToFit(__FILE__, __LINE__, __PRETTY_FUNCTION__, "UnableToFit-GaussFitter", "Could not fit the gaussian to the data");
 
194
      }
 
195
 
 
196
      // write the result in a GaussFitResult struct
 
197
      GaussFitResult result;
 
198
      result.A = gsl_vector_get(s->x, 0);
 
199
      result.x0 = gsl_vector_get(s->x, 1);
 
200
      result.sigma = gsl_vector_get(s->x, 2);
 
201
 
 
202
      // build a formula with the fitted parameters for gnuplot
 
203
      stringstream formula;
 
204
      formula << "f(x)=" << result.A << " * exp(-(x - " << result.x0 << ") ** 2 / 2 / (" << result.sigma << ") ** 2)";
 
205
      gnuplot_formula_ = formula.str();
 
206
#ifdef GAUSS_FITTER_VERBOSE
 
207
      cout << gnuplot_formula_ << endl;
 
208
#endif
 
209
 
 
210
      gsl_rng_free(r);
 
211
      gsl_multifit_fdfsolver_free(s);
 
212
 
 
213
      return result;
 
214
    }
 
215
 
 
216
  }   //namespace Math
209
217
} // namespace OpenMS
210