Parallelization, fft convolve and binning for FFMI

src/openms/include/OpenMS/ANALYSIS/OPENSWATH/MRMTransitionGroupPicker.h

@@ -46,6 +46,7 @@
 
 #include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
 #include <OpenMS/CONCEPT/LogStream.h>
+#include <omp.h>
 
 // Cross-correlation
 #include <OpenMS/OPENSWATHALGO/ALGO/Scoring.h>
@@ -115,45 +116,57 @@ namespace OpenMS
       OPENMS_PRECONDITION(transition_group.isInternallyConsistent(), "Consistent state required")
       OPENMS_PRECONDITION(transition_group.chromatogramIdsMatch(), "Chromatogram native IDs need to match keys in transition group")
 
-      std::vector<MSChromatogram > picked_chroms;
-      std::vector<MSChromatogram > smoothed_chroms;
+      std::vector<MSChromatogram> picked_chroms;
+      std::vector<MSChromatogram> smoothed_chroms;
 
-      // Pick fragment ion chromatograms
-      for (Size k = 0; k < transition_group.getChromatograms().size(); k++)
+      picked_chroms.resize(transition_group.getChromatograms().size());
+      smoothed_chroms.resize(transition_group.getChromatograms().size());
+
+      // Pick precursor chromatograms
+      if (use_precursors_)
       {
-        MSChromatogram& chromatogram = transition_group.getChromatograms()[k];
-        String native_id = chromatogram.getNativeID();
-
-        // only pick detecting transitions (skip all others)
-        if (transition_group.getTransitions().size() > 0 && 
-            transition_group.hasTransition(native_id)  && 
-            !transition_group.getTransition(native_id).isDetectingTransition() )
+        #pragma omp parallel num_threads(transition_group.getPrecursorChromatograms().size())
         {
-          continue;
-        }
-
-        MSChromatogram picked_chrom, smoothed_chrom;
-        smoothed_chrom.setNativeID(native_id);
-        picker_.pickChromatogram(chromatogram, picked_chrom, smoothed_chrom);
-        picked_chrom.sortByIntensity();
-        picked_chroms.push_back(std::move(picked_chrom));
-        smoothed_chroms.push_back(std::move(smoothed_chrom));
-      }
+          PeakPickerMRM picker_temp(picker_);
+          #pragma omp for
+          for (Size k = 0; k < transition_group.getPrecursorChromatograms().size(); k++)
+          {
+            SpectrumT& picked_chrom = picked_chroms[k];
+            SpectrumT& smoothed_chrom = smoothed_chroms[k];
+            SpectrumT& chromatogram = transition_group.getPrecursorChromatograms()[k];
 
-      // Pick precursor chromatograms
-      if (use_precursors_)
+            picker_temp.pickChromatogram(chromatogram, picked_chrom, smoothed_chrom);
+            picked_chrom.sortByIntensity();
+          }
+        }
+      } 
+      else 
       {
-        for (Size k = 0; k < transition_group.getPrecursorChromatograms().size(); k++)
+        #pragma omp parallel num_threads(transition_group.getChromatograms().size())
         {
-          SpectrumT picked_chrom, smoothed_chrom;
-          SpectrumT& chromatogram = transition_group.getPrecursorChromatograms()[k];
+          PeakPickerMRM picker_temp(picker_);
+          #pragma omp for
+          for (Size k = 0; k < transition_group.getChromatograms().size(); k++) // Pick fragment ion chromatograms
+          {
+            MSChromatogram& chromatogram = transition_group.getChromatograms()[k];
+            String native_id = chromatogram.getNativeID();
 
-          picker_.pickChromatogram(chromatogram, picked_chrom, smoothed_chrom);
-          picked_chrom.sortByIntensity();
-          picked_chroms.push_back(picked_chrom);
-          smoothed_chroms.push_back(smoothed_chrom);
+            // only pick detecting transitions (skip all others)
+            if (transition_group.getTransitions().size() > 0 && 
+              transition_group.hasTransition(native_id)  && 
+              !transition_group.getTransition(native_id).isDetectingTransition() )
+            {
+            continue;
+            }
+
+            MSChromatogram& picked_chrom = picked_chroms[k];
+            MSChromatogram& smoothed_chrom = smoothed_chroms[k];
+            smoothed_chrom.setNativeID(native_id);
+            picker_temp.pickChromatogram(chromatogram, picked_chrom, smoothed_chrom);
+            picked_chrom.sortByIntensity();
+          }
         }
-      }
+      } 
 
       // Find features (peak groups) in this group of transitions.
       // While there are still peaks left, one will be picked and used to create

