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