00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_GAUSSTRACEFITTER_H
00036 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_GAUSSTRACEFITTER_H
00037
00038 #include <sstream>
00039
00040 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/TraceFitter.h>
00041
00042 namespace OpenMS
00043 {
00044
00052 template <typename PeakType>
00053 class GaussTraceFitter :
00054 public TraceFitter<PeakType>
00055 {
00056 public:
00057 GaussTraceFitter()
00058 {
00059
00060 }
00061
00062 GaussTraceFitter(const GaussTraceFitter & other) :
00063 TraceFitter<PeakType>(other)
00064 {
00065 this->height_ = other.height_;
00066 this->x0_ = other.x0_;
00067 this->sigma_ = other.sigma_;
00068
00069 updateMembers_();
00070 }
00071
00072 GaussTraceFitter & operator=(const GaussTraceFitter & source)
00073 {
00074 TraceFitter<PeakType>::operator=(source);
00075
00076 this->height_ = source.height_;
00077 this->x0_ = source.x0_;
00078 this->sigma_ = source.sigma_;
00079
00080 updateMembers_();
00081
00082 return *this;
00083 }
00084
00085 virtual ~GaussTraceFitter()
00086 {
00087 }
00088
00089
00090 void fit(FeatureFinderAlgorithmPickedHelperStructs::MassTraces<PeakType> & traces)
00091 {
00092 LOG_DEBUG << "Traces length: " << traces.size() << std::endl;
00093 setInitialParameters_(traces);
00094
00095 double x_init[NUM_PARAMS_] = {height_, x0_, sigma_};
00096
00097 Size num_params = NUM_PARAMS_;
00098
00099 TraceFitter<PeakType>::optimize_(traces, num_params, x_init,
00100 &(GaussTraceFitter<PeakType>::residual_),
00101 &(GaussTraceFitter<PeakType>::jacobian_),
00102 &(GaussTraceFitter<PeakType>::evaluate_));
00103 }
00104
00105 DoubleReal getLowerRTBound() const
00106 {
00107 return x0_ - 2.5 * sigma_;
00108 }
00109
00110 DoubleReal getUpperRTBound() const
00111 {
00112 return x0_ + 2.5 * sigma_;
00113 }
00114
00115 DoubleReal getHeight() const
00116 {
00117 return height_;
00118 }
00119
00120 DoubleReal getCenter() const
00121 {
00122 return x0_;
00123 }
00124
00125 DoubleReal getFWHM() const
00126 {
00127 return 2.0 * sigma_;
00128 }
00129
00133 DoubleReal getSigma() const
00134 {
00135 return sigma_;
00136 }
00137
00138 bool checkMaximalRTSpan(const DoubleReal max_rt_span)
00139 {
00140 return 5.0 * sigma_ > max_rt_span * region_rt_span_;
00141 }
00142
00143 bool checkMinimalRTSpan(const std::pair<DoubleReal, DoubleReal> & rt_bounds, const DoubleReal min_rt_span)
00144 {
00145 return (rt_bounds.second - rt_bounds.first) < (min_rt_span * 5.0 * sigma_);
00146 }
00147
00148 DoubleReal computeTheoretical(const FeatureFinderAlgorithmPickedHelperStructs::MassTrace<PeakType> & trace, Size k)
00149 {
00150 return trace.theoretical_int * height_ * exp(-0.5 * pow(trace.peaks[k].first - x0_, 2) / pow(sigma_, 2));
00151 }
00152
00153 DoubleReal getFeatureIntensityContribution()
00154 {
00155 return 2.5 * height_ * sigma_;
00156 }
00157
00158 String getGnuplotFormula(FeatureFinderAlgorithmPickedHelperStructs::MassTrace<PeakType> const & trace, const char function_name, const DoubleReal baseline, const DoubleReal rt_shift)
00159 {
00160 std::stringstream s;
00161 s << String(function_name) << "(x)= " << baseline << " + ";
00162 s << (trace.theoretical_int * height_) << " * exp(-0.5*(x-" << (rt_shift + x0_) << ")**2/(" << sigma_ << ")**2)";
00163 return String(s.str());
00164 }
00165
00166 protected:
00167 DoubleReal sigma_;
00168 DoubleReal x0_;
00169 DoubleReal height_;
00170 DoubleReal region_rt_span_;
00171
00172 static const Size NUM_PARAMS_ = 3;
00173
00174 void getOptimizedParameters_(gsl_multifit_fdfsolver * s)
00175 {
00176 height_ = gsl_vector_get(s->x, 0);
00177 x0_ = gsl_vector_get(s->x, 1);
00178 sigma_ = std::fabs(gsl_vector_get(s->x, 2));
00179 }
00180
00181 static Int residual_(const gsl_vector * param, void * data, gsl_vector * f)
00182 {
00183 FeatureFinderAlgorithmPickedHelperStructs::MassTraces<PeakType> * traces = static_cast<FeatureFinderAlgorithmPickedHelperStructs::MassTraces<PeakType> *>(data);
00184 double height = gsl_vector_get(param, 0);
00185 double x0 = gsl_vector_get(param, 1);
00186 double sig = gsl_vector_get(param, 2);
00187 double c_fac = -0.5 / pow(sig, 2);
00188
00189 Size count = 0;
00190 for (Size t = 0; t < traces->size(); ++t)
00191 {
00192 const FeatureFinderAlgorithmPickedHelperStructs::MassTrace<PeakType> & trace = (*traces)[t];
00193 for (Size i = 0; i < trace.peaks.size(); ++i)
00194 {
00195 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());
00196 ++count;
00197 }
00198 }
00199 return GSL_SUCCESS;
00200 }
00201
00202 static Int jacobian_(const gsl_vector * param, void * data, gsl_matrix * J)
00203 {
00204 FeatureFinderAlgorithmPickedHelperStructs::MassTraces<PeakType> * traces = static_cast<FeatureFinderAlgorithmPickedHelperStructs::MassTraces<PeakType> *>(data);
00205 double height = gsl_vector_get(param, 0);
00206 double x0 = gsl_vector_get(param, 1);
00207 double sig = gsl_vector_get(param, 2);
00208 double sig_sq = pow(sig, 2);
00209 double sig_3 = pow(sig, 3);
00210 double c_fac = -0.5 / sig_sq;
00211
00212 Size count = 0;
00213 for (Size t = 0; t < traces->size(); ++t)
00214 {
00215 const FeatureFinderAlgorithmPickedHelperStructs::MassTrace<PeakType> & trace = (*traces)[t];
00216 for (Size i = 0; i < trace.peaks.size(); ++i)
00217 {
00218 DoubleReal rt = trace.peaks[i].first;
00219 DoubleReal e = exp(c_fac * pow(rt - x0, 2));
00220 gsl_matrix_set(J, count, 0, trace.theoretical_int * e);
00221 gsl_matrix_set(J, count, 1, trace.theoretical_int * height * e * (rt - x0) / sig_sq);
00222 gsl_matrix_set(J, count, 2, 0.125 * trace.theoretical_int * height * e * pow(rt - x0, 2) / sig_3);
00223 ++count;
00224 }
00225 }
00226 return GSL_SUCCESS;
00227 }
00228
00229 static Int evaluate_(const gsl_vector * param, void * data, gsl_vector * f, gsl_matrix * J)
00230 {
00231 residual_(param, data, f);
00232 jacobian_(param, data, J);
00233 return GSL_SUCCESS;
00234 }
00235
00236 void setInitialParameters_(FeatureFinderAlgorithmPickedHelperStructs::MassTraces<PeakType> & traces)
00237 {
00238 LOG_DEBUG << "GaussTraceFitter->setInitialParameters(..)" << std::endl;
00239 LOG_DEBUG << "Traces length: " << traces.size() << std::endl;
00240 LOG_DEBUG << "Max trace: " << traces.max_trace << std::endl;
00241
00242
00243 height_ = traces[traces.max_trace].max_peak->getIntensity() - traces.baseline;
00244 LOG_DEBUG << "height: " << height_ << std::endl;
00245 x0_ = traces[traces.max_trace].max_rt;
00246 LOG_DEBUG << "x0: " << x0_ << std::endl;
00247 region_rt_span_ = traces[traces.max_trace].peaks.back().first - traces[traces.max_trace].peaks[0].first;
00248 LOG_DEBUG << "region_rt_span_: " << region_rt_span_ << std::endl;
00249 sigma_ = region_rt_span_ / 20.0;
00250 LOG_DEBUG << "sigma_: " << sigma_ << std::endl;
00251 }
00252
00253 virtual void updateMembers_()
00254 {
00255 TraceFitter<PeakType>::updateMembers_();
00256 }
00257
00258 void printState_(SignedSize iter, gsl_multifit_fdfsolver * s)
00259 {
00260 LOG_DEBUG << "iter " << iter << ": " <<
00261 "height: " << gsl_vector_get(s->x, 0) << " " <<
00262 "x0: " << gsl_vector_get(s->x, 1) << " " <<
00263 "sigma: " << std::fabs(gsl_vector_get(s->x, 2)) << " " <<
00264 "|f(x)| = " << gsl_blas_dnrm2(s->f) << std::endl;
00265 }
00266
00267 };
00268
00269 }
00270
00271 #endif // #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKEDTRACEFITTERGAUSS_H