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