Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
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.
 
template<class T >
double operator() (const T &getChi2)
 Fit.
 

Public Attributes

size_t Nmax
 maximum number of iterations
 
size_t Nextra
 maximum number of extra steps
 
double epsilon
 epsilon
 
int debug
 debug
 
size_t numberOfIterations
 

Private Member Functions

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

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
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
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.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
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: