Logo ROOT   6.08/07
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 #ifndef ROOT_Math_Error
19 #include "Math/Error.h"
20 #endif
21 
22 #include <iostream>
23 
24 namespace ROOT {
25 namespace Math {
26 
27 namespace BrentMethods {
28 
29  double MinimStep(const IGenFunction* function, int type, double &xmin, double &xmax, double fy, int npx,bool logStep)
30 {
31  // Grid search implementation, used to bracket the minimum and later
32  // use Brent's method with the bracketed interval
33  // The step of the search is set to (xmax-xmin)/npx
34  // type: 0-returns MinimumX
35  // 1-returns Minimum
36  // 2-returns MaximumX
37  // 3-returns Maximum
38  // 4-returns X corresponding to fy
39 
40  if (logStep) {
41  xmin = std::log(xmin);
42  xmax = std::log(xmax);
43  }
44 
45 
46  if (npx < 2) return 0.5*(xmax-xmin); // no bracketing - return just mid-point
47  double dx = (xmax-xmin)/(npx-1);
48  double xxmin = (logStep) ? std::exp(xmin) : xmin;
49  double yymin;
50  if (type < 2)
51  yymin = (*function)(xxmin);
52  else if (type < 4)
53  yymin = -(*function)(xxmin);
54  else
55  yymin = std::fabs((*function)(xxmin)-fy);
56 
57  for (int i=1; i<=npx-1; i++) {
58  double x = xmin + i*dx;
59  if (logStep) x = std::exp(x);
60  double y = 0;
61  if (type < 2)
62  y = (*function)(x);
63  else if (type < 4)
64  y = -(*function)(x);
65  else
66  y = std::fabs((*function)(x)-fy);
67  if (y < yymin) {
68  xxmin = x;
69  yymin = y;
70  }
71  }
72 
73  if (logStep) {
74  xmin = std::exp(xmin);
75  xmax = std::exp(xmax);
76  }
77 
78 
79  xmin = std::max(xmin,xxmin-dx);
80  xmax = std::min(xmax,xxmin+dx);
81 
82  return std::min(xxmin, xmax);
83 }
84 
85  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)
86 {
87  //Finds a minimum of a function, if the function is unimodal between xmin and xmax
88  //This method uses a combination of golden section search and parabolic interpolation
89  //Details about convergence and properties of this algorithm can be
90  //found in the book by R.P.Brent "Algorithms for Minimization Without Derivatives"
91  //or in the "Numerical Recipes", chapter 10.2
92  // convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs)
93  //
94  //type: 0-returns MinimumX
95  // 1-returns Minimum
96  // 2-returns MaximumX
97  // 3-returns Maximum
98  // 4-returns X corresponding to fy
99  //if ok=true the method has converged
100 
101  const double c = 3.81966011250105097e-01; // (3.-std::sqrt(5.))/2. (comes from golden section)
102  double u, v, w, x, fv, fu, fw, fx, e, p, q, r, t2, d=0, m, tol;
103  v = w = x = xmiddle;
104  e=0;
105 
106  double a=xmin;
107  double b=xmax;
108  if (type < 2)
109  fv = fw = fx = (*function)(x);
110  else if (type < 4)
111  fv = fw = fx = -(*function)(x);
112  else
113  fv = fw = fx = std::fabs((*function)(x)-fy);
114 
115  for (int i=0; i<itermax; i++){
116  m=0.5*(a + b);
117  tol = epsrel*(std::fabs(x))+epsabs;
118  t2 = 2*tol;
119 
120  if (std::fabs(x-m) <= (t2-0.5*(b-a))) {
121  //converged, return x
122  ok=true;
123  niter = i-1;
124  if (type==1)
125  return fx;
126  else if (type==3)
127  return -fx;
128  else
129  return x;
130  }
131 
132  if (std::fabs(e)>tol){
133  //fit parabola
134  r = (x-w)*(fx-fv);
135  q = (x-v)*(fx-fw);
136  p = (x-v)*q - (x-w)*r;
137  q = 2*(q-r);
138  if (q>0) p=-p;
139  else q=-q;
140  r=e; // Deltax before last
141  e=d; // last delta x
142  // current Deltax = p/q
143  // take a parabolic step only if:
144  // Deltax < 0.5* (DeltaX before last) && Deltax > a && Deltax < b
145  // (a BUG in testing this condition is fixed 11/3/2010 (with revision 32544)
146  if (std::fabs(p) >= std::fabs(0.5*q*r) || p <= q*(a-x) || p >= q*(b-x)) {
147  // condition fails - do not take parabolic step
148  e=(x>=m ? a-x : b-x);
149  d = c*e;
150  } else {
151  // take a parabolic interpolation step
152  d = p/q;
153  u = x+d;
154  if (u-a < t2 || b-u < t2)
155  //d=TMath::Sign(tol, m-x);
156  d=(m-x >= 0) ? std::fabs(tol) : -std::fabs(tol);
157  }
158  } else {
159  e=(x>=m ? a-x : b-x);
160  d = c*e;
161  }
162  u = (std::fabs(d)>=tol ? x+d : x+ ((d >= 0) ? std::fabs(tol) : -std::fabs(tol)) );
163  if (type < 2)
164  fu = (*function)(u);
165  else if (type < 4)
166  fu = -(*function)(u);
167  else
168  fu = std::fabs((*function)(u)-fy);
169  //update a, b, v, w and x
170  if (fu<=fx){
171  if (u<x) b=x;
172  else a=x;
173  v=w; fv=fw; w=x; fw=fx; x=u; fx=fu;
174  } else {
175  if (u<x) a=u;
176  else b=u;
177  if (fu<=fw || w==x){
178  v=w; fv=fw; w=u; fw=fu;
179  }
180  else if (fu<=fv || v==x || v==w){
181  v=u; fv=fu;
182  }
183  }
184  }
185  //didn't converge
186  ok = false;
187  xmin = a;
188  xmax = b;
189  niter = itermax;
190  return x;
191 
192 }
193 
194 } // end namespace BrentMethods
195 } // end namespace Math
196 } // 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:133
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
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...
return c
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
TRandom2 r(17)
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
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
double exp(double)
float * q
Definition: THbookFile.cxx:87
double log(double)