JFIT::JGradient Struct Reference

Conjugate gradient fit. More...

#include <JGradient.hh>

Inheritance diagram for JFIT::JGradient:

Public Member Functions
	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)
	Constructor. More...
template<class T >
double	operator() (const T &getChi2)
	Fit. More...

Public Attributes
size_t	Nmax
	maximum number of iterations More...
size_t	Nextra
	maximum number of extra steps More...
double	epsilon
	epsilon More...
int	debug
	debug More...

Private Member Functions
template<class T >
double	evaluate (const T &getChi2)
	Evaluate gradient. More...
void	move (const double factor)
	Move. More...

Private Attributes
std::vector< double >	gradient

Detailed Description

Conjugate gradient fit.

Definition at line 73 of file JGradient.hh.

Constructor & Destructor Documentation

JFIT::JGradient::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  )

inline

Constructor.

The number of iterations and epsilon refer to the number of steps and the distance to the minimum, respectively.
The number of extra steps can be used to overcome a possible hurdle on the way.

Parameters

Nmax	maximum number of iterations
Nextra	maximum number of extra steps
epsilon	epsilon
debug	debug

Definition at line 88 of file JGradient.hh.

                                         :
      Nmax   (Nmax),
      Nextra (Nextra),
      epsilon(epsilon),
      debug  (debug)
    {}

Member Function Documentation

template<class T >

double JFIT::JGradient::operator() ( const T &  getChi2 )

inline

Fit.

The template parameter should provide for the following function operator.

   double operator()(int option);

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.

Parameters

getChi2

chi2 function

Returns

chi2

Definition at line 115 of file JGradient.hh.

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

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

        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) << number_of_iterations << ' ' << FIXED(12,5) << chi2[0] << endl);

      return chi2[0];
    }

template<class T >

double JFIT::JGradient::evaluate ( const T &  getChi2 )

inlineprivate

Evaluate gradient.

Returns

chi2

Definition at line 280 of file JGradient.hh.

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

void JFIT::JGradient::move ( const double  factor )

inlineprivate

Move.

Parameters

factor

factor

Definition at line 343 of file JGradient.hh.

    {
      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);
        }
      }
    }

Member Data Documentation

size_t JFIT::JGradient::Nmax

maximum number of iterations

Definition at line 268 of file JGradient.hh.

size_t JFIT::JGradient::Nextra

maximum number of extra steps

Definition at line 269 of file JGradient.hh.

double JFIT::JGradient::epsilon

epsilon

Definition at line 270 of file JGradient.hh.

int JFIT::JGradient::debug

debug

Definition at line 271 of file JGradient.hh.

std::vector<double> JFIT::JGradient::gradient

private

Definition at line 356 of file JGradient.hh.

---

The documentation for this struct was generated from the following file:

• JGradient.hh