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 "JMath/JZero.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.\n
   * This data structure should have arithmetic capabilities.
   *
   * The data member JGandalf::value corresponds to the start c.q.\ final value of
   * the model of the fit procedure and JGandalf::error to the uncertainties.\n
   * The co-variance matrix is stored in data member JGandalf::H.\n
   *
   * The data member JGandalf::parameters constitutes a list of those parameters of the model that should actually be fitted.\n
   * For this, the model should contain the type definition for <tt>parameter_type</tt>.\n
   * Normally, this type definition corresponds to a pointer to a data member of the model.\n
   * Alternatively, the type definition can be <tt>size_t</tt> or <tt>int</tt>.\n
   * In that case, the model class should provide for the element access <tt>operator[]</tt>.\n
   *
   * The template fit function in JGandalf::operator() should provide for an implementation of the function.\n
   * This operator should return the data type JGandalf::result_type which is composed of
   * the value of the chi2 and the gradient of a data point, respectively.\n
   * The function operator returns the chi2 of the overall 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  __begin        begin of data
     * \param  __end          end   of data
     * \return                chi2
     */
    template<class JFunction_t, class T, class ...Args>
    result_type operator()(const JFunction_t& fit, T __begin, T __end, Args ...args)
    {
      using namespace std;
      using namespace JPP;

      const int N = parameters.size();

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

      previous  = value;

      successor.gradient = value;
      successor.gradient = zero;

      error     = value;
      error     = zero;

      lambda    = LAMBDA_MIN;

      result_type precessor = result_type(numeric_limits<double>::max(), value);

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

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

        reset();

        __evaluate__(fit, __begin, __end, args...);

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

        if (successor.chi2 < precessor.chi2) {

          if (numberOfIterations != 0) {

            if (fabs(precessor.chi2 - successor.chi2) < EPSILON*fabs(precessor.chi2)) {
              return successor;
            }

            if (lambda > LAMBDA_MIN) {
              lambda /= LAMBDA_DOWN;
            }
          }

          precessor = successor;
          previous  = value;

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

          reset();
          
          __evaluate__(fit, __begin, __end, args...);
        }

        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 precessor;
        }
        
        
        for (int i = 0; i != N; ++i) {

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

          for (int j = 0; j != N; ++j) {
            alias(value, parameters[i]) -= H(i,j) * alias(successor.gradient, parameters[j]) * h[i] * h[j];
          }
          
          DEBUG(' ' << FIXED(15,5) << alias(value, parameters[i]) << noshowpos << endl);

          alias(error, parameters[i]) = h[i];
        }
      }

      return precessor;
    }


    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()
    {
      using namespace std;
      using namespace JPP;

      successor.chi2     = 0.0;
      successor.gradient = zero;

      H.reset();          
    }


    /**
     * Recursive method for evaluation of fit.
     *
     * \param  fit            fit function
     * \param  __begin        begin of data
     * \param  __end          end   of data
     */
    template<class JFunction_t, class T, class ...Args>
    inline void __evaluate__(const JFunction_t& fit, T __begin, T __end, Args ...args)
    {
      using namespace std;

      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) += alias(result.gradient, parameters[i]) * alias(result.gradient, parameters[j]);
          }
        }
      }

      __evaluate__(fit, args...);
    }


    /**
     * Termination method for evaluation of fit.
     *
     * \param  fit            fit function
     */
    template<class JFunction_t>
    inline void __evaluate__(const JFunction_t& fit)
    {}


    /**
     * Read/write access to parameter value by data member.
     *
     * \param  model                    model
     * \param  parameter                parameter
     * \return                          value
     */
    template<class T>
    static inline double alias(const T& model, const typename T::parameter_type parameter)
    {
      return model.*parameter;
    }


    /**
     * Read/write access to parameter value by data member.
     *
     * \param  model                    model
     * \param  parameter                parameter
     * \return                          value
     */
    template<class T>
    static inline double& alias(T& model, const typename T::parameter_type parameter)
    {
      return model.*parameter;
    }


    /**
     * Read/write access to parameter value by index.
     *
     * \param  model                    model
     * \param  index                    index
     * \return                          value
     */
    template<class T>
    static inline double alias(const T& model, const size_t index)
    {
      return model[index];
    }


    /**
     * Read/write access to parameter value by index.
     *
     * \param  model                    model
     * \param  index                    index
     * \return                          value
     */
    template<class T>
    static inline double& alias(T& model, const size_t index)
    {
      return model[index];
    }


    /**
     * Read/write access to parameter value by index.
     *
     * \param  model                    model
     * \param  index                    index
     * \return                          value
     */
    template<class T>
    static inline double alias(const T& model, const int index)
    {
      return model[index];
    }


    /**
     * Read/write access to parameter value by index.
     *
     * \param  model                    model
     * \param  index                    index
     * \return                          value
     */
    template<class T>
    static inline double& alias(T& model, const int index)
    {
      return model[index];
    }


    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