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 <type_traits>
 
#include "Jeep/JMessage.hh"
#include "JMath/JVectorND.hh"
#include "JMath/JMatrixNS.hh"
#include "JMath/JZero.hh"
#include "JLang/JManip.hh"
#include "JLang/JException.hh"
 
#include "Jeep/JMessage.hh"
 
/**
 * \author mdejong
 */
 
namespace JFIT {}
namespace JPP { using namespace JFIT; }
 
namespace JFIT {
 
  using JEEP::JMessage;
  using JLANG::JException;
 
  namespace JFIT_LOCAL {
 
    template<class T>
    class JTypedef {
      template<class U>   static auto parameter_type(U*)   -> decltype(std::declval<typename U::parameter_type>());
      template<typename>  static auto parameter_type(...)  -> std::false_type;
 
    public:
      static const bool has_parameter_type  = !std::is_same<std::false_type, decltype(parameter_type<T> (0))>::value;
    };
 
    template<class T, bool has_parameter_type = JTypedef<T>::has_parameter_type>
    struct JTypedef_t;
 
    template<class T> struct JTypedef_t<T, true>  { typedef typename T::parameter_type  parameter_type;  };
    template<class T> struct JTypedef_t<T, false> { typedef double   T::*parameter_type;  };
  }
  
 
  /**
   * Auxiliary function to constrain model during fit.
   *
   * \param  value          model (I/O)  
   */
  template<class JModel_t>
  inline void model(JModel_t& value)
  {}
 
 
  /**
   * 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::V.\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
   * If not defined, the parameters are assumed to be data members of type <tt>double</tt>.\n
   * Alternatively, the type definition can be <tt>size_t</tt> or <tt>int</tt>.\n
   * In that case, the model should provide for the element access <tt>operator[]</tt>.\n
   *
   * The first template parameter in the function operator should provide for an implementation of the actual fit function.\n
   * This function should return the data type JGandalf::result_type.\n
   * This data structure comprises the values of the chi2 and the gradient for a given data point.\n
   * The function operator returns the minimal chi2 and summed gradient of all data points.
   */
 
  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 JFIT_LOCAL::JTypedef_t<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;              //!< partial derivatives of chi2
    };
 
 
    /**
     * Default constructor.
     */
    JGandalf()
    {}
 
 
    /**
     * Multi-dimensional fit of multiple data sets.
     *
     * The fit function should return the chi2 as well as the partial derivatives
     * for the current value of the model and a given data point.
     *
     * \param  fit            fit function
     * \param  __begin        begin of data
     * \param  __end          end   of data
     * \param  args           optional data
     * \return                chi2 and gradient
     */
    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;
 
      // note that all model values should be assigned to the start value of the model before use
      // because the actual list of model parameters can vary from fit to fit
      // (e.g. if model consists of a container).
 
      const size_t N = parameters.size();
 
      V.resize(N);
      h.resize(N);
      x.resize(N);
 
      previous.result.chi2 = numeric_limits<double>::max();
 
      current.result.chi2     = numeric_limits<double>::max();
      current.result.gradient = value;
      current.result.gradient = zero;
 
      error     = value;
      error     = zero;
 
      lambda    = LAMBDA_MIN;
 
      for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
 
        DEBUG("step:     " << numberOfIterations << endl);
 
        reset();
 
        update(fit, __begin, __end, args...);
 
        DEBUG("lambda:   " << FIXED(12,5) << lambda              << endl);
        DEBUG("chi2:     " << FIXED(12,5) << current.result.chi2 << endl);
 
        if (current.result.chi2 < previous.result.chi2) {
 
          if (numberOfIterations != 0) {
 
            const double tolerance = EPSILON * (EPSILON_ABSOLUTE ? 1.0 : fabs(previous.result.chi2)); 
 
            if (fabs(previous.result.chi2 - current.result.chi2) <= tolerance) {
 
              // normal end
 
              const result_type result = current.result;
 
              reset();
 
              update(fit, __begin, __end, args...);
 
              try {
                V.invert();
              }
              catch (const exception& error) {}
 
              for (size_t i = 0; i != N; ++i) {
                get(error, parameters[i]) = sqrt(V(i,i));
              }
 
              return result;
            }
 
            if (lambda > LAMBDA_MIN) {
              lambda /= LAMBDA_DOWN;
            }
          }
 
          // store current values
 
          previous.value  = value;
          previous.result = current.result;
 
        } else {
 
          value   = previous.value;                                 // restore value
 
          lambda *= LAMBDA_UP;
          
          if (lambda > LAMBDA_MAX) {
            break;
          }
 
          reset();
 
          update(fit, __begin, __end, args...);
        }
 
        DEBUG("Hesse matrix:" << endl << V << endl);
 
        // force definite positiveness
 
        for (size_t i = 0; i != N; ++i) {
 
          if (V(i,i) < PIVOT) {
            V(i,i) = PIVOT;
          }
 
          h[i] = 1.0 / sqrt(V(i,i));
        }
 
        // normalisation
        
        for (size_t row = 0; row != N; ++row) {
          for (size_t col = 0; col != row; ++col) {
            V(row,col) *= h[row] * h[col];
            V(col,row)  = V(row,col);
          }
        }
        
        for (size_t i = 0; i != N; ++i) {
          V(i,i) = 1.0 + lambda;
        }
 
        // solve A x = b
 
        for (size_t col = 0; col != N; ++col) {
          x[col] = h[col] * get(current.result.gradient, parameters[col]);
        }
 
        try {
          V.solve(x);
        }
        catch (const exception& error) {
 
          ERROR("JGandalf: " << error.what() << endl << V << endl);
 
          break;
        }
 
        // update value
 
        for (size_t row = 0; row != N; ++row) {
 
          DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(15,5) << get(value, parameters[row]));
 
          get(value, parameters[row]) -= h[row] * x[row];
 
          DEBUG(" -> " << FIXED(15,5) << get(value, parameters[row]) << noshowpos << endl);
        }
 
        model(value);
      }
 
      // abnormal end
 
      const result_type result = previous.result;
 
      value   = previous.value;                                     // restore value
 
      reset();
 
      update(fit, __begin, __end, args...);
 
      try {
        V.invert();
      }
      catch (const exception& error) {}
 
      for (size_t i = 0; i != N; ++i) {
        get(error, parameters[i]) = sqrt(V(i,i));
      }
 
      return result;
    }
 
 
    static int     MAXIMUM_ITERATIONS;                         //!< maximal number of iterations
    static double  EPSILON;                                    //!< maximal distance to minimum
    static bool    EPSILON_ABSOLUTE;                           //!< set epsilon to absolute difference instead of relative
    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 Hesse matrix
 
    std::vector<parameter_type> parameters;                    //!< fit parameters
    int                         numberOfIterations;            //!< number of iterations
    double                      lambda;                        //!< control parameter
    JModel_t                    value;                         //!< value
    JModel_t                    error;                         //!< error
    JMATH::JMatrixNS            V;                             //!< Hesse matrix
 
  private:
    /**
     * Reset current parameters.
     */
    void reset()
    {
      using namespace JPP;
 
      current.result.chi2     = 0.0;
      current.result.gradient = zero;
 
      V.reset();          
    }
 
 
    /**
     * Recursive method to update current parameters.
     *
     * \param  fit            fit function
     * \param  __begin        begin of data
     * \param  __end          end   of data
     * \param  args           optional data
     */
    template<class JFunction_t, class T, class ...Args>
    inline void update(const JFunction_t& fit, T __begin, T __end, Args ...args)
    {
      for (T i = __begin; i != __end; ++i) {
 
        const result_type& result = fit(value, *i);
 
        current.result.chi2     += result.chi2;
        current.result.gradient += result.gradient;
 
        for (size_t row = 0; row != parameters.size(); ++row) {
          for (size_t col = row; col != parameters.size(); ++col) {
            V(row,col) += get(result.gradient, parameters[row]) * get(result.gradient, parameters[col]);
          }
        }
      }
 
      update(fit, args...);
    }
 
 
    /**
     * Termination method to update current parameters.
     *
     * \param  fit            fit function
     */
    template<class JFunction_t>
    inline void update(const JFunction_t& fit)
    {
      for (size_t row = 0; row != parameters.size(); ++row) {
        for (size_t col = 0; col != row; ++col) {
          V(row,col) = V(col,row);
        }
      }
    }
 
 
    /**
     * Read/write access to parameter value by data member.
     *
     * \param  model                    model
     * \param  parameter                parameter
     * \return                          value
     */
    static inline double get(const JModel_t& model, double JModel_t::*parameter)
    {
      return model.*parameter;
    }
 
 
    /**
     * Read/write access to parameter value by data member.
     *
     * \param  model                    model
     * \param  parameter                parameter
     * \return                          value
     */
    static inline double& get(JModel_t& model, double JModel_t::*parameter)
    {
      return model.*parameter;
    }
 
 
    /**
     * Read/write access to parameter value by index.
     *
     * \param  model                    model
     * \param  index                    index
     * \return                          value
     */
    static inline double get(const JModel_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
     */
    static inline double& get(JModel_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
     */
    static inline double get(const JModel_t& model, const int index)
    {
      return model[index];
    }
 
 
    /**
     * Read/write access to parameter value by index.
     *
     * \param  model                    model
     * \param  index                    index
     * \return                          value
     */
    static inline double& get(JModel_t& model, const int index)
    {
      return model[index];
    }
 
    std::vector<double>         h;                             // normalisation vector
    JMATH::JVectorND            x;
 
    struct {
      result_type               result;                        // result
    } current;
 
    struct {
      JModel_t                  value;                         // value
      result_type               result;                        // result
    } previous;
  };
 
 
  /**
   * 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;
  
  /**
   * set epsilon to absolute difference instead of relative.
   */
  template<class JModel_t>
  bool JGandalf<JModel_t>::EPSILON_ABSOLUTE  = false;
 
  /**
   * 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       =  10.0;
 
 
  /**
   * multiplication factor control parameter
   */
  template<class JModel_t>
  double JGandalf<JModel_t>::LAMBDA_DOWN     =  10.0;
 
 
  /**
   * minimal value diagonal element of matrix
   */
  template<class JModel_t>
  double JGandalf<JModel_t>::PIVOT           =   std::numeric_limits<double>::epsilon();
}
 
#endif