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: cadff589b17e
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: d78a0ac6e028
Choose a head ref
  • 8 commits
  • 3 files changed
  • 2 contributors

Commits on Jul 17, 2020

  1. Implementation skeleton.

    pleroy committed Jul 17, 2020
    Copy the full SHA
    ecbbab4 View commit details
  2. Complete algorithm.

    pleroy committed Jul 17, 2020
    Copy the full SHA
    18b6365 View commit details
  3. Test passing.

    pleroy committed Jul 17, 2020
    Copy the full SHA
    27af8ee View commit details
  4. Maximum.

    pleroy committed Jul 17, 2020
    Copy the full SHA
    6a51791 View commit details
  5. Comment.

    pleroy committed Jul 17, 2020
    Copy the full SHA
    789ef18 View commit details
  6. After egg's review.

    pleroy committed Jul 17, 2020
    Copy the full SHA
    25b05e3 View commit details
  7. A new test.

    pleroy committed Jul 17, 2020
    Copy the full SHA
    3115f2a View commit details
  8. Merge pull request #2652 from pleroy/GoldenSectionSearch

    Golden section search
    pleroy authored Jul 17, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    d78a0ac View commit details
Showing with 109 additions and 0 deletions.
  1. +15 −0 numerics/root_finders.hpp
  2. +53 −0 numerics/root_finders_body.hpp
  3. +41 −0 numerics/root_finders_test.cpp
15 changes: 15 additions & 0 deletions numerics/root_finders.hpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@

#pragma once

#include <functional>
#include <vector>

#include "base/array.hpp"
@@ -18,11 +19,24 @@ 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,
Argument const& upper_bound);

// Performs a golden-section search to find a minimum 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>()))>>
Argument GoldenSectionSearch(Function f,
Argument const& lower_bound,
Argument const& upper_bound,
Compare comp = Compare());

// Returns the solutions of the quadratic equation:
// a2 * (x - origin)^2 + a1 * (x - origin) + a0 == 0
// The result may have 0, 1 or 2 values and is sorted.
@@ -36,6 +50,7 @@ BoundedArray<Argument, 2> SolveQuadraticEquation(
} // namespace internal_root_finders

using internal_root_finders::Bisect;
using internal_root_finders::GoldenSectionSearch;
using internal_root_finders::SolveQuadraticEquation;

} // namespace numerics
53 changes: 53 additions & 0 deletions numerics/root_finders_body.hpp
Original file line number Diff line number Diff line change
@@ -3,6 +3,7 @@

#include "root_finders.hpp"

#include <algorithm>
#include <vector>

#include "geometry/barycentre_calculator.hpp"
@@ -56,6 +57,58 @@ Argument Bisect(Function f,
}
}

// See https://en.wikipedia.org/wiki/Golden-section_search for a description of
// this algorithm.
template<typename Argument, typename Function, typename Compare>
Argument GoldenSectionSearch(Function f,
Argument const& lower_bound,
Argument const& upper_bound,
Compare const comp) {
static constexpr double lower_interior_ratio = 2 - φ;
static constexpr double upper_interior_ratio = φ - 1;
using Value = decltype(f(lower_bound));

Argument upper = upper_bound;
Value f_upper = f(upper);

Argument lower = lower_bound;
Value f_lower = f(lower);

Argument lower_interior = Barycentre<Argument, double>(
{lower, upper}, {upper_interior_ratio, lower_interior_ratio});
Value f_lower_interior = f(lower_interior);

Argument upper_interior = Barycentre<Argument, double>(
{lower, upper}, {lower_interior_ratio, upper_interior_ratio});
Value f_upper_interior = f(upper_interior);

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)) {
upper = upper_interior;
f_upper = f_upper_interior;
upper_interior = lower_interior;
f_upper_interior = f_lower_interior;
lower_interior = Barycentre<Argument, double>(
{lower, upper}, {upper_interior_ratio, lower_interior_ratio});
f_lower_interior = f(lower_interior);
} else {
lower = lower_interior;
f_lower = f_lower_interior;
lower_interior = upper_interior;
f_lower_interior = f_upper_interior;
upper_interior = Barycentre<Argument, double>(
{lower, upper}, {lower_interior_ratio, upper_interior_ratio});
f_upper_interior = f(upper_interior);
}
}

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

template<typename Argument, typename Value>
BoundedArray<Argument, 2> SolveQuadraticEquation(
Argument const& origin,
41 changes: 41 additions & 0 deletions numerics/root_finders_test.cpp
Original file line number Diff line number Diff line change
@@ -1,11 +1,13 @@

#include "numerics/root_finders.hpp"

#include <functional>
#include <vector>

#include "geometry/named_quantities.hpp"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "quantities/elementary_functions.hpp"
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/almost_equals.hpp"
@@ -16,9 +18,11 @@ using geometry::Instant;
using quantities::Acceleration;
using quantities::Length;
using quantities::Pow;
using quantities::Sin;
using quantities::Sqrt;
using quantities::Time;
using quantities::si::Metre;
using quantities::si::Radian;
using quantities::si::Second;
using testing_utilities::AlmostEquals;
using ::testing::AllOf;
@@ -53,6 +57,43 @@ TEST_F(RootFindersTest, SquareRoots) {
}
}

TEST_F(RootFindersTest, GoldenSectionSearch) {
Instant const t_0;
auto sin = [t_0](Instant const& t) {
return Sin((t - t_0) * Radian / Second);
};

// Minimum.
for (int l = 16; l <= 47; ++l) {
for (int u = 48; u <= 62; ++u) {
// The result is not overly precise because near its maximum a the
// function is:
// f(a) + fʺ(a) (x - a) / 2 + o((x - a)²)
// The second order term vanishes when x and a match on the first 26
// leading bits (roughly).
EXPECT_THAT(
GoldenSectionSearch(sin,
t_0 + l * 0.1 * Second,
t_0 + u * 0.1 * Second),
AlmostEquals(t_0 + 3 * π / 2 * Second, 11'863'280, 11'863'284));
}
}

// Maximum.
EXPECT_THAT(
(GoldenSectionSearch<Instant, decltype(sin), std::greater<double>>)(
sin,
t_0 + 1.5 * Second,
t_0 + 1.6 * Second),
AlmostEquals(t_0 + π / 2 * Second, 47453132));

// A big interval will yield a semi-random minimum.
EXPECT_THAT(GoldenSectionSearch(sin,
t_0 - 100 * Second,
t_0 + 666 * Second),
AlmostEquals(t_0 + 119 * π / 2 * Second, 370'728));
}

TEST_F(RootFindersTest, QuadraticEquations) {
// Golden ratio.
auto const s1 = SolveQuadraticEquation(0.0, -1.0, -1.0, 1.0);