Logo ROOT   6.16/01
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
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 "RooFit.h"
28
29#include "RooBrentRootFinder.h"
30#include "RooBrentRootFinder.h"
31#include "RooAbsFunc.h"
32#include <math.h>
33#include "Riostream.h"
34#include "RooMsgService.h"
35
36using namespace std;
37
39;
40
41
42////////////////////////////////////////////////////////////////////////////////
43/// Constructor taking function binding as input
44
47 _tol(2.2204460492503131e-16)
48{
49}
50
51
52
53////////////////////////////////////////////////////////////////////////////////
54/// Do the root finding using the Brent-Decker method. Returns a boolean status and
55/// loads 'result' with our best guess at the root if true.
56/// Prints a warning if the initial interval does not bracket a single
57/// root or if the root is not found after a fixed number of iterations.
58
60{
62
63 Double_t a(xlo),b(xhi);
64 Double_t fa= (*_function)(&a) - value;
65 Double_t fb= (*_function)(&b) - value;
66 if(fb*fa > 0) {
67 oocxcoutD((TObject*)0,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): initial interval does not bracket a root: ("
68 << a << "," << b << "), value = " << value << " f[xlo] = " << fa << " f[xhi] = " << fb << endl;
69 return kFALSE;
70 }
71
72 Bool_t ac_equal(kFALSE);
73 Double_t fc= fb;
74 Double_t c(0),d(0),e(0);
75 for(Int_t iter= 0; iter <= MaxIterations; iter++) {
76
77 if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
78 // Rename a,b,c and adjust bounding interval d
79 ac_equal = kTRUE;
80 c = a;
81 fc = fa;
82 d = b - a;
83 e = b - a;
84 }
85
86 if (fabs (fc) < fabs (fb)) {
87 ac_equal = kTRUE;
88 a = b;
89 b = c;
90 c = a;
91 fa = fb;
92 fb = fc;
93 fc = fa;
94 }
95
96 Double_t tol = 0.5 * _tol * fabs(b);
97 Double_t m = 0.5 * (c - b);
98
99
100 if (fb == 0 || fabs(m) <= tol) {
101 //cout << "RooBrentRootFinder: iter = " << iter << " m = " << m << " tol = " << tol << endl ;
102 result= b;
104 return kTRUE;
105 }
106
107 if (fabs (e) < tol || fabs (fa) <= fabs (fb)) {
108 // Bounds decreasing too slowly: use bisection
109 d = m;
110 e = m;
111 }
112 else {
113 // Attempt inverse cubic interpolation
114 Double_t p, q, r;
115 Double_t s = fb / fa;
116
117 if (ac_equal) {
118 p = 2 * m * s;
119 q = 1 - s;
120 }
121 else {
122 q = fa / fc;
123 r = fb / fc;
124 p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
125 q = (q - 1) * (r - 1) * (s - 1);
126 }
127 // Check whether we are in bounds
128 if (p > 0) {
129 q = -q;
130 }
131 else {
132 p = -p;
133 }
134
135 Double_t min1= 3 * m * q - fabs (tol * q);
136 Double_t min2= fabs (e * q);
137 if (2 * p < (min1 < min2 ? min1 : min2)) {
138 // Accept the interpolation
139 e = d;
140 d = p / q;
141 }
142 else {
143 // Interpolation failed: use bisection.
144 d = m;
145 e = m;
146 }
147 }
148 // Move last best guess to a
149 a = b;
150 fa = fb;
151 // Evaluate new trial root
152 if (fabs (d) > tol) {
153 b += d;
154 }
155 else {
156 b += (m > 0 ? +tol : -tol);
157 }
158 fb= (*_function)(&b) - value;
159
160 }
161 // Return our best guess if we run out of iterations
162 oocoutE((TObject*)0,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): maximum iterations exceeded." << endl;
163 result= b;
164
166
167 return kFALSE;
168}
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
#define oocxcoutD(o, a)
Definition: RooMsgService.h:81
#define oocoutE(o, a)
Definition: RooMsgService.h:47
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:363
float * q
Definition: THbookFile.cxx:87
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
virtual void saveXVec() const
Definition: RooAbsFunc.h:51
virtual void restoreXVec() const
Definition: RooAbsFunc.h:54
virtual const char * getName() const
Definition: RooAbsFunc.h:59
RooAbsRootFinder is the abstract interface for finding roots of real-valued 1-dimensional function th...
const RooAbsFunc * _function
Implement the abstract 1-dimensional root finding interface using the Brent-Decker method.
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.
RooBrentRootFinder(const RooAbsFunc &function)
Constructor taking function binding as input.
Mother of all ROOT objects.
Definition: TObject.h:37
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
static constexpr double s
STL namespace.
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12