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