1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file solver1d.hpp
21 \brief Abstract 1-D solver class
22*/
23
24#ifndef quantlib_solver1d_hpp
25#define quantlib_solver1d_hpp
26
27#include <ql/math/comparison.hpp>
28#include <ql/utilities/null.hpp>
29#include <ql/patterns/curiouslyrecurring.hpp>
30#include <ql/errors.hpp>
31#include <algorithm>
32#include <iomanip>
33
34namespace QuantLib {
35
36 #define MAX_FUNCTION_EVALUATIONS 100
37
38 //! Base class for 1-D solvers
39 /*! The implementation of this class uses the so-called
40 "Barton-Nackman trick", also known as "the curiously recurring
41 template pattern". Concrete solvers will be declared as:
42 \code
43 class Foo : public Solver1D<Foo> {
44 public:
45 ...
46 template <class F>
47 Real solveImpl(const F& f, Real accuracy) const {
48 ...
49 }
50 };
51 \endcode
52 Before calling <tt>solveImpl</tt>, the base class will set its
53 protected data members so that:
54 - <tt>xMin_</tt> and <tt>xMax_</tt> form a valid bracket;
55 - <tt>fxMin_</tt> and <tt>fxMax_</tt> contain the values of
56 the function in <tt>xMin_</tt> and <tt>xMax_</tt>;
57 - <tt>root_</tt> is a valid initial guess.
58 The implementation of <tt>solveImpl</tt> can safely assume all
59 of the above.
60
61 \todo
62 - clean up the interface so that it is clear whether the
63 accuracy is specified for \f$ x \f$ or \f$ f(x) \f$.
64 - add target value (now the target value is 0.0)
65 */
66 template <class Impl>
67 class Solver1D : public CuriouslyRecurringTemplate<Impl> {
68 public:
69 Solver1D() = default;
70 //! \name Modifiers
71 //@{
72 /*! This method returns the zero of the function \f$ f \f$,
73 determined with the given accuracy \f$ \epsilon \f$;
74 depending on the particular solver, this might mean that
75 the returned \f$ x \f$ is such that \f$ |f(x)| < \epsilon
76 \f$, or that \f$ |x-\xi| < \epsilon \f$ where \f$ \xi \f$
77 is the real zero.
78
79 This method contains a bracketing routine to which an
80 initial guess must be supplied as well as a step used to
81 scan the range of the possible bracketing values.
82 */
83 template <class F>
84 Real solve(const F& f,
85 Real accuracy,
86 Real guess,
87 Real step) const {
88
89 QL_REQUIRE(accuracy>0.0,
90 "accuracy (" << accuracy << ") must be positive");
91 // check whether we really want to use epsilon
92 accuracy = std::max(a: accuracy, QL_EPSILON);
93
94 const Real growthFactor = 1.6;
95 Integer flipflop = -1;
96
97 root_ = guess;
98 fxMax_ = f(root_);
99
100 // monotonically crescent bias, as in optionValue(volatility)
101 if (close(x: fxMax_,y: 0.0))
102 return root_;
103 else if (fxMax_ > 0.0) {
104 xMin_ = enforceBounds_(x: root_ - step);
105 fxMin_ = f(xMin_);
106 xMax_ = root_;
107 } else {
108 xMin_ = root_;
109 fxMin_ = fxMax_;
110 xMax_ = enforceBounds_(x: root_+step);
111 fxMax_ = f(xMax_);
112 }
113
114 evaluationNumber_ = 2;
115 while (evaluationNumber_ <= maxEvaluations_) {
116 if (fxMin_*fxMax_ <= 0.0) {
117 if (close(x: fxMin_, y: 0.0))
118 return xMin_;
119 if (close(x: fxMax_, y: 0.0))
120 return xMax_;
121 root_ = (xMax_+xMin_)/2.0;
122 return this->impl().solveImpl(f, accuracy);
123 }
124 if (std::fabs(x: fxMin_) < std::fabs(x: fxMax_)) {
125 xMin_ = enforceBounds_(x: xMin_+growthFactor*(xMin_ - xMax_));
126 fxMin_= f(xMin_);
127 } else if (std::fabs(x: fxMin_) > std::fabs(x: fxMax_)) {
128 xMax_ = enforceBounds_(x: xMax_+growthFactor*(xMax_ - xMin_));
129 fxMax_= f(xMax_);
130 } else if (flipflop == -1) {
131 xMin_ = enforceBounds_(x: xMin_+growthFactor*(xMin_ - xMax_));
132 fxMin_= f(xMin_);
133 evaluationNumber_++;
134 flipflop = 1;
135 } else if (flipflop == 1) {
136 xMax_ = enforceBounds_(x: xMax_+growthFactor*(xMax_ - xMin_));
137 fxMax_= f(xMax_);
138 flipflop = -1;
139 }
140 evaluationNumber_++;
141 }
142
143 QL_FAIL("unable to bracket root in " << maxEvaluations_
144 << " function evaluations (last bracket attempt: "
145 << "f[" << xMin_ << "," << xMax_ << "] "
146 << "-> [" << fxMin_ << "," << fxMax_ << "])");
147 }
148 /*! This method returns the zero of the function \f$ f \f$,
149 determined with the given accuracy \f$ \epsilon \f$;
150 depending on the particular solver, this might mean that
151 the returned \f$ x \f$ is such that \f$ |f(x)| < \epsilon
152 \f$, or that \f$ |x-\xi| < \epsilon \f$ where \f$ \xi \f$
153 is the real zero.
154
155 An initial guess must be supplied, as well as two values
156 \f$ x_\mathrm{min} \f$ and \f$ x_\mathrm{max} \f$ which
157 must bracket the zero (i.e., either \f$ f(x_\mathrm{min})
158 \leq 0 \leq f(x_\mathrm{max}) \f$, or \f$
159 f(x_\mathrm{max}) \leq 0 \leq f(x_\mathrm{min}) \f$ must
160 be true).
161 */
162 template <class F>
163 Real solve(const F& f,
164 Real accuracy,
165 Real guess,
166 Real xMin,
167 Real xMax) const {
168
169 QL_REQUIRE(accuracy>0.0,
170 "accuracy (" << accuracy << ") must be positive");
171 // check whether we really want to use epsilon
172 accuracy = std::max(a: accuracy, QL_EPSILON);
173
174 xMin_ = xMin;
175 xMax_ = xMax;
176
177 QL_REQUIRE(xMin_ < xMax_,
178 "invalid range: xMin_ (" << xMin_
179 << ") >= xMax_ (" << xMax_ << ")");
180 QL_REQUIRE(!lowerBoundEnforced_ || xMin_ >= lowerBound_,
181 "xMin_ (" << xMin_
182 << ") < enforced low bound (" << lowerBound_ << ")");
183 QL_REQUIRE(!upperBoundEnforced_ || xMax_ <= upperBound_,
184 "xMax_ (" << xMax_
185 << ") > enforced hi bound (" << upperBound_ << ")");
186
187 fxMin_ = f(xMin_);
188 if (close(x: fxMin_, y: 0.0))
189 return xMin_;
190
191 fxMax_ = f(xMax_);
192 if (close(x: fxMax_, y: 0.0))
193 return xMax_;
194
195 evaluationNumber_ = 2;
196
197 QL_REQUIRE(fxMin_*fxMax_ < 0.0,
198 "root not bracketed: f["
199 << xMin_ << "," << xMax_ << "] -> ["
200 << std::scientific
201 << fxMin_ << "," << fxMax_ << "]");
202
203 QL_REQUIRE(guess > xMin_,
204 "guess (" << guess << ") < xMin_ (" << xMin_ << ")");
205 QL_REQUIRE(guess < xMax_,
206 "guess (" << guess << ") > xMax_ (" << xMax_ << ")");
207
208 root_ = guess;
209
210 return this->impl().solveImpl(f, accuracy);
211 }
212
213 /*! This method sets the maximum number of function
214 evaluations for the bracketing routine. An error is thrown
215 if a bracket is not found after this number of
216 evaluations.
217 */
218 void setMaxEvaluations(Size evaluations);
219 //! sets the lower bound for the function domain
220 void setLowerBound(Real lowerBound);
221 //! sets the upper bound for the function domain
222 void setUpperBound(Real upperBound);
223 //@}
224 protected:
225 mutable Real root_, xMin_, xMax_, fxMin_, fxMax_;
226 Size maxEvaluations_ = MAX_FUNCTION_EVALUATIONS;
227 mutable Size evaluationNumber_;
228 private:
229 Real enforceBounds_(Real x) const;
230 Real lowerBound_, upperBound_;
231 bool lowerBoundEnforced_ = false, upperBoundEnforced_ = false;
232 };
233
234
235 // inline definitions
236
237 template <class T>
238 inline void Solver1D<T>::setMaxEvaluations(Size evaluations) {
239 maxEvaluations_ = evaluations;
240 }
241
242 template <class T>
243 inline void Solver1D<T>::setLowerBound(Real lowerBound) {
244 lowerBound_ = lowerBound;
245 lowerBoundEnforced_ = true;
246 }
247
248 template <class T>
249 inline void Solver1D<T>::setUpperBound(Real upperBound) {
250 upperBound_ = upperBound;
251 upperBoundEnforced_ = true;
252 }
253
254 template <class T>
255 inline Real Solver1D<T>::enforceBounds_(Real x) const {
256 if (lowerBoundEnforced_ && x < lowerBound_)
257 return lowerBound_;
258 if (upperBoundEnforced_ && x > upperBound_)
259 return upperBound_;
260 return x;
261 }
262
263}
264
265#endif
266

source code of quantlib/ql/math/solver1d.hpp