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/Error.h"
19
20#include <iostream>
21
22namespace ROOT {
23namespace Math {
24
25namespace 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) {
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
#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
ROOT::R::TRInterface & r
Definition: Object.C:4
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
int type
Definition: TGX11.cxx:120
float xmin
Definition: THbookFile.cxx:93
float * q
Definition: THbookFile.cxx:87
float xmax
Definition: THbookFile.cxx:93
double exp(double)
double log(double)
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
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...
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:151
VSD Structures.
Definition: StringConv.hxx:21
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12