TwoDOptimization.h

Go to the documentation of this file.

// Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: $
// --------------------------------------------------------------------------
 
#pragma once
 
//#define DEBUG_2D
#undef DEBUG_2D
 
#ifdef DEBUG_2D
#include <iostream>
#include <fstream>
#endif
 
#include <vector>
#include <utility>
#include <cmath>
#include <set>
 
#include <OpenMS/KERNEL/StandardTypes.h>
#include <OpenMS/TRANSFORMATIONS/RAW2PEAK/PeakShape.h>
#include <OpenMS/TRANSFORMATIONS/RAW2PEAK/OptimizePeakDeconvolution.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/KERNEL/MSSpectrum.h>
#include <OpenMS/KERNEL/PeakIndex.h>
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/DATASTRUCTURES/IsotopeCluster.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/MATH/MISC/MathFunctions.h>
 
 
#ifndef OPENMS_SYSTEM_STOPWATCH_H
#endif
 
#include <OpenMS/TRANSFORMATIONS/RAW2PEAK/OptimizePick.h>
#include <OpenMS/TRANSFORMATIONS/RAW2PEAK/PeakShape.h>
 
namespace OpenMS
{
 
  class OPENMS_DLLAPI TwoDOptimization :
    public DefaultParamHandler
  {
public:
 
    typedef MSExperiment::const_iterator InputSpectrumIterator;
 
    TwoDOptimization();
 
    TwoDOptimization(const TwoDOptimization& opt);
 
    ~TwoDOptimization() override{}
 
    TwoDOptimization& operator=(const TwoDOptimization& opt);
 
 
    inline double getMZTolerance() const {return tolerance_mz_; }
    inline void setMZTolerance(double tolerance_mz)
    {
      tolerance_mz_ = tolerance_mz;
      param_.setValue("2d:tolerance_mz", tolerance_mz);
    }
 
    inline double getMaxPeakDistance() const {return max_peak_distance_; }
    inline void setMaxPeakDistance(double max_peak_distance)
    {
      max_peak_distance_ = max_peak_distance;
      param_.setValue("2d:max_peak_distance", max_peak_distance);
    }
 
    inline UInt getMaxIterations() const {return max_iteration_; }
    inline void setMaxIterations(UInt max_iteration)
    {
      max_iteration_ = max_iteration;
      param_.setValue("iterations", max_iteration);
    }
 
    inline const OptimizationFunctions::PenaltyFactorsIntensity& getPenalties() const {return penalties_; }
    inline void setPenalties(OptimizationFunctions::PenaltyFactorsIntensity& penalties)
    {
      penalties_ = penalties;
      param_.setValue("penalties:position", penalties.pos);
      param_.setValue("penalties:height", penalties.height);
      param_.setValue("penalties:left_width", penalties.lWidth);
      param_.setValue("penalties:right_width", penalties.rWidth);
    }
 
    void optimize(InputSpectrumIterator first,
                  InputSpectrumIterator last,
                  PeakMap& ms_exp, bool real2D = true);
 
 
protected:
    struct Data
    {
      std::vector<std::pair<SignedSize, SignedSize> > signal2D;
      std::multimap<double, IsotopeCluster>::iterator iso_map_iter;
      Size total_nr_peaks;
      std::map<Int, std::vector<PeakIndex> > matching_peaks;
      PeakMap picked_peaks;
      PeakMap::ConstIterator raw_data_first;
      OptimizationFunctions::PenaltyFactorsIntensity penalties;
      std::vector<double> positions;
      std::vector<double> signal;
    };
 
    class OPENMS_DLLAPI TwoDOptFunctor
    {
    public:
      int inputs() const { return m_inputs; }
      int values() const { return m_values; }
 
      TwoDOptFunctor(unsigned dimensions, unsigned num_data_points, const TwoDOptimization::Data* data)
      : m_inputs(dimensions), m_values(num_data_points), m_data(data) {}
 
      int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec);
      // compute Jacobian matrix for the different parameters
      int df(const Eigen::VectorXd &x, Eigen::MatrixXd &J);
 
    private:
      const int m_inputs, m_values;
      const TwoDOptimization::Data* m_data;
    };
 
 
    std::multimap<double, IsotopeCluster> iso_map_;
 
    std::multimap<double, IsotopeCluster>::const_iterator curr_region_;
 
    double max_peak_distance_;
 
    double tolerance_mz_;
 
    //  std::map<Int, std::vector<PeakMap::SpectrumType::Iterator > > matching_peaks_;
    std::map<Int, std::vector<PeakIndex> > matching_peaks_;
 
 
    UInt max_iteration_;
 
    bool real_2D_;
 
 
    OptimizationFunctions::PenaltyFactorsIntensity penalties_;
 
 
    std::vector<double>::iterator searchInScan_(std::vector<double>::iterator scan_begin,
                                                    std::vector<double>::iterator scan_end,
                                                    double current_mz);
 
    void optimizeRegions_(InputSpectrumIterator& first,
                          InputSpectrumIterator& last,
                          PeakMap& ms_exp);
 
    void optimizeRegionsScanwise_(InputSpectrumIterator& first,
                                  InputSpectrumIterator& last,
                                  PeakMap& ms_exp);
 
 
    void getRegionEndpoints_(PeakMap& exp,
                             InputSpectrumIterator& first,
                             InputSpectrumIterator& last,
                             Size iso_map_idx,
                             double noise_level,
                             TwoDOptimization::Data& d);
 
    void findMatchingPeaks_(std::multimap<double, IsotopeCluster>::iterator& it,
                            PeakMap& ms_exp);
 
 
    void updateMembers_() override;
  };
 
}