Logo ROOT  
Reference Guide
RooBrentRootFinder.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitCore *
4  * @(#)root/roofitcore:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16 
17 /**
18 \file RooBrentRootFinder.cxx
19 \class RooBrentRootFinder
20 \ingroup Roofitcore
21 
22 Implement the abstract 1-dimensional root finding interface using
23 the Brent-Decker method. This implementation is based on the one
24 in the GNU scientific library (v0.99).
25 **/
26 
27 #include "RooFit.h"
28 
29 #include "RooBrentRootFinder.h"
30 #include "RooAbsFunc.h"
31 #include <math.h>
32 #include "Riostream.h"
33 #include "RooMsgService.h"
34 
35 using namespace std;
36 
38 ;
39 
40 
41 ////////////////////////////////////////////////////////////////////////////////
42 /// Constructor taking function binding as input
43 
46  _tol(2.2204460492503131e-16)
47 {
48 }
49 
50 
51 
52 ////////////////////////////////////////////////////////////////////////////////
53 /// Do the root finding using the Brent-Decker method. Returns a boolean status and
54 /// loads 'result' with our best guess at the root if true.
55 /// Prints a warning if the initial interval does not bracket a single
56 /// root or if the root is not found after a fixed number of iterations.
57 
59 {
60  _function->saveXVec() ;
61 
62  Double_t a(xlo),b(xhi);
63  Double_t fa= (*_function)(&a) - value;
64  Double_t fb= (*_function)(&b) - value;
65  if(fb*fa > 0) {
66  oocxcoutD((TObject*)0,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): initial interval does not bracket a root: ("
67  << a << "," << b << "), value = " << value << " f[xlo] = " << fa << " f[xhi] = " << fb << endl;
68  return kFALSE;
69  }
70 
71  Bool_t ac_equal(kFALSE);
72  Double_t fc= fb;
73  Double_t c(0),d(0),e(0);
74  for(Int_t iter= 0; iter <= MaxIterations; iter++) {
75 
76  if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
77  // Rename a,b,c and adjust bounding interval d
78  ac_equal = kTRUE;
79  c = a;
80  fc = fa;
81  d = b - a;
82  e = b - a;
83  }
84 
85  if (fabs (fc) < fabs (fb)) {
86  ac_equal = kTRUE;
87  a = b;
88  b = c;
89  c = a;
90  fa = fb;
91  fb = fc;
92  fc = fa;
93  }
94 
95  Double_t tol = 0.5 * _tol * fabs(b);
96  Double_t m = 0.5 * (c - b);
97 
98 
99  if (fb == 0 || fabs(m) <= tol) {
100  //cout << "RooBrentRootFinder: iter = " << iter << " m = " << m << " tol = " << tol << endl ;
101  result= b;
102  _function->restoreXVec() ;
103  return kTRUE;
104  }
105 
106  if (fabs (e) < tol || fabs (fa) <= fabs (fb)) {
107  // Bounds decreasing too slowly: use bisection
108  d = m;
109  e = m;
110  }
111  else {
112  // Attempt inverse cubic interpolation
113  Double_t p, q, r;
114  Double_t s = fb / fa;
115 
116  if (ac_equal) {
117  p = 2 * m * s;
118  q = 1 - s;
119  }
120  else {
121  q = fa / fc;
122  r = fb / fc;
123  p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
124  q = (q - 1) * (r - 1) * (s - 1);
125  }
126  // Check whether we are in bounds
127  if (p > 0) {
128  q = -q;
129  }
130  else {
131  p = -p;
132  }
133 
134  Double_t min1= 3 * m * q - fabs (tol * q);
135  Double_t min2= fabs (e * q);
136  if (2 * p < (min1 < min2 ? min1 : min2)) {
137  // Accept the interpolation
138  e = d;
139  d = p / q;
140  }
141  else {
142  // Interpolation failed: use bisection.
143  d = m;
144  e = m;
145  }
146  }
147  // Move last best guess to a
148  a = b;
149  fa = fb;
150  // Evaluate new trial root
151  if (fabs (d) > tol) {
152  b += d;
153  }
154  else {
155  b += (m > 0 ? +tol : -tol);
156  }
157  fb= (*_function)(&b) - value;
158 
159  }
160  // Return our best guess if we run out of iterations
161  oocoutE((TObject*)0,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): maximum iterations exceeded." << endl;
162  result= b;
163 
165 
166  return kFALSE;
167 }
c
#define c(i)
Definition: RSha256.hxx:119
m
auto * m
Definition: textangle.C:8
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
e
#define e(i)
Definition: RSha256.hxx:121
RooMsgService.h
RooFit.h
fc
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
RooAbsFunc::restoreXVec
virtual void restoreXVec() const
Definition: RooAbsFunc.h:54
r
ROOT::R::TRInterface & r
Definition: Object.C:4
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:168
RooAbsFunc.h
RooAbsFunc::saveXVec
virtual void saveXVec() const
Definition: RooAbsFunc.h:51
oocoutE
#define oocoutE(o, a)
Definition: RooMsgService.h:48
b
#define b(i)
Definition: RSha256.hxx:118
bool
RooBrentRootFinder::MaxIterations
@ MaxIterations
Definition: RooBrentRootFinder.h:47
RooAbsFunc
Definition: RooAbsFunc.h:23
q
float * q
Definition: THbookFile.cxx:89
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
oocxcoutD
#define oocxcoutD(o, a)
Definition: RooMsgService.h:83
a
auto * a
Definition: textangle.C:12
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
ROOT::R::function
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:151
RooBrentRootFinder
Definition: RooBrentRootFinder.h:21
RooBrentRootFinder::_tol
Double_t _tol
Definition: RooBrentRootFinder.h:49
RooBrentRootFinder::RooBrentRootFinder
RooBrentRootFinder(const RooAbsFunc &function)
Constructor taking function binding as input.
Definition: RooBrentRootFinder.cxx:44
Double_t
double Double_t
Definition: RtypesCore.h:59
RooBrentRootFinder.h
TObject
Definition: TObject.h:37
d
#define d(i)
Definition: RSha256.hxx:120
RooAbsFunc::getName
virtual const char * getName() const
Definition: RooAbsFunc.h:59
RooFit::Eval
@ Eval
Definition: RooGlobalFunc.h:68
RooAbsRootFinder
Definition: RooAbsRootFinder.h:23
Riostream.h
RooBrentRootFinder::findRoot
virtual Bool_t findRoot(Double_t &result, Double_t xlo, Double_t xhi, Double_t value=0) const
Do the root finding using the Brent-Decker method.
Definition: RooBrentRootFinder.cxx:58
RooAbsRootFinder::_function
const RooAbsFunc * _function
Definition: RooAbsRootFinder.h:31
int