Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
22Implement the abstract 1-dimensional root finding interface using
23the Brent-Decker method. This implementation is based on the one
24in the GNU scientific library (v0.99).
25**/
26
27#include "RooBrentRootFinder.h"
28#include "RooAbsFunc.h"
29#include <cmath>
30#include "Riostream.h"
31#include "RooMsgService.h"
32
33using namespace std;
34
36
37
38////////////////////////////////////////////////////////////////////////////////
39/// Constructor taking function binding as input
40
42 _function(&function), _valid(function.isValid()),
43 _tol(2.2204460492503131e-16)
44{
45 if(_function->getDimension() != 1) {
46 oocoutE(nullptr,Eval) << "RooBrentRootFinder:: cannot find roots for function of dimension "
47 << _function->getDimension() << endl;
48 _valid= false;
49 }
50}
51
52
53
54////////////////////////////////////////////////////////////////////////////////
55/// Do the root finding using the Brent-Decker method. Returns a boolean status and
56/// loads 'result' with our best guess at the root if true.
57/// Prints a warning if the initial interval does not bracket a single
58/// root or if the root is not found after a fixed number of iterations.
59
60bool RooBrentRootFinder::findRoot(double &result, double xlo, double xhi, double value) const
61{
63
64 double a(xlo);
65 double b(xhi);
66 double fa= (*_function)(&a) - value;
67 double fb= (*_function)(&b) - value;
68 if(fb*fa > 0) {
69 oocxcoutD((TObject*)nullptr,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): initial interval does not bracket a root: ("
70 << a << "," << b << "), value = " << value << " f[xlo] = " << fa << " f[xhi] = " << fb << endl;
71 return false;
72 }
73
74 bool ac_equal(false);
75 double fc= fb;
76 double c(0);
77 double d(0);
78 double e(0);
79 for(Int_t iter= 0; iter <= MaxIterations; iter++) {
80
81 if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
82 // Rename a,b,c and adjust bounding interval d
83 ac_equal = true;
84 c = a;
85 fc = fa;
86 d = b - a;
87 e = b - a;
88 }
89
90 if (std::abs(fc) < std::abs(fb)) {
91 ac_equal = true;
92 a = b;
93 b = c;
94 c = a;
95 fa = fb;
96 fb = fc;
97 fc = fa;
98 }
99
100 double tol = 0.5 * _tol * std::abs(b);
101 double m = 0.5 * (c - b);
102
103
104 if (fb == 0 || std::abs(m) <= tol) {
105 //cout << "RooBrentRootFinder: iter = " << iter << " m = " << m << " tol = " << tol << endl ;
106 result= b;
108 return true;
109 }
110
111 if (std::abs(e) < tol || std::abs(fa) <= std::abs(fb)) {
112 // Bounds decreasing too slowly: use bisection
113 d = m;
114 e = m;
115 }
116 else {
117 // Attempt inverse cubic interpolation
118 double p;
119 double q;
120 double r;
121 double s = fb / fa;
122
123 if (ac_equal) {
124 p = 2 * m * s;
125 q = 1 - s;
126 }
127 else {
128 q = fa / fc;
129 r = fb / fc;
130 p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
131 q = (q - 1) * (r - 1) * (s - 1);
132 }
133 // Check whether we are in bounds
134 if (p > 0) {
135 q = -q;
136 }
137 else {
138 p = -p;
139 }
140
141 double min1= 3 * m * q - std::abs(tol * q);
142 double min2= std::abs(e * q);
143 if (2 * p < (min1 < min2 ? min1 : min2)) {
144 // Accept the interpolation
145 e = d;
146 d = p / q;
147 }
148 else {
149 // Interpolation failed: use bisection.
150 d = m;
151 e = m;
152 }
153 }
154 // Move last best guess to a
155 a = b;
156 fa = fb;
157 // Evaluate new trial root
158 if (std::abs(d) > tol) {
159 b += d;
160 }
161 else {
162 b += (m > 0 ? +tol : -tol);
163 }
164 fb= (*_function)(&b) - value;
165
166 }
167 // Return our best guess if we run out of iterations
168 oocoutE(nullptr,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): maximum iterations exceeded." << endl;
169 result= b;
170
172
173 return false;
174}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
RooAbsReal & function()
#define oocxcoutD(o, a)
#define oocoutE(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
float * q
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition RooAbsFunc.h:27
virtual void saveXVec() const
Definition RooAbsFunc.h:56
virtual void restoreXVec() const
Definition RooAbsFunc.h:59
UInt_t getDimension() const
Definition RooAbsFunc.h:33
virtual const char * getName() const
Name of function binding.
Definition RooAbsFunc.h:65
Implement the abstract 1-dimensional root finding interface using the Brent-Decker method.
const RooAbsFunc * _function
Pointer to input function.
bool findRoot(double &result, double xlo, double xhi, double value=0) const
Do the root finding using the Brent-Decker method.
static constexpr int MaxIterations
RooBrentRootFinder(const RooAbsFunc &function)
Constructor taking function binding as input.
bool _valid
True if current state is valid.
Mother of all ROOT objects.
Definition TObject.h:41
TMarker m
Definition textangle.C:8