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.
template<class T >
double	operator() (const T &getChi2)
	Fit.
template<class T >
void	operator() (std::ostream &out, const T &getChi2, const std::vector< double > &xs) const
	Make scan at current position.

Public Attributes
size_t	Nmax
	maximum number of iterations
size_t	Nextra
	maximum number of extra steps
double	epsilon
	epsilon
int	debug
	debug
bool	absolute = false
	absolute chi2
bool	normalise = false
	normalise gradient
size_t	numberOfIterations

Private Member Functions
template<class T >
double	evaluate (const T &getChi2)
	Evaluate gradient.
void	move (const double factor)
	Move.
template<class T >
void	scan (std::ostream &out, const T &getChi2, const std::vector< double > &xs, std::vector< double > &ys) const
	Scan.

Private Attributes
std::vector< double >	gradient

Detailed Description

Conjugate gradient fit.

Definition at line 74 of file JGradient.hh.

Constructor & Destructor Documentation

JGradient()

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 89 of file JGradient.hh.

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

Member Function Documentation

operator()() [1/2]

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 116 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];
      }
 
      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 * (absolute ? 1.0 : 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];
    }

operator()() [2/2]

template<class T >

void JFIT::JGradient::operator() ( std::ostream & out, const T & getChi2, const std::vector< double > & xs ) const

inline

Make scan at current position.

Parameters

out	output stream
getChi2	chi2 function
xs	set of steps

Definition at line 277 of file JGradient.hh.

    {
      using namespace std;
 
      out << "scan: ";
 
      for (const auto& i : *this) { 
        out << i.name << ' ';
      }
 
      out << "chi2" << endl;
 
      std::vector<double> ys;
 
      scan(out, getChi2, xs, ys);
    }

evaluate()

template<class T >

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

inlineprivate

Evaluate gradient.

Returns

chi2

Definition at line 312 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];
      }
 
      V = sqrt(V);
 
      DEBUG(setw(width) << left << "|gradient|" << right << ' ' << FIXED(12,5) << V << endl);
 
      if (normalise) {
        for (size_t i = 0; i != N; ++i) {
          gradient[i] /= V;
        }
      }
 
      return x0;
    }

move()

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

inlineprivate

Move.

Parameters

factor

factor

Definition at line 383 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);
        }
      }
    }

scan()

template<class T >

void JFIT::JGradient::scan ( std::ostream & out, const T & getChi2, const std::vector< double > & xs, std::vector< double > & ys ) const

inlineprivate

Scan.

Parameters

out	output stream
getChi2	chi2 function
xs	set of steps
ys	set of current steps

Definition at line 406 of file JGradient.hh.

    {
      using namespace std;
 
      if (ys.size() == this->size()) {
 
        for (size_t i = 0; i != this->size(); ++i) {
          (*this)[ i ]->apply(+ys[ i ]);
        }
 
        for (const double y : ys) {
          out << FIXED(12,5) << y << ' ';
        }
 
        out << FIXED(12,5) << getChi2(1) << endl; 
 
        for (size_t i = this->size(); i != 0; --i) {
          (*this)[i-1]->apply(-ys[i-1]);
        }
 
      } else {
 
        for (const double x : xs) {
 
          ys.push_back(x * (*this)[ys.size()].value);
 
          scan(out, getChi2, xs, ys);
 
          ys.pop_back();
        }
      }
    }

Member Data Documentation

Nmax

size_t JFIT::JGradient::Nmax

maximum number of iterations

Definition at line 295 of file JGradient.hh.

Nextra

size_t JFIT::JGradient::Nextra

maximum number of extra steps

Definition at line 296 of file JGradient.hh.

epsilon

double JFIT::JGradient::epsilon

epsilon

Definition at line 297 of file JGradient.hh.

debug

int JFIT::JGradient::debug

debug

Definition at line 298 of file JGradient.hh.

absolute

bool JFIT::JGradient::absolute = false

absolute chi2

Definition at line 300 of file JGradient.hh.

normalise

bool JFIT::JGradient::normalise = false

normalise gradient

Definition at line 301 of file JGradient.hh.

numberOfIterations

size_t JFIT::JGradient::numberOfIterations

Definition at line 303 of file JGradient.hh.

gradient

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

private

Definition at line 439 of file JGradient.hh.

---

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

• JGradient.hh