Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
24 namespace ROOT {
25 
26  namespace Minuit2 {
27 
28 
29 
30 MnCross 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 
178 L300:
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 
235 L460:
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 
342 L500:
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