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: 46eb77706670
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 38b0917cd2e9
Choose a head ref

Commits on Oct 9, 2020

  1. Adjust tolerances.

    pleroy committed Oct 9, 2020
    Copy the full SHA
    0c57138 View commit details
  2. Fix one test.

    pleroy committed Oct 9, 2020
    Copy the full SHA
    d258b34 View commit details

Commits on Oct 11, 2020

  1. Remove an include.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    0f602a1 View commit details
  2. Lots of traces.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    45523e7 View commit details
  3. Eliminate some tests.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    de40d85 View commit details
  4. Copy the full SHA
    64d2328 View commit details
  5. Better name for a new test.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    6cb519a View commit details
  6. Copy the full SHA
    52ae001 View commit details
  7. Copy the full SHA
    ad7e8c9 View commit details
  8. Quadrature cleanup.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    b8583b7 View commit details
  9. Clean up frequency analysis.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    e4fa7a7 View commit details
  10. Copy the full SHA
    01917ef View commit details
  11. A test that tests something.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    d710915 View commit details
  12. Running the physics test.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    7e5813c View commit details
  13. Double tolerance.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    daa6f9f View commit details
  14. The right parameters.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    61cc456 View commit details
  15. Cleanup.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    b782525 View commit details
  16. Proper testing.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    e807d31 View commit details
  17. Test tolerances.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    1ddafc0 View commit details
  18. Copy the full SHA
    da07284 View commit details
  19. Cleanup quadrature.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    7bcca3a View commit details
  20. More cleanup.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    8ba4a3f View commit details
  21. Revert project changes.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    b8da338 View commit details
  22. Tolerances and lint.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    994b0cb View commit details
  23. Tolerances.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    e61a7ee View commit details
  24. Fix Sin tests.

    pleroy committed Oct 11, 2020
    Copy the full SHA
    8771ec9 View commit details

Commits on Oct 12, 2020

  1. Tolerances.

    pleroy committed Oct 12, 2020
    Copy the full SHA
    9287670 View commit details
  2. Merge pull request #2759 from pleroy/Tests

    Improved exit condition for the Clenshaw-Curtis quadrature
    pleroy authored Oct 12, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    38b0917 View commit details
1 change: 0 additions & 1 deletion numerics/finite_difference_test.cpp
Original file line number Diff line number Diff line change
@@ -6,7 +6,6 @@
#include "base/file.hpp"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "mathematica/mathematica.hpp"
#include "numerics/cbrt.hpp"
#include "testing_utilities/almost_equals.hpp"
#include "testing_utilities/is_near.hpp"
29 changes: 13 additions & 16 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -55,6 +55,7 @@ using testing_utilities::IsNear;
using testing_utilities::RelativeErrorFrom;
using testing_utilities::operator""_⑴;
using ::testing::AllOf;
using ::testing::Ge;
using ::testing::Gt;
using ::testing::Lt;
namespace si = quantities::si;
@@ -211,7 +212,6 @@ TEST_F(FrequencyAnalysisTest, PreciseModeVector) {
EXPECT_THAT(precise_mode, RelativeErrorFrom(ω, IsNear(4.2e-11_⑴)));
}

#if 0
TEST_F(FrequencyAnalysisTest, PoissonSeriesScalarProjection) {
AngularFrequency const ω = 666.543 * π * Radian / Second;
std::mt19937_64 random(42);
@@ -234,7 +234,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesScalarProjection) {
t_min, t_max);
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(projection4(t0_ + i * Radian / ω),
AlmostEquals(series(t0_ + i * Radian / ω), 0, 256));
AlmostEquals(series(t0_ + i * Radian / ω), 0, 768));
}

// Projection on a 5th degree basis is also accurate.
@@ -245,7 +245,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesScalarProjection) {
t_min, t_max);
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(projection5(t0_ + i * Radian / ω),
AlmostEquals(series(t0_ + i * Radian / ω), 0, 256));
AlmostEquals(series(t0_ + i * Radian / ω), 0, 768));
}

// Projection on a 3rd degree basis introduces significant errors.
@@ -312,7 +312,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesVectorProjection) {
t_min, t_max);
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(projection4(t0_ + i * Radian / ω),
AlmostEquals(series(t0_ + i * Radian / ω), 0, 384));
AlmostEquals(series(t0_ + i * Radian / ω), 0, 640));
}

