ROOT  6.06/09
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 
14 #include <cmath>
15 #include <algorithm>
16 
17 #ifndef ROOT_Math_Error
18 #include "Math/Error.h"
19 #endif
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-returns MinimumX
34  // 1-returns Minimum
35  // 2-returns MaximumX
36  // 3-returns Maximum
37  // 4-returns X corresponding to fy
38 
39  if (logStep) {
40  xmin = std::log(xmin);
41  xmax = std::log(xmax);
42  }
43 
44 
45  if (npx < 2) return 0.5*(xmax-xmin); // no bracketing - return just mid-point
46  double dx = (xmax-xmin)/(npx-1);
47  double xxmin = (logStep) ? std::exp(xmin) : xmin;
48  double yymin;
49  if (type < 2)
50  yymin = (*function)(xxmin);
51  else if (type < 4)
52  yymin = -(*function)(xxmin);
53  else
54  yymin = std::fabs((*function)(xxmin)-fy);
55 
56  for (int i=1; i<=npx-1; i++) {
57  double x = xmin + i*dx;
58  if (logStep) x = std::exp(x);
59  double y = 0;
60  if (type < 2)
61  y = (*function)(x);
62  else if (type < 4)
63  y = -(*function)(x);
64  else
65  y = std::fabs((*function)(x)-fy);
66  if (y < yymin) {
67  xxmin = x;
68  yymin = y;
69  }
70  }
71 
72  if (logStep) {
73  xmin = std::exp(xmin);
74  xmax = std::exp(xmax);
75  }
76 
77 
78  xmin = std::max(xmin,xxmin-dx);
79  xmax = std::min(xmax,xxmin+dx);
80 
81  return std::min(xxmin, xmax);
82 }
83 
84  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)
85 {
86  //Finds a minimum of a function, if the function is unimodal between xmin and xmax
87  //This method uses a combination of golden section search and parabolic interpolation
88  //Details about convergence and properties of this algorithm can be
89  //found in the book by R.P.Brent "Algorithms for Minimization Without Derivatives"
90  //or in the "Numerical Recipes", chapter 10.2
91  // convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs)
92  //
93  //type: 0-returns MinimumX
94  // 1-returns Minimum
95  // 2-returns MaximumX
96  // 3-returns Maximum
97  // 4-returns X corresponding to fy
98  //if ok=true the method has converged
99 
100  const double c = 3.81966011250105097e-01; // (3.-std::sqrt(5.))/2. (comes from golden section)
101  double u, v, w, x, fv, fu, fw, fx, e, p, q, r, t2, d=0, m, tol;
102  v = w = x = xmiddle;
103  e=0;
104 
105  double a=xmin;
106  double b=xmax;
107  if (type < 2)
108  fv = fw = fx = (*function)(x);
109  else if (type < 4)
110  fv = fw = fx = -(*function)(x);
111  else
112  fv = fw = fx = std::fabs((*function)(x)-fy);
113 
114  for (int i=0; i<itermax; i++){
115  m=0.5*(a + b);
116  tol = epsrel*(std::fabs(x))+epsabs;
117  t2 = 2*tol;
118 
119  if (std::fabs(x-m) <= (t2-0.5*(b-a))) {
120  //converged, return x
121  ok=true;
122  niter = i-1;
123  if (type==1)
124  return fx;
125  else if (type==3)
126  return -fx;
127  else
128  return x;
129  }
130 
131  if (std::fabs(e)>tol){
132  //fit parabola
133  r = (x-w)*(fx-fv);
134  q = (x-v)*(fx-fw);
135  p = (x-v)*q - (x-w)*r;
136  q = 2*(q-r);
137  if (q>0) p=-p;
138  else q=-q;
139  r=e; // Deltax before last
140  e=d; // last delta x
141  // current Deltax = p/q
142  // take a parabolic step only if:
143  // Deltax < 0.5* (DeltaX before last) && Deltax > a && Deltax < b
144  // (a BUG in testing this condition is fixed 11/3/2010 (with revision 32544)
145  if (std::fabs(p) >= std::fabs(0.5*q*r) || p <= q*(a-x) || p >= q*(b-x)) {
146  // condition fails - do not take parabolic step
147  e=(x>=m ? a-x : b-x);
148  d = c*e;
149  } else {
150  // take a parabolic interpolation step
151  d = p/q;
152  u = x+d;
153  if (u-a < t2 || b-u < t2)
154  //d=TMath::Sign(tol, m-x);
155  d=(m-x >= 0) ? std::fabs(tol) : -std::fabs(tol);
156  }
157  } else {
158  e=(x>=m ? a-x : b-x);
159  d = c*e;
160  }
161  u = (std::fabs(d)>=tol ? x+d : x+ ((d >= 0) ? std::fabs(tol) : -std::fabs(tol)) );
162  if (type < 2)
163  fu = (*function)(u);
164  else if (type < 4)
165  fu = -(*function)(u);
166  else
167  fu = std::fabs((*function)(u)-fy);
168  //update a, b, v, w and x
169  if (fu<=fx){
170  if (u<x) b=x;
171  else a=x;
172  v=w; fv=fw; w=x; fw=fx; x=u; fx=fu;
173  } else {
174  if (u<x) a=u;
175  else b=u;
176  if (fu<=fw || w==x){
177  v=w; fv=fw; w=u; fw=fu;
178  }
179  else if (fu<=fv || v==x || v==w){
180  v=u; fv=fu;
181  }
182  }
183  }
184  //didn't converge
185  ok = false;
186  xmin = a;
187  xmax = b;
188  niter = itermax;
189  return x;
190 
191 }
192 
193 } // end namespace BrentMethods
194 } // end namespace Math
195 } // ned namespace ROOT
float xmin
Definition: THbookFile.cxx:93
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
Definition: vector.h:433
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:133
static float fu(float x)
Definition: main.cpp:53
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
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...
TArc * a
Definition: textangle.C:12
Double_t x[n]
Definition: legend1.C:17
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
const double tol
ROOT::R::TRInterface & r
Definition: Object.C:4
SVector< double, 2 > v
Definition: Dict.h:5
TMarker * m
Definition: textangle.C:8
float xmax
Definition: THbookFile.cxx:93
int type
Definition: TGX11.cxx:120
Double_t y[n]
Definition: legend1.C:17
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
Namespace for new Math classes and functions.
double exp(double)
float * q
Definition: THbookFile.cxx:87
double log(double)