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 double 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 >
void	evaluate (const T &getChi2)
	Evaluate gradient. More...
void	move (const double factor)
	Move. More...

Private Attributes
std::vector< double >	gradient
double	chi2 [5]

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 double  debug = 3  )

inline

Constructor.

Parameters

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

Definition at line 84 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 syntax.

   double getChi2(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

Definition at line 110 of file JGradient.hh.

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

template<class T >

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

inlineprivate

Evaluate gradient.

Definition at line 271 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;
      }
        
      chi2[0] = getChi2(0);

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

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

inlineprivate

Move.

Parameters

factor

factor

Definition at line 332 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 261 of file JGradient.hh.

size_t JFIT::JGradient::Nextra

maximum number of extra steps

Definition at line 262 of file JGradient.hh.

double JFIT::JGradient::epsilon

epsilon

Definition at line 263 of file JGradient.hh.

int JFIT::JGradient::debug

debug

Definition at line 264 of file JGradient.hh.

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

private

Definition at line 345 of file JGradient.hh.

double JFIT::JGradient::chi2[5]

private

Definition at line 346 of file JGradient.hh.

---

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

• JGradient.hh