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-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: 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 
59 
60 #ifndef OPENMS_SYSTEM_STOPWATCH_H
61 #endif
62 
63 #include <boost/math/special_functions/acosh.hpp>
64 #include <gsl/gsl_vector.h>
65 #include <gsl/gsl_multifit_nlin.h>
66 #include <gsl/gsl_blas.h>
69 
70 namespace OpenMS
71 {
86  class OPENMS_DLLAPI TwoDOptimization :
87  public DefaultParamHandler
88  {
89 public:
90 
93 
96 
98  virtual ~TwoDOptimization(){}
99 
101  TwoDOptimization& operator=(const TwoDOptimization& opt);
102 
103 
105  inline DoubleReal getMZTolerance() const {return tolerance_mz_; }
107  inline void setMZTolerance(DoubleReal tolerance_mz)
108  {
109  tolerance_mz_ = tolerance_mz;
110  param_.setValue("2d:tolerance_mz", tolerance_mz);
111  }
112 
114  inline DoubleReal getMaxPeakDistance() const {return max_peak_distance_; }
116  inline void setMaxPeakDistance(DoubleReal max_peak_distance)
117  {
118  max_peak_distance_ = max_peak_distance;
119  param_.setValue("2d:max_peak_distance", max_peak_distance);
120  }
121 
123  inline DoubleReal getMaxAbsError() const {return eps_abs_; }
125  inline void setMaxAbsError(DoubleReal eps_abs)
126  {
127  eps_abs_ = eps_abs;
128  param_.setValue("delta_abs_error", eps_abs);
129  }
130 
132  inline DoubleReal getMaxRelError() const {return eps_rel_; }
134  inline void setMaxRelError(DoubleReal eps_rel)
135  {
136  eps_rel_ = eps_rel;
137  param_.setValue("delta_rel_error", eps_rel);
138  }
139 
141  inline UInt getMaxIterations() const {return max_iteration_; }
143  inline void setMaxIterations(UInt max_iteration)
144  {
145  max_iteration_ = max_iteration;
146  param_.setValue("iterations", max_iteration);
147  }
148 
150  inline const OptimizationFunctions::PenaltyFactorsIntensity& getPenalties() const {return penalties_; }
153  {
154  penalties_ = penalties;
155  param_.setValue("penalties:position", penalties.pos);
156  param_.setValue("penalties:height", penalties.height);
157  param_.setValue("penalties:left_width", penalties.lWidth);
158  param_.setValue("penalties:right_width", penalties.rWidth);
159  }
160 
177  template <typename InputSpectrumIterator, typename OutputPeakType>
178  void optimize(InputSpectrumIterator first,
179  InputSpectrumIterator last,
180  MSExperiment<OutputPeakType>& ms_exp, bool real2D = true);
181 
182 
183 protected:
185  struct Data
186  {
187  std::vector<std::pair<SignedSize, SignedSize> > signal2D;
188  std::multimap<DoubleReal, IsotopeCluster>::iterator iso_map_iter;
190  std::map<Int, std::vector<PeakIndex> > matching_peaks;
194  std::vector<DoubleReal> positions;
195  std::vector<DoubleReal> signal;
196  };
197 
199  std::multimap<DoubleReal, IsotopeCluster> iso_map_;
200 
202  std::multimap<DoubleReal, IsotopeCluster>::const_iterator curr_region_;
203 
206 
209 
211  // std::map<Int, std::vector<MSExperiment<>::SpectrumType::Iterator > > matching_peaks_;
212  std::map<Int, std::vector<PeakIndex> > matching_peaks_;
213 
214 
217 
220 
223 
225  bool real_2D_;
226 
227 
230 
231 
236  static Int residual2D_(const gsl_vector* x, void* params, gsl_vector* f);
239  static Int jacobian2D_(const gsl_vector* x, void* params, gsl_matrix* J);
241  static Int evaluate2D_(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* J);
242 
243 
248  std::vector<DoubleReal>::iterator searchInScan_(std::vector<DoubleReal>::iterator scan_begin,
249  std::vector<DoubleReal>::iterator scan_end,
250  DoubleReal current_mz);
251 
253  template <typename InputSpectrumIterator, typename OutputPeakType>
254  void optimizeRegions_(InputSpectrumIterator& first,
255  InputSpectrumIterator& last,
257 
259  template <typename InputSpectrumIterator, typename OutputPeakType>
260  void optimizeRegionsScanwise_(InputSpectrumIterator& first,
261  InputSpectrumIterator& last,
263 
264 
266  template <typename InputSpectrumIterator, typename OutputPeakType>
267  void getRegionEndpoints_(MSExperiment<OutputPeakType>& exp,
268  InputSpectrumIterator& first,
269  InputSpectrumIterator& last,
270  Size iso_map_idx,
271  DoubleReal noise_level,
273 
275  void findMatchingPeaks_(std::multimap<DoubleReal, IsotopeCluster>::iterator& it,
276  MSExperiment<>& ms_exp);
277 
279 
281  void updateMembers_();
282  };
283 
284 
285  template <typename InputSpectrumIterator, typename OutputPeakType>
286  void TwoDOptimization::optimize(InputSpectrumIterator first, InputSpectrumIterator last, MSExperiment<OutputPeakType>& ms_exp, bool real2D)
287  {
288  //#define DEBUG_2D
289  //check if the input maps have the same number of spectra
290  if ((UInt)distance(first, last) != ms_exp.size())
291  {
292  throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__, "Error in Two2Optimization: Raw and peak map do not have the same number of spectra");
293  }
294  //do nothing if there are no scans
295  if (ms_exp.empty())
296  {
297  return;
298  }
299  //check if required meta data arrays are present (for each scan)
300  for (Size i = 0; i < ms_exp.size(); ++i)
301  {
302  //check if enough meta data arrays are present
303  if (ms_exp[i].getFloatDataArrays().size() < 6)
304  {
305  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)");
306  }
307  bool area = ms_exp[i].getFloatDataArrays()[1].getName() == "maximumIntensity";
308  bool wleft = ms_exp[i].getFloatDataArrays()[3].getName() == "leftWidth";
309  bool wright = ms_exp[i].getFloatDataArrays()[4].getName() == "rightWidth";
310  bool shape = ms_exp[i].getFloatDataArrays()[5].getName() == "peakShape";
311 
312  if (!area || !wleft || !wright || !shape)
313  {
314  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)");
315  }
316  }
317  real_2D_ = real2D;
318  typedef typename InputSpectrumIterator::value_type InputSpectrumType;
319  typedef typename InputSpectrumType::value_type PeakType;
320  typedef MSSpectrum<PeakType> SpectrumType;
321 
322  typename MSExperiment<OutputPeakType>::Iterator ms_exp_it = ms_exp.begin();
323  typename MSExperiment<OutputPeakType>::Iterator ms_exp_it_end = ms_exp.end();
324  if (ms_exp.empty())
325  {
326  std::cout << "empty experiment" << std::endl;
327  return;
328  }
329  // stores the monoisotopic peaks of isotopic clusters
330  std::vector<DoubleReal> iso_last_scan;
331  std::vector<DoubleReal> iso_curr_scan;
332  std::vector<std::multimap<DoubleReal, IsotopeCluster>::iterator> clusters_last_scan;
333  std::vector<std::multimap<DoubleReal, IsotopeCluster>::iterator> clusters_curr_scan;
334  std::multimap<DoubleReal, IsotopeCluster>::iterator cluster_iter;
335  DoubleReal current_rt = ms_exp_it->getRT(), last_rt = 0;
336 
337  // retrieve values for accepted peaks distances
338  max_peak_distance_ = param_.getValue("2d:max_peak_distance");
339  DoubleReal tolerance_mz = param_.getValue("2d:tolerance_mz");
340 
341  UInt current_charge = 0; // charge state of the current isotopic cluster
342  DoubleReal mz_in_hash = 0; // used as reference to the current isotopic peak
343 
344  // sweep through scans
345  for (UInt curr_scan = 0; ms_exp_it + curr_scan != ms_exp_it_end; ++curr_scan)
346  {
347  Size nr_peaks_in_scan = (ms_exp_it + curr_scan)->size();
348  if (nr_peaks_in_scan == 0)
349  continue;
350 
351  //last_rt = current_rt;
352  current_rt = (ms_exp_it + curr_scan)->getRT();
353  typename MSExperiment<OutputPeakType>::SpectrumType::Iterator peak_it = (ms_exp_it + curr_scan)->begin();
354 
355  // copy cluster information of least scan
356  iso_last_scan = iso_curr_scan;
357  iso_curr_scan.clear();
358  clusters_last_scan = clusters_curr_scan;
359  clusters_curr_scan.clear();
360 
361 #ifdef DEBUG_2D
362  std::cout << "Next scan with rt: " << current_rt << std::endl;
363  std::cout << "Next scan, rt = " << current_rt << " last_rt: " << last_rt << std::endl;
364  std::cout << "---------------------------------------------------------------------------" << std::endl;
365 #endif
367  s.setRT(current_rt);
368  // check if there were scans in between
369  if (last_rt == 0 || // are we in the first scan
370  ((lower_bound(first, last, s, typename SpectrumType::RTLess()) - 1)->getRT() == last_rt))
371  {
372 
373 
374  for (UInt curr_peak = 0; curr_peak < (ms_exp_it + curr_scan)->size() - 1; ++curr_peak)
375  {
376 
377  // store the m/z of the current peak
378  DoubleReal curr_mz = (peak_it + curr_peak)->getMZ();
379  DoubleReal dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - curr_mz;
380 
381  if (dist2nextpeak <= max_peak_distance_) // one single peak without neighbors isn't optimized
382  {
383 #ifdef DEBUG_2D
384  std::cout << "Isotopic pattern found ! " << std::endl;
385  std::cout << "We are at: " << (peak_it + curr_peak)->getMZ() << " " << curr_mz << std::endl;
386 #endif
387  if (!iso_last_scan.empty()) // Did we find any isotopic cluster in the last scan?
388  {
389  std::sort(iso_last_scan.begin(), iso_last_scan.end());
390  // there were some isotopic clustures in the last scan...
391  std::vector<DoubleReal>::iterator it =
392  searchInScan_(iso_last_scan.begin(), iso_last_scan.end(), curr_mz);
393 
394  DoubleReal delta_mz = fabs(*it - curr_mz);
395  //std::cout << delta_mz << " "<< tolerance_mz << std::endl;
396  if (delta_mz > tolerance_mz) // check if first peak of last cluster is close enough
397  {
398  mz_in_hash = curr_mz; // update current hash key
399 
400  // create new isotopic cluster
401 // #ifdef DEBUG_2D
402 // std::cout << "Last peak cluster too far, creating new cluster at "<<curr_mz << std::endl;
403 // #endif
404  IsotopeCluster new_cluster;
405  new_cluster.peaks.charge = current_charge;
406  new_cluster.scans.push_back(curr_scan);
407  cluster_iter = iso_map_.insert(std::pair<DoubleReal, IsotopeCluster>(mz_in_hash, new_cluster));
408 
409  }
410  else
411  {
412 // //#ifdef DEBUG_2D
413 // std::cout << "Found neighbouring peak with distance (m/z) " << delta_mz << std::endl;
414 // //#endif
415  cluster_iter = clusters_last_scan[distance(iso_last_scan.begin(), it)];
416 
417  // check whether this scan is already contained
418  if (find(cluster_iter->second.scans.begin(), cluster_iter->second.scans.end(), curr_scan)
419  == cluster_iter->second.scans.end())
420  {
421  cluster_iter->second.scans.push_back(curr_scan);
422  }
423 
424 // //#ifdef DEBUG_2D
425 // std::cout << "Cluster with " << cluster_iter->second.peaks.size()
426 // << " peaks retrieved." << std::endl;
427 // //#endif
428  }
429 
430  }
431  else // last scan did not contain any isotopic cluster
432  {
433 // //#ifdef DEBUG_2D
434 // std::cout << "Last scan was empty => creating new cluster." << std::endl;
435 // std::cout << "Creating new cluster at m/z: " << curr_mz << std::endl;
436 // //#endif
437 
438  mz_in_hash = curr_mz; // update current hash key
439 
440  // create new isotopic cluster
441  IsotopeCluster new_cluster;
442  new_cluster.peaks.charge = current_charge;
443  new_cluster.scans.push_back(curr_scan);
444  cluster_iter = iso_map_.insert(std::pair<DoubleReal, IsotopeCluster>(mz_in_hash, new_cluster));
445 
446  }
447 
448 // //#ifdef DEBUG_2D
449 // std::cout << "Storing found peak in current isotopic cluster" << std::endl;
450 // //#endif
451 
452 
453 
454  cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
455 
456  iso_curr_scan.push_back(mz_in_hash);
457  clusters_curr_scan.push_back(cluster_iter);
458  ++curr_peak;
459 
460  cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
461  iso_curr_scan.push_back((peak_it + curr_peak)->getMZ());
462  clusters_curr_scan.push_back(cluster_iter);
463 
464  // check distance to next peak
465  if ((curr_peak + 1) >= nr_peaks_in_scan)
466  break;
467  dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - (peak_it + curr_peak)->getMZ();
468 
469 
470  // loop until end of isotopic pattern in this scan
471  while (dist2nextpeak <= max_peak_distance_
472  && curr_peak < (nr_peaks_in_scan - 1))
473  {
474  cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak + 1)); // save peak in cluster
475  iso_curr_scan.push_back((peak_it + curr_peak + 1)->getMZ());
476  clusters_curr_scan.push_back(cluster_iter);
477  // std::cout << "new enter'd: "<<(peak_it+curr_peak+1)->getMZ()<<" im while"<<std::endl;
478  ++curr_peak;
479  if (curr_peak >= nr_peaks_in_scan - 1)
480  break;
481  dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - (peak_it + curr_peak)->getMZ(); // get distance to next peak
482 
483 
484  } // end while(...)
485 
486 
487 
488  } // end of if (dist2nextpeak <= max_peak_distance_)
489  else
490  {
491  if (!iso_last_scan.empty()) // Did we find any isotopic cluster in the last scan?
492  {
493  std::sort(iso_last_scan.begin(), iso_last_scan.end());
494  // there were some isotopic clusters in the last scan...
495  std::vector<DoubleReal>::iterator it =
496  searchInScan_(iso_last_scan.begin(), iso_last_scan.end(), curr_mz);
497 
498  DoubleReal delta_mz = fabs(*it - curr_mz);
499  // std::cout << delta_mz << " "<< tolerance_mz << std::endl;
500  if (delta_mz > tolerance_mz) // check if first peak of last cluster is close enough
501  {
502  mz_in_hash = curr_mz; // update current hash key
503 
504  // create new isotopic cluster
505 // //#ifdef DEBUG_2D
506 // std::cout << "Last peak cluster too far, creating new cluster at "<<curr_mz << std::endl;
507 // //#endif
508  IsotopeCluster new_cluster;
509  new_cluster.peaks.charge = current_charge;
510  new_cluster.scans.push_back(curr_scan);
511  cluster_iter = iso_map_.insert(std::pair<DoubleReal, IsotopeCluster>(mz_in_hash, new_cluster));
512 
513  }
514  else
515  {
516 // //#ifdef DEBUG_2D
517 // std::cout << "Found neighbouring peak with distance (m/z) " << delta_mz << std::endl;
518 // //#endif
519  cluster_iter = clusters_last_scan[distance(iso_last_scan.begin(), it)];
520 
521  // check whether this scan is already contained
522  if (find(cluster_iter->second.scans.begin(), cluster_iter->second.scans.end(), curr_scan)
523  == cluster_iter->second.scans.end())
524  {
525  cluster_iter->second.scans.push_back(curr_scan);
526  }
527 
528 // //#ifdef DEBUG_2D
529 // std::cout << "Cluster with " << cluster_iter->second.peaks.size()
530 // << " peaks retrieved." << std::endl;
531 // //#endif
532  }
533 
534  }
535  else // last scan did not contain any isotopic cluster
536  {
537 // //#ifdef DEBUG_2D
538 // std::cout << "Last scan was empty => creating new cluster." << std::endl;
539 // std::cout << "Creating new cluster at m/z: " << curr_mz << std::endl;
540 // //#endif
541 
542  mz_in_hash = curr_mz; // update current hash key
543 
544  // create new isotopic cluster
545  IsotopeCluster new_cluster;
546  new_cluster.peaks.charge = current_charge;
547  new_cluster.scans.push_back(curr_scan);
548  cluster_iter = iso_map_.insert(std::pair<DoubleReal, IsotopeCluster>(mz_in_hash, new_cluster));
549 
550  }
551 
552 // //#ifdef DEBUG_2D
553 // std::cout << "Storing found peak in current isotopic cluster" << std::endl;
554 // //#endif
555 
556 
557 
558  cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
559 
560  iso_curr_scan.push_back(mz_in_hash);
561  clusters_curr_scan.push_back(cluster_iter);
562 
563 
564  }
565 
566  current_charge = 0; // reset charge
567  } // end for (...)
568  }
569  last_rt = current_rt;
570  }
571  curr_region_ = iso_map_.begin();
572 #ifdef DEBUG_2D
573  std::cout << iso_map_.size() << " isotopic clusters were found ! " << std::endl;
574 #endif
575 
576  if (real_2D_)
577  optimizeRegions_(first, last, ms_exp);
578  else
579  optimizeRegionsScanwise_(first, last, ms_exp);
580  //#undef DEBUG_2D
581  }
582 
583  template <typename InputSpectrumIterator, typename OutputPeakType>
584  void TwoDOptimization::optimizeRegions_(InputSpectrumIterator& first,
585  InputSpectrumIterator& last,
587  {
588  Int counter = 0;
589  // go through the clusters
590  for (std::multimap<DoubleReal, IsotopeCluster>::iterator it = iso_map_.begin();
591  it != iso_map_.end();
592  ++it)
593  {
594 #ifdef DEBUG_2D
595  std::cout << "element: " << counter << std::endl;
596  std::cout << "mz: " << it->first << std::endl << "rts: ";
597 // for(Size i=0;i<it->second.scans.size();++i) std::cout << it->second.scans[i] << "\n";
598  std::cout << std::endl << "peaks: ";
599  IsotopeCluster::IndexSet::const_iterator iter = it->second.peaks.begin();
600  for (; iter != it->second.peaks.end(); ++iter)
601  std::cout << ms_exp[iter->first].getRT() << " " << (ms_exp[iter->first][iter->second]).getMZ() << std::endl;
602 
603 //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;
604  std::cout << std::endl << std::endl;
605 
606 #endif
607 
608  // prepare for optimization:
609  // determine the matching peaks
610  matching_peaks_.clear();
611  findMatchingPeaks_(it, ms_exp);
613  d.penalties = penalties_;
615  // and the endpoints of each isotope pattern in the cluster
616  getRegionEndpoints_(ms_exp, first, last, counter, 400, d);
617 
618  // peaks have to be stored globally
619  d.iso_map_iter = it;
620 
621  d.picked_peaks = ms_exp;
622  d.raw_data_first = first;
623 
624  Size nr_diff_peaks = matching_peaks_.size();
625  d.total_nr_peaks = it->second.peaks.size();
626 
627  Size nr_parameters = nr_diff_peaks * 3 + d.total_nr_peaks;
628 
629  gsl_vector* start_value = gsl_vector_alloc(nr_parameters);
630  gsl_vector_set_zero(start_value);
631 
632  // initialize parameters for optimization
633  std::map<Int, std::vector<PeakIndex> >::iterator m_peaks_it = d.matching_peaks.begin();
634  DoubleReal av_mz = 0, av_lw = 0, av_rw = 0, avr_height = 0, height;
635  Int peak_counter = 0;
636  Int diff_peak_counter = 0;
637  // go through the matching peaks
638  for (; m_peaks_it != d.matching_peaks.end(); ++m_peaks_it)
639  {
640  av_mz = 0, av_lw = 0, av_rw = 0, avr_height = 0;
641  std::vector<PeakIndex>::iterator iter_iter = (m_peaks_it)->second.begin();
642  for (; iter_iter != m_peaks_it->second.end(); ++iter_iter)
643  {
644  height = ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[1][(iter_iter)->peak]; //(iter_iter)->getPeak(ms_exp).getIntensity();
645  avr_height += height;
646  av_mz += (iter_iter)->getPeak(ms_exp).getMZ() * height;
647  av_lw += ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[3][(iter_iter)->peak] * height; //left width
648  av_rw += ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[4][(iter_iter)->peak] * height; //right width
649  gsl_vector_set(start_value, peak_counter, height);
650  ++peak_counter;
651  }
652  gsl_vector_set(start_value, d.total_nr_peaks + 3 * diff_peak_counter, av_mz / avr_height);
653  gsl_vector_set(start_value, d.total_nr_peaks + 3 * diff_peak_counter + 1, av_lw / avr_height);
654  gsl_vector_set(start_value, d.total_nr_peaks + 3 * diff_peak_counter + 2, av_rw / avr_height);
655  ++diff_peak_counter;
656  }
657 
658 #ifdef DEBUG_2D
659  std::cout << "----------------------------\n\nstart_value: " << std::endl;
660  for (Size k = 0; k < start_value->size; ++k)
661  {
662  std::cout << gsl_vector_get(start_value, k) << std::endl;
663  }
664 #endif
665  Int num_positions = 0;
666  for (Size i = 0; i < d.signal2D.size(); i += 2)
667  {
668  num_positions += (d.signal2D[i + 1].second - d.signal2D[i].second + 1);
669 #ifdef DEBUG_2D
670  std::cout << d.signal2D[i + 1].second << " - " << d.signal2D[i].second << " +1 " << std::endl;
671 #endif
672 
673  }
674 #ifdef DEBUG_2D
675  std::cout << "num_positions : " << num_positions << std::endl;
676 #endif
677  // The gsl algorithms require us to provide function pointers for the evaluation of
678  // the target function.
679  gsl_multifit_function_fdf fit_function;
680  fit_function.f = (Int (*)(const gsl_vector* x, void* params, gsl_vector* f)) & OpenMS::TwoDOptimization::residual2D_;
681  fit_function.df = (Int (*)(const gsl_vector* x, void* params, gsl_matrix* J)) & OpenMS::TwoDOptimization::jacobian2D_;
682  fit_function.fdf = (Int (*)(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* J)) & OpenMS::TwoDOptimization::evaluate2D_;
683  // gsl crashes when n is smaller than p!
684  fit_function.n = std::max(num_positions + 1, (Int)(nr_parameters));
685  fit_function.p = nr_parameters;
686  fit_function.params = &d;
687 #ifdef DEBUG_2D
688  std::cout << "fit_function.n " << fit_function.n
689  << "\tfit_function.p " << fit_function.p << std::endl;
690 #endif
691  const gsl_multifit_fdfsolver_type* type = gsl_multifit_fdfsolver_lmsder;
692 
693  gsl_multifit_fdfsolver* fit = gsl_multifit_fdfsolver_alloc(type,
694  std::max(num_positions + 1, (Int)(nr_parameters)),
695  nr_parameters);
696 
697  gsl_multifit_fdfsolver_set(fit, &fit_function, start_value);
698 
699 
700 
701  // initial norm
702 #ifdef DEBUG_2D
703  std::cout << "Before optimization: ||f|| = " << gsl_blas_dnrm2(fit->f) << std::endl;
704 #endif
705  // Iteration
706  UInt iteration = 0;
707  Int status;
708 
709  do
710  {
711  iteration++;
712  status = gsl_multifit_fdfsolver_iterate(fit);
713 #ifdef DEBUG_2D
714  std::cout << "Iteration " << iteration << "; Status " << gsl_strerror(status) << "; " << std::endl;
715  std::cout << "||f|| = " << gsl_blas_dnrm2(fit->f) << std::endl;
716  std::cout << "Number of parms: " << nr_parameters << std::endl;
717  std::cout << "Delta: " << gsl_blas_dnrm2(fit->dx) << std::endl;
718 #endif
719 
720  status = gsl_multifit_test_delta(fit->dx, fit->x, eps_abs_, eps_rel_);
721  if (status != GSL_CONTINUE)
722  break;
723 
724  }
725  while (status == GSL_CONTINUE && iteration < max_iteration_);
726 
727 #ifdef DEBUG_2D
728  std::cout << "Finished! No. of iterations" << iteration << std::endl;
729  std::cout << "Delta: " << gsl_blas_dnrm2(fit->dx) << std::endl;
730  DoubleReal chi = gsl_blas_dnrm2(fit->f);
731  std::cout << "After optimization: || f || = " << gsl_blas_dnrm2(fit->f) << std::endl;
732  std::cout << "chisq/dof = " << pow(chi, 2.0) / (num_positions - nr_parameters);
733 
734 
735  std::cout << "----------------------------------------------\n\nnachher" << std::endl;
736  for (Size k = 0; k < fit->x->size; ++k)
737  {
738  std::cout << gsl_vector_get(fit->x, k) << std::endl;
739  }
740 #endif
741  Int peak_idx = 0;
742  std::map<Int, std::vector<PeakIndex> >::iterator itv
743  = d.matching_peaks.begin();
744  for (; itv != d.matching_peaks.end(); ++itv)
745  {
746  Int i = distance(d.matching_peaks.begin(), itv);
747  for (Size j = 0; j < itv->second.size(); ++j)
748  {
749 
750 #ifdef DEBUG_2D
751  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()
752  << "\nlw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[3][itv->second[j].peak]
753  << "\nrw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[4][itv->second[j].peak] << "\n";
754 
755 #endif
756  DoubleReal mz = gsl_vector_get(fit->x, d.total_nr_peaks + 3 * i);
757  ms_exp[itv->second[j].spectrum][itv->second[j].peak].setMZ(mz);
758  DoubleReal height = (gsl_vector_get(fit->x, peak_idx));
759  ms_exp[itv->second[j].spectrum].getFloatDataArrays()[1][itv->second[j].peak] = height;
760  DoubleReal left_width = gsl_vector_get(fit->x, d.total_nr_peaks + 3 * i + 1);
761  ms_exp[itv->second[j].spectrum].getFloatDataArrays()[3][itv->second[j].peak] = left_width;
762  DoubleReal right_width = gsl_vector_get(fit->x, d.total_nr_peaks + 3 * i + 2);
763  ms_exp[itv->second[j].spectrum].getFloatDataArrays()[4][itv->second[j].peak] = right_width;
764  // calculate area
765  if ((PeakShape::Type)(Int)ms_exp[itv->second[j].spectrum].getFloatDataArrays()[5][itv->second[j].peak] == PeakShape::LORENTZ_PEAK)
766  {
767  DoubleReal x_left_endpoint = mz - 1 / left_width* sqrt(height / 1 - 1);
768  DoubleReal x_rigth_endpoint = mz + 1 / right_width* sqrt(height / 1 - 1);
769  DoubleReal area_left = -height / left_width* atan(left_width * (x_left_endpoint - mz));
770  DoubleReal area_right = -height / right_width* atan(right_width * (mz - x_rigth_endpoint));
771  ms_exp[itv->second[j].spectrum][itv->second[j].peak].setIntensity(area_left + area_right);
772  }
773  else // it's a sech peak
774  {
775  DoubleReal x_left_endpoint = mz - 1 / left_width* boost::math::acosh(sqrt(height / 0.001));
776  DoubleReal x_rigth_endpoint = mz + 1 / right_width* boost::math::acosh(sqrt(height / 0.001));
777  DoubleReal area_left = -height / left_width * (sinh(left_width * (mz - x_left_endpoint)) / cosh(left_width * (mz - x_left_endpoint)));
778  DoubleReal area_right = -height / right_width * (sinh(right_width * (mz - x_rigth_endpoint)) / cosh(right_width * (mz - x_rigth_endpoint)));
779  ms_exp[itv->second[j].spectrum][itv->second[j].peak].setIntensity(area_left + area_right);
780  }
781 
782 
783 #ifdef DEBUG_2D
784  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()
785  << "\nlw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[3][itv->second[j].peak]
786  << "\nrw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[4][itv->second[j].peak] << "\n";
787 
788 // std::cout << "pos: "<<itv->second[j]->getMZ()<<"\nint: "<<itv->second[j]->getIntensity()
789 // <<"\nlw: "<<itv->second[j]->getLeftWidthParameter()
790 // <<"\nrw: "<<itv->second[j]->getRightWidthParameter() << "\n";
791 
792 #endif
793 
794  ++peak_idx;
795 
796 
797  }
798  }
799 
800  gsl_multifit_fdfsolver_free(fit);
801  gsl_vector_free(start_value);
802  ++counter;
803  } // end for
804  //#undef DEBUG_2D
805  }
806 
807  template <typename InputSpectrumIterator, typename OutputPeakType>
808  void TwoDOptimization::optimizeRegionsScanwise_(InputSpectrumIterator& first,
809  InputSpectrumIterator& last,
811  {
812  Int counter = 0;
814  d.picked_peaks = ms_exp;
815  d.raw_data_first = first;
816 
817  //std::cout << "richtig hier" << std::endl;
819 
820 
821  DataValue dv = param_.getValue("penalties:position");
822  if (dv.isEmpty() || dv.toString() == "")
823  penalties.pos = 0.;
824  else
825  penalties.pos = (float)dv;
826 
827  dv = param_.getValue("penalties:left_width");
828  if (dv.isEmpty() || dv.toString() == "")
829  penalties.lWidth = 1.;
830  else
831  penalties.lWidth = (float)dv;
832 
833  dv = param_.getValue("penalties:right_width");
834  if (dv.isEmpty() || dv.toString() == "")
835  penalties.rWidth = 1.;
836  else
837  penalties.rWidth = (float)dv;
838 #ifdef DEBUG_2D
839  std::cout << penalties.pos << " "
840  << penalties.rWidth << " "
841  << penalties.lWidth << std::endl;
842 #endif
843 // MSExperiment<Peak1D >::const_iterator help = first;
844 // // std::cout << "\n\n\n\n---------------------------------------------------------------";
845 // while(help!=last)
846 // {
847 // // std::cout<<help->getRT()<<std::endl;
848 // ++help;
849 // }
850  // std::cout << "---------------------------------------------------------------\n\n\n\n";
851 
852  UInt max_iteration;
853  dv = param_.getValue("iterations");
854  if (dv.isEmpty() || dv.toString() == "")
855  max_iteration = 15;
856  else
857  max_iteration = (UInt)dv;
858 
859  DoubleReal eps_abs;
860  dv = param_.getValue("delta_abs_error");
861  if (dv.isEmpty() || dv.toString() == "")
862  eps_abs = 1e-04f;
863  else
864  eps_abs = (DoubleReal)dv;
865 
866  DoubleReal eps_rel;
867  dv = param_.getValue("delta_rel_error");
868  if (dv.isEmpty() || dv.toString() == "")
869  eps_rel = 1e-04f;
870  else
871  eps_rel = (DoubleReal)dv;
872 
873  std::vector<PeakShape> peak_shapes;
874 
875 
876  // go through the clusters
877  for (std::multimap<DoubleReal, IsotopeCluster>::iterator it = iso_map_.begin();
878  it != iso_map_.end();
879  ++it)
880  {
881  d.iso_map_iter = it;
882 #ifdef DEBUG_2D
883  std::cerr << "element: " << counter << std::endl;
884  std::cerr << "mz: " << it->first << std::endl << "rts: ";
885  for (Size i = 0; i < it->second.scans.size(); ++i)
886  std::cerr << it->second.scans[i] << "\n";
887  std::cerr << std::endl << "peaks: ";
888  IsotopeCluster::IndexSet::const_iterator iter = it->second.peaks.begin();
889  for (; iter != it->second.peaks.end(); ++iter)
890  std::cerr << ms_exp[iter->first].getRT() << " " << (ms_exp[iter->first][iter->second]).getMZ() << std::endl;
891  //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;
892  std::cerr << std::endl << std::endl;
893 
894 #endif
895  // prepare for optimization:
896  // determine the matching peaks
897  // and the endpoints of each isotope pattern in the cluster
898 
899  getRegionEndpoints_(ms_exp, first, last, counter, 400, d);
900  OptimizePick::Data data;
901 
902 
903  Size idx = 0;
904  for (Size i = 0; i < d.signal2D.size() / 2; ++i)
905  {
906  data.positions.clear();
907  data.signal.clear();
908 
910  (d.raw_data_first + d.signal2D[2 * i].first)->begin() + d.signal2D[2 * i].second;
911  Int size = distance(ms_it, (d.raw_data_first + d.signal2D[2 * i].first)->begin() + d.signal2D[2 * i + 1].second);
912  data.positions.reserve(size);
913  data.signal.reserve(size);
914 
915  while (ms_it != (d.raw_data_first + d.signal2D[2 * i].first)->begin() + d.signal2D[2 * i + 1].second)
916  {
917  data.positions.push_back(ms_it->getMZ());
918  data.signal.push_back(ms_it->getIntensity());
919  ++ms_it;
920  }
921 
922 
924  pair.first = d.iso_map_iter->second.peaks.begin()->first + idx;
925 
926  IsotopeCluster::IndexSet::const_iterator set_iter = lower_bound(d.iso_map_iter->second.peaks.begin(),
927  d.iso_map_iter->second.peaks.end(),
929 
930 
931  // find the last entry with this rt-value
932  ++pair.first;
933  IsotopeCluster::IndexSet::const_iterator set_iter2 = lower_bound(d.iso_map_iter->second.peaks.begin(),
934  d.iso_map_iter->second.peaks.end(),
936 
937  while (set_iter != set_iter2)
938  {
939  const Size peak_index = set_iter->second;
940  const MSSpectrum<>& spec = ms_exp[set_iter->first];
941  PeakShape shape(spec.getFloatDataArrays()[1][peak_index], //intensity
942  spec[peak_index].getMZ(),
943  spec.getFloatDataArrays()[3][peak_index], //left width
944  spec.getFloatDataArrays()[4][peak_index], //right width
945  spec[peak_index].getIntensity(), //area is stored in peak intensity
946  std::vector<Peak1D>::iterator(),
947  std::vector<Peak1D>::iterator(),
948  PeakShape::Type(Int(spec.getFloatDataArrays()[5][peak_index]))); //shape
949  peak_shapes.push_back(shape);
950  ++set_iter;
951  }
952 #ifdef DEBUG_2D
953  std::cout << "rt "
954  << (d.raw_data_first + d.signal2D[2 * i].first)->getRT()
955  << "\n";
956 #endif
957  OptimizePick opt(penalties, max_iteration, eps_abs, eps_rel);
958 #ifdef DEBUG_2D
959  std::cout << "vorher\n";
960 
961  for (Size p = 0; p < peak_shapes.size(); ++p)
962  {
963  std::cout << peak_shapes[p].mz_position << "\t" << peak_shapes[p].height
964  << "\t" << peak_shapes[p].left_width << "\t" << peak_shapes[p].right_width << std::endl;
965  }
966 #endif
967  opt.optimize(peak_shapes, data);
968 #ifdef DEBUG_2D
969  std::cout << "nachher\n";
970  for (Size p = 0; p < peak_shapes.size(); ++p)
971  {
972  std::cout << peak_shapes[p].mz_position << "\t" << peak_shapes[p].height
973  << "\t" << peak_shapes[p].left_width << "\t" << peak_shapes[p].right_width << std::endl;
974  }
975 #endif
976  std::sort(peak_shapes.begin(), peak_shapes.end(), PeakShape::PositionLess());
977  pair.first = d.iso_map_iter->second.peaks.begin()->first + idx;
978 
979  set_iter = lower_bound(d.iso_map_iter->second.peaks.begin(),
980  d.iso_map_iter->second.peaks.end(),
982  Size p = 0;
983  while (p < peak_shapes.size())
984  {
985  MSSpectrum<>& spec = ms_exp[set_iter->first];
986  spec[set_iter->second].setMZ(peak_shapes[p].mz_position);
987  spec.getFloatDataArrays()[3][set_iter->second] = peak_shapes[p].left_width;
988  spec.getFloatDataArrays()[4][set_iter->second] = peak_shapes[p].right_width;
989  spec.getFloatDataArrays()[1][set_iter->second] = peak_shapes[p].height; // maximum intensity
990  // calculate area
991  if (peak_shapes[p].type == PeakShape::LORENTZ_PEAK)
992  {
993  PeakShape& ps = peak_shapes[p];
994  double x_left_endpoint = ps.mz_position - 1 / ps.left_width* sqrt(ps.height / 1 - 1);
995  double x_rigth_endpoint = ps.mz_position + 1 / ps.right_width* sqrt(ps.height / 1 - 1);
996  double area_left = -ps.height / ps.left_width* atan(ps.left_width * (x_left_endpoint - ps.mz_position));
997  double area_right = -ps.height / ps.right_width* atan(ps.right_width * (ps.mz_position - x_rigth_endpoint));
998  spec[set_iter->second].setIntensity(area_left + area_right); // area is stored as peak intensity
999  }
1000  else //It's a Sech - Peak
1001  {
1002  PeakShape& ps = peak_shapes[p];
1003  double x_left_endpoint = ps.mz_position - 1 / ps.left_width* boost::math::acosh(sqrt(ps.height / 0.001));
1004  double x_rigth_endpoint = ps.mz_position + 1 / ps.right_width* boost::math::acosh(sqrt(ps.height / 0.001));
1005  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)));
1006  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)));
1007  spec[set_iter->second].setIntensity(area_left + area_right); // area is stored as peak intensity
1008  }
1009  ++set_iter;
1010  ++p;
1011  }
1012  ++idx;
1013  peak_shapes.clear();
1014  }
1015 
1016  ++counter;
1017  }
1018  }
1019 
1020  template <typename InputSpectrumIterator, typename OutputPeakType>
1022  InputSpectrumIterator& first,
1023  InputSpectrumIterator& last,
1024  Size iso_map_idx,
1025  DoubleReal noise_level,
1027  {
1028  d.signal2D.clear();
1029  typedef typename InputSpectrumIterator::value_type InputExperimentType;
1030  typedef typename InputExperimentType::value_type InputPeakType;
1031  typedef std::multimap<DoubleReal, IsotopeCluster> MapType;
1032 
1033  DoubleReal rt, first_peak_mz, last_peak_mz;
1034 
1035  //MSSpectrum<InputPeakType> spec;
1037  InputPeakType peak;
1038 
1039  MapType::iterator iso_map_iter = iso_map_.begin();
1040  for (Size i = 0; i < iso_map_idx; ++i)
1041  ++iso_map_iter;
1042 
1043 #ifdef DEBUG2D
1044  std::cout << "rt begin: " << exp[iso_map_iter->second.scans[0]].getRT()
1045  << "\trt end: " << exp[iso_map_iter->second.scans[iso_map_iter->second.scans.size() - 1]].getRT()
1046  << " \t" << iso_map_iter->second.scans.size() << " scans"
1047  << std::endl;
1048 #endif
1049 
1050  // get left and right endpoint for all scans in the current cluster
1051  for (Size i = 0; i < iso_map_iter->second.scans.size(); ++i)
1052  {
1054 
1055  // first the right scan through binary search
1056  rt = exp[iso_map_iter->second.scans[i]].getRT();
1057  spec.setRT(rt);
1058  InputSpectrumIterator iter = lower_bound(first, last, spec, typename MSSpectrum<InputPeakType>::RTLess());
1059  // if(iter->getRT() != rt) --iter;
1060  exp_it = exp.RTBegin(rt);
1061 #ifdef DEBUG2D
1062  std::cout << exp_it->getRT() << " vs " << iter->getRT() << std::endl;
1063 #endif
1064  // now the right mz
1066  pair.first = iso_map_iter->second.peaks.begin()->first + i;
1067  // get iterator in peaks-set that points to the first peak in the current scan
1068  IsotopeCluster::IndexSet::const_iterator set_iter = lower_bound(iso_map_iter->second.peaks.begin(),
1069  iso_map_iter->second.peaks.end(),
1071 
1072  // consider a bit more of the signal to the left
1073  first_peak_mz = (exp_it->begin() + set_iter->second)->getMZ() - 1;
1074 
1075  // find the last entry with this rt-value
1076  ++pair.first;
1077  IsotopeCluster::IndexSet::const_iterator set_iter2 = lower_bound(iso_map_iter->second.peaks.begin(),
1078  iso_map_iter->second.peaks.end(),
1080 
1081  if (i == iso_map_iter->second.scans.size() - 1)
1082  {
1083  set_iter2 = iso_map_iter->second.peaks.end();
1084  --set_iter2;
1085  }
1086  else if (set_iter2 != iso_map_iter->second.peaks.begin())
1087  --set_iter2;
1088 
1089  last_peak_mz = (exp_it->begin() + set_iter2->second)->getMZ() + 1;
1090 
1091  //std::cout << rt<<": first peak mz "<<first_peak_mz << "\tlast peak mz "<<last_peak_mz <<std::endl;
1092  peak.setPosition(first_peak_mz);
1094  = lower_bound(iter->begin(), iter->end(), peak, typename InputPeakType::PositionLess());
1095  if (raw_data_iter != iter->begin())
1096  {
1097  --raw_data_iter;
1098  }
1099  DoubleReal intensity = raw_data_iter->getIntensity();
1100  // while the intensity is falling go to the left
1101  while (raw_data_iter != iter->begin() && (raw_data_iter - 1)->getIntensity() < intensity &&
1102  (raw_data_iter - 1)->getIntensity() > noise_level)
1103  {
1104  --raw_data_iter;
1105  intensity = raw_data_iter->getIntensity();
1106  }
1107  ++raw_data_iter;
1108  IsotopeCluster::IndexPair left, right;
1109  left.first = distance(first, iter);
1110  left.second = raw_data_iter - iter->begin();
1111 #ifdef DEBUG2D
1112  std::cout << "left: " << iter->getRT() << "\t" << raw_data_iter->getMZ() << std::endl;
1113 #endif
1114  // consider a bit more of the signal to the right
1115  peak.setPosition(last_peak_mz + 1);
1116  raw_data_iter
1117  = upper_bound(iter->begin(), iter->end(), peak, typename InputPeakType::PositionLess());
1118  if (raw_data_iter == iter->end())
1119  --raw_data_iter;
1120  intensity = raw_data_iter->getIntensity();
1121  // while the intensity is falling go to the right
1122  while (raw_data_iter + 1 != iter->end() && (raw_data_iter + 1)->getIntensity() < intensity)
1123  {
1124  ++raw_data_iter;
1125  intensity = raw_data_iter->getIntensity();
1126  if ((raw_data_iter + 1 != iter->end()) && (raw_data_iter + 1)->getIntensity() > noise_level)
1127  break;
1128  }
1129  right.first = left.first;
1130  right.second = raw_data_iter - iter->begin();
1131 #ifdef DEBUG2D
1132  std::cout << "right: " << iter->getRT() << "\t" << raw_data_iter->getMZ() << std::endl;
1133 #endif
1134  // region endpoints are stored in global vector
1135  d.signal2D.push_back(left);
1136  d.signal2D.push_back(right);
1137  }
1138 #ifdef DEBUG2D
1139  //std::cout << "fertig"<< std::endl;
1140  std::cout << first_peak_mz << "\t" << last_peak_mz << std::endl;
1141 #endif
1142  }
1143 
1144 }
1145 
1146 #endif //OPENMS_TRANSFORMATIONS_RAW2PEAK_TWODOPTIMIZATION_H
std::vector< double > positions
Positions and intensity values of the raw data.
Definition: OptimizePick.h:123
const double k
Internal representation of a peak shape (used by the PeakPickerCWT)
Definition: PeakShape.h:51
void getRegionEndpoints_(MSExperiment< OutputPeakType > &exp, InputSpectrumIterator &first, InputSpectrumIterator &last, Size iso_map_idx, DoubleReal noise_level, TwoDOptimization::Data &d)
Get the indices of the first and last raw data point of this region.
Definition: TwoDOptimization.h:1021
Definition: OptimizePick.h:120
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:212
Size size() const
Definition: MSExperiment.h:117
void setMZTolerance(DoubleReal tolerance_mz)
Mutable access to the matching epsilon.
Definition: TwoDOptimization.h:107
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:150
DoubleReal getMZTolerance() const
Non-mutable access to the matching epsilon.
Definition: TwoDOptimization.h:105
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:148
Comparison of mz_positions.
Definition: PeakShape.h:154
DoubleReal height
Maximum intensity of the peak shape.
Definition: PeakShape.h:131
DoubleReal eps_rel_
Convergence Parameter: Maximal relative error.
Definition: TwoDOptimization.h:219
Comparator for the retention time.
Definition: MSSpectrum.h:94
unsigned int UInt
Unsigned integer type.
Definition: Types.h:92
static Int residual2D_(const gsl_vector *x, void *params, gsl_vector *f)
Function computing estimated signal and its deviation to the experimental signal*/.
Peak2D PeakType
Definition: MassTrace.h:49
Helper struct (contains the size of an area and a raw data container)
Definition: TwoDOptimization.h:185
std::vector< DoubleReal >::iterator searchInScan_(std::vector< DoubleReal >::iterator scan_begin, std::vector< DoubleReal >::iterator scan_end, DoubleReal current_mz)
Iterator begin()
Definition: MSExperiment.h:147
virtual ~TwoDOptimization()
Destructor.
Definition: TwoDOptimization.h:98
std::multimap< DoubleReal, IsotopeCluster > iso_map_
stores the retention time of each isotopic cluster
Definition: TwoDOptimization.h:199
ContainerType::iterator Iterator
Mutable iterator.
Definition: MSSpectrum.h:123
double lWidth
Penalty factor for the peak shape&#39;s left width parameter.
Definition: OptimizePick.h:91
DoubleReal tolerance_mz_
threshold for the difference in the peak position of two matching peaks
Definition: TwoDOptimization.h:208
void setMaxRelError(DoubleReal eps_rel)
Mutable access to the maximal relative error.
Definition: TwoDOptimization.h:134
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:584
std::multimap< DoubleReal, IsotopeCluster >::iterator iso_map_iter
Definition: TwoDOptimization.h:188
OptimizationFunctions::PenaltyFactorsIntensity penalties
Definition: TwoDOptimization.h:193
DoubleReal eps_abs_
Convergence Parameter: Maximal absolute error.
Definition: TwoDOptimization.h:216
DoubleReal getMaxPeakDistance() const
Non-mutable access to the maximal peak distance in a cluster.
Definition: TwoDOptimization.h:114
Iterator end()
Definition: MSExperiment.h:157
DoubleReal getMaxRelError() const
Non-mutable access to the maximal relative error.
Definition: TwoDOptimization.h:132
ConstIterator RTBegin(CoordinateType rt) const
Fast search for spectrum range begin.
Definition: MSExperiment.h:349
void setMaxAbsError(DoubleReal eps_abs)
Mutable access to the maximal absolute error.
Definition: TwoDOptimization.h:125
Base::iterator iterator
Definition: MSExperiment.h:114
void findMatchingPeaks_(std::multimap< DoubleReal, IsotopeCluster >::iterator &it, MSExperiment<> &ms_exp)
Identify matching peak in a peak cluster.
static Int jacobian2D_(const gsl_vector *x, void *params, gsl_matrix *J)
Function computing the Jacobian */.
Type
Peak shape type (asymmetric lorentzian or asymmetric hyperbolic secans squared).
Definition: PeakShape.h:68
void setPenalties(OptimizationFunctions::PenaltyFactorsIntensity &penalties)
Mutable access to the minimal number of adjacent scans.
Definition: TwoDOptimization.h:152
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.
void optimize(std::vector< PeakShape > &peaks, Data &data)
Start the optimization of the peak shapes peaks. The original peak shapes will be subsituted by the o...
bool real_2D_
Optimization considering all scans of a cluster or optimization of each scan separately.
Definition: TwoDOptimization.h:225
DoubleReal max_peak_distance_
upper bound for distance between two peaks belonging to the same region
Definition: TwoDOptimization.h:205
Definition: PeakShape.h:70
Size total_nr_peaks
Definition: TwoDOptimization.h:189
std::vector< DoubleReal > signal
Definition: TwoDOptimization.h:195
void setMaxPeakDistance(DoubleReal max_peak_distance)
Mutable access to the maximal peak distance in a cluster.
Definition: TwoDOptimization.h:116
OptimizationFunctions::PenaltyFactorsIntensity penalties_
Penalty factors for some parameters in the optimization.
Definition: TwoDOptimization.h:229
std::vector< std::pair< SignedSize, SignedSize > > signal2D
Definition: TwoDOptimization.h:187
DoubleReal right_width
Right width parameter.
Definition: PeakShape.h:137
std::vector< DoubleReal > positions
Definition: TwoDOptimization.h:194
MSExperiment< Peak1D >::ConstIterator raw_data_first
Definition: TwoDOptimization.h:192
double pos
Penalty factor for the peak shape&#39;s position.
Definition: OptimizePick.h:89
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
std::multimap< DoubleReal, IsotopeCluster >::const_iterator curr_region_
Pointer to the current region.
Definition: TwoDOptimization.h:202
MSExperiment picked_peaks
Definition: TwoDOptimization.h:191
double rWidth
Penalty factor for the peak shape&#39;s right width parameter.
Definition: OptimizePick.h:93
void setRT(DoubleReal rt)
Sets the absolute retention time (is seconds)
Definition: MSSpectrum.h:221
std::vector< Size > scans
the scans of this cluster
Definition: IsotopeCluster.h:75
DoubleReal left_width
Left width parameter.
Definition: PeakShape.h:135
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
bool empty() const
Definition: MSExperiment.h:127
void setMaxIterations(UInt max_iteration)
Mutable access to the maximal number of iterations.
Definition: TwoDOptimization.h:143
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:71
Int charge
charge estimate (convention: zero means &quot;no charge estimate&quot;)
Definition: IsotopeCluster.h:62
UInt max_iteration_
Convergence Parameter: Maximal number of iterations.
Definition: TwoDOptimization.h:222
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:90
std::vector< double > signal
Definition: OptimizePick.h:124
UInt getMaxIterations() const
Non-mutable access to the maximal number of iterations.
Definition: TwoDOptimization.h:141
DoubleReal getMaxAbsError() const
Non-mutable access to the maximal absolute error.
Definition: TwoDOptimization.h:123
std::map< Int, std::vector< PeakIndex > > matching_peaks
Definition: TwoDOptimization.h:190
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:270
void optimizeRegionsScanwise_(InputSpectrumIterator &first, InputSpectrumIterator &last, MSExperiment< OutputPeakType > &ms_exp)
Definition: TwoDOptimization.h:808
DoubleReal height
Definition: OptimizePeakDeconvolution.h:85
int Int
Signed integer type.
Definition: Types.h:100
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:286
static Int evaluate2D_(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J)
Function that calls residual2D and jacobian2D*/.
This class provides the non-linear optimization of the peak parameters.
Definition: OptimizePick.h:116
This class provides the two-dimensional optimization of the picked peak parameters.
Definition: TwoDOptimization.h:86
Class for the penalty factors used during the optimization.
Definition: OptimizePeakDeconvolution.h:66
double DoubleReal
Double-precision real type.
Definition: Types.h:118
DoubleReal mz_position
Centroid position.
Definition: PeakShape.h:133

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