Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Declares function RootFinder::newton_raphson
       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 "Utilities/ErrorHandling/Exceptions.hpp"
      15             : #include "Utilities/MakeString.hpp"
      16             : 
      17             : namespace RootFinder {
      18             : /*!
      19             :  * \ingroup NumericalAlgorithmsGroup
      20             :  * \brief Finds the root of the function `f` with the Newton-Raphson method.
      21             :  *
      22             :  * `f` is a unary invokable that takes a `double` which is the current value at
      23             :  * which to evaluate `f`. `f` must return a `std::pair<double, double>` where
      24             :  * the first element is the function value and the second element is the
      25             :  * derivative of the function.  An example is below.
      26             :  *
      27             :  * \snippet Test_NewtonRaphson.cpp double_newton_raphson_root_find
      28             :  *
      29             :  * See the [Boost](http://www.boost.org/) documentation for more details.
      30             :  *
      31             :  * \requires Function `f` is invokable with a `double`
      32             :  * \note The parameter `digits` specifies the precision of the result in its
      33             :  * desired number of base-10 digits.
      34             :  *
      35             :  * \throws `convergence_error` if the requested precision is not met after
      36             :  * `max_iterations` iterations.
      37             :  */
      38             : template <typename Function>
      39           1 : double newton_raphson(const Function& f, const double initial_guess,
      40             :                       const double lower_bound, const double upper_bound,
      41             :                       const size_t digits, const size_t max_iterations = 50) {
      42             :   ASSERT(digits < std::numeric_limits<double>::digits10,
      43             :          "The desired accuracy of " << digits
      44             :                                     << " base-10 digits must be smaller than "
      45             :                                        "the machine numeric limit of "
      46             :                                     << std::numeric_limits<double>::digits10
      47             :                                     << " base-10 digits.");
      48             : 
      49             :   boost::uintmax_t max_iters = max_iterations;
      50             :   // clang-tidy: internal boost warning, can't fix it.
      51             :   const auto result = boost::math::tools::newton_raphson_iterate(  // NOLINT
      52             :       f, initial_guess, lower_bound, upper_bound,
      53             :       std::round(std::log2(std::pow(10, digits))), max_iters);
      54             :   if (max_iters >= max_iterations) {
      55             :     throw convergence_error(MakeString{}
      56             :                             << "newton_raphson reached max iterations of "
      57             :                             << max_iterations
      58             :                             << " without converging. Best result is: " << result
      59             :                             << " with residual " << f(result).first);
      60             :   }
      61             :   return result;
      62             : }
      63             : 
      64             : /*!
      65             :  * \ingroup NumericalAlgorithmsGroup
      66             :  * \brief Finds the root of the function `f` with the Newton-Raphson method on
      67             :  * each element in a `DataVector`.
      68             :  *
      69             :  * `f` is a binary invokable that takes a `double` as its first argument and a
      70             :  * `size_t` as its second. The `double` is the current value at which to
      71             :  * evaluate `f`, and the `size_t` is the current index into the `DataVector`s.
      72             :  *  `f` must return a `std::pair<double, double>` where the first element is
      73             :  * the function value and the second element is the derivative of the function.
      74             :  * Below is an example of how to root find different functions by indexing into
      75             :  * a lambda-captured `DataVector` using the `size_t` passed to `f`.
      76             :  *
      77             :  * \snippet Test_NewtonRaphson.cpp datavector_newton_raphson_root_find
      78             :  *
      79             :  * See the [Boost](http://www.boost.org/) documentation for more details.
      80             :  *
      81             :  * \requires Function `f` be callable with a `double` and a `size_t`
      82             :  * \note The parameter `digits` specifies the precision of the result in its
      83             :  * desired number of base-10 digits.
      84             :  *
      85             :  * \throws `convergence_error` if, for any index, the requested precision is not
      86             :  * met after `max_iterations` iterations.
      87             :  */
      88             : template <typename Function>
      89           1 : DataVector newton_raphson(const Function& f, const DataVector& initial_guess,
      90             :                           const DataVector& lower_bound,
      91             :                           const DataVector& upper_bound, const size_t digits,
      92             :                           const size_t max_iterations = 50) {
      93             :   ASSERT(digits < std::numeric_limits<double>::digits10,
      94             :          "The desired accuracy of " << digits
      95             :                                     << " base-10 digits must be smaller than "
      96             :                                        "the machine numeric limit of "
      97             :                                     << std::numeric_limits<double>::digits10
      98             :                                     << " base-10 digits.");
      99             :   const auto digits_binary = std::round(std::log2(std::pow(10, digits)));
     100             : 
     101             :   DataVector result_vector{lower_bound.size()};
     102             :   for (size_t i = 0; i < result_vector.size(); ++i) {
     103             :     boost::uintmax_t max_iters = max_iterations;
     104             :     // clang-tidy: internal boost warning, can't fix it.
     105             :     result_vector[i] = boost::math::tools::newton_raphson_iterate(  // NOLINT
     106             :         [&f, i](double x) { return f(x, i); }, initial_guess[i], lower_bound[i],
     107             :         upper_bound[i], digits_binary, max_iters);
     108             :     if (max_iters >= max_iterations) {
     109             :       throw convergence_error(MakeString{}
     110             :                               << "newton_raphson reached max iterations of "
     111             :                               << max_iterations
     112             :                               << " without converging. Best result is: "
     113             :                               << result_vector[i] << " with residual "
     114             :                               << f(result_vector[i], i).first);
     115             :     }
     116             :   }
     117             :   return result_vector;
     118             : }
     119             : 
     120             : }  // namespace RootFinder