SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/RootFinding - RootBracketing.hpp Hit Total Coverage
Commit: cd74d65bdc718fd7e344eaec61dc6334dd4d366b Lines: 2 5 40.0 %
Date: 2022-08-12 23:56:47
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <optional>
       7             : #include <sstream>
       8             : #include <stdexcept>
       9             : #include <tuple>
      10             : #include <vector>
      11             : 
      12             : #include "DataStructures/DataVector.hpp"
      13             : #include "Utilities/ErrorHandling/Error.hpp"
      14             : #include "Utilities/Gsl.hpp"
      15             : 
      16             : namespace RootFinder {
      17             : 
      18             : namespace bracketing_detail {
      19             : // Brackets a root, given a functor f(x) that returns a
      20             : // std::optional<double> and given two arrays x and y (with y=f(x))
      21             : // containing points that have already been tried for bracketing.
      22             : //
      23             : // Returns a std::tuple<double, double, double, double>
      24             : // containing {{x1,x2,y1,y2}} where
      25             : // x1 and x2 bracket the root, and y1=f(x1) and y2=f(x2).
      26             : //
      27             : // Note that y might be undefined (i.e. an invalid std::optional)
      28             : // for values of x at or near the endpoints of the interval.  We
      29             : // assume that if f(x1) and f(x2) are valid for some x1 and x2, then
      30             : // f(x) is valid for all x between x1 and x2.
      31             : //
      32             : // Assumes that there is a root between the first and last points of
      33             : // the array.
      34             : // We also assume that there is only one root.
      35             : //
      36             : // So this means we have only 2 possibilities for the validity of the
      37             : // points in the input x,y arrays:
      38             : // 1) All points are invalid, e.g. "X X X X X X X X".
      39             : //    (here X represents an invalid point)
      40             : // 2) All valid points are adjacent, with the same sign, e.g. "X X o o X X"
      41             : //    or "o o X X X" or "X X X o o".
      42             : //    (here o represents a valid point)
      43             : // Note that we assume that all valid points have the same sign; otherwise
      44             : // the caller would have known that the root was bracketed and the caller would
      45             : // not have called bracket_by_contracting.
      46             : //
      47             : // Also note that we exclude the case "o o o o o" (no roots), since the
      48             : // caller would have known that too.  If such a case is found, an error is
      49             : // thrown.  An error is also thrown if the size of the region where the
      50             : // sign changes is so small that the number of iterations is exceeded.
      51             : //
      52             : // For case 1) above, we bisect each pair of points, and call
      53             : // bracket_by_contracting recursively until we find a valid point.
      54             : // For case 2) above, it is sufficent to check for a bracket only
      55             : // between valid and invalid points.  That is, for "X X + + X X" we
      56             : // check only between points 1 and 2 and between points 3 and 4 (where
      57             : // points are numbered starting from zero).  For "+ + X X X" we check
      58             : // only between points 1 and 2.
      59             : template <typename Functor>
      60             : auto bracket_by_contracting(const std::vector<double>& x,
      61             :                             const std::vector<std::optional<double>>& y,
      62             :                             const Functor& f, const size_t level = 0)
      63             :     -> std::tuple<double, double, double, double> {
      64             :   constexpr size_t max_level = 27;
      65             :   if (level > max_level) {
      66             :     std::stringstream ss;
      67             :     ss << "Too many iterations in bracket_by_contracting.  Either the "
      68             :           "region where the root changes sign is so small that we cannot "
      69             :           "find it, or the given interval does not actually bracket a "
      70             :           "root.  The points we are checking are "
      71             :        << x.size() << " almost-evenly-spaced points from " << x.front()
      72             :        << " to " << x.back() << " with difference " << x.back()-x.front();
      73             :     throw std::runtime_error(ss.str());
      74             :   }
      75             : 
      76             :   // First check if we have any valid points.
      77             :   size_t last_valid_index = y.size();
      78             :   for (size_t i = y.size(); i >= 1; --i) {
      79             :     if (y[i - 1].has_value()) {
      80             :       last_valid_index = i - 1;
      81             :       break;
      82             :     }
      83             :   }
      84             : 
      85             :   if (last_valid_index == y.size()) {
      86             :     // No valid points!
      87             : 
      88             :     // Create larger arrays with one point between each of the already
      89             :     // computed points.
      90             :     std::vector<double> bisected_x(x.size() * 2 - 1);
      91             :     std::vector<std::optional<double>> bisected_y(y.size() * 2 - 1);
      92             : 
      93             :     // Copy all even-numbered points in the range.
      94             :     for (size_t i = 0; i < x.size(); ++i) {
      95             :       bisected_x[2 * i] = x[i];
      96             :       bisected_y[2 * i] = y[i];
      97             :     }
      98             : 
      99             :     // Fill midpoints and check for bracket on each one.
     100             :     for (size_t i = 0; i < x.size() - 1; ++i) {
     101             :       bisected_x[2 * i + 1] = x[i] + 0.5 * (x[i + 1] - x[i]);
     102             :       bisected_y[2 * i + 1] = f(bisected_x[2 * i + 1]);
     103             :       if (bisected_y[2 * i + 1].has_value()) {
     104             :         // Valid point! We know that all the other points are
     105             :         // invalid, so we need to check only 3 points in the next
     106             :         // iteration: the new valid point and its neighbors.
     107             :         return bracket_by_contracting({{x[i], bisected_x[2 * i + 1], x[i + 1]}},
     108             :                                       {{y[i], bisected_y[2 * i + 1], y[i + 1]}},
     109             :                                       f, level + 1);
     110             :       }
     111             :     }
     112             :     // We still have no valid points. So recurse, using all points.
     113             :     // The next iteration will bisect all the points.
     114             :     return bracket_by_contracting(bisected_x, bisected_y, f, level + 1);
     115             :   }
     116             : 
     117             :   // If we get here, we have found a valid point; in particular we have
     118             :   // found the last valid point in the array.
     119             : 
     120             :   // Find the first valid point in the array.
     121             :   size_t first_valid_index = 0;
     122             :   for (size_t i = 0; i < y.size(); ++i) {
     123             :     if (y[i].has_value()) {
     124             :       first_valid_index = i;
     125             :       break;
     126             :     }
     127             :   }
     128             : 
     129             :   // Make a new set of points that includes only the points that
     130             :   // neighbor the boundary between valid and invalid points.
     131             :   std::vector<double> x_near_valid_point;
     132             :   std::vector<std::optional<double>> y_near_valid_point;
     133             : 
     134             :   if (first_valid_index == 0 and last_valid_index == y.size() - 1) {
     135             :     ERROR(
     136             :         "bracket_while_contracting: found a case where all points are valid,"
     137             :         "which should not happen under our assumptions.");
     138             :   }
     139             : 
     140             :   if (first_valid_index > 0) {
     141             :     // Check for a root between first_valid_index-1 and first_valid_index.
     142             :     const double x_test =
     143             :         x[first_valid_index - 1] +
     144             :         0.5 * (x[first_valid_index] - x[first_valid_index - 1]);
     145             :     const auto y_test = f(x_test);
     146             :     if (y_test.has_value() and
     147             :         y[first_valid_index].value() * y_test.value() <= 0.0) {
     148             :       // Bracketed!
     149             :       return {
     150             :           x_test, x[first_valid_index], y_test.value(),
     151             :           y[first_valid_index].value()};
     152             :     } else {
     153             :       x_near_valid_point.push_back(x[first_valid_index - 1]);
     154             :       y_near_valid_point.push_back(y[first_valid_index - 1]);
     155             :       x_near_valid_point.push_back(x_test);
     156             :       y_near_valid_point.push_back(y_test);
     157             :       x_near_valid_point.push_back(x[first_valid_index]);
     158             :       y_near_valid_point.push_back(y[first_valid_index]);
     159             :     }
     160             :   }
     161             :   if (last_valid_index < y.size() - 1) {
     162             :     // Check for a root between last_valid_index and last_valid_index+1.
     163             :     const double x_test = x[last_valid_index] +
     164             :                           0.5 * (x[last_valid_index + 1] - x[last_valid_index]);
     165             :     const auto y_test = f(x_test);
     166             :     if (y_test.has_value() and
     167             :         y[last_valid_index].value() * y_test.value() <= 0.0) {
     168             :       // Bracketed!
     169             :       return {x[last_valid_index], x_test, y[last_valid_index].value(),
     170             :               y_test.value()};
     171             :     } else {
     172             :       if (first_valid_index != last_valid_index or first_valid_index == 0) {
     173             :         x_near_valid_point.push_back(x[last_valid_index]);
     174             :         y_near_valid_point.push_back(y[last_valid_index]);
     175             :       }  // else we already pushed back last_valid_index (==first_valid_index).
     176             :       x_near_valid_point.push_back(x_test);
     177             :       y_near_valid_point.push_back(y_test);
     178             :       x_near_valid_point.push_back(x[last_valid_index + 1]);
     179             :       y_near_valid_point.push_back(y[last_valid_index + 1]);
     180             :     }
     181             :   }
     182             : 
     183             :   // We have one or more valid points but we didn't find a bracket.
     184             :   // That is, we have something like "X X o o X X" or "X X o o" or "o o X X".
     185             :   // So recurse, zooming in to the boundary (either one boundary or two
     186             :   // boundaries) between valid and invalid points.
     187             :   // Note that "o o o o" is prohibited by our assumptions, and checked for
     188             :   // above just in case it occurs by mistake.
     189             :   return bracket_by_contracting(x_near_valid_point, y_near_valid_point, f,
     190             :                                 level + 1);
     191             : }
     192             : }  // namespace bracketing_detail
     193             : 
     194             : /*!
     195             :  * \ingroup NumericalAlgorithmsGroup
     196             :  * \brief Brackets the root of the function `f`, assuming a single
     197             :  * root in a given interval \f$f[x_\mathrm{lo},x_\mathrm{up}]\f$
     198             :  * and assuming that `f` is defined only in an unknown smaller
     199             :  * interval \f$f[x_a,x_b]\f$ where
     200             :  * \f$x_\mathrm{lo} \leq x_a \leq x_b \leq x_\mathrm{hi}\f$.
     201             :  *
     202             :  * `f` is a unary invokable that takes a `double` which is the current value at
     203             :  * which to evaluate `f`.  `f` returns a `std::optional<double>` which
     204             :  * evaluates to false if the function is undefined at the supplied point.
     205             :  *
     206             :  * Assumes that there is only one root in the interval.
     207             :  *
     208             :  * Assumes that if \f$f(x_1)\f$ and \f$f(x_2)\f$ are both defined for
     209             :  * some \f$(x_1,x_2)\f$, then \f$f(x)\f$ is defined for all \f$x\f$
     210             :  * between \f$x_1\f$ and \f$x_2\f$.
     211             :  *
     212             :  * On input, assumes that the root lies in the interval
     213             :  * [`lower_bound`,`upper_bound`].  Optionally takes a `guess` for the
     214             :  * location of the root.  If `guess` is supplied, then evaluates the
     215             :  * function first at `guess` and `upper_bound` before trying
     216             :  * `lower_bound`: this means that it would be optimal if `guess`
     217             :  * underestimates the actual root and if `upper_bound` was less likely
     218             :  * to be undefined than `lower_bound`.
     219             :  *
     220             :  * On return, `lower_bound` and `upper_bound` are replaced with values that
     221             :  * bracket the root and for which the function is defined, and
     222             :  * `f_at_lower_bound` and `f_at_upper_bound` are replaced with
     223             :  * `f` evaluated at those bracketing points.
     224             :  *
     225             :  * `bracket_possibly_undefined_function_in_interval` throws an error if
     226             :  *  all points are valid but of the same sign (because that would indicate
     227             :  *  multiple roots but we assume only one root), if no root exists, or
     228             :  *  if the range of a sign change is sufficently small relative to the
     229             :  *  given interval that the number of iterations to find the root is exceeded.
     230             :  *
     231             :  */
     232             : template <typename Functor>
     233           1 : void bracket_possibly_undefined_function_in_interval(
     234             :     const gsl::not_null<double*> lower_bound,
     235             :     const gsl::not_null<double*> upper_bound,
     236             :     const gsl::not_null<double*> f_at_lower_bound,
     237             :     const gsl::not_null<double*> f_at_upper_bound, const Functor& f,
     238             :     const double guess) {
     239             :   // Initial values of x1,x2,y1,y2.  Use `guess` and `upper_bound`,
     240             :   // because in typical usage `guess` underestimates the actual
     241             :   // root, and `lower_bound` is more likely than `upper_bound` to be
     242             :   // invalid.
     243             :   double x1 = guess;
     244             :   double x2 = *upper_bound;
     245             :   auto y1 = f(x1);
     246             :   auto y2 = f(x2);
     247             :   const bool y1_defined = y1.has_value();
     248             :   const bool y2_defined = y2.has_value();
     249             :   if (not(y1_defined and y2_defined and y1.value() * y2.value() <= 0.0)) {
     250             :     // Root is not bracketed.
     251             :     // Before moving to the general algorithm, try the remaining
     252             :     // input point that was supplied.
     253             :     const double x3 = *lower_bound;
     254             :     const auto y3 = f(x3);
     255             :     const bool y3_defined = y3.has_value();
     256             :     if (y1_defined and y3_defined and y1.value() * y3.value() <= 0.0) {
     257             :       // Bracketed! Throw out x2,y2.  Rename variables to keep x1 < x2.
     258             :       x2 = x1;
     259             :       y2 = y1;
     260             :       x1 = x3;
     261             :       y1 = y3;
     262             :     } else {
     263             :       // Our simple checks didn't work, so call the more general method.
     264             :       // There are 8 cases:
     265             :       //
     266             :       // y3 y1 y2
     267             :       // --------
     268             :       // X  X  X
     269             :       // o  X  X
     270             :       // X  o  X
     271             :       // o  o  X
     272             :       // X  o  o
     273             :       // o  o  o
     274             :       // X  X  o
     275             :       // o  X  o
     276             :       //
     277             :       // where X means an invalid point, o means a valid point.
     278             :       // All valid points have the same sign, or we would have found a
     279             :       // bracket already.
     280             :       //
     281             :       // Before calling the general case, error on "o o o" and "o X o".
     282             :       // Both of these are prohibited by our assumptions (we
     283             :       // assume the root is in the interval so no "o o o", and we
     284             :       // assume that all invalid points are at the end of interval, so no
     285             :       // "o X o").
     286             :       if (y2_defined and y3_defined) {
     287             :         ERROR(
     288             :             "bracket_possibly_undefined_function_in_interval: found "
     289             :             "case that should not happen under our assumptions.");
     290             :       }
     291             :       try {
     292             :         std::tie(x1, x2, y1, y2) = bracketing_detail::bracket_by_contracting(
     293             :             {{x3, x1, x2}}, {{y3, y1, y2}}, f);
     294             :       } catch (std::runtime_error& e) {
     295             :         std::stringstream ss;
     296             :         ss << "bracket_by_contracting: Cannot bracket root between "
     297             :            << *lower_bound << " and " << *upper_bound
     298             :            << ". Internal error message is: '" << e.what() << "'";
     299             :         throw std::runtime_error(ss.str());
     300             :       }
     301             :     }
     302             :   }
     303             :   *f_at_lower_bound = y1.value();
     304             :   *f_at_upper_bound = y2.value();
     305             :   *lower_bound = x1;
     306             :   *upper_bound = x2;
     307             : }
     308             : 
     309             : /*!
     310             :  * \ingroup NumericalAlgorithmsGroup
     311             :  * \brief Brackets the single root of the
     312             :  * function `f` for each element in a `DataVector`, assuming the root
     313             :  * lies in the given interval and that `f` may be undefined at some
     314             :  * points in the interval.
     315             :  *
     316             :  * `f` is a binary invokable that takes a `double` and a `size_t` as
     317             :  * arguments.  The `double` is the current value at which to evaluate
     318             :  * `f`, and the `size_t` is the index into the `DataVector`s.  `f`
     319             :  * returns a `std::optional<double>` which evaluates to false if the
     320             :  * function is undefined at the supplied point.
     321             :  *
     322             :  * Assumes that there is only one root in the interval.
     323             :  *
     324             :  * Assumes that if \f$f(x_1)\f$ and \f$f(x_2)\f$ are both defined for
     325             :  * some \f$(x_1,x_2)\f$, then \f$f(x)\f$ is defined for all \f$x\f$
     326             :  * between \f$x_1\f$ and \f$x_2\f$.
     327             :  *
     328             :  * On input, assumes that the root lies in the interval
     329             :  * [`lower_bound`,`upper_bound`].  Optionally takes a `guess` for the
     330             :  * location of the root.
     331             :  *
     332             :  * On return, `lower_bound` and `upper_bound` are replaced with values that
     333             :  * bracket the root and for which the function is defined, and
     334             :  * `f_at_lower_bound` and `f_at_upper_bound` are replaced with
     335             :  * `f` evaluated at those bracketing points.
     336             :  *
     337             :  */
     338             : template <typename Functor>
     339           1 : void bracket_possibly_undefined_function_in_interval(
     340             :     const gsl::not_null<DataVector*> lower_bound,
     341             :     const gsl::not_null<DataVector*> upper_bound,
     342             :     const gsl::not_null<DataVector*> f_at_lower_bound,
     343             :     const gsl::not_null<DataVector*> f_at_upper_bound, const Functor& f,
     344             :     const DataVector& guess) {
     345             :   for (size_t s = 0; s < lower_bound->size(); ++s) {
     346             :     bracket_possibly_undefined_function_in_interval(
     347             :         &((*lower_bound)[s]), &((*upper_bound)[s]), &((*f_at_lower_bound)[s]),
     348             :         &((*f_at_upper_bound)[s]), [&f, &s](const double x) { return f(x, s); },
     349             :         guess[s]);
     350             :   }
     351             : }
     352             : 
     353             : /*
     354             :  * Version of `bracket_possibly_undefined_function_in_interval`
     355             :  * without a supplied initial guess; uses the mean of `lower_bound` and
     356             :  * `upper_bound` as the guess.
     357             :  */
     358             : template <typename Functor>
     359           0 : void bracket_possibly_undefined_function_in_interval(
     360             :     const gsl::not_null<double*> lower_bound,
     361             :     const gsl::not_null<double*> upper_bound,
     362             :     const gsl::not_null<double*> f_at_lower_bound,
     363             :     const gsl::not_null<double*> f_at_upper_bound, const Functor& f) {
     364             :   bracket_possibly_undefined_function_in_interval(
     365             :       lower_bound, upper_bound, f_at_lower_bound, f_at_upper_bound, f,
     366             :       *lower_bound + 0.5 * (*upper_bound - *lower_bound));
     367             : }
     368             : 
     369             : /*
     370             :  * Version of `bracket_possibly_undefined_function_in_interval`
     371             :  * without a supplied initial guess; uses the mean of `lower_bound` and
     372             :  * `upper_bound` as the guess.
     373             :  */
     374             : template <typename Functor>
     375           0 : void bracket_possibly_undefined_function_in_interval(
     376             :     const gsl::not_null<DataVector*> lower_bound,
     377             :     const gsl::not_null<DataVector*> upper_bound,
     378             :     const gsl::not_null<DataVector*> f_at_lower_bound,
     379             :     const gsl::not_null<DataVector*> f_at_upper_bound, const Functor& f) {
     380             :   bracket_possibly_undefined_function_in_interval(
     381             :       lower_bound, upper_bound, f_at_lower_bound, f_at_upper_bound, f,
     382             :       *lower_bound + 0.5 * (*upper_bound - *lower_bound));
     383             : }
     384             : }  // namespace RootFinder

Generated by: LCOV version 1.14