Logo 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 
23 namespace ROOT {
24 namespace Math {
25 
26 namespace 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) {
38  xmin = std::log(xmin);
39  xmax = std::log(xmax);
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
c
#define c(i)
Definition: RSha256.hxx:119
m
auto * m
Definition: textangle.C:8
Util.h
e
#define e(i)
Definition: RSha256.hxx:121
IFunction.h
r
ROOT::R::TRInterface & r
Definition: Object.C:4
xmax
float xmax
Definition: THbookFile.cxx:95
exp
double exp(double)
log
double log(double)
IFunctionfwd.h
x
Double_t x[n]
Definition: legend1.C:17
v
@ v
Definition: rootcling_impl.cxx:3635
b
#define b(i)
Definition: RSha256.hxx:118
q
float * q
Definition: THbookFile.cxx:89
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Math::BrentMethods::MinimStep
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...
Definition: BrentMethods.cxx:43
xmin
float xmin
Definition: THbookFile.cxx:95
Error.h
a
auto * a
Definition: textangle.C:12
ROOT::Math::Util::ToString
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition: Util.h:50
y
Double_t y[n]
Definition: legend1.C:17
MATH_INFO_MSG
#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:77
BrentMethods.h
ROOT::Math::IGenFunction
IBaseFunctionOneDim IGenFunction
Definition: IFunctionfwd.h:47
d
#define d(i)
Definition: RSha256.hxx:120
type
int type
Definition: TGX11.cxx:121
ROOT
VSD Structures.
Definition: StringConv.hxx:21
Math
ROOT::Math::BrentMethods::MinimBrent
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...
Definition: BrentMethods.cxx:145