Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
TwoDOptimization.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-2015.
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: Alexandra Zerck $
32 // $Authors: $
33 // --------------------------------------------------------------------------
34 
35 #ifndef OPENMS_TRANSFORMATIONS_RAW2PEAK_TWODOPTIMIZATION_H
36 #define OPENMS_TRANSFORMATIONS_RAW2PEAK_TWODOPTIMIZATION_H
37 
38 //#define DEBUG_2D
39 #undef DEBUG_2D
40 
41 #ifdef DEBUG_2D
42 #include <iostream>
43 #include <fstream>
44 #endif
45 
46 #include <vector>
47 #include <utility>
48 #include <cmath>
49 #include <set>
50 
60 
61 #include <Eigen/Core>
62 #include <unsupported/Eigen/NonLinearOptimization>
63 
64 #ifndef OPENMS_SYSTEM_STOPWATCH_H
65 #endif
66 
67 #include <boost/math/special_functions/acosh.hpp>
70 
71 namespace OpenMS
72 {
73 
88  class OPENMS_DLLAPI TwoDOptimization :
89  public DefaultParamHandler
90  {
91 public:
92 
95 
98 
100  virtual ~TwoDOptimization(){}
101 
103  TwoDOptimization& operator=(const TwoDOptimization& opt);
104 
105 
107  inline double getMZTolerance() const {return tolerance_mz_; }
109  inline void setMZTolerance(double tolerance_mz)
110  {
111  tolerance_mz_ = tolerance_mz;
112  param_.setValue("2d:tolerance_mz", tolerance_mz);
113  }
114 
116  inline double getMaxPeakDistance() const {return max_peak_distance_; }
118  inline void setMaxPeakDistance(double max_peak_distance)
119  {
120  max_peak_distance_ = max_peak_distance;
121  param_.setValue("2d:max_peak_distance", max_peak_distance);
122  }
123 
125  inline UInt getMaxIterations() const {return max_iteration_; }
127  inline void setMaxIterations(UInt max_iteration)
128  {
129  max_iteration_ = max_iteration;
130  param_.setValue("iterations", max_iteration);
131  }
132 
134  inline const OptimizationFunctions::PenaltyFactorsIntensity& getPenalties() const {return penalties_; }
137  {
138  penalties_ = penalties;
139  param_.setValue("penalties:position", penalties.pos);
140  param_.setValue("penalties:height", penalties.height);
141  param_.setValue("penalties:left_width", penalties.lWidth);
142  param_.setValue("penalties:right_width", penalties.rWidth);
143  }
144 
161  template <typename InputSpectrumIterator, typename OutputPeakType>
162  void optimize(InputSpectrumIterator first,
163  InputSpectrumIterator last,
164  MSExperiment<OutputPeakType>& ms_exp, bool real2D = true);
165 
166 
167 protected:
169  struct Data
170  {
171  std::vector<std::pair<SignedSize, SignedSize> > signal2D;
172  std::multimap<double, IsotopeCluster>::iterator iso_map_iter;
174  std::map<Int, std::vector<PeakIndex> > matching_peaks;
178  std::vector<double> positions;
179  std::vector<double> signal;
180  };
181 
182  class OPENMS_DLLAPI TwoDOptFunctor
183  {
184  public:
185  int inputs() const { return m_inputs; }
186  int values() const { return m_values; }
187 
188  TwoDOptFunctor(unsigned dimensions, unsigned num_data_points, const TwoDOptimization::Data* data)
189  : m_inputs(dimensions), m_values(num_data_points), m_data(data) {}
190 
191  int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec);
192  // compute Jacobian matrix for the different parameters
193  int df(const Eigen::VectorXd &x, Eigen::MatrixXd &J);
194 
195  private:
196  const int m_inputs, m_values;
198  };
199 
200 
202  std::multimap<double, IsotopeCluster> iso_map_;
203 
205  std::multimap<double, IsotopeCluster>::const_iterator curr_region_;
206 
209 
212 
214  // std::map<Int, std::vector<MSExperiment<>::SpectrumType::Iterator > > matching_peaks_;
215  std::map<Int, std::vector<PeakIndex> > matching_peaks_;
216 
217 
220 
222  bool real_2D_;
223 
224 
227 
228 
233  std::vector<double>::iterator searchInScan_(std::vector<double>::iterator scan_begin,
234  std::vector<double>::iterator scan_end,
235  double current_mz);
236 
238  template <typename InputSpectrumIterator, typename OutputPeakType>
239  void optimizeRegions_(InputSpectrumIterator& first,
240  InputSpectrumIterator& last,
242 
244  template <typename InputSpectrumIterator, typename OutputPeakType>
245  void optimizeRegionsScanwise_(InputSpectrumIterator& first,
246  InputSpectrumIterator& last,
248 
249 
251  template <typename InputSpectrumIterator, typename OutputPeakType>
252  void getRegionEndpoints_(MSExperiment<OutputPeakType>& exp,
253  InputSpectrumIterator& first,
254  InputSpectrumIterator& last,
255  Size iso_map_idx,
256  double noise_level,
258 
260  void findMatchingPeaks_(std::multimap<double, IsotopeCluster>::iterator& it,
261  MSExperiment<>& ms_exp);
262 
264 
266  void updateMembers_();
267  };
268 
269 
270  template <typename InputSpectrumIterator, typename OutputPeakType>
271  void TwoDOptimization::optimize(InputSpectrumIterator first, InputSpectrumIterator last, MSExperiment<OutputPeakType>& ms_exp, bool real2D)
272  {
273  //#define DEBUG_2D
274  //check if the input maps have the same number of spectra
275  if ((UInt)distance(first, last) != ms_exp.size())
276  {
277  throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__, "Error in Two2Optimization: Raw and peak map do not have the same number of spectra");
278  }
279  //do nothing if there are no scans
280  if (ms_exp.empty())
281  {
282  return;
283  }
284  //check if required meta data arrays are present (for each scan)
285  for (Size i = 0; i < ms_exp.size(); ++i)
286  {
287  //check if enough meta data arrays are present
288  if (ms_exp[i].getFloatDataArrays().size() < 6)
289  {
290  throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__, "Error in Two2Optimization: Not enough meta data arrays present (1:area, 5:shape, 3:left width, 4:right width)");
291  }
292  bool area = ms_exp[i].getFloatDataArrays()[1].getName() == "maximumIntensity";
293  bool wleft = ms_exp[i].getFloatDataArrays()[3].getName() == "leftWidth";
294  bool wright = ms_exp[i].getFloatDataArrays()[4].getName() == "rightWidth";
295  bool shape = ms_exp[i].getFloatDataArrays()[5].getName() == "peakShape";
296 
297  if (!area || !wleft || !wright || !shape)
298  {
299  throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__, "Error in Two2Optimization: One or several meta data arrays missing (1:intensity, 5:shape, 3:left width, 4:right width)");
300  }
301  }
302  real_2D_ = real2D;
303  typedef typename InputSpectrumIterator::value_type InputSpectrumType;
304  typedef typename InputSpectrumType::value_type PeakType;
305  typedef MSSpectrum<PeakType> SpectrumType;
306 
307  typename MSExperiment<OutputPeakType>::Iterator ms_exp_it = ms_exp.begin();
308  typename MSExperiment<OutputPeakType>::Iterator ms_exp_it_end = ms_exp.end();
309  if (ms_exp.empty())
310  {
311  return;
312  }
313  // stores the monoisotopic peaks of isotopic clusters
314  std::vector<double> iso_last_scan;
315  std::vector<double> iso_curr_scan;
316  std::vector<std::multimap<double, IsotopeCluster>::iterator> clusters_last_scan;
317  std::vector<std::multimap<double, IsotopeCluster>::iterator> clusters_curr_scan;
318  std::multimap<double, IsotopeCluster>::iterator cluster_iter;
319  double current_rt = ms_exp_it->getRT(), last_rt = 0;
320 
321  // retrieve values for accepted peaks distances
322  max_peak_distance_ = param_.getValue("2d:max_peak_distance");
323  double tolerance_mz = param_.getValue("2d:tolerance_mz");
324 
325  UInt current_charge = 0; // charge state of the current isotopic cluster
326  double mz_in_hash = 0; // used as reference to the current isotopic peak
327 
328  // sweep through scans
329  for (UInt curr_scan = 0; ms_exp_it + curr_scan != ms_exp_it_end; ++curr_scan)
330  {
331  Size nr_peaks_in_scan = (ms_exp_it + curr_scan)->size();
332  if (nr_peaks_in_scan == 0)
333  continue;
334 
335  //last_rt = current_rt;
336  current_rt = (ms_exp_it + curr_scan)->getRT();
337  typename MSExperiment<OutputPeakType>::SpectrumType::Iterator peak_it = (ms_exp_it + curr_scan)->begin();
338 
339  // copy cluster information of least scan
340  iso_last_scan = iso_curr_scan;
341  iso_curr_scan.clear();
342  clusters_last_scan = clusters_curr_scan;
343  clusters_curr_scan.clear();
344 
345 #ifdef DEBUG_2D
346  std::cout << "Next scan with rt: " << current_rt << std::endl;
347  std::cout << "Next scan, rt = " << current_rt << " last_rt: " << last_rt << std::endl;
348  std::cout << "---------------------------------------------------------------------------" << std::endl;
349 #endif
351  s.setRT(current_rt);
352  // check if there were scans in between
353  if (last_rt == 0 || // are we in the first scan
354  ((lower_bound(first, last, s, typename SpectrumType::RTLess()) - 1)->getRT() == last_rt))
355  {
356 
357 
358  for (UInt curr_peak = 0; curr_peak < (ms_exp_it + curr_scan)->size() - 1; ++curr_peak)
359  {
360 
361  // store the m/z of the current peak
362  double curr_mz = (peak_it + curr_peak)->getMZ();
363  double dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - curr_mz;
364 
365  if (dist2nextpeak <= max_peak_distance_) // one single peak without neighbors isn't optimized
366  {
367 #ifdef DEBUG_2D
368  std::cout << "Isotopic pattern found ! " << std::endl;
369  std::cout << "We are at: " << (peak_it + curr_peak)->getMZ() << " " << curr_mz << std::endl;
370 #endif
371  if (!iso_last_scan.empty()) // Did we find any isotopic cluster in the last scan?
372  {
373  std::sort(iso_last_scan.begin(), iso_last_scan.end());
374  // there were some isotopic clusters in the last scan...
375  std::vector<double>::iterator it =
376  searchInScan_(iso_last_scan.begin(), iso_last_scan.end(), curr_mz);
377 
378  double delta_mz = fabs(*it - curr_mz);
379  //std::cout << delta_mz << " "<< tolerance_mz << std::endl;
380  if (delta_mz > tolerance_mz) // check if first peak of last cluster is close enough
381  {
382  mz_in_hash = curr_mz; // update current hash key
383 
384  // create new isotopic cluster
385 // #ifdef DEBUG_2D
386 // std::cout << "Last peak cluster too far, creating new cluster at "<<curr_mz << std::endl;
387 // #endif
388  IsotopeCluster new_cluster;
389  new_cluster.peaks.charge = current_charge;
390  new_cluster.scans.push_back(curr_scan);
391  cluster_iter = iso_map_.insert(std::pair<double, IsotopeCluster>(mz_in_hash, new_cluster));
392 
393  }
394  else
395  {
396 // //#ifdef DEBUG_2D
397 // std::cout << "Found neighbouring peak with distance (m/z) " << delta_mz << std::endl;
398 // //#endif
399  cluster_iter = clusters_last_scan[distance(iso_last_scan.begin(), it)];
400 
401  // check whether this scan is already contained
402  if (find(cluster_iter->second.scans.begin(), cluster_iter->second.scans.end(), curr_scan)
403  == cluster_iter->second.scans.end())
404  {
405  cluster_iter->second.scans.push_back(curr_scan);
406  }
407 
408 // //#ifdef DEBUG_2D
409 // std::cout << "Cluster with " << cluster_iter->second.peaks.size()
410 // << " peaks retrieved." << std::endl;
411 // //#endif
412  }
413 
414  }
415  else // last scan did not contain any isotopic cluster
416  {
417 // //#ifdef DEBUG_2D
418 // std::cout << "Last scan was empty => creating new cluster." << std::endl;
419 // std::cout << "Creating new cluster at m/z: " << curr_mz << std::endl;
420 // //#endif
421 
422  mz_in_hash = curr_mz; // update current hash key
423 
424  // create new isotopic cluster
425  IsotopeCluster new_cluster;
426  new_cluster.peaks.charge = current_charge;
427  new_cluster.scans.push_back(curr_scan);
428  cluster_iter = iso_map_.insert(std::pair<double, IsotopeCluster>(mz_in_hash, new_cluster));
429 
430  }
431 
432 // //#ifdef DEBUG_2D
433 // std::cout << "Storing found peak in current isotopic cluster" << std::endl;
434 // //#endif
435 
436 
437 
438  cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
439 
440  iso_curr_scan.push_back(mz_in_hash);
441  clusters_curr_scan.push_back(cluster_iter);
442  ++curr_peak;
443 
444  cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
445  iso_curr_scan.push_back((peak_it + curr_peak)->getMZ());
446  clusters_curr_scan.push_back(cluster_iter);
447 
448  // check distance to next peak
449  if ((curr_peak + 1) >= nr_peaks_in_scan)
450  break;
451  dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - (peak_it + curr_peak)->getMZ();
452 
453 
454  // loop until end of isotopic pattern in this scan
455  while (dist2nextpeak <= max_peak_distance_
456  && curr_peak < (nr_peaks_in_scan - 1))
457  {
458  cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak + 1)); // save peak in cluster
459  iso_curr_scan.push_back((peak_it + curr_peak + 1)->getMZ());
460  clusters_curr_scan.push_back(cluster_iter);
461  // std::cout << "new entered: "<<(peak_it+curr_peak+1)->getMZ()<<" im while"<<std::endl;
462  ++curr_peak;
463  if (curr_peak >= nr_peaks_in_scan - 1)
464  break;
465  dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - (peak_it + curr_peak)->getMZ(); // get distance to next peak
466 
467 
468  } // end while(...)
469 
470 
471 
472  } // end of if (dist2nextpeak <= max_peak_distance_)
473  else
474  {
475  if (!iso_last_scan.empty()) // Did we find any isotopic cluster in the last scan?
476  {
477  std::sort(iso_last_scan.begin(), iso_last_scan.end());
478  // there were some isotopic clusters in the last scan...
479  std::vector<double>::iterator it =
480  searchInScan_(iso_last_scan.begin(), iso_last_scan.end(), curr_mz);
481 
482  double delta_mz = fabs(*it - curr_mz);
483  // std::cout << delta_mz << " "<< tolerance_mz << std::endl;
484  if (delta_mz > tolerance_mz) // check if first peak of last cluster is close enough
485  {
486  mz_in_hash = curr_mz; // update current hash key
487 
488  // create new isotopic cluster
489 // //#ifdef DEBUG_2D
490 // std::cout << "Last peak cluster too far, creating new cluster at "<<curr_mz << std::endl;
491 // //#endif
492  IsotopeCluster new_cluster;
493  new_cluster.peaks.charge = current_charge;
494  new_cluster.scans.push_back(curr_scan);
495  cluster_iter = iso_map_.insert(std::pair<double, IsotopeCluster>(mz_in_hash, new_cluster));
496 
497  }
498  else
499  {
500 // //#ifdef DEBUG_2D
501 // std::cout << "Found neighbouring peak with distance (m/z) " << delta_mz << std::endl;
502 // //#endif
503  cluster_iter = clusters_last_scan[distance(iso_last_scan.begin(), it)];
504 
505  // check whether this scan is already contained
506  if (find(cluster_iter->second.scans.begin(), cluster_iter->second.scans.end(), curr_scan)
507  == cluster_iter->second.scans.end())
508  {
509  cluster_iter->second.scans.push_back(curr_scan);
510  }
511 
512 // //#ifdef DEBUG_2D
513 // std::cout << "Cluster with " << cluster_iter->second.peaks.size()
514 // << " peaks retrieved." << std::endl;
515 // //#endif
516  }
517 
518  }
519  else // last scan did not contain any isotopic cluster
520  {
521 // //#ifdef DEBUG_2D
522 // std::cout << "Last scan was empty => creating new cluster." << std::endl;
523 // std::cout << "Creating new cluster at m/z: " << curr_mz << std::endl;
524 // //#endif
525 
526  mz_in_hash = curr_mz; // update current hash key
527 
528  // create new isotopic cluster
529  IsotopeCluster new_cluster;
530  new_cluster.peaks.charge = current_charge;
531  new_cluster.scans.push_back(curr_scan);
532  cluster_iter = iso_map_.insert(std::pair<double, IsotopeCluster>(mz_in_hash, new_cluster));
533 
534  }
535 
536 // //#ifdef DEBUG_2D
537 // std::cout << "Storing found peak in current isotopic cluster" << std::endl;
538 // //#endif
539 
540 
541 
542  cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
543 
544  iso_curr_scan.push_back(mz_in_hash);
545  clusters_curr_scan.push_back(cluster_iter);
546 
547 
548  }
549 
550  current_charge = 0; // reset charge
551  } // end for (...)
552  }
553  last_rt = current_rt;
554  }
555  curr_region_ = iso_map_.begin();
556 #ifdef DEBUG_2D
557  std::cout << iso_map_.size() << " isotopic clusters were found ! " << std::endl;
558 #endif
559 
560  if (real_2D_)
561  optimizeRegions_(first, last, ms_exp);
562  else
563  optimizeRegionsScanwise_(first, last, ms_exp);
564  //#undef DEBUG_2D
565  }
566 
567 
568 
569  template <typename InputSpectrumIterator, typename OutputPeakType>
570  void TwoDOptimization::optimizeRegions_(InputSpectrumIterator& first,
571  InputSpectrumIterator& last,
573  {
574  Int counter = 0;
575  // go through the clusters
576  for (std::multimap<double, IsotopeCluster>::iterator it = iso_map_.begin();
577  it != iso_map_.end();
578  ++it)
579  {
580 #ifdef DEBUG_2D
581  std::cout << "element: " << counter << std::endl;
582  std::cout << "mz: " << it->first << std::endl << "rts: ";
583 // for(Size i=0;i<it->second.scans.size();++i) std::cout << it->second.scans[i] << "\n";
584  std::cout << std::endl << "peaks: ";
585  IsotopeCluster::IndexSet::const_iterator iter = it->second.peaks.begin();
586  for (; iter != it->second.peaks.end(); ++iter)
587  std::cout << ms_exp[iter->first].getRT() << " " << (ms_exp[iter->first][iter->second]).getMZ() << std::endl;
588 
589 //for(Size i=0;i<it->second.peaks.size();++i) std::cout << ms_exp[it->first].getRT() << " "<<(ms_exp[it->first][it->second]).getMZ()<<std::endl;
590  std::cout << std::endl << std::endl;
591 
592 #endif
593 
594  // prepare for optimization:
595  // determine the matching peaks
596  matching_peaks_.clear();
597  findMatchingPeaks_(it, ms_exp);
598  TwoDOptimization::Data twoD_data;
599  twoD_data.penalties = penalties_;
600  twoD_data.matching_peaks = matching_peaks_;
601  // and the endpoints of each isotope pattern in the cluster
602  getRegionEndpoints_(ms_exp, first, last, counter, 400, twoD_data);
603 
604  // peaks have to be stored globally
605  twoD_data.iso_map_iter = it;
606 
607  twoD_data.picked_peaks = ms_exp;
608  twoD_data.raw_data_first = first;
609 
610  Size nr_diff_peaks = matching_peaks_.size();
611  twoD_data.total_nr_peaks = it->second.peaks.size();
612 
613  Size nr_parameters = nr_diff_peaks * 3 + twoD_data.total_nr_peaks;
614 
615  // initialize and set parameters for optimization
616  Eigen::VectorXd x_init (nr_parameters);
617  x_init.setZero();
618 
619  std::map<Int, std::vector<PeakIndex> >::iterator m_peaks_it = twoD_data.matching_peaks.begin();
620  Int peak_counter = 0;
621  Int diff_peak_counter = 0;
622  // go through the matching peaks
623  for (; m_peaks_it != twoD_data.matching_peaks.end(); ++m_peaks_it)
624  {
625  double av_mz = 0, av_lw = 0, av_rw = 0, avr_height = 0, height;
626  std::vector<PeakIndex>::iterator iter_iter = (m_peaks_it)->second.begin();
627  for (; iter_iter != m_peaks_it->second.end(); ++iter_iter)
628  {
629  height = ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[1][(iter_iter)->peak]; //(iter_iter)->getPeak(ms_exp).getIntensity();
630  avr_height += height;
631  av_mz += (iter_iter)->getPeak(ms_exp).getMZ() * height;
632  av_lw += ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[3][(iter_iter)->peak] * height; //left width
633  av_rw += ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[4][(iter_iter)->peak] * height; //right width
634  x_init(peak_counter) = height;
635  ++peak_counter;
636  }
637  x_init(twoD_data.total_nr_peaks + 3 * diff_peak_counter) = av_mz / avr_height;
638  x_init(twoD_data.total_nr_peaks + 3 * diff_peak_counter + 1) = av_lw / avr_height;
639  x_init(twoD_data.total_nr_peaks + 3 * diff_peak_counter + 2) = av_rw / avr_height;
640  ++diff_peak_counter;
641  }
642 
643 #ifdef DEBUG_2D
644  std::cout << "----------------------------\n\nstart_value: " << std::endl;
645  for (Size k = 0; k < start_value->size; ++k)
646  {
647  std::cout << x_init(k) << std::endl;
648  }
649 #endif
650  Int num_positions = 0;
651  for (Size i = 0; i < twoD_data.signal2D.size(); i += 2)
652  {
653  num_positions += (twoD_data.signal2D[i + 1].second - twoD_data.signal2D[i].second + 1);
654 #ifdef DEBUG_2D
655  std::cout << twoD_data.signal2D[i + 1].second << " - " << twoD_data.signal2D[i].second << " +1 " << std::endl;
656 #endif
657 
658  }
659 #ifdef DEBUG_2D
660  std::cout << "num_positions : " << num_positions << std::endl;
661 #endif
662 
663  TwoDOptFunctor functor (nr_parameters, std::max(num_positions + 1, (Int)(nr_parameters)), &twoD_data);
664  Eigen::LevenbergMarquardt<TwoDOptFunctor> lmSolver (functor);
665  Eigen::LevenbergMarquardtSpace::Status status = lmSolver.minimize(x_init);
666 
667  //the states are poorly documented. after checking the source, we believe that
668  //all states except NotStarted, Running and ImproperInputParameters are good
669  //termination states.
670  if (status <= Eigen::LevenbergMarquardtSpace::ImproperInputParameters)
671  {
672  throw Exception::UnableToFit(__FILE__, __LINE__, __PRETTY_FUNCTION__, "UnableToFit-TwoDOptimization:", "Could not fit the data: Error " + String(status));
673  }
674 
675  Int peak_idx = 0;
676  std::map<Int, std::vector<PeakIndex> >::iterator itv
677  = twoD_data.matching_peaks.begin();
678  for (; itv != twoD_data.matching_peaks.end(); ++itv)
679  {
680  Int i = distance(twoD_data.matching_peaks.begin(), itv);
681  for (Size j = 0; j < itv->second.size(); ++j)
682  {
683 
684 #ifdef DEBUG_2D
685  std::cout << "pos: " << itv->second[j].getPeak(ms_exp).getMZ() << "\nint: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[1][itv->second[j].peak] //itv->second[j].getPeak(ms_exp).getIntensity()
686  << "\nlw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[3][itv->second[j].peak]
687  << "\nrw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[4][itv->second[j].peak] << "\n";
688 
689 #endif
690  double mz = x_init(twoD_data.total_nr_peaks + 3 * i);
691  ms_exp[itv->second[j].spectrum][itv->second[j].peak].setMZ(mz);
692  double height = x_init(peak_idx);
693  ms_exp[itv->second[j].spectrum].getFloatDataArrays()[1][itv->second[j].peak] = height;
694  double left_width = x_init(twoD_data.total_nr_peaks + 3 * i + 1);
695  ms_exp[itv->second[j].spectrum].getFloatDataArrays()[3][itv->second[j].peak] = left_width;
696  double right_width = x_init(twoD_data.total_nr_peaks + 3 * i + 2);
697 
698  ms_exp[itv->second[j].spectrum].getFloatDataArrays()[4][itv->second[j].peak] = right_width;
699  // calculate area
700  if ((PeakShape::Type)(Int)ms_exp[itv->second[j].spectrum].getFloatDataArrays()[5][itv->second[j].peak] == PeakShape::LORENTZ_PEAK)
701  {
702  double x_left_endpoint = mz - 1 / left_width* sqrt(height / 1 - 1);
703  double x_rigth_endpoint = mz + 1 / right_width* sqrt(height / 1 - 1);
704  double area_left = -height / left_width* atan(left_width * (x_left_endpoint - mz));
705  double area_right = -height / right_width* atan(right_width * (mz - x_rigth_endpoint));
706  ms_exp[itv->second[j].spectrum][itv->second[j].peak].setIntensity(area_left + area_right);
707  }
708  else // it's a sech peak
709  {
710  double x_left_endpoint = mz - 1 / left_width* boost::math::acosh(sqrt(height / 0.001));
711  double x_rigth_endpoint = mz + 1 / right_width* boost::math::acosh(sqrt(height / 0.001));
712  double area_left = -height / left_width * (sinh(left_width * (mz - x_left_endpoint)) / cosh(left_width * (mz - x_left_endpoint)));
713  double area_right = -height / right_width * (sinh(right_width * (mz - x_rigth_endpoint)) / cosh(right_width * (mz - x_rigth_endpoint)));
714  ms_exp[itv->second[j].spectrum][itv->second[j].peak].setIntensity(area_left + area_right);
715  }
716 
717 
718 #ifdef DEBUG_2D
719  std::cout << "pos: " << itv->second[j].getPeak(ms_exp).getMZ() << "\nint: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[1][itv->second[j].peak] //itv->second[j].getPeak(ms_exp).getIntensity()
720  << "\nlw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[3][itv->second[j].peak]
721  << "\nrw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[4][itv->second[j].peak] << "\n";
722 #endif
723 
724  ++peak_idx;
725 
726 
727  }
728  }
729 
730  ++counter;
731  } // end for
732  //#undef DEBUG_2D
733  }
734 
735  template <typename InputSpectrumIterator, typename OutputPeakType>
736  void TwoDOptimization::optimizeRegionsScanwise_(InputSpectrumIterator& first,
737  InputSpectrumIterator& last,
739  {
740  Int counter = 0;
742  d.picked_peaks = ms_exp;
743  d.raw_data_first = first;
744 
745  //std::cout << "richtig hier" << std::endl;
747 
748 
749  DataValue dv = param_.getValue("penalties:position");
750  if (dv.isEmpty() || dv.toString() == "")
751  penalties.pos = 0.;
752  else
753  penalties.pos = (float)dv;
754 
755  dv = param_.getValue("penalties:left_width");
756  if (dv.isEmpty() || dv.toString() == "")
757  penalties.lWidth = 1.;
758  else
759  penalties.lWidth = (float)dv;
760 
761  dv = param_.getValue("penalties:right_width");
762  if (dv.isEmpty() || dv.toString() == "")
763  penalties.rWidth = 1.;
764  else
765  penalties.rWidth = (float)dv;
766 #ifdef DEBUG_2D
767  std::cout << penalties.pos << " "
768  << penalties.rWidth << " "
769  << penalties.lWidth << std::endl;
770 #endif
771 // MSExperiment<Peak1D >::const_iterator help = first;
772 // // std::cout << "\n\n\n\n---------------------------------------------------------------";
773 // while(help!=last)
774 // {
775 // // std::cout<<help->getRT()<<std::endl;
776 // ++help;
777 // }
778  // std::cout << "---------------------------------------------------------------\n\n\n\n";
779 
780  UInt max_iteration;
781  dv = param_.getValue("iterations");
782  if (dv.isEmpty() || dv.toString() == "")
783  max_iteration = 15;
784  else
785  max_iteration = (UInt)dv;
786 
787  std::vector<PeakShape> peak_shapes;
788 
789 
790  // go through the clusters
791  for (std::multimap<double, IsotopeCluster>::iterator it = iso_map_.begin();
792  it != iso_map_.end();
793  ++it)
794  {
795  d.iso_map_iter = it;
796 #ifdef DEBUG_2D
797  std::cerr << "element: " << counter << std::endl;
798  std::cerr << "mz: " << it->first << std::endl << "rts: ";
799  for (Size i = 0; i < it->second.scans.size(); ++i)
800  std::cerr << it->second.scans[i] << "\n";
801  std::cerr << std::endl << "peaks: ";
802  IsotopeCluster::IndexSet::const_iterator iter = it->second.peaks.begin();
803  for (; iter != it->second.peaks.end(); ++iter)
804  std::cerr << ms_exp[iter->first].getRT() << " " << (ms_exp[iter->first][iter->second]).getMZ() << std::endl;
805  //for(Size i=0;i<it->second.peaks_.size();++i) std::cout << ms_exp[it->first].getRT() << " "<<(ms_exp[it->first][it->second]).getMZ()<<std::endl;
806  std::cerr << std::endl << std::endl;
807 
808 #endif
809  // prepare for optimization:
810  // determine the matching peaks
811  // and the endpoints of each isotope pattern in the cluster
812 
813  getRegionEndpoints_(ms_exp, first, last, counter, 400, d);
814  OptimizePick::Data data;
815 
816 
817  Size idx = 0;
818  for (Size i = 0; i < d.signal2D.size() / 2; ++i)
819  {
820  data.positions.clear();
821  data.signal.clear();
822 
824  (d.raw_data_first + d.signal2D[2 * i].first)->begin() + d.signal2D[2 * i].second;
825  Int size = distance(ms_it, (d.raw_data_first + d.signal2D[2 * i].first)->begin() + d.signal2D[2 * i + 1].second);
826  data.positions.reserve(size);
827  data.signal.reserve(size);
828 
829  while (ms_it != (d.raw_data_first + d.signal2D[2 * i].first)->begin() + d.signal2D[2 * i + 1].second)
830  {
831  data.positions.push_back(ms_it->getMZ());
832  data.signal.push_back(ms_it->getIntensity());
833  ++ms_it;
834  }
835 
836 
838  pair.first = d.iso_map_iter->second.peaks.begin()->first + idx;
839 
840  IsotopeCluster::IndexSet::const_iterator set_iter = lower_bound(d.iso_map_iter->second.peaks.begin(),
841  d.iso_map_iter->second.peaks.end(),
843 
844 
845  // find the last entry with this rt-value
846  ++pair.first;
847  IsotopeCluster::IndexSet::const_iterator set_iter2 = lower_bound(d.iso_map_iter->second.peaks.begin(),
848  d.iso_map_iter->second.peaks.end(),
850 
851  while (set_iter != set_iter2)
852  {
853  const Size peak_index = set_iter->second;
854  const MSSpectrum<>& spec = ms_exp[set_iter->first];
855  PeakShape shape(spec.getFloatDataArrays()[1][peak_index], //intensity
856  spec[peak_index].getMZ(),
857  spec.getFloatDataArrays()[3][peak_index], //left width
858  spec.getFloatDataArrays()[4][peak_index], //right width
859  spec[peak_index].getIntensity(), //area is stored in peak intensity
860  PeakShape::Type(Int(spec.getFloatDataArrays()[5][peak_index]))); //shape
861  peak_shapes.push_back(shape);
862  ++set_iter;
863  }
864 #ifdef DEBUG_2D
865  std::cout << "rt "
866  << (d.raw_data_first + d.signal2D[2 * i].first)->getRT()
867  << "\n";
868 #endif
869  OptimizePick opt(penalties, max_iteration);
870 #ifdef DEBUG_2D
871  std::cout << "vorher\n";
872 
873  for (Size p = 0; p < peak_shapes.size(); ++p)
874  {
875  std::cout << peak_shapes[p].mz_position << "\t" << peak_shapes[p].height
876  << "\t" << peak_shapes[p].left_width << "\t" << peak_shapes[p].right_width << std::endl;
877  }
878 #endif
879  opt.optimize(peak_shapes, data);
880 #ifdef DEBUG_2D
881  std::cout << "nachher\n";
882  for (Size p = 0; p < peak_shapes.size(); ++p)
883  {
884  std::cout << peak_shapes[p].mz_position << "\t" << peak_shapes[p].height
885  << "\t" << peak_shapes[p].left_width << "\t" << peak_shapes[p].right_width << std::endl;
886  }
887 #endif
888  std::sort(peak_shapes.begin(), peak_shapes.end(), PeakShape::PositionLess());
889  pair.first = d.iso_map_iter->second.peaks.begin()->first + idx;
890 
891  set_iter = lower_bound(d.iso_map_iter->second.peaks.begin(),
892  d.iso_map_iter->second.peaks.end(),
894  Size p = 0;
895  while (p < peak_shapes.size())
896  {
897  MSSpectrum<>& spec = ms_exp[set_iter->first];
898  spec[set_iter->second].setMZ(peak_shapes[p].mz_position);
899  spec.getFloatDataArrays()[3][set_iter->second] = peak_shapes[p].left_width;
900  spec.getFloatDataArrays()[4][set_iter->second] = peak_shapes[p].right_width;
901  spec.getFloatDataArrays()[1][set_iter->second] = peak_shapes[p].height; // maximum intensity
902  // calculate area
903  if (peak_shapes[p].type == PeakShape::LORENTZ_PEAK)
904  {
905  PeakShape& ps = peak_shapes[p];
906  double x_left_endpoint = ps.mz_position - 1 / ps.left_width* sqrt(ps.height / 1 - 1);
907  double x_rigth_endpoint = ps.mz_position + 1 / ps.right_width* sqrt(ps.height / 1 - 1);
908  double area_left = -ps.height / ps.left_width* atan(ps.left_width * (x_left_endpoint - ps.mz_position));
909  double area_right = -ps.height / ps.right_width* atan(ps.right_width * (ps.mz_position - x_rigth_endpoint));
910  spec[set_iter->second].setIntensity(area_left + area_right); // area is stored as peak intensity
911  }
912  else //It's a Sech - Peak
913  {
914  PeakShape& ps = peak_shapes[p];
915  double x_left_endpoint = ps.mz_position - 1 / ps.left_width* boost::math::acosh(sqrt(ps.height / 0.001));
916  double x_rigth_endpoint = ps.mz_position + 1 / ps.right_width* boost::math::acosh(sqrt(ps.height / 0.001));
917  double area_left = ps.height / ps.left_width * (sinh(ps.left_width * (ps.mz_position - x_left_endpoint)) / cosh(ps.left_width * (ps.mz_position - x_left_endpoint)));
918  double area_right = -ps.height / ps.right_width * (sinh(ps.right_width * (ps.mz_position - x_rigth_endpoint)) / cosh(ps.right_width * (ps.mz_position - x_rigth_endpoint)));
919  spec[set_iter->second].setIntensity(area_left + area_right); // area is stored as peak intensity
920  }
921  ++set_iter;
922  ++p;
923  }
924  ++idx;
925  peak_shapes.clear();
926  }
927 
928  ++counter;
929  }
930  }
931 
932  template <typename InputSpectrumIterator, typename OutputPeakType>
934  InputSpectrumIterator& first,
935  InputSpectrumIterator& last,
936  Size iso_map_idx,
937  double noise_level,
939  {
940  d.signal2D.clear();
941  typedef typename InputSpectrumIterator::value_type InputExperimentType;
942  typedef typename InputExperimentType::value_type InputPeakType;
943  typedef std::multimap<double, IsotopeCluster> MapType;
944 
945  double rt, first_peak_mz, last_peak_mz;
946 
947  //MSSpectrum<InputPeakType> spec;
949  InputPeakType peak;
950 
951  MapType::iterator iso_map_iter = iso_map_.begin();
952  for (Size i = 0; i < iso_map_idx; ++i)
953  ++iso_map_iter;
954 
955 #ifdef DEBUG2D
956  std::cout << "rt begin: " << exp[iso_map_iter->second.scans[0]].getRT()
957  << "\trt end: " << exp[iso_map_iter->second.scans[iso_map_iter->second.scans.size() - 1]].getRT()
958  << " \t" << iso_map_iter->second.scans.size() << " scans"
959  << std::endl;
960 #endif
961 
962  // get left and right endpoint for all scans in the current cluster
963  for (Size i = 0; i < iso_map_iter->second.scans.size(); ++i)
964  {
966 
967  // first the right scan through binary search
968  rt = exp[iso_map_iter->second.scans[i]].getRT();
969  spec.setRT(rt);
970  InputSpectrumIterator iter = lower_bound(first, last, spec, typename MSSpectrum<InputPeakType>::RTLess());
971  // if(iter->getRT() != rt) --iter;
972  exp_it = exp.RTBegin(rt);
973 #ifdef DEBUG2D
974  std::cout << exp_it->getRT() << " vs " << iter->getRT() << std::endl;
975 #endif
976  // now the right mz
978  pair.first = iso_map_iter->second.peaks.begin()->first + i;
979  // get iterator in peaks-set that points to the first peak in the current scan
980  IsotopeCluster::IndexSet::const_iterator set_iter = lower_bound(iso_map_iter->second.peaks.begin(),
981  iso_map_iter->second.peaks.end(),
983 
984  // consider a bit more of the signal to the left
985  first_peak_mz = (exp_it->begin() + set_iter->second)->getMZ() - 1;
986 
987  // find the last entry with this rt-value
988  ++pair.first;
989  IsotopeCluster::IndexSet::const_iterator set_iter2 = lower_bound(iso_map_iter->second.peaks.begin(),
990  iso_map_iter->second.peaks.end(),
992 
993  if (i == iso_map_iter->second.scans.size() - 1)
994  {
995  set_iter2 = iso_map_iter->second.peaks.end();
996  --set_iter2;
997  }
998  else if (set_iter2 != iso_map_iter->second.peaks.begin())
999  --set_iter2;
1000 
1001  last_peak_mz = (exp_it->begin() + set_iter2->second)->getMZ() + 1;
1002 
1003  //std::cout << rt<<": first peak mz "<<first_peak_mz << "\tlast peak mz "<<last_peak_mz <<std::endl;
1004  peak.setPosition(first_peak_mz);
1006  = lower_bound(iter->begin(), iter->end(), peak, typename InputPeakType::PositionLess());
1007  if (raw_data_iter != iter->begin())
1008  {
1009  --raw_data_iter;
1010  }
1011  double intensity = raw_data_iter->getIntensity();
1012  // while the intensity is falling go to the left
1013  while (raw_data_iter != iter->begin() && (raw_data_iter - 1)->getIntensity() < intensity &&
1014  (raw_data_iter - 1)->getIntensity() > noise_level)
1015  {
1016  --raw_data_iter;
1017  intensity = raw_data_iter->getIntensity();
1018  }
1019  ++raw_data_iter;
1020  IsotopeCluster::IndexPair left, right;
1021  left.first = distance(first, iter);
1022  left.second = raw_data_iter - iter->begin();
1023 #ifdef DEBUG2D
1024  std::cout << "left: " << iter->getRT() << "\t" << raw_data_iter->getMZ() << std::endl;
1025 #endif
1026  // consider a bit more of the signal to the right
1027  peak.setPosition(last_peak_mz + 1);
1028  raw_data_iter
1029  = upper_bound(iter->begin(), iter->end(), peak, typename InputPeakType::PositionLess());
1030  if (raw_data_iter == iter->end())
1031  --raw_data_iter;
1032  intensity = raw_data_iter->getIntensity();
1033  // while the intensity is falling go to the right
1034  while (raw_data_iter + 1 != iter->end() && (raw_data_iter + 1)->getIntensity() < intensity)
1035  {
1036  ++raw_data_iter;
1037  intensity = raw_data_iter->getIntensity();
1038  if ((raw_data_iter + 1 != iter->end()) && (raw_data_iter + 1)->getIntensity() > noise_level)
1039  break;
1040  }
1041  right.first = left.first;
1042  right.second = raw_data_iter - iter->begin();
1043 #ifdef DEBUG2D
1044  std::cout << "right: " << iter->getRT() << "\t" << raw_data_iter->getMZ() << std::endl;
1045 #endif
1046  // region endpoints are stored in global vector
1047  d.signal2D.push_back(left);
1048  d.signal2D.push_back(right);
1049  }
1050 #ifdef DEBUG2D
1051  //std::cout << "fertig"<< std::endl;
1052  std::cout << first_peak_mz << "\t" << last_peak_mz << std::endl;
1053 #endif
1054  }
1055 
1056 }
1057 
1058 #endif //OPENMS_TRANSFORMATIONS_RAW2PEAK_TWODOPTIMIZATION_H
std::vector< double > positions
Positions and intensity values of the raw data.
Definition: OptimizePick.h:103
const int m_values
Definition: TwoDOptimization.h:196
const double k
double left_width
Left width parameter.
Definition: PeakShape.h:122
int inputs() const
Definition: TwoDOptimization.h:185
Definition: OptimizePick.h:100
void setMaxPeakDistance(double max_peak_distance)
Mutable access to the maximal peak distance in a cluster.
Definition: TwoDOptimization.h:118
A more convenient string class.
Definition: String.h:57
std::map< Int, std::vector< PeakIndex > > matching_peaks_
Indices of peaks in the adjacent scans matching peaks in the scan with no. ref_scan.
Definition: TwoDOptimization.h:215
Size size() const
Definition: MSExperiment.h:117
std::pair< Size, Size > IndexPair
An index e.g. in an MSExperiment.
Definition: IsotopeCluster.h:48
std::vector< SpectrumType >::iterator Iterator
Mutable iterator.
Definition: MSExperiment.h:101
const OptimizationFunctions::PenaltyFactorsIntensity & getPenalties() const
Non-mutable access to the minimal number of adjacent scans.
Definition: TwoDOptimization.h:134
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:150
Type
Peak shape type (asymmetric Lorentzian or asymmetric hyperbolic secans squared).
Definition: PeakShape.h:69
Comparison of mz_positions.
Definition: PeakShape.h:141
void setMZTolerance(double tolerance_mz)
Mutable access to the matching epsilon.
Definition: TwoDOptimization.h:109
Comparator for the retention time.
Definition: MSSpectrum.h:92
unsigned int UInt
Unsigned integer type.
Definition: Types.h:88
double getMZTolerance() const
Non-mutable access to the matching epsilon.
Definition: TwoDOptimization.h:107
Peak2D PeakType
Definition: MassTrace.h:48
Helper struct (contains the size of an area and a raw data container)
Definition: TwoDOptimization.h:169
Iterator begin()
Definition: MSExperiment.h:147
virtual ~TwoDOptimization()
Destructor.
Definition: TwoDOptimization.h:100
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
std::multimap< double, IsotopeCluster >::iterator iso_map_iter
Definition: TwoDOptimization.h:172
ContainerType::iterator Iterator
Mutable iterator.
Definition: MSSpectrum.h:121
double lWidth
Penalty factor for the peak shape's left width parameter.
Definition: OptimizePick.h:83
const TwoDOptimization::Data * m_data
Definition: TwoDOptimization.h:197
void getRegionEndpoints_(MSExperiment< OutputPeakType > &exp, InputSpectrumIterator &first, InputSpectrumIterator &last, Size iso_map_idx, double noise_level, TwoDOptimization::Data &d)
Get the indices of the first and last raw data point of this region.
Definition: TwoDOptimization.h:933
Class to hold strings, numeric values, lists of strings and lists of numeric values.
Definition: DataValue.h:57
void optimizeRegions_(InputSpectrumIterator &first, InputSpectrumIterator &last, MSExperiment< OutputPeakType > &ms_exp)
Definition: TwoDOptimization.h:570
OptimizationFunctions::PenaltyFactorsIntensity penalties
Definition: TwoDOptimization.h:177
Iterator end()
Definition: MSExperiment.h:157
MSExperiment< Peak1D > MapType
Definition: PeakPickerIterative.cpp:87
ConstIterator RTBegin(CoordinateType rt) const
Fast search for spectrum range begin.
Definition: MSExperiment.h:374
Base::iterator iterator
Definition: MSExperiment.h:114
void findMatchingPeaks_(std::multimap< double, IsotopeCluster >::iterator &it, MSExperiment<> &ms_exp)
Identify matching peak in a peak cluster.
Internal representation of a peak shape (used by the PeakPickerCWT)
Definition: PeakShape.h:51
void setPenalties(OptimizationFunctions::PenaltyFactorsIntensity &penalties)
Mutable access to the minimal number of adjacent scans.
Definition: TwoDOptimization.h:136
A method or algorithm argument contains illegal values.
Definition: Exception.h:634
ChargedIndexSet peaks
peaks in this cluster
Definition: IsotopeCluster.h:72
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
double max_peak_distance_
upper bound for distance between two peaks belonging to the same region
Definition: TwoDOptimization.h:208
TwoDOptFunctor(unsigned dimensions, unsigned num_data_points, const TwoDOptimization::Data *data)
Definition: TwoDOptimization.h:188
std::vector< double >::iterator searchInScan_(std::vector< double >::iterator scan_begin, std::vector< double >::iterator scan_end, double current_mz)
void optimize(std::vector< PeakShape > &peaks, Data &data)
Start the optimization of the peak shapes peaks. The original peak shapes will be substituted by the ...
bool real_2D_
Optimization considering all scans of a cluster or optimization of each scan separately.
Definition: TwoDOptimization.h:222
double height
Maximum intensity of the peak shape.
Definition: PeakShape.h:118
double tolerance_mz_
threshold for the difference in the peak position of two matching peaks
Definition: TwoDOptimization.h:211
Size total_nr_peaks
Definition: TwoDOptimization.h:173
std::multimap< double, IsotopeCluster >::const_iterator curr_region_
Pointer to the current region.
Definition: TwoDOptimization.h:205
OptimizationFunctions::PenaltyFactorsIntensity penalties_
Penalty factors for some parameters in the optimization.
Definition: TwoDOptimization.h:226
std::vector< std::pair< SignedSize, SignedSize > > signal2D
Definition: TwoDOptimization.h:171
void setRT(double rt)
Sets the absolute retention time (is seconds)
Definition: MSSpectrum.h:249
MSExperiment< Peak1D >::ConstIterator raw_data_first
Definition: TwoDOptimization.h:176
double pos
Penalty factor for the peak shape's position.
Definition: OptimizePick.h:81
In-Memory representation of a mass spectrometry experiment.
Definition: MSExperiment.h:69
MSExperiment picked_peaks
Definition: TwoDOptimization.h:175
Exception used if an error occurred while fitting a model to a given dataset.
Definition: Exception.h:662
std::vector< double > positions
Definition: TwoDOptimization.h:178
double rWidth
Penalty factor for the peak shape's right width parameter.
Definition: OptimizePick.h:85
std::vector< Size > scans
the scans of this cluster
Definition: IsotopeCluster.h:75
double getMaxPeakDistance() const
Non-mutable access to the maximal peak distance in a cluster.
Definition: TwoDOptimization.h:116
std::vector< double > signal
Definition: TwoDOptimization.h:179
bool empty() const
Definition: MSExperiment.h:127
void setMaxIterations(UInt max_iteration)
Mutable access to the maximal number of iterations.
Definition: TwoDOptimization.h:127
double right_width
Right width parameter.
Definition: PeakShape.h:124
Stores information about an isotopic cluster (i.e. potential peptide charge variants) ...
Definition: IsotopeCluster.h:45
Class for the penalty factors used during the optimization.
Definition: OptimizePick.h:63
Int charge
charge estimate (convention: zero means "no charge estimate")
Definition: IsotopeCluster.h:62
std::multimap< double, IsotopeCluster > iso_map_
stores the retention time of each isotopic cluster
Definition: TwoDOptimization.h:202
UInt max_iteration_
Convergence Parameter: Maximal number of iterations.
Definition: TwoDOptimization.h:219
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:92
std::vector< double > signal
Definition: OptimizePick.h:104
UInt getMaxIterations() const
Non-mutable access to the maximal number of iterations.
Definition: TwoDOptimization.h:125
std::map< Int, std::vector< PeakIndex > > matching_peaks
Definition: TwoDOptimization.h:174
Class for comparison of std::pair using first ONLY e.g. for use with std::sort.
Definition: ComparatorUtils.h:326
const FloatDataArrays & getFloatDataArrays() const
Returns a const reference to the float meta data arrays.
Definition: MSSpectrum.h:298
void optimizeRegionsScanwise_(InputSpectrumIterator &first, InputSpectrumIterator &last, MSExperiment< OutputPeakType > &ms_exp)
Definition: TwoDOptimization.h:736
Definition: PeakShape.h:71
int Int
Signed integer type.
Definition: Types.h:96
void optimize(InputSpectrumIterator first, InputSpectrumIterator last, MSExperiment< OutputPeakType > &ms_exp, bool real2D=true)
Find two dimensional peak clusters and optimize their peak parameters.
Definition: TwoDOptimization.h:271
double mz_position
Centroid position.
Definition: PeakShape.h:120
double height
Definition: OptimizePeakDeconvolution.h:77
Definition: TwoDOptimization.h:182
This class provides the non-linear optimization of the peak parameters.
Definition: OptimizePick.h:96
This class provides the two-dimensional optimization of the picked peak parameters.
Definition: TwoDOptimization.h:88
int values() const
Definition: TwoDOptimization.h:186
Class for the penalty factors used during the optimization.
Definition: OptimizePeakDeconvolution.h:58

OpenMS / TOPP release 2.0.0 Documentation generated on Thu Aug 20 2015 01:44:30 using doxygen 1.8.9.1