OpenMS  2.6.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 
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>
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;
86  };
87 
88  typedef std::multimap<UInt, BoxElement> Box;
89 
90 
95  {
97 
98 public:
99 
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 
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 
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 
480  std::vector<double> c_mzs_, c_spacings_, psi_, prod_, xs_;
481  std::vector<double> interpol_xs_, interpol_ys_;
482 
485  bool hr_data_;
488  std::vector<int> indices_;
489 
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>
502  {
503  return a.getIntensity() < b.getIntensity();
504  }
505 
506  template <typename PeakType>
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>
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);
554  }
555 
556  template <typename PeakType>
558  {
559  delete (tmp_boxes_);
560  }
561 
562  template <typename PeakType>
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>
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>
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>
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>
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());
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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);
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::IsotopeWaveletTransform::tmp_boxes_
std::vector< std::multimap< double, Box > > * tmp_boxes_
Definition: IsotopeWaveletTransform.h:477
OpenMS::Constants::IW_HALF_NEUTRON_MASS
const double IW_HALF_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:55
OpenMS::Constants::DEFAULT_NUM_OF_INTERPOLATION_POINTS
const unsigned int DEFAULT_NUM_OF_INTERPOLATION_POINTS
Definition: IsotopeWaveletConstants.h:43
OpenMS::IsotopeWaveletTransform::BoxElement::ref_intens
double ref_intens
Definition: IsotopeWaveletTransform.h:81
OpenMS::IsotopeWaveletTransform::TransSpectrum::getMZ
double getMZ(const UInt i) const
Definition: IsotopeWaveletTransform.h:134
double
OpenMS::IsotopeWaveletTransform::TransSpectrum::end
MSSpectrum::const_iterator end() const
Definition: IsotopeWaveletTransform.h:191
OpenMS::ConvexHull2D::addPoints
void addPoints(const PointArrayType &points)
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
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::intensityAscendingComparator
bool intensityAscendingComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:501
OpenMS::IsotopeWaveletTransform::av_MZ_spacing_
double av_MZ_spacing_
Definition: IsotopeWaveletTransform.h:479
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::Peak2D
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:54
OpenMS::IsotopeWaveletTransform::TransSpectrum::reference_
const MSSpectrum * reference_
The reference spectrum.
Definition: IsotopeWaveletTransform.h:205
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::String
A more convenient string class.
Definition: String.h:59
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
OpenMS::MSExperiment::begin
Iterator begin()
Definition: MSExperiment.h:157
IsotopeWavelet.h
OpenMS::IsotopeWaveletTransform::getAvIntens_
double getAvIntens_(const TransSpectrum &scan)
Computes the average (transformed) intensity (neglecting negative values) of scan.
Definition: IsotopeWaveletTransform.h:1199
OpenMS::MSExperiment
In-Memory representation of a mass spectrometry experiment.
Definition: MSExperiment.h:77
OpenMS::Size
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
OpenMS::Constants::c
const double c
OpenMS::MSSpectrum::getRT
double getRT() const
OpenMS::ConstRefVector::ConstRefVectorIterator
Mutable iterator for the ConstRefVector.
Definition: ConstRefVector.h:236
OpenMS::IsotopeWaveletTransform::getSigma
double getSigma() const
Definition: IsotopeWaveletTransform.h:300
OpenMS::IsotopeWaveletTransform::interpol_xs_
std::vector< double > interpol_xs_
Definition: IsotopeWaveletTransform.h:481
OpenMS::IsotopeWaveletTransform::c_spacings_
std::vector< double > c_spacings_
Definition: IsotopeWaveletTransform.h:480
OpenMS::Peak2D::getIntensity
IntensityType getIntensity() const
Definition: Peak2D.h:166
OpenMS::IsotopeWaveletTransform::c_mzs_
std::vector< double > c_mzs_
Definition: IsotopeWaveletTransform.h:480
OpenMS::IsotopeWaveletTransform::front_boxes_
std::multimap< double, Box > front_boxes_
Definition: IsotopeWaveletTransform.h:476
OpenMS::IsotopeWaveletTransform::max_mz_cutoff_
double max_mz_cutoff_
Definition: IsotopeWaveletTransform.h:490
OpenMS::IsotopeWaveletTransform::TransSpectrum::~TransSpectrum
virtual ~TransSpectrum()
Definition: IsotopeWaveletTransform.h:114
OpenMS::IsotopeWaveletTransform::IsotopeWaveletTransform
IsotopeWaveletTransform()
Default Constructor.
Definition: IsotopeWaveletTransform.h:519
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::TransSpectrum::MZEnd
MSSpectrum::const_iterator MZEnd(const double mz) const
Definition: IsotopeWaveletTransform.h:184
OpenMS::IsotopeWaveletTransform::TransSpectrum::getRefSpectrum
const MSSpectrum * getRefSpectrum() const
Definition: IsotopeWaveletTransform.h:170
OpenMS::intensityPointerComparator
bool intensityPointerComparator(PeakType *a, PeakType *b)
Definition: IsotopeWaveletTransform.h:507
OpenMS::Peak2D::setIntensity
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:172
OpenMS::Constants::IW_PROTON_MASS
const double IW_PROTON_MASS
Definition: IsotopeWaveletConstants.h:68
OpenMS::IsotopeWaveletTransform::BoxElement::c
UInt c
Note, this is not the charge (it is charge-1!!!)
Definition: IsotopeWaveletTransform.h:78
OpenMS::MSExperiment::size
Size size() const
Definition: MSExperiment.h:127
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::getMinSpacing
double getMinSpacing() const
Definition: IsotopeWaveletTransform.h:312
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::xs_
std::vector< double > xs_
Definition: IsotopeWaveletTransform.h:480
LinearRegression.h
OpenMS::Int
int Int
Signed integer type.
Definition: Types.h:102
OpenMS::ConvexHull2D
A 2-dimensional hull representation in [counter]clockwise direction - depending on axis labelling.
Definition: ConvexHull2D.h:72
OpenMS::Peak2D::setRT
void setRT(CoordinateType coordinate)
Mutable access to the RT coordinate (index 0)
Definition: Peak2D.h:214
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
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
Exception.h
OpenMS::IsotopeWaveletTransform::indices_
std::vector< int > indices_
Definition: IsotopeWaveletTransform.h:488
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::Feature::setConvexHulls
void setConvexHulls(const std::vector< ConvexHull2D > &hulls)
Set the convex hulls of single mass traces.
OpenMS::IsotopeWaveletTransform::psi_
std::vector< double > psi_
Definition: IsotopeWaveletTransform.h:480
OpenMS::IsotopeWaveletTransform::from_max_to_right_
Int from_max_to_right_
Definition: IsotopeWaveletTransform.h:487
OpenMS::IsotopeWaveletTransform::TransSpectrum::MZBegin
MSSpectrum::const_iterator MZBegin(const double mz) const
Definition: IsotopeWaveletTransform.h:177
OpenMS::IsotopeWaveletTransform::TransSpectrum::getRefSpectrum
const MSSpectrum * getRefSpectrum()
Definition: IsotopeWaveletTransform.h:164
OpenMS::IsotopeWaveletTransform::TransSpectrum::getRefIntensity
double getRefIntensity(const UInt i) const
Definition: IsotopeWaveletTransform.h:140
int
OpenMS::IsotopeWaveletTransform::scores_
std::vector< float > scores_
Definition: IsotopeWaveletTransform.h:491
OpenMS::IsotopeWaveletTransform::getClosedBoxes
virtual std::multimap< double, Box > getClosedBoxes()
Returns the closed boxes.
Definition: IsotopeWaveletTransform.h:276
OpenMS::IsotopeWaveletTransform::data_length_
UInt data_length_
Definition: IsotopeWaveletTransform.h:484
OpenMS::IsotopeWaveletTransform
A class implementing the isotope wavelet transform. If you just want to find features using the isoto...
Definition: IsotopeWaveletTransform.h:69
FeatureMap.h
OpenMS::IsotopeWaveletTransform::TransSpectrum::getRT
double getRT() const
Definition: IsotopeWaveletTransform.h:128
IsotopeWaveletConstants.h
OpenMS::IsotopeWaveletTransform::BoxElement::intens
double intens
The transformed intensity at the monoisotopic mass.
Definition: IsotopeWaveletTransform.h:80
OpenMS::IsotopeWaveletTransform::initializeScan
virtual void initializeScan(const MSSpectrum &c_ref, const UInt c=0)
Definition: IsotopeWaveletTransform.h:644
OpenMS::MSSpectrum::MZEnd
Iterator MZEnd(CoordinateType mz)
Binary search for peak range end (returns the past-the-end iterator)
OpenMS::IsotopeWaveletTransform::getAvMZSpacing_
double getAvMZSpacing_(const MSSpectrum &scan)
Computes the average MZ spacing of scan.
Definition: IsotopeWaveletTransform.h:1180
OpenMS::IsotopeWaveletTransform::TransSpectrum::TransSpectrum
TransSpectrum(const MSSpectrum *reference)
Definition: IsotopeWaveletTransform.h:107
OpenMS::MSSpectrum::MZBegin
Iterator MZBegin(CoordinateType mz)
Binary search for peak range begin.
OpenMS::IsotopeWaveletTransform::BoxElement::score
double score
The associated score.
Definition: IsotopeWaveletTransform.h:79
OpenMS::IsotopeWaveletTransform::end_boxes_
std::multimap< double, Box > end_boxes_
Definition: IsotopeWaveletTransform.h:476
OpenMS::Peak2D::setMZ
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:202
OpenMS::Feature::setOverallQuality
void setOverallQuality(QualityType q)
Set the overall quality.
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::TransSpectrum::begin
MSSpectrum::const_iterator begin() const
Definition: IsotopeWaveletTransform.h:198
OpenMS::IsotopeWaveletTransform::prod_
std::vector< double > prod_
Definition: IsotopeWaveletTransform.h:480
OpenMS::IsotopeWavelet::getMzPeakCutOffAtMonoPos
static UInt getMzPeakCutOffAtMonoPos(const double mass, const UInt z)
OpenMS::IsotopeWaveletTransform::from_max_to_left_
Int from_max_to_left_
Definition: IsotopeWaveletTransform.h:487
OpenMS::IsotopeWaveletTransform::~IsotopeWaveletTransform
virtual ~IsotopeWaveletTransform()
Destructor.
Definition: IsotopeWaveletTransform.h:557
OpenMS::IsotopeWavelet::getNumPeakCutOff
static UInt getNumPeakCutOff(const double mass, const UInt z)
OpenMS::IsotopeWaveletTransform::hr_data_
bool hr_data_
Definition: IsotopeWaveletTransform.h:485
OpenMS::IsotopeWaveletTransform::TransSpectrum::destroy
virtual void destroy()
Definition: IsotopeWaveletTransform.h:119
OpenMS::DPosition< 2 >
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::UInt
unsigned int UInt
Unsigned integer type.
Definition: Types.h:94
OpenMS::IsotopeWaveletTransform::TransSpectrum::getTransIntensity
double getTransIntensity(const UInt i) const
Definition: IsotopeWaveletTransform.h:146
OpenMS::IsotopeWaveletTransform::intenstype_
String intenstype_
Definition: IsotopeWaveletTransform.h:486
OpenMS::IsotopeWaveletTransform::checkPPMTheoModel_
virtual std::pair< double, double > checkPPMTheoModel_(const MSSpectrum &ref, const double c_mz, const UInt c)
Definition: IsotopeWaveletTransform.h:2149
MSExperiment.h
OpenMS::SignedSize
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:134
OpenMS::IsotopeWaveletTransform::TransSpectrum::setTransIntensity
void setTransIntensity(const UInt i, const double intens)
Definition: IsotopeWaveletTransform.h:152
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::FeatureMap
A container for features.
Definition: FeatureMap.h:97
OpenMS::ConstRefVector
This vector holds pointer to the elements of another container.
Definition: ConstRefVector.h:70
OpenMS::IsotopeWaveletTransform::TransSpectrum::TransSpectrum
TransSpectrum()
Definition: IsotopeWaveletTransform.h:101
OpenMS::intensityComparator
bool intensityComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:495
OpenMS::Feature
An LC-MS feature.
Definition: Feature.h:70
OpenMS::IsotopeWaveletTransform::sigma_
double sigma_
Definition: IsotopeWaveletTransform.h:479
OpenMS::IsotopeWavelet::init
static IsotopeWavelet * init(const double max_m, const UInt max_charge)
OpenMS::Math::mean
static double mean(IteratorType begin, IteratorType end)
Calculates the mean of a range of values.
Definition: StatisticFunctions.h:133
OpenMS::IsotopeWaveletTransform::closed_boxes_
std::multimap< double, Box > closed_boxes_
Definition: IsotopeWaveletTransform.h:476
OpenMS::ConstRefVector::sortByIntensity
void sortByIntensity(bool reverse=false)
Sorting.
Definition: ConstRefVector.h:741
ConstRefVector.h
OpenMS::IsotopeWaveletTransform::BoxElement::MZ_end
UInt MZ_end
Index.
Definition: IsotopeWaveletTransform.h:85
OpenMS::Constants::PEPTIDE_MASS_RULE_BOUND
const double PEPTIDE_MASS_RULE_BOUND
Definition: IsotopeWaveletConstants.h:50
OpenMS::Constants::IW_QUARTER_NEUTRON_MASS
const double IW_QUARTER_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:56
OpenMS::IsotopeWaveletTransform::max_charge_
UInt max_charge_
Definition: IsotopeWaveletTransform.h:484
OpenMS::IsotopeWaveletTransform::BoxElement::RT
double RT
The elution time (not the scan index)
Definition: IsotopeWaveletTransform.h:82
OpenMS::IsotopeWaveletTransform::setSigma
void setSigma(const double sigma)
Definition: IsotopeWaveletTransform.h:305
OpenMS::IsotopeWaveletTransform::TransSpectrum::size
Size size() const
Definition: IsotopeWaveletTransform.h:158
OpenMS::IsotopeWaveletTransform::getMaxScanSize
Size getMaxScanSize() const
Definition: IsotopeWaveletTransform.h:317
OpenMS::IsotopeWaveletTransform::BoxElement::MZ_begin
UInt MZ_begin
Index.
Definition: IsotopeWaveletTransform.h:84
OpenMS::Constants::IW_NEUTRON_MASS
const double IW_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:54
OpenMS::MSExperiment::clear
void clear(bool clear_meta_data)
Clears all data and meta data.
OpenMS::IsotopeWaveletTransform::sampleTheCMarrWavelet_
void sampleTheCMarrWavelet_(const MSSpectrum &scan, const Int wavelet_length, const Int mz_index, const UInt charge)
OpenMS::IsotopeWaveletTransform::BoxElement
Internally used data structure.
Definition: IsotopeWaveletTransform.h:75
OpenMS::IsotopeWaveletTransform::interpol_ys_
std::vector< double > interpol_ys_
Definition: IsotopeWaveletTransform.h:481
OpenMS::Peak2D::getMZ
CoordinateType getMZ() const
Returns the m/z coordinate (index 1)
Definition: Peak2D.h:196
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,...
OpenMS::IsotopeWaveletTransform::min_spacing_
double min_spacing_
Definition: IsotopeWaveletTransform.h:490
OpenMS::IsotopeWaveletTransform::Box
std::multimap< UInt, BoxElement > Box
Key: RT index, value: BoxElement.
Definition: IsotopeWaveletTransform.h:88
OpenMS::IsotopeWaveletTransform::TransSpectrum
Internally (only by GPUs) used data structure . It allows efficient data exchange between CPU and GPU...
Definition: IsotopeWaveletTransform.h:94
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::IsotopeWaveletTransform::max_scan_size_
Size max_scan_size_
Definition: IsotopeWaveletTransform.h:483
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::MSSpectrum
The representation of a 1D spectrum.
Definition: MSSpectrum.h:67
OpenMS::IsotopeWaveletTransform::max_num_peaks_per_pattern_
UInt max_num_peaks_per_pattern_
Definition: IsotopeWaveletTransform.h:484
OpenMS::IsotopeWaveletTransform::open_boxes_
std::multimap< double, Box > open_boxes_
Definition: IsotopeWaveletTransform.h:476
OpenMS::positionComparator
bool positionComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:513
StandardTypes.h
OpenMS::IsotopeWaveletTransform::computeMinSpacing
virtual void computeMinSpacing(const MSSpectrum &c_ref)
Definition: IsotopeWaveletTransform.h:686
OpenMS::IsotopeWaveletTransform::BoxElement::RT_index
UInt RT_index
The elution time (map) index.
Definition: IsotopeWaveletTransform.h:83
OpenMS::BaseFeature::setCharge
void setCharge(const ChargeType &ch)
Set charge state.
OpenMS::IsotopeWaveletTransform::TransSpectrum::trans_intens_
std::vector< float > * trans_intens_
The intensities of the transform.
Definition: IsotopeWaveletTransform.h:206
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::Constants::PEPTIDE_MASS_RULE_FACTOR
const double PEPTIDE_MASS_RULE_FACTOR
Definition: IsotopeWaveletConstants.h:49
OpenMS::Constants::PEPTIDE_MASS_RULE_THEO_PPM_BOUND
const double PEPTIDE_MASS_RULE_THEO_PPM_BOUND
Definition: IsotopeWaveletConstants.h:51
MSSpectrum.h
OpenMS::IsotopeWaveletTransform::BoxElement::mz
double mz
The monoisotopic position.
Definition: IsotopeWaveletTransform.h:77
OpenMS::IsotopeWavelet::getLambdaL
static double getLambdaL(const double m)
Returns the mass-parameter lambda (linear fit).
OpenMS::IsotopeWaveletTransform::zeros_
std::vector< float > zeros_
Definition: IsotopeWaveletTransform.h:491