ROOT   Reference Guide
MnFunctionCross.cxx
Go to the documentation of this file.
1// @(#)root/minuit2:$Id$
2// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7 * *
8 **********************************************************************/
9
13#include "Minuit2/FCNBase.h"
14#include "Minuit2/MnParabola.h"
17#include "Minuit2/MnCross.h"
19
20//#define DEBUG
21#include "Minuit2/MnPrint.h"
22
23
24namespace ROOT {
25
26 namespace Minuit2 {
27
28
29
30MnCross MnFunctionCross::operator()(const std::vector<unsigned int>& par, const std::vector<double>& pmid,
31 const std::vector<double>& pdir, double tlr, unsigned int maxcalls) const {
32 // evaluate crossing point where function is equal to MIN + UP,
33 // with direction pdir from values pmid
34 // tlr indicate tolerance and maxcalls maximum number of calls
35
36// double edmmax = 0.5*0.001*toler*fFCN.Up();
37
38
39 unsigned int npar = par.size();
40 unsigned int nfcn = 0;
41 const MnMachinePrecision& prec = fState.Precision();
42 // tolerance used when calling Migrad
43 double mgr_tlr = 0.5 * tlr; // to be consistent with F77 version (for default values of tlr which is 0.1)
44 // other olerance values are fixed at 0.01
45 tlr = 0.01;
46 // convergence when F is within tlf of aim and next prediction
47 // of aopt is within tla of previous value of aopt
48 double up = fFCN.Up();
49 // for finding the point :
50 double tlf = tlr*up;
51 double tla = tlr;
52 unsigned int maxitr = 15;
53 unsigned int ipt = 0;
54 double aminsv = fFval;
55 double aim = aminsv + up;
56 //std::cout<<"aim= "<<aim<<std::endl;
57 double aopt = 0.;
58 bool limset = false;
59 std::vector<double> alsb(3, 0.), flsb(3, 0.);
60
61 int printLevel = MnPrint::Level();
62#ifdef DEBUG
63 printLevel = 3;
64#endif
65
66#ifdef DEBUG
67 for (unsigned int i = 0; i < par.size(); ++i)
68 std::cout << "MnFunctionCross for parameter " << par[i] << " value " << pmid[i] << " dir " << pdir[i]
69 << " function min = " << aminsv << " contur value aim = (fmin + up) = " << aim << std::endl;
70#endif
71
72
73 // find the largest allowed aulim
74
75 double aulim = 100.;
76 for(unsigned int i = 0; i < par.size(); i++) {
77 unsigned int kex = par[i];
78 if(fState.Parameter(kex).HasLimits()) {
79 double zmid = pmid[i];
80 double zdir = pdir[i];
81 // double zlim = 0.;
82 if(zdir > 0. && fState.Parameter(kex).HasUpperLimit()) {
83 double zlim = fState.Parameter(kex).UpperLimit();
84 if (fabs(zdir) < fState.Precision().Eps()) {
85 // we have a limit
86 if (fabs(zlim-zmid) < fState.Precision().Eps() ) limset = true;
87 continue;
88 }
89 aulim = std::min(aulim, (zlim-zmid)/zdir);
90 }
91 else if(zdir < 0. && fState.Parameter(kex).HasLowerLimit()) {
92 double zlim = fState.Parameter(kex).LowerLimit();
93 if (fabs(zdir) < fState.Precision().Eps()) {
94 // we have a limit
95 if (fabs(zlim-zmid) < fState.Precision().Eps() ) limset = true;
96 continue;
97 }
98 aulim = std::min(aulim, (zlim-zmid)/zdir);
99 }
100 }
101 }
102
103#ifdef DEBUG
104 std::cout<<"Largest allowed aulim "<< aulim << std::endl;
105#endif
106
107
108 // case of a single parameter and we are at limit
109 if (limset && npar == 1) {
110 if (printLevel > 1)
111 std::cout << "MnFunctionCross: parameter is at limit " << pmid[0] << " delta "
112 << pdir[0] << std::endl;
113 return MnCross(fState, nfcn, MnCross::CrossParLimit());
114 }
115
116 if (aulim < aopt + tla)
117 limset = true;
118
120
121 if (printLevel > 2) {
122 std::cout << "MnFunctionCross: Run Migrad fixing parameters :";
123 if (npar > 1) std::cout << std::endl;
124 }
125 for (unsigned int i = 0; i < npar; i++) {
126
128
129 if (printLevel > 2) {
130 std::cout << "\t" << fState.Name(par[i]) << " #" << par[i] << " to " << pmid[i] << std::endl;
131 }
132 }
133 // find minimum with respect all the other parameters (n- npar) (npar are the fixed ones)
134
135 FunctionMinimum min0 = migrad(maxcalls, mgr_tlr);
136 nfcn += min0.NFcn();
137
138 if (printLevel > 2) {
139 MnPrint::PrintState(std::cout, min0.State(), "MnFunctionCross: Result after Migrad");
140 std::cout << min0.UserState().Parameters() << std::endl;
141 }
142
143 // case a new minimum is found
144 if (min0.Fval() < fFval - tlf) {
145 // case of new minimum is found
146 if (printLevel > 1) {
147 std::cout << "MnFunctionCross: A new minimum is found when scanning parameter " << par.front() << " new value = "
148 << min0.Fval() << " old value : " << fFval << std::endl;
149 }
150 return MnCross(min0.UserState(), nfcn, MnCross::CrossNewMin());
151 }
152 if(min0.HasReachedCallLimit())
153 return MnCross(min0.UserState(), nfcn, MnCross::CrossFcnLimit());
154 if(!min0.IsValid()) return MnCross(fState, nfcn);
155 if(limset == true && min0.Fval() < aim)
156 return MnCross(min0.UserState(), nfcn, MnCross::CrossParLimit());
157
158 ipt++;
159 alsb[0] = 0.;
160 flsb[0] = min0.Fval();
161 flsb[0] = std::max(flsb[0], aminsv + 0.1*up);
162 aopt = sqrt(up/(flsb[0]-aminsv)) - 1.;
163 if(fabs(flsb[0] - aim) < tlf) return MnCross(aopt, min0.UserState(), nfcn);
164
165 if(aopt > 1.) aopt = 1.;
166 if(aopt < -0.5) aopt = -0.5;
167 limset = false;
168 if(aopt > aulim) {
169 aopt = aulim;
170 limset = true;
171 }
172#ifdef DEBUG
173 std::cout << "MnFunctionCross: flsb[0] = " << flsb[0] << " aopt = " << aopt << std::endl;
174#endif
175
176 if (printLevel > 2) {
177 std::cout << "MnFunctionCross: Run Migrad again (#2) : ";
178 if (npar > 1) std::cout << std::endl;
179 }
180 for (unsigned int i = 0; i < npar; i++) {
182 if (printLevel > 2) {
183 std::cout << "\t - parameter " << i << " fixed to " << pmid[i] + (aopt)*pdir[i] << std::endl;
184 }
185 }
186
187 FunctionMinimum min1 = migrad(maxcalls, mgr_tlr);
188 nfcn += min1.NFcn();
189
190 if (printLevel > 2) {
191 MnPrint::PrintState(std::cout, min1.State(), "MnFunctionCross: Result after 2nd Migrad");
192 std::cout << min1.UserState().Parameters() << std::endl;
193 }
194
195 if (min1.Fval() < fFval - tlf) // case of new minimum found
196 return MnCross(min1.UserState(), nfcn, MnCross::CrossNewMin());
197 if(min1.HasReachedCallLimit())
198 return MnCross(min1.UserState(), nfcn, MnCross::CrossFcnLimit());
199 if(!min1.IsValid()) return MnCross(fState, nfcn);
200 if(limset == true && min1.Fval() < aim)
201 return MnCross(min1.UserState(), nfcn, MnCross::CrossParLimit());
202
203 ipt++;
204 alsb[1] = aopt;
205 flsb[1] = min1.Fval();
206 double dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
207
208#ifdef DEBUG
209 std::cout << "aopt = " << aopt << " min1Val = " << flsb[1] << " dfda = " << dfda << std::endl;
210#endif
211
212
213L300:
214 if(dfda < 0.) {
215 // looking for slope of the right sign
216#ifdef DEBUG
217 std::cout << "MnFunctionCross: dfda < 0 - iterate from " << ipt << " to max of " << maxitr << std::endl;
218#endif
219 // iterate (max times is maxitr) incrementing aopt
220
221 unsigned int maxlk = maxitr - ipt;
222 for(unsigned int it = 0; it < maxlk; it++) {
223 alsb[0] = alsb[1];
224 flsb[0] = flsb[1];
225 // LM: Add + 1, looking at Fortran code it starts from 1 ( see bug #8396)
226 aopt = alsb[0] + 0.2*(it+1);
227 limset = false;
228 if(aopt > aulim) {
229 aopt = aulim;
230 limset = true;
231 }
232 if (printLevel > 2) {
233 std::cout << "MnFunctionCross: Run Migrad again (iteration " << it << " ) : ";
234 if (npar > 1) std::cout << std::endl;
235 }
236 for(unsigned int i = 0; i < npar; i++) {
238 if (printLevel > 2) {
239 std::cout << "\t - parameter " << i << " fixed to " << pmid[i] + (aopt)*pdir[i] << std::endl;
240 }
241 }
243 nfcn += min1.NFcn();
244
245 if (printLevel > 2) {
246 MnPrint::PrintState(std::cout, min1.State(), "MnFunctionCross: Result after Migrad");
247 std::cout << min1.UserState().Parameters() << std::endl;
248 }
249
250 if (min1.Fval() < fFval - tlf) // case of new minimum found
251 return MnCross(min1.UserState(), nfcn, MnCross::CrossNewMin());
252 if(min1.HasReachedCallLimit())
253 return MnCross(min1.UserState(), nfcn, MnCross::CrossFcnLimit());
254 if(!min1.IsValid()) return MnCross(fState, nfcn);
255 if(limset == true && min1.Fval() < aim)
256 return MnCross(min1.UserState(), nfcn, MnCross::CrossParLimit());
257 ipt++;
258 alsb[1] = aopt;
259 flsb[1] = min1.Fval();
260 dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
261 // if(dfda > 0.) goto L460;
262
263#ifdef DEBUG
264 std::cout << "aopt = " << aopt << " min1Val = " << flsb[1] << " dfda = " << dfda << std::endl;
265#endif
266
267 if(dfda > 0.) break;
268 }
269 if(ipt > maxitr) return MnCross(fState, nfcn);
270 } //if(dfda < 0.)
271
272L460:
273
274 // dfda > 0: we have two points with the right slope
275
276 aopt = alsb[1] + (aim-flsb[1])/dfda;
277
278#ifdef DEBUG
279 std::cout << "MnFunctionCross: dfda > 0 : aopt = " << aopt << std::endl;
280#endif
281
282 double fdist = std::min(fabs(aim - flsb[0]), fabs(aim - flsb[1]));
283 double adist = std::min(fabs(aopt - alsb[0]), fabs(aopt - alsb[1]));
284 tla = tlr;
285 if(fabs(aopt) > 1.) tla = tlr*fabs(aopt);
286 if(adist < tla && fdist < tlf) return MnCross(aopt, min1.UserState(), nfcn);
287 if(ipt > maxitr) return MnCross(fState, nfcn);
288 double bmin = std::min(alsb[0], alsb[1]) - 1.;
289 if(aopt < bmin) aopt = bmin;
290 double bmax = std::max(alsb[0], alsb[1]) + 1.;
291 if(aopt > bmax) aopt = bmax;
292
293 limset = false;
294 if(aopt > aulim) {
295 aopt = aulim;
296 limset = true;
297 }
298
299 if (printLevel > 2) {
300 std::cout << "MnFunctionCross: Run Migrad again (#3) : ";
301 if (npar > 1) std::cout << std::endl;
302 }
303 for(unsigned int i = 0; i < npar; i++) {
305 if (printLevel > 2) {
306 std::cout << "\t : parameter " << i << " fixed to " << pmid[i] + (aopt)*pdir[i] << std::endl;
307 }
308 }
309 FunctionMinimum min2 = migrad(maxcalls, mgr_tlr);
310 nfcn += min2.NFcn();
311
312 if (printLevel > 2) {
313 MnPrint::PrintState(std::cout, min2.State(), "MnFunctionCross: Result after Migrad (#3): ");
314 std::cout << min2.UserState().Parameters() << std::endl;
315 }
316
317 if (min2.Fval() < fFval - tlf) // case of new minimum found
318 return MnCross(min2.UserState(), nfcn, MnCross::CrossNewMin());
319 if(min2.HasReachedCallLimit())
320 return MnCross(min2.UserState(), nfcn, MnCross::CrossFcnLimit());
321 if(!min2.IsValid()) return MnCross(fState, nfcn);
322 if(limset == true && min2.Fval() < aim)
323 return MnCross(min2.UserState(), nfcn, MnCross::CrossParLimit());
324
325 ipt++;
326 alsb[2] = aopt;
327 flsb[2] = min2.Fval();
328
329 // now we have three points, ask how many < AIM
330
331 double ecarmn = fabs(flsb[2] - aim);
332 double ecarmx = 0.;
333 unsigned int ibest = 2;
334 unsigned int iworst = 0;
335 unsigned int noless = 0;
336
337 for(unsigned int i = 0; i < 3; i++) {
338 double ecart = fabs(flsb[i] - aim);
339 if(ecart > ecarmx) {
340 ecarmx = ecart;
341 iworst = i;
342 }
343 if(ecart < ecarmn) {
344 ecarmn = ecart;
345 ibest = i;
346 }
347 if(flsb[i] < aim) noless++;
348 }
349
350#ifdef DEBUG
351 std::cout << "MnFunctionCross: have three points : nless < aim = " << noless << " ibest = " << ibest << " iworst = " << iworst << std::endl;
352#endif
353
354 //std::cout<<"480"<<std::endl;
355
356 // at least one on each side of AIM (contour), fit a parabola
357 if(noless == 1 || noless == 2) goto L500;
358 // if all three are above AIM, third point must be the closest to AIM, return it
359 if(noless == 0 && ibest != 2) return MnCross(fState, nfcn);
360 // if all three below and third is not best then the slope has again gone negative,
361 // re-iterate and look for positive slope
362 if(noless == 3 && ibest != 2) {
363 alsb[1] = alsb[2];
364 flsb[1] = flsb[2];
365#ifdef DEBUG
366 std::cout << "MnFunctionCross: all three points below - look again fir positive slope " << std::endl;
367#endif
368 goto L300;
369 }
370
371 // in other case new straight line thru first two points
372
373 flsb[iworst] = flsb[2];
374 alsb[iworst] = alsb[2];
375 dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
376#ifdef DEBUG
377 std::cout << "MnFunctionCross: new straight line using point 1-2 - dfda = " << dfda << std::endl;
378#endif
379 goto L460;
380
381L500:
382
383 do {
384 // do parabola fit
385 MnParabola parbol = MnParabolaFactory()(MnParabolaPoint(alsb[0], flsb[0]), MnParabolaPoint(alsb[1], flsb[1]), MnParabolaPoint(alsb[2], flsb[2]));
386 // aopt = parbol.X_pos(aim);
387 //std::cout<<"alsb1,2,3= "<<alsb[0]<<", "<<alsb[1]<<", "<<alsb[2]<<std::endl;
388 //std::cout<<"flsb1,2,3= "<<flsb[0]<<", "<<flsb[1]<<", "<<flsb[2]<<std::endl;
389
390#ifdef DEBUG
391 std::cout << "MnFunctionCross: parabola fit: iteration " << ipt << std::endl;
392#endif
393
394 double coeff1 = parbol.C();
395 double coeff2 = parbol.B();
396 double coeff3 = parbol.A();
397 double determ = coeff2*coeff2 - 4.*coeff3*(coeff1 - aim);
398
399#ifdef DEBUG
400 std::cout << "MnFunctionCross: parabola fit: a = " << coeff3 << " b = "
401 << coeff2 << " c = " << coeff1 << " determ = " << determ << std::endl;
402#endif
403 // curvature is negative
404 if(determ < prec.Eps()) return MnCross(fState, nfcn);
405 double rt = sqrt(determ);
406 double x1 = (-coeff2 + rt)/(2.*coeff3);
407 double x2 = (-coeff2 - rt)/(2.*coeff3);
408 double s1 = coeff2 + 2.*x1*coeff3;
409 double s2 = coeff2 + 2.*x2*coeff3;
410
411#ifdef DEBUG
412 std::cout << "MnFunctionCross: parabola fit: x1 = " << x1 << " x2 = "
413 << x2 << " s1 = " << s1 << " s2 = " << s2 << std::endl;
414#endif
415
416#ifdef WARNINGMSG
417 if(s1*s2 > 0.) MN_INFO_MSG("MnFunctionCross problem 1");
418#endif
419 // find with root is the right one
420 aopt = x1;
421 double slope = s1;
422 if(s2 > 0.) {
423 aopt = x2;
424 slope = s2;
425 }
426#ifdef DEBUG
427 std::cout << "MnFunctionCross: parabola fit: aopt = " << aopt << " slope = "
428 << slope << std::endl;
429#endif
430
432 tla = tlr;
433 if(fabs(aopt) > 1.) tla = tlr*fabs(aopt);
434
435#ifdef DEBUG
436 std::cout << "MnFunctionCross: Delta(aopt) = " << fabs(aopt - alsb[ibest]) << " tla = "
437 << tla << "Delta(F) = " << fabs(flsb[ibest] - aim) << " tlf = " << tlf << std::endl;
438#endif
439
440
441 if(fabs(aopt - alsb[ibest]) < tla && fabs(flsb[ibest] - aim) < tlf)
442 return MnCross(aopt, min2.UserState(), nfcn);
443
444 // if(ipt > maxitr) return MnCross();
445
446 // see if proposed point is in acceptable zone between L and R
447 // first find ileft, iright, iout and ibest
448
449 unsigned int ileft = 3;
450 unsigned int iright = 3;
451 unsigned int iout = 3;
452 ibest = 0;
453 ecarmx = 0.;
454 ecarmn = fabs(aim-flsb[0]);
455 for(unsigned int i = 0; i < 3; i++) {
456 double ecart = fabs(flsb[i] - aim);
457 if(ecart < ecarmn) {
458 ecarmn = ecart;
459 ibest = i;
460 }
461 if(ecart > ecarmx) ecarmx = ecart;
462 if(flsb[i] > aim) {
463 if(iright == 3) iright = i;
464 else if(flsb[i] > flsb[iright]) iout = i;
465 else {
466 iout = iright;
467 iright = i;
468 }
469 } else if(ileft == 3) ileft = i;
470 else if(flsb[i] < flsb[ileft]) iout = i;
471 else {
472 iout = ileft;
473 ileft = i;
474 }
475 }
476
477#ifdef DEBUG
478 std::cout << "MnFunctionCross: ileft = " << ileft << " iright = "
479 << iright << " iout = " << iout << " ibest = " << ibest << std::endl;
480#endif
481
482 // avoid keeping a bad point nest time around
483
484 if(ecarmx > 10.*fabs(flsb[iout] - aim))
485 aopt = 0.5*(aopt + 0.5*(alsb[iright] + alsb[ileft]));
486
487 // knowing ileft and iright, get acceptable window
488 double smalla = 0.1*tla;
489 if(slope*smalla > tlf) smalla = tlf/slope;
490 double aleft = alsb[ileft] + smalla;
491 double aright = alsb[iright] - smalla;
492
493 // move proposed point AOPT into window if necessary
494 if(aopt < aleft) aopt = aleft;
495 if(aopt > aright) aopt = aright;
496 if(aleft > aright) aopt = 0.5*(aleft + aright);
497
498 // see if proposed point outside limits (should be impossible)
499 limset = false;
500 if(aopt > aulim) {
501 aopt = aulim;
502 limset = true;
503 }
504
505 // evaluate at new point aopt
506 if (printLevel > 2) {
507 std::cout << "MnFunctionCross: Run Migrad again at new point (#iter = " << ipt << " ): ";
508 if (npar > 1) std::cout << std::endl;
509 }
510 for(unsigned int i = 0; i < npar; i++) {
512 if (printLevel > 2) {
513 std::cout << "\t : parameter " << i << " fixed to " << pmid[i] + (aopt)*pdir[i] << std::endl;
514 }
515 }
517 nfcn += min2.NFcn();
518
519 if (printLevel > 2) {
520 MnPrint::PrintState(std::cout, min2.State(), "MnFunctionCross: Result after new Migrad : ");
521 std::cout << min2.UserState().Parameters() << std::endl;
522 }
523
524 if (min2.Fval() < fFval - tlf) // case of new minimum found
525 return MnCross(min2.UserState(), nfcn, MnCross::CrossNewMin());
526 if(min2.HasReachedCallLimit())
527 return MnCross(min2.UserState(), nfcn, MnCross::CrossFcnLimit());
528 if(!min2.IsValid()) return MnCross(fState, nfcn);
529 if(limset == true && min2.Fval() < aim)
530 return MnCross(min2.UserState(), nfcn, MnCross::CrossParLimit());
531
532 ipt++;
533 // replace odd point with new one (which is the best of three)
534 alsb[iout] = aopt;
535 flsb[iout] = min2.Fval();
536 ibest = iout;
537 } while(ipt < maxitr);
538
539 // goto L500;
540
541 return MnCross(fState, nfcn);
542}
543
544 } // namespace Minuit2
545
546} // namespace ROOT
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
#define s1(x)
Definition: RSha256.hxx:91
static const double x2[5]
static const double x1[5]
double sqrt(double)
virtual double Up() const =0
Error definition of the function.
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
const MnUserParameterState & UserState() const
const MinimumState & State() const
void SetValue(unsigned int, double)
const MnUserParameterState & fState
MnCross operator()(const std::vector< unsigned int > &, const std::vector< double > &, const std::vector< double > &, double, unsigned int) const
Sets the relative floating point (double) arithmetic precision.
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
A point of a parabola.
This class defines a parabola of the form a*x*x + b*x + c.
Definition: MnParabola.h:31
double B() const
Accessor to the coefficient of the linear term.
Definition: MnParabola.h:149
double C() const
Accessor to the coefficient of the constant term.
Definition: MnParabola.h:160
double A() const
Accessor to the coefficient of the quadratic term.
Definition: MnParabola.h:138
static void PrintState(std::ostream &os, const MinimumState &state, const char *msg, int iter=-1)
Definition: MnPrint.cxx:58
static int Level()
Definition: MnPrint.cxx:47
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
unsigned int Strategy() const
Definition: MnStrategy.h:39
const MnMachinePrecision & Precision() const
const MnUserParameters & Parameters() const
const MinuitParameter & Parameter(unsigned int i) const
const char * Name(unsigned int) const
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: StringConv.hxx:21