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
     */
    template<class T>
    double operator()(const T& getChi2)
    {
      using namespace std;
      using namespace JPP;

      if (this->empty()) {
        return numeric_limits<double>::max();
      }

      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];
      }

      size_t number_of_iterations = 0;

      for ( ; number_of_iterations != Nmax; ++number_of_iterations) {

        DEBUG("chi2[0]  " << setw(4) << number_of_iterations << ' ' << 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) << number_of_iterations << ' ' << 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;
        }

        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) << number_of_iterations << ' ' << 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

  private:
    /**
     * Evaluate gradient.
     */
    template<class T>
    void 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;
      }
        
      chi2[0] = 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);

          chi2[1] = getChi2(2);

          (*this)[i]->apply(-0.5 * (*this)[i].value);
          (*this)[i]->apply(-0.5 * (*this)[i].value);

          chi2[2] = getChi2(2);

          gradient[i] = chi2[1] - chi2[2];
        
          (*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);
    }
    

    /**
     * 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;
    double              chi2[5];
  };
}

#endif