Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
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-2013.
6 //
7 // This software is released under a three-clause BSD license:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of any author or any participating institution
14 // may be used to endorse or promote products derived from this software
15 // without specific prior written permission.
16 // For a full list of authors, refer to the file AUTHORS.
17 // --------------------------------------------------------------------------
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
22 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
25 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
27 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
28 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Rene Hussong$
32 // $Authors: $
33 // --------------------------------------------------------------------------
34 
35 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_ISOTOPEWAVELETTRANSFORM_H
36 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_ISOTOPEWAVELETTRANSFORM_H
37 
46 #include <cmath>
47 #include <math.h>
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 #ifdef OPENMS_HAS_CUDA
56 #include <cuda.h>
58 #endif
59 
60 // we are not yet sure if we really want to drag in cutil.h and the CUDA_SAFE_CALL definitions...
61 #ifndef CUDA_SAFE_CALL
62 #define CUDA_SAFE_CALL(call) call;
63 #endif
64 
65 namespace OpenMS
66 {
67 
69  class cudaHelp
70  {
71 public:
72  float getMZ()
73  { return mz; }
74  float getIntensity()
75  { return intens; }
76 
77  float mz;
78  float intens;
79  float score;
80  };
81 
82 
86  template <typename PeakType>
88  {
89 public:
90 
91 
93  struct BoxElement
94  {
95  DoubleReal mz; //<The monoisotopic position
96  UInt c; //<Note, this is not the charge (it is charge-1!!!)
97  DoubleReal score; //<The associated score
98  DoubleReal intens; //<The transformed intensity at the monoisotopic mass
100  DoubleReal RT; //<The elution time (not the scan index)
101  UInt RT_index; //<The elution time (map) index
102  UInt MZ_begin; //<Index
103  UInt MZ_end; //<Index
104  };
105 
106  typedef std::multimap<UInt, BoxElement> Box;
107 
108 
113  {
115 
116 public:
117 
121  {
122  }
123 
126  reference_(reference)
127  {
128  trans_intens_ = new std::vector<float>(reference_->size(), 0.0);
129  }
130 
132  virtual ~TransSpectrum()
133  {
134  delete (trans_intens_);
135  }
136 
137  virtual void destroy()
138  {
139  delete (trans_intens_);
141  delete (reference_);
142  reference_ = NULL;
143  }
144 
146  inline DoubleReal getRT() const
147  {
148  return reference_->getRT();
149  }
150 
152  inline DoubleReal getMZ(const UInt i) const
153  {
154  return (*reference_)[i].getMZ();
155  }
156 
158  inline DoubleReal getRefIntensity(const UInt i) const
159  {
160  return (*reference_)[i].getIntensity();
161  }
162 
164  inline DoubleReal getTransIntensity(const UInt i) const
165  {
166  return (*trans_intens_)[i];
167  }
168 
170  inline void setTransIntensity(const UInt i, const DoubleReal intens)
171  {
172  (*trans_intens_)[i] = intens;
173  }
174 
176  inline Size size() const
177  {
178  return trans_intens_->size();
179  }
180 
183  {
184  return reference_;
185  }
186 
188  inline const MSSpectrum<PeakType>* getRefSpectrum() const
189  {
190  return reference_;
191  }
192 
196  {
197  return reference_->MZBegin(mz);
198  }
199 
202  inline typename MSSpectrum<PeakType>::const_iterator MZEnd(const DoubleReal mz) const
203  {
204  return reference_->MZEnd(mz);
205  }
206 
210  {
211  return reference_->end();
212  }
213 
217  {
218  return reference_->begin();
219  }
220 
221 protected:
222 
223  const MSSpectrum<PeakType>* reference_; //<The reference spectrum
224  std::vector<float>* trans_intens_; //<The intensities of the transform
225 
226  };
227 
228 
229 
235  IsotopeWaveletTransform(const DoubleReal min_mz, const DoubleReal max_mz, const UInt max_charge, const Size max_scan_size = 0, const bool use_cuda = false, const bool hr_data = false, const String intenstype = "ref");
236 
238  virtual ~IsotopeWaveletTransform();
239 
240 
245  virtual void getTransform(MSSpectrum<PeakType>& c_trans, const MSSpectrum<PeakType>& c_ref, const UInt c);
246 
251  virtual void getTransformHighRes(MSSpectrum<PeakType>& c_trans, const MSSpectrum<PeakType>& c_ref, const UInt c);
252 
272  virtual void identifyCharge(const MSSpectrum<PeakType>& candidates, const MSSpectrum<PeakType>& ref, const UInt scan_index, const UInt c,
273  const DoubleReal ampl_cutoff, const bool check_PPMs);
274 
275  virtual void initializeScan(const MSSpectrum<PeakType>& c_ref, const UInt c = 0);
276 
277 #ifdef OPENMS_HAS_CUDA
278 
280  virtual int initializeScanCuda(const MSSpectrum<PeakType>& scan, const UInt c = 0);
281 
283  virtual void finalizeScanCuda();
284 
288  virtual void getTransformCuda(TransSpectrum& c_trans, const UInt c);
289 
305  virtual void identifyChargeCuda(const TransSpectrum& candidates, const UInt scan_index, const UInt c,
306  const DoubleReal ampl_cutoff, const bool check_PPMs);
307 
310  virtual int sortCuda(MSSpectrum<PeakType>& sorted);
311 #endif
312 
313 
320  void updateBoxStates(const MSExperiment<PeakType>& map, const Size scan_index, const UInt RT_interleave,
321  const UInt RT_votes_cutoff, const Int front_bound = -1, const Int end_bound = -1);
322 
323 
324  void mergeFeatures(IsotopeWaveletTransform<PeakType>* later_iwt, const UInt RT_interleave, const UInt RT_votes_cutoff);
325 
326 
331  FeatureMap<Feature> mapSeeds2Features(const MSExperiment<PeakType>& map, const UInt RT_votes_cutoff);
332 
334  virtual std::multimap<DoubleReal, Box> getClosedBoxes()
335  { return closed_boxes_; }
336 
337 
342  inline DoubleReal getLinearInterpolation(const typename MSSpectrum<PeakType>::const_iterator& left_iter, const DoubleReal mz_pos, const typename MSSpectrum<PeakType>::const_iterator& right_iter)
343  {
344  return left_iter->getIntensity() + (right_iter->getIntensity() - left_iter->getIntensity()) / (right_iter->getMZ() - left_iter->getMZ()) * (mz_pos - left_iter->getMZ());
345  }
346 
353  inline DoubleReal getLinearInterpolation(const DoubleReal mz_a, const DoubleReal intens_a, const DoubleReal mz_pos, const DoubleReal mz_b, const DoubleReal intens_b)
354  {
355  return intens_a + (intens_b - intens_a) / (mz_b - mz_a) * (mz_pos - mz_a);
356  }
357 
358  inline DoubleReal getSigma() const
359  {
360  return sigma_;
361  }
362 
363  inline void setSigma(const DoubleReal sigma)
364  {
365  sigma_ = sigma;
366  }
367 
368  virtual void computeMinSpacing(const MSSpectrum<PeakType>& c_ref);
369 
370  inline DoubleReal getMinSpacing() const
371  {
372  return min_spacing_;
373  }
374 
375  inline Size getMaxScanSize() const
376  {
377  return max_scan_size_;
378  }
379 
380 protected:
381 
382 
386 
387 
388  inline void sampleTheCMarrWavelet_(const MSSpectrum<PeakType>& scan, const Int wavelet_length, const Int mz_index, const UInt charge);
389 
390 
397  virtual DoubleReal scoreThis_(const TransSpectrum& candidate, const UInt peak_cutoff,
398  const DoubleReal seed_mz, const UInt c, const DoubleReal ampl_cutoff);
399 
406  virtual DoubleReal scoreThis_(const MSSpectrum<PeakType>& candidate, const UInt peak_cutoff,
407  const DoubleReal seed_mz, const UInt c, const DoubleReal ampl_cutoff);
408 
409 
417  virtual bool checkPositionForPlausibility_(const TransSpectrum& candidate, const MSSpectrum<PeakType>& ref, const DoubleReal seed_mz,
418  const UInt c, const UInt scan_index, const bool check_PPMs, const DoubleReal transintens, const DoubleReal prev_score);
419 
427  virtual bool checkPositionForPlausibility_(const MSSpectrum<PeakType>& candidate, const MSSpectrum<PeakType>& ref, const DoubleReal seed_mz,
428  const UInt c, const UInt scan_index, const bool check_PPMs, const DoubleReal transintens, const DoubleReal prev_score);
429 
430  virtual std::pair<DoubleReal, DoubleReal> checkPPMTheoModel_(const MSSpectrum<PeakType>& ref, const DoubleReal c_mz, const UInt c);
431 
432 
434  inline DoubleReal getAvIntens_(const TransSpectrum& scan);
436  inline DoubleReal getAvIntens_(const MSSpectrum<PeakType>& scan);
437 
439  inline DoubleReal getSdIntens_(const TransSpectrum& scan, const DoubleReal mean);
441  inline DoubleReal getSdIntens_(const MSSpectrum<PeakType>& scan, const DoubleReal mean);
442 
453  virtual void push2Box_(const DoubleReal mz, const UInt scan, UInt c, const DoubleReal score,
454  const DoubleReal intens, const DoubleReal rt, const UInt MZ_begin, const UInt MZ_end, const DoubleReal ref_intens);
455 
471  virtual void push2TmpBox_(const DoubleReal mz, const UInt scan, UInt charge, const DoubleReal score,
472  const DoubleReal intens, const DoubleReal rt, const UInt MZ_begin, const UInt MZ_end);
473 
479  inline DoubleReal getAvMZSpacing_(const MSSpectrum<PeakType>& scan);
480 
481 
486  void clusterSeeds_(const TransSpectrum& candidates, const MSSpectrum<PeakType>& ref,
487  const UInt scan_index, const UInt c, const bool check_PPMs);
488 
493  virtual void clusterSeeds_(const MSSpectrum<PeakType>& candidates, const MSSpectrum<PeakType>& ref,
494  const UInt scan_index, const UInt c, const bool check_PPMs);
495 
496 
501  void extendBox_(const MSExperiment<PeakType>& map, const Box box);
502 
505  inline DoubleReal peptideMassRule_(const DoubleReal c_mass) const
506  {
507  DoubleReal correction_fac = c_mass / Constants::PEPTIDE_MASS_RULE_BOUND;
508  DoubleReal old_frac_mass = c_mass - (Int)(c_mass);
509  DoubleReal new_mass = ((Int)(c_mass)) * (1. + Constants::PEPTIDE_MASS_RULE_FACTOR) - (Int)(correction_fac);
510  DoubleReal new_frac_mass = new_mass - (Int)(new_mass);
511 
512  if (new_frac_mass - old_frac_mass > 0.5)
513  {
514  new_mass -= 1.;
515  }
516 
517  if (new_frac_mass - old_frac_mass < -0.5)
518  {
519  new_mass += 1.;
520  }
521 
522  return new_mass;
523  }
524 
528  inline DoubleReal getPPMs_(const DoubleReal mass_a, const DoubleReal mass_b) const
529  {
530  return fabs(mass_a - mass_b) / (0.5 * (mass_a + mass_b)) * 1e6;
531  }
532 
533  //internally used data structures for the sweep line algorithm
534  std::multimap<DoubleReal, Box> open_boxes_, closed_boxes_, end_boxes_, front_boxes_; //DoubleReal = average m/z position
535  std::vector<std::multimap<DoubleReal, Box> >* tmp_boxes_; //for each charge we need a separate container
536 
538  std::vector<DoubleReal> c_mzs_, c_spacings_, psi_, prod_, xs_;
539  std::vector<DoubleReal> interpol_xs_, interpol_ys_;
540 
543  bool hr_data_;
546  std::vector<int> indices_;
547 
550  std::vector<float> scores_, zeros_;
551 
552 #ifdef OPENMS_HAS_CUDA
553  float* h_data_;
554  int* h_pos_;
555  UInt largest_array_size_, overall_size_, block_size_, to_load_, to_compute_;
556  Int num_elements_;
557  void* cuda_device_intens_;
558  void* cuda_device_pos_;
559  void* cuda_device_trans_intens_;
560  void* cuda_device_fwd2_;
561  void* cuda_device_posindices_sorted_;
562  void* cuda_device_trans_intens_sorted_;
563  void* cuda_device_scores_;
564  std::vector<float> cuda_positions_, cuda_intensities_;
565  dim3 dimGrid_, dimBlock_;
566 #endif
567  };
568 
569 
570  bool myCudaComparator(const cudaHelp& a, const cudaHelp& b);
571 
572  template <typename PeakType>
573  bool intensityComparator(const PeakType& a, const PeakType& b)
574  {
575  return a.getIntensity() > b.getIntensity();
576  }
577 
578  template <typename PeakType>
580  {
581  return a.getIntensity() < b.getIntensity();
582  }
583 
584  template <typename PeakType>
586  {
587  return a->getIntensity() > b->getIntensity();
588  }
589 
590  template <typename PeakType>
591  bool positionComparator(const PeakType& a, const PeakType& b)
592  {
593  return a.getMZ() < b.getMZ();
594  }
595 
596  template <typename PeakType>
598  {
599  tmp_boxes_ = new std::vector<std::multimap<DoubleReal, Box> >(1);
600  av_MZ_spacing_ = 1;
601  max_scan_size_ = 0;
602  max_mz_cutoff_ = 3;
603  max_num_peaks_per_pattern_ = 3;
604  hr_data_ = false;
605  intenstype_ = "ref";
606 #ifdef OPENMS_HAS_CUDA
607  largest_array_size_ = 0;
608  num_elements_ = 0;
609  h_data_ = NULL;
610  h_pos_ = NULL;
611 #endif
612  }
613 
614  template <typename PeakType>
615  IsotopeWaveletTransform<PeakType>::IsotopeWaveletTransform(const DoubleReal min_mz, const DoubleReal max_mz, const UInt max_charge, const Size max_scan_size, const bool use_cuda, const bool hr_data, String intenstype)
616  {
617  max_charge_ = max_charge;
618  max_scan_size_ = max_scan_size;
619  hr_data_ = hr_data;
620  intenstype_ = intenstype;
621  tmp_boxes_ = new std::vector<std::multimap<DoubleReal, Box> >(max_charge);
622  if (max_scan_size <= 0) //only important for the CPU
623  {
624  IsotopeWavelet::init(max_mz, max_charge);
625  }
626 
627  av_MZ_spacing_ = 1;
628  max_mz_cutoff_ = IsotopeWavelet::getMzPeakCutOffAtMonoPos(max_mz, max_charge);
629  max_num_peaks_per_pattern_ = IsotopeWavelet::getNumPeakCutOff(max_mz, max_charge);
630 
631 #ifdef OPENMS_HAS_CUDA
632  if (use_cuda) //only important for the GPU
633  {
634  if (hr_data_)
635  {
636  max_scan_size_ = max_scan_size * (4 * (max_mz_cutoff_ - 1) - 1);
637  }
638  largest_array_size_ = pow(2, ceil(log(max_scan_size_ + Constants::CUDA_EXTENDED_BLOCK_SIZE_MAX) / log(2.0)));
639 
640  cuda_positions_.reserve(largest_array_size_);
641  cuda_intensities_.reserve(largest_array_size_);
642  indices_.resize(largest_array_size_);
643  for (UInt q = 0; q < largest_array_size_; ++q)
644  {
645  indices_[q] = q;
646  }
647 
648  h_data_ = (float*) malloc(largest_array_size_ * sizeof(float));
649  h_pos_ = (int*) malloc(largest_array_size_ * sizeof(int));
650  }
651  else
652  {
653  h_data_ = NULL;
654  h_pos_ = NULL;
655  Int size_estimate((Int)ceil(max_scan_size_ / (max_mz - min_mz)));
656  Int to_reserve((Int)ceil(size_estimate * max_num_peaks_per_pattern_ * Constants::IW_NEUTRON_MASS));
657  psi_.reserve(to_reserve); //The wavelet
658  prod_.reserve(to_reserve);
659  xs_.reserve(to_reserve);
662  }
663 #else
664  if (!use_cuda)
665  {
666  Int size_estimate((Int)ceil(max_scan_size_ / (max_mz - min_mz)));
667  Int to_reserve((Int)ceil(size_estimate * max_num_peaks_per_pattern_ * Constants::IW_NEUTRON_MASS));
668  psi_.reserve(to_reserve); //The wavelet
669  prod_.reserve(to_reserve);
670  xs_.reserve(to_reserve);
673  }
674 #endif
675  }
676 
677  template <typename PeakType>
679  {
680 #ifdef OPENMS_HAS_CUDA
681  if (h_data_ != NULL)
682  free(h_data_);
683  if (h_pos_ != NULL)
684  free(h_pos_);
685  h_data_ = NULL;
686  h_pos_ = NULL;
687 #endif
688 
689  delete (tmp_boxes_);
690  }
691 
692  template <typename PeakType>
694  {
695  Int spec_size((Int)c_ref.size());
696  //in the very unlikely case that size_t will not fit to int anymore this will be a problem of course
697  //for the sake of simplicity (we need here a signed int) we do not cast at every following comparison individually
698  UInt charge = c + 1;
699  DoubleReal value, T_boundary_left, T_boundary_right, old, c_diff, current, old_pos, my_local_MZ, my_local_lambda, origin, c_mz;
700 
701  for (Int my_local_pos = 0; my_local_pos < spec_size; ++my_local_pos)
702  {
703  value = 0; T_boundary_left = 0, T_boundary_right = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_ref[my_local_pos].getMZ(), charge) / (DoubleReal)charge;
704  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_;
705  my_local_MZ = c_ref[my_local_pos].getMZ(); my_local_lambda = IsotopeWavelet::getLambdaL(my_local_MZ * charge);
706  c_diff = 0;
707  origin = -my_local_MZ + Constants::IW_QUARTER_NEUTRON_MASS / (DoubleReal)charge;
708 
709  for (Int current_conv_pos = std::max(0, my_local_pos - from_max_to_left_); c_diff < T_boundary_right; ++current_conv_pos)
710  {
711  if (current_conv_pos >= spec_size)
712  {
713  value += 0.5 * old * min_spacing_;
714  break;
715  }
716 
717  c_mz = c_ref[current_conv_pos].getMZ();
718  c_diff = c_mz + origin;
719 
720  //Attention! The +1. has nothing to do with the charge, it is caused by the wavelet's formula (tz1).
721  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;
722 
723  value += 0.5 * (current + old) * (c_mz - old_pos);
724 
725  old = current;
726  old_pos = c_mz;
727  }
728 
729 
730 
731  c_trans[my_local_pos].setIntensity(value);
732  }
733  }
734 
735  template <typename PeakType>
737  {
738  Int spec_size((Int)c_ref.size());
739  //in the very unlikely case that size_t will not fit to int anymore this will be a problem of course
740  //for the sake of simplicity (we need here a signed int) we do not cast at every following comparison individually
741  UInt charge = c + 1;
742  DoubleReal value, T_boundary_left, T_boundary_right, c_diff, current, my_local_MZ, my_local_lambda, origin, c_mz;
743 
744  for (Int my_local_pos = 0; my_local_pos < spec_size; ++my_local_pos)
745  {
746  value = 0; T_boundary_left = 0, T_boundary_right = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_ref[my_local_pos].getMZ(), charge) / (DoubleReal)charge;
747 
748 
749  my_local_MZ = c_ref[my_local_pos].getMZ(); my_local_lambda = IsotopeWavelet::getLambdaL(my_local_MZ * charge);
750  c_diff = 0;
751  origin = -my_local_MZ + Constants::IW_QUARTER_NEUTRON_MASS / (DoubleReal)charge;
752 
753  for (Int current_conv_pos = std::max(0, my_local_pos - from_max_to_left_); c_diff < T_boundary_right; ++current_conv_pos)
754  {
755  if (current_conv_pos >= spec_size)
756  {
757  break;
758  }
759 
760  c_mz = c_ref[current_conv_pos].getMZ();
761  c_diff = c_mz + origin;
762 
763  //Attention! The +1. has nothing to do with the charge, it is caused by the wavelet's formula (tz1).
764  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;
765 
766  value += current;
767  }
768 
769  c_trans[my_local_pos].setIntensity(value);
770  }
771  }
772 
773  template <typename PeakType>
775  {
776  data_length_ = (UInt) c_ref.size();
777  computeMinSpacing(c_ref);
778  Int wavelet_length = 0, quarter_length = 0;
779 
780  if (hr_data_) //We have to check this separately, because the simply estimation for LowRes data is destroyed by large gaps
781  {
782  UInt c_mz_cutoff;
783  typename MSSpectrum<PeakType>::const_iterator start_iter, end_iter;
784  for (UInt i = 0; i < data_length_; ++i)
785  {
786  c_mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_ref[i].getMZ(), c + 1);
787  start_iter = c_ref.MZEnd(c_ref[i].getMZ());
788  end_iter = c_ref.MZBegin(c_ref[i].getMZ() + c_mz_cutoff);
789  wavelet_length = std::max((SignedSize) wavelet_length, distance(start_iter, end_iter) + 1);
790  end_iter = c_ref.MZEnd(c_ref[i].getMZ() - Constants::IW_QUARTER_NEUTRON_MASS / DoubleReal(c + 1.));
791  quarter_length = std::max((SignedSize) quarter_length, distance(end_iter, start_iter) + 1);
792  }
793  }
794  else
795  {
796  //CHANGED
797  max_mz_cutoff_ = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_ref[data_length_ - 1].getMZ(), max_charge_);
798  wavelet_length = (UInt) ceil(max_mz_cutoff_ / min_spacing_);
799  }
800  //... done
801 
802 
803  if (wavelet_length > (Int) c_ref.size())
804  {
805  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;
806  std::cout << "Minimal spacing: " << min_spacing_ << std::endl;
807  std::cout << "Warning/Error generated at scan with RT " << c_ref.getRT() << "." << std::endl;
808  }
809 
810  Int max_index = (UInt) (Constants::IW_QUARTER_NEUTRON_MASS / min_spacing_);
811  from_max_to_left_ = max_index;
812  from_max_to_right_ = wavelet_length - 1 - from_max_to_left_;
813  }
814 
815  template <typename PeakType>
817  {
818  min_spacing_ = INT_MAX;
819  for (UInt c_conv_pos = 1; c_conv_pos < c_ref.size(); ++c_conv_pos)
820  {
821  min_spacing_ = std::min(min_spacing_, c_ref[c_conv_pos].getMZ() - c_ref[c_conv_pos - 1].getMZ());
822  }
823  }
824 
825 #ifdef OPENMS_HAS_CUDA
826  template <typename PeakType>
828  {
829  (cudaFree(cuda_device_pos_));
830  (cudaFree(cuda_device_intens_));
831  (cudaFree(cuda_device_trans_intens_));
832  (cudaFree(cuda_device_fwd2_));
833  (cudaFree(cuda_device_trans_intens_sorted_));
834  (cudaFree(cuda_device_posindices_sorted_));
835  (cudaFree(cuda_device_scores_));
836  }
837 
838  template <typename PeakType>
839  int IsotopeWaveletTransform<PeakType>::initializeScanCuda(const MSSpectrum<PeakType>& scan, const UInt c)
840  {
841  data_length_ = scan.size();
842 
843  std::vector<float> pre_positions(data_length_), pre_intensities(data_length_);
844  float c_spacing;
845  min_spacing_ = INT_MAX;
846  pre_positions[0] = scan[0].getMZ();
847  pre_intensities[0] = scan[0].getIntensity();
848 
849  for (UInt i = 1; i < data_length_; ++i)
850  {
851  pre_positions[i] = scan[i].getMZ();
852  c_spacing = pre_positions[i] - pre_positions[i - 1];
853  if (c_spacing < min_spacing_)
854  {
855  min_spacing_ = c_spacing;
856  }
857  pre_intensities[i] = scan[i].getIntensity();
858  }
859  if (min_spacing_ == INT_MAX) //spectrum consists of a single data point
860  {
861  std::cout << "Scan consits of a single point. Unable to compute transform." << std::endl;
863  }
864 
865  //Estimating the wave_length ...
866  UInt wavelet_length = 0, quarter_length = 0;
867  if (hr_data_) //We have to check this separately, because the simply estimation for LowRes data is destroyed by large gaps
868  {
869  UInt c_mz_cutoff;
870  typename MSSpectrum<PeakType>::const_iterator start_iter, end_iter;
871  for (UInt i = 0; i < data_length_; ++i)
872  {
873  c_mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(scan[i].getMZ(), c + 1);
874  start_iter = scan.MZEnd(scan[i].getMZ());
875  end_iter = scan.MZBegin(scan[i].getMZ() + c_mz_cutoff);
876  wavelet_length = std::max((long int) wavelet_length, distance(start_iter, end_iter) + 1);
877  end_iter = scan.MZEnd(scan[i].getMZ() - Constants::IW_QUARTER_NEUTRON_MASS / DoubleReal(c + 1.));
878  quarter_length = std::max((long int) quarter_length, distance(end_iter, start_iter) + 1);
879  }
880  }
881  else
882  {
883  //CHANGED
884  max_mz_cutoff_ = IsotopeWavelet::getMzPeakCutOffAtMonoPos(scan[data_length_ - 1].getMZ(), max_charge_);
885  wavelet_length = (UInt) ceil(max_mz_cutoff_ / min_spacing_);
886  }
887  //... done
888 
889  if (wavelet_length > data_length_ || wavelet_length == 1) //==1, because of 'ceil'
890  {
891  std::cout << "Warning: the extremal length of the wavelet is larger (" << wavelet_length << ") than the number of data points (" << data_length_ << "). This might (!) severely affect the transform." << std::endl;
892  std::cout << "Minimal spacing: " << min_spacing_ << std::endl;
893  std::cout << "Warning/Error generated at scan with RT " << scan.getRT() << "." << std::endl;
894  }
895 
896  UInt max_index;
897  if (hr_data_)
898  {
899  max_index = quarter_length;
900  }
901  else
902  {
903  max_index = (UInt) (Constants::IW_QUARTER_NEUTRON_MASS / min_spacing_);
904  }
905  from_max_to_left_ = max_index;
906  from_max_to_right_ = wavelet_length - 1 - from_max_to_left_;
907 
908  Int problem_size = Constants::CUDA_BLOCK_SIZE_MAX;
909  to_load_ = problem_size + from_max_to_left_ + from_max_to_right_;
910 
911  UInt missing_points = problem_size - (data_length_ % problem_size);
912  overall_size_ = wavelet_length - 1 + data_length_ + missing_points;
913 
914  num_elements_ = overall_size_;
915  Int dev_num_elements = 1, tmp = overall_size_ >> 1;
916 
917  //Get power of 2 elements (necessary for the sorting algorithm)
918  while (tmp)
919  {
920  dev_num_elements <<= 1;
921  tmp >>= 1;
922  }
923 
924  if (num_elements_ > dev_num_elements)
925  {
926  dev_num_elements <<= 1;
927  }
928 
929  if (dev_num_elements < Constants::CUDA_MIN_SORT_SIZE)
930  {
931  dev_num_elements = Constants::CUDA_MIN_SORT_SIZE;
932  }
933 
934  overall_size_ = dev_num_elements;
935 
936  cuda_intensities_.resize(overall_size_, 0); cuda_positions_.resize(overall_size_, 0);
937  //Pad the values to the left; the positions should not matter if the values are zero
938  float first_pos = pre_positions[0];
939  for (Int i = 0; i < from_max_to_left_; ++i)
940  {
941  cuda_positions_[i] = first_pos - (from_max_to_left_ - i) * min_spacing_;
942  }
943 
944  for (UInt i = 0; i < data_length_; ++i)
945  {
946  cuda_positions_[from_max_to_left_ + i] = pre_positions[i];
947  cuda_intensities_[from_max_to_left_ + i] = pre_intensities[i];
948  }
949 
950  float last_pos = pre_positions[pre_positions.size() - 1];
951  for (UInt i = 0; i < missing_points + from_max_to_right_ + dev_num_elements - num_elements_; ++i)
952  {
953  cuda_positions_[from_max_to_left_ + data_length_ + i] = last_pos + (i + 1) * min_spacing_;
954  }
955 
956 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
957  std::stringstream name; name << "cuda_input_" << scan.getRT() << ".out\0";
958  std::fstream outfile(name.str().c_str(), std::ios::out);
959  for (size_t i = 0; i < overall_size_; ++i)
960  outfile << cuda_positions_[i] << " " << cuda_intensities_[i] << std::endl;
961  outfile.close();
962 #endif
963 
964 
965  dimBlock_ = dim3(Constants::CUDA_BLOCK_SIZE_MAX);
966  to_compute_ = problem_size;
967  dimGrid_ = dim3((data_length_ + missing_points) / problem_size);
968 
969  (cudaMalloc(&cuda_device_posindices_sorted_, overall_size_ * sizeof(int)));
970  (cudaMalloc(&cuda_device_pos_, overall_size_ * sizeof(float)));
971  (cudaMemcpy(cuda_device_pos_, &(cuda_positions_[0]), overall_size_ * sizeof(float), cudaMemcpyHostToDevice));
972  (cudaMalloc(&cuda_device_intens_, overall_size_ * sizeof(float)));
973  (cudaMemcpy(cuda_device_intens_, &(cuda_intensities_[0]), overall_size_ * sizeof(float), cudaMemcpyHostToDevice));
974  (cudaMalloc(&cuda_device_trans_intens_, overall_size_ * sizeof(float)));
975  (cudaMalloc(&cuda_device_fwd2_, overall_size_ * sizeof(float)));
976  (cudaMalloc(&cuda_device_trans_intens_sorted_, overall_size_ * sizeof(float)));
977 
978  c_sorted_candidate_.resize(overall_size_);
979  scores_.resize(data_length_);
980  zeros_.resize(overall_size_);
981  memset(&zeros_[0], 0., overall_size_ * sizeof(float));
982 
983  (cudaMalloc(&cuda_device_scores_, overall_size_ * sizeof(float)));
984 
986  }
987 
988  template <typename PeakType>
989  void IsotopeWaveletTransform<PeakType>::getTransformCuda(TransSpectrum& c_trans, const UInt c)
990  {
991 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
992  std::cout << "res in vector" << std::endl;
993  std::vector<float> res(overall_size_, 0);
994 #endif
995  (cudaMemcpy(cuda_device_trans_intens_, &zeros_[0], overall_size_ * sizeof(float), cudaMemcpyHostToDevice));
996 
997  (cudaMemcpy(cuda_device_fwd2_, &zeros_[0], overall_size_ * sizeof(float), cudaMemcpyHostToDevice));
998  getExternalCudaTransforms(dimGrid_, dimBlock_, (float*)cuda_device_pos_, (float*)cuda_device_intens_, from_max_to_left_, from_max_to_right_, (float*)cuda_device_trans_intens_,
999  c + 1, to_load_, to_compute_, data_length_, (float*)cuda_device_fwd2_, hr_data_);
1000 
1001  (cudaMemcpy(cuda_device_trans_intens_sorted_, cuda_device_fwd2_, overall_size_ * sizeof(float), cudaMemcpyDeviceToDevice));
1002 
1003  (cudaMemcpy(&((*c_trans.trans_intens_)[0]), (float*)cuda_device_trans_intens_ + from_max_to_left_, data_length_ * sizeof(float), cudaMemcpyDeviceToHost));
1004 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1005  (cudaMemcpy(&(res[0]), (float*)cuda_device_trans_intens_ + from_max_to_left_, data_length_ * sizeof(float), cudaMemcpyDeviceToHost));
1006  for (UInt i = 0; i < data_length_; ++i)
1007  {
1008  c_trans.setTransIntensity(i, res[i]);
1009  }
1010 #endif
1011  }
1012 
1013  template <typename PeakType>
1014  int IsotopeWaveletTransform<PeakType>::sortCuda(MSSpectrum<PeakType>& sorted)
1015  {
1016  (cudaMemcpy(cuda_device_posindices_sorted_, &indices_[0], overall_size_ * sizeof(int), cudaMemcpyHostToDevice));
1017  Int gpu_index = sortOnDevice((float*)cuda_device_trans_intens_sorted_, (int*) cuda_device_posindices_sorted_, overall_size_, 0);
1018 
1019  if (gpu_index < 0) //i.e., there is no positive intensity value at all
1020  {
1021  return gpu_index;
1022  }
1023 
1024  (cudaMemcpy(h_data_, (float*)cuda_device_trans_intens_sorted_ + gpu_index, sizeof(float) * (overall_size_ - gpu_index), cudaMemcpyDeviceToHost));
1025  (cudaMemcpy(h_pos_, (int*)cuda_device_posindices_sorted_ + gpu_index, sizeof(int) * (overall_size_ - gpu_index), cudaMemcpyDeviceToHost));
1026 
1027  for (UInt i = 0; i < (overall_size_ - gpu_index); ++i)
1028  {
1029  sorted[i].setIntensity(h_data_[i]);
1030  sorted[i].setMZ(cuda_positions_[h_pos_[i]]);
1031  }
1032 
1033  return gpu_index;
1034  }
1035 
1036  template <typename PeakType>
1037  void IsotopeWaveletTransform<PeakType>::identifyChargeCuda(const TransSpectrum& candidates,
1038  const UInt scan_index, const UInt c, const DoubleReal ampl_cutoff, const bool check_PPMs)
1039  {
1040  const MSSpectrum<PeakType>& ref(*candidates.getRefSpectrum());
1041  UInt index, MZ_start, MZ_end;
1042  typename MSSpectrum<PeakType>::iterator iter, bound_iter;
1043  typename MSSpectrum<PeakType>::const_iterator iter_start, iter_end, iter_p, iter2, seed_iter;
1044  DoubleReal mz_cutoff, seed_mz, c_av_intens = 0, c_score = 0, c_sd_intens = 0, threshold = 0, help_mz;
1045 
1046  Int gpu_index = sortCuda(c_sorted_candidate_), c_index;
1047  if (gpu_index < 0) //the transform produced non-exploitable data
1048  {
1049  return;
1050  }
1051 
1052  std::vector<UInt> processed(data_length_, 0);
1053  if (ampl_cutoff < 0)
1054  {
1055  threshold = 0;
1056  }
1057  else
1058  {
1059  c_av_intens = getAvIntens_(candidates);
1060  c_sd_intens = getSdIntens_(candidates, c_av_intens);
1061  threshold = ampl_cutoff * c_sd_intens + c_av_intens;
1062  }
1063 
1064  Int num_of_scores = overall_size_ - gpu_index;
1065 
1066  (cudaMemcpy(cuda_device_scores_, &zeros_[0], num_of_scores * sizeof(float), cudaMemcpyHostToDevice));
1067 
1068  scoreOnDevice((int*)cuda_device_posindices_sorted_, (float*)cuda_device_trans_intens_, (float*)cuda_device_pos_, (float*)cuda_device_scores_,
1069  c, num_of_scores, overall_size_, max_num_peaks_per_pattern_, threshold);
1070 
1071  (cudaMemcpy(&scores_[0], cuda_device_scores_, num_of_scores * sizeof(float), cudaMemcpyDeviceToHost));
1072 
1073  std::vector<float>::iterator score_iter;
1074 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1075  std::stringstream stream;
1076  stream << "sorted_gpu_" << candidates.getRT() << "_" << c + 1 << ".trans\0";
1077  std::ofstream ofile(stream.str().c_str());
1078  for (c_index = overall_size_ - gpu_index - 1, score_iter = scores_.begin() + num_of_scores - 1; c_index >= 0; --c_index, --score_iter)
1079  {
1080  ofile << c_sorted_candidate_[c_index].getMZ() << "\t" << c_sorted_candidate_[c_index].getIntensity() << "\t" << *score_iter << std::endl;
1081  }
1082  ofile.close();
1083 #endif
1084 
1085  for (c_index = overall_size_ - gpu_index - 1, score_iter = scores_.begin() + num_of_scores - 1; c_index >= 0; --c_index, --score_iter)
1086  {
1087  seed_mz = c_sorted_candidate_[c_index].getMZ();
1088 
1089  //We can replace the following two lines ...
1090  //seed_iter = ref.MZBegin(seed_mz);
1091  //index = distance(ref.begin(), seed_iter);
1092  //... with:
1093  index = h_pos_[c_index] - from_max_to_left_;
1094  seed_iter = ref.begin() + index;
1095 
1096  if (seed_iter == ref.end() || processed[distance(ref.begin(), seed_iter)] || index <= 0)
1097  {
1098  continue;
1099  }
1100 
1101  mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(seed_mz, c + 1);
1102  //Mark the region as processed
1103  //Do not move this further down, since we have to mark this as processed in any case,
1104  //even when score <=0; otherwise we would look around the maximum's position unless
1105  //any significant point is found
1106  iter_start = ref.MZBegin(ref.begin(), seed_mz - Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.), seed_iter);
1107  iter_end = ref.MZEnd(seed_iter, seed_mz + mz_cutoff / (c + 1.), ref.end());
1108 
1109  if (iter_end == ref.end())
1110  {
1111  --iter_end;
1112  }
1113 
1114  MZ_start = distance(ref.begin(), iter_start);
1115  MZ_end = distance(ref.begin(), iter_end);
1116 
1117  memset(&(processed[MZ_start]), 1, sizeof(UInt) * (MZ_end - MZ_start + 1));
1118 
1119  c_score = *score_iter;
1120 
1121  if (c_score <= 0 && c_score != -1000)
1122  {
1123  continue;
1124  }
1125 
1126  //Push the seed into its corresponding box (or create a new one, if necessary)
1127  //Do ***NOT*** move this further down!
1128 
1129  push2TmpBox_(seed_mz, scan_index, c, c_score, c_sorted_candidate_[c_index].getIntensity(), ref.getRT(), MZ_start, MZ_end);
1130 
1131  //Push neighboring peaks to compute finally a derivative over the isotope pattern envelope
1132  help_mz = seed_mz - Constants::IW_NEUTRON_MASS / (c + 1.);
1133  iter2 = candidates.MZBegin(help_mz);
1134 
1135  if (iter2 == candidates.end() || iter2 == candidates.begin())
1136  {
1137  continue;
1138  }
1139 
1140  if (fabs(iter2->getMZ() - seed_mz) > 0.5 * Constants::IW_NEUTRON_MASS / (c + 1.))
1141  {
1142  //In the other case, we are too close to the peak, leading to incorrect derivatives.
1143  if (iter2 != candidates.end())
1144  {
1145  UInt dist = distance(candidates.begin(), iter2);
1146  push2TmpBox_(iter2->getMZ(), scan_index, c, 0,
1147  getLinearInterpolation((iter2 - 1)->getMZ(), candidates.getTransIntensity(dist - 1), help_mz, iter2->getMZ(), candidates.getTransIntensity(dist)),
1148  candidates.getRT(), MZ_start, MZ_end);
1149  }
1150  }
1151 
1152  help_mz = seed_mz + Constants::IW_NEUTRON_MASS / (c + 1.);
1153  iter2 = candidates.MZBegin(help_mz);
1154 
1155  if (iter2 == candidates.end() || iter2 == candidates.begin())
1156  {
1157  continue;
1158  }
1159 
1160  if (fabs(iter2->getMZ() - seed_mz) > 0.5 * Constants::IW_NEUTRON_MASS / (c + 1.))
1161  {
1162  //In the other case, we are too close to the peak, leading to incorrect derivatives.
1163  if (iter2 != candidates.end())
1164  {
1165  UInt dist = distance(candidates.begin(), iter2);
1166  push2TmpBox_(iter2->getMZ(), scan_index, c, 0,
1167  getLinearInterpolation((iter2 - 1)->getMZ(), candidates.getTransIntensity(dist - 1), help_mz, iter2->getMZ(), candidates.getTransIntensity(dist)),
1168  candidates.getRT(), MZ_start, MZ_end);
1169  }
1170  }
1171  }
1172 
1173  clusterSeeds_(candidates, ref, scan_index, c, check_PPMs);
1174  }
1175 
1176 #endif
1177 
1178 
1179  template <typename PeakType>
1181  const MSSpectrum<PeakType>& ref, const UInt scan_index, const UInt c, const DoubleReal ampl_cutoff, const bool check_PPMs)
1182  {
1183  Size scan_size(candidates.size());
1184  typename ConstRefVector<MSSpectrum<PeakType> >::iterator iter;
1185  typename MSSpectrum<PeakType>::const_iterator iter_start, iter_end, iter_p, seed_iter, iter2;
1186  DoubleReal mz_cutoff, seed_mz, c_av_intens = 0, c_score = 0, c_sd_intens = 0, threshold = 0, help_mz, share, share_pos, bwd, fwd;
1187  UInt MZ_start, MZ_end;
1188 
1189  MSSpectrum<PeakType> diffed(candidates);
1190  diffed[0].setIntensity(0); diffed[scan_size - 1].setIntensity(0);
1191 
1192 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1193  std::stringstream stream;
1194  stream << "diffed_" << ref.getRT() << "_" << c + 1 << ".trans\0";
1195  std::ofstream ofile(stream.str().c_str());
1196 #endif
1197 
1198  if (!hr_data_) //LowRes data
1199  {
1200  for (UInt i = 0; i < scan_size - 2; ++i)
1201  {
1202  share = candidates[i + 1].getIntensity(), share_pos = candidates[i + 1].getMZ();
1203  bwd = (share - candidates[i].getIntensity()) / (share_pos - candidates[i].getMZ());
1204  fwd = (candidates[i + 2].getIntensity() - share) / (candidates[i + 2].getMZ() - share_pos);
1205 
1206  if (!(bwd >= 0 && fwd <= 0) || share > ref[i + 1].getIntensity())
1207  {
1208  diffed[i + 1].setIntensity(0);
1209  }
1210 
1211 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1212  ofile << diffed[i + 1].getMZ() << "\t" << diffed[i + 1].getIntensity() << std::endl;
1213 #endif
1214  }
1215  }
1216  else //HighRes data
1217  {
1218  for (UInt i = 0; i < scan_size - 2; ++i)
1219  {
1220  share = candidates[i + 1].getIntensity(), share_pos = candidates[i + 1].getMZ();
1221  bwd = (share - candidates[i].getIntensity()) / (share_pos - candidates[i].getMZ());
1222  fwd = (candidates[i + 2].getIntensity() - share) / (candidates[i + 2].getMZ() - share_pos);
1223 
1224  if (!(bwd >= 0 && fwd <= 0))
1225  {
1226  diffed[i + 1].setIntensity(0);
1227  }
1228 
1229 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1230  ofile << diffed[i + 1].getMZ() << "\t" << diffed[i + 1].getIntensity() << std::endl;
1231 #endif
1232  }
1233  }
1234 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1235  ofile.close();
1236 #endif
1237 
1238  ConstRefVector<MSSpectrum<PeakType> > c_sorted_candidate_(diffed.begin(), diffed.end());
1239 
1240  //Sort the transform in descending order according to the intensities present in the transform
1241  c_sorted_candidate_.sortByIntensity();
1242 
1243 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1244  std::stringstream stream2;
1245  stream2 << "sorted_cpu_" << candidates.getRT() << "_" << c + 1 << ".trans\0";
1246  std::ofstream ofile2(stream2.str().c_str());
1247  for (iter = c_sorted_candidate_.end() - 1; iter != c_sorted_candidate_.begin(); --iter)
1248  {
1249  ofile2 << iter->getMZ() << "\t" << iter->getIntensity() << std::endl;
1250  }
1251  ofile2.close();
1252 #endif
1253 
1254  std::vector<UInt> processed(scan_size, 0);
1255 
1256  if (ampl_cutoff < 0)
1257  {
1258  threshold = 0;
1259  }
1260  else
1261  {
1262  c_av_intens = getAvIntens_(candidates);
1263  c_sd_intens = getSdIntens_(candidates, c_av_intens);
1264  threshold = ampl_cutoff * c_sd_intens + c_av_intens;
1265  }
1266 
1267  for (iter = c_sorted_candidate_.end() - 1; iter != c_sorted_candidate_.begin(); --iter)
1268  {
1269  if (iter->getIntensity() <= 0)
1270  {
1271  break;
1272  }
1273 
1274  seed_mz = iter->getMZ();
1275  seed_iter = ref.MZBegin(seed_mz);
1276 
1277  if (seed_iter == ref.end() || processed[distance(ref.begin(), seed_iter)])
1278  {
1279  continue;
1280  }
1281 
1282  mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(seed_mz, c + 1);
1283  //Mark the region as processed
1284  //Do not move this further down, since we have to mark this as processed in any case,
1285  //even when score <=0; otherwise we would look around the maximum's position unless
1286  //any significant point is found
1287  iter_start = ref.MZBegin(ref.begin(), seed_mz - Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.), seed_iter);
1288  iter_end = ref.MZEnd(seed_iter, seed_mz + mz_cutoff / (c + 1.), ref.end());
1289  if (iter_end == ref.end())
1290  {
1291  --iter_end;
1292  }
1293 
1294  MZ_start = distance(ref.begin(), iter_start);
1295  MZ_end = distance(ref.begin(), iter_end);
1296 
1297  memset(&(processed[MZ_start]), 1, sizeof(UInt) * (MZ_end - MZ_start + 1));
1298 
1299  c_score = scoreThis_(candidates, IsotopeWavelet::getNumPeakCutOff(seed_mz * (c + 1.)), seed_mz, c, threshold);
1300 
1301 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1302  if (trunc(seed_mz) == 874)
1303  std::cout << seed_mz << "\t" << c_score << std::endl;
1304 #endif
1305 
1306  if (c_score <= 0 && c_score != -1000)
1307  {
1308  continue;
1309  }
1310 
1311  //Push the seed into its corresponding box (or create a new one, if necessary)
1312  //Do ***NOT*** move this further down!
1313  push2TmpBox_(seed_mz, scan_index, c, c_score, iter->getIntensity(), ref.getRT(), MZ_start, MZ_end);
1314 
1315  help_mz = seed_mz - Constants::IW_NEUTRON_MASS / (c + 1.);
1316  iter2 = candidates.MZBegin(help_mz);
1317  if (iter2 == candidates.end() || iter2 == candidates.begin())
1318  {
1319  continue;
1320  }
1321 
1322  if (fabs(iter2->getMZ() - seed_mz) > 0.5 * Constants::IW_NEUTRON_MASS / (c + 1.))
1323  {
1324  //In the other case, we are too close to the peak, leading to incorrect derivatives.
1325  if (iter2 != candidates.end())
1326  {
1327  push2TmpBox_(iter2->getMZ(), scan_index, c, 0, getLinearInterpolation(iter2 - 1, help_mz, iter2), ref.getRT(), MZ_start, MZ_end);
1328  }
1329  }
1330 
1331  help_mz = seed_mz + Constants::IW_NEUTRON_MASS / (c + 1.);
1332  iter2 = candidates.MZBegin(help_mz);
1333  if (iter2 == candidates.end() || iter2 == candidates.begin())
1334  {
1335  continue;
1336  }
1337 
1338  if (fabs(iter2->getMZ() - seed_mz) > 0.5 * Constants::IW_NEUTRON_MASS / (c + 1.))
1339  {
1340  //In the other case, we are too close to the peak, leading to incorrect derivatives.
1341  if (iter2 != candidates.end())
1342  {
1343  push2TmpBox_(iter2->getMZ(), scan_index, c, 0, getLinearInterpolation(iter2 - 1, help_mz, iter2), ref.getRT(), MZ_start, MZ_end);
1344  }
1345  }
1346  }
1347 
1348  clusterSeeds_(candidates, ref, scan_index, c, check_PPMs);
1349  }
1350 
1351 #if defined(OPENMS_HAS_TBB) && defined(OPENMS_HAS_CUDA)
1352  template <typename PeakType>
1353  void IsotopeWaveletTransform<PeakType>::mergeFeatures(IsotopeWaveletTransform<PeakType>* later_iwt, const UInt RT_interleave, const UInt RT_votes_cutoff)
1354  {
1355  typename std::multimap<DoubleReal, Box>::iterator front_iter, end_iter, best_match, help_iter;
1356 
1357  //First of all do the trivial part of the merge
1358  for (end_iter = later_iwt->closed_boxes_.begin(); end_iter != later_iwt->closed_boxes_.end(); ++end_iter)
1359  {
1360  closed_boxes_.insert(*end_iter);
1361  }
1362 
1363  typename std::multimap<DoubleReal, Box>& end_container(this->end_boxes_);
1364  typename std::multimap<DoubleReal, Box>& front_container(later_iwt->front_boxes_);
1365 
1366 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1367  std::cout << "FontBox: " << front_container.size() << std::endl;
1368  for (front_iter = front_container.begin(); front_iter != front_container.end(); ++front_iter)
1369  std::cout << front_iter->first << "\t" << front_iter->second.size() << std::endl;
1370 
1371  std::cout << "EndBox: " << end_container.size() << std::endl;
1372  for (front_iter = end_container.begin(); front_iter != end_container.end(); ++front_iter)
1373  std::cout << front_iter->first << "\t" << front_iter->second.size() << std::endl;
1374 #endif
1375 
1376  typename std::multimap<UInt, BoxElement>::iterator biter;
1377 
1378  DoubleReal best_dist, c_dist; UInt c;
1379  //Now, try to find matching boxes for the rest
1380  for (front_iter = front_container.begin(); front_iter != front_container.end(); )
1381  {
1382  best_match = end_container.end(); best_dist = INT_MAX;
1383  //This is everything else than efficient, but both containers should be very small in size
1384  for (end_iter = end_container.begin(); end_iter != end_container.end(); ++end_iter)
1385  {
1386  c = 0;
1387  for (biter = front_iter->second.begin(); biter != front_iter->second.end(); ++biter)
1388  {
1389  c = std::max(c, biter->second.c);
1390  }
1391 
1392 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1393  std::cout << "Trying to match: " << end_iter->first << " to front " << front_iter->first << std::endl;
1394 #endif
1395 
1396  c_dist = fabs(end_iter->first - front_iter->first);
1397  if (c_dist < Constants::IW_HALF_NEUTRON_MASS / (c + 1.) && c_dist < best_dist)
1398  {
1399 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1400  std::cout << "best match " << front_iter->second.begin()->first << "\t" << (--(end_iter->second.end()))->first << std::endl;
1401 #endif
1402  if ((front_iter->second.begin()->first - (--(end_iter->second.end()))->first) <= RT_interleave + 1)
1403  {
1404  //otherwise, there are too many blank scans in between
1405  best_match = end_iter;
1406  best_dist = c_dist;
1407  }
1408  }
1409  }
1410  if (best_match == end_container.end()) //No matching pair found
1411  {
1412  if (front_iter->second.size() >= RT_votes_cutoff)
1413  {
1414  closed_boxes_.insert(*front_iter);
1415  //extendBox_ (map, front_iter->second);
1416  }
1417  ++front_iter;
1418  }
1419  else //That's the funny part
1420  {
1421 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1422  std::cout << "Merging the boxes: " << front_iter->first << "\t" << best_match->first << std::endl;
1423 #endif
1424 
1425  front_iter->second.insert(best_match->second.begin(), best_match->second.end());
1426  Box replacement(front_iter->second);
1427 
1428  //We cannot divide both m/z by 2, since we already inserted some m/zs whose weight would be lowered.
1429  DoubleReal c_mz = front_iter->first * (front_iter->second.size() - best_match->second.size()) + best_match->first * best_match->second.size();
1430  c_mz /= ((DoubleReal) (front_iter->second.size()));
1431 
1432  help_iter = front_iter;
1433  ++help_iter;
1434  std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
1435  closed_boxes_.insert(help3);
1436  //extendBox_ (map, help3.second);
1437  front_container.erase(front_iter);
1438  end_container.erase(best_match);
1439  front_iter = help_iter;
1440  }
1441  }
1442 
1443  //Merge the rest in end_container
1444  for (end_iter = end_container.begin(); end_iter != end_container.end(); ++end_iter)
1445  {
1446  if (end_iter->second.size() >= RT_votes_cutoff)
1447  {
1448  closed_boxes_.insert(*end_iter);
1449  //extendBox_ (map, end_iter->second);
1450  }
1451  }
1452  }
1453 
1454 #endif
1455 
1456 
1457  template <typename PeakType>
1459  const UInt peak_cutoff, const DoubleReal seed_mz, const UInt c, const DoubleReal ampl_cutoff)
1460  {
1461  DoubleReal c_score = 0, c_val;
1462  typename MSSpectrum<PeakType>::const_iterator c_left_iter2, c_right_iter2;
1463  Int signal_size((Int)candidate.size());
1464  //in the very unlikely case that size_t will not fit to int anymore this will be a problem of course
1465  //for the sake of simplicity (we need here a signed int) we do not cast at every following comparison individually
1466 
1467  //p_h_ind indicates if we are looking for a whole or a peak
1468  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
1469 
1470  std::vector<DoubleReal> positions(end);
1471  for (Int i = 0; i < end; ++i)
1472  {
1473  positions[i] = seed_mz - ((peak_cutoff - 1) * Constants::IW_NEUTRON_MASS - (i + 1) * Constants::IW_HALF_NEUTRON_MASS) / ((DoubleReal)c + 1);
1474  }
1475 
1476  DoubleReal l_score = 0, mid_val = 0;
1477  Int start_index = distance(candidate.begin(), candidate.MZBegin(positions[0])) - 1;
1478  for (Int v = 1; v <= end; ++v, ++p_h_ind)
1479  {
1480  do
1481  {
1482  if (start_index < signal_size - 1)
1483  ++start_index;
1484  else
1485  break;
1486  }
1487  while (candidate[start_index].getMZ() < positions[v - 1]);
1488 
1489  if (start_index <= 0 || start_index >= signal_size - 1) //unable to interpolate
1490  {
1491  continue;
1492  }
1493 
1494  c_left_iter2 = candidate.begin() + start_index - 1;
1495  c_right_iter2 = c_left_iter2 + 1;
1496 
1497  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());
1498 
1499  if (v == (int)(ceil(end / 2.)))
1500  {
1501  l_score = c_score;
1502  mid_val = c_val;
1503  }
1504 
1505  if (p_h_ind % 2 == 1) //I.e. a whole
1506  {
1507  c_score -= c_val;
1508 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1509  if (trunc(seed_mz) == 874)
1510  std::cout << -c_val << std::endl;
1511 #endif
1512  }
1513  else
1514  {
1515  c_score += c_val;
1516 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1517  if (trunc(seed_mz) == 874)
1518  std::cout << c_val << std::endl;
1519 #endif
1520  }
1521 
1522 
1523  start_index = distance(candidate.begin(), c_left_iter2);
1524  }
1525 
1526 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1527  std::ofstream ofile_score("scores.dat", ios::app);
1528  std::ofstream ofile_check_score("check_scores.dat", ios::app);
1529  ofile_score.close();
1530  ofile_check_score.close();
1531 #endif
1532 
1533 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1534  if (trunc(seed_mz) == 874)
1535  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;
1536 #endif
1537 
1538  if (c_score - mid_val <= 0)
1539  {
1540  return 0;
1541  }
1542 
1543  if (c_score - mid_val <= ampl_cutoff)
1544  {
1545  return -1000;
1546  }
1547 
1548  if (l_score <= 0 || c_score - l_score - mid_val <= 0)
1549  {
1550  return 0;
1551  }
1552 
1553  return c_score;
1554  }
1555 
1556  template <typename PeakType>
1558  const UInt peak_cutoff, const DoubleReal seed_mz, const UInt c, const DoubleReal ampl_cutoff)
1559  {
1560  DoubleReal c_score = 0, c_val;
1561  typename MSSpectrum<PeakType>::const_iterator c_left_iter2, c_right_iter2;
1562  Int signal_size((Int)candidate.size());
1563 
1564  //p_h_ind indicates if we are looking for a whole or a peak
1565  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
1566 
1567  std::vector<DoubleReal> positions(end);
1568  for (Int i = 0; i < end; ++i)
1569  {
1570  positions[i] = seed_mz - ((peak_cutoff - 1) * Constants::IW_NEUTRON_MASS - (i + 1) * Constants::IW_HALF_NEUTRON_MASS) / ((DoubleReal)c + 1);
1571  }
1572 
1573  DoubleReal l_score = 0, mid_val = 0;
1574  Int start_index = distance(candidate.begin(), candidate.MZBegin(positions[0])) - 1;
1575  for (Int v = 1; v <= end; ++v, ++p_h_ind)
1576  {
1577  do
1578  {
1579  if (start_index < signal_size - 1)
1580  ++start_index;
1581  else
1582  break;
1583  }
1584  while (candidate.getMZ(start_index) < positions[v - 1]);
1585 
1586  if (start_index <= 0 || start_index >= signal_size - 1) //unable to interpolate
1587  {
1588  continue;
1589  }
1590 
1591  c_left_iter2 = candidate.begin() + start_index - 1;
1592  c_right_iter2 = c_left_iter2 + 1;
1593 
1594  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());
1595  if (v == (int)(ceil(end / 2.)))
1596  {
1597  l_score = c_score;
1598  mid_val = c_val;
1599  }
1600 
1601  if (p_h_ind % 2 == 1) //I.e. a whole
1602  {
1603  c_score -= c_val;
1604  }
1605  else
1606  {
1607  c_score += c_val;
1608  }
1609 
1610  start_index = distance(candidate.begin(), c_left_iter2);
1611  }
1612 
1613 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1614  std::ofstream ofile_score("scores.dat", ios::app);
1615  std::ofstream ofile_check_score("check_scores.dat", ios::app);
1616  ofile_score << c_check_point << "\t" << c_score << std::endl;
1617  ofile_score.close();
1618  ofile_check_score.close();
1619 #endif
1620 
1621  if (l_score <= 0 || c_score - l_score - mid_val <= 0 || c_score - mid_val <= ampl_cutoff)
1622  {
1623  return 0;
1624  }
1625 
1626  return c_score;
1627  }
1628 
1629  template <typename PeakType>
1631  const MSSpectrum<PeakType>& ref, const UInt scan_index, const UInt c, const bool check_PPMs)
1632  {
1633  typename std::multimap<DoubleReal, Box>::iterator iter;
1634  typename Box::iterator box_iter;
1635  std::vector<BoxElement> final_box;
1636  DoubleReal c_mz, av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
1637  DoubleReal virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
1638 
1639  typename std::pair<DoubleReal, DoubleReal> c_extend;
1640  for (iter = tmp_boxes_->at(c).begin(); iter != tmp_boxes_->at(c).end(); ++iter)
1641  {
1642 
1643  Box& c_box = iter->second;
1644  av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
1645  virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
1646 
1647  //Now, let's get the RT boundaries for the box
1648  for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1649  {
1650  if (box_iter->second.score == 0) //virtual helping point
1651  {
1652  if (count != 0)
1653  continue; //it is in any way not pure virtual
1654 
1655  c_mz = box_iter->second.mz;
1656  virtual_av_intens += box_iter->second.intens;
1657  virtual_av_abs_intens += fabs(box_iter->second.intens);
1658  virtual_av_mz += c_mz * fabs(box_iter->second.intens);
1659  ++virtual_count;
1660  }
1661  else
1662  {
1663  c_mz = box_iter->second.mz;
1664  av_score += box_iter->second.score;
1665  av_intens += box_iter->second.intens;
1666  av_abs_intens += fabs(box_iter->second.intens);
1667  av_mz += c_mz * fabs(box_iter->second.intens);
1668  ++count;
1669  }
1670  }
1671 
1672  if (count == 0) //pure virtual helping box
1673  {
1674  av_intens = virtual_av_intens / virtual_count;
1675  av_score = 0;
1676  av_mz = virtual_av_mz / virtual_av_abs_intens;
1677  }
1678  else
1679  {
1680  av_intens /= count;
1681  av_score /= count;
1682  av_mz /= av_abs_intens;
1683  }
1684 
1685  BoxElement c_box_element;
1686  c_box_element.mz = av_mz;
1687  c_box_element.c = c;
1688  c_box_element.score = av_score;
1689  c_box_element.intens = av_intens;
1690 
1691  c_box_element.RT = c_box.begin()->second.RT;
1692  final_box.push_back(c_box_element);
1693  }
1694 
1695  Size num_o_feature = final_box.size();
1696  if (num_o_feature == 0)
1697  {
1698  tmp_boxes_->at(c).clear();
1699  return;
1700  }
1701 
1702  //Computing the derivatives
1703  std::vector<DoubleReal> bwd_diffs(num_o_feature, 0);
1704 
1705  bwd_diffs[0] = 0;
1706  for (Size i = 1; i < num_o_feature; ++i)
1707  {
1708  bwd_diffs[i] = (final_box[i].intens - final_box[i - 1].intens) / (final_box[i].mz - final_box[i - 1].mz);
1709  }
1710 
1711 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1712  std::ofstream ofile_bwd("bwd_cpu.dat");
1713  for (Size i = 0; i < num_o_feature; ++i)
1714  {
1715  ofile_bwd << final_box[i].mz << "\t" << bwd_diffs[i] << std::endl;
1716  }
1717  ofile_bwd.close();
1718 #endif
1719 
1720 
1721  for (Size i = 0; i < num_o_feature - 1; ++i)
1722  {
1723  while (i < num_o_feature - 2)
1724  {
1725  if (final_box[i].score > 0 || final_box[i].score == -1000) //this has been an helping point
1726  break;
1727  ++i;
1728  }
1729 
1730  if (bwd_diffs[i] > 0 && bwd_diffs[i + 1] < 0)
1731  {
1732  checkPositionForPlausibility_(candidate, ref, final_box[i].mz, final_box[i].c, scan_index, check_PPMs, final_box[i].intens, final_box[i].score);
1733  continue;
1734  }
1735  }
1736 
1737  tmp_boxes_->at(c).clear();
1738  }
1739 
1740  template <typename PeakType>
1742  {
1743  DoubleReal av_intens = 0;
1744  for (UInt i = 0; i < scan.size(); ++i)
1745  {
1746  if (scan[i].getIntensity() >= 0)
1747  {
1748  av_intens += scan[i].getIntensity();
1749  }
1750  }
1751  return av_intens / (double)scan.size();
1752  }
1753 
1754  template <typename PeakType>
1756  {
1757  DoubleReal res = 0, intens;
1758  for (UInt i = 0; i < scan.size(); ++i)
1759  {
1760  if (scan[i].getIntensity() >= 0)
1761  {
1762  intens = scan[i].getIntensity();
1763  res += (intens - mean) * (intens - mean);
1764  }
1765  }
1766  return sqrt(res / (double)(scan.size() - 1));
1767  }
1768 
1769  template <typename PeakType>
1771  {
1772  std::vector<DoubleReal> diffs(scan.size() - 1, 0);
1773  for (UInt i = 0; i < scan.size() - 1; ++i)
1774  {
1775  diffs[i] = scan[i + 1].getMZ() - scan[i].getMZ();
1776  }
1777 
1778  sort(diffs.begin(), diffs.end());
1779  DoubleReal av_MZ_spacing = 0;
1780  for (UInt i = 0; i < diffs.size() / 2; ++i)
1781  {
1782  av_MZ_spacing += diffs[i];
1783  }
1784 
1785  return av_MZ_spacing / (diffs.size() / 2);
1786  }
1787 
1788  template <typename PeakType>
1790  {
1791  DoubleReal av_intens = 0;
1792  for (UInt i = 0; i < scan.size(); ++i)
1793  {
1794  if (scan.getTransIntensity(i) >= 0)
1795  {
1796  av_intens += scan.getTransIntensity(i);
1797  }
1798  }
1799  return av_intens / (double)scan.size();
1800  }
1801 
1802  template <typename PeakType>
1804  {
1805  DoubleReal res = 0, intens;
1806  for (UInt i = 0; i < scan.size(); ++i)
1807  {
1808  if (scan.getTransIntensity(i) >= 0)
1809  {
1810  intens = scan.getTransIntensity(i);
1811  res += (intens - mean) * (intens - mean);
1812  }
1813  }
1814  return sqrt(res / (double)(scan.size() - 1));
1815  }
1816 
1817  template <typename PeakType>
1819  const DoubleReal score, const DoubleReal intens, const DoubleReal rt, const UInt MZ_begin, const UInt MZ_end, DoubleReal ref_intens)
1820  {
1821  const DoubleReal dist_constraint(Constants::IW_HALF_NEUTRON_MASS / (DoubleReal)max_charge_);
1822 
1823  typename std::multimap<DoubleReal, Box>::iterator upper_iter(open_boxes_.upper_bound(mz));
1824  typename std::multimap<DoubleReal, Box>::iterator lower_iter(open_boxes_.lower_bound(mz));
1825 
1826  if (lower_iter != open_boxes_.end())
1827  {
1828  //Ugly, but necessary due to the implementation of STL lower_bound
1829  if (mz != lower_iter->first && lower_iter != open_boxes_.begin())
1830  {
1831  --lower_iter;
1832  }
1833  }
1834 
1835  typename std::multimap<DoubleReal, Box>::iterator insert_iter;
1836  bool create_new_box = true;
1837  if (lower_iter == open_boxes_.end()) //I.e. there is no open Box for that mz position
1838  {
1839  //There is another special case to be considered here:
1840  //Assume that the current box contains only a single element that is (slightly) smaller than the new mz value,
1841  //then the lower bound for the new mz value is box.end and this would usually force a new entry
1842  if (!open_boxes_.empty())
1843  {
1844  if (fabs((--lower_iter)->first - mz) < dist_constraint) //matching box
1845  {
1846  create_new_box = false;
1847  insert_iter = lower_iter;
1848  }
1849  }
1850  else
1851  {
1852  create_new_box = true;
1853  }
1854  }
1855  else
1856  {
1857  if (upper_iter == open_boxes_.end() && fabs(lower_iter->first - mz) < dist_constraint) //Found matching Box
1858  {
1859  insert_iter = lower_iter;
1860  create_new_box = false;
1861  }
1862  else
1863  {
1864  create_new_box = true;
1865  }
1866  }
1867 
1868 
1869  if (upper_iter != open_boxes_.end() && lower_iter != open_boxes_.end())
1870  {
1871  //Here is the question if you should figure out the smallest charge .... and then
1872 
1873  //Figure out which entry is closer to m/z
1874  DoubleReal dist_lower = fabs(lower_iter->first - mz);
1875  DoubleReal dist_upper = fabs(upper_iter->first - mz);
1876  dist_lower = (dist_lower < dist_constraint) ? dist_lower : INT_MAX;
1877  dist_upper = (dist_upper < dist_constraint) ? dist_upper : INT_MAX;
1878 
1879  if (dist_lower >= dist_constraint && dist_upper >= dist_constraint) // they are both too far away
1880  {
1881  create_new_box = true;
1882  }
1883  else
1884  {
1885  insert_iter = (dist_lower < dist_upper) ? lower_iter : upper_iter;
1886  create_new_box = false;
1887  }
1888  }
1889 
1890  BoxElement element;
1891  element.c = c; element.mz = mz; element.score = score; element.RT = rt; element.intens = intens; element.ref_intens = ref_intens;
1892  element.RT_index = scan; element.MZ_begin = MZ_begin; element.MZ_end = MZ_end;
1893 
1894 
1895  if (create_new_box == false)
1896  {
1897  std::pair<UInt, BoxElement> help2(scan, element);
1898  insert_iter->second.insert(help2);
1899 
1900  //Unfortunately, we need to change the m/z key to the average of all keys inserted in that box.
1901  Box replacement(insert_iter->second);
1902 
1903  //We cannot divide both m/z by 2, since we already inserted some m/zs whose weight would be lowered.
1904  //Also note that we already inserted the new entry, leading to size-1.
1905  DoubleReal c_mz = insert_iter->first * (insert_iter->second.size() - 1) + mz;
1906  c_mz /= ((DoubleReal) insert_iter->second.size());
1907 
1908  //Now let's remove the old and insert the new one
1909  open_boxes_.erase(insert_iter);
1910  std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
1911  open_boxes_.insert(help3);
1912  }
1913  else
1914  {
1915  std::pair<UInt, BoxElement> help2(scan, element);
1916  std::multimap<UInt, BoxElement> help3;
1917  help3.insert(help2);
1918  std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help4(mz, help3);
1919  open_boxes_.insert(help4);
1920  }
1921  }
1922 
1923  template <typename PeakType>
1925  const DoubleReal score, const DoubleReal intens, const DoubleReal rt, const UInt MZ_begin, const UInt MZ_end)
1926  {
1927  const DoubleReal dist_constraint(Constants::IW_HALF_NEUTRON_MASS / (DoubleReal)max_charge_);
1928 
1929  std::multimap<DoubleReal, Box>& tmp_box(tmp_boxes_->at(c));
1930  typename std::multimap<DoubleReal, Box>::iterator upper_iter(tmp_box.upper_bound(mz));
1931  typename std::multimap<DoubleReal, Box>::iterator lower_iter(tmp_box.lower_bound(mz));
1932 
1933  if (lower_iter != tmp_box.end())
1934  {
1935  //Ugly, but necessary due to the implementation of STL lower_bound
1936  if (mz != lower_iter->first && lower_iter != tmp_box.begin())
1937  {
1938  --lower_iter;
1939  }
1940  }
1941 
1942  typename std::multimap<DoubleReal, Box>::iterator insert_iter;
1943  bool create_new_box = true;
1944  if (lower_iter == tmp_box.end()) //I.e. there is no tmp Box for that mz position
1945  {
1946  //There is another special case to be considered here:
1947  //Assume that the current box contains only a single element that is (slightly) smaller than the new mz value,
1948  //then the lower bound for the new mz value is box.end and this would usually force a new entry
1949  if (!tmp_box.empty())
1950  {
1951  if (fabs((--lower_iter)->first - mz) < dist_constraint) //matching box
1952  {
1953  create_new_box = false;
1954  insert_iter = lower_iter;
1955  }
1956  }
1957  else
1958  {
1959  create_new_box = true;
1960  }
1961  }
1962  else
1963  {
1964  if (upper_iter == tmp_box.end() && fabs(lower_iter->first - mz) < dist_constraint) //Found matching Box
1965  {
1966  insert_iter = lower_iter;
1967  create_new_box = false;
1968  }
1969  else
1970  {
1971  create_new_box = true;
1972  }
1973  }
1974 
1975 
1976  if (upper_iter != tmp_box.end() && lower_iter != tmp_box.end())
1977  {
1978  //Figure out which entry is closer to m/z
1979  DoubleReal dist_lower = fabs(lower_iter->first - mz);
1980  DoubleReal dist_upper = fabs(upper_iter->first - mz);
1981  dist_lower = (dist_lower < dist_constraint) ? dist_lower : INT_MAX;
1982  dist_upper = (dist_upper < dist_constraint) ? dist_upper : INT_MAX;
1983 
1984  if (dist_lower >= dist_constraint && dist_upper >= dist_constraint) // they are both too far away
1985  {
1986  create_new_box = true;
1987  }
1988  else
1989  {
1990  insert_iter = (dist_lower < dist_upper) ? lower_iter : upper_iter;
1991  create_new_box = false;
1992  }
1993  }
1994 
1995  BoxElement element;
1996  element.c = c; element.mz = mz; element.score = score; element.RT = rt; element.intens = intens; element.ref_intens = -1000;
1997  element.RT_index = scan; element.MZ_begin = MZ_begin; element.MZ_end = MZ_end;
1998 
1999  if (create_new_box == false)
2000  {
2001  std::pair<UInt, BoxElement> help2(scan, element);
2002  insert_iter->second.insert(help2);
2003 
2004  //Unfortunately, we need to change the m/z key to the average of all keys inserted in that box.
2005  Box replacement(insert_iter->second);
2006 
2007  //We cannot divide both m/z by 2, since we already inserted some m/zs whose weight would be lowered.
2008  //Also note that we already inserted the new entry, leading to size-1.
2009  DoubleReal c_mz = insert_iter->first * (insert_iter->second.size() - 1) + mz;
2010  c_mz /= ((DoubleReal) insert_iter->second.size());
2011 
2012  //Now let's remove the old and insert the new one
2013  tmp_box.erase(insert_iter);
2014  std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
2015  tmp_box.insert(help3);
2016  }
2017  else
2018  {
2019  std::pair<UInt, BoxElement> help2(scan, element);
2020  std::multimap<UInt, BoxElement> help3;
2021  help3.insert(help2);
2022 
2023  std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help4(mz, help3);
2024  tmp_box.insert(help4);
2025  }
2026  }
2027 
2028  template <typename PeakType>
2029  void IsotopeWaveletTransform<PeakType>::updateBoxStates(const MSExperiment<PeakType>& map, const Size scan_index, const UInt RT_interleave,
2030  const UInt RT_votes_cutoff, const Int front_bound, const Int end_bound)
2031  {
2032  typename std::multimap<DoubleReal, Box>::iterator iter, iter2;
2033 
2034  if ((Int)scan_index == end_bound && end_bound != (Int)map.size() - 1)
2035  {
2036  for (iter = open_boxes_.begin(); iter != open_boxes_.end(); ++iter)
2037  {
2038 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2039  std::cout << "LOW THREAD insert in end_box " << iter->first << std::endl;
2040  typename Box::iterator dings;
2041  for (dings = iter->second.begin(); dings != iter->second.end(); ++dings)
2042  std::cout << map[dings->first].getRT() << "\t" << dings->second.c + 1 << std::endl;
2043 #endif
2044  end_boxes_.insert(*iter);
2045  }
2046  open_boxes_.clear();
2047  return;
2048  }
2049 
2050  for (iter = open_boxes_.begin(); iter != open_boxes_.end(); )
2051  {
2052 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2053  if (front_bound > 0)
2054  {
2055  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;
2056  }
2057 #endif
2058 
2059  //For each Box we need to figure out, if and when the last RT value has been inserted
2060  UInt lastScan = (--(iter->second.end()))->first;
2061  if (scan_index - lastScan > RT_interleave + 1 || scan_index == map.size() - 1) //I.e. close the box!
2062  {
2063  if (iter->second.begin()->first - front_bound <= RT_interleave + 1 && front_bound > 0)
2064  {
2065  iter2 = iter;
2066  ++iter2;
2067 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2068  std::cout << "HIGH THREAD insert in front_box " << iter->first << std::endl;
2069 #endif
2070  front_boxes_.insert(*iter);
2071  open_boxes_.erase(iter);
2072  iter = iter2;
2073  continue;
2074  }
2075 
2076  iter2 = iter;
2077  ++iter2;
2078  //Please do **NOT** simplify the upcoming lines.
2079  //The 'obvious' overhead is necessary since the object represented by iter might be erased
2080  //by push2Box which might be called by extendBox_.
2081  if (iter->second.size() >= RT_votes_cutoff)
2082  {
2083  //extendBox_ (map, iter->second);
2084  iter = iter2;
2085  closed_boxes_.insert(*(--iter));
2086  }
2087  open_boxes_.erase(iter);
2088  iter = iter2;
2089  }
2090  else
2091  {
2092  ++iter;
2093  }
2094  }
2095  }
2096 
2097  template <typename PeakType>
2099  {
2100 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2101  std::cout << "**** CHECKING FOR BOX EXTENSIONS ****" << std::endl;
2102 #endif
2103 
2104  //Determining the elution profile
2105  typename Box::const_iterator iter;
2106  std::vector<DoubleReal> elution_profile(box.size());
2107  UInt index = 0;
2108  for (iter = box.begin(); iter != box.end(); ++iter, ++index)
2109  {
2110  for (Size i = iter->second.MZ_begin; i != iter->second.MZ_end; ++i)
2111  {
2112  elution_profile[index] += map[iter->second.RT_index][i].getIntensity();
2113  }
2114  elution_profile[index] /= iter->second.MZ_end - iter->second.MZ_begin + 1.;
2115  }
2116 
2117  DoubleReal max = 0;
2118  Int max_index = INT_MIN;
2119  for (Size i = 0; i < elution_profile.size(); ++i)
2120  {
2121  if (elution_profile[i] > max)
2122  {
2123  max_index = i;
2124  max = elution_profile[i];
2125  }
2126  }
2127 
2128  Int max_extension = (Int)(elution_profile.size()) - 2 * max_index;
2129 
2130  DoubleReal av_elution = 0;
2131  for (Size i = 0; i < elution_profile.size(); ++i)
2132  {
2133  av_elution += elution_profile[i];
2134  }
2135  av_elution /= (DoubleReal)elution_profile.size();
2136 
2137  DoubleReal sd_elution = 0;
2138  for (Size i = 0; i < elution_profile.size(); ++i)
2139  {
2140  sd_elution += (av_elution - elution_profile[i]) * (av_elution - elution_profile[i]);
2141  }
2142  sd_elution /= (DoubleReal)(elution_profile.size() - 1);
2143  sd_elution = sqrt(sd_elution);
2144 
2145  //Determine average m/z monoisotopic pos
2146  DoubleReal av_mz = 0;
2147  for (iter = box.begin(); iter != box.end(); ++iter, ++index)
2148  {
2149  av_mz += iter->second.mz;
2150 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2151  std::cout << iter->second.RT << "\t" << iter->second.mz << "\t" << iter->second.c + 1 << std::endl;
2152 #endif
2153  }
2154  av_mz /= (DoubleReal)box.size();
2155 
2156 
2157  //Boundary check
2158  if ((Int)(box.begin()->second.RT_index) - 1 < 0)
2159  {
2160  return;
2161  }
2162 
2163  UInt pre_index = box.begin()->second.RT_index - 1;
2164  typename MSSpectrum<PeakType>::const_iterator c_iter = map[pre_index].MZBegin(av_mz);
2165  DoubleReal pre_elution = 0;
2166 
2167  DoubleReal mz_start = map[pre_index + 1][box.begin()->second.MZ_begin].getMZ();
2168  DoubleReal mz_end = map[pre_index + 1][box.begin()->second.MZ_end].getMZ();
2169 
2170  typename MSSpectrum<PeakType>::const_iterator mz_start_iter = map[pre_index].MZBegin(mz_start), mz_end_iter = map[pre_index].MZBegin(mz_end);
2171  for (typename MSSpectrum<PeakType>::const_iterator mz_iter = mz_start_iter; mz_iter != mz_end_iter; ++mz_iter)
2172  {
2173  pre_elution += mz_iter->getIntensity();
2174  }
2175 
2176 
2177  //Do we need to extend at all?
2178  if (pre_elution <= av_elution - 2 * sd_elution)
2179  {
2180  return;
2181  }
2182 
2183  Int c_index = max_extension;
2184  Int first_index = box.begin()->second.RT_index;
2185  for (Int i = 1; i < max_extension; ++i)
2186  {
2187  c_index = first_index - i;
2188  if (c_index < 0)
2189  {
2190  break;
2191  }
2192 
2193  //CHECK Majority vote for charge???????????????
2194 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2195  std::cout << box.begin()->second.RT << "\t" << av_mz << "\t" << box.begin()->second.c + 1 << "\t" << " extending the box " << std::endl;
2196 #endif
2197 
2198  push2Box_(av_mz, c_index, box.begin()->second.c, box.begin()->second.score, c_iter->getIntensity(),
2199  map[c_index].getRT(), box.begin()->second.MZ_begin, box.begin()->second.MZ_end);
2200  }
2201  }
2202 
2203  template <typename PeakType>
2205  const MSSpectrum<PeakType>& ref, const UInt scan_index, const UInt c, const bool check_PPMs)
2206  {
2207  typename std::multimap<DoubleReal, Box>::iterator iter;
2208  typename Box::iterator box_iter;
2209  std::vector<BoxElement> final_box;
2210  DoubleReal c_mz, av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
2211  DoubleReal virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
2212 
2213  typename std::pair<DoubleReal, DoubleReal> c_extend;
2214  for (iter = tmp_boxes_->at(c).begin(); iter != tmp_boxes_->at(c).end(); ++iter)
2215  {
2216  Box& c_box = iter->second;
2217  av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
2218  virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
2219 
2220  for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
2221  {
2222  if (box_iter->second.score == 0) //virtual helping point
2223  {
2224  if (count != 0)
2225  continue; //it is in any way not pure virtual
2226 
2227  c_mz = box_iter->second.mz;
2228  virtual_av_intens += box_iter->second.intens;
2229  virtual_av_abs_intens += fabs(box_iter->second.intens);
2230  virtual_av_mz += c_mz * fabs(box_iter->second.intens);
2231  ++virtual_count;
2232  }
2233  else
2234  {
2235  c_mz = box_iter->second.mz;
2236  av_score += box_iter->second.score;
2237  av_intens += box_iter->second.intens;
2238  av_abs_intens += fabs(box_iter->second.intens);
2239  av_mz += c_mz * fabs(box_iter->second.intens);
2240  ++count;
2241  }
2242  }
2243 
2244  if (count == 0) //pure virtual helping box
2245  {
2246  av_intens = virtual_av_intens / virtual_count;
2247  av_score = 0;
2248  av_mz = virtual_av_mz / virtual_av_abs_intens;
2249  }
2250  else
2251  {
2252  av_intens /= count;
2253  av_score /= count;
2254  av_mz /= av_abs_intens;
2255  }
2256 
2257  BoxElement c_box_element;
2258  c_box_element.mz = av_mz;
2259  c_box_element.c = c;
2260  c_box_element.score = av_score;
2261  c_box_element.intens = av_intens;
2262 
2263  c_box_element.RT = c_box.begin()->second.RT;
2264 
2265  final_box.push_back(c_box_element);
2266  }
2267 
2268  UInt num_o_feature = final_box.size();
2269  if (num_o_feature == 0)
2270  {
2271  tmp_boxes_->at(c).clear();
2272  return;
2273  }
2274 
2275  //Computing the derivatives
2276  std::vector<DoubleReal> bwd_diffs(num_o_feature, 0);
2277 
2278  bwd_diffs[0] = 0;
2279  for (UInt i = 1; i < num_o_feature; ++i)
2280  {
2281  bwd_diffs[i] = (final_box[i].intens - final_box[i - 1].intens) / (final_box[i].mz - final_box[i - 1].mz);
2282  }
2283 
2284 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2285  std::ofstream ofile_bwd("bwd_gpu.dat");
2286  for (UInt i = 0; i < num_o_feature; ++i)
2287  {
2288  ofile_bwd << final_box[i].mz << "\t" << bwd_diffs[i] << std::endl;
2289  }
2290  ofile_bwd.close();
2291 #endif
2292 
2293  for (UInt i = 0; i < num_o_feature - 1; ++i)
2294  {
2295  while (i < num_o_feature - 2)
2296  {
2297  if (final_box[i].score > 0 || final_box[i].score == -1000) //this has been an helping point
2298  break;
2299  ++i;
2300  }
2301 
2302  if (bwd_diffs[i] > 0 && bwd_diffs[i + 1] < 0)
2303  {
2304  checkPositionForPlausibility_(candidates, ref, final_box[i].mz, final_box[i].c, scan_index, check_PPMs, final_box[i].intens, final_box[i].score);
2305  continue;
2306  }
2307  }
2308 
2309  tmp_boxes_->at(c).clear();
2310  }
2311 
2312  template <typename PeakType>
2314  {
2315  FeatureMap<Feature> feature_map;
2316  typename std::multimap<DoubleReal, Box>::iterator iter;
2317  typename Box::iterator box_iter;
2318  UInt best_charge_index; DoubleReal best_charge_score, c_mz, c_RT; UInt c_charge;
2319  DoubleReal av_intens = 0, av_ref_intens = 0, av_score = 0, av_mz = 0, av_RT = 0, mz_cutoff, sum_of_ref_intenses_g;
2320  bool restart = false;
2321 
2322  typename std::pair<DoubleReal, DoubleReal> c_extend;
2323  for (iter = closed_boxes_.begin(); iter != closed_boxes_.end(); ++iter)
2324  {
2325  sum_of_ref_intenses_g = 0;
2326  Box& c_box = iter->second;
2327  std::vector<DoubleReal> charge_votes(max_charge_, 0), charge_binary_votes(max_charge_, 0);
2328  restart = false;
2329 
2330  //Let's first determine the charge
2331  //Therefor, we can use two types of votes: qualitative ones (charge_binary_votes) or quantitative ones (charge_votes)
2332  for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
2333  {
2334 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2335  if (trunc(box_iter->second.mz) == 874)
2336  std::cout << box_iter->second.c << "\t" << box_iter->second.intens << "\t" << box_iter->second.score << std::endl;
2337 #endif
2338 
2339  if (box_iter->second.score == -1000)
2340  {
2341  restart = true;
2342  break;
2343  }
2344 
2345  charge_votes[box_iter->second.c] += box_iter->second.intens; //score; Do not use score, can get problematic for charge state 2 vs 4
2346  ++charge_binary_votes[box_iter->second.c];
2347  }
2348 
2349  if (restart)
2350  {
2351  continue;
2352  }
2353 
2354  //... determining the best fitting charge
2355  best_charge_index = 0; best_charge_score = 0;
2356  for (UInt i = 0; i < max_charge_; ++i)
2357  {
2358  if (charge_votes[i] > best_charge_score)
2359  {
2360  best_charge_index = i;
2361  best_charge_score = charge_votes[i];
2362  }
2363  }
2364 
2365  //Pattern found in too few RT scan
2366  if (charge_binary_votes[best_charge_index] < RT_votes_cutoff && RT_votes_cutoff <= map.size())
2367  {
2368  continue;
2369  }
2370 
2371  c_charge = best_charge_index + 1; //that's the finally predicted charge state for the pattern
2372 
2373  av_intens = 0, av_ref_intens = 0, av_score = 0, av_mz = 0, av_RT = 0;
2374  //Now, let's get the RT boundaries for the box
2375  std::vector<DPosition<2> > point_set;
2376  DoubleReal sum_of_ref_intenses_l;
2377  for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
2378  {
2379  sum_of_ref_intenses_l = 0;
2380  c_mz = box_iter->second.mz;
2381  c_RT = box_iter->second.RT;
2382 
2383  mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(c_mz, c_charge);
2384 
2385  point_set.push_back(DPosition<2>(c_RT, c_mz - Constants::IW_QUARTER_NEUTRON_MASS / (DoubleReal)c_charge));
2386  //-1 since we are already at the first peak and +0.75, since this includes the last peak of the wavelet as a whole
2387  point_set.push_back(DPosition<2>(c_RT, c_mz + mz_cutoff / (DoubleReal)c_charge));
2388 
2389 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2390  std::cout << "Intenstype: " << intenstype_ << std::endl;
2391 #endif
2392  if (intenstype_ == "ref")
2393  {
2394  //Find monoisotopic max
2395  const MSSpectrum<PeakType>& c_spec(map[box_iter->second.RT_index]);
2396  //'Correct' possible shift
2397  for (unsigned int i = 0; i < mz_cutoff; ++i)
2398  {
2399  typename MSSpectrum<PeakType>::const_iterator h_iter = c_spec.MZBegin(c_mz + i * Constants::IW_NEUTRON_MASS / c_charge + Constants::IW_QUARTER_NEUTRON_MASS / (DoubleReal)c_charge), hc_iter = c_spec.MZBegin(c_mz + i * Constants::IW_NEUTRON_MASS / c_charge);
2400 
2401  hc_iter = c_spec.MZBegin(c_mz + i * Constants::IW_NEUTRON_MASS / c_charge);
2402 
2403  while (h_iter != c_spec.begin())
2404  {
2405 
2406 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2407  if (trunc(c_mz) == 874)
2408  {
2409  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;
2410  }
2411 #endif
2412 
2413  --h_iter;
2414  if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2415  {
2416  hc_iter = h_iter;
2417  }
2418 
2419  if (c_mz + i * Constants::IW_NEUTRON_MASS / c_charge - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (DoubleReal)c_charge)
2420  {
2421  break;
2422  }
2423  }
2424 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2425  if (trunc(c_mz) == 874)
2426  {
2427  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";
2428  }
2429 #endif
2430  sum_of_ref_intenses_l += hc_iter->getIntensity();
2431 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2432  if (trunc(c_mz) == 874)
2433  {
2434  std::cout << sum_of_ref_intenses_l << "********" << std::endl;
2435  }
2436 #endif
2437  }
2438  }
2439 
2440  if (best_charge_index == box_iter->second.c)
2441  {
2442  av_score += box_iter->second.score;
2443  av_intens += box_iter->second.intens;
2444  av_ref_intens += box_iter->second.ref_intens;
2445  sum_of_ref_intenses_g += sum_of_ref_intenses_l;
2446  av_mz += c_mz * box_iter->second.intens;
2447  }
2448  av_RT += c_RT;
2449  }
2450 
2451  av_mz /= av_intens;
2452  av_ref_intens /= (DoubleReal)charge_binary_votes[best_charge_index];
2453  av_score /= (DoubleReal)charge_binary_votes[best_charge_index];
2454  av_RT /= (DoubleReal)c_box.size();
2455 
2456 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2457  if (trunc(av_mz) == 874)
2458  std::cout << av_mz << "\t" << best_charge_index << "\t" << best_charge_score << std::endl;
2459 #endif
2460 
2461  Feature c_feature;
2462  ConvexHull2D c_conv_hull;
2463  c_conv_hull.addPoints(point_set);
2464  c_feature.setCharge(c_charge);
2465  c_feature.setConvexHulls(std::vector<ConvexHull2D>(1, c_conv_hull));
2466 
2467  //This makes the intensity value independent of the m/z (the lambda) value (Skellam distribution)
2468  if (intenstype_ == "corrected")
2469  {
2470  DoubleReal lambda = IsotopeWavelet::getLambdaL(av_mz * c_charge);
2471  av_intens /= exp(-2 * lambda) * boost::math::cyl_bessel_i(0, 2 * lambda);
2472  }
2473  if (intenstype_ == "ref")
2474  {
2475 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2476  if (trunc(c_mz) == 874)
2477  {
2478  std::cout << sum_of_ref_intenses_g << "####" << std::endl;
2479  }
2480 #endif
2481 
2482  av_intens = sum_of_ref_intenses_g;
2483  }
2484 
2485  c_feature.setMZ(av_mz);
2486  c_feature.setIntensity(av_intens);
2487  c_feature.setRT(av_RT);
2488  c_feature.setOverallQuality(av_score);
2489  feature_map.push_back(c_feature);
2490  }
2491 
2492  return feature_map;
2493  }
2494 
2495  template <typename PeakType>
2497  const MSSpectrum<PeakType>& ref, const DoubleReal seed_mz, const UInt c, const UInt scan_index, const bool check_PPMs, const DoubleReal transintens, const DoubleReal prev_score)
2498  {
2499  typename MSSpectrum<PeakType>::const_iterator iter, ref_iter;
2500  UInt peak_cutoff;
2501  peak_cutoff = IsotopeWavelet::getNumPeakCutOff(seed_mz, c + 1);
2502 
2503  iter = candidate.MZBegin(seed_mz);
2504  //we can ignore those cases
2505  if (iter == candidate.begin() || iter == candidate.end())
2506  {
2507  return false;
2508  }
2509 
2510  std::pair<DoubleReal, DoubleReal> reals;
2511  ref_iter = ref.MZBegin(seed_mz);
2512  //Correct the position
2513  DoubleReal real_mz, real_intens;
2514  if (check_PPMs)
2515  {
2516  reals = checkPPMTheoModel_(ref, iter->getMZ(), c);
2517  real_mz = reals.first, real_intens = reals.second;
2518  //if (real_mz <= 0 || real_intens <= 0)
2519  //{
2520  typename MSSpectrum<PeakType>::const_iterator h_iter = ref_iter, hc_iter = ref_iter;
2521  while (h_iter != ref.begin())
2522  {
2523  --h_iter;
2524  if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2525  {
2526  hc_iter = h_iter;
2527  }
2528  else
2529  {
2530  break;
2531  }
2532 
2533  if (seed_mz - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.))
2534  {
2535  return false;
2536  }
2537  }
2538  reals = checkPPMTheoModel_(ref, h_iter->getMZ(), c);
2539  real_mz = reals.first, real_intens = reals.second;
2540  if (real_mz <= 0 || real_intens <= 0)
2541  {
2542  return false;
2543  }
2544  real_mz = h_iter->getMZ();
2545  real_intens = h_iter->getIntensity();
2546  //}
2547  }
2548  else
2549  {
2550  reals = std::pair<DoubleReal, DoubleReal>(seed_mz, ref_iter->getIntensity());
2551  real_mz = reals.first, real_intens = reals.second;
2552 
2553  if (real_mz <= 0 || real_intens <= 0)
2554  {
2555  typename MSSpectrum<PeakType>::const_iterator h_iter = ref_iter, hc_iter = ref_iter;
2556  while (h_iter != ref.begin())
2557  {
2558  --h_iter;
2559  if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2560  {
2561  hc_iter = h_iter;
2562  }
2563  else
2564  {
2565  break;
2566  }
2567 
2568  if (seed_mz - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.))
2569  {
2570  return false;
2571  }
2572  }
2573  real_mz = h_iter->getMZ(), real_intens = h_iter->getIntensity();
2574  if (real_mz <= 0 || real_intens <= 0)
2575  {
2576  return false;
2577  }
2578  real_mz = h_iter->getMZ();
2579  real_intens = h_iter->getIntensity();
2580  }
2581  }
2582 
2583  DoubleReal c_score = scoreThis_(candidate, peak_cutoff, real_mz, c, 0);
2584 
2585  if (c_score <= 0)
2586  {
2587  return false;
2588  }
2589 
2590  DoubleReal mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(real_mz, c + 1);
2591  typename MSSpectrum<PeakType>::const_iterator real_l_MZ_iter = ref.MZBegin(real_mz - Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.));
2592  typename MSSpectrum<PeakType>::const_iterator real_r_MZ_iter = ref.MZBegin(real_l_MZ_iter, real_mz + mz_cutoff / (c + 1.), ref.end());
2593  if (real_r_MZ_iter == ref.end())
2594  {
2595  --real_r_MZ_iter;
2596  }
2597 
2598 
2599  UInt real_mz_begin = distance(ref.begin(), real_l_MZ_iter);
2600  UInt real_mz_end = distance(ref.begin(), real_r_MZ_iter);
2601 
2602  if (prev_score == -1000)
2603  {
2604  push2Box_(real_mz, scan_index, c, prev_score, transintens, ref.getRT(), real_mz_begin, real_mz_end, real_intens);
2605  }
2606  else
2607  {
2608  push2Box_(real_mz, scan_index, c, c_score, transintens, ref.getRT(), real_mz_begin, real_mz_end, real_intens);
2609  }
2610  return true;
2611  }
2612 
2613  template <typename PeakType>
2615  const MSSpectrum<PeakType>& ref, const DoubleReal seed_mz, const UInt c, const UInt scan_index, const bool check_PPMs, const DoubleReal transintens, const DoubleReal prev_score)
2616  {
2617  typename MSSpectrum<PeakType>::const_iterator iter, ref_iter;
2618  UInt peak_cutoff;
2619  peak_cutoff = IsotopeWavelet::getNumPeakCutOff(seed_mz, c + 1);
2620 
2621  iter = candidate.MZBegin(seed_mz);
2622  //we can ignore those cases
2623  if (iter == candidate.begin() || iter == candidate.end())
2624  {
2625  return false;
2626  }
2627 
2628  std::pair<DoubleReal, DoubleReal> reals;
2629  ref_iter = ref.MZBegin(seed_mz);
2630  //Correct the position
2631  DoubleReal real_mz, real_intens;
2632  if (check_PPMs)
2633  {
2634  reals = checkPPMTheoModel_(ref, iter->getMZ(), c);
2635  real_mz = reals.first, real_intens = reals.second;
2636  //if (real_mz <= 0 || real_intens <= 0)
2637  //{
2638  typename MSSpectrum<PeakType>::const_iterator h_iter = ref_iter, hc_iter = ref_iter;
2639  while (h_iter != ref.begin())
2640  {
2641  --h_iter;
2642  if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2643  {
2644  hc_iter = h_iter;
2645  }
2646  else
2647  {
2648  break;
2649  }
2650 
2651  if (seed_mz - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.))
2652  {
2653  return false;
2654  }
2655  }
2656  ++h_iter;
2657  reals = checkPPMTheoModel_(ref, h_iter->getMZ(), c);
2658  real_mz = reals.first, real_intens = reals.second;
2659 
2660 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2661  std::cout << "Plausibility check old_mz: " << iter->getMZ() << "\t" << real_mz << std::endl;
2662 #endif
2663 
2664  if (real_mz <= 0 || real_intens <= 0)
2665  {
2666  return false;
2667  }
2668  real_mz = h_iter->getMZ();
2669  real_intens = h_iter->getIntensity();
2670  //}
2671  }
2672  else
2673  {
2674  reals = std::pair<DoubleReal, DoubleReal>(seed_mz, ref_iter->getIntensity());
2675  real_mz = reals.first, real_intens = reals.second;
2676 
2677  if (real_mz <= 0 || real_intens <= 0)
2678  {
2679  typename MSSpectrum<PeakType>::const_iterator h_iter = ref_iter, hc_iter = ref_iter;
2680  while (h_iter != ref.begin())
2681  {
2682  --h_iter;
2683  if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2684  {
2685  hc_iter = h_iter;
2686  }
2687  else
2688  {
2689  break;
2690  }
2691 
2692  if (seed_mz - h_iter->getMZ() > Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.))
2693  {
2694  return false;
2695  }
2696  }
2697  real_mz = h_iter->getMZ(), real_intens = h_iter->getIntensity();
2698  if (real_mz <= 0 || real_intens <= 0)
2699  {
2700  return false;
2701  }
2702  real_mz = h_iter->getMZ();
2703  real_intens = h_iter->getIntensity();
2704  }
2705  }
2706 
2707  DoubleReal c_score = scoreThis_(candidate, peak_cutoff, real_mz, c, 0);
2708 
2709  if (c_score <= 0)
2710  {
2711  return false;
2712  }
2713 
2714  DoubleReal mz_cutoff = IsotopeWavelet::getMzPeakCutOffAtMonoPos(real_mz, c + 1);
2715  typename MSSpectrum<PeakType>::const_iterator real_l_MZ_iter = ref.MZBegin(real_mz - Constants::IW_QUARTER_NEUTRON_MASS / (c + 1.));
2716  typename MSSpectrum<PeakType>::const_iterator real_r_MZ_iter = ref.MZBegin(real_l_MZ_iter, real_mz + mz_cutoff / (c + 1.), ref.end());
2717  if (real_r_MZ_iter == ref.end())
2718  {
2719  --real_r_MZ_iter;
2720  }
2721 
2722 
2723  UInt real_mz_begin = distance(ref.begin(), real_l_MZ_iter);
2724  UInt real_mz_end = distance(ref.begin(), real_r_MZ_iter);
2725 
2726  if (prev_score == -1000)
2727  {
2728  push2Box_(real_mz, scan_index, c, prev_score, transintens, ref.getRT(), real_mz_begin, real_mz_end, real_intens);
2729  }
2730  else
2731  {
2732  push2Box_(real_mz, scan_index, c, c_score, transintens, ref.getRT(), real_mz_begin, real_mz_end, real_intens);
2733  }
2734 
2735  return true;
2736  }
2737 
2738  template <typename PeakType>
2739  std::pair<DoubleReal, DoubleReal> IsotopeWaveletTransform<PeakType>::checkPPMTheoModel_(const MSSpectrum<PeakType>& ref, const DoubleReal c_mz, const UInt c)
2740  {
2741  DoubleReal mass = c_mz * (c + 1) - Constants::IW_PROTON_MASS * (c);
2742  DoubleReal ppms = getPPMs_(peptideMassRule_(mass), mass);
2744  {
2745 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2746  std::cout << ::std::setprecision(8) << std::fixed << c_mz << "\t =(" << "ISO_WAVE" << ")> " << "REJECT \t" << ppms << " (rule: " << peptideMassRule_(mass) << " got: " << mass << ")" << std::endl;
2747 #endif
2748  return std::pair<DoubleReal, DoubleReal>(-1, -1);
2749  }
2750 
2751 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2752  std::cout << ::std::setprecision(8) << std::fixed << c_mz << "\t =(" << "ISO_WAVE" << ")> " << "ACCEPT \t" << ppms << " (rule: " << peptideMassRule_(mass) << " got: " << mass << ")" << std::endl;
2753 #endif
2754  return std::pair<DoubleReal, DoubleReal>(c_mz, ref.MZBegin(c_mz)->getIntensity());
2755  }
2756 
2757 } //namespace
2758 
2759 #endif
static DoubleReal getLambdaL(const DoubleReal m)
Returns the mass-parameter lambda (linear fit).
const int CUDA_INIT_SUCCESS
Definition: IsotopeWaveletConstants.h:92
float getMZ()
Definition: IsotopeWaveletTransform.h:72
UInt MZ_end
Definition: IsotopeWaveletTransform.h:103
FeatureMap< Feature > mapSeeds2Features(const MSExperiment< PeakType > &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:2313
virtual void getTransform(MSSpectrum< PeakType > &c_trans, const MSSpectrum< PeakType > &c_ref, const UInt c)
Computes the isotope wavelet transform of charge state c.
Definition: IsotopeWaveletTransform.h:693
const MSSpectrum< PeakType > * getRefSpectrum() const
Definition: IsotopeWaveletTransform.h:188
void mergeFeatures(IsotopeWaveletTransform< PeakType > *later_iwt, const UInt RT_interleave, const UInt RT_votes_cutoff)
DoubleReal getLinearInterpolation(const DoubleReal mz_a, const DoubleReal intens_a, const DoubleReal mz_pos, const DoubleReal mz_b, const DoubleReal intens_b)
Computes a linear (intensity) interpolation.
Definition: IsotopeWaveletTransform.h:353
void setTransIntensity(const UInt i, const DoubleReal intens)
Definition: IsotopeWaveletTransform.h:170
virtual void push2Box_(const DoubleReal mz, const UInt scan, UInt c, const DoubleReal score, const DoubleReal intens, const DoubleReal rt, const UInt MZ_begin, const UInt MZ_end, const DoubleReal ref_intens)
Inserts a potential isotopic pattern into an open box or - if no such box exists - creates a new one...
Definition: IsotopeWaveletTransform.h:1818
UInt max_num_peaks_per_pattern_
Definition: IsotopeWaveletTransform.h:542
virtual void destroy()
Definition: IsotopeWaveletTransform.h:137
A more convenient string class.
Definition: String.h:56
TransSpectrum(const MSSpectrum< PeakType > *reference)
Definition: IsotopeWaveletTransform.h:125
bool intensityPointerComparator(PeakType *a, PeakType *b)
Definition: IsotopeWaveletTransform.h:585
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:55
Size size() const
Definition: MSExperiment.h:117
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:197
DoubleReal getPPMs_(const DoubleReal mass_a, const DoubleReal mass_b) const
Returns the parts-per-million deviation of the masses.
Definition: IsotopeWaveletTransform.h:528
bool hr_data_
Definition: IsotopeWaveletTransform.h:543
std::vector< DoubleReal > interpol_xs_
Definition: IsotopeWaveletTransform.h:539
DoubleReal score
Definition: IsotopeWaveletTransform.h:97
static UInt getMzPeakCutOffAtMonoPos(const DoubleReal mass, const UInt z)
DoubleReal mz
Definition: IsotopeWaveletTransform.h:95
MSSpectrum< PeakType >::const_iterator end() const
Definition: IsotopeWaveletTransform.h:209
DoubleReal av_MZ_spacing_
Definition: IsotopeWaveletTransform.h:537
IntensityType getIntensity() const
Definition: Peak2D.h:161
void clusterSeeds_(const TransSpectrum &candidates, const MSSpectrum< PeakType > &ref, const UInt scan_index, const UInt c, const bool check_PPMs)
Clusters the seeds stored by push2TmpBox_.
Definition: IsotopeWaveletTransform.h:2204
bool myCudaComparator(const cudaHelp &a, const cudaHelp &b)
const double IW_PROTON_MASS
Definition: IsotopeWaveletConstants.h:70
std::multimap< DoubleReal, Box > end_boxes_
Definition: IsotopeWaveletTransform.h:534
DoubleReal getSigma() const
Definition: IsotopeWaveletTransform.h:358
const MSSpectrum< PeakType > * getRefSpectrum()
Definition: IsotopeWaveletTransform.h:182
A container for features.
Definition: FeatureMap.h:111
void sortByIntensity(bool reverse=false)
Sorting.
Definition: ConstRefVector.h:753
unsigned int UInt
Unsigned integer type.
Definition: Types.h:92
#define NULL
Definition: IsotopeWaveletParallelFor.h:41
DoubleReal getSdIntens_(const TransSpectrum &scan, const DoubleReal mean)
Computes the standard deviation (neglecting negative values) of the (transformed) intensities of scan...
Definition: IsotopeWaveletTransform.h:1803
DoubleReal getRT() const
Definition: IsotopeWaveletTransform.h:146
std::vector< DoubleReal > xs_
Definition: IsotopeWaveletTransform.h:538
const int CUDA_MIN_SORT_SIZE
Definition: IsotopeWaveletConstants.h:108
std::vector< DoubleReal > prod_
Definition: IsotopeWaveletTransform.h:538
Iterator begin()
Definition: MSExperiment.h:147
DoubleReal max_mz_cutoff_
Definition: IsotopeWaveletTransform.h:549
Iterator MZEnd(CoordinateType mz)
Binary search for peak range end (returns the past-the-end iterator)
Definition: MSSpectrum.h:530
DoubleReal getAvIntens_(const TransSpectrum &scan)
Computes the average (transformed) intensity (neglecting negative values) of scan.
Definition: IsotopeWaveletTransform.h:1789
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:151
DoubleReal intens
Definition: IsotopeWaveletTransform.h:98
std::vector< std::multimap< DoubleReal, Box > > * tmp_boxes_
Definition: IsotopeWaveletTransform.h:535
static IsotopeWavelet * init(const DoubleReal max_m, const UInt max_charge)
virtual void push2TmpBox_(const DoubleReal mz, const UInt scan, UInt charge, const DoubleReal score, const DoubleReal intens, const DoubleReal rt, const UInt MZ_begin, const UInt MZ_end)
Essentially the same function as.
Definition: IsotopeWaveletTransform.h:1924
const double c
virtual bool checkPositionForPlausibility_(const TransSpectrum &candidate, const MSSpectrum< PeakType > &ref, const DoubleReal seed_mz, const UInt c, const UInt scan_index, const bool check_PPMs, const DoubleReal transintens, const DoubleReal prev_score)
A ugly but necessary function to handle &quot;off-by-1-Dalton predictions&quot; due to idiosyncrasies of the da...
Definition: IsotopeWaveletTransform.h:2614
DoubleReal RT
Definition: IsotopeWaveletTransform.h:100
A 2-dimensional hull representation in [counter]clockwise direction - depending on axis labelling...
Definition: ConvexHull2D.h:75
Internally (only by GPUs) used data structure . It allows efficient data exchange between CPU and GPU...
Definition: IsotopeWaveletTransform.h:112
Iterator MZBegin(CoordinateType mz)
Binary search for peak range begin.
Definition: MSSpectrum.h:506
bool intensityComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:573
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:167
DoubleReal getLinearInterpolation(const typename MSSpectrum< PeakType >::const_iterator &left_iter, const DoubleReal mz_pos, const typename MSSpectrum< PeakType >::const_iterator &right_iter)
Computes a linear (intensity) interpolation.
Definition: IsotopeWaveletTransform.h:342
virtual void getTransformHighRes(MSSpectrum< PeakType > &c_trans, const MSSpectrum< PeakType > &c_ref, const UInt c)
Computes the isotope wavelet transform of charge state c.
Definition: IsotopeWaveletTransform.h:736
void sampleTheCMarrWavelet_(const MSSpectrum< PeakType > &scan, const Int wavelet_length, const Int mz_index, const UInt charge)
DoubleReal min_spacing_
Definition: IsotopeWaveletTransform.h:549
std::vector< DoubleReal > c_mzs_
Definition: IsotopeWaveletTransform.h:538
DoubleReal peptideMassRule_(const DoubleReal c_mass) const
Returns the monoisotopic mass (with corresponding decimal values) we would expect at c_mass...
Definition: IsotopeWaveletTransform.h:505
Int from_max_to_right_
Definition: IsotopeWaveletTransform.h:545
DoubleReal getMZ(const UInt i) const
Definition: IsotopeWaveletTransform.h:152
void setRT(CoordinateType coordinate)
Mutable access to the RT coordinate (index 0)
Definition: Peak2D.h:209
DoubleReal sigma_
Definition: IsotopeWaveletTransform.h:537
UInt MZ_begin
Definition: IsotopeWaveletTransform.h:102
const int CUDA_INIT_FAIL
Definition: IsotopeWaveletConstants.h:93
virtual DoubleReal scoreThis_(const TransSpectrum &candidate, const UInt peak_cutoff, const DoubleReal seed_mz, const UInt c, const DoubleReal ampl_cutoff)
Given a candidate for an isotopic pattern, this function computes the corresponding score...
Definition: IsotopeWaveletTransform.h:1557
virtual void initializeScan(const MSSpectrum< PeakType > &c_ref, const UInt c=0)
Definition: IsotopeWaveletTransform.h:774
A class implementing the isotope wavelet transform. If you just want to find features using the isoto...
Definition: IsotopeWaveletTransform.h:87
static DoubleReal mean(IteratorType begin, IteratorType end)
Calculates the mean of a range of values.
Definition: StatisticFunctions.h:69
DoubleReal getMinSpacing() const
Definition: IsotopeWaveletTransform.h:370
An LC-MS feature.
Definition: Feature.h:66
int sortOnDevice(float *array, int *pos_indices, int numElements, int padding)
Size getMaxScanSize() const
Definition: IsotopeWaveletTransform.h:375
void updateBoxStates(const MSExperiment< PeakType > &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:2029
String intenstype_
Definition: IsotopeWaveletTransform.h:544
virtual void identifyCharge(const MSSpectrum< PeakType > &candidates, const MSSpectrum< PeakType > &ref, const UInt scan_index, const UInt c, const DoubleReal ampl_cutoff, const bool check_PPMs)
Given an isotope wavelet transformed spectrum candidates, this function assigns to every significant ...
Definition: IsotopeWaveletTransform.h:1180
std::vector< DoubleReal > interpol_ys_
Definition: IsotopeWaveletTransform.h:539
IsotopeWaveletTransform()
Default Constructor.
Definition: IsotopeWaveletTransform.h:597
bool positionComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:591
const double IW_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:56
Int from_max_to_left_
Definition: IsotopeWaveletTransform.h:545
CoordinateType getMZ() const
Returns the m/z coordinate (index 1)
Definition: Peak2D.h:191
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
DoubleReal getTransIntensity(const UInt i) const
Definition: IsotopeWaveletTransform.h:164
DoubleReal getAvMZSpacing_(const MSSpectrum< PeakType > &scan)
Computes the average MZ spacing of scan.
Definition: IsotopeWaveletTransform.h:1770
const MSSpectrum< PeakType > * reference_
Definition: IsotopeWaveletTransform.h:223
void setConvexHulls(const std::vector< ConvexHull2D > &hulls)
Set the convex hulls of single mass traces.
MSSpectrum< PeakType >::const_iterator begin() const
Definition: IsotopeWaveletTransform.h:216
const double PEPTIDE_MASS_RULE_BOUND
Definition: IsotopeWaveletConstants.h:52
void setOverallQuality(QualityType q)
Set the overall quality.
virtual void computeMinSpacing(const MSSpectrum< PeakType > &c_ref)
Definition: IsotopeWaveletTransform.h:816
std::vector< float > scores_
Definition: IsotopeWaveletTransform.h:550
MSSpectrum< PeakType >::const_iterator MZEnd(const DoubleReal mz) const
Definition: IsotopeWaveletTransform.h:202
std::vector< DoubleReal > psi_
Definition: IsotopeWaveletTransform.h:538
MSSpectrum< PeakType > c_sorted_candidate_
Definition: IsotopeWaveletTransform.h:548
UInt data_length_
Definition: IsotopeWaveletTransform.h:542
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
bool intensityAscendingComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:579
const int CUDA_EXTENDED_BLOCK_SIZE_MAX
Definition: IsotopeWaveletConstants.h:97
TransSpectrum()
Definition: IsotopeWaveletTransform.h:119
virtual ~IsotopeWaveletTransform()
Destructor.
Definition: IsotopeWaveletTransform.h:678
void clear(bool clear_meta_data)
Clears all data and meta data.
Definition: MSExperiment.h:809
float mz
Definition: IsotopeWaveletTransform.h:77
UInt c
Definition: IsotopeWaveletTransform.h:96
float score
Definition: IsotopeWaveletTransform.h:79
std::vector< int > indices_
Definition: IsotopeWaveletTransform.h:546
MSSpectrum< PeakType >::const_iterator MZBegin(const DoubleReal mz) const
Definition: IsotopeWaveletTransform.h:195
void setCharge(const ChargeType &ch)
Set charge state.
void addPoints(const PointArrayType &points)
std::multimap< DoubleReal, Box > closed_boxes_
Definition: IsotopeWaveletTransform.h:534
const int CUDA_BLOCK_SIZE_MAX
Definition: IsotopeWaveletConstants.h:96
std::vector< float > * trans_intens_
Definition: IsotopeWaveletTransform.h:224
Size size() const
Definition: IsotopeWaveletTransform.h:176
virtual ~TransSpectrum()
Definition: IsotopeWaveletTransform.h:132
UInt max_charge_
Definition: IsotopeWaveletTransform.h:542
std::multimap< UInt, BoxElement > Box
Key: RT index, value: BoxElement.
Definition: IsotopeWaveletTransform.h:106
void extendBox_(const MSExperiment< PeakType > &map, const Box box)
A currently still necessary function that extends the box box in order to capture also signals whose ...
Definition: IsotopeWaveletTransform.h:2098
DoubleReal getRT() const
Definition: MSSpectrum.h:215
This vector holds pointer to the elements of another container.
Definition: ConstRefVector.h:72
const double IW_HALF_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:57
float intens
Definition: IsotopeWaveletTransform.h:78
UInt RT_index
Definition: IsotopeWaveletTransform.h:101
const double PEPTIDE_MASS_RULE_FACTOR
Definition: IsotopeWaveletConstants.h:51
float getIntensity()
Definition: IsotopeWaveletTransform.h:74
DoubleReal getRefIntensity(const UInt i) const
Definition: IsotopeWaveletTransform.h:158
Size max_scan_size_
Definition: IsotopeWaveletTransform.h:541
virtual std::multimap< DoubleReal, Box > getClosedBoxes()
Returns the closed boxes.
Definition: IsotopeWaveletTransform.h:334
virtual std::pair< DoubleReal, DoubleReal > checkPPMTheoModel_(const MSSpectrum< PeakType > &ref, const DoubleReal c_mz, const UInt c)
Definition: IsotopeWaveletTransform.h:2739
void scoreOnDevice(int *sorted_positions_indices, float *trans_intensities, float *pos, float *scores, const int c, const int num_of_scores, const int overall_size, const unsigned int max_peak_cutoff, const float ampl_cutoff)
int Int
Signed integer type.
Definition: Types.h:100
const double IW_QUARTER_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:58
void setSigma(const DoubleReal sigma)
Definition: IsotopeWaveletTransform.h:363
const unsigned int DEFAULT_NUM_OF_INTERPOLATION_POINTS
Definition: IsotopeWaveletConstants.h:45
void getExternalCudaTransforms(dim3 dimGrid, dim3 dimBlock, float *positions_dev, float *intensities_dev, int from_max_to_left, int from_max_to_right, float *result_dev, const int charge, const int to_load, const int to_compute, const int size, float *fwd2, bool highres=false)
static UInt getNumPeakCutOff(const DoubleReal mass, const UInt z)
std::vector< float > zeros_
Definition: IsotopeWaveletTransform.h:550
std::multimap< DoubleReal, Box > front_boxes_
Definition: IsotopeWaveletTransform.h:534
DoubleReal ref_intens
Definition: IsotopeWaveletTransform.h:99
static DoubleReal getValueByLambda(const DoubleReal lambda, const DoubleReal 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.
std::multimap< DoubleReal, Box > open_boxes_
Definition: IsotopeWaveletTransform.h:534
Internally used data structure.
Definition: IsotopeWaveletTransform.h:93
double DoubleReal
Double-precision real type.
Definition: Types.h:118
std::vector< DoubleReal > c_spacings_
Definition: IsotopeWaveletTransform.h:538
An internally used class, subsuming several variables.
Definition: IsotopeWaveletTransform.h:69
const double PEPTIDE_MASS_RULE_THEO_PPM_BOUND
Definition: IsotopeWaveletConstants.h:53

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