JGandalf.hh

Go to the documentation of this file.

#ifndef __JFIT__JGANDALF__
#define __JFIT__JGANDALF__

#include <limits>
#include <vector>
#include <cmath>
#include <ostream>
#include <iomanip>

#include "Jeep/JPrint.hh"
#include "Jeep/JMessage.hh"
#include "JMath/JMatrixNS.hh"
#include "JLang/JException.hh"


/**
 * \author mdejong
 */

namespace JFIT {}
namespace JPP { using namespace JFIT; }

namespace JFIT {

  using JEEP::JMessage;
  using JLANG::JException;

  
  /**
   * Fit method based on the Levenberg-Marquardt method.
   *
   * The template argument refers to the model that should be fitted to the data.
   * This data structure should have arithmetic capabalities.
   *
   * The data member JGandalf::value corresponds to the start or final value of
   * the model of the fit procedure and JGandalf::error to the uncertainties.
   * The co-variance matrix is stored in data member JGandalf::H.
   * The data member JGandalf::parameters is a list of pointers to those 
   * data members of the model that should be fitted.
   * The template fit function should return the data type JGandalf::result_type
   * which is composed of the values of the chi2 and gradient of a data point, respectively.
   * The function operator returns the chi2 of the fit.
   */
  template<class JModel_t>
  class JGandalf :
    public JMessage< JGandalf<JModel_t> >
  {
  public:

    using JMessage< JGandalf<JModel_t> >::debug;


    /**
     * Data type of fit parameter.
     */
    typedef typename JModel_t::parameter_type  parameter_type;


    /**
     * Data structure for return value of fit function.
     */
    struct result_type {
      /**
       * Default constructor.
       */
      result_type() :
        chi2(0.0),
        gradient()
      {}
      

      /**
       * Constructor.
       *
       * \param  __chi2       chi2
       * \param  __model      gradient
       */
      result_type(const double   __chi2,
                  const JModel_t __model) :
        chi2(__chi2),
        gradient(__model)
      {}
      

      /**
       * Type conversion operator.
       *
       * \return              chi2
       */
      operator double() const 
      {
        return chi2;
      }


      double   chi2;                  //!<   chi2 of data point
      JModel_t gradient;              //!< d(chi2)/d(...)
    };


    /**
     * Default constructor.
     */
    JGandalf()
    {}


    /**
     * Multi-dimensional fit of two data sets.
     *
     * The fit function should return the equivalent of chi2 for the current value
     * of the model and the given data point as well as the partial derivatives.
     *
     * \param  fit            fit function
     * \param  __begin1       begin of first  data set
     * \param  __end1         end   of first  data set
     * \param  __begin2       begin of second data set
     * \param  __end2         end   of second data set
     * \return                chi2
     */
    template<class JFunction_t, class T1, class T2>
    double operator()(const JFunction_t& fit,
                      T1 __begin1, T1 __end1,
                      T2 __begin2, T2 __end2)
    {
      using namespace std;

      const int N = parameters.size();

      H.resize(N);
      h.resize(N);

      previous  = value;
      error     = JModel_t();
      lambda    = LAMBDA_MIN;

      double chi2_old = numeric_limits<double>::max();

      for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {

        DEBUG("step:     " << numberOfIterations << endl);

        reset();

        evaluate(fit, __begin1, __end1);
        evaluate(fit, __begin2, __end2);

        DEBUG("lambda:   " << SCIENTIFIC(12,5) << lambda << endl);
        DEBUG("chi2:     " << FIXED     (12,5) << chi2   << endl);

        if (chi2 < chi2_old) {

          if (numberOfIterations != 0) {

            if (fabs(chi2_old - chi2) < EPSILON*fabs(chi2_old)) {
              return chi2;
            }

            lambda /= LAMBDA_DOWN;
          }

          chi2_old = chi2;
          previous = value;

        } else {

          value   = previous;
          lambda *= LAMBDA_UP;
            
          if (lambda > LAMBDA_MAX) {
            return chi2_old;                                   // no improvement found
          }

          reset();
          
          evaluate(fit, __begin1, __end1);
          evaluate(fit, __begin2, __end2);
        }

        DEBUG("Hesse matrix:" << endl);
        DEBUG(H << endl);


        // force definite positiveness

        for (int i = 0; i != N; ++i) {

          if (H(i,i) < PIVOT) {
            H(i,i) = PIVOT;
          }

          h[i] = 1.0 / sqrt(H(i,i));
        }


        // normalisation

        for (int i = 0; i != N; ++i) {
          for (int j = 0; j != i; ++j) {
            H(j,i) *= h[i] * h[j];
            H(i,j)  = H(j,i);
          }
        }
        
        for (int i = 0; i != N; ++i) {
          H(i,i) = 1.0 + lambda;
        }

        
        try {
          H.invert();
        }
        catch (const JException& error) {
          ERROR("JGandalf: " << error.what() << endl);
          return chi2_old;
        }
        
        
        for (int i = 0; i != N; ++i) {

          DEBUG("u[" << noshowpos << i << "] = " << showpos << FIXED(15,5) << value.*parameters[i]);

          for (int j = 0; j != N; ++j) {
            value.*parameters[i] -= H(i,j) * gradient.*parameters[j] * h[i] * h[j];
          }

          DEBUG(' ' << FIXED(15,5) << value.*parameters[i] << noshowpos << endl);

          error.*parameters[i] = h[i];
        }
      }

      return chi2_old;
    }


    /**
     * Multi-dimensional fit of one data set.
     *
     * The fit function should return the equivalent of chi2 for the current value
     * of the model and the given data point as well as the partial derivatives.
     *
     * \param  fit            fit function
     * \param  __begin        begin of data
     * \param  __end          end   of data
     * \return                chi2
     */
    template<class JFunction_t, class T>
    double operator()(const JFunction_t& fit, T __begin, T __end)
    {
      return (*this)(fit, __begin, __end, __end, __end);
    }


    static int     MAXIMUM_ITERATIONS;                         //!< maximal number of iterations
    static double  EPSILON;                                    //!< maximal distance to minimum
    static double  LAMBDA_MIN;                                 //!< minimal value control parameter
    static double  LAMBDA_MAX;                                 //!< maximal value control parameter
    static double  LAMBDA_UP;                                  //!< multiplication factor control parameter
    static double  LAMBDA_DOWN;                                //!< multiplication factor control parameter
    static double  PIVOT;                                      //!< minimal value diagonal element of matrix


    double                      lambda;
    JModel_t                    value;
    JModel_t                    error;
    std::vector<parameter_type> parameters;
    int                         numberOfIterations;
    JMATH::JMatrixNS            H;

  private:
    /**
     * Reset.
     */
    void reset()
    {
      chi2     = 0.0;      
      gradient = JModel_t();

      H.reset();          
    }


    /**
     * Evaluate fit for given data set.
     *
     * \param  fit            fit function
     * \param  __begin        begin of data
     * \param  __end          end   of data
     */
    template<class JFunction_t, class T>
    inline void evaluate(const JFunction_t& fit, T __begin, T __end)
    {
      for (T hit = __begin; hit != __end; ++hit) {
          
        const result_type& result = fit(value, *hit);

        chi2     += result.chi2;
        gradient += result.gradient;
        
        for (unsigned int i = 0; i != parameters.size(); ++i) {
          for (unsigned int j = i; j != parameters.size(); ++j) {
            H(i,j) += result.gradient.*parameters[i] * result.gradient.*parameters[j];
          }
        }
      }
    }

    double                      chi2;
    JModel_t                    gradient;
    JModel_t                    previous;
    std::vector<double>         h;                             // normalisation vector
  };


  /**
   * maximal number of iterations.
   */
  template<class JModel_t>
  int JGandalf<JModel_t>::MAXIMUM_ITERATIONS = 1000;


  /**
   * maximal distance to minimum.
   */
  template<class JModel_t>
  double JGandalf<JModel_t>::EPSILON         = 1.0e-4;


  /**
   * minimal value control parameter
   */
  template<class JModel_t>
  double JGandalf<JModel_t>::LAMBDA_MIN      =   0.01;


  /**
   * maximal value control parameter
   */
  template<class JModel_t>
  double JGandalf<JModel_t>::LAMBDA_MAX      = 100.0;


  /**
   * multiplication factor control parameter
   */
  template<class JModel_t>
  double JGandalf<JModel_t>::LAMBDA_UP       =   9.0;


  /**
   * multiplication factor control parameter
   */
  template<class JModel_t>
  double JGandalf<JModel_t>::LAMBDA_DOWN     =  11.0;


  /**
   * minimal value diagonal element of matrix
   */
  template<class JModel_t>
  double JGandalf<JModel_t>::PIVOT           =   1.0e-3;
}

#endif