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