29 #include "Math/Minimizer1D.h" 65 std::cout<<
"gdel= "<<gdel<<std::endl;
66 std::cout<<
"step= "<<step<<std::endl;
69 double overal = 1000.;
70 double undral = -100.;
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;
85 if(
fabs(slamin) < prec.
Eps()) slamin = prec.
Eps();
86 slamin *= prec.
Eps2();
88 double f0 = st.
Fval();
89 double f1 = fcn(st.
Vec()+step);
91 double fvmin = st.
Fval();
98 double toler8 = toler;
99 double slamax = slambg;
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;
126 double denom = 2.*(flast-f0-gdel*slam)/(slam*slam);
127 if (debug) std::cout<<
"denom= "<<denom<<std::endl;
134 if (debug) std::cout<<
"new slam= "<<slam<<std::endl;
137 #ifdef TRY_OPPOSIT_DIR 139 bool slamIsNeg =
false;
144 if (debug) std::cout<<
"slam is negative- set to " << slamax << std::endl;
145 #ifdef TRY_OPPOSITE_DIR 148 if (debug) std::cout<<
"slam is negative- compare values between "<< slamNeg <<
" and " << slamax << std::endl;
155 if (debug) std::cout <<
"slam larger than mac value - set to " << slamax << std::endl;
159 if (debug) std::cout <<
"slam too small - set to " << toler8 << std::endl;
164 if (debug) std::cout <<
"slam smaller than " << slamin <<
" return " << std::endl;
170 if(
fabs(slam - 1.) < toler8 && p1.
Y() < p0.
Y()) {
176 if(
fabs(slam - 1.) < toler8) slam = 1. + toler8;
183 f2 = fcn(st.
Vec() + slam*step);
187 #ifdef TRY_OPPOSITE_DIR 190 double f2alt = fcn(st.
Vec() + slamNeg*step);
208 overal = slam - toler8;
213 }
while(iterate && niter < maxiter);
214 if(niter >= maxiter) {
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;
229 slamax = std::max(slamax, alpha*
fabs(xvmin));
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;
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;
247 if(slam > xvmin + slamax) slam = xvmin + slamax;
248 if(slam < xvmin - slamax) slam = xvmin - slamax;
251 if(slam > overal) slam = overal;
253 if(slam < undral) slam = undral;
257 std::cout<<
" slam= "<<slam<<
" undral " << undral <<
" overal " << overal << std::endl;
263 if (debug) std::cout <<
" iterate on f3- slam " << niter <<
" slam = " << slam <<
" xvmin " << xvmin << std::endl;
266 double toler9 = std::max(toler8,
fabs(toler8*slam));
268 if(
fabs(p0.
X() - slam) < toler9 ||
269 fabs(p1.X() - slam) < toler9 ||
270 fabs(p2.
X() - slam) < toler9) {
280 f3 = fcn(st.
Vec() + slam*step);
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;
286 if(f3 > p0.
Y() && f3 > p1.Y() && f3 > p2.
Y()) {
288 std::cout<<
"f3 worse than all three previous"<<std::endl;
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);
294 std::cout<<
"new slam= "<<slam<<std::endl;
299 }
while(iterate && niter < maxiter);
300 if(niter >= maxiter) {
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;
310 if (debug) std::cout <<
" f3 " << f3 <<
" fvmin " << fvmin <<
" xvmin " << xvmin << std::endl;
315 if(slam > xvmin) overal = std::min(overal, slam-toler8);
316 if(slam < xvmin) undral = std::max(undral, slam+toler8);
320 }
while(niter < maxiter);
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;
342 std::cout<<
"gdel= "<<gdel<<std::endl;
343 std::cout<<
"g2del= "<<g2del<<std::endl;
344 std::cout<<
"step= "<<step<<std::endl;
348 double overal = 100.;
349 double undral = -100.;
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;
361 if(
fabs(slamin) < prec.
Eps()) slamin = prec.
Eps();
362 slamin *= prec.
Eps2();
364 double f0 = st.
Fval();
365 double f1 = fcn(st.
Vec()+step);
366 double fvmin = st.
Fval();
368 if (debug) std::cout <<
"f0 " << f0 <<
" f1 " << f1 << std::endl;
374 double toler8 = toler;
375 double slamax = slambg;
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;
406 if (debug) std::cout <<
"Vec:\n " << bVec << std::endl;
409 if (!cubicMatrix.
Invert() ) {
410 std::cout <<
"Inversion failed - return " << std::endl;
414 cubicCoeff = cubicMatrix * bVec;
415 if (debug) std::cout <<
"Cubic:\n " << cubicCoeff << std::endl;
418 double ddd = cubicCoeff(1)*cubicCoeff(1)-4*cubicCoeff(0)*cubicCoeff(2);
419 double slam1, slam2 = 0;
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) );
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) );
430 slam1 = - gdel/g2del;
434 std::cout <<
"slam1 & slam 2 " << slam1 <<
" " << slam2 << std::endl;
438 if (slam2 < undral) slam2 = undral;
439 if (slam2 > overal) slam2 = overal;
443 if (
std::fabs(slam2) < toler) slam2 = ( slam2 >=0 ) ? slamax : -slamax;
446 std::cout <<
"try with slam 2 " << slam2 << std::endl;
448 double f2 = fcn(st.
Vec()+slam2*step);
450 std::cout <<
"try with slam 2 " << slam2 <<
" f2 " << f2 << std::endl;
463 if (slam1 < undral) slam1 = undral;
464 if (slam1 > overal) slam1 = overal;
467 if (
std::fabs(slam1) < toler) slam1 = ( slam1 >=0 ) ? -slamax : slamax;
469 double f3 = fcn(st.
Vec()+slam1*step);
471 std::cout <<
"try with slam 1 " << slam1 <<
" f3 " << f3 << std::endl;
490 bool iterate2 =
false;
498 std::cout <<
"\n iter" << niter <<
" test approx deriv ad second deriv at " << slam <<
" fp " << fp << std::endl;
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);
507 std::cout <<
"deriv: " << df <<
" , " << df2 << std::endl;
512 slam = ( slam >=0 ) ? -slamax : slamax;
514 fp = fcn(st.
Vec()+slam*step);
534 }
while (niter < maxiter);
558 double DoEval(
double x)
const {
559 return fFcn(fPar.Vec() + x*fStep);
571 std::cout<<
"gdel= "<<gdel<<std::endl;
572 std::cout<<
"g2del= "<<g2del<<std::endl;
573 for (
unsigned int i = 0; i < step.
size(); ++i) {
575 std::cout<<
"step(i)= "<<step(i)<<std::endl;
576 std::cout<<
"par(i) " <<st.
Vec()(i) <<std::endl;
582 ProjectedFcn func(fcn, st, step);
587 double f0 = st.
Fval();
588 double f1 = fcn(st.
Vec()+step);
591 double fvmin = st.
Fval();
593 if (debug) std::cout <<
"f0 " << f0 <<
" f1 " << f1 << std::endl;
604 double undral = -1000;
605 double overal = 1000;
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;
636 if (debug) std::cout <<
"Vec:\n " << bVec << std::endl;
639 if (!cubicMatrix.
Invert() ) {
640 std::cout <<
"Inversion failed - return " << std::endl;
644 cubicCoeff = cubicMatrix * bVec;
645 if (debug) std::cout <<
"Cubic:\n " << cubicCoeff << std::endl;
648 double ddd = cubicCoeff(1)*cubicCoeff(1)-4*cubicCoeff(0)*cubicCoeff(2);
649 double slam1, slam2 = 0;
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) );
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) );
660 slam1 = - gdel/g2del;
664 if (slam1 < undral) slam1 = undral;
665 if (slam1 > overal) slam1 = overal;
667 if (slam2 < undral) slam2 = undral;
668 if (slam2 > overal) slam2 = overal;
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;
695 double a = x0 -astart;
696 double b = x0 + astart;
700 double relTol = 0.05;
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;
717 double delta = delta0;
719 bool scanmin =
false;
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) {
741 if (
std::fabs ( (fa -f2 )/(a-x2) ) > toler ) {
750 if (nreset == 0) dir = -1;
752 x0 = x0 + dir * delta;
755 if (
std::fabs (( f0 -fa )/(x0 -a) ) < toler ) {
756 delta = 10 * delta0/(10.*(nreset+1));
763 std::cout <<
" A: Done a reset - scan in opposite direction ! " << std::endl;
766 if (x0 < b && x0 > a )
783 if (
std::fabs ( (fb -f2 )/(x2-b) ) > toler ) {
792 if (nreset == 0) dir = 1;
794 x0 = x0 + dir * delta;
797 if (
std::fabs ( ( f0 -fb )/(x0-b) ) < toler ) {
798 delta = 10 * delta0/(10.*(nreset+1));
805 std::cout <<
" B: Done a reset - scan in opposite direction ! " << std::endl;
808 if (x0 > a && x0 < b)
824 std::cout <<
" new x0 " << x0 <<
" " << f0 << std::endl;
827 iterate = scanmin || enlarge;
834 }
while (iterate && iter < maxiter );
836 if ( f0 > fa || f0 > fb)
844 std::cout <<
"1D minimization using " << minType << std::endl;
846 ROOT::Math::Minimizer1D min(minType);
848 min.SetFunction(func,x0, a, b);
849 int ret = min.Minimize(maxiter, absTol, relTol);
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;
double Y() const
Accessor to the y (second) coordinate.
int iterate(rng_state_t *X)
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Namespace for new ROOT classes and functions.
const MnAlgebraicVector & Vec() const
double Min() const
Calculates the x coordinate of the Minimum of the parabola.
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
double A() const
Accessor to the coefficient of the quadratic term.
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
determines the relative floating point arithmetic precision.
SMatrix: a generic fixed size D1 x D2 Matrix class.
double C() const
Accessor to the coefficient of the constant term.
This class defines a parabola of the form a*x*x + b*x + c.
Wrapper class to FCNBase interface used internally by Minuit.
double B() const
Accessor to the coefficient of the linear term.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
double X() const
Accessor to the x (first) coordinate.
unsigned int size() const
* x
Deprecated and error prone model selection interface.
double Eps2() const
eps2 returns 2*sqrt(eps)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Type
Enumeration with One Dimensional Minimizer Algorithms.
MnParabolaPoint operator()(const MnFcn &, const MinimumParameters &, const MnAlgebraicVector &, double, const MnMachinePrecision &, bool debug=false) const
Perform a line search from position defined by the vector st along the direction step, where the length of vector step gives the expected position of Minimum.
SVector: a generic fixed size Vector class.