TwoDOptimization.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2022.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $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;
  };

}