TOMS748.hpp
Go to the documentation of this file.
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 DataVectors.
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