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
      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>
    result_type 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;
 
      result_type precursor = result_type(numeric_limits<double>::max(), JModel_t());
 
      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) << successor.chi2 << endl);
 
        if (successor.chi2 < precursor.chi2) {
 
          if (numberOfIterations != 0) {
 
            if (fabs(precursor.chi2 - successor.chi2) < EPSILON*fabs(precursor.chi2)) {
              return successor;
            }
 
            if (lambda > LAMBDA_MIN) {
              lambda /= LAMBDA_DOWN;
            }
          }
 
          precursor = successor;
          previous  = value;
 
        } else {
          
          value   = previous;
          lambda *= LAMBDA_UP;
          
          if (lambda > LAMBDA_MAX) {
            return precursor;                                   // 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 precursor;
        }
        
        
        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) * successor.gradient.*parameters[j] * h[i] * h[j];
          }
          
          DEBUG(' ' << FIXED(15,5) << value.*parameters[i] << noshowpos << endl);
 
          error.*parameters[i] = h[i];
        }
 
      }
 
      return precursor;
    }
 
 
    /**
     * 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>
    result_type 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()
    {
      successor = result_type();
 
      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);
 
        successor.chi2     += result.chi2;
        successor.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];
          }
        }
      }
    }
 
    result_type                 successor;
    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-3;
 
 
  /**
   * 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