| OpenMS 2.5.0

IsotopeWaveletTransform.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 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/IsotopeWaveletConstants.h>
38 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/IsotopeWavelet.h>
39 #include <OpenMS/KERNEL/StandardTypes.h>
40 #include <OpenMS/KERNEL/FeatureMap.h>
41 #include <OpenMS/KERNEL/MSExperiment.h>
42 #include <OpenMS/KERNEL/MSSpectrum.h>
43 #include <OpenMS/CONCEPT/Exception.h>
44 #include <OpenMS/MATH/STATISTICS/LinearRegression.h>
45 #include <OpenMS/DATASTRUCTURES/ConstRefVector.h>
46 #include <cmath>
47 #include <cmath>
48 #include <boost/math/special_functions/bessel.hpp>
49 #include <vector>
50 #include <map>
51 #include <sstream>
52 #include <fstream>
53 #include <iomanip>
54
55 // This code has quite a few strange things in it triggering warnings which
56 // clutters the rest of the diagnostics
57 #pragma clang diagnostic push
58 #pragma clang diagnostic ignored "-Wfloat-equal"
59 #pragma clang diagnostic ignored "-Wconversion"
60 #pragma clang diagnostic ignored "-Wshorten-64-to-32"
61
62 namespace OpenMS
63 {
64
68 template <typename PeakType>
69 class IsotopeWaveletTransform
70 {
71 public:
72
73
75 struct BoxElement
76 {
77 double mz;
78 UInt c;
79 double score;
80 double intens;
81 double ref_intens;
82 double RT;
83 UInt RT_index;
84 UInt MZ_begin;
85 UInt MZ_end;
86 };
87
88 typedef std::multimap<UInt, BoxElement> Box;
89
90
94 class TransSpectrum
95 {
96 friend class IsotopeWaveletTransform;
97
98 public:
99
101 TransSpectrum() :
102 reference_(nullptr), trans_intens_(nullptr)
103 {
104 }
105
107 TransSpectrum(const MSSpectrum* reference) :
108 reference_(reference)
109 {
110 trans_intens_ = new std::vector<float>(reference_->size(), 0.0);
111 }
112
114 virtual ~TransSpectrum()
115 {
116 delete (trans_intens_);
117 }
118
119 virtual void destroy()
120 {
121 delete (trans_intens_);
122 trans_intens_ = NULL;
123 delete (reference_);
124 reference_ = NULL;
125 }
126
128 inline double getRT() const
129 {
130 return reference_->getRT();
131 }
132
134 inline double getMZ(const UInt i) const
135 {
136 return (*reference_)[i].getMZ();
137 }
138
140 inline double getRefIntensity(const UInt i) const
141 {
142 return (*reference_)[i].getIntensity();
143 }
144
146 inline double getTransIntensity(const UInt i) const
147 {
148 return (*trans_intens_)[i];
149 }
150
152 inline void setTransIntensity(const UInt i, const double intens)
153 {
154 (*trans_intens_)[i] = intens;
155 }
156
158 inline Size size() const
159 {
160 return trans_intens_->size();
161 }
162
164 inline const MSSpectrum* getRefSpectrum()
165 {
166 return reference_;
167 }
168
170 inline const MSSpectrum* getRefSpectrum() const
171 {
172 return reference_;
173 }
174
177 inline typename MSSpectrum::const_iterator MZBegin(const double mz) const
178 {
179 return reference_->MZBegin(mz);
180 }
181
184 inline typename MSSpectrum::const_iterator MZEnd(const double mz) const
185 {
186 return reference_->MZEnd(mz);
187 }
188
191 inline typename MSSpectrum::const_iterator end() const
192 {
193 return reference_->end();
194 }
195
198 inline typename MSSpectrum::const_iterator begin() const
199 {
200 return reference_->begin();
201 }
202
203 protected:
204
205 const MSSpectrum* reference_;
206 std::vector<float>* trans_intens_;
207
208 };
209
210
211
217 IsotopeWaveletTransform(const double min_mz, const double max_mz, const UInt max_charge, const Size max_scan_size = 0, const bool hr_data = false, const String intenstype = "ref");
218
220 virtual ~IsotopeWaveletTransform();
221
222
227 virtual void getTransform(MSSpectrum& c_trans, const MSSpectrum& c_ref, const UInt c);
228
233 virtual void getTransformHighRes(MSSpectrum& c_trans, const MSSpectrum& c_ref, const UInt c);
234
254 virtual void identifyCharge(const MSSpectrum& candidates, const MSSpectrum& ref, const UInt scan_index, const UInt c,
255 const double ampl_cutoff, const bool check_PPMs);
256
257 virtual void initializeScan(const MSSpectrum& c_ref, const UInt c = 0);
258
265 void updateBoxStates(const PeakMap& map, const Size scan_index, const UInt RT_interleave,
266 const UInt RT_votes_cutoff, const Int front_bound = -1, const Int end_bound = -1);
267
268
273 FeatureMap mapSeeds2Features(const PeakMap& map, const UInt RT_votes_cutoff);
274
276 virtual std::multimap<double, Box> getClosedBoxes()
277 { return closed_boxes_; }
278
279
284 inline double getLinearInterpolation(const typename MSSpectrum::const_iterator& left_iter, const double mz_pos, const typename MSSpectrum::const_iterator& right_iter)
285 {
286 return left_iter->getIntensity() + (right_iter->getIntensity() - left_iter->getIntensity()) / (right_iter->getMZ() - left_iter->getMZ()) * (mz_pos - left_iter->getMZ());
287 }
288
295 inline double getLinearInterpolation(const double mz_a, const double intens_a, const double mz_pos, const double mz_b, const double intens_b)
296 {
297 return intens_a + (intens_b - intens_a) / (mz_b - mz_a) * (mz_pos - mz_a);
298 }
299
300 inline double getSigma() const
301 {
302 return sigma_;
303 }
304
305 inline void setSigma(const double sigma)
306 {
307 sigma_ = sigma;
308 }
309
310 virtual void computeMinSpacing(const MSSpectrum& c_ref);
311
312 inline double getMinSpacing() const
313 {
314 return min_spacing_;
315 }
316
317 inline Size getMaxScanSize() const
318 {
319 return max_scan_size_;
320 }
321
322 protected:
323
324
327 IsotopeWaveletTransform();
328
329
330 inline void sampleTheCMarrWavelet_(const MSSpectrum& scan, const Int wavelet_length, const Int mz_index, const UInt charge);
331
332
339 virtual double scoreThis_(const TransSpectrum& candidate, const UInt peak_cutoff,
340 const double seed_mz, const UInt c, const double ampl_cutoff);
341
348 virtual double scoreThis_(const MSSpectrum& candidate, const UInt peak_cutoff,
349 const double seed_mz, const UInt c, const double ampl_cutoff);
350
351
359 virtual bool checkPositionForPlausibility_(const TransSpectrum& candidate, const MSSpectrum& ref, const double seed_mz,
360 const UInt c, const UInt scan_index, const bool check_PPMs, const double transintens, const double prev_score);
361
369 virtual bool checkPositionForPlausibility_(const MSSpectrum& candidate, const MSSpectrum& ref, const double seed_mz,
370 const UInt c, const UInt scan_index, const bool check_PPMs, const double transintens, const double prev_score);
371
372 virtual std::pair<double, double> checkPPMTheoModel_(const MSSpectrum& ref, const double c_mz, const UInt c);
373
374
376 inline double getAvIntens_(const TransSpectrum& scan);
378 inline double getAvIntens_(const MSSpectrum& scan);
379
381 inline double getSdIntens_(const TransSpectrum& scan, const double mean);
383 inline double getSdIntens_(const MSSpectrum& scan, const double mean);
384
395 virtual void push2Box_(const double mz, const UInt scan, UInt c, const double score,
396 const double intens, const double rt, const UInt MZ_begin, const UInt MZ_end, const double ref_intens);
397
413 virtual void push2TmpBox_(const double mz, const UInt scan, UInt charge, const double score,
414 const double intens, const double rt, const UInt MZ_begin, const UInt MZ_end);
415
421 inline double getAvMZSpacing_(const MSSpectrum& scan);
422
423
428 void clusterSeeds_(const TransSpectrum& candidates, const MSSpectrum& ref,
429 const UInt scan_index, const UInt c, const bool check_PPMs);
430
435 virtual void clusterSeeds_(const MSSpectrum& candidates, const MSSpectrum& ref,
436 const UInt scan_index, const UInt c, const bool check_PPMs);
437
438
443 void extendBox_(const PeakMap& map, const Box& box);
444
447 inline double peptideMassRule_(const double c_mass) const
448 {
449 double correction_fac = c_mass / Constants::PEPTIDE_MASS_RULE_BOUND;
450 double old_frac_mass = c_mass - (Int)(c_mass);
451 double new_mass = ((Int)(c_mass)) * (1. + Constants::PEPTIDE_MASS_RULE_FACTOR) - (Int)(correction_fac);
452 double new_frac_mass = new_mass - (Int)(new_mass);
453
454 if (new_frac_mass - old_frac_mass > 0.5)
455 {
456 new_mass -= 1.;
457 }
458
459 if (new_frac_mass - old_frac_mass < -0.5)
460 {
461 new_mass += 1.;
462 }
463
464 return new_mass;
465 }
466
470 inline double getPPMs_(const double mass_a, const double mass_b) const
471 {
472 return fabs(mass_a - mass_b) / (0.5 * (mass_a + mass_b)) * 1e6;
473 }
474
475 //internally used data structures for the sweep line algorithm
476 std::multimap<double, Box> open_boxes_, closed_boxes_, end_boxes_, front_boxes_; //double = average m/z position
477 std::vector<std::multimap<double, Box> >* tmp_boxes_; //for each charge we need a separate container
478
479 double av_MZ_spacing_, sigma_;
480 std::vector<double> c_mzs_, c_spacings_, psi_, prod_, xs_;
481 std::vector<double> interpol_xs_, interpol_ys_;
482
483 Size max_scan_size_;
484 UInt max_num_peaks_per_pattern_, max_charge_, data_length_;
485 bool hr_data_;
486 String intenstype_;
487 Int from_max_to_left_, from_max_to_right_;
488 std::vector<int> indices_;
489
490 double min_spacing_, max_mz_cutoff_;
491 std::vector<float> scores_, zeros_;
492 };
493
494 template <typename PeakType>
495 bool intensityComparator(const PeakType& a, const PeakType& b)
496 {
497 return a.getIntensity() > b.getIntensity();
498 }
499
500 template <typename PeakType>
501 bool intensityAscendingComparator(const PeakType& a, const PeakType& b)
502 {
503 return a.getIntensity() < b.getIntensity();
504 }
505
506 template <typename PeakType>
507 bool intensityPointerComparator(PeakType* a, PeakType* b)
508 {
509 return a->getIntensity() > b->getIntensity();
510 }
511
512 template <typename PeakType>
513 bool positionComparator(const PeakType& a, const PeakType& b)
514 {
515 return a.getMZ() < b.getMZ();
516 }
517
518 template <typename PeakType>
519 IsotopeWaveletTransform<PeakType>::IsotopeWaveletTransform()
520 {
521 tmp_boxes_ = new std::vector<std::multimap<double, Box> >(1);
522 av_MZ_spacing_ = 1;
523 max_scan_size_ = 0;
524 max_mz_cutoff_ = 3;
525 max_num_peaks_per_pattern_ = 3;
526 hr_data_ = false;
527 intenstype_ = "ref";
528 }
529
530 template <typename PeakType>
531 IsotopeWaveletTransform<PeakType>::IsotopeWaveletTransform(const double min_mz, const double max_mz, const UInt max_charge, const Size max_scan_size, const bool hr_data, String intenstype)
532 {
533 max_charge_ = max_charge;
534 max_scan_size_ = max_scan_size;
535 hr_data_ = hr_data;
536 intenstype_ = intenstype;
537 tmp_boxes_ = new std::vector<std::multimap<double, Box> >(max_charge);
538 if (max_scan_size <= 0) //only important for the CPU
539 {
540 IsotopeWavelet::init(max_mz, max_charge);
541 }
542
543 av_MZ_spacing_ = 1;
544 max_mz_cutoff_ = IsotopeWavelet::getMzPeakCutOffAtMonoPos(max_mz, max_charge);
545 max_num_peaks_per_pattern_ = IsotopeWavelet::getNumPeakCutOff(max_mz, max_charge);
546
547 Int size_estimate((Int)ceil(max_scan_size_ / (max_mz - min_mz)));
548 Int to_reserve((Int)ceil(size_estimate * max_num_peaks_per_pattern_ * Constants::IW_NEUTRON_MASS));
549 psi_.reserve(to_reserve); //The wavelet
550 prod_.reserve(to_reserve);
551 xs_.reserve(to_reserve);
552 interpol_xs_.resize(Constants::DEFAULT_NUM_OF_INTERPOLATION_POINTS);
553 interpol_ys_.resize(Constants::DEFAULT_NUM_OF_INTERPOLATION_POINTS);
554 }
555
556 template <typename PeakType>
557 IsotopeWaveletTransform<PeakType>::~IsotopeWaveletTransform()
558 {
559 delete (tmp_boxes_);
560 }
561
562 template <typename PeakType>
563 void IsotopeWaveletTransform<PeakType>::getTransform(MSSpectrum& c_trans, const MSSpectrum& c_ref, const UInt c)
564 {
565 Int spec_size((Int)c_ref.size());
566 //in the very unlikely case that size_t will not fit to int anymore this will be a problem of course
567 //for the sake of simplicity (we need here a signed int) we do not cast at every following comparison individually
568 UInt charge = c + 1;
569 double value, T_boundary_left, T_boundary_right, old, c_diff, current, old_pos, my_local_MZ, my_local_lambda, origin, c_mz;
570
571 for (Int my_local_pos = 0; my_local_pos < spec_size; ++my_local_pos)
572 {
573 value = 0; T_boundary_left = 0; T_boundary_right = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_ref[my_local_pos].getMZ(), charge) / (double)charge;
574 old = 0; old_pos = (my_local_pos - from_max_to_left_ - 1 >= 0) ? c_ref[my_local_pos - from_max_to_left_ - 1].getMZ() : c_ref[0].getMZ() - min_spacing_;
575 my_local_MZ = c_ref[my_local_pos].getMZ(); my_local_lambda = IsotopeWavelet::getLambdaL(my_local_MZ * charge);
576 c_diff = 0;
577 origin = -my_local_MZ + Constants::IW_QUARTER_NEUTRON_MASS / (double)charge;
578
579 for (Int current_conv_pos = std::max(0, my_local_pos - from_max_to_left_); c_diff < T_boundary_right; ++current_conv_pos)
580 {
581 if (current_conv_pos >= spec_size)
582 {
583 value += 0.5 * old * min_spacing_;
584 break;
585 }
586
587 c_mz = c_ref[current_conv_pos].getMZ();
588 c_diff = c_mz + origin;
589
590 //Attention! The +1. has nothing to do with the charge, it is caused by the wavelet's formula (tz1).
591 current = c_diff > T_boundary_left && c_diff <= T_boundary_right ? IsotopeWavelet::getValueByLambda(my_local_lambda, c_diff * charge + 1.) * c_ref[current_conv_pos].getIntensity() : 0;
592
593 value += 0.5 * (current + old) * (c_mz - old_pos);
594
595 old = current;
596 old_pos = c_mz;
597 }
598
599
600
601 c_trans[my_local_pos].setIntensity(value);
602 }
603 }
604
605 template <typename PeakType>
606 void IsotopeWaveletTransform<PeakType>::getTransformHighRes(MSSpectrum& c_trans, const MSSpectrum& c_ref, const UInt c)
607 {
608 Int spec_size((Int)c_ref.size());
609 //in the very unlikely case that size_t will not fit to int anymore this will be a problem of course
610 //for the sake of simplicity (we need here a signed int) we do not cast at every following comparison individually
611 UInt charge = c + 1;
612 double value, T_boundary_left, T_boundary_right, c_diff, current, my_local_MZ, my_local_lambda, origin, c_mz;
613
614 for (Int my_local_pos = 0; my_local_pos < spec_size; ++my_local_pos)
615 {
616 value = 0; T_boundary_left = 0; T_boundary_right = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_ref[my_local_pos].getMZ(), charge) / (double)charge;
617
618
619 my_local_MZ = c_ref[my_local_pos].getMZ(); my_local_lambda = IsotopeWavelet::getLambdaL(my_local_MZ * charge);
620 c_diff = 0;
621 origin = -my_local_MZ + Constants::IW_QUARTER_NEUTRON_MASS / (double)charge;
622
623 for (Int current_conv_pos = std::max(0, my_local_pos - from_max_to_left_); c_diff < T_boundary_right; ++current_conv_pos)
624 {
625 if (current_conv_pos >= spec_size)
626 {
627 break;
628 }
629
630 c_mz = c_ref[current_conv_pos].getMZ();
631 c_diff = c_mz + origin;
632
633 //Attention! The +1. has nothing to do with the charge, it is caused by the wavelet's formula (tz1).
634 current = c_diff > T_boundary_left && c_diff <= T_boundary_right ? IsotopeWavelet::getValueByLambda(my_local_lambda, c_diff * charge + 1.) * c_ref[current_conv_pos].getIntensity() : 0;
635
636 value += current;
637 }
638
639 c_trans[my_local_pos].setIntensity(value);
640 }
641 }
642
643 template <typename PeakType>
644 void IsotopeWaveletTransform<PeakType>::initializeScan(const MSSpectrum& c_ref, const UInt c)
645 {
646 data_length_ = (UInt) c_ref.size();
647 computeMinSpacing(c_ref);
648 Int wavelet_length = 0, quarter_length = 0;
649
650 if (hr_data_) //We have to check this separately, because the simply estimation for LowRes data is destroyed by large gaps
651 {
652 UInt c_mz_cutoff;
653 typename MSSpectrum::const_iterator start_iter, end_iter;
654 for (UInt i = 0; i < data_length_; ++i)
655 {
656 c_mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_ref[i].getMZ(), c + 1);
657 start_iter = c_ref.MZEnd(c_ref[i].getMZ());
658 end_iter = c_ref.MZBegin(c_ref[i].getMZ() + c_mz_cutoff);
659 wavelet_length = std::max((SignedSize) wavelet_length, distance(start_iter, end_iter) + 1);
660 end_iter = c_ref.MZEnd(c_ref[i].getMZ() - Constants::IW_QUARTER_NEUTRON_MASS / double(c + 1.));
661 quarter_length = std::max((SignedSize) quarter_length, distance(end_iter, start_iter) + 1);
662 }
663 }
664 else
665 {
666 //CHANGED
667 max_mz_cutoff_ = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_ref[data_length_ - 1].getMZ(), max_charge_);
668 wavelet_length = (UInt) ceil(max_mz_cutoff_ / min_spacing_);
669 }
670 //... done
671
672
673 if (wavelet_length > (Int) c_ref.size())
674 {
675 std::cout << "Warning: the extremal length of the wavelet is larger (" << wavelet_length << ") than the number of data points (" << c_ref.size() << "). This might (!) severely affect the transform." << std::endl;
676 std::cout << "Minimal spacing: " << min_spacing_ << std::endl;
677 std::cout << "Warning/Error generated at scan with RT " << c_ref.getRT() << "." << std::endl;
678 }
679
680 Int max_index = (UInt) (Constants::IW_QUARTER_NEUTRON_MASS / min_spacing_);
681 from_max_to_left_ = max_index;
682 from_max_to_right_ = wavelet_length - 1 - from_max_to_left_;
683 }
684
685 template <typename PeakType>
686 void IsotopeWaveletTransform<PeakType>::computeMinSpacing(const MSSpectrum& c_ref)
687 {
688 min_spacing_ = INT_MAX;
689 for (UInt c_conv_pos = 1; c_conv_pos < c_ref.size(); ++c_conv_pos)
690 {
691 min_spacing_ = std::min(min_spacing_, c_ref[c_conv_pos].getMZ() - c_ref[c_conv_pos - 1].getMZ());
692 }
693 }
694
695 template <typename PeakType>
696 void IsotopeWaveletTransform<PeakType>::identifyCharge(const MSSpectrum& candidates,
697 const MSSpectrum& ref, const UInt scan_index, const UInt c, const double ampl_cutoff, const bool check_PPMs)
698 {
699 Size scan_size(candidates.size());
700 typename ConstRefVector<MSSpectrum >::iterator iter;
701 typename MSSpectrum::const_iterator iter_start, iter_end, iter_p, seed_iter, iter2;
702 double mz_cutoff, seed_mz, c_av_intens = 0, c_score = 0, c_sd_intens = 0, threshold = 0, help_mz, share, share_pos, bwd, fwd;
703 UInt MZ_start, MZ_end;
704
705 MSSpectrum diffed(candidates);
706 diffed[0].setIntensity(0); diffed[scan_size - 1].setIntensity(0);
707
708 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
709 std::stringstream stream;
710 stream << "diffed_" << ref.getRT() << "_" << c + 1 << ".trans\0";
711 std::ofstream ofile(stream.str().c_str());
712 #endif
713
714 if (!hr_data_) //LowRes data
715 {
716 for (UInt i = 0; i < scan_size - 2; ++i)
717 {
718 share = candidates[i + 1].getIntensity(); share_pos = candidates[i + 1].getMZ();
719 bwd = (share - candidates[i].getIntensity()) / (share_pos - candidates[i].getMZ());
720 fwd = (candidates[i + 2].getIntensity() - share) / (candidates[i + 2].getMZ() - share_pos);
721
722 if (!(bwd >= 0 && fwd <= 0) || share > ref[i + 1].getIntensity())
723 {
724 diffed[i + 1].setIntensity(0);
725 }
726
727 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
728 ofile << diffed[i + 1].getMZ() << "\t" << diffed[i + 1].getIntensity() << std::endl;
729 #endif
730 }
731 }
732 else //HighRes data
733 {
734 for (UInt i = 0; i < scan_size - 2; ++i)
735 {
736 share = candidates[i + 1].getIntensity(); share_pos = candidates[i + 1].getMZ();
737 bwd = (share - candidates[i].getIntensity()) / (share_pos - candidates[i].getMZ());
738 fwd = (candidates[i + 2].getIntensity() - share) / (candidates[i + 2].getMZ() - share_pos);
739
740 if (!(bwd >= 0 && fwd <= 0))
741 {
742 diffed[i + 1].setIntensity(0);
743 }
744
745 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
746 ofile << diffed[i + 1].getMZ() << "\t" << diffed[i + 1].getIntensity() << std::endl;
747 #endif
748 }
749 }
750 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
751 ofile.close();
752 #endif
753
754 ConstRefVector<MSSpectrum > c_sorted_candidate(diffed.begin(), diffed.end());
755
756 //Sort the transform in descending order according to the intensities present in the transform
757 c_sorted_candidate.sortByIntensity();
758
759 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
760 std::stringstream stream2;
761 stream2 << "sorted_cpu_" << candidates.getRT() << "_" << c + 1 << ".trans\0";
762 std::ofstream ofile2(stream2.str().c_str());
763 for (iter = c_sorted_candidate.end() - 1; iter != c_sorted_candidate.begin(); --iter)
764 {
765 ofile2 << iter->getMZ() << "\t" << iter->getIntensity() << std::endl;
766 }
767 ofile2.close();
768 #endif
769
770 std::vector<UInt> processed(scan_size, 0);
771
772 if (ampl_cutoff < 0)
773 {
774 threshold = 0;
775 }
776 else
777 {
778 c_av_intens = getAvIntens_(candidates);
779 c_sd_intens = getSdIntens_(candidates, c_av_intens);
780 threshold = ampl_cutoff * c_sd_intens + c_av_intens;
781 }
782
783 for (iter = c_sorted_candidate.end() - 1; iter != c_sorted_candidate.begin(); --iter)
784 {
785 if (iter->getIntensity() <= 0)
786 {
787 break;
788 }
789
790 seed_mz = iter->getMZ();
791 seed_iter = ref.MZBegin(seed_mz);
792
793 if (seed_iter == ref.end() || processed[distance(ref.begin(), seed_iter)])
794 {
795 continue;
796 }
797
798 mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(seed_mz, c + 1);
799 //Mark the region as processed
800 //Do not move this further down, since we have to mark this as processed in any case,
801 //even when score <=0; otherwise we would look around the maximum's position unless
802 //any significant point is found
803 iter_start = ref.MZBegin(ref.begin(), seed_mz - Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.), seed_iter);
804 iter_end = ref.MZEnd(seed_iter, seed_mz + mz_cutoff / (c + 1.), ref.end());
805 if (iter_end == ref.end())
806 {
807 --iter_end;
808 }
809
810 MZ_start = distance(ref.begin(), iter_start);
811 MZ_end = distance(ref.begin(), iter_end);
812
813 memset(&(processed[MZ_start]), 1, sizeof(UInt) * (MZ_end - MZ_start + 1));
814
815 c_score = scoreThis_(candidates, IsotopeWavelet::getNumPeakCutOff(seed_mz * (c + 1.)), seed_mz, c, threshold);
816
817 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
818 if (trunc(seed_mz) == 874)
819 std::cout << seed_mz << "\t" << c_score << std::endl;
820 #endif
821
822 if (c_score <= 0 && c_score != -1000)
823 {
824 continue;
825 }
826
827 //Push the seed into its corresponding box (or create a new one, if necessary)
828 //Do ***NOT*** move this further down!
829 push2TmpBox_(seed_mz, scan_index, c, c_score, iter->getIntensity(), ref.getRT(), MZ_start, MZ_end);
830
831 help_mz = seed_mz - Constants::IW_NEUTRON_MASS / (c + 1.);
832 iter2 = candidates.MZBegin(help_mz);
833 if (iter2 == candidates.end() || iter2 == candidates.begin())
834 {
835 continue;
836 }
837
838 if (fabs(iter2->getMZ() - seed_mz) > 0.5 * Constants::IW_NEUTRON_MASS / (c + 1.))
839 {
840 //In the other case, we are too close to the peak, leading to incorrect derivatives.
841 if (iter2 != candidates.end())
842 {
843 push2TmpBox_(iter2->getMZ(), scan_index, c, 0, getLinearInterpolation(iter2 - 1, help_mz, iter2), ref.getRT(), MZ_start, MZ_end);
844 }
845 }
846
847 help_mz = seed_mz + Constants::IW_NEUTRON_MASS / (c + 1.);
848 iter2 = candidates.MZBegin(help_mz);
849 if (iter2 == candidates.end() || iter2 == candidates.begin())
850 {
851 continue;
852 }
853
854 if (fabs(iter2->getMZ() - seed_mz) > 0.5 * Constants::IW_NEUTRON_MASS / (c + 1.))
855 {
856 //In the other case, we are too close to the peak, leading to incorrect derivatives.
857 if (iter2 != candidates.end())
858 {
859 push2TmpBox_(iter2->getMZ(), scan_index, c, 0, getLinearInterpolation(iter2 - 1, help_mz, iter2), ref.getRT(), MZ_start, MZ_end);
860 }
861 }
862 }
863
864 clusterSeeds_(candidates, ref, scan_index, c, check_PPMs);
865 }
866
867 template <typename PeakType>
868 double IsotopeWaveletTransform<PeakType>::scoreThis_(const MSSpectrum& candidate,
869 const UInt peak_cutoff, const double seed_mz, const UInt c, const double ampl_cutoff)
870 {
871 double c_score = 0, c_val;
872 typename MSSpectrum::const_iterator c_left_iter2, c_right_iter2;
873 Int signal_size((Int)candidate.size());
874 //in the very unlikely case that size_t will not fit to int anymore this will be a problem of course
875 //for the sake of simplicity (we need here a signed int) we do not cast at every following comparison individually
876
877 //p_h_ind indicates if we are looking for a whole or a peak
878 Int p_h_ind = 1, end = 4 * (peak_cutoff - 1) - 1; //4 times and not 2 times, since we move by 0.5 m/z entities
879
880 std::vector<double> positions(end);
881 for (Int i = 0; i < end; ++i)
882 {
883 positions[i] = seed_mz - ((peak_cutoff - 1) * Constants::IW_NEUTRON_MASS - (i + 1) * Constants::IW_HALF_NEUTRON_MASS) / ((double)c + 1);
884 }
885
886 double l_score = 0, mid_val = 0;
887 Int start_index = distance(candidate.begin(), candidate.MZBegin(positions[0])) - 1;
888 for (Int v = 1; v <= end; ++v, ++p_h_ind)
889 {
890 do
891 {
892 if (start_index < signal_size - 1)
893 ++start_index;
894 else
895 break;
896 }
897 while (candidate[start_index].getMZ() < positions[v - 1]);
898
899 if (start_index <= 0 || start_index >= signal_size - 1) //unable to interpolate
900 {
901 continue;
902 }
903
904 c_left_iter2 = candidate.begin() + start_index - 1;
905 c_right_iter2 = c_left_iter2 + 1;
906
907 c_val = c_left_iter2->getIntensity() + (c_right_iter2->getIntensity() - c_left_iter2->getIntensity()) / (c_right_iter2->getMZ() - c_left_iter2->getMZ()) * (positions[v - 1] - c_left_iter2->getMZ());
908
909 if (v == (int)(ceil(end / 2.)))
910 {
911 l_score = c_score;
912 mid_val = c_val;
913 }
914
915 if (p_h_ind % 2 == 1) //I.e. a whole
916 {
917 c_score -= c_val;
918 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
919 if (trunc(seed_mz) == 874)
920 std::cout << -c_val << std::endl;
921 #endif
922 }
923 else
924 {
925 c_score += c_val;
926 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
927 if (trunc(seed_mz) == 874)
928 std::cout << c_val << std::endl;
929 #endif
930 }
931
932
933 start_index = distance(candidate.begin(), c_left_iter2);
934 }
935
936 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
937 std::ofstream ofile_score("scores.dat", ios::app);
938 std::ofstream ofile_check_score("check_scores.dat", ios::app);
939 ofile_score.close();
940 ofile_check_score.close();
941 #endif
942
943 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
944 if (trunc(seed_mz) == 874)
945 std::cout << "final_score: " << seed_mz << "\t" << c_score << "\t l_score: " << l_score << "\t" << c_score - l_score - mid_val << "\t" << c_score - mid_val << "\t" << ampl_cutoff << std::endl;
946 #endif
947
948 if (c_score - mid_val <= 0)
949 {
950 return 0;
951 }
952
953 if (c_score - mid_val <= ampl_cutoff)
954 {
955 return -1000;
956 }
957
958 if (l_score <= 0 || c_score - l_score - mid_val <= 0)
959 {
960 return 0;
961 }
962
963 return c_score;
964 }
965
966 template <typename PeakType>
967 double IsotopeWaveletTransform<PeakType>::scoreThis_(const TransSpectrum& candidate,
968 const UInt peak_cutoff, const double seed_mz, const UInt c, const double ampl_cutoff)
969 {
970 double c_score = 0, c_val;
971 typename MSSpectrum::const_iterator c_left_iter2, c_right_iter2;
972 Int signal_size((Int)candidate.size());
973
974 //p_h_ind indicates if we are looking for a whole or a peak
975 Int p_h_ind = 1, end = 4 * (peak_cutoff - 1) - 1; //4 times and not 2 times, since we move by 0.5 m/z entities
976
977 std::vector<double> positions(end);
978 for (Int i = 0; i < end; ++i)
979 {
980 positions[i] = seed_mz - ((peak_cutoff - 1) * Constants::IW_NEUTRON_MASS - (i + 1) * Constants::IW_HALF_NEUTRON_MASS) / ((double)c + 1);
981 }
982
983 double l_score = 0, mid_val = 0;
984 Int start_index = distance(candidate.begin(), candidate.MZBegin(positions[0])) - 1;
985 for (Int v = 1; v <= end; ++v, ++p_h_ind)
986 {
987 do
988 {
989 if (start_index < signal_size - 1)
990 ++start_index;
991 else
992 break;
993 }
994 while (candidate.getMZ(start_index) < positions[v - 1]);
995
996 if (start_index <= 0 || start_index >= signal_size - 1) //unable to interpolate
997 {
998 continue;
999 }
1000
1001 c_left_iter2 = candidate.begin() + start_index - 1;
1002 c_right_iter2 = c_left_iter2 + 1;
1003
1004 c_val = candidate.getTransIntensity(start_index - 1) + (candidate.getTransIntensity(start_index) - candidate.getTransIntensity(start_index - 1)) / (c_right_iter2->getMZ() - c_left_iter2->getMZ()) * (positions[v - 1] - c_left_iter2->getMZ());
1005 if (v == (int)(ceil(end / 2.)))
1006 {
1007 l_score = c_score;
1008 mid_val = c_val;
1009 }
1010
1011 if (p_h_ind % 2 == 1) //I.e. a whole
1012 {
1013 c_score -= c_val;
1014 }
1015 else
1016 {
1017 c_score += c_val;
1018 }
1019
1020 start_index = distance(candidate.begin(), c_left_iter2);
1021 }
1022
1023 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1024 std::ofstream ofile_score("scores.dat", ios::app);
1025 std::ofstream ofile_check_score("check_scores.dat", ios::app);
1026 ofile_score << c_check_point << "\t" << c_score << std::endl;
1027 ofile_score.close();
1028 ofile_check_score.close();
1029 #endif
1030
1031 if (l_score <= 0 || c_score - l_score - mid_val <= 0 || c_score - mid_val <= ampl_cutoff)
1032 {
1033 return 0;
1034 }
1035
1036 return c_score;
1037 }
1038
1039 template <typename PeakType>
1040 void IsotopeWaveletTransform<PeakType>::clusterSeeds_(const MSSpectrum& candidate,
1041 const MSSpectrum& ref, const UInt scan_index, const UInt c, const bool check_PPMs)
1042 {
1043 typename std::multimap<double, Box>::iterator iter;
1044 typename Box::iterator box_iter;
1045 std::vector<BoxElement> final_box;
1046 double c_mz, av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
1047 double virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
1048
1049 typename std::pair<double, double> c_extend;
1050 for (iter = tmp_boxes_->at(c).begin(); iter != tmp_boxes_->at(c).end(); ++iter)
1051 {
1052
1053 Box& c_box = iter->second;
1054 av_score = 0; av_mz = 0; av_intens = 0; av_abs_intens = 0; count = 0;
1055 virtual_av_mz = 0; virtual_av_intens = 0; virtual_av_abs_intens = 0; virtual_count = 0;
1056
1057 //Now, let's get the RT boundaries for the box
1058 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1059 {
1060 if (box_iter->second.score == 0) //virtual helping point
1061 {
1062 if (count != 0)
1063 continue; //it is in any way not pure virtual
1064
1065 c_mz = box_iter->second.mz;
1066 virtual_av_intens += box_iter->second.intens;
1067 virtual_av_abs_intens += fabs(box_iter->second.intens);
1068 virtual_av_mz += c_mz * fabs(box_iter->second.intens);
1069 ++virtual_count;
1070 }
1071 else
1072 {
1073 c_mz = box_iter->second.mz;
1074 av_score += box_iter->second.score;
1075 av_intens += box_iter->second.intens;
1076 av_abs_intens += fabs(box_iter->second.intens);
1077 av_mz += c_mz * fabs(box_iter->second.intens);
1078 ++count;
1079 }
1080 }
1081
1082 if (count == 0) //pure virtual helping box
1083 {
1084 av_intens = virtual_av_intens / virtual_count;
1085 av_score = 0;
1086 av_mz = virtual_av_mz / virtual_av_abs_intens;
1087 }
1088 else
1089 {
1090 av_intens /= count;
1091 av_score /= count;
1092 av_mz /= av_abs_intens;
1093 }
1094
1095 BoxElement c_box_element;
1096 c_box_element.mz = av_mz;
1097 c_box_element.c = c;
1098 c_box_element.score = av_score;
1099 c_box_element.intens = av_intens;
1100
1101 c_box_element.RT = c_box.begin()->second.RT;
1102 final_box.push_back(c_box_element);
1103 }
1104
1105 Size num_o_feature = final_box.size();
1106 if (num_o_feature == 0)
1107 {
1108 tmp_boxes_->at(c).clear();
1109 return;
1110 }
1111
1112 //Computing the derivatives
1113 std::vector<double> bwd_diffs(num_o_feature, 0);
1114
1115 bwd_diffs[0] = 0;
1116 for (Size i = 1; i < num_o_feature; ++i)
1117 {
1118 bwd_diffs[i] = (final_box[i].intens - final_box[i - 1].intens) / (final_box[i].mz - final_box[i - 1].mz);
1119 }
1120
1121 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1122 std::ofstream ofile_bwd("bwd_cpu.dat");
1123 for (Size i = 0; i < num_o_feature; ++i)
1124 {
1125 ofile_bwd << final_box[i].mz << "\t" << bwd_diffs[i] << std::endl;
1126 }
1127 ofile_bwd.close();
1128 #endif
1129
1130
1131 for (Size i = 0; i < num_o_feature - 1; ++i)
1132 {
1133 while (i < num_o_feature - 2)
1134 {
1135 if (final_box[i].score > 0 || final_box[i].score == -1000) //this has been an helping point
1136 break;
1137 ++i;
1138 }
1139
1140 if (bwd_diffs[i] > 0 && bwd_diffs[i + 1] < 0)
1141 {
1142 checkPositionForPlausibility_(candidate, ref, final_box[i].mz, final_box[i].c, scan_index, check_PPMs, final_box[i].intens, final_box[i].score);
1143 continue;
1144 }
1145 }
1146
1147 tmp_boxes_->at(c).clear();
1148 }
1149
1150 template <typename PeakType>
1151 double IsotopeWaveletTransform<PeakType>::getAvIntens_(const MSSpectrum& scan)
1152 {
1153 double av_intens = 0;
1154 for (UInt i = 0; i < scan.size(); ++i)
1155 {
1156 if (scan[i].getIntensity() >= 0)
1157 {
1158 av_intens += scan[i].getIntensity();
1159 }
1160 }
1161 return av_intens / (double)scan.size();
1162 }
1163
1164 template <typename PeakType>
1165 double IsotopeWaveletTransform<PeakType>::getSdIntens_(const MSSpectrum& scan, const double mean)
1166 {
1167 double res = 0, intens;
1168 for (UInt i = 0; i < scan.size(); ++i)
1169 {
1170 if (scan[i].getIntensity() >= 0)
1171 {
1172 intens = scan[i].getIntensity();
1173 res += (intens - mean) * (intens - mean);
1174 }
1175 }
1176 return sqrt(res / (double)(scan.size() - 1));
1177 }
1178
1179 template <typename PeakType>
1180 double IsotopeWaveletTransform<PeakType>::getAvMZSpacing_(const MSSpectrum& scan) //, Int start_index, Int end_index)
1181 {
1182 std::vector<double> diffs(scan.size() - 1, 0);
1183 for (UInt i = 0; i < scan.size() - 1; ++i)
1184 {
1185 diffs[i] = scan[i + 1].getMZ() - scan[i].getMZ();
1186 }
1187
1188 sort(diffs.begin(), diffs.end());
1189 double av_MZ_spacing = 0;
1190 for (UInt i = 0; i < diffs.size() / 2; ++i)
1191 {
1192 av_MZ_spacing += diffs[i];
1193 }
1194
1195 return av_MZ_spacing / (diffs.size() / 2);
1196 }
1197
1198 template <typename PeakType>
1199 double IsotopeWaveletTransform<PeakType>::getAvIntens_(const TransSpectrum& scan)
1200 {
1201 double av_intens = 0;
1202 for (UInt i = 0; i < scan.size(); ++i)
1203 {
1204 if (scan.getTransIntensity(i) >= 0)
1205 {
1206 av_intens += scan.getTransIntensity(i);
1207 }
1208 }
1209 return av_intens / (double)scan.size();
1210 }
1211
1212 template <typename PeakType>
1213 double IsotopeWaveletTransform<PeakType>::getSdIntens_(const TransSpectrum& scan, const double mean)
1214 {
1215 double res = 0, intens;
1216 for (UInt i = 0; i < scan.size(); ++i)
1217 {
1218 if (scan.getTransIntensity(i) >= 0)
1219 {
1220 intens = scan.getTransIntensity(i);
1221 res += (intens - mean) * (intens - mean);
1222 }
1223 }
1224 return sqrt(res / (double)(scan.size() - 1));
1225 }
1226
1227 template <typename PeakType>
1228 void IsotopeWaveletTransform<PeakType>::push2Box_(const double mz, const UInt scan, UInt c,
1229 const double score, const double intens, const double rt, const UInt MZ_begin, const UInt MZ_end, double ref_intens)
1230 {
1231 const double dist_constraint(Constants::IW_HALF_NEUTRON_MASS / (double)max_charge_);
1232
1233 typename std::multimap<double, Box>::iterator upper_iter(open_boxes_.upper_bound(mz));
1234 typename std::multimap<double, Box>::iterator lower_iter(open_boxes_.lower_bound(mz));
1235
1236 if (lower_iter != open_boxes_.end())
1237 {
1238 //Ugly, but necessary due to the implementation of STL lower_bound
1239 if (mz != lower_iter->first && lower_iter != open_boxes_.begin())
1240 {
1241 --lower_iter;
1242 }
1243 }
1244
1245 typename std::multimap<double, Box>::iterator insert_iter;
1246 bool create_new_box = true;
1247 if (lower_iter == open_boxes_.end()) //I.e. there is no open Box for that mz position
1248 {
1249 //There is another special case to be considered here:
1250 //Assume that the current box contains only a single element that is (slightly) smaller than the new mz value,
1251 //then the lower bound for the new mz value is box.end and this would usually force a new entry
1252 if (!open_boxes_.empty())
1253 {
1254 if (fabs((--lower_iter)->first - mz) < dist_constraint) //matching box
1255 {
1256 create_new_box = false;
1257 insert_iter = lower_iter;
1258 }
1259 }
1260 else
1261 {
1262 create_new_box = true;
1263 }
1264 }
1265 else
1266 {
1267 if (upper_iter == open_boxes_.end() && fabs(lower_iter->first - mz) < dist_constraint) //Found matching Box
1268 {
1269 insert_iter = lower_iter;
1270 create_new_box = false;
1271 }
1272 else
1273 {
1274 create_new_box = true;
1275 }
1276 }
1277
1278
1279 if (upper_iter != open_boxes_.end() && lower_iter != open_boxes_.end())
1280 {
1281 //Here is the question if you should figure out the smallest charge .... and then
1282
1283 //Figure out which entry is closer to m/z
1284 double dist_lower = fabs(lower_iter->first - mz);
1285 double dist_upper = fabs(upper_iter->first - mz);
1286 dist_lower = (dist_lower < dist_constraint) ? dist_lower : INT_MAX;
1287 dist_upper = (dist_upper < dist_constraint) ? dist_upper : INT_MAX;
1288
1289 if (dist_lower >= dist_constraint && dist_upper >= dist_constraint) // they are both too far away
1290 {
1291 create_new_box = true;
1292 }
1293 else
1294 {
1295 insert_iter = (dist_lower < dist_upper) ? lower_iter : upper_iter;
1296 create_new_box = false;
1297 }
1298 }
1299
1300 BoxElement element;
1301 element.c = c; element.mz = mz; element.score = score; element.RT = rt; element.intens = intens; element.ref_intens = ref_intens;
1302 element.RT_index = scan; element.MZ_begin = MZ_begin; element.MZ_end = MZ_end;
1303
1304
1305 if (create_new_box == false)
1306 {
1307 std::pair<UInt, BoxElement> help2(scan, element);
1308 insert_iter->second.insert(help2);
1309
1310 //Unfortunately, we need to change the m/z key to the average of all keys inserted in that box.
1311 Box replacement(insert_iter->second);
1312
1313 //We cannot divide both m/z by 2, since we already inserted some m/zs whose weight would be lowered.
1314 //Also note that we already inserted the new entry, leading to size-1.
1315 double c_mz = insert_iter->first * (insert_iter->second.size() - 1) + mz;
1316 c_mz /= ((double) insert_iter->second.size());
1317
1318 //Now let's remove the old and insert the new one
1319 open_boxes_.erase(insert_iter);
1320 std::pair<double, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
1321 open_boxes_.insert(help3);
1322 }
1323 else
1324 {
1325 std::pair<UInt, BoxElement> help2(scan, element);
1326 std::multimap<UInt, BoxElement> help3;
1327 help3.insert(help2);
1328 std::pair<double, std::multimap<UInt, BoxElement> > help4(mz, help3);
1329 open_boxes_.insert(help4);
1330 }
1331 }
1332
1333 template <typename PeakType>
1334 void IsotopeWaveletTransform<PeakType>::push2TmpBox_(const double mz, const UInt scan, UInt c,
1335 const double score, const double intens, const double rt, const UInt MZ_begin, const UInt MZ_end)
1336 {
1337 const double dist_constraint(Constants::IW_HALF_NEUTRON_MASS / (double)max_charge_);
1338
1339 std::multimap<double, Box>& tmp_box(tmp_boxes_->at(c));
1340 typename std::multimap<double, Box>::iterator upper_iter(tmp_box.upper_bound(mz));
1341 typename std::multimap<double, Box>::iterator lower_iter(tmp_box.lower_bound(mz));
1342
1343 if (lower_iter != tmp_box.end())
1344 {
1345 //Ugly, but necessary due to the implementation of STL lower_bound
1346 if (mz != lower_iter->first && lower_iter != tmp_box.begin())
1347 {
1348 --lower_iter;
1349 }
1350 }
1351
1352 typename std::multimap<double, Box>::iterator insert_iter;
1353 bool create_new_box = true;
1354 if (lower_iter == tmp_box.end()) //I.e. there is no tmp Box for that mz position
1355 {
1356 //There is another special case to be considered here:
1357 //Assume that the current box contains only a single element that is (slightly) smaller than the new mz value,
1358 //then the lower bound for the new mz value is box.end and this would usually force a new entry
1359 if (!tmp_box.empty())
1360 {
1361 if (fabs((--lower_iter)->first - mz) < dist_constraint) //matching box
1362 {
1363 create_new_box = false;
1364 insert_iter = lower_iter;
1365 }
1366 }
1367 else
1368 {
1369 create_new_box = true;
1370 }
1371 }
1372 else
1373 {
1374 if (upper_iter == tmp_box.end() && fabs(lower_iter->first - mz) < dist_constraint) //Found matching Box
1375 {
1376 insert_iter = lower_iter;
1377 create_new_box = false;
1378 }
1379 else
1380 {
1381 create_new_box = true;
1382 }
1383 }
1384
1385
1386 if (upper_iter != tmp_box.end() && lower_iter != tmp_box.end())
1387 {
1388 //Figure out which entry is closer to m/z
1389 double dist_lower = fabs(lower_iter->first - mz);
1390 double dist_upper = fabs(upper_iter->first - mz);
1391 dist_lower = (dist_lower < dist_constraint) ? dist_lower : INT_MAX;
1392 dist_upper = (dist_upper < dist_constraint) ? dist_upper : INT_MAX;
1393
1394 if (dist_lower >= dist_constraint && dist_upper >= dist_constraint) // they are both too far away
1395 {
1396 create_new_box = true;
1397 }
1398 else
1399 {
1400 insert_iter = (dist_lower < dist_upper) ? lower_iter : upper_iter;
1401 create_new_box = false;
1402 }
1403 }
1404
1405 BoxElement element;
1406 element.c = c; element.mz = mz; element.score = score; element.RT = rt; element.intens = intens; element.ref_intens = -1000;
1407 element.RT_index = scan; element.MZ_begin = MZ_begin; element.MZ_end = MZ_end;
1408
1409 if (create_new_box == false)
1410 {
1411 std::pair<UInt, BoxElement> help2(scan, element);
1412 insert_iter->second.insert(help2);
1413
1414 //Unfortunately, we need to change the m/z key to the average of all keys inserted in that box.
1415 Box replacement(insert_iter->second);
1416
1417 //We cannot divide both m/z by 2, since we already inserted some m/zs whose weight would be lowered.
1418 //Also note that we already inserted the new entry, leading to size-1.
1419 double c_mz = insert_iter->first * (insert_iter->second.size() - 1) + mz;
1420 c_mz /= ((double) insert_iter->second.size());
1421
1422 //Now let's remove the old and insert the new one
1423 tmp_box.erase(insert_iter);
1424 std::pair<double, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
1425 tmp_box.insert(help3);
1426 }
1427 else
1428 {
1429 std::pair<UInt, BoxElement> help2(scan, element);
1430 std::multimap<UInt, BoxElement> help3;
1431 help3.insert(help2);
1432
1433 std::pair<double, std::multimap<UInt, BoxElement> > help4(mz, help3);
1434 tmp_box.insert(help4);
1435 }
1436 }
1437
1438 template <typename PeakType>
1439 void IsotopeWaveletTransform<PeakType>::updateBoxStates(const PeakMap& map, const Size scan_index, const UInt RT_interleave,
1440 const UInt RT_votes_cutoff, const Int front_bound, const Int end_bound)
1441 {
1442 typename std::multimap<double, Box>::iterator iter, iter2;
1443
1444 if ((Int)scan_index == end_bound && end_bound != (Int)map.size() - 1)
1445 {
1446 for (iter = open_boxes_.begin(); iter != open_boxes_.end(); ++iter)
1447 {
1448 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1449 std::cout << "LOW THREAD insert in end_box " << iter->first << std::endl;
1450 typename Box::iterator dings;
1451 for (dings = iter->second.begin(); dings != iter->second.end(); ++dings)
1452 std::cout << map[dings->first].getRT() << "\t" << dings->second.c + 1 << std::endl;
1453 #endif
1454 end_boxes_.insert(*iter);
1455 }
1456 open_boxes_.clear();
1457 return;
1458 }
1459
1460 for (iter = open_boxes_.begin(); iter != open_boxes_.end(); )
1461 {
1462 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1463 if (front_bound > 0)
1464 {
1465 std::cout << "HIGH THREAD open box. " << iter->first << "\t current scan index " << scan_index << "\t" << ((iter->second.begin()))->first << "\t of last scan " << map.size() - 1 << "\t" << front_bound << std::endl;
1466 }
1467 #endif
1468
1469 //For each Box we need to figure out, if and when the last RT value has been inserted
1470 UInt lastScan = (--(iter->second.end()))->first;
1471 if (scan_index - lastScan > RT_interleave + 1 || scan_index == map.size() - 1) //I.e. close the box!
1472 {
1473 if (iter->second.begin()->first - front_bound <= RT_interleave + 1 && front_bound > 0)
1474 {
1475 iter2 = iter;
1476 ++iter2;
1477 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1478 std::cout << "HIGH THREAD insert in front_box " << iter->first << std::endl;
1479 #endif
1480 front_boxes_.insert(*iter);
1481 open_boxes_.erase(iter);
1482 iter = iter2;
1483 continue;
1484 }
1485
1486 iter2 = iter;
1487 ++iter2;
1488 //Please do **NOT** simplify the upcoming lines.
1489 //The 'obvious' overhead is necessary since the object represented by iter might be erased
1490 //by push2Box which might be called by extendBox_.
1491 if (iter->second.size() >= RT_votes_cutoff)
1492 {
1493 //extendBox_ (map, iter->second);
1494 iter = iter2;
1495 closed_boxes_.insert(*(--iter));
1496 }
1497 open_boxes_.erase(iter);
1498 iter = iter2;
1499 }
1500 else
1501 {
1502 ++iter;
1503 }
1504 }
1505 }
1506
1507 template <typename PeakType>
1508 void IsotopeWaveletTransform<PeakType>::extendBox_(const PeakMap& map, const Box& box)
1509 {
1510 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1511 std::cout << "**** CHECKING FOR BOX EXTENSIONS ****" << std::endl;
1512 #endif
1513
1514 //Determining the elution profile
1515 typename Box::const_iterator iter;
1516 std::vector<double> elution_profile(box.size());
1517 UInt index = 0;
1518 for (iter = box.begin(); iter != box.end(); ++iter, ++index)
1519 {
1520 for (Size i = iter->second.MZ_begin; i != iter->second.MZ_end; ++i)
1521 {
1522 elution_profile[index] += map[iter->second.RT_index][i].getIntensity();
1523 }
1524 elution_profile[index] /= iter->second.MZ_end - iter->second.MZ_begin + 1.;
1525 }
1526
1527 double max = 0;
1528 Int max_index = INT_MIN;
1529 for (Size i = 0; i < elution_profile.size(); ++i)
1530 {
1531 if (elution_profile[i] > max)
1532 {
1533 max_index = i;
1534 max = elution_profile[i];
1535 }
1536 }
1537
1538 Int max_extension = (Int)(elution_profile.size()) - 2 * max_index;
1539
1540 double av_elution = 0;
1541 for (Size i = 0; i < elution_profile.size(); ++i)
1542 {
1543 av_elution += elution_profile[i];
1544 }
1545 av_elution /= (double)elution_profile.size();
1546
1547 double sd_elution = 0;
1548 for (Size i = 0; i < elution_profile.size(); ++i)
1549 {
1550 sd_elution += (av_elution - elution_profile[i]) * (av_elution - elution_profile[i]);
1551 }
1552 sd_elution /= (double)(elution_profile.size() - 1);
1553 sd_elution = sqrt(sd_elution);
1554
1555 //Determine average m/z monoisotopic pos
1556 double av_mz = 0;
1557 for (iter = box.begin(); iter != box.end(); ++iter, ++index)
1558 {
1559 av_mz += iter->second.mz;
1560 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1561 std::cout << iter->second.RT << "\t" << iter->second.mz << "\t" << iter->second.c + 1 << std::endl;
1562 #endif
1563 }
1564 av_mz /= (double)box.size();
1565
1566
1567 //Boundary check
1568 if ((Int)(box.begin()->second.RT_index) - 1 < 0)
1569 {
1570 return;
1571 }
1572
1573 UInt pre_index = box.begin()->second.RT_index - 1;
1574 typename MSSpectrum::const_iterator c_iter = map[pre_index].MZBegin(av_mz);
1575 double pre_elution = 0;
1576
1577 double mz_start = map[pre_index + 1][box.begin()->second.MZ_begin].getMZ();
1578 double mz_end = map[pre_index + 1][box.begin()->second.MZ_end].getMZ();
1579
1580 typename MSSpectrum::const_iterator mz_start_iter = map[pre_index].MZBegin(mz_start), mz_end_iter = map[pre_index].MZBegin(mz_end);
1581 for (typename MSSpectrum::const_iterator mz_iter = mz_start_iter; mz_iter != mz_end_iter; ++mz_iter)
1582 {
1583 pre_elution += mz_iter->getIntensity();
1584 }
1585
1586
1587 //Do we need to extend at all?
1588 if (pre_elution <= av_elution - 2 * sd_elution)
1589 {
1590 return;
1591 }
1592
1593 Int c_index = max_extension;
1594 Int first_index = box.begin()->second.RT_index;
1595 for (Int i = 1; i < max_extension; ++i)
1596 {
1597 c_index = first_index - i;
1598 if (c_index < 0)
1599 {
1600 break;
1601 }
1602
1603 //CHECK Majority vote for charge???????????????
1604 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1605 std::cout << box.begin()->second.RT << "\t" << av_mz << "\t" << box.begin()->second.c + 1 << "\t" << " extending the box " << std::endl;
1606 #endif
1607
1608 push2Box_(av_mz, c_index, box.begin()->second.c, box.begin()->second.score, c_iter->getIntensity(),
1609 map[c_index].getRT(), box.begin()->second.MZ_begin, box.begin()->second.MZ_end);
1610 }
1611 }
1612
1613 template <typename PeakType>
1614 void IsotopeWaveletTransform<PeakType>::clusterSeeds_(const TransSpectrum& candidates,
1615 const MSSpectrum& ref, const UInt scan_index, const UInt c, const bool check_PPMs)
1616 {
1617 typename std::multimap<double, Box>::iterator iter;
1618 typename Box::iterator box_iter;
1619 std::vector<BoxElement> final_box;
1620 double c_mz, av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
1621 double virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
1622
1623 typename std::pair<double, double> c_extend;
1624 for (iter = tmp_boxes_->at(c).begin(); iter != tmp_boxes_->at(c).end(); ++iter)
1625 {
1626 Box& c_box = iter->second;
1627 av_score = 0; av_mz = 0; av_intens = 0; av_abs_intens = 0; count = 0;
1628 virtual_av_mz = 0; virtual_av_intens = 0; virtual_av_abs_intens = 0; virtual_count = 0;
1629
1630 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1631 {
1632 if (box_iter->second.score == 0) //virtual helping point
1633 {
1634 if (count != 0)
1635 continue; //it is in any way not pure virtual
1636
1637 c_mz = box_iter->second.mz;
1638 virtual_av_intens += box_iter->second.intens;
1639 virtual_av_abs_intens += fabs(box_iter->second.intens);
1640 virtual_av_mz += c_mz * fabs(box_iter->second.intens);
1641 ++virtual_count;
1642 }
1643 else
1644 {
1645 c_mz = box_iter->second.mz;
1646 av_score += box_iter->second.score;
1647 av_intens += box_iter->second.intens;
1648 av_abs_intens += fabs(box_iter->second.intens);
1649 av_mz += c_mz * fabs(box_iter->second.intens);
1650 ++count;
1651 }
1652 }
1653
1654 if (count == 0) //pure virtual helping box
1655 {
1656 av_intens = virtual_av_intens / virtual_count;
1657 av_score = 0;
1658 av_mz = virtual_av_mz / virtual_av_abs_intens;
1659 }
1660 else
1661 {
1662 av_intens /= count;
1663 av_score /= count;
1664 av_mz /= av_abs_intens;
1665 }
1666
1667 BoxElement c_box_element;
1668 c_box_element.mz = av_mz;
1669 c_box_element.c = c;
1670 c_box_element.score = av_score;
1671 c_box_element.intens = av_intens;
1672
1673 c_box_element.RT = c_box.begin()->second.RT;
1674
1675 final_box.push_back(c_box_element);
1676 }
1677
1678 UInt num_o_feature = final_box.size();
1679 if (num_o_feature == 0)
1680 {
1681 tmp_boxes_->at(c).clear();
1682 return;
1683 }
1684
1685 //Computing the derivatives
1686 std::vector<double> bwd_diffs(num_o_feature, 0);
1687
1688 bwd_diffs[0] = 0;
1689 for (UInt i = 1; i < num_o_feature; ++i)
1690 {
1691 bwd_diffs[i] = (final_box[i].intens - final_box[i - 1].intens) / (final_box[i].mz - final_box[i - 1].mz);
1692 }
1693
1694 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1695 std::ofstream ofile_bwd("bwd_gpu.dat");
1696 for (UInt i = 0; i < num_o_feature; ++i)
1697 {
1698 ofile_bwd << final_box[i].mz << "\t" << bwd_diffs[i] << std::endl;
1699 }
1700 ofile_bwd.close();
1701 #endif
1702
1703 for (UInt i = 0; i < num_o_feature - 1; ++i)
1704 {
1705 while (i < num_o_feature - 2)
1706 {
1707 if (final_box[i].score > 0 || final_box[i].score == -1000) //this has been an helping point
1708 break;
1709 ++i;
1710 }
1711
1712 if (bwd_diffs[i] > 0 && bwd_diffs[i + 1] < 0)
1713 {
1714 checkPositionForPlausibility_(candidates, ref, final_box[i].mz, final_box[i].c, scan_index, check_PPMs, final_box[i].intens, final_box[i].score);
1715 continue;
1716 }
1717 }
1718
1719 tmp_boxes_->at(c).clear();
1720 }
1721
1722 template <typename PeakType>
1723 FeatureMap IsotopeWaveletTransform<PeakType>::mapSeeds2Features(const PeakMap& map, const UInt RT_votes_cutoff)
1724 {
1725 FeatureMap feature_map;
1726 typename std::multimap<double, Box>::iterator iter;
1727 typename Box::iterator box_iter;
1728 UInt best_charge_index; double best_charge_score, c_mz, c_RT; UInt c_charge;
1729 double av_intens = 0, av_ref_intens = 0, av_score = 0, av_mz = 0, av_RT = 0, mz_cutoff, sum_of_ref_intenses_g;
1730 bool restart = false;
1731
1732 typename std::pair<double, double> c_extend;
1733 for (iter = closed_boxes_.begin(); iter != closed_boxes_.end(); ++iter)
1734 {
1735 sum_of_ref_intenses_g = 0;
1736 Box& c_box = iter->second;
1737 std::vector<double> charge_votes(max_charge_, 0), charge_binary_votes(max_charge_, 0);
1738 restart = false;
1739
1740 //Let's first determine the charge
1741 //Therefor, we can use two types of votes: qualitative ones (charge_binary_votes) or quantitative ones (charge_votes)
1742 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1743 {
1744 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1745 if (trunc(box_iter->second.mz) == 874)
1746 std::cout << box_iter->second.c << "\t" << box_iter->second.intens << "\t" << box_iter->second.score << std::endl;
1747 #endif
1748
1749 if (box_iter->second.score == -1000)
1750 {
1751 restart = true;
1752 break;
1753 }
1754
1755 charge_votes[box_iter->second.c] += box_iter->second.intens; //score; Do not use score, can get problematic for charge state 2 vs 4
1756 ++charge_binary_votes[box_iter->second.c];
1757 }
1758
1759 if (restart)
1760 {
1761 continue;
1762 }
1763
1764 //... determining the best fitting charge
1765 best_charge_index = 0; best_charge_score = 0;
1766 for (UInt i = 0; i < max_charge_; ++i)
1767 {
1768 if (charge_votes[i] > best_charge_score)
1769 {
1770 best_charge_index = i;
1771 best_charge_score = charge_votes[i];
1772 }
1773 }
1774
1775 //Pattern found in too few RT scan
1776 if (charge_binary_votes[best_charge_index] < RT_votes_cutoff && RT_votes_cutoff <= map.size())
1777 {
1778 continue;
1779 }
1780
1781 c_charge = best_charge_index + 1; //that's the finally predicted charge state for the pattern
1782
1783 av_intens = 0; av_ref_intens = 0; av_score = 0; av_mz = 0; av_RT = 0;
1784 //Now, let's get the RT boundaries for the box
1785 std::vector<DPosition<2> > point_set;
1786 double sum_of_ref_intenses_l;
1787 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1788 {
1789 sum_of_ref_intenses_l = 0;
1790 c_mz = box_iter->second.mz;
1791 c_RT = box_iter->second.RT;
1792
1793 mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_mz, c_charge);
1794
1795 point_set.push_back(DPosition<2>(c_RT, c_mz - Constants::IW_QUARTER_NEUTRON_MASS / (double)c_charge));
1796 //-1 since we are already at the first peak and +0.75, since this includes the last peak of the wavelet as a whole
1797 point_set.push_back(DPosition<2>(c_RT, c_mz + mz_cutoff / (double)c_charge));
1798
1799 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1800 std::cout << "Intenstype: " << intenstype_ << std::endl;
1801 #endif
1802 if (intenstype_ == "ref")
1803 {
1804 //Find monoisotopic max
1805 const MSSpectrum& c_spec(map[box_iter->second.RT_index]);
1806 //'Correct' possible shift
1807 for (unsigned int i = 0; i < mz_cutoff; ++i)
1808 {
1809 typename MSSpectrum::const_iterator h_iter = c_spec.MZBegin(c_mz + i * Constants::IW_NEUTRON_MASS / c_charge + Constants::IW_QUARTER_NEUTRON_MASS / (double)c_charge), hc_iter = c_spec.MZBegin(c_mz + i * Constants::IW_NEUTRON_MASS / c_charge);
1810
1811 hc_iter = c_spec.MZBegin(c_mz + i * Constants::IW_NEUTRON_MASS / c_charge);
1812
1813 while (h_iter != c_spec.begin())
1814 {
1815
1816 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1817 if (trunc(c_mz) == 874)
1818 {
1819 std::cout << "cmz: " << c_mz + i * Constants::IW_NEUTRON_MASS / c_charge << "\t" << hc_iter->getMZ() << "\t" << hc_iter->getIntensity() << "\t" << h_iter->getMZ() << "\t" << h_iter->getIntensity() << std::endl;
1820 }
1821 #endif
1822
1823 --h_iter;
1824 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
1825 {
1826 hc_iter = h_iter;
1827 }
1828
1829 if (c_mz + i * Constants::IW_NEUTRON_MASS / c_charge - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (double)c_charge)
1830 {
1831 break;
1832 }
1833 }
1834 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1835 if (trunc(c_mz) == 874)
1836 {
1837 std::cout << "c_mz: " << c_mz + i * Constants::IW_NEUTRON_MASS / c_charge << "\t" << hc_iter->getMZ() << "\t" << hc_iter->getIntensity() << "\t" << i * Constants::IW_NEUTRON_MASS / c_charge << "\t";
1838 }
1839 #endif
1840 sum_of_ref_intenses_l += hc_iter->getIntensity();
1841 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1842 if (trunc(c_mz) == 874)
1843 {
1844 std::cout << sum_of_ref_intenses_l << "********" << std::endl;
1845 }
1846 #endif
1847 }
1848 }
1849
1850 if (best_charge_index == box_iter->second.c)
1851 {
1852 av_score += box_iter->second.score;
1853 av_intens += box_iter->second.intens;
1854 av_ref_intens += box_iter->second.ref_intens;
1855 sum_of_ref_intenses_g += sum_of_ref_intenses_l;
1856 av_mz += c_mz * box_iter->second.intens;
1857 }
1858 av_RT += c_RT;
1859 }
1860
1861 av_mz /= av_intens;
1862 av_ref_intens /= (double)charge_binary_votes[best_charge_index];
1863 av_score /= (double)charge_binary_votes[best_charge_index];
1864 av_RT /= (double)c_box.size();
1865
1866 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1867 if (trunc(av_mz) == 874)
1868 std::cout << av_mz << "\t" << best_charge_index << "\t" << best_charge_score << std::endl;
1869 #endif
1870
1871 Feature c_feature;
1872 ConvexHull2D c_conv_hull;
1873 c_conv_hull.addPoints(point_set);
1874 c_feature.setCharge(c_charge);
1875 c_feature.setConvexHulls(std::vector<ConvexHull2D>(1, c_conv_hull));
1876
1877 //This makes the intensity value independent of the m/z (the lambda) value (Skellam distribution)
1878 if (intenstype_ == "corrected")
1879 {
1880 double lambda = IsotopeWavelet::getLambdaL(av_mz * c_charge);
1881 av_intens /= exp(-2 * lambda) * boost::math::cyl_bessel_i(0, 2 * lambda);
1882 }
1883 if (intenstype_ == "ref")
1884 {
1885 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1886 if (trunc(c_mz) == 874)
1887 {
1888 std::cout << sum_of_ref_intenses_g << "####" << std::endl;
1889 }
1890 #endif
1891
1892 av_intens = sum_of_ref_intenses_g;
1893 }
1894
1895 c_feature.setMZ(av_mz);
1896 c_feature.setIntensity(av_intens);
1897 c_feature.setRT(av_RT);
1898 c_feature.setOverallQuality(av_score);
1899 feature_map.push_back(c_feature);
1900 }
1901
1902 return feature_map;
1903 }
1904
1905 template <typename PeakType>
1906 bool IsotopeWaveletTransform<PeakType>::checkPositionForPlausibility_(const MSSpectrum& candidate,
1907 const MSSpectrum& ref, const double seed_mz, const UInt c, const UInt scan_index, const bool check_PPMs, const double transintens, const double prev_score)
1908 {
1909 typename MSSpectrum::const_iterator iter, ref_iter;
1910 UInt peak_cutoff;
1911 peak_cutoff = IsotopeWavelet::getNumPeakCutOff(seed_mz, c + 1);
1912
1913 iter = candidate.MZBegin(seed_mz);
1914 //we can ignore those cases
1915 if (iter == candidate.begin() || iter == candidate.end())
1916 {
1917 return false;
1918 }
1919
1920 std::pair<double, double> reals;
1921 ref_iter = ref.MZBegin(seed_mz);
1922 //Correct the position
1923 double real_mz, real_intens;
1924 if (check_PPMs)
1925 {
1926 reals = checkPPMTheoModel_(ref, iter->getMZ(), c);
1927 real_mz = reals.first; real_intens = reals.second;
1928 //if (real_mz <= 0 || real_intens <= 0)
1929 //{
1930 typename MSSpectrum::const_iterator h_iter = ref_iter, hc_iter = ref_iter;
1931 while (h_iter != ref.begin())
1932 {
1933 --h_iter;
1934 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
1935 {
1936 hc_iter = h_iter;
1937 }
1938 else
1939 {
1940 break;
1941 }
1942
1943 if (seed_mz - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.))
1944 {
1945 return false;
1946 }
1947 }
1948 reals = checkPPMTheoModel_(ref, h_iter->getMZ(), c);
1949 real_mz = reals.first; real_intens = reals.second;
1950 if (real_mz <= 0 || real_intens <= 0)
1951 {
1952 return false;
1953 }
1954 real_mz = h_iter->getMZ();
1955 real_intens = h_iter->getIntensity();
1956 //}
1957 }
1958 else
1959 {
1960 reals = std::pair<double, double>(seed_mz, ref_iter->getIntensity());
1961 real_mz = reals.first; real_intens = reals.second;
1962
1963 if (real_mz <= 0 || real_intens <= 0)
1964 {
1965 typename MSSpectrum::const_iterator h_iter = ref_iter, hc_iter = ref_iter;
1966 while (h_iter != ref.begin())
1967 {
1968 --h_iter;
1969 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
1970 {
1971 hc_iter = h_iter;
1972 }
1973 else
1974 {
1975 break;
1976 }
1977
1978 if (seed_mz - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.))
1979 {
1980 return false;
1981 }
1982 }
1983 real_mz = h_iter->getMZ(); real_intens = h_iter->getIntensity();
1984 if (real_mz <= 0 || real_intens <= 0)
1985 {
1986 return false;
1987 }
1988 real_mz = h_iter->getMZ();
1989 real_intens = h_iter->getIntensity();
1990 }
1991 }
1992
1993 double c_score = scoreThis_(candidate, peak_cutoff, real_mz, c, 0);
1994
1995 if (c_score <= 0)
1996 {
1997 return false;
1998 }
1999
2000 double mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(real_mz, c + 1);
2001 typename MSSpectrum::const_iterator real_l_MZ_iter = ref.MZBegin(real_mz - Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.));
2002 typename MSSpectrum::const_iterator real_r_MZ_iter = ref.MZBegin(real_l_MZ_iter, real_mz + mz_cutoff / (c + 1.), ref.end());
2003 if (real_r_MZ_iter == ref.end())
2004 {
2005 --real_r_MZ_iter;
2006 }
2007
2008
2009 UInt real_mz_begin = distance(ref.begin(), real_l_MZ_iter);
2010 UInt real_mz_end = distance(ref.begin(), real_r_MZ_iter);
2011
2012 if (prev_score == -1000)
2013 {
2014 push2Box_(real_mz, scan_index, c, prev_score, transintens, ref.getRT(), real_mz_begin, real_mz_end, real_intens);
2015 }
2016 else
2017 {
2018 push2Box_(real_mz, scan_index, c, c_score, transintens, ref.getRT(), real_mz_begin, real_mz_end, real_intens);
2019 }
2020 return true;
2021 }
2022
2023 template <typename PeakType>
2024 bool IsotopeWaveletTransform<PeakType>::checkPositionForPlausibility_(const TransSpectrum& candidate,
2025 const MSSpectrum& ref, const double seed_mz, const UInt c, const UInt scan_index, const bool check_PPMs, const double transintens, const double prev_score)
2026 {
2027 typename MSSpectrum::const_iterator iter, ref_iter;
2028 UInt peak_cutoff;
2029 peak_cutoff = IsotopeWavelet::getNumPeakCutOff(seed_mz, c + 1);
2030
2031 iter = candidate.MZBegin(seed_mz);
2032 //we can ignore those cases
2033 if (iter == candidate.begin() || iter == candidate.end())
2034 {
2035 return false;
2036 }
2037
2038 std::pair<double, double> reals;
2039 ref_iter = ref.MZBegin(seed_mz);
2040 //Correct the position
2041 double real_mz, real_intens;
2042 if (check_PPMs)
2043 {
2044 reals = checkPPMTheoModel_(ref, iter->getMZ(), c);
2045 real_mz = reals.first; real_intens = reals.second;
2046 //if (real_mz <= 0 || real_intens <= 0)
2047 //{
2048 auto h_iter = ref_iter, hc_iter = ref_iter;
2049 while (h_iter != ref.begin())
2050 {
2051 --h_iter;
2052 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2053 {
2054 hc_iter = h_iter;
2055 }
2056 else
2057 {
2058 break;
2059 }
2060
2061 if (seed_mz - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.))
2062 {
2063 return false;
2064 }
2065 }
2066 ++h_iter;
2067 reals = checkPPMTheoModel_(ref, h_iter->getMZ(), c);
2068 real_mz = reals.first; real_intens = reals.second;
2069
2070 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2071 std::cout << "Plausibility check old_mz: " << iter->getMZ() << "\t" << real_mz << std::endl;
2072 #endif
2073
2074 if (real_mz <= 0 || real_intens <= 0)
2075 {
2076 return false;
2077 }
2078 real_mz = h_iter->getMZ();
2079 real_intens = h_iter->getIntensity();
2080 //}
2081 }
2082 else
2083 {
2084 reals = std::pair<double, double>(seed_mz, ref_iter->getIntensity());
2085 real_mz = reals.first; real_intens = reals.second;
2086
2087 if (real_mz <= 0 || real_intens <= 0)
2088 {
2089 auto h_iter = ref_iter, hc_iter = ref_iter;
2090 while (h_iter != ref.begin())
2091 {
2092 --h_iter;
2093 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2094 {
2095 hc_iter = h_iter;
2096 }
2097 else
2098 {
2099 break;
2100 }
2101
2102 if (seed_mz - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.))
2103 {
2104 return false;
2105 }
2106 }
2107 real_mz = h_iter->getMZ(); real_intens = h_iter->getIntensity();
2108 if (real_mz <= 0 || real_intens <= 0)
2109 {
2110 return false;
2111 }
2112 real_mz = h_iter->getMZ();
2113 real_intens = h_iter->getIntensity();
2114 }
2115 }
2116
2117 double c_score = scoreThis_(candidate, peak_cutoff, real_mz, c, 0);
2118
2119 if (c_score <= 0)
2120 {
2121 return false;
2122 }
2123
2124 double mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(real_mz, c + 1);
2125 typename MSSpectrum::const_iterator real_l_MZ_iter = ref.MZBegin(real_mz - Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.));
2126 typename MSSpectrum::const_iterator real_r_MZ_iter = ref.MZBegin(real_l_MZ_iter, real_mz + mz_cutoff / (c + 1.), ref.end());
2127 if (real_r_MZ_iter == ref.end())
2128 {
2129 --real_r_MZ_iter;
2130 }
2131
2132
2133 UInt real_mz_begin = distance(ref.begin(), real_l_MZ_iter);
2134 UInt real_mz_end = distance(ref.begin(), real_r_MZ_iter);
2135
2136 if (prev_score == -1000)
2137 {
2138 push2Box_(real_mz, scan_index, c, prev_score, transintens, ref.getRT(), real_mz_begin, real_mz_end, real_intens);
2139 }
2140 else
2141 {
2142 push2Box_(real_mz, scan_index, c, c_score, transintens, ref.getRT(), real_mz_begin, real_mz_end, real_intens);
2143 }
2144
2145 return true;
2146 }
2147
2148 template <typename PeakType>
2149 std::pair<double, double> IsotopeWaveletTransform<PeakType>::checkPPMTheoModel_(const MSSpectrum& ref, const double c_mz, const UInt c)
2150 {
2151 double mass = c_mz * (c + 1) - Constants::IW_PROTON_MASS * (c);
2152 double ppms = getPPMs_(peptideMassRule_(mass), mass);
2153 if (ppms >= Constants::PEPTIDE_MASS_RULE_THEO_PPM_BOUND)
2154 {
2155 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2156 std::cout << ::std::setprecision(8) << std::fixed << c_mz << "\t =(" << "ISO_WAVE" << ")> " << "REJECT \t" << ppms << " (rule: " << peptideMassRule_(mass) << " got: " << mass << ")" << std::endl;
2157 #endif
2158 return std::pair<double, double>(-1, -1);
2159 }
2160
2161 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2162 std::cout << ::std::setprecision(8) << std::fixed << c_mz << "\t =(" << "ISO_WAVE" << ")> " << "ACCEPT \t" << ppms << " (rule: " << peptideMassRule_(mass) << " got: " << mass << ")" << std::endl;
2163 #endif
2164 return std::pair<double, double>(c_mz, ref.MZBegin(c_mz)->getIntensity());
2165 }
2166
2167 } //namespace
2168
2169 #pragma clang diagnostic pop
2170
2171
OpenMS::Feature
An LC-MS feature.
Definition: Feature.h:70
OpenMS::Peak2D::setRT
void setRT(CoordinateType coordinate)
Mutable access to the RT coordinate (index 0)
Definition: Peak2D.h:214
OpenMS::IsotopeWaveletTransform::peptideMassRule_
double peptideMassRule_(const double c_mass) const
Returns the monoisotopic mass (with corresponding decimal values) we would expect at c_mass...
Definition: IsotopeWaveletTransform.h:447
OpenMS::SignedSize
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:134
OpenMS::positionComparator
bool positionComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:513
OpenMS::Constants::DEFAULT_NUM_OF_INTERPOLATION_POINTS
const unsigned int DEFAULT_NUM_OF_INTERPOLATION_POINTS
Definition: IsotopeWaveletConstants.h:43
OpenMS::String
A more convenient string class.
Definition: String.h:58
OpenMS::IsotopeWaveletTransform::initializeScan
virtual void initializeScan(const MSSpectrum &c_ref, const UInt c=0)
Definition: IsotopeWaveletTransform.h:644
OpenMS::Feature::setOverallQuality
void setOverallQuality(QualityType q)
Set the overall quality.
OpenMS::IsotopeWaveletTransform::TransSpectrum::end
MSSpectrum::const_iterator end() const
Definition: IsotopeWaveletTransform.h:191
OpenMS::IsotopeWaveletTransform::c_mzs_
std::vector< double > c_mzs_
Definition: IsotopeWaveletTransform.h:480
LinearRegression.h
OpenMS::Peak2D::getIntensity
IntensityType getIntensity() const
Definition: Peak2D.h:166
OpenMS::IsotopeWaveletTransform::intenstype_
String intenstype_
Definition: IsotopeWaveletTransform.h:486
OpenMS::IsotopeWaveletTransform::hr_data_
bool hr_data_
Definition: IsotopeWaveletTransform.h:485
OpenMS::IsotopeWaveletTransform::getAvMZSpacing_
double getAvMZSpacing_(const MSSpectrum &scan)
Computes the average MZ spacing of scan.
Definition: IsotopeWaveletTransform.h:1180
OpenMS::Peak2D
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:54
OpenMS::IsotopeWavelet::init
static IsotopeWavelet * init(const double max_m, const UInt max_charge)
OpenMS::IsotopeWaveletTransform::getMaxScanSize
Size getMaxScanSize() const
Definition: IsotopeWaveletTransform.h:317
OpenMS::MSSpectrum::getRT
double getRT() const
OpenMS::IsotopeWaveletTransform::TransSpectrum::getRT
double getRT() const
Definition: IsotopeWaveletTransform.h:128
OpenMS::IsotopeWaveletTransform::BoxElement::mz
double mz
The monoisotopic position.
Definition: IsotopeWaveletTransform.h:77
OpenMS::IsotopeWaveletTransform::TransSpectrum::size
Size size() const
Definition: IsotopeWaveletTransform.h:158
OpenMS::IsotopeWaveletTransform::data_length_
UInt data_length_
Definition: IsotopeWaveletTransform.h:484
OpenMS::IsotopeWaveletTransform::max_num_peaks_per_pattern_
UInt max_num_peaks_per_pattern_
Definition: IsotopeWaveletTransform.h:484
OpenMS::Peak2D::setMZ
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:202
OpenMS::Constants::IW_PROTON_MASS
const double IW_PROTON_MASS
Definition: IsotopeWaveletConstants.h:68
OpenMS::IsotopeWaveletTransform::TransSpectrum::TransSpectrum
TransSpectrum()
Definition: IsotopeWaveletTransform.h:101
OpenMS::ConstRefVector::sortByIntensity
void sortByIntensity(bool reverse=false)
Sorting.
Definition: ConstRefVector.h:741
OpenMS::IsotopeWaveletTransform::push2TmpBox_
virtual void push2TmpBox_(const double mz, const UInt scan, UInt charge, const double score, const double intens, const double rt, const UInt MZ_begin, const UInt MZ_end)
Essentially the same function as.
Definition: IsotopeWaveletTransform.h:1334
OpenMS::Constants::PEPTIDE_MASS_RULE_BOUND
const double PEPTIDE_MASS_RULE_BOUND
Definition: IsotopeWaveletConstants.h:50
OpenMS::IsotopeWaveletTransform::max_mz_cutoff_
double max_mz_cutoff_
Definition: IsotopeWaveletTransform.h:490
OpenMS::MSExperiment
In-Memory representation of a mass spectrometry experiment.
Definition: MSExperiment.h:77
OpenMS::IsotopeWaveletTransform::Box
std::multimap< UInt, BoxElement > Box
Key: RT index, value: BoxElement.
Definition: IsotopeWaveletTransform.h:88
OpenMS::IsotopeWaveletTransform::from_max_to_right_
Int from_max_to_right_
Definition: IsotopeWaveletTransform.h:487
OpenMS::IsotopeWaveletTransform::push2Box_
virtual void push2Box_(const double mz, const UInt scan, UInt c, const double score, const double intens, const double rt, const UInt MZ_begin, const UInt MZ_end, const double ref_intens)
Inserts a potential isotopic pattern into an open box or - if no such box exists - creates a new one...
Definition: IsotopeWaveletTransform.h:1228
OpenMS::IsotopeWaveletTransform::max_charge_
UInt max_charge_
Definition: IsotopeWaveletTransform.h:484
Exception.h
OpenMS::FeatureMap
A container for features.
Definition: FeatureMap.h:95
OpenMS::IsotopeWaveletTransform::updateBoxStates
void updateBoxStates(const PeakMap &map, const Size scan_index, const UInt RT_interleave, const UInt RT_votes_cutoff, const Int front_bound=-1, const Int end_bound=-1)
A function keeping track of currently open and closed sweep line boxes. This function is used by the ...
Definition: IsotopeWaveletTransform.h:1439
OpenMS::IsotopeWaveletTransform::scores_
std::vector< float > scores_
Definition: IsotopeWaveletTransform.h:491
OpenMS::Constants::IW_NEUTRON_MASS
const double IW_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:54
OpenMS::IsotopeWaveletTransform::getClosedBoxes
virtual std::multimap< double, Box > getClosedBoxes()
Returns the closed boxes.
Definition: IsotopeWaveletTransform.h:276
OpenMS::IsotopeWaveletTransform::getLinearInterpolation
double getLinearInterpolation(const double mz_a, const double intens_a, const double mz_pos, const double mz_b, const double intens_b)
Computes a linear (intensity) interpolation.
Definition: IsotopeWaveletTransform.h:295
OpenMS::IsotopeWaveletTransform::TransSpectrum::getTransIntensity
double getTransIntensity(const UInt i) const
Definition: IsotopeWaveletTransform.h:146
StandardTypes.h
OpenMS::IsotopeWaveletTransform::from_max_to_left_
Int from_max_to_left_
Definition: IsotopeWaveletTransform.h:487
OpenMS::Constants::IW_QUARTER_NEUTRON_MASS
const double IW_QUARTER_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:56
OpenMS::IsotopeWaveletTransform::av_MZ_spacing_
double av_MZ_spacing_
Definition: IsotopeWaveletTransform.h:479
OpenMS::IsotopeWaveletTransform::BoxElement
Internally used data structure.
Definition: IsotopeWaveletTransform.h:75
OpenMS::IsotopeWaveletTransform::BoxElement::intens
double intens
The transformed intensity at the monoisotopic mass.
Definition: IsotopeWaveletTransform.h:80
OpenMS::IsotopeWaveletTransform
A class implementing the isotope wavelet transform. If you just want to find features using the isoto...
Definition: IsotopeWaveletTransform.h:69
OpenMS::IsotopeWaveletTransform::TransSpectrum::getRefSpectrum
const MSSpectrum * getRefSpectrum() const
Definition: IsotopeWaveletTransform.h:170
OpenMS::Math::mean
static double mean(IteratorType begin, IteratorType end)
Calculates the mean of a range of values.
Definition: StatisticFunctions.h:133
OpenMS::MSSpectrum::MZBegin
Iterator MZBegin(CoordinateType mz)
Binary search for peak range begin.
OpenMS::IsotopeWaveletTransform::sampleTheCMarrWavelet_
void sampleTheCMarrWavelet_(const MSSpectrum &scan, const Int wavelet_length, const Int mz_index, const UInt charge)
OpenMS::intensityAscendingComparator
bool intensityAscendingComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:501
IsotopeWaveletConstants.h
OpenMS::IsotopeWaveletTransform::prod_
std::vector< double > prod_
Definition: IsotopeWaveletTransform.h:480
OpenMS::ConvexHull2D::addPoints
void addPoints(const PointArrayType &points)
OpenMS::IsotopeWaveletTransform::TransSpectrum::MZBegin
MSSpectrum::const_iterator MZBegin(const double mz) const
Definition: IsotopeWaveletTransform.h:177
OpenMS::Feature::setConvexHulls
void setConvexHulls(const std::vector< ConvexHull2D > &hulls)
Set the convex hulls of single mass traces.
OpenMS::IsotopeWaveletTransform::TransSpectrum::getRefSpectrum
const MSSpectrum * getRefSpectrum()
Definition: IsotopeWaveletTransform.h:164
OpenMS::Int
int Int
Signed integer type.
Definition: Types.h:102
IsotopeWavelet.h
OpenMS::IsotopeWaveletTransform::TransSpectrum::MZEnd
MSSpectrum::const_iterator MZEnd(const double mz) const
Definition: IsotopeWaveletTransform.h:184
OpenMS::IsotopeWaveletTransform::TransSpectrum
Internally (only by GPUs) used data structure . It allows efficient data exchange between CPU and GPU...
Definition: IsotopeWaveletTransform.h:94
OpenMS::MSSpectrum::MZEnd
Iterator MZEnd(CoordinateType mz)
Binary search for peak range end (returns the past-the-end iterator)
OpenMS::MSSpectrum
The representation of a 1D spectrum.
Definition: MSSpectrum.h:67
OpenMS::intensityComparator
bool intensityComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:495
OpenMS::IsotopeWaveletTransform::BoxElement::c
UInt c
Note, this is not the charge (it is charge-1!!!)
Definition: IsotopeWaveletTransform.h:78
FeatureMap.h
OpenMS::IsotopeWaveletTransform::getTransform
virtual void getTransform(MSSpectrum &c_trans, const MSSpectrum &c_ref, const UInt c)
Computes the isotope wavelet transform of charge state c.
Definition: IsotopeWaveletTransform.h:563
OpenMS::IsotopeWavelet::getLambdaL
static double getLambdaL(const double m)
Returns the mass-parameter lambda (linear fit).
OpenMS::BaseFeature::setCharge
void setCharge(const ChargeType &ch)
Set charge state.
OpenMS::IsotopeWavelet::getMzPeakCutOffAtMonoPos
static UInt getMzPeakCutOffAtMonoPos(const double mass, const UInt z)
OpenMS::IsotopeWaveletTransform::max_scan_size_
Size max_scan_size_
Definition: IsotopeWaveletTransform.h:483
OpenMS::IsotopeWaveletTransform::TransSpectrum::trans_intens_
std::vector< float > * trans_intens_
The intensities of the transform.
Definition: IsotopeWaveletTransform.h:206
OpenMS::MSExperiment::clear
void clear(bool clear_meta_data)
Clears all data and meta data.
OpenMS::Constants::PEPTIDE_MASS_RULE_THEO_PPM_BOUND
const double PEPTIDE_MASS_RULE_THEO_PPM_BOUND
Definition: IsotopeWaveletConstants.h:51
OpenMS::IsotopeWaveletTransform::BoxElement::RT_index
UInt RT_index
The elution time (map) index.
Definition: IsotopeWaveletTransform.h:83
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
OpenMS::IsotopeWaveletTransform::interpol_xs_
std::vector< double > interpol_xs_
Definition: IsotopeWaveletTransform.h:481
OpenMS::IsotopeWaveletTransform::getPPMs_
double getPPMs_(const double mass_a, const double mass_b) const
Returns the parts-per-million deviation of the masses.
Definition: IsotopeWaveletTransform.h:470
OpenMS::IsotopeWaveletTransform::clusterSeeds_
void clusterSeeds_(const TransSpectrum &candidates, const MSSpectrum &ref, const UInt scan_index, const UInt c, const bool check_PPMs)
Clusters the seeds stored by push2TmpBox_.
Definition: IsotopeWaveletTransform.h:1614
OpenMS::IsotopeWaveletTransform::interpol_ys_
std::vector< double > interpol_ys_
Definition: IsotopeWaveletTransform.h:481
MSExperiment.h
OpenMS::IsotopeWaveletTransform::IsotopeWaveletTransform
IsotopeWaveletTransform()
Default Constructor.
Definition: IsotopeWaveletTransform.h:519
OpenMS::IsotopeWaveletTransform::BoxElement::MZ_end
UInt MZ_end
Index.
Definition: IsotopeWaveletTransform.h:85
OpenMS::IsotopeWaveletTransform::TransSpectrum::setTransIntensity
void setTransIntensity(const UInt i, const double intens)
Definition: IsotopeWaveletTransform.h:152
OpenMS::ConstRefVector
This vector holds pointer to the elements of another container.
Definition: ConstRefVector.h:70
OpenMS::IsotopeWaveletTransform::xs_
std::vector< double > xs_
Definition: IsotopeWaveletTransform.h:480
OpenMS::IsotopeWaveletTransform::extendBox_
void extendBox_(const PeakMap &map, const Box &box)
A currently still necessary function that extends the box box in order to capture also signals whose ...
Definition: IsotopeWaveletTransform.h:1508
OpenMS::IsotopeWaveletTransform::TransSpectrum::getMZ
double getMZ(const UInt i) const
Definition: IsotopeWaveletTransform.h:134
OpenMS::IsotopeWavelet::getValueByLambda
static double getValueByLambda(const double lambda, const double tz1)
Returns the value of the isotope wavelet at position t via a fast table lookup. Usually, you do not need to call this function. Please use.
MSSpectrum.h
OpenMS::UInt
unsigned int UInt
Unsigned integer type.
Definition: Types.h:94
OpenMS::Constants::IW_HALF_NEUTRON_MASS
const double IW_HALF_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:55
OpenMS::IsotopeWaveletTransform::BoxElement::RT
double RT
The elution time (not the scan index)
Definition: IsotopeWaveletTransform.h:82
OpenMS::IsotopeWaveletTransform::c_spacings_
std::vector< double > c_spacings_
Definition: IsotopeWaveletTransform.h:480
OpenMS::DPosition< 2 >
OpenMS::IsotopeWaveletTransform::BoxElement::score
double score
The associated score.
Definition: IsotopeWaveletTransform.h:79
OpenMS::IsotopeWaveletTransform::zeros_
std::vector< float > zeros_
Definition: IsotopeWaveletTransform.h:491
OpenMS::Peak2D::getMZ
CoordinateType getMZ() const
Returns the m/z coordinate (index 1)
Definition: Peak2D.h:196
double
OpenMS::Constants::c
const double c
OpenMS::MSExperiment::size
Size size() const
Definition: MSExperiment.h:127
OpenMS::ConstRefVector::ConstRefVectorIterator
Mutable iterator for the ConstRefVector.
Definition: ConstRefVector.h:236
OpenMS::IsotopeWaveletTransform::sigma_
double sigma_
Definition: IsotopeWaveletTransform.h:479
OpenMS::IsotopeWaveletTransform::getTransformHighRes
virtual void getTransformHighRes(MSSpectrum &c_trans, const MSSpectrum &c_ref, const UInt c)
Computes the isotope wavelet transform of charge state c.
Definition: IsotopeWaveletTransform.h:606
int
OpenMS::MSExperiment::begin
Iterator begin()
Definition: MSExperiment.h:157
OpenMS::IsotopeWaveletTransform::BoxElement::MZ_begin
UInt MZ_begin
Index.
Definition: IsotopeWaveletTransform.h:84
OpenMS::IsotopeWaveletTransform::getSigma
double getSigma() const
Definition: IsotopeWaveletTransform.h:300
OpenMS::IsotopeWaveletTransform::front_boxes_
std::multimap< double, Box > front_boxes_
Definition: IsotopeWaveletTransform.h:476
OpenMS::IsotopeWaveletTransform::TransSpectrum::destroy
virtual void destroy()
Definition: IsotopeWaveletTransform.h:119
OpenMS::IsotopeWaveletTransform::~IsotopeWaveletTransform
virtual ~IsotopeWaveletTransform()
Destructor.
Definition: IsotopeWaveletTransform.h:557
OpenMS::IsotopeWaveletTransform::TransSpectrum::~TransSpectrum
virtual ~TransSpectrum()
Definition: IsotopeWaveletTransform.h:114
OpenMS::IsotopeWaveletTransform::open_boxes_
std::multimap< double, Box > open_boxes_
Definition: IsotopeWaveletTransform.h:476
OpenMS::IsotopeWaveletTransform::getSdIntens_
double getSdIntens_(const TransSpectrum &scan, const double mean)
Computes the standard deviation (neglecting negative values) of the (transformed) intensities of scan...
Definition: IsotopeWaveletTransform.h:1213
OpenMS::Constants::PEPTIDE_MASS_RULE_FACTOR
const double PEPTIDE_MASS_RULE_FACTOR
Definition: IsotopeWaveletConstants.h:49
OpenMS::IsotopeWaveletTransform::computeMinSpacing
virtual void computeMinSpacing(const MSSpectrum &c_ref)
Definition: IsotopeWaveletTransform.h:686
OpenMS::IsotopeWaveletTransform::checkPPMTheoModel_
virtual std::pair< double, double > checkPPMTheoModel_(const MSSpectrum &ref, const double c_mz, const UInt c)
Definition: IsotopeWaveletTransform.h:2149
OpenMS::IsotopeWaveletTransform::checkPositionForPlausibility_
virtual bool checkPositionForPlausibility_(const TransSpectrum &candidate, const MSSpectrum &ref, const double seed_mz, const UInt c, const UInt scan_index, const bool check_PPMs, const double transintens, const double prev_score)
A ugly but necessary function to handle "off-by-1-Dalton predictions" due to idiosyncrasies of the da...
Definition: IsotopeWaveletTransform.h:2024
OpenMS::IsotopeWaveletTransform::TransSpectrum::TransSpectrum
TransSpectrum(const MSSpectrum *reference)
Definition: IsotopeWaveletTransform.h:107
OpenMS::IsotopeWaveletTransform::TransSpectrum::begin
MSSpectrum::const_iterator begin() const
Definition: IsotopeWaveletTransform.h:198
OpenMS::IsotopeWaveletTransform::end_boxes_
std::multimap< double, Box > end_boxes_
Definition: IsotopeWaveletTransform.h:476
OpenMS::IsotopeWaveletTransform::getAvIntens_
double getAvIntens_(const TransSpectrum &scan)
Computes the average (transformed) intensity (neglecting negative values) of scan.
Definition: IsotopeWaveletTransform.h:1199
OpenMS::IsotopeWaveletTransform::TransSpectrum::getRefIntensity
double getRefIntensity(const UInt i) const
Definition: IsotopeWaveletTransform.h:140
OpenMS::IsotopeWaveletTransform::TransSpectrum::reference_
const MSSpectrum * reference_
The reference spectrum.
Definition: IsotopeWaveletTransform.h:205
OpenMS::intensityPointerComparator
bool intensityPointerComparator(PeakType *a, PeakType *b)
Definition: IsotopeWaveletTransform.h:507
OpenMS::IsotopeWaveletTransform::mapSeeds2Features
FeatureMap mapSeeds2Features(const PeakMap &map, const UInt RT_votes_cutoff)
Filters the candidates further more and maps the internally used data structures to the OpenMS framew...
Definition: IsotopeWaveletTransform.h:1723
OpenMS::IsotopeWaveletTransform::min_spacing_
double min_spacing_
Definition: IsotopeWaveletTransform.h:490
OpenMS::IsotopeWaveletTransform::closed_boxes_
std::multimap< double, Box > closed_boxes_
Definition: IsotopeWaveletTransform.h:476
OpenMS::IsotopeWavelet::getNumPeakCutOff
static UInt getNumPeakCutOff(const double mass, const UInt z)
OpenMS::IsotopeWaveletTransform::psi_
std::vector< double > psi_
Definition: IsotopeWaveletTransform.h:480
OpenMS::Peak2D::setIntensity
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.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::IsotopeWaveletTransform::getMinSpacing
double getMinSpacing() const
Definition: IsotopeWaveletTransform.h:312
OpenMS::IsotopeWaveletTransform::setSigma
void setSigma(const double sigma)
Definition: IsotopeWaveletTransform.h:305
OpenMS::ConvexHull2D
A 2-dimensional hull representation in [counter]clockwise direction - depending on axis labelling...
Definition: ConvexHull2D.h:72
OpenMS::IsotopeWaveletTransform::tmp_boxes_
std::vector< std::multimap< double, Box > > * tmp_boxes_
Definition: IsotopeWaveletTransform.h:477
ConstRefVector.h
OpenMS::IsotopeWaveletTransform::scoreThis_
virtual double scoreThis_(const TransSpectrum &candidate, const UInt peak_cutoff, const double seed_mz, const UInt c, const double ampl_cutoff)
Given a candidate for an isotopic pattern, this function computes the corresponding score...
Definition: IsotopeWaveletTransform.h:967
OpenMS::IsotopeWaveletTransform::identifyCharge
virtual void identifyCharge(const MSSpectrum &candidates, const MSSpectrum &ref, const UInt scan_index, const UInt c, const double ampl_cutoff, const bool check_PPMs)
Given an isotope wavelet transformed spectrum candidates, this function assigns to every significant ...
Definition: IsotopeWaveletTransform.h:696
OpenMS::IsotopeWaveletTransform::BoxElement::ref_intens
double ref_intens
Definition: IsotopeWaveletTransform.h:81
OpenMS::IsotopeWaveletTransform::indices_
std::vector< int > indices_
Definition: IsotopeWaveletTransform.h:488
OpenMS::IsotopeWaveletTransform::getLinearInterpolation
double getLinearInterpolation(const typename MSSpectrum::const_iterator &left_iter, const double mz_pos, const typename MSSpectrum::const_iterator &right_iter)
Computes a linear (intensity) interpolation.
Definition: IsotopeWaveletTransform.h:284