// @(#)root/mathcore:\$Id: BrentMethods.cxx 36905 2010-11-24 15:44:34Z moneta \$ // Authors: David Gonzalez Maline 01/2008 (MinimStep & MinimBrent) // Amnon Harel 04/2012 (BracketStep, GridBracketStep & BrentMethod) /********************************************************************** * * * Copyright (c) 2006 , LCG ROOT MathLib Team * * * * * **********************************************************************/ #include "ModBrentMethods.h" #include "Math/IFunction.h" #include #include #ifndef ROOT_Math_Error #include "Math/Error.h" #endif #include using std::cout; using std::cerr; using std::endl; namespace ModBrentMethods { double GridBracketStep(const ROOT::Math::IGenFunction* function, double &xmin, double &xmax, double fy, int npx,bool logStep) { // Grid search implementation, used to bracket the solution f(x)=fy // and later use Brent's method with the bracketed interval // The step of the search is set to (xmax-xmin)/npx // If no sign changes (of g(x) = f(x)-fy) are found, optimistically returns // double-step sized bracket around the point of closest approach (a-la MinimStep). // Caller must ensure xmax > xmin cout<<"GBstep [ "<= zmin) { zmin = z_opp - dz; zmid = z_opp - 0.5*dz; zmax = z_opp; } else { zmin = std::max(zmin,z_closest-1.021*dz); // a bit different than dz so we'll get more new grid points zmid = z_closest; zmax = std::min(zmax,z_closest+1.031*dz); // a bit different than dz so we'll get more new grid points } xmin = logStep ? std::exp(zmin) : zmin; xmax = logStep ? std::exp(zmax) : zmax; double xmid = logStep ? std::exp(zmid) : zmid; cout<<"GBstep gave "< itermax) { return -2; } } return 0; } double BrentMethod(const ROOT::Math::IGenFunction* function, double &xmin, double &xmax, double fy, bool &ok, int &niter, double epsabs, double epsrel, int itermax) { //Finds x in bracket [xmin, xmax] such that function(x)=fy and x. //The input xmin and xmax must define a valid bracket, at one func>=fy and at the other func<=fy. //x is found using Brent's method, as a root of g(x)=function(x)-fy. //Brent's method uses either bisection, linear interpolation, or inverse quadratic //interpolation in each iteration. //Details about convergence and properties of this algorithm can be //found in the book by R.P.Brent "Algorithms for Minimization Without Derivatives" //or in the "Numerical Recipes", chapter 10.2 // //The output root position should be accurate to 4*epsrel*abs(x)+epsabs //Brent had esprel as the relative machine precision defined as the smallest //representable number such that 1.+macheps .gt. 1. // //if ok=true the method has converged // Implementation based on Brent's Fortran code (http://www.netlib.org/go/zeroin.f) // Hence this is more or less assembly code with comments // Here f(x) = function(x) - fy double a=xmin, b=xmax; double fa = (*function)(xmin)-fy; double fb = (*function)(xmax)-fy; if( fa * fb > 0 ) { ok = false; niter = -1; return xmin - std::fabs( xmin ) - 1.; }; // mark illegal inputs double c = a; double fc = fa; double step = b-a; // proposed step size double old_step = step; // An old step size to be used in comparisons. // After a bisection, this is the previous step. // After an interpolation, the next-to-previous step. double tol; // the tolerance for the current b double xm; // bracket half width, signed double s,p,r,q; // temporaries. p and q may be meaningful as commented below. niter = 0; while( true ) { // at this point the bracket edges are b and c ++niter; // ensure b has smaller |f| than c. If not so originally, also set a=c (thus forgeting // about the previous guess, probably because in this case b isn't between a and c) if (std::fabs(fc) < std::fabs(fb)) { a = b; b = c; c = a; fa = fb; fb = fc; fc = fa; } // at this point b is the current best guess. a contains the previous estimate, which may be equal to c. tol = 2*epsrel*(std::fabs(b))+0.5*epsabs; xm = 0.5*(c - b); if (std::fabs(xm) <= tol || fb==0) { // convergence and succesful exit ok = true; return b; } if( niter > itermax ) { // didn't converge ok = false; xmin = std::min( b, c ); xmax = std::max( b, c ); return b; } // Choose Step size if( std::fabs(old_step) < tol || // the old step was suspicously small std::fabs(fa) <= std::fabs(fb) // haven't improved the guess in the previous step ) { // use bisection. old_step = step = xm; } else { // try to use an interpolation s = fb/fa; // no division by zero since |fa|=(?)=|fc|>|fb|>0 // the interpolations calculate a denominator, q, and an enumerator, p, for d if (a == c) { // use linear interpolation, since only two points are in play p = 2.*xm*s; q = 1.-s; } else { // use inverse quadratic interpolation q = fa/fc; // here used as a temp variable, like s and r r = fb/fc; p = s*(2.*xm*q*(q-r) - (b-a)*(r-1.)); q = (q-1.)*(r-1.)*(s-1.); } if( p > 0 ) { // ensure positive p and flip the sign of p/q q *= -1; } else { p *= -1; } // Note: slight simplification of zeroin.f, as "s" isn't really used if ( ( 2.*p < 3.*xm*q-std::fabs(tol*q)) && // ensures the |b-new_b| interval is reduced (from |b-c|) by atleast 25% + tolerance ( p < std::fabs(0.5*old_step*q)) // proposed interpolation step smaller than half of "the old step" ) { // going ahead with the interpolation old_step = step; step = p/q; } else { // giving up on the interpolation and using bisection instead old_step = step = xm; } } // Prepare the new interval: // 1) first update b ( and store the old guess in a ) a = b; fa = fb; if (std::fabs(step) > tol) { b = b + step; } else { if (xm > 0.) { b += tol; } else { b -= tol; } } fb = (*function)(b)-fy; // 2) If needed, update c, to ensure b and c have opposite signs in the next iteration if (fb * fc > 0.) { c = a; fc = fa; old_step = step = b-a; } } } double MinimStep(const ROOT::Math::IGenFunction* function, int type, double &xmin, double &xmax, double fy, int npx,bool logStep) { // Grid search implementation, used to bracket the minimum and later // use Brent's method with the bracketed interval // The step of the search is set to (xmax-xmin)/npx // type: 0-returns MinimumX // 1-returns Minimum // 2-returns MaximumX // 3-returns Maximum // 4-returns X corresponding to fy if (logStep) { xmin = std::log(xmin); xmax = std::log(xmax); } //cout<<"AH MinimStep type: "< [ "<tol){ //fit parabola r = (x-w)*(fx-fv); q = (x-v)*(fx-fw); p = (x-v)*q - (x-w)*r; q = 2*(q-r); if (q>0) p=-p; else q=-q; r=e; // Deltax before last e=d; // last delta x // current Deltax = p/q // take a parabolic step only if: // Deltax < 0.5* (DeltaX before last) && Deltax > a && Deltax < b // (a BUG in testing this condition is fixed 11/3/2010 (with revision 32544) if (std::fabs(p) >= std::fabs(0.5*q*r) || p <= q*(a-x) || p >= q*(b-x)) { // condition fails - do not take parabolic step e=(x>=m ? a-x : b-x); d = c*e; } else { // take a parabolic interpolation step d = p/q; u = x+d; if (u-a < t2 || b-u < t2) //d=TMath::Sign(tol, m-x); d=(m-x >= 0) ? std::fabs(tol) : -std::fabs(tol); } } else { e=(x>=m ? a-x : b-x); d = c*e; } u = (std::fabs(d)>=tol ? x+d : x+ ((d >= 0) ? std::fabs(tol) : -std::fabs(tol)) ); if (type < 2) fu = (*function)(u); else if (type < 4) fu = -(*function)(u); else fu = std::fabs((*function)(u)-fy); //update a, b, v, w and x if (fu<=fx){ if (u