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