// Projection on a 5th degree basis is also accurate.
@@ -323,7 +323,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesVectorProjection) {
t_min, t_max);
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(projection5(t0_ + i * Radian / ω),
AlmostEquals(series(t0_ + i * Radian / ω), 0, 384));
AlmostEquals(series(t0_ + i * Radian / ω), 0, 640));
}

// Projection on a 3rd degree basis introduces significant errors.
@@ -371,8 +371,8 @@ TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {
Instant const t_min = piecewise_series.t_min();
Instant const t_max = piecewise_series.t_max();

// Projection on a 4th degree basis. The approximation is only as good as the
// Gauss-Legendre integration.
// Projection on a 4th degree basis. The approximation is reasonably
// accurate.
auto const projection4 =
Projection<4>(piecewise_series,
ω,
@@ -382,12 +382,10 @@ TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {
EXPECT_THAT(
projection4(t_min + i * (t_max - t_min) / 100),
RelativeErrorFrom(series(t0_ + i * (t_max - t_min) / 100),
AllOf(Gt(2.1e-10), Lt(9.0e-4))));
AllOf(Gt(1.4e-9), Lt(9.9e-5))));
}
}

// TODO(phl): This test produces NaNs on casanova. Revisit once we have better
// integration.
TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionNoSecular) {
std::mt19937_64 random(42);
std::uniform_real_distribution<> frequency_distribution(2000.0, 3000.0);
@@ -428,7 +426,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionNoSecular) {
? AllOf(Gt(6.7e-2 * Metre), Lt(7.9 * Metre))
: ω_index == 2
? AllOf(Gt(1.1e-4 * Metre), Lt(9.7e-1 * Metre))
: AllOf(Gt(2.8e-11 * Metre), Lt(1.2e-6 * Metre)))
: AllOf(Gt(1.4e-16 * Metre), Lt(2.8e-13 * Metre)))
<< ω_index;
}
if (ω_index == ωs.size()) {
@@ -449,7 +447,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionNoSecular) {
EXPECT_THAT(
projection4(t_min + i * (t_max - t_min) / 100),
RelativeErrorFrom(series.value()(t_min + i * (t_max - t_min) / 100),
AllOf(Gt(5.9e-12), Lt(1.9e-6))));
AllOf(Ge(0), Lt(1.4e-11))));
}
}

