JGradient.hh

Go to the documentation of this file.

#ifndef __JFIT__JGRADIENT__
#define __JFIT__JGRADIENT__
 
#include <limits>
#include <vector>
#include <cmath>
#include <memory>
#include <ostream>
#include <iomanip>
 
#include "JLang/JManip.hh"
#include "JLang/JException.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JFIT {}
namespace JPP { using namespace JFIT; }
 
namespace JFIT {
 
  /**
   * Auxiliary data structure for fit parameter.
   */
  struct JParameter_t {
    /**
     * Virtual destructor.
     */
    virtual ~JParameter_t()
    {}
 
    /**
     * Apply step.
     *
     * \param  step              step
     */
    virtual void apply(const double step) = 0;
  };
 
 
  /**
   * Auxiliary data structure for editable parameter.
   */
  struct JModifier_t :
    public std::shared_ptr<JParameter_t>
  {
    /**
     * Constructor.
     *
     * \param  name              name
     * \param  parameter         parameter
     * \param  value             value
     */
    JModifier_t(const std::string& name,
                JParameter_t*      parameter,
                const double       value) :
      std::shared_ptr<JParameter_t>(parameter),
      name (name),
      value(value)
    {}
 
    std::string name;
    double      value;
  };
  
 
  /**
   * Conjugate gradient fit.
   */
  struct JGradient :
    public std::vector<JModifier_t>
  {
    /**
     * Constructor.
     *
     * The number of iterations and epsilon refer to the number of steps and 
     * the distance to the minimum, respectively.\n
     * The number of extra steps can be used to overcome a possible hurdle on the way.
     *
     * \param  Nmax              maximum number of iterations
     * \param  Nextra            maximum number of extra steps
     * \param  epsilon           epsilon
     * \param  debug             debug
     */
    JGradient(const size_t Nmax     = std::numeric_limits<size_t>::max(),
              const size_t Nextra   = 0,
              const double epsilon  = 1.0e-4,
              const int    debug    = 3) :
      Nmax   (Nmax),
      Nextra (Nextra),
      epsilon(epsilon),
      debug  (debug)
    {}
 
 
    /**
     * Fit.
     *
     * The template parameter should provide for the following function operator.
     * <pre>
     *    double operator()(int option);
     * </pre>
     * The value of the option corresponds to the following cases.
     *  - 0 => step wise improvement of the chi2;
     *  - 1 => evaluation of the chi2 before the determination of the gradient of the chi2; and
     *  - 2 => evaluation of the derivative of the chi2 to each fit parameter.
     *
     * \param  getChi2           chi2 function
     * \return                   chi2
     */
    template<class T>
    double operator()(const T& getChi2)
    {
      using namespace std;
      using namespace JPP;
 
      vector<double> chi2(5, numeric_limits<double>::max());
 
      if (this->empty()) {
        return numeric_limits<double>::max();
      }
 
      chi2[0] = this->evaluate(getChi2);
      
      const size_t N = this->size();
 
      vector<double> H(N);
      vector<double> G(N);
 
      for (size_t i = 0; i != N; ++i) {
        G[i] = -1.0 * gradient[i];
        H[i] = G[i];
        gradient[i] = H[i];
      }
 
      numberOfIterations = 0;
 
      for ( ; numberOfIterations != Nmax; ++numberOfIterations) {
 
        DEBUG("chi2[0]  " << setw(4) << numberOfIterations << ' ' << FIXED(12,5) << chi2[0] << endl);
 
        // minimise chi2 in direction of gradient
 
        chi2[1] = chi2[0];
        chi2[2] = chi2[1];
 
        size_t m = 0;
        
        for (double ds = 1.0; ds > 1.0e-3; ) {
 
          this->move(+1.0 * ds);
 
          chi2[3] = getChi2(0);
 
          DEBUG("chi2[3]  " << setw(4) << m << ' ' << FIXED(12,5) << chi2[3] << ' ' << FIXED(12,5) << ds << endl);
 
          if (chi2[3] < chi2[2]) {
 
            chi2[1] = chi2[2];
            chi2[2] = chi2[3];
 
            m = 0;
            
            continue;
          }
 
          if (ds == 1.0) {
 
            if (m == 0) {
              chi2[4] = chi2[3];
            }
 
            if (m != Nextra) {
 
              ++m;
 
              continue;
 
            } else {
 
              for ( ; m != 0; --m) {
                this->move(-1.0 * ds);
              }
 
              chi2[3] = chi2[4];
            }
          }
 
          this->move(-1.0 * ds);
 
          if (chi2[2] < chi2[3]) {
 
            // final step based on parabolic interpolation through following points
            //
            // x1 = -1 * ds  ->  chi2[1]
            // x2 =  0 * ds  ->  chi2[2]
            // x3 = +1 * ds  ->  chi2[3]
 
            const double f21 = chi2[2] - chi2[1];   // f(x2) - f(x1)
            const double f23 = chi2[2] - chi2[3];   // f(x2) - f(x3)
 
            const double xs  =  0.5 * (f21 - f23) / (f23 + f21);
 
            this->move(+1.0 * xs * ds);
 
            chi2[3] = getChi2(0);
            
            if (chi2[3] < chi2[2]) {
 
              chi2[2] = chi2[3];
 
            } else {
 
              this->move(-1.0 * xs * ds);
 
              chi2[2] = getChi2(0);
            }
 
            DEBUG("chi2[2]  " << setw(4) << numberOfIterations << ' ' << FIXED(12,5) << chi2[2] << ' ' << SCIENTIFIC(12,5) << ds << endl);
 
            break;
 
          } else {
 
            ds *= 0.5;
          }
        }
 
        if (fabs(chi2[2] - chi2[0]) < epsilon * 0.5 * (fabs(chi2[0]) + fabs(chi2[2]))) {
 
          chi2[0] = chi2[2];
 
          break;
        }
 
        chi2[0] = this->evaluate(getChi2);
 
        double gg  = 0.0;
        double dgg = 0.0;
 
        for (size_t i = 0; i != N; ++i){
          gg  +=  G[i]*G[i];
          dgg += (gradient[i] + G[i]) * gradient[i];
        }
 
        if (gg == 0.0) {
          break;
        }
 
        dgg /= gg;
 
        for (size_t i = 0; i != N; ++i){
          G[i] = -1.0 * gradient[i];
          H[i] =  G[i] + dgg * H[i];
          gradient[i] = H[i];
        }
      }
 
      DEBUG("chi2[0]  " << setw(4) << numberOfIterations << ' ' << FIXED(12,5) << chi2[0] << endl);
 
      return chi2[0];
    }
 
 
    size_t  Nmax;              //!< maximum number of iterations
    size_t  Nextra;            //!< maximum number of extra steps
    double  epsilon;           //!< epsilon
    int     debug;             //!< debug
 
    size_t  numberOfIterations;
 
  private:
    /**
     * Evaluate gradient.
     *
     * \return                   chi2
     */
    template<class T>
    double evaluate(const T& getChi2)
    {
      using namespace std;
      using namespace JPP;
 
      const size_t N = this->size();
 
      gradient.resize(N);
 
      for (std::vector<double>::iterator i = gradient.begin(); i != gradient.end(); ++i) {
        *i = 0.0;
      }
      
      const double x0 = getChi2(1);
 
      size_t width = 1;
 
      for (size_t i = 0; i != N; ++i) {
        if ((*this)[i].name.size() > width) {
          width = (*this)[i].name.size();
        }
      }
 
      double V = 0.0;
 
      for (size_t i = 0; i != N; ++i) {
 
        if ((*this)[i].value != 0.0) {
 
          (*this)[i]->apply(+0.5 * (*this)[i].value);
 
          const double x1 = getChi2(2);
 
          (*this)[i]->apply(-0.5 * (*this)[i].value);
          (*this)[i]->apply(-0.5 * (*this)[i].value);
 
          const double x2 = getChi2(2);
 
          gradient[i] = x1 - x2;
        
          (*this)[i]->apply(+0.5 * (*this)[i].value);
 
          DEBUG(setw(width) << left << (*this)[i].name << right << ' ' << FIXED(12,5) << (*this)[i].value << ' ' << FIXED(12,5) << gradient[i] << endl);
          
        } else {
 
          gradient[i] = 0.0;
        }
 
        V += gradient[i] * gradient[i];
      }
 
      DEBUG(setw(width) << left << "|gradient|" << right << ' ' << FIXED(12,5) << sqrt(V) << endl);
 
      return x0;
    }
    
 
    /**
     * Move.
     *
     * \param  factor            factor
     */
    void move(const double factor)
    {
      if        (factor > 0.0) {
        for (size_t i = 0; i != this->size(); ++i) {
          (*this)[ i ]->apply((*this)[ i ].value * gradient[ i ] * factor);
        }
      } else if (factor < 0.0) {
        for (size_t i = this->size(); i != 0; --i) {
          (*this)[i-1]->apply((*this)[i-1].value * gradient[i-1] * factor);
        }
      }
    }
    
    std::vector<double> gradient;
  };
}
 
#endif