src/openms/include/OpenMS/FILTERING/SMOOTHING/GaussFilterAlgorithm.h

@@ -37,8 +37,11 @@
 #include <OpenMS/CONCEPT/Types.h>
 #include <OpenMS/CONCEPT/Constants.h>
 #include <OpenMS/INTERFACES/DataStructures.h>
+#include <Evergreen/evergreen.hpp>
 #include <cmath>
 #include <vector>
+#include <numeric>      // Reihenfolge TODO
+#include <stdio.h>
 
 namespace OpenMS
 {
@@ -126,7 +129,7 @@ namespace OpenMS
     /**
       @brief Smoothes an two data arrays containing data.
 
-      Convolutes the filter and the profile data and writes the results into the output iterators mz_out and int_out. 
+      Convolutes the filter and the profile data and writes the results into the output iterators mz_out and int_out.
     */
     template <typename ConstIterT, typename IterT>
     bool filter(
@@ -137,56 +140,150 @@ namespace OpenMS
         IterT int_out)
     {
       bool found_signal = false;
-
-      ConstIterT mz_it = mz_in_start;
+
       ConstIterT int_it = int_in_start;
-      for (; mz_it != mz_in_end; mz_it++, int_it++)
-      {
-        // if ppm tolerance is used, calculate a reasonable width value for this m/z
+      ConstIterT mz_it = mz_in_start;
+      // if ppm tolerance is used, calculate a reasonable width value for this m/z
         if (use_ppm_tolerance_)
-        {
-          initialize((*mz_it) * ppm_tolerance_ * 10e-6, spacing_, ppm_tolerance_, use_ppm_tolerance_ );
+        { 
+          for (; mz_it != mz_in_end; mz_it++, int_it++)
+          {
+            initialize((*mz_it) * ppm_tolerance_ * 10e-6, spacing_, ppm_tolerance_, use_ppm_tolerance_ ); 
+
+            double new_int = integrate_(mz_it, int_it, mz_in_start, mz_in_end);
+
+            // store new intensity and m/z into output iterator
+            *mz_out = *mz_it;
+            *int_out = new_int;
+            ++mz_out;
+            ++int_out;
+
+            if (new_int != 0) found_signal = true;
+          }
         }
+        else 
+        {
+          evergreen::Tensor<double> gaussian_data({coeffs_.size()});
 
-        double new_int = integrate_(mz_it, int_it, mz_in_start, mz_in_end);
-
-        // store new intensity and m/z into output iterator
-        *mz_out = *mz_it;
-        *int_out = new_int;
-        ++mz_out;
-        ++int_out;
+          // normalize the gaussian coefficients to scale the intensities correctly
+          double sum = std::accumulate(coeffs_.begin(), coeffs_.end(), 0.0);
+          if (sum != 1.0)
+          {
+            for (long unsigned i = 0; i < coeffs_.size(); i++)
+            {
+              gaussian_data[i] = coeffs_[i] / sum;
+            }
+          } 
+          else
+          {
+            for (long unsigned i = 0; i < coeffs_.size(); i++)
+            {
+              gaussian_data[i] = coeffs_[i];
+            }
+          }
+
+          // saving the spacing_ nearest mz values in the same bin (binning)
+          std::vector<int> bin_numbers = {};
+          evergreen::Tensor<double> bins({(double(*(mz_in_end - 1) - *mz_in_start) / spacing_) + 1.0});
+          binning(bins, bin_numbers, mz_in_start, mz_in_end, int_it);
+
+          //Fast Fourier Transformation from evergreen
+          evergreen::Tensor<double> out_fft = evergreen::fft_convolve(bins, gaussian_data);
+
+          mz_it = mz_in_start;
+          int_it = int_in_start;
 
-        if (fabs(new_int) > 0) found_signal = true;
+          out_fft = debinning(out_fft, bin_numbers, mz_in_start, mz_in_end, int_it);
+
+          mz_it = mz_in_start;
+          int_it = int_in_start;
+
+          //cutting off the access data that was calculated by convolve
+          for (int i = ((out_fft.flat_size() - bins.flat_size())/ 2); mz_it != mz_in_end; mz_it++, int_it++, i++)
+          {
+            *mz_out = *mz_it;
+            *int_out = out_fft[i];
+            ++mz_out;
+            ++int_out;
+
+            if (out_fft[i] != 0) found_signal = true;
+          }
       }
       return found_signal;
     }
 
+
     void initialize(double gaussian_width, double spacing, double ppm_tolerance, bool use_ppm_tolerance);
 
-protected:
+    template <typename ConstIterT, typename DataPackage, typename BinPosPackage>
+    void binning(DataPackage& bins, BinPosPackage& bin_numbers, ConstIterT mz_it, ConstIterT mz_end, ConstIterT int_it) const
+    {
+      double start = *mz_it;
+      unsigned j = 0;
+      for (unsigned i = 0; mz_it != mz_end; j++)
+      {
+        if (*mz_it < (start + (i + 1) * spacing_))
+        {
+          bins[i] += *int_it;
+          int_it++;
+          mz_it++;
+          bin_numbers.push_back(1);
 
-    ///Coefficients
-    std::vector<double> coeffs_;
-    /// The standard derivation  \f$ \sigma \f$.
-    double sigma_;
-    /// The spacing of the pre-tabulated kernel coefficients
-    double spacing_;
+          continue;
+        }
 
-    // tolerance in ppm
-    bool use_ppm_tolerance_;
-    double ppm_tolerance_;
+        i++;
+        bin_numbers.push_back(0);
+      }
+    }
+
+
+    template <typename ConstIterT, typename DataPackage, typename BinPosPackage>
+    DataPackage debinning(DataPackage& out_fft, BinPosPackage& bin_numbers, ConstIterT mz_it, ConstIterT mz_end, ConstIterT int_it) const
+    {
+      unsigned start_bin = 0;
+      unsigned end_bin = 0;
+      for (unsigned i = 0; i < out_fft.flat_size() - 1; i++)
+      {
+        // search for bin with original data
+        if (bin_numbers[i] != 0)
+        {
+          end_bin = i;
+        }
+        else 
+        {
+          end_bin++;
+          continue;
+        }
+
+        // adding the values in the [start_bin, end_bin] interval (while considering the distance to the start_bin and end_bin peaks)
+        unsigned j = start_bin + 1;
+        while (j != end_bin)
+        {
+          out_fft[end_bin] += out_fft[j] * ((1.0 - double(end_bin - j)) / double(end_bin - start_bin));
+          out_fft[start_bin] += out_fft[j] * ((1.0 - double(j - start_bin)) / double(end_bin - start_bin));
+
+          out_fft[j] = 0.0;
+          j++;
+        }
+        start_bin = end_bin;
+      }
 
+      return out_fft;
+    }
+
+protected:
     /// Computes the convolution of the raw data at position x and the gaussian kernel
     template <typename InputPeakIterator>
-    double integrate_(InputPeakIterator x /* mz */, InputPeakIterator y /* int */, InputPeakIterator first, InputPeakIterator last)
+    double integrate_(InputPeakIterator x /* mz */, InputPeakIterator y /* int */, InputPeakIterator first, InputPeakIterator last) const
     {
       double v = 0.;
       // norm the gaussian kernel area to one
       double norm = 0.;
       Size middle = coeffs_.size();
 
-      double start_pos = (( (*x) - (middle * spacing_)) > (*first)) ? ((*x) - (middle * spacing_)) : (*first);
-      double end_pos = (( (*x) + (middle * spacing_)) < (*(last - 1))) ? ((*x) + (middle * spacing_)) : (*(last - 1));
+      double start_pos = (((*x) - (middle * spacing_)) > (*first)) ? ((*x) - (middle * spacing_)) : (*first);
+      double end_pos = (((*x) + (middle * spacing_)) < (*(last - 1))) ? ((*x) + (middle * spacing_)) : (*(last - 1));
 
       InputPeakIterator help_x = x;
       InputPeakIterator help_y = y;
@@ -203,7 +300,7 @@ namespace OpenMS
         Size left_position = (Size)floor(distance_in_gaussian / spacing_);
 
         // search for the true left adjacent data point (because of rounding errors)
-        for (int j = 0; ((j < 3) &&  (distance(first, help_x - j) >= 0)); ++j)
+        for (UInt j = 0; ((j < 3) && (distance(first, help_x - j) >= 0)); ++j)
         {
           if (((left_position - j) * spacing_ <= distance_in_gaussian) && ((left_position - j + 1) * spacing_ >= distance_in_gaussian))
           {
@@ -234,7 +331,6 @@ namespace OpenMS
         std::cout << "interpolated value left " << coeffs_right << std::endl;
 #endif
 
-
         // search for the corresponding datapoint for (help-1) in the gaussian (take the left most adjacent point)
         distance_in_gaussian = fabs((*x) - (*(help_x - 1)));
         left_position = (Size)floor(distance_in_gaussian / spacing_);
@@ -273,7 +369,6 @@ namespace OpenMS
                   << std::endl;
 #endif
 
-
         norm += fabs((*(help_x - 1)) - (*help_x)) / 2. * (coeffs_left + coeffs_right);
 
         v += fabs((*(help_x - 1)) - (*help_x)) / 2. * (*(help_y - 1) * coeffs_left + (*help_y) * coeffs_right);
@@ -364,6 +459,7 @@ namespace OpenMS
                   << "* " << coeffs_right
                   << std::endl;
 #endif
+
         norm += fabs((*help_x) - (*(help_x + 1)) ) / 2. * (coeffs_left + coeffs_right);
 
         v += fabs((*help_x) - (*(help_x + 1)) ) / 2. * ((*help_y) * coeffs_left + (*(help_y + 1)) * coeffs_right);
@@ -381,6 +477,17 @@ namespace OpenMS
       }
     }
 
+    /// Coefficients
+    std::vector<double> coeffs_;          // Gaussian numbers for FFT
+    /// The standard derivation \f$ \sigma \f$.
+    double sigma_;
+    /// The spacing of the pre-tabulated kernel coefficients
+    double spacing_;
+
+    // tolerance in ppm
+    bool use_ppm_tolerance_;
+    double ppm_tolerance_;
+
   };
 
 } // namespace OpenMS

src/openms/include/OpenMS/FILTERING/SMOOTHING/SavitzkyGolayFilter.h

@@ -112,7 +112,7 @@ namespace OpenMS
     // low level template to filters spectra and chromatograms
     // raw data and meta data needs to be copied to the output container before calling this function
     template<class InputIt, class OutputIt>
-    void filter(InputIt first, InputIt last, OutputIt d_first)
+    void filter(InputIt first, InputIt last, OutputIt d_first) const
     {
       size_t n = std::distance(first, last);
 
@@ -198,7 +198,7 @@ namespace OpenMS
     /**
       @brief Removed the noise from an MSChromatogram
     */
-    void filter(MSChromatogram & chromatogram)
+    void filter(MSChromatogram & chromatogram) const
     {
       // copy the data AND META DATA to the output container
       MSChromatogram output = chromatogram;

src/openms/source/ANALYSIS/OPENSWATH/PeakPickerMRM.cpp

@@ -125,6 +125,8 @@ namespace OpenMS
       gauss_.filter(smoothed_chrom);
     }
 
+    // läuft top
+
     // Find initial seeds (peak picking)
     pp_.pick(smoothed_chrom, picked_chrom);
     OPENMS_LOG_DEBUG << "Picked " << picked_chrom.size() << " chromatographic peaks." << std::endl;
@@ -165,9 +167,9 @@ namespace OpenMS
     picked_chrom.getFloatDataArrays()[IDX_RIGHTBORDER].reserve(picked_chrom.size());
     for (Size i = 0; i < picked_chrom.size(); i++)
     {
-      picked_chrom.getFloatDataArrays()[IDX_ABUNDANCE].push_back(integrated_intensities_[i]);
-      picked_chrom.getFloatDataArrays()[IDX_LEFTBORDER].push_back((float)chromatogram[left_width_[i]].getRT());
-      picked_chrom.getFloatDataArrays()[IDX_RIGHTBORDER].push_back((float)chromatogram[right_width_[i]].getRT());
+      picked_chrom.getFloatDataArrays()[IDX_ABUNDANCE][i] = integrated_intensities_[i];
+      picked_chrom.getFloatDataArrays()[IDX_LEFTBORDER][i] = (float)chromatogram[left_width_[i]].getRT();
+      picked_chrom.getFloatDataArrays()[IDX_RIGHTBORDER][i] = (float)chromatogram[right_width_[i]].getRT();
     }
   }