Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
MnLineSearch.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/MnLineSearch.h"
11 #include "Minuit2/MnFcn.h"
12 #include "Minuit2/MinimumParameters.h"
13 #include "Minuit2/MnMachinePrecision.h"
14 #include "Minuit2/MnParabola.h"
15 #include "Minuit2/MnParabolaPoint.h"
16 #include "Minuit2/MnParabolaFactory.h"
17 #include "Minuit2/LaSum.h"
18 
19 #include <iostream>
20 #include "Minuit2/MnPrint.h"
21 #include <algorithm>
22 
23 #ifdef USE_OTHER_LS
24 
25 #include "Math/SMatrix.h"
26 #include "Math/SVector.h"
27 
28 #include "Math/IFunction.h"
29 #include "Math/Minimizer1D.h"
30 
31 #endif
32 
33 //#define DEBUG
34 
35 namespace ROOT {
36 
37  namespace Minuit2 {
38 
39 
40 /** Perform a line search from position defined by the vector st
41  along the direction step, where the length of vector step
42  gives the expected position of Minimum.
43  fcn is Value of function at the starting position ,
44  gdel (if non-zero) is df/dx along step at st.
45  Return a parabola point containing Minimum x position and y (function Value)
46  - add a falg to control the debug
47 */
48 
49 MnParabolaPoint MnLineSearch::operator()(const MnFcn& fcn, const MinimumParameters& st, const MnAlgebraicVector& step, double gdel, const MnMachinePrecision& prec, bool debug) const {
50 
51  //*-*-*-*-*-*-*-*-*-*Perform a line search from position st along step *-*-*-*-*-*-*-*
52  //*-* =========================================
53  //*-* SLAMBG and ALPHA control the maximum individual steps allowed.
54  //*-* The first step is always =1. The max length of second step is SLAMBG.
55  //*-* The max size of subsequent steps is the maximum previous successful
56  //*-* step multiplied by ALPHA + the size of most recent successful step,
57  //*-* but cannot be smaller than SLAMBG.
58  //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
59 
60 #ifdef DEBUG
61  debug = true;
62 #endif
63 
64  if (debug) {
65  std::cout<<"gdel= "<<gdel<<std::endl;
66  std::cout<<"step= "<<step<<std::endl;
67  }
68 
69  double overal = 1000.;
70  double undral = -100.;
71  double toler = 0.05;
72  double slamin = 0.;
73  double slambg = 5.;
74  double alpha = 2.;
75  int maxiter = 12;
76  // start as in Fortran from 1 and count all the time we evaluate the function
77  int niter = 1;
78 
79  for(unsigned int i = 0; i < step.size(); i++) {
80  if(step(i) == 0 ) continue;
81  double ratio = fabs(st.Vec()(i)/step(i));
82  if( slamin == 0) slamin = ratio;
83  if(ratio < slamin) slamin = ratio;
84  }
85  if(fabs(slamin) < prec.Eps()) slamin = prec.Eps();
86  slamin *= prec.Eps2();
87 
88  double f0 = st.Fval();
89  double f1 = fcn(st.Vec()+step);
90  niter++;
91  double fvmin = st.Fval();
92  double xvmin = 0.;
93 
94  if(f1 < f0) {
95  fvmin = f1;
96  xvmin = 1.;
97  }
98  double toler8 = toler;
99  double slamax = slambg;
100  double flast = f1;
101  double slam = 1.;
102 
103  bool iterate = false;
104  MnParabolaPoint p0(0., f0);
105  MnParabolaPoint p1(slam, flast);
106  double f2 = 0.;
107  // quadratic interpolation using the two points p0,p1 and the slope at p0
108  do {
109  // cut toler8 as function goes up
110  iterate = false;
111  //MnParabola pb = MnParabolaFactory()(p0, gdel, p1);
112 
113  if (debug) {
114  // std::cout<<"pb.Min() = "<<pb.Min()<<std::endl;
115  std::cout<<"flast, f0= "<<flast<<", "<<f0<<std::endl;
116  std::cout<<"flast-f0= "<<flast-f0<<std::endl;
117  std::cout<<"slam= "<<slam<<std::endl;
118  }
119  // double df = flast-f0;
120  // if(fabs(df) < prec.Eps2()) {
121  // if(flast-f0 < 0.) df = -prec.Eps2();
122  // else df = prec.Eps2();
123  // }
124  // std::cout<<"df= "<<df<<std::endl;
125  // double denom = 2.*(df-gdel*slam)/(slam*slam);
126  double denom = 2.*(flast-f0-gdel*slam)/(slam*slam);
127  if (debug) std::cout<<"denom= "<<denom<<std::endl;
128  if(denom != 0) {
129  slam = - gdel/denom;
130  } else {
131  denom = -0.1*gdel;
132  slam = 1.;
133  }
134  if (debug) std::cout<<"new slam= "<<slam<<std::endl;
135 
136 
137 #ifdef TRY_OPPOSIT_DIR
138  // if is less than zero indicates maximum position. Use then slamax or x = x0 - 2slam + 1
139  bool slamIsNeg = false;
140  double slamNeg = 0;
141 #endif
142 
143  if(slam < 0.) {
144  if (debug) std::cout<<"slam is negative- set to " << slamax << std::endl;
145 #ifdef TRY_OPPOSITE_DIR
146  slamNeg = 2*slam -1;
147  slamIsNeg = true;
148  if (debug) std::cout<<"slam is negative- compare values between "<< slamNeg << " and " << slamax << std::endl;
149 #endif
150  slam = slamax;
151  }
152  // std::cout<<"slam= "<<slam<<std::endl;
153  if(slam > slamax) {
154  slam = slamax;
155  if (debug) std::cout << "slam larger than mac value - set to " << slamax << std::endl;
156  }
157 
158  if(slam < toler8) {
159  if (debug) std::cout << "slam too small - set to " << toler8 << std::endl;
160  slam = toler8;
161  }
162  // std::cout<<"slam= "<<slam<<std::endl;
163  if(slam < slamin) {
164  if (debug) std::cout << "slam smaller than " << slamin << " return " << std::endl;
165  // std::cout<<"f1, f2= "<<p0.Y()<<", "<<p1.Y()<<std::endl;
166  // std::cout<<"x1, x2= "<<p0.X()<<", "<<p1.X()<<std::endl;
167  // std::cout<<"x, f= "<<xvmin<<", "<<fvmin<<std::endl;
168  return MnParabolaPoint(xvmin, fvmin);
169  }
170  if(fabs(slam - 1.) < toler8 && p1.Y() < p0.Y()) {
171  // std::cout<<"f1, f2= "<<p0.Y()<<", "<<p1.Y()<<std::endl;
172  // std::cout<<"x1, x2= "<<p0.X()<<", "<<p1.X()<<std::endl;
173  // std::cout<<"x, f= "<<xvmin<<", "<<fvmin<<std::endl;
174  return MnParabolaPoint(xvmin, fvmin);
175  }
176  if(fabs(slam - 1.) < toler8) slam = 1. + toler8;
177 
178  // if(fabs(gdel) < prec.Eps2() && fabs(denom) < prec.Eps2())
179  // slam = 1000.;
180  // MnAlgebraicVector tmp = step;
181  // tmp *= slam;
182  // f2 = fcn(st.Vec()+tmp);
183  f2 = fcn(st.Vec() + slam*step);
184 
185  niter++; // do as in Minuit (count all func evalu)
186 
187 #ifdef TRY_OPPOSITE_DIR
188  if (slamIsNeg) {
189  // try alternative in the opposite direction
190  double f2alt = fcn(st.Vec() + slamNeg*step);
191  if (f2alt < f2) {
192  slam = slamNeg;
193  f2 = f2alt;
194  undral += slam;
195  }
196  }
197 #endif
198  if(f2 < fvmin) {
199  fvmin = f2;
200  xvmin = slam;
201  }
202  // LM : correct a bug using precision
203  if (fabs( p0.Y() - fvmin) < fabs(fvmin)*prec.Eps() ) {
204  // if(p0.Y()-prec.Eps() < fvmin && fvmin < p0.Y()+prec.Eps()) {
205  iterate = true;
206  flast = f2;
207  toler8 = toler*slam;
208  overal = slam - toler8;
209  slamax = overal;
210  p1 = MnParabolaPoint(slam, flast);
211  //niter++;
212  }
213  } while(iterate && niter < maxiter);
214  if(niter >= maxiter) {
215  // exhausted max number of iterations
216  return MnParabolaPoint(xvmin, fvmin);
217  }
218 
219  if (debug){
220  std::cout<<"after initial 2-point iter: "<<std::endl;
221  std::cout<<"x0, x1, x2= "<<p0.X()<<", "<<p1.X()<<", "<<slam<<std::endl;
222  std::cout<<"f0, f1, f2= "<<p0.Y()<<", "<<p1.Y()<<", "<<f2<<std::endl;
223  }
224 
225  MnParabolaPoint p2(slam, f2);
226 
227  // do now the quadratic interpolation with 3 points
228  do {
229  slamax = std::max(slamax, alpha*fabs(xvmin));
230  MnParabola pb = MnParabolaFactory()(p0, p1, p2);
231  if (debug) {
232  std::cout << "\nLS Iteration " << niter << std::endl;
233  std::cout<<"x0, x1, x2= "<<p0.X()<<", "<<p1.X()<<", "<<p2.X() <<std::endl;
234  std::cout<<"f0, f1, f2= "<<p0.Y()<<", "<<p1.Y()<<", "<<p2.Y() <<std::endl;
235  std::cout << "slamax = " << slamax << std::endl;
236  std::cout<<"p2-p0,p1: "<<p2.Y() - p0.Y()<<", "<<p2.Y() - p1.Y()<<std::endl;
237  std::cout<<"a, b, c= "<<pb.A()<<" "<<pb.B()<<" "<<pb.C()<<std::endl;
238  }
239  if(pb.A() < prec.Eps2()) {
240  double slopem = 2.*pb.A()*xvmin + pb.B();
241  if(slopem < 0.) slam = xvmin + slamax;
242  else slam = xvmin - slamax;
243  if (debug) std::cout << "xvmin = " << xvmin << " slopem " << slopem << " slam " << slam << std::endl;
244  } else {
245  slam = pb.Min();
246  // std::cout<<"pb.Min() slam= "<<slam<<std::endl;
247  if(slam > xvmin + slamax) slam = xvmin + slamax;
248  if(slam < xvmin - slamax) slam = xvmin - slamax;
249  }
250  if(slam > 0.) {
251  if(slam > overal) slam = overal;
252  } else {
253  if(slam < undral) slam = undral;
254  }
255 
256  if (debug) {
257  std::cout<<" slam= "<<slam<< " undral " << undral << " overal " << overal << std::endl;
258  }
259 
260  double f3 = 0.;
261  do {
262 
263  if (debug) std::cout << " iterate on f3- slam " << niter << " slam = " << slam << " xvmin " << xvmin << std::endl;
264 
265  iterate = false;
266  double toler9 = std::max(toler8, fabs(toler8*slam));
267  // min. of parabola at one point
268  if(fabs(p0.X() - slam) < toler9 ||
269  fabs(p1.X() - slam) < toler9 ||
270  fabs(p2.X() - slam) < toler9) {
271  // std::cout<<"f1, f2, f3= "<<p0.Y()<<", "<<p1.Y()<<", "<<p2.Y()<<std::endl;
272  // std::cout<<"x1, x2, x3= "<<p0.X()<<", "<<p1.X()<<", "<<p2.X()<<std::endl;
273  // std::cout<<"x, f= "<<xvmin<<", "<<fvmin<<std::endl;
274  return MnParabolaPoint(xvmin, fvmin);
275  }
276 
277  // take the step
278  // MnAlgebraicVector tmp = step;
279  // tmp *= slam;
280  f3 = fcn(st.Vec() + slam*step);
281  if (debug) {
282  std::cout<<"f3= "<<f3<<std::endl;
283  std::cout<<"f3-p(2-0).Y()= "<<f3-p2.Y()<<" "<<f3-p1.Y()<<" "<<f3-p0.Y()<<std::endl;
284  }
285  // if latest point worse than all three previous, cut step
286  if(f3 > p0.Y() && f3 > p1.Y() && f3 > p2.Y()) {
287  if (debug) {
288  std::cout<<"f3 worse than all three previous"<<std::endl;
289  }
290  if(slam > xvmin) overal = std::min(overal, slam-toler8);
291  if(slam < xvmin) undral = std::max(undral, slam+toler8);
292  slam = 0.5*(slam + xvmin);
293  if (debug) {
294  std::cout<<"new slam= "<<slam<<std::endl;
295  }
296  iterate = true;
297  niter++;
298  }
299  } while(iterate && niter < maxiter);
300  if(niter >= maxiter) {
301  // exhausted max number of iterations
302  return MnParabolaPoint(xvmin, fvmin);
303  }
304 
305  // find worst previous point out of three and replace
306  MnParabolaPoint p3(slam, f3);
307  if(p0.Y() > p1.Y() && p0.Y() > p2.Y()) p0 = p3;
308  else if(p1.Y() > p0.Y() && p1.Y() > p2.Y()) p1 = p3;
309  else p2 = p3;
310  if (debug) std::cout << " f3 " << f3 << " fvmin " << fvmin << " xvmin " << xvmin << std::endl;
311  if(f3 < fvmin) {
312  fvmin = f3;
313  xvmin = slam;
314  } else {
315  if(slam > xvmin) overal = std::min(overal, slam-toler8);
316  if(slam < xvmin) undral = std::max(undral, slam+toler8);
317  }
318 
319  niter++;
320  } while(niter < maxiter);
321 
322  if (debug) {
323  std::cout<<"f1, f2= "<<p0.Y()<<", "<<p1.Y()<<std::endl;
324  std::cout<<"x1, x2= "<<p0.X()<<", "<<p1.X()<<std::endl;
325  std::cout<<"x, f= "<<xvmin<<", "<<fvmin<<std::endl;
326  }
327  return MnParabolaPoint(xvmin, fvmin);
328 
329 }
330 
331 #ifdef USE_OTHER_LS
332 
333 /** Perform a line search using a cubic interpolation using x0, x1 , df/dx(x0) and d2/dx(x0) (second derivative)
334  This is used at the beginning when the second derivative is known to be negative
335 */
336 
337 MnParabolaPoint MnLineSearch::CubicSearch(const MnFcn& fcn, const MinimumParameters& st, const MnAlgebraicVector& step, double gdel, double g2del, const MnMachinePrecision& prec, bool debug) const {
338 
339 
340 
341  if (debug) {
342  std::cout<<"gdel= "<<gdel<<std::endl;
343  std::cout<<"g2del= "<<g2del<<std::endl;
344  std::cout<<"step= "<<step<<std::endl;
345  }
346 
347  // change ot large values
348  double overal = 100.;
349  double undral = -100.;
350  double toler = 0.05;
351  double slamin = 0.;
352  double slambg = 5.;
353  double alpha = 2.;
354 
355  for(unsigned int i = 0; i < step.size(); i++) {
356  if(step(i) == 0 ) continue;
357  double ratio = fabs(st.Vec()(i)/step(i));
358  if( slamin == 0) slamin = ratio;
359  if(ratio < slamin) slamin = ratio;
360  }
361  if(fabs(slamin) < prec.Eps()) slamin = prec.Eps();
362  slamin *= prec.Eps2();
363 
364  double f0 = st.Fval();
365  double f1 = fcn(st.Vec()+step);
366  double fvmin = st.Fval();
367  double xvmin = 0.;
368  if (debug) std::cout << "f0 " << f0 << " f1 " << f1 << std::endl;
369 
370  if(f1 < f0) {
371  fvmin = f1;
372  xvmin = 1.;
373  }
374  double toler8 = toler;
375  double slamax = slambg;
376  double flast = f1;
377  double slam = 1.;
378 
379 // MnParabolaPoint p0(0., f0);
380 // MnParabolaPoint p1(slam, flast);
381 
382  ROOT::Math::SMatrix<double, 3> cubicMatrix;
383  ROOT::Math::SVector<double, 3> cubicCoeff; // cubic coefficients to be found
384  ROOT::Math::SVector<double, 3> bVec; // cubic coefficients to be found
385  double x0 = 0;
386 
387  // cubic interpolation using the two points p0,p1 and the slope at p0 and the second derivative at p0
388 
389  // cut toler8 as function goes up
390  double x1 = slam;
391  cubicMatrix(0,0) = (x0*x0*x0 - x1*x1*x1)/3.;
392  cubicMatrix(0,1) = (x0*x0 - x1*x1)/2.;
393  cubicMatrix(0,2) = (x0 - x1);
394  cubicMatrix(1,0) = x0*x0;
395  cubicMatrix(1,1) = x0;
396  cubicMatrix(1,2) = 1;
397  cubicMatrix(2,0) = 2.*x0;
398  cubicMatrix(2,1) = 1;
399  cubicMatrix(2,2) = 0;
400 
401  bVec(0) = f0-f1;
402  bVec(1) = gdel;
403  bVec(2) = g2del;
404 
405  //if (debug) std::cout << "Matrix:\n " << cubicMatrix << std::endl;
406  if (debug) std::cout << "Vec:\n " << bVec << std::endl;
407 
408  // find the minimum need to invert the matrix
409  if (!cubicMatrix.Invert() ) {
410  std::cout << "Inversion failed - return " << std::endl;
411  return MnParabolaPoint(xvmin, fvmin);
412  }
413 
414  cubicCoeff = cubicMatrix * bVec;
415  if (debug) std::cout << "Cubic:\n " << cubicCoeff << std::endl;
416 
417 
418  double ddd = cubicCoeff(1)*cubicCoeff(1)-4*cubicCoeff(0)*cubicCoeff(2); // b**2 - 4ac
419  double slam1, slam2 = 0;
420  // slam1 should be minimum and slam2 the maximum
421  if (cubicCoeff(0) > 0) {
422  slam1 = (- cubicCoeff(1) - std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
423  slam2 = (- cubicCoeff(1) + std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
424  }
425  else if (cubicCoeff(0) < 0) {
426  slam1 = (- cubicCoeff(1) + std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
427  slam2 = (- cubicCoeff(1) - std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
428  }
429  else { // case == 0 (-b/c)
430  slam1 = - gdel/g2del;
431  slam2 = slam1;
432  }
433 
434  std::cout << "slam1 & slam 2 " << slam1 << " " << slam2 << std::endl;
435 
436  // this should be the minimum otherwise inversion failed and I should do something else
437 
438  if (slam2 < undral) slam2 = undral;
439  if (slam2 > overal) slam2 = overal;
440 
441 
442  // I am stack somewhere - take a large step
443  if (std::fabs(slam2) < toler) slam2 = ( slam2 >=0 ) ? slamax : -slamax;
444 
445 
446  std::cout << "try with slam 2 " << slam2 << std::endl;
447 
448  double f2 = fcn(st.Vec()+slam2*step);
449 
450  std::cout << "try with slam 2 " << slam2 << " f2 " << f2 << std::endl;
451 
452  double fp;
453  // use this as new minimum
454  //bool noImpr = false;
455  if (f2 < fvmin) {
456  slam = slam2;
457  xvmin = slam;
458  fvmin = f2;
459  fp = fvmin;
460  } else {
461  // try with slam2 if it is better
462 
463  if (slam1 < undral) slam1 = undral;
464  if (slam1 > overal) slam1 = overal;
465 
466 
467  if (std::fabs(slam1) < toler) slam1 = ( slam1 >=0 ) ? -slamax : slamax;
468 
469  double f3 = fcn(st.Vec()+slam1*step);
470 
471  std::cout << "try with slam 1 " << slam1 << " f3 " << f3 << std::endl;
472 
473  if (f3 < fvmin) {
474  slam = slam1;
475  fp = fvmin;
476  xvmin = slam;
477  fvmin = f3;
478  } else {
479  // case both f2 and f3 did not produce a better result
480  if (f2 < f3) {
481  slam = slam1;
482  fp = f2;
483  } else {
484  slam = slam2;
485  fp = f3;
486  }
487  }
488  }
489 
490  bool iterate2 = false;
491  int niter = 0;
492 
493  int maxiter = 10;
494 
495  do {
496  iterate2 = false;
497 
498  std::cout << "\n iter" << niter << " test approx deriv ad second deriv at " << slam << " fp " << fp << std::endl;
499 
500  // estimate grad and second derivative at new point taking a step of 10-3
501  double h = 0.001*slam;
502  double fh = fcn(st.Vec()+(slam+h)*step);
503  double fl = fcn(st.Vec()+(slam-h)*step);
504  double df = (fh-fl)/(2.*h);
505  double df2 = (fh+fl-2.*fp )/(h*h);
506 
507  std::cout << "deriv: " << df << " , " << df2 << std::endl;
508 
509  // if I am in a point of still negative derivative
510  if ( std::fabs(df ) < prec.Eps() && std::fabs( df2 ) < prec.Eps() ) {
511  // try in opposite direction with an opposite value
512  slam = ( slam >=0 ) ? -slamax : slamax;
513  slamax *= 10;
514  fp = fcn(st.Vec()+slam*step);
515  }
516  else if ( std::fabs(df2 ) <= 0) { // gradient is significative different than zero then redo a cubic interpolation
517  // from new point
518 
519  return MnParabolaPoint(slam, fp); // should redo a cubic interpol. ??
520 // niter ++;
521 // if (niter > maxiter) break;
522 
523 // MinimumParameters pa = MinimumParameters(st.Vec() + stepNew, fp);
524 // gdel = stepNew(i)
525 // MnParabolaPoint pp = CubicSearch(fcn, st, stepNew, df, df2
526 
527 
528  }
529 
530  else
531  return MnParabolaPoint(slam, fp);
532 
533  niter++;
534  } while (niter < maxiter);
535 
536  return MnParabolaPoint(xvmin, fvmin);
537 
538 
539 }
540 
541 
542 // class describing Fcn function in one dimension
543  class ProjectedFcn : public ROOT::Math::IGenFunction {
544 
545  public:
546 
547  ProjectedFcn(const MnFcn & fcn, const MinimumParameters & pa, const MnAlgebraicVector & step) :
548  fFcn(fcn),
549  fPar(pa),
550  fStep(step)
551  {}
552 
553 
554  ROOT::Math::IGenFunction * Clone() const { return new ProjectedFcn(*this); }
555 
556  private:
557 
558  double DoEval(double x) const {
559  return fFcn(fPar.Vec() + x*fStep);
560  }
561 
562  const MnFcn & fFcn;
563  const MinimumParameters & fPar;
564  const MnAlgebraicVector & fStep;
565  };
566 
567 
568 MnParabolaPoint MnLineSearch::BrentSearch(const MnFcn& fcn, const MinimumParameters& st, const MnAlgebraicVector& step, double gdel, double g2del, const MnMachinePrecision& prec, bool debug) const {
569 
570  if (debug) {
571  std::cout<<"gdel= "<<gdel<<std::endl;
572  std::cout<<"g2del= "<<g2del<<std::endl;
573  for (unsigned int i = 0; i < step.size(); ++i) {
574  if (step(i) != 0) {
575  std::cout<<"step(i)= "<<step(i)<<std::endl;
576  std::cout<<"par(i) " <<st.Vec()(i) <<std::endl;
577  break;
578  }
579  }
580  }
581 
582  ProjectedFcn func(fcn, st, step);
583 
584 
585 // do first a cubic interpolation
586 
587  double f0 = st.Fval();
588  double f1 = fcn(st.Vec()+step);
589  f0 = func(0.0);
590  f1 = func(1.);
591  double fvmin = st.Fval();
592  double xvmin = 0.;
593  if (debug) std::cout << "f0 " << f0 << " f1 " << f1 << std::endl;
594 
595  if(f1 < f0) {
596  fvmin = f1;
597  xvmin = 1.;
598  }
599 // double toler8 = toler;
600 // double slamax = slambg;
601 // double flast = f1;
602  double slam = 1.;
603 
604  double undral = -1000;
605  double overal = 1000;
606 
607  double x0 = 0;
608 
609 // MnParabolaPoint p0(0., f0);
610 // MnParabolaPoint p1(slam, flast);
611 #ifdef USE_CUBIC
612 
613  ROOT::Math::SMatrix<double, 3> cubicMatrix;
614  ROOT::Math::SVector<double, 3> cubicCoeff; // cubic coefficients to be found
615  ROOT::Math::SVector<double, 3> bVec; // cubic coefficients to be found
616 
617  // cubic interpolation using the two points p0,p1 and the slope at p0 and the second derivative at p0
618 
619  // cut toler8 as function goes up
620  double x1 = slam;
621  cubicMatrix(0,0) = (x0*x0*x0 - x1*x1*x1)/3.;
622  cubicMatrix(0,1) = (x0*x0 - x1*x1)/2.;
623  cubicMatrix(0,2) = (x0 - x1);
624  cubicMatrix(1,0) = x0*x0;
625  cubicMatrix(1,1) = x0;
626  cubicMatrix(1,2) = 1;
627  cubicMatrix(2,0) = 2.*x0;
628  cubicMatrix(2,1) = 1;
629  cubicMatrix(2,2) = 0;
630 
631  bVec(0) = f0-f1;
632  bVec(1) = gdel;
633  bVec(2) = g2del;
634 
635  //if (debug) std::cout << "Matrix:\n " << cubicMatrix << std::endl;
636  if (debug) std::cout << "Vec:\n " << bVec << std::endl;
637 
638  // find the minimum need to invert the matrix
639  if (!cubicMatrix.Invert() ) {
640  std::cout << "Inversion failed - return " << std::endl;
641  return MnParabolaPoint(xvmin, fvmin);
642  }
643 
644  cubicCoeff = cubicMatrix * bVec;
645  if (debug) std::cout << "Cubic:\n " << cubicCoeff << std::endl;
646 
647 
648  double ddd = cubicCoeff(1)*cubicCoeff(1)-4*cubicCoeff(0)*cubicCoeff(2); // b**2 - 4ac
649  double slam1, slam2 = 0;
650  // slam1 should be minimum and slam2 the maximum
651  if (cubicCoeff(0) > 0) {
652  slam1 = (- cubicCoeff(1) - std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
653  slam2 = (- cubicCoeff(1) + std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
654  }
655  else if (cubicCoeff(0) < 0) {
656  slam1 = (- cubicCoeff(1) + std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
657  slam2 = (- cubicCoeff(1) - std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
658  }
659  else { // case == 0 (-b/c)
660  slam1 = - gdel/g2del;
661  slam2 = slam1;
662  }
663 
664  if (slam1 < undral) slam1 = undral;
665  if (slam1 > overal) slam1 = overal;
666 
667  if (slam2 < undral) slam2 = undral;
668  if (slam2 > overal) slam2 = overal;
669 
670 
671  double fs1 = func(slam1);
672  double fs2 = func(slam2);
673  std::cout << "slam1 & slam 2 " << slam1 << " " << slam2 << std::endl;
674  std::cout << "f(slam1) & f(slam 2) " << fs1 << " " << fs2 << std::endl;
675 
676  if (fs1 < fs2) {
677  x0 = slam1;
678  f0 = fs1;
679  }
680  else {
681  x0 = slam2;
682  f0 = fs2;
683  }
684 
685 #else
686  x0 = xvmin;
687  f0 = fvmin;
688 #endif
689 
690 
691  double astart = 100;
692 
693 
694  // do a Brent search in a large interval
695  double a = x0 -astart;
696  double b = x0 + astart;
697  //double x0 = 1;
698  int maxiter = 20;
699  double absTol = 0.5;
700  double relTol = 0.05;
701 
702  ROOT::Math::Minim1D::Type minType = ROOT::Math::Minim1D::BRENT;
703 
704  bool iterate = false;
705  int iter = 0;
706 
707  std::cout << "f(0)" << func(0.) << std::endl;
708  std::cout << "f(1)" << func(1.0) << std::endl;
709  std::cout << "f(3)" << func(3.0) << std::endl;
710  std::cout << "f(5)" << func(5.0) << std::endl;
711 
712  double fa = func(a);
713  double fb = func(b);
714  //double f0 = func(x0);
715  double toler = 0.01;
716  double delta0 = 5;
717  double delta = delta0;
718  bool enlarge = true;
719  bool scanmin = false;
720  double x2, f2 = 0;
721  double dir = 1;
722  int nreset = 0;
723 
724  do {
725 
726  std::cout << " iter " << iter << " a , b , x0 " << a << " " << b << " x0 " << x0 << std::endl;
727  std::cout << " fa , fb , f0 " << fa << " " << fb << " f0 " << f0 << std::endl;
728  if (fa <= f0 || fb <= f0) {
729  scanmin = false;
730  if (fa < fb) {
731  if (fa < fvmin) {
732  fvmin = fa;
733  xvmin = a;
734  }
735 // double f2 = fa;
736 // double x2 = a;
737  if (enlarge) {
738  x2 = a - 1000; // move lower
739  f2 = func(x2);
740  }
741  if ( std::fabs ( (fa -f2 )/(a-x2) ) > toler ) { // significant change in f continue to enlarge interval
742  x0 = a;
743  f0 = fa;
744  a = x2;
745  fa = f2;
746  enlarge = true;
747  } else {
748  // move direction of [a
749  // values change a little, start from central point try with x0 = x0 - delta
750  if (nreset == 0) dir = -1;
751  enlarge = false;
752  x0 = x0 + dir * delta;
753  f0 = func(x0);
754  // if reached limit try opposite direction , direction b]
755  if ( std::fabs (( f0 -fa )/(x0 -a) ) < toler ) {
756  delta = 10 * delta0/(10.*(nreset+1)); // decrease the delta if still not change observed
757  a = x0;
758  fa = f0;
759  x0 = delta;
760  f0 = func(x0);
761  dir *= -1;
762  nreset++;
763  std::cout << " A: Done a reset - scan in opposite direction ! " << std::endl;
764  }
765  delta *= 2; // double delta at next iteration
766  if (x0 < b && x0 > a )
767  scanmin = true;
768  else {
769  x0 = 0;
770  f0 = st.Fval();
771  }
772  }
773  }
774  else { // fb < fa
775  if (fb < fvmin) {
776  fvmin = fb;
777  xvmin = b;
778  }
779  if (enlarge) {
780  x2 = b + 1000; // move upper
781  f2 = func(x2);
782  }
783  if ( std::fabs ( (fb -f2 )/(x2-b) ) > toler ) { // significant change in f continue to enlarge interval
784  f0 = fb;
785  x0 = b;
786  b = x2; //
787  fb = f2;
788  enlarge = true;
789  } else {
790  // here move direction b
791  // values change a little, try with x0 = fa + delta
792  if (nreset == 0) dir = 1;
793  enlarge = false;
794  x0 = x0 + dir * delta;
795  f0 = func(x0);
796  // if reached limit try other side
797  if ( std::fabs ( ( f0 -fb )/(x0-b) ) < toler ) {
798  delta = 10 * delta0/(10.*(nreset+1)); // decrease the delta if still not change observed
799  b = x0;
800  fb = f0;
801  x0 = -delta;
802  f0 = func(x0);
803  dir *= -1;
804  nreset++;
805  std::cout << " B: Done a reset - scan in opposite direction ! " << std::endl;
806  }
807  delta *= 2; // double delta at next iteration
808  if (x0 > a && x0 < b)
809  scanmin = true;
810  else {
811  x0 = 0;
812  f0 = st.Fval();
813  }
814 
815  }
816 
817  }
818 
819  if ( f0 < fvmin) {
820  x0 = xvmin;
821  fvmin = f0;
822  }
823 
824  std::cout << " new x0 " << x0 << " " << f0 << std::endl;
825 
826  // use golden section
827  iterate = scanmin || enlarge;
828 
829  } else { // x0 is a min of [a,b]
830  iterate = false;
831  }
832 
833  iter++;
834  } while (iterate && iter < maxiter );
835 
836  if ( f0 > fa || f0 > fb)
837  // skip minim 1d try Minuit LS
838  // return (*this) (fcn, st, step, gdel, prec, debug);
839  return MnParabolaPoint(xvmin, fvmin);
840 
841 
842 
843 
844  std::cout << "1D minimization using " << minType << std::endl;
845 
846  ROOT::Math::Minimizer1D min(minType);
847 
848  min.SetFunction(func,x0, a, b);
849  int ret = min.Minimize(maxiter, absTol, relTol);
850 
851  std::cout << "result of GSL 1D minimization: " << ret << " niter " << min.Iterations() << std::endl;
852  std::cout << " xmin = " << min.XMinimum() << " fmin " << min.FValMinimum() << std::endl;
853 
854  return MnParabolaPoint(min.XMinimum(), min.FValMinimum());
855 
856 }
857 
858 #endif
859 
860 
861 
862 } // namespace Minuit2
863 
864 } // namespace ROOT