TraceFitter.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-2012.
//
// 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: Stephan Aiche$
// $Authors: Stephan Aiche, Marc Sturm $
// --------------------------------------------------------------------------

#ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_TRACEFITTER_H
#define OPENMS_TRANSFORMATIONS_FEATUREFINDER_TRACEFITTER_H

#include <OpenMS/CONCEPT/LogStream.h>

#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/FeatureFinderAlgorithmPickedHelperStructs.h>

#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>

namespace OpenMS
{

  template <class PeakType>
  class TraceFitter :
    public DefaultParamHandler
  {

public:
    TraceFitter() :
      DefaultParamHandler("TraceFitter")
    {
      this->defaults_.setValue("max_iteration", 500, "Maximum number of iterations using by Levenberg-Marquardt algorithm.", StringList::create("advanced"));
      this->defaults_.setValue("epsilon_abs", 0.0001, "Absolute error used by the Levenberg-Marquardt algorithm.", StringList::create("advanced"));
      this->defaults_.setValue("epsilon_rel", 0.0001, "Relative error used by the Levenberg-Marquardt algorithm.", StringList::create("advanced"));
    }

    TraceFitter(const TraceFitter & source) :
      DefaultParamHandler(source),
      epsilon_abs_(source.epsilon_abs_),
      epsilon_rel_(source.epsilon_rel_),
      max_iterations_(source.max_iterations_)
    {
    }

    virtual TraceFitter & operator=(const TraceFitter & source)
    {
      DefaultParamHandler::operator=(source);
      max_iterations_ = source.max_iterations_;
      epsilon_abs_  = source.epsilon_abs_;
      epsilon_rel_ = source.epsilon_rel_;
      updateMembers_();

      return *this;
    }

    virtual ~TraceFitter()
    {
    }

    virtual void fit(FeatureFinderAlgorithmPickedHelperStructs::MassTraces<PeakType> & traces) = 0;

    virtual DoubleReal getLowerRTBound() const = 0;

    virtual DoubleReal getUpperRTBound() const = 0;

    virtual DoubleReal getHeight() const = 0;

    virtual DoubleReal getCenter() const = 0;

    virtual DoubleReal getFWHM() const = 0;

    virtual DoubleReal computeTheoretical(const FeatureFinderAlgorithmPickedHelperStructs::MassTrace<PeakType> & trace, Size k) = 0;

    virtual bool checkMinimalRTSpan(const std::pair<DoubleReal, DoubleReal> & rt_bounds, const DoubleReal min_rt_span) = 0;

    virtual bool checkMaximalRTSpan(const DoubleReal max_rt_span) = 0;

    virtual DoubleReal getFeatureIntensityContribution() = 0;

    virtual String getGnuplotFormula(FeatureFinderAlgorithmPickedHelperStructs::MassTrace<PeakType> const & trace, const char function_name, const DoubleReal baseline, const DoubleReal rt_shift) = 0;

protected:

    virtual void printState_(SignedSize iter, gsl_multifit_fdfsolver * s) = 0;

    virtual void updateMembers_()
    {
      max_iterations_ = this->param_.getValue("max_iteration");
      epsilon_abs_ = this->param_.getValue("epsilon_abs");
      epsilon_rel_ = this->param_.getValue("epsilon_rel");
    }

    virtual void getOptimizedParameters_(gsl_multifit_fdfsolver * s) = 0;

    void optimize_(FeatureFinderAlgorithmPickedHelperStructs::MassTraces<PeakType> & traces, const Size num_params, double x_init[],
                   Int (* residual)(const gsl_vector * x, void * params, gsl_vector * f),
                   Int (* jacobian)(const gsl_vector * x, void * params, gsl_matrix * J),
                   Int (* evaluate)(const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix * J))
    {
      const gsl_multifit_fdfsolver_type * T;
      gsl_multifit_fdfsolver * s;

      const size_t data_count = traces.getPeakCount();

      // gsl always expects N>=p or default gsl error handler invoked,
      // cause Jacobian be rectangular M x N with M>=N
      if (data_count < num_params) throw Exception::UnableToFit(__FILE__, __LINE__, __PRETTY_FUNCTION__, "UnableToFit-FinalSet", "Skipping feature, gsl always expects N>=p");

      gsl_multifit_function_fdf func;
      gsl_vector_view x = gsl_vector_view_array(x_init, num_params);
      gsl_rng_env_setup();
      func.f = (residual);
      func.df = (jacobian);
      func.fdf = (evaluate);
      func.n = data_count;
      func.p = num_params;
      func.params = &traces;
      T = gsl_multifit_fdfsolver_lmsder;
      s = gsl_multifit_fdfsolver_alloc(T, data_count, num_params);
      gsl_multifit_fdfsolver_set(s, &func, &x.vector);
      SignedSize iter = 0;
      Int gsl_status_;
      do
      {
        iter++;
        gsl_status_ = gsl_multifit_fdfsolver_iterate(s);
        printState_(iter, s);
        if (gsl_status_) break;
        gsl_status_ = gsl_multifit_test_delta(s->dx, s->x, epsilon_abs_, epsilon_rel_);
      }
      while (gsl_status_ == GSL_CONTINUE && iter < max_iterations_);

      // get the parameters out of the fdfsolver
      getOptimizedParameters_(s);

      gsl_multifit_fdfsolver_free(s);
    }


    DoubleReal epsilon_abs_;
    DoubleReal epsilon_rel_;
    SignedSize max_iterations_;

  };

}

#endif // #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKED_RTFITTING_H