| OpenMS 2.5.0

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