Jpp  19.1.0-rc.1
the software that should make you happy
Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
JFIT::JGradient Struct Reference

Conjugate gradient fit. More...

#include <JGradient.hh>

Inheritance diagram for JFIT::JGradient:
std::vector< JModifier_t >

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...
 
size_t numberOfIterations
 

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

◆ 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
Nmaxmaximum number of iterations
Nextramaximum number of extra steps
epsilonepsilon
debugdebug

Definition at line 88 of file JGradient.hh.

91  :
92  Nmax (Nmax),
93  Nextra (Nextra),
95  debug (debug)
96  {}
size_t Nmax
maximum number of iterations
Definition: JGradient.hh:268
int debug
debug
Definition: JGradient.hh:271
double epsilon
epsilon
Definition: JGradient.hh:270
size_t Nextra
maximum number of extra steps
Definition: JGradient.hh:269

Member Function Documentation

◆ operator()()

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
getChi2chi2 function
Returns
chi2

Definition at line 115 of file JGradient.hh.

116  {
117  using namespace std;
118  using namespace JPP;
119 
120  vector<double> chi2(5, numeric_limits<double>::max());
121 
122  if (this->empty()) {
123  return numeric_limits<double>::max();
124  }
125 
126  chi2[0] = this->evaluate(getChi2);
127 
128  const size_t N = this->size();
129 
130  vector<double> H(N);
131  vector<double> G(N);
132 
133  for (size_t i = 0; i != N; ++i) {
134  G[i] = -1.0 * gradient[i];
135  H[i] = G[i];
136  gradient[i] = H[i];
137  }
138 
139  numberOfIterations = 0;
140 
142 
143  DEBUG("chi2[0] " << setw(4) << numberOfIterations << ' ' << FIXED(12,5) << chi2[0] << endl);
144 
145  // minimise chi2 in direction of gradient
146 
147  chi2[1] = chi2[0];
148  chi2[2] = chi2[1];
149 
150  size_t m = 0;
151 
152  for (double ds = 1.0; ds > 1.0e-3; ) {
153 
154  this->move(+1.0 * ds);
155 
156  chi2[3] = getChi2(0);
157 
158  DEBUG("chi2[3] " << setw(4) << m << ' ' << FIXED(12,5) << chi2[3] << ' ' << FIXED(12,5) << ds << endl);
159 
160  if (chi2[3] < chi2[2]) {
161 
162  chi2[1] = chi2[2];
163  chi2[2] = chi2[3];
164 
165  m = 0;
166 
167  continue;
168  }
169 
170  if (ds == 1.0) {
171 
172  if (m == 0) {
173  chi2[4] = chi2[3];
174  }
175 
176  if (m != Nextra) {
177 
178  ++m;
179 
180  continue;
181 
182  } else {
183 
184  for ( ; m != 0; --m) {
185  this->move(-1.0 * ds);
186  }
187 
188  chi2[3] = chi2[4];
189  }
190  }
191 
192  this->move(-1.0 * ds);
193 
194  if (chi2[2] < chi2[3]) {
195 
196  // final step based on parabolic interpolation through following points
197  //
198  // x1 = -1 * ds -> chi2[1]
199  // x2 = 0 * ds -> chi2[2]
200  // x3 = +1 * ds -> chi2[3]
201 
202  const double f21 = chi2[2] - chi2[1]; // f(x2) - f(x1)
203  const double f23 = chi2[2] - chi2[3]; // f(x2) - f(x3)
204 
205  const double xs = 0.5 * (f21 - f23) / (f23 + f21);
206 
207  this->move(+1.0 * xs * ds);
208 
209  chi2[3] = getChi2(0);
210 
211  if (chi2[3] < chi2[2]) {
212 
213  chi2[2] = chi2[3];
214 
215  } else {
216 
217  this->move(-1.0 * xs * ds);
218 
219  chi2[2] = getChi2(0);
220  }
221 
222  DEBUG("chi2[2] " << setw(4) << numberOfIterations << ' ' << FIXED(12,5) << chi2[2] << ' ' << SCIENTIFIC(12,5) << ds << endl);
223 
224  break;
225 
226  } else {
227 
228  ds *= 0.5;
229  }
230  }
231 
232  if (fabs(chi2[2] - chi2[0]) < epsilon * 0.5 * (fabs(chi2[0]) + fabs(chi2[2]))) {
233 
234  chi2[0] = chi2[2];
235 
236  break;
237  }
238 
239  chi2[0] = this->evaluate(getChi2);
240 
241  double gg = 0.0;
242  double dgg = 0.0;
243 
244  for (size_t i = 0; i != N; ++i){
245  gg += G[i]*G[i];
246  dgg += (gradient[i] + G[i]) * gradient[i];
247  }
248 
249  if (gg == 0.0) {
250  break;
251  }
252 
253  dgg /= gg;
254 
255  for (size_t i = 0; i != N; ++i){
256  G[i] = -1.0 * gradient[i];
257  H[i] = G[i] + dgg * H[i];
258  gradient[i] = H[i];
259  }
260  }
261 
262  DEBUG("chi2[0] " << setw(4) << numberOfIterations << ' ' << FIXED(12,5) << chi2[0] << endl);
263 
264  return chi2[0];
265  }
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:56
static const double H
Planck constant [eV s].
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JSTDTypes.hh:14
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
void move(const double factor)
Move.
Definition: JGradient.hh:345
std::vector< double > gradient
Definition: JGradient.hh:358
double evaluate(const T &getChi2)
Evaluate gradient.
Definition: JGradient.hh:282
size_t numberOfIterations
Definition: JGradient.hh:273
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:488

