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
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 "RooAbsFunc.h"
31#include <math.h>
32#include "Riostream.h"
33#include "RooMsgService.h"
34
35using 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{
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;
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}
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:83
#define oocoutE(o, a)
Definition: RooMsgService.h:48
int Int_t
Definition: RtypesCore.h:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
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
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12