Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 27435934e0dd
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 97ccf5692d7c
Choose a head ref
  • 6 commits
  • 3 files changed
  • 1 contributor

Commits on Sep 23, 2020

  1. no defaults

    eggrobin committed Sep 23, 2020
    Copy the full SHA
    224e9a1 View commit details

Commits on Sep 24, 2020

  1. localmin

    eggrobin committed Sep 24, 2020
    Copy the full SHA
    a1edcbb View commit details
  2. after pleroy’s review

    eggrobin committed Sep 24, 2020
    Copy the full SHA
    3a63f7c View commit details
  3. Copy the full SHA
    4437ad8 View commit details
  4. lint

    eggrobin committed Sep 24, 2020
    Copy the full SHA
    7385605 View commit details
  5. Merge pull request #2732 from eggrobin/good-minimizers

    they're good minimization algorithms Brent
    eggrobin authored Sep 24, 2020
    Copy the full SHA
    97ccf56 View commit details
Showing with 302 additions and 22 deletions.
  1. +19 −8 numerics/root_finders.hpp
  2. +134 −7 numerics/root_finders_body.hpp
  3. +149 −7 numerics/root_finders_test.cpp
27 changes: 19 additions & 8 deletions numerics/root_finders.hpp
Original file line number Diff line number Diff line change
@@ -19,7 +19,6 @@ using quantities::Derivative;
// function agreeing with |f| on the values of |Argument|.
// If |f(lower_bound)| and |f(upper_bound)| are both nonzero, they must be of
// opposite signs.
// TODO(phl): Use Brent's algorithm.
template<typename Argument, typename Function>
Argument Bisect(Function f,
Argument const& lower_bound,
@@ -32,17 +31,29 @@ Argument Brent(Function f,
Argument const& lower_bound,
Argument const& upper_bound);

// Performs a golden-section search to find a minimum of |f| between
// Performs a golden-section search to find a local extremum of |f| between
// |lower_bound| and |upper_bound|.
// TODO(phl): Use Brent's algorithm.
template<typename Argument,
typename Function,
typename Compare = std::less<
decltype(std::declval<Function>()(std::declval<Argument>()))>>
// The function searches for a minimum if compare is <, and a maximum if compare
// is >. Arbitrary order relations are allowed; in general, this function
// searches for a value x such that compare(y, f(x)) is false for all y in some
// neighbourhood of x.
template<typename Argument, typename Function, typename Compare>
Argument GoldenSectionSearch(Function f,
Argument const& lower_bound,
Argument const& upper_bound,
Compare comp = Compare());
Compare compare);

// Performs Brent’s procedure |localmin| from [Bre73], chapter 5, with an
// absolute tolerance t=0.
// The function searches for a minimum if compare is <, and a maximum if compare
// is >. No values of Compare other than std::less<> and std::greater<> are
// allowed.
template<typename Argument, typename Function, typename Compare>
Argument Brent(Function f,
Argument const& lower_bound,
Argument const& upper_bound,
Compare compare,
double eps);

// Returns the solutions of the quadratic equation:
// a2 * (x - origin)^2 + a1 * (x - origin) + a0 == 0
141 changes: 134 additions & 7 deletions numerics/root_finders_body.hpp
Original file line number Diff line number Diff line change
@@ -4,8 +4,9 @@
#include "numerics/root_finders.hpp"

#include <algorithm>
#include <vector>
#include <functional>
#include <limits>
#include <vector>

#include "geometry/barycentre_calculator.hpp"
#include "geometry/sign.hpp"
@@ -21,6 +22,7 @@ using geometry::Barycentre;
using geometry::Sign;
using quantities::Abs;
using quantities::Difference;
using quantities::Product;
using quantities::Sqrt;
using quantities::Square;