◆ evaluate()

template<class T >
double JFIT::JGradient::evaluate ( const T &  getChi2)
inlineprivate

Evaluate gradient.

Returns
chi2

Definition at line 282 of file JGradient.hh.

283  {
284  using namespace std;
285  using namespace JPP;
286 
287  const size_t N = this->size();
288 
289  gradient.resize(N);
290 
291  for (std::vector<double>::iterator i = gradient.begin(); i != gradient.end(); ++i) {
292  *i = 0.0;
293  }
294 
295  const double x0 = getChi2(1);
296 
297  size_t width = 1;
298 
299  for (size_t i = 0; i != N; ++i) {
300  if ((*this)[i].name.size() > width) {
301  width = (*this)[i].name.size();
302  }
303  }
304 
305  double V = 0.0;
306 
307  for (size_t i = 0; i != N; ++i) {
308 
309  if ((*this)[i].value != 0.0) {
310 
311  (*this)[i]->apply(+0.5 * (*this)[i].value);
312 
313  const double x1 = getChi2(2);
314 
315  (*this)[i]->apply(-0.5 * (*this)[i].value);
316  (*this)[i]->apply(-0.5 * (*this)[i].value);
317 
318  const double x2 = getChi2(2);
319 
320  gradient[i] = x1 - x2;
321 
322  (*this)[i]->apply(+0.5 * (*this)[i].value);
323 
324  DEBUG(setw(width) << left << (*this)[i].name << right << ' ' << FIXED(12,5) << (*this)[i].value << ' ' << FIXED(12,5) << gradient[i] << endl);
325 
326  } else {
327 
328  gradient[i] = 0.0;
329  }
330 
331  V += gradient[i] * gradient[i];
332  }
333 
334  DEBUG(setw(width) << left << "|gradient|" << right << ' ' << FIXED(12,5) << sqrt(V) << endl);
335 
336  return x0;
337  }

◆ move()

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

Move.

Parameters
factorfactor

Definition at line 345 of file JGradient.hh.

346  {
347  if (factor > 0.0) {
348  for (size_t i = 0; i != this->size(); ++i) {
349  (*this)[ i ]->apply((*this)[ i ].value * gradient[ i ] * factor);
350  }
351  } else if (factor < 0.0) {
352  for (size_t i = this->size(); i != 0; --i) {
353  (*this)[i-1]->apply((*this)[i-1].value * gradient[i-1] * factor);
354  }
355  }
356  }

Member Data Documentation

◆ Nmax

size_t JFIT::JGradient::Nmax

maximum number of iterations

Definition at line 268 of file JGradient.hh.

◆ Nextra

size_t JFIT::JGradient::Nextra

maximum number of extra steps

Definition at line 269 of file JGradient.hh.

◆ epsilon

double JFIT::JGradient::epsilon

epsilon

Definition at line 270 of file JGradient.hh.

◆ debug

int JFIT::JGradient::debug

debug

Definition at line 271 of file JGradient.hh.

◆ numberOfIterations

size_t JFIT::JGradient::numberOfIterations

Definition at line 273 of file JGradient.hh.

◆ gradient

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

Definition at line 358 of file JGradient.hh.


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