ROOT   Reference Guide
BrentMethods.cxx
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Authors: David Gonzalez Maline 01/2008
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 , LCG ROOT MathLib Team *
7 * *
8 * *
9 **********************************************************************/
10
11#include "Math/BrentMethods.h"
12#include "Math/IFunction.h"
13#include "Math/IFunctionfwd.h"
14
15#include <cmath>
16#include <algorithm>
17
18#include "Math/Util.h"
19#include "Math/Error.h"
20
21#include <iostream>
22
23namespace ROOT {
24namespace Math {
25
26namespace BrentMethods {
27
28 double MinimStep(const IGenFunction* function, int type, double &xmin, double &xmax, double fy, int npx,bool logStep)
29{
30 // Grid search implementation, used to bracket the minimum and later
31 // use Brent's method with the bracketed interval
32 // The step of the search is set to (xmax-xmin)/npx
33 // type: 0,1 returns MinimumX
34 // 2,3 returns MaximumX
35 // 4-returns best X corresponding to fy
36
37 if (logStep) {
40 }
41
42
43 if (npx < 2) return 0.5*(xmax-xmin); // no bracketing - return just mid-point
44 double dx = (xmax-xmin)/(npx-1);
45 double xxmin = (logStep) ? std::exp(xmin) : xmin;
46 double yymin;
47 double xmiddle = 0;
48 // found an interval containing root (for type=4).
49 // For type !=4 an interval is always found
50 bool foundInterval = (type != 4);
51 if (type < 2)
52 yymin = (*function)(xxmin);
53 else if (type < 4)
54 yymin = -(*function)(xxmin);
55 else {
56 yymin = (*function)(xxmin)-fy;
57
58 // case root is at the interval boundaries
59 if (yymin==0) {
60 xmin = xxmin;
61 xmax = xxmin;
62 return xxmin;
63 }
64 }
65 for (int i=1; i<=npx-1; i++) {
66 double x = xmin + i*dx;
67 if (logStep) x = std::exp(x);
68 double y = 0;
69 if (type < 2)
70 y = (*function)(x);
71 else if (type < 4)
72 y = -(*function)(x);
73 else
74 y = (*function)(x)-fy;
75
76 if ( type < 4 && y < yymin) {
77 xxmin = x;
78 yymin = y;
79 }
80 // when looking for root break at first instance
81 if (type == 4 ) {
82 // if root is at interval boundaries
83 if (y == 0) {
84 xmin = x;
85 xmax = x;
86 xmiddle = x;
87 foundInterval = true;
88 break;
89 }
90 // found good interval if sign product is negative or equal zero to
91 if (std::copysign(1.,y)*std::copysign(1.,yymin) < 0 ) {
92 xmin = xxmin; // previous value
93 xmax = x; // current value
94 xmiddle = 0.5*(xmax+xmin);
95 foundInterval = true;
96 break;
97 }
98 // continue bracketing
99 xxmin = x;
100 yymin = y;
101 }
102 }
103
104 if (type < 4 ) {
105 if (logStep) {
106 xmin = std::exp(xmin);
107 xmax = std::exp(xmax);
108 }
109 xmin = std::max(xmin,xxmin-dx);
110 xmax = std::min(xmax,xxmin+dx);
111 xmiddle = std::min(xxmin,xmax);
112 }
113
114 //std::cout << "bracketing result " << xmin << " " << xmax << "middle value " << xmiddle << std::endl;
115 if (!foundInterval) {
116 MATH_INFO_MSG("BrentMethods::MinimStep", "Grid search failed to find a root in the interval ");
117 std::string msg = "xmin = ";
119 msg += std::string(" xmax = ");
121 msg += std::string(" npts = ");
122 msg += ROOT::Math::Util::ToString(npx);
123 MATH_INFO_MSG("BrentMethods::MinimStep", msg.c_str());
124 xmin = 1;
125 xmax = 0;
126 }
127
128 return xmiddle;
129}
130 double MinimBrent(const IGenFunction* function, int type, double &xmin, double &xmax, double xmiddle, double fy, bool &ok, int &niter, double epsabs, double epsrel, int itermax)
131{
132 //Finds a minimum of a function, if the function is unimodal between xmin and xmax
133 //This method uses a combination of golden section search and parabolic interpolation
134 //Details about convergence and properties of this algorithm can be
135 //found in the book by R.P.Brent "Algorithms for Minimization Without Derivatives"
136 //or in the "Numerical Recipes", chapter 10.2
137 // convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs)
138 //
139 //type: 0-returns MinimumX
140 // 1-returns Minimum
141 // 2-returns MaximumX
142 // 3-returns Maximum
143 // 4-returns X corresponding to fy
144 //if ok=true the method has converged
145
146 const double c = 3.81966011250105097e-01; // (3.-std::sqrt(5.))/2. (comes from golden section)
147 double u, v, w, x, fv, fu, fw, fx, e, p, q, r, t2, d=0, m, tol;
148 v = w = x = xmiddle;
149 e=0;
150
151 double a=xmin;
152 double b=xmax;
153 if (type < 2)
154 fv = fw = fx = (*function)(x);
155 else if (type < 4)
156 fv = fw = fx = -(*function)(x);
157 else
158 fv = fw = fx = std::fabs((*function)(x)-fy);
159
160 for (int i=0; i<itermax; i++){
161 m=0.5*(a + b);
162 tol = epsrel*(std::fabs(x))+epsabs;
163 t2 = 2*tol;
164
165 if (std::fabs(x-m) <= (t2-0.5*(b-a))) {
166 //converged, return x
167 ok=true;
168 niter = i-1;
169 if (type==1)
170 return fx;
171 else if (type==3)
172 return -fx;
173 else
174 return x;
175 }
176
177 if (std::fabs(e)>tol){
178 //fit parabola
179 r = (x-w)*(fx-fv);
180 q = (x-v)*(fx-fw);
181 p = (x-v)*q - (x-w)*r;
182 q = 2*(q-r);
183 if (q>0) p=-p;
184 else q=-q;
185 r=e; // Deltax before last
186 e=d; // last delta x
187 // current Deltax = p/q
188 // take a parabolic step only if:
189 // Deltax < 0.5* (DeltaX before last) && Deltax > a && Deltax < b
190 // (a BUG in testing this condition is fixed 11/3/2010 (with revision 32544)
191 if (std::fabs(p) >= std::fabs(0.5*q*r) || p <= q*(a-x) || p >= q*(b-x)) {
192 // condition fails - do not take parabolic step
193 e=(x>=m ? a-x : b-x);
194 d = c*e;
195 } else {
196 // take a parabolic interpolation step
197 d = p/q;
198 u = x+d;
199 if (u-a < t2 || b-u < t2)
200 //d=TMath::Sign(tol, m-x);
201 d=(m-x >= 0) ? std::fabs(tol) : -std::fabs(tol);
202 }
203 } else {
204 e=(x>=m ? a-x : b-x);
205 d = c*e;
206 }
207 u = (std::fabs(d)>=tol ? x+d : x+ ((d >= 0) ? std::fabs(tol) : -std::fabs(tol)) );
208 if (type < 2)
209 fu = (*function)(u);
210 else if (type < 4)
211 fu = -(*function)(u);
212 else
213 fu = std::fabs((*function)(u)-fy);
214 //update a, b, v, w and x
215 if (fu<=fx){
216 if (u<x) b=x;
217 else a=x;
218 v=w; fv=fw; w=x; fw=fx; x=u; fx=fu;
219 } else {
220 if (u<x) a=u;
221 else b=u;
222 if (fu<=fw || w==x){
223 v=w; fv=fw; w=u; fw=fu;
224 }
225 else if (fu<=fv || v==x || v==w){
226 v=u; fv=fu;
227 }
228 }
229 }
230 //didn't converge
231 ok = false;
232 xmin = a;
233 xmax = b;
234 niter = itermax;
235 return x;
236
237}
238
239} // end namespace BrentMethods
240} // end namespace Math
241} // ned namespace ROOT
#define MATH_INFO_MSG(loc, str)
Pre-processor macro to report messages which can be configured to use ROOT error or simply an std::io...
Definition: Error.h:76
ROOT::R::TRInterface & r
Definition: Object.C:4
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
int type
Definition: TGX11.cxx:120
float xmin
Definition: THbookFile.cxx:93
float * q
Definition: THbookFile.cxx:87
float xmax
Definition: THbookFile.cxx:93
double exp(double)
double log(double)
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
double MinimStep(const IGenFunction *f, int type, double &xmin, double &xmax, double fy, int npx=100, bool useLog=false)
Grid search implementation, used to bracket the minimum and later use Brent's method with the bracket...
double MinimBrent(const IGenFunction *f, int type, double &xmin, double &xmax, double xmiddle, double fy, bool &ok, int &niter, double epsabs=1.E-8, double epsrel=1.E-10, int maxiter=100)
Finds a minimum of a function, if the function is unimodal between xmin and xmax This method uses a c...
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition: Util.h:50
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:151
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: StringConv.hxx:21
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12