@@ -71,7 +73,6 @@ Argument Brent(Function f,
using Value = decltype(f(lower_bound));
Value const zero{};


// We do not use |std::numeric_limits<double>::epsilon()|, because it is 2ϵ in
// Brent’s notation: Brent uses ϵ = β^(1-τ) / 2 for rounded arithmetic, see
// (2.9).
@@ -161,7 +162,7 @@ template<typename Argument, typename Function, typename Compare>
Argument GoldenSectionSearch(Function f,
Argument const& lower_bound,
Argument const& upper_bound,
Compare const comp) {
Compare const compare) {
static constexpr double lower_interior_ratio = 2 - φ;
static constexpr double upper_interior_ratio = φ - 1;
using Value = decltype(f(lower_bound));
@@ -183,9 +184,9 @@ Argument GoldenSectionSearch(Function f,
while (lower < lower_interior &&
lower_interior < upper_interior &&
upper_interior < upper) {
Value const f_lower_min = std::min(f_lower, f_lower_interior, comp);
Value const f_upper_min = std::min(f_upper_interior, f_upper, comp);
if (comp(f_lower_min, f_upper_min)) {
Value const f_lower_min = std::min(f_lower, f_lower_interior, compare);
Value const f_upper_min = std::min(f_upper_interior, f_upper, compare);
if (compare(f_lower_min, f_upper_min)) {
upper = upper_interior;
f_upper = f_upper_interior;
upper_interior = lower_interior;
@@ -203,10 +204,136 @@ Argument GoldenSectionSearch(Function f,
f_upper_interior = f(upper_interior);
}
}

return Barycentre<Argument, double>({lower, upper}, {1, 1});
}

// The implementation is translated from the ALGOL 60 in [Bre73], chapter 5,
// section 8.
template<typename Argument, typename Function, typename Compare>
Argument Brent(Function f,
Argument const& lower_bound,
Argument const& upper_bound,
Compare compare,
double eps) {
using Value = decltype(f(lower_bound));

static_assert(std::is_same_v<Compare, std::less<>> ||
std::is_same_v<Compare, std::greater<>>,
"Brent’s method relies on the consistency of the order whose "
"extremum is sought with the arithmetic operations. For "
"arbitrary order relations, use golden section search.");

// The code from [Bre73] looks for a minimum; for a maximum, we look for a
// minimum of the opposite.
auto const minimized_f = [&f](Argument const& x) {
if constexpr (std::is_same_v<Compare, std::greater<>>) {
return -f(x);
} else {
return f(x);
}
};
{
auto const& f = minimized_f;

// We do not use |std::numeric_limits<double>::epsilon()|, because it is 2ϵ
// in Brent’s notation: Brent uses ϵ = β^(1-τ) / 2 for rounded arithmetic,
// see [Bre73], chapter 4, section 2, (2.9).
constexpr double ϵ = ScaleB(0.5, 1 - std::numeric_limits<double>::digits);
// In order to ensure convergence, eps should be no smaller than 2ϵ, see
// [Bre73] chapter 5, section 5.
eps = std::max(eps, 2 * ϵ);

Argument a = lower_bound;
Argument b = upper_bound;
constexpr double c = 2 - φ;
Difference<Argument> d;
Argument u;
Argument v;
Argument w;
Argument x;
Value f_u;
Value f_v;
Value f_w;
Value f_x;

v = w = x = a + c * (b - a);
Difference<Argument> e{};
f_v = f_w = f_x = f(x);
for (;;) {
Argument const m = Barycentre<Argument, double>({a, b}, {1, 1});
Difference<Argument> const tol = eps * Abs(x - Argument{});
Difference<Argument> const t2 = 2 * tol;
// Check stopping criterion.
if (Abs(x - m) <= t2 - 0.5 * (b - a)) {
return x;
}
// p = q = r = 0;
Product<Square<Difference<Argument>>, Difference<Value>> p{};
Product<Difference<Argument>, Difference<Value>> q{};
if (Abs(e) > tol) {
// Fit parabola.
auto const r₁ = (x - w) * (f_x - f_v);
auto const r₂ = (x - v) * (f_x - f_w);
p = (x - v) * r₂ - (x - w) * r₁;
q = 2 * (r₂ - r₁);
if (Sign(q).is_positive()) {
p = -p;
} else {
q = -q;
}
}
// The second clause is incorrectly p < q * (a - x) in [Bre73] p.80, see
// the errata.
if (Abs(p) < Abs(0.5 * q * e) && p > q * (a - x) && p < q * (b - x)) {
e = d;
// A “parabolic interpolation” step.
d = p / q;
u = x + d;
// f must not be evaluated too close to a or b.
if (u - a < t2 || b - u < t2) {
d = x < m ? tol : -tol;
}
} else {
// A “golden section” step.
e = (x < m ? b : a) - x;
d = c * e;
}
// f must not be evaluated too close to x.
u = x + (Abs(d) > tol ? d : tol * Sign(d));
f_u = f(u);
// Update a, b, v, w, and x.
if (f_u <= f_x) {
if (u < x) {
b = x;
} else {
a = x;
}
v = w;
f_v = f_w;
w = x;
f_w = f_x;
x = u;
f_x = f_u;
} else {
if (u < x) {
a = u;
} else {
b = u;
}
if (f_u <= f_w || w == x) {
v = w;
f_v = f_w;
w = u;
f_w = f_u;
} else if (f_u <= f_v || v == x || v == w) {
v = u;
f_v = f_u;
}
}
}
}
}

template<typename Argument, typename Value>
BoundedArray<Argument, 2> SolveQuadraticEquation(
Argument const& origin,
Loading