TOMS748.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /// \file
5 /// Declares function RootFinder::toms748
6 
7 #pragma once
8 
9 #include <boost/math/tools/roots.hpp>
10 #include <functional>
11 #include <limits>
12 
13 #include "DataStructures/DataVector.hpp"
14 #include "ErrorHandling/Exceptions.hpp"
15 
16 namespace RootFinder {
17 
18 /*!
19  * \ingroup NumericalAlgorithmsGroup
20  * \brief Finds the root of the function `f` with the TOMS_748 method.
21  *
22  * `f` is a unary invokable that takes a `double` which is the current value at
23  * which to evaluate `f`. An example is below.
24  *
25  * \snippet Test_TOMS748.cpp double_root_find
26  *
27  * The TOMS_748 algorithm searches for a root in the interval [`lower_bound`,
28  * `upper_bound`], and will throw if this interval does not bracket a root,
29  * i.e. if `f(lower_bound) * f(upper_bound) > 0`.
30  *
31  * See the [Boost](http://www.boost.org/) documentation for more details.
32  *
33  * \requires Function `f` is invokable with a `double`
34  *
35  * \throws `std::domain_error` if the bounds do not bracket a root.
36  * \throws `convergence_error` if the requested tolerance is not met after
37  * `max_iterations` iterations.
38  */
39 template <typename Function>
40 double toms748(const Function& f, const double lower_bound,
41  const double upper_bound, const double absolute_tolerance,
42  const double relative_tolerance,
43  const size_t max_iterations = 100) {
44  ASSERT(relative_tolerance > std::numeric_limits<double>::epsilon(),
45  "The relative tolerance is too small.");
46 
47  boost::uintmax_t max_iters = max_iterations;
48 
49  // This solver requires tol to be passed as a termination condition. This
50  // termination condition is equivalent to the convergence criteria used by the
51  // GSL
52  auto tol = [absolute_tolerance, relative_tolerance](double lhs, double rhs) {
53  return (fabs(lhs - rhs) <=
54  absolute_tolerance +
55  relative_tolerance * fmin(fabs(lhs), fabs(rhs)));
56  };
57  // clang-tidy: internal boost warning, can't fix it.
58  auto result = boost::math::tools::toms748_solve( // NOLINT
59  f, lower_bound, upper_bound, tol, max_iters);
60  if (max_iters >= max_iterations) {
61  throw convergence_error(
62  "toms748 reached max iterations without converging");
63  }
64  return result.first + 0.5 * (result.second - result.first);
65 }
66 
67 /*!
68  * \ingroup NumericalAlgorithmsGroup
69  * \brief Finds the root of the function `f` with the TOMS_748 method on each
70  * element in a `DataVector`.
71  *
72  * `f` is a binary invokable that takes a `double` as its first argument and a
73  * `size_t` as its second. The `double` is the current value at which to
74  * evaluate `f`, and the `size_t` is the current index into the `DataVector`s.
75  * Below is an example of how to root find different functions by indexing into
76  * a lambda-captured `DataVector` using the `size_t` passed to `f`.
77  *
78  * \snippet Test_TOMS748.cpp datavector_root_find
79  *
80  * For each index `i` into the DataVector, the TOMS_748 algorithm searches for a
81  * root in the interval [`lower_bound[i]`, `upper_bound[i]`], and will throw if
82  * this interval does not bracket a root,
83  * i.e. if `f(lower_bound[i], i) * f(upper_bound[i], i) > 0`.
84  *
85  * See the [Boost](http://www.boost.org/) documentation for more details.
86  *
87  * \requires Function `f` be callable with a `double` and a `size_t`
88  *
89  * \throws `std::domain_error` if, for any index, the bounds do not bracket a
90  * root.
91  * \throws `convergence_error` if, for any index, the requested tolerance is not
92  * met after `max_iterations` iterations.
93  */
94 template <typename Function>
95 DataVector toms748(const Function& f, const DataVector& lower_bound,
96  const DataVector& upper_bound,
97  const double absolute_tolerance,
98  const double relative_tolerance,
99  const size_t max_iterations = 100) {
100  ASSERT(relative_tolerance > std::numeric_limits<double>::epsilon(),
101  "The relative tolerance is too small.");
102  // This solver requires tol to be passed as a termination condition. This
103  // termination condition is equivalent to the convergence criteria used by the
104  // GSL
105  auto tol = [absolute_tolerance, relative_tolerance](const double lhs,
106  const double rhs) {
107  return (fabs(lhs - rhs) <=
108  absolute_tolerance +
109  relative_tolerance * fmin(fabs(lhs), fabs(rhs)));
110  };
111  DataVector result_vector{lower_bound.size()};
112  for (size_t i = 0; i < result_vector.size(); ++i) {
113  // toms748_solver modifies the max_iters after the root is found to the
114  // number of iterations that it took to find the root, so we reset it to
115  // max_iterations after each root find.
116  boost::uintmax_t max_iters = max_iterations;
117  // clang-tidy: internal boost warning, can't fix it.
118  auto result = boost::math::tools::toms748_solve( // NOLINT
119  [&f, i](double x) { return f(x, i); }, lower_bound[i], upper_bound[i],
120  tol, max_iters);
121  if (max_iters >= max_iterations) {
122  throw convergence_error(
123  "toms748 reached max iterations without converging");
124  }
125  result_vector[i] = result.first + 0.5 * (result.second - result.first);
126  }
127  return result_vector;
128 }
129 
130 } // namespace RootFinder
Definition: GslMultiRoot.cpp:10
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Exception indicating convergence failure.
Definition: Exceptions.hpp:10
DataVector toms748(const Function &f, const DataVector &lower_bound, const DataVector &upper_bound, const double absolute_tolerance, const double relative_tolerance, const size_t max_iterations=100)
Finds the root of the function f with the TOMS_748 method on each element in a DataVector.
Definition: TOMS748.hpp:95
Stores a collection of function values.
Definition: DataVector.hpp:46