Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
GaussTraceFitter.h
Go to the documentation of this file.
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.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Stephan Aiche$
32 // $Authors: Stephan Aiche, Marc Sturm$
33 // --------------------------------------------------------------------------
34 
35 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_GAUSSTRACEFITTER_H
36 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_GAUSSTRACEFITTER_H
37 
38 #include <sstream>
39 
41 
42 namespace OpenMS
43 {
44 
52  template <typename PeakType>
54  public TraceFitter<PeakType>
55  {
56 public:
58  {
59  //setName("GaussTraceFitter");
60  }
61 
63  TraceFitter<PeakType>(other)
64  {
65  this->height_ = other.height_;
66  this->x0_ = other.x0_;
67  this->sigma_ = other.sigma_;
68 
70  }
71 
73  {
75 
76  this->height_ = source.height_;
77  this->x0_ = source.x0_;
78  this->sigma_ = source.sigma_;
79 
81 
82  return *this;
83  }
84 
86  {
87  }
88 
89  // override important methods
91  {
92  LOG_DEBUG << "Traces length: " << traces.size() << std::endl;
93  setInitialParameters_(traces);
94 
95  double x_init[NUM_PARAMS_] = {height_, x0_, sigma_};
96 
97  Size num_params = NUM_PARAMS_;
98 
99  TraceFitter<PeakType>::optimize_(traces, num_params, x_init,
103  }
104 
106  {
107  return x0_ - 2.5 * sigma_;
108  }
109 
111  {
112  return x0_ + 2.5 * sigma_;
113  }
114 
116  {
117  return height_;
118  }
119 
121  {
122  return x0_;
123  }
124 
126  {
127  return 2.0 * sigma_;
128  }
129 
134  {
135  return sigma_;
136  }
137 
138  bool checkMaximalRTSpan(const DoubleReal max_rt_span)
139  {
140  return 5.0 * sigma_ > max_rt_span * region_rt_span_;
141  }
142 
143  bool checkMinimalRTSpan(const std::pair<DoubleReal, DoubleReal> & rt_bounds, const DoubleReal min_rt_span)
144  {
145  return (rt_bounds.second - rt_bounds.first) < (min_rt_span * 5.0 * sigma_);
146  }
147 
149  {
150  return trace.theoretical_int * height_ * exp(-0.5 * pow(trace.peaks[k].first - x0_, 2) / pow(sigma_, 2));
151  }
152 
154  {
155  return 2.5 * height_ * sigma_;
156  }
157 
158  String getGnuplotFormula(FeatureFinderAlgorithmPickedHelperStructs::MassTrace<PeakType> const & trace, const char function_name, const DoubleReal baseline, const DoubleReal rt_shift)
159  {
160  std::stringstream s;
161  s << String(function_name) << "(x)= " << baseline << " + ";
162  s << (trace.theoretical_int * height_) << " * exp(-0.5*(x-" << (rt_shift + x0_) << ")**2/(" << sigma_ << ")**2)";
163  return String(s.str());
164  }
165 
166 protected:
171 
172  static const Size NUM_PARAMS_ = 3;
173 
174  void getOptimizedParameters_(gsl_multifit_fdfsolver * s)
175  {
176  height_ = gsl_vector_get(s->x, 0);
177  x0_ = gsl_vector_get(s->x, 1);
178  sigma_ = std::fabs(gsl_vector_get(s->x, 2));
179  }
180 
181  static Int residual_(const gsl_vector * param, void * data, gsl_vector * f)
182  {
184  double height = gsl_vector_get(param, 0);
185  double x0 = gsl_vector_get(param, 1);
186  double sig = gsl_vector_get(param, 2);
187  double c_fac = -0.5 / pow(sig, 2);
188 
189  Size count = 0;
190  for (Size t = 0; t < traces->size(); ++t)
191  {
193  for (Size i = 0; i < trace.peaks.size(); ++i)
194  {
195  gsl_vector_set(f, count, traces->baseline + trace.theoretical_int * height * exp(c_fac * pow(trace.peaks[i].first - x0, 2)) - trace.peaks[i].second->getIntensity());
196  ++count;
197  }
198  }
199  return GSL_SUCCESS;
200  }
201 
202  static Int jacobian_(const gsl_vector * param, void * data, gsl_matrix * J)
203  {
205  double height = gsl_vector_get(param, 0);
206  double x0 = gsl_vector_get(param, 1);
207  double sig = gsl_vector_get(param, 2);
208  double sig_sq = pow(sig, 2);
209  double sig_3 = pow(sig, 3);
210  double c_fac = -0.5 / sig_sq;
211 
212  Size count = 0;
213  for (Size t = 0; t < traces->size(); ++t)
214  {
216  for (Size i = 0; i < trace.peaks.size(); ++i)
217  {
218  DoubleReal rt = trace.peaks[i].first;
219  DoubleReal e = exp(c_fac * pow(rt - x0, 2));
220  gsl_matrix_set(J, count, 0, trace.theoretical_int * e);
221  gsl_matrix_set(J, count, 1, trace.theoretical_int * height * e * (rt - x0) / sig_sq);
222  gsl_matrix_set(J, count, 2, 0.125 * trace.theoretical_int * height * e * pow(rt - x0, 2) / sig_3);
223  ++count;
224  }
225  }
226  return GSL_SUCCESS;
227  }
228 
229  static Int evaluate_(const gsl_vector * param, void * data, gsl_vector * f, gsl_matrix * J)
230  {
231  residual_(param, data, f);
232  jacobian_(param, data, J);
233  return GSL_SUCCESS;
234  }
235 
237  {
238  LOG_DEBUG << "GaussTraceFitter->setInitialParameters(..)" << std::endl;
239  LOG_DEBUG << "Traces length: " << traces.size() << std::endl;
240  LOG_DEBUG << "Max trace: " << traces.max_trace << std::endl;
241 
242  // initial values for externals
243  height_ = traces[traces.max_trace].max_peak->getIntensity() - traces.baseline;
244  LOG_DEBUG << "height: " << height_ << std::endl;
245  x0_ = traces[traces.max_trace].max_rt;
246  LOG_DEBUG << "x0: " << x0_ << std::endl;
247  region_rt_span_ = traces[traces.max_trace].peaks.back().first - traces[traces.max_trace].peaks[0].first;
248  LOG_DEBUG << "region_rt_span_: " << region_rt_span_ << std::endl;
249  sigma_ = region_rt_span_ / 20.0;
250  LOG_DEBUG << "sigma_: " << sigma_ << std::endl;
251  }
252 
253  virtual void updateMembers_()
254  {
256  }
257 
258  void printState_(SignedSize iter, gsl_multifit_fdfsolver * s)
259  {
260  LOG_DEBUG << "iter " << iter << ": " <<
261  "height: " << gsl_vector_get(s->x, 0) << " " <<
262  "x0: " << gsl_vector_get(s->x, 1) << " " <<
263  "sigma: " << std::fabs(gsl_vector_get(s->x, 2)) << " " <<
264  "|f(x)| = " << gsl_blas_dnrm2(s->f) << std::endl;
265  }
266 
267  };
268 
269 } // namespace OpenMS
270 
271 #endif // #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKEDTRACEFITTERGAUSS_H
const double k
DoubleReal getHeight() const
Definition: GaussTraceFitter.h:115
A more convenient string class.
Definition: String.h:56
void getOptimizedParameters_(gsl_multifit_fdfsolver *s)
Definition: GaussTraceFitter.h:174
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:55
Size max_trace
Maximum intensity trace.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:269
void optimize_(FeatureFinderAlgorithmPickedHelperStructs::MassTraces< PeakType > &traces, const Size num_params, double x_init[], Int(*residual)(const gsl_vector *x, void *params, gsl_vector *f), Int(*jacobian)(const gsl_vector *x, void *params, gsl_matrix *J), Int(*evaluate)(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J))
Definition: TraceFitter.h:199
virtual void updateMembers_()
This method is used to update extra member variables at the end of the setParameters() method...
Definition: TraceFitter.h:182
DoubleReal x0_
Definition: GaussTraceFitter.h:168
Helper struct for a collection of mass traces used in FeatureFinderAlgorithmPicked.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:153
Abstract fitter for RT profile fitting.
Definition: TraceFitter.h:62
virtual void updateMembers_()
This method is used to update extra member variables at the end of the setParameters() method...
Definition: GaussTraceFitter.h:253
bool checkMinimalRTSpan(const std::pair< DoubleReal, DoubleReal > &rt_bounds, const DoubleReal min_rt_span)
Definition: GaussTraceFitter.h:143
GaussTraceFitter(const GaussTraceFitter &other)
Definition: GaussTraceFitter.h:62
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:151
DoubleReal sigma_
Definition: GaussTraceFitter.h:167
GaussTraceFitter()
Definition: GaussTraceFitter.h:57
#define LOG_DEBUG
Macro for general debugging information.
Definition: LogStream.h:459
DoubleReal computeTheoretical(const FeatureFinderAlgorithmPickedHelperStructs::MassTrace< PeakType > &trace, Size k)
Definition: GaussTraceFitter.h:148
Fitter for RT profiles using a gaussian background model.
Definition: GaussTraceFitter.h:53
virtual TraceFitter & operator=(const TraceFitter &source)
assignment operator
Definition: TraceFitter.h:86
DoubleReal getSigma() const
Returns the sigma of the fitted gaussian model.
Definition: GaussTraceFitter.h:133
DoubleReal getCenter() const
Definition: GaussTraceFitter.h:120
DoubleReal theoretical_int
Theoretical intensity value (scaled to [0,1])
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:91
String getGnuplotFormula(FeatureFinderAlgorithmPickedHelperStructs::MassTrace< PeakType > const &trace, const char function_name, const DoubleReal baseline, const DoubleReal rt_shift)
Definition: GaussTraceFitter.h:158
DoubleReal getFeatureIntensityContribution()
Definition: GaussTraceFitter.h:153
GaussTraceFitter & operator=(const GaussTraceFitter &source)
Definition: GaussTraceFitter.h:72
void printState_(SignedSize iter, gsl_multifit_fdfsolver *s)
Definition: GaussTraceFitter.h:258
DoubleReal baseline
Estimated baseline in the region of the feature (used for the fit)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:271
DoubleReal getFWHM() const
Definition: GaussTraceFitter.h:125
void fit(FeatureFinderAlgorithmPickedHelperStructs::MassTraces< PeakType > &traces)
Definition: GaussTraceFitter.h:90
virtual ~GaussTraceFitter()
Definition: GaussTraceFitter.h:85
DoubleReal height_
Definition: GaussTraceFitter.h:169
DoubleReal region_rt_span_
Definition: GaussTraceFitter.h:170
void setInitialParameters_(FeatureFinderAlgorithmPickedHelperStructs::MassTraces< PeakType > &traces)
Definition: GaussTraceFitter.h:236
static Int evaluate_(const gsl_vector *param, void *data, gsl_vector *f, gsl_matrix *J)
Definition: GaussTraceFitter.h:229
static Int residual_(const gsl_vector *param, void *data, gsl_vector *f)
Definition: GaussTraceFitter.h:181
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
std::vector< std::pair< DoubleReal, const PeakType * > > peaks
Contained peaks (pair of RT and pointer to peak)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:94
DoubleReal getUpperRTBound() const
Definition: GaussTraceFitter.h:110
static const Size NUM_PARAMS_
Definition: GaussTraceFitter.h:172
DoubleReal getLowerRTBound() const
Definition: GaussTraceFitter.h:105
bool checkMaximalRTSpan(const DoubleReal max_rt_span)
Definition: GaussTraceFitter.h:138
int Int
Signed integer type.
Definition: Types.h:100
static Int jacobian_(const gsl_vector *param, void *data, gsl_matrix *J)
Definition: GaussTraceFitter.h:202
Helper struct for mass traces used in FeatureFinderAlgorithmPicked.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:83

OpenMS / TOPP release 1.11.1 Documentation generated on Thu Nov 14 2013 11:19:14 using doxygen 1.8.5