35 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_EGHTRACEFITTER_H
36 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_EGHTRACEFITTER_H
57 template <
class PeakType>
159 double rt = trace.
peaks[
k].first;
160 double t_diff, t_diff2, denominator = 0.0;
164 t_diff2 = t_diff * t_diff;
168 if (denominator > 0.0)
185 return bounds.second - bounds.first;
191 s <<
String(function_name) <<
"(x)= " << baseline <<
" + ";
222 std::pair<DoubleReal, DoubleReal> bounds;
229 s1 = (-1 * (L *
tau_) / 2) + s;
230 s2 = (-1 * (L *
tau_) / 2) - s;
233 bounds.first =
apex_rt_ + std::min(s1, s2);
235 bounds.second =
apex_rt_ + std::max(s1, s2);
242 height_ = gsl_vector_get(fdfsolver->x, 0);
243 apex_rt_ = gsl_vector_get(fdfsolver->x, 1);
245 tau_ = gsl_vector_get(fdfsolver->x, 3);
254 static Int residual_(
const gsl_vector * param,
void * data, gsl_vector * f)
258 double H = gsl_vector_get(param, 0);
259 double tR = gsl_vector_get(param, 1);
260 double sigma_square = gsl_vector_get(param, 2);
261 double tau = gsl_vector_get(param, 3);
263 double t_diff, t_diff2, denominator = 0.0;
268 for (
Size t = 0; t < traces->size(); ++t)
271 for (
Size i = 0; i < trace.
peaks.size(); ++i)
276 t_diff2 = t_diff * t_diff;
278 denominator = 2 * sigma_square + tau * t_diff;
280 if (denominator > 0.0)
289 gsl_vector_set(f, count, (fegh - trace.
peaks[i].second->getIntensity()));
296 static Int jacobian_(
const gsl_vector * param,
void * data, gsl_matrix * J)
300 double H = gsl_vector_get(param, 0);
301 double tR = gsl_vector_get(param, 1);
302 double sigma_square = gsl_vector_get(param, 2);
303 double tau = gsl_vector_get(param, 3);
305 double derivative_H, derivative_tR, derivative_sigma_square, derivative_tau = 0.0;
306 double t_diff, t_diff2, exp1, denominator = 0.0;
309 for (
Size t = 0; t < traces->size(); ++t)
312 for (
Size i = 0; i < trace.
peaks.size(); ++i)
317 t_diff2 = t_diff * t_diff;
319 denominator = 2 * sigma_square + tau * t_diff;
323 exp1 = exp(-t_diff2 / denominator);
329 derivative_tR = trace.
theoretical_int * H * exp1 * (((4 * sigma_square + tau * t_diff) * t_diff) / (denominator * denominator));
332 derivative_sigma_square = trace.
theoretical_int * H * exp1 * ((2 * t_diff2) / (denominator * denominator));
335 derivative_tau = trace.
theoretical_int * H * exp1 * ((t_diff * t_diff2) / (denominator * denominator));
341 derivative_sigma_square = 0.0;
342 derivative_tau = 0.0;
346 gsl_matrix_set(J, count, 0, derivative_H);
347 gsl_matrix_set(J, count, 1, derivative_tR);
348 gsl_matrix_set(J, count, 2, derivative_sigma_square);
349 gsl_matrix_set(J, count, 3, derivative_tau);
357 static Int evaluate_(
const gsl_vector * param,
void * data, gsl_vector * f, gsl_matrix * J)
366 LOG_DEBUG <<
"EGHTraceFitter->setInitialParameters(..)" << std::endl;
367 LOG_DEBUG <<
"Traces length: " << traces.size() << std::endl;
381 for (
Size i = 1; i < traces[traces.
max_trace].peaks.size(); ++i)
385 max_peak = traces[traces.
max_trace].peaks[i].second;
391 LOG_DEBUG <<
"max_pos: " << max_pos << std::endl;
392 if (traces[traces.
max_trace].peaks.size() < 3)
398 Size filter_max_pos = traces[traces.
max_trace].peaks.size() - 2;
403 if ((max_pos < 2) || (max_pos + 2 >= traces[traces.
max_trace].peaks.size()))
406 smoothed_height = traces[traces.
max_trace].peaks[max_pos].second->getIntensity();
411 smoothed_height = (traces[traces.
max_trace].peaks[max_pos - 2].second->getIntensity()
412 + traces[traces.
max_trace].peaks[max_pos - 1].second->getIntensity()
413 + traces[traces.
max_trace].peaks[max_pos].second->getIntensity()
414 + traces[traces.
max_trace].peaks[max_pos + 1].second->getIntensity()
415 + traces[traces.
max_trace].peaks[max_pos + 2].second->getIntensity()) / 5.0;
421 while (i > 2 && i < filter_max_pos)
425 + traces[traces.
max_trace].peaks[i - 1].second->getIntensity()
426 + traces[traces.
max_trace].peaks[i].second->getIntensity()
427 + traces[traces.
max_trace].peaks[i + 1].second->getIntensity()
428 + traces[traces.
max_trace].peaks[i + 2].second->getIntensity()) / 5.0;
430 if (smoothed / smoothed_height < 0.5)
break;
433 LOG_DEBUG <<
"Left alpha at " << i <<
" with " << traces[traces.
max_trace].peaks[i].first << std::endl;
437 while (i < filter_max_pos && i > 2)
440 + traces[traces.
max_trace].peaks[i - 1].second->getIntensity()
441 + traces[traces.
max_trace].peaks[i].second->getIntensity()
442 + traces[traces.
max_trace].peaks[i + 1].second->getIntensity()
443 + traces[traces.
max_trace].peaks[i + 2].second->getIntensity()) / 5.0;
445 if (smoothed / smoothed_height < 0.5)
break;
448 LOG_DEBUG <<
"Right alpha at " << i <<
" with " << traces[traces.
max_trace].peaks[i].first << std::endl;
455 double log_alpha = log(0.5);
457 tau_ = (-1 / log_alpha) * (B - A);
471 <<
"height: " << gsl_vector_get(s->x, 0) <<
" "
472 <<
"apex_rt: " << gsl_vector_get(s->x, 1) <<
" "
473 <<
"sigma_square: " << gsl_vector_get(s->x, 2) <<
" "
474 <<
"tau: " << gsl_vector_get(s->x, 3) <<
" "
475 <<
"|f(x)| = " << gsl_blas_dnrm2(s->f) << std::endl;
482 #endif // #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKEDTRACEFITTERGAUSS_H
DoubleReal computeTheoretical(const FeatureFinderAlgorithmPickedHelperStructs::MassTrace< PeakType > &trace, Size k)
Definition: EGHTraceFitter.h:157
A more convenient string class.
Definition: String.h:56
EGHTraceFitter()
Definition: EGHTraceFitter.h:62
DoubleReal tau_
Definition: EGHTraceFitter.h:206
void getOptimizedParameters_(gsl_multifit_fdfsolver *fdfsolver)
Definition: EGHTraceFitter.h:240
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:55
DoubleReal getTau() const
Definition: EGHTraceFitter.h:122
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
void fit(FeatureFinderAlgorithmPickedHelperStructs::MassTraces< PeakType > &traces)
Definition: EGHTraceFitter.h:103
IntensityType getIntensity() const
Definition: Peak2D.h:161
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
static const Size NUM_PARAMS_
Definition: EGHTraceFitter.h:213
DoubleReal getFWHM() const
Definition: EGHTraceFitter.h:181
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:151
virtual DoubleReal getFeatureIntensityContribution()
Definition: EGHTraceFitter.h:176
std::pair< DoubleReal, DoubleReal > getAlphaBoundaries_(const DoubleReal alpha) const
Return an ordered pair of the positions where the EGH reaches a height of alpha * height of the EGH...
Definition: EGHTraceFitter.h:220
DoubleReal getSigmaSquare() const
Definition: EGHTraceFitter.h:137
#define LOG_DEBUG
Macro for general debugging information.
Definition: LogStream.h:459
DoubleReal height_
Definition: EGHTraceFitter.h:203
static Int jacobian_(const gsl_vector *param, void *data, gsl_matrix *J)
Definition: EGHTraceFitter.h:296
virtual TraceFitter & operator=(const TraceFitter &source)
assignment operator
Definition: TraceFitter.h:86
virtual ~EGHTraceFitter()
Definition: EGHTraceFitter.h:98
DoubleReal theoretical_int
Theoretical intensity value (scaled to [0,1])
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:91
std::pair< DoubleReal, DoubleReal > sigma_5_bound_
Definition: EGHTraceFitter.h:208
EGHTraceFitter & operator=(const EGHTraceFitter &source)
Definition: EGHTraceFitter.h:81
DoubleReal getUpperRTBound() const
Definition: EGHTraceFitter.h:127
std::pair< DoubleReal, DoubleReal > fwhm_bound_
Definition: EGHTraceFitter.h:209
DoubleReal getCenter() const
Definition: EGHTraceFitter.h:142
virtual void updateMembers_()
This method is used to update extra member variables at the end of the setParameters() method...
Definition: EGHTraceFitter.h:463
DoubleReal baseline
Estimated baseline in the region of the feature (used for the fit)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:271
virtual bool checkMinimalRTSpan(const std::pair< DoubleReal, DoubleReal > &rt_bounds, const DoubleReal min_rt_span)
Definition: EGHTraceFitter.h:152
DoubleReal region_rt_span_
Definition: EGHTraceFitter.h:211
DoubleReal apex_rt_
Definition: EGHTraceFitter.h:202
void printState_(SignedSize iter, gsl_multifit_fdfsolver *s)
Definition: EGHTraceFitter.h:468
DoubleReal getHeight() const
Definition: EGHTraceFitter.h:132
DoubleReal sigma_square_
Definition: EGHTraceFitter.h:205
static Int residual_(const gsl_vector *param, void *data, gsl_vector *f)
Definition: EGHTraceFitter.h:254
bool checkMaximalRTSpan(const DoubleReal max_rt_span)
Definition: EGHTraceFitter.h:147
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
EGHTraceFitter(const EGHTraceFitter &other)
Definition: EGHTraceFitter.h:67
virtual String getGnuplotFormula(FeatureFinderAlgorithmPickedHelperStructs::MassTrace< PeakType > const &trace, const char function_name, const DoubleReal baseline, const DoubleReal rt_shift)
Definition: EGHTraceFitter.h:188
A RT Profile fitter using an Exponential Gaussian Hybrid background model.
Definition: EGHTraceFitter.h:58
void setInitialParameters_(FeatureFinderAlgorithmPickedHelperStructs::MassTraces< PeakType > &traces)
Definition: EGHTraceFitter.h:364
DoubleReal getLowerRTBound() const
Definition: EGHTraceFitter.h:117
int Int
Signed integer type.
Definition: Types.h:100
static Int evaluate_(const gsl_vector *param, void *data, gsl_vector *f, gsl_matrix *J)
Definition: EGHTraceFitter.h:357
Helper struct for mass traces used in FeatureFinderAlgorithmPicked.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:83