@@ -492,8 +490,8 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionSecular) {
? AllOf(Gt(3.3e-2 * Metre), Lt(3.6 * Metre))
: ω_index == 3
? AllOf(Gt(7.5e-3 * Metre), Lt(5.4 * Metre))
: AllOf(Gt(3.7e-15 * Metre),
Lt(4.4e-11 * Metre)))
: AllOf(Gt(1.0e-16 * Metre),
Lt(6.0e-14 * Metre)))
<< ω_index;
}
if (ω_index == ωs.size()) {
@@ -514,10 +512,9 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionSecular) {
EXPECT_THAT(
projection4(t_min + i * (t_max - t_min) / 100),
RelativeErrorFrom(series(t_min + i * (t_max - t_min) / 100),
AllOf(Gt(1.9e-16), Lt(3.5e-12))));
AllOf(Ge(0), Lt(5.9e-15))));
}
}
#endif

} // namespace frequency_analysis
} // namespace numerics
9 changes: 9 additions & 0 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -189,6 +189,15 @@ class PoissonSeries {
template<typename V, int d, template<typename, typename, int> class E>
friend std::ostream& operator<<(std::ostream& out,
PoissonSeries<V, d, E> const& series);
template<typename L, typename R,
int l, int r, int w,
template<typename, typename, int> class E>
friend typename Hilbert<L, R>::InnerProductType InnerProduct(
PoissonSeries<L, l, E> const& left,
PoissonSeries<R, r, E> const& right,
PoissonSeries<double, w, E> const& weight,
Instant const& t_min,
Instant const& t_max);
template<typename V, int d,
template<typename, typename, int> class E,
typename O>
46 changes: 44 additions & 2 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -31,6 +31,12 @@ using quantities::Variation;
using quantities::si::Radian;
namespace si = quantities::si;

// These parameters have been tuned for approximation of the Moon over 3 months
// with 10 periods.
constexpr int clenshaw_curtis_min_points_overall = 33;
constexpr int clenshaw_curtis_point_per_period = 4;
constexpr double clenshaw_curtis_relative_error = 0x1p-32;

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
typename PoissonSeries<Primitive<Value, Time>, degree_ + 1, Evaluator>::
@@ -345,10 +351,27 @@ PoissonSeries<Value, degree_, Evaluator>::Norm(
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) const {
AngularFrequency const max_ω =
(periodic_.empty() ? AngularFrequency{}
: 2 * periodic_.back().first) +
(weight.periodic_.empty() ? AngularFrequency{}
: weight.periodic_.back().first);
std::optional<int> max_points =
max_ω == AngularFrequency()
? std::optional<int>{}
: std::max(
clenshaw_curtis_min_points_overall,
static_cast<int>(clenshaw_curtis_point_per_period *
(t_max - t_min) * max_ω / (2 * π * Radian)));

auto integrand = [this, &weight](Instant const& t) {
return Hilbert<Value>::InnerProduct((*this)(t), (*this)(t)) * weight(t);
};
return Sqrt(quadrature::AutomaticClenshawCurtis(integrand, t_min, t_max) /
return Sqrt(quadrature::AutomaticClenshawCurtis(
integrand,
t_min, t_max,
/*max_relative_error=*/clenshaw_curtis_relative_error,
/*max_points=*/max_points) /
(t_max - t_min));
}

@@ -646,10 +669,29 @@ typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
AngularFrequency const max_ω =
(left.periodic_.empty() ? AngularFrequency{}
: left.periodic_.back().first) +
(right.periodic_.empty() ? AngularFrequency{}
: right.periodic_.back().first) +
(weight.periodic_.empty() ? AngularFrequency{}
: weight.periodic_.back().first);
std::optional<int> max_points =
max_ω == AngularFrequency()
? std::optional<int>{}
: std::max(
clenshaw_curtis_min_points_overall,
static_cast<int>(clenshaw_curtis_point_per_period *
(t_max - t_min) * max_ω / (2 * π * Radian)));

auto integrand = [&left, &right, &weight](Instant const& t) {
return Hilbert<LValue, RValue>::InnerProduct(left(t), right(t)) * weight(t);
};
return quadrature::AutomaticClenshawCurtis(integrand, t_min, t_max) /
return quadrature::AutomaticClenshawCurtis(
integrand,
t_min, t_max,
/*max_relative_error=*/clenshaw_curtis_relative_error,
/*max_points=*/max_points) /
(t_max - t_min);
}

54 changes: 53 additions & 1 deletion numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -36,6 +36,7 @@ using geometry::Velocity;
using quantities::Acceleration;
using quantities::AngularFrequency;
using quantities::Cos;
using quantities::Length;
using quantities::Pow;
using quantities::Sin;
using quantities::Sqrt;
@@ -312,7 +313,58 @@ TEST_F(PoissonSeriesTest, InnerProduct) {
apodization::Hann<HornerEvaluator>(t_min, t_max),
t_min,
t_max),
AlmostEquals(-381.25522770148542400, 7));
AlmostEquals(-381.25522770148542400, 4, 7));
}

TEST_F(PoissonSeriesTest, PoorlyConditionedInnerProduct) {
using Degree4 = PoissonSeries<Length, 4, HornerEvaluator>;
using Degree5 = PoissonSeries<Length, 5, HornerEvaluator>;
Instant const t_min = t0_;
Instant const t_max = t0_ + 4.77553415434249021e-02 * Second;
AngularFrequency const ω = 2.09400659210170170e+03 * Radian / Second;
Degree4 const f(Degree4::Polynomial({}, t0_),
{{ω,
{/*sin=*/Degree4::Polynomial(
{+5.10311065909077932e+00 * Metre,
+2.78062787709394854e+00 * Metre / Second,
+5.04290401496053242e+00 * Metre / Pow<2>(Second),
-7.27454632735125806e+00 * Metre / Pow<3>(Second),
+8.06537932856756612e+00 * Metre / Pow<4>(Second)},
t0_),
/*cos=*/Degree4::Polynomial(
{-8.11863376474325804e+00 * Metre,
+1.49140608216528037e+00 * Metre / Second,
-2.54224601087630298e+00 * Metre / Pow<2>(Second),
-4.52251796525658367e+00 * Metre / Pow<3>(Second),
-2.19458237171412751e+00 * Metre / Pow<4>(Second)},
t0_)}}});
Degree5 const q(Degree5::Polynomial({}, t0_),
{{ω,
{/*sin=*/Degree5::Polynomial(
{-4.41249783881549433e+01 * Metre,
+1.50208859053174347e+04 * Metre / Second,
-1.70674564621978020e+06 * Metre / Pow<2>(Second),
+8.52015772027946562e+07 * Metre / Pow<3>(Second),
-1.92799129073151779e+09 * Metre / Pow<4>(Second),
+1.61514557548221931e+10 * Metre / Pow<5>(Second)},
t0_),
/*cos=*/Degree5::Polynomial(
{-1.00752842659088765e-01 * Metre,
+2.25402995957193006e+01 * Metre / Second,
-1.66819064858902379e+03 * Metre / Pow<2>(Second),
+4.98682536071893774e+04 * Metre / Pow<3>(Second),
-5.18229522289936838e+05 * Metre / Pow<4>(Second),
0 * Metre / Pow<5>(Second)},
t0_)}}});

// The integral is very small compared to the functions, so we end up in the
// numerical noise, and adding more points would not help much.
auto const product = InnerProduct(
f, q, apodization::Hann<HornerEvaluator>(t_min, t_max), t_min, t_max);
// Exact result obtained using Mathematica.
EXPECT_THAT(product,
RelativeErrorFrom(-4.848079980325297e-13 * Metre * Metre,
IsNear(0.05_⑴)));
}

TEST_F(PoissonSeriesTest, Output) {
12 changes: 7 additions & 5 deletions numerics/quadrature.hpp
Original file line number Diff line number Diff line change
@@ -1,17 +1,16 @@

#pragma once

#include <optional>
#include <type_traits>

#include "geometry/hilbert.hpp"
#include "quantities/named_quantities.hpp"

namespace principia {
namespace numerics {
namespace quadrature {
namespace internal_quadrature {

using geometry::Hilbert;
using quantities::Primitive;

template<int points, typename Argument, typename Function>
@@ -21,15 +20,18 @@ Primitive<std::invoke_result_t<Function, Argument>, Argument> GaussLegendre(
Argument const& upper_bound);

// Computes a Clenshaw-Curtis quadrature on 2ᵖ + 1 points for successive p until
// the tolerance is satisfied.
// TODO(egg): The semantics of the |relative_tolerance| are unclear.
// the tolerance is satisfied. The client controls the accuracy of the result
// using |max_relative_error| (returns when the relative error on the result is
// estimated to be less that the specified value) and |max_points| (returns when
// the number of points would exceed the specified value).
template<int initial_points = 2, typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument>
AutomaticClenshawCurtis(
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
double relative_tolerance = 0x1p-16);
std::optional<double> max_relative_error,
std::optional<int> max_points);

// |points| must be of the form 2ᵖ + 1 for some p ∈ ℕ. Returns the
// Clenshaw-Curtis quadrature of f with the given number of points.
31 changes: 19 additions & 12 deletions numerics/quadrature_body.hpp
Original file line number Diff line number Diff line change
@@ -6,6 +6,7 @@
#include <memory>
#include <vector>

#include "geometry/hilbert.hpp"
#include "numerics/fast_fourier_transform.hpp"
#include "numerics/gauss_legendre_weights.mathematica.h"
#include "numerics/legendre_roots.mathematica.h"
@@ -17,6 +18,7 @@ namespace quadrature {
namespace internal_quadrature {

using base::FloorLog2;
using geometry::Hilbert;
using quantities::Angle;
using quantities::Cos;
using quantities::Difference;
@@ -59,26 +61,30 @@ AutomaticClenshawCurtisImplementation(
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
double const relative_tolerance,
std::optional<double> const max_relative_error,
std::optional<int> const max_points,
Primitive<std::invoke_result_t<Function, Argument>, Argument> const
previous_estimate) {
using Result = Primitive<std::invoke_result_t<Function, Argument>, Argument>;
Result const estimate = ClenshawCurtis<points>(f, lower_bound, upper_bound);

// This is the naïve estimate mentioned in [Gen72b], p. 339.
double const relative_error_estimate =
Hilbert<Result>::Norm(previous_estimate - estimate) /
Hilbert<Result>::Norm(estimate);
// We look for an estimated relative error smaller than
// |relative_tolerance * points|: since the integral is computed from |points|
// evaluations of |f|, it will necessarily carry a relative error proportional
// to |points|, so it makes no sense to look for convergence beyond that.
if (relative_error_estimate > relative_tolerance * points) {
auto const absolute_error_estimate =
Hilbert<Result>::Norm(previous_estimate - estimate);

if ((!max_relative_error.has_value() ||
absolute_error_estimate >
max_relative_error.value() * Hilbert<Result>::Norm(estimate)) &&
(!max_points.has_value() || points <= max_points.value())) {
if constexpr (points > 1 << 24) {
LOG(FATAL) << "Too many refinements while integrating from "
<< lower_bound << " to " << upper_bound;
} else {
return AutomaticClenshawCurtisImplementation<points + points - 1>(
f, lower_bound, upper_bound, relative_tolerance, estimate);
f,
lower_bound, upper_bound,
max_relative_error, max_points,
estimate);
}
}
return estimate;
@@ -90,14 +96,15 @@ AutomaticClenshawCurtis(
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
double const relative_tolerance) {
std::optional<double> const max_relative_error,
std::optional<int> const max_points) {
using Result = Primitive<std::invoke_result_t<Function, Argument>, Argument>;
// TODO(egg): factor the evaluations of f.
Result const estimate =
ClenshawCurtis<initial_points>(f, lower_bound, upper_bound);
return AutomaticClenshawCurtisImplementation<
initial_points + initial_points - 1>(
f, lower_bound, upper_bound, relative_tolerance, estimate);
f, lower_bound, upper_bound, max_relative_error, max_points, estimate);
}

template<int points, typename Argument, typename Function>
27 changes: 16 additions & 11 deletions numerics/quadrature_test.cpp
Original file line number Diff line number Diff line change
@@ -22,6 +22,7 @@ using testing_utilities::AlmostEquals;
using testing_utilities::IsNear;
using testing_utilities::RelativeErrorFrom;
using testing_utilities::operator""_⑴;
using ::testing::AnyOf;
using ::testing::Eq;

class QuadratureTest : public ::testing::Test {};
@@ -65,11 +66,13 @@ TEST_F(QuadratureTest, Sin) {
AlmostEquals(ʃf, 3, 4));

evaluations = 0;
EXPECT_THAT(AutomaticClenshawCurtis(f,
-2.0 * Radian,
5.0 * Radian,
std::numeric_limits<double>::epsilon()),
AlmostEquals(ʃf, 4));
EXPECT_THAT(AutomaticClenshawCurtis(
f,
-2.0 * Radian,
5.0 * Radian,
/*max_relative_error=*/std::numeric_limits<double>::epsilon(),
/*max_points=*/std::nullopt),
AlmostEquals(ʃf, 3, 4));
EXPECT_THAT(evaluations, Eq(134));
}

@@ -117,12 +120,14 @@ TEST_F(QuadratureTest, Sin10) {
return Sin(10 * x);
};
auto const ʃf = (Cos(20 * Radian) - Cos(50 * Radian)) / 10 * Radian;
EXPECT_THAT(AutomaticClenshawCurtis(f,
-2.0 * Radian,
5.0 * Radian,
std::numeric_limits<double>::epsilon()),
AlmostEquals(ʃf, 196));
EXPECT_THAT(evaluations, Eq(263));
EXPECT_THAT(AutomaticClenshawCurtis(
f,
-2.0 * Radian,
5.0 * Radian,
/*max_relative_error=*/std::numeric_limits<double>::epsilon(),
/*max_points=*/std::nullopt),
AlmostEquals(ʃf, 9, 277));
EXPECT_THAT(evaluations, AnyOf(Eq(131088), Eq(524306)));
}

} // namespace quadrature