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

Commits on Oct 6, 2020

  1. we are going to need a DCT

    eggrobin committed Oct 6, 2020
    Copy the full SHA
    3cea830 View commit details
  2. another tangent

    eggrobin committed Oct 6, 2020
    Copy the full SHA
    0b6851d View commit details
  3. merge

    eggrobin committed Oct 6, 2020
    Copy the full SHA
    2ebdbaa View commit details
  4. weigh the scalars

    eggrobin committed Oct 6, 2020
    Copy the full SHA
    4dda81b View commit details
  5. no time to lose

    eggrobin committed Oct 6, 2020
    Copy the full SHA
    a39c0c0 View commit details
  6. Copy the full SHA
    cb77bf0 View commit details

Commits on Oct 7, 2020

  1. it sort of works

    eggrobin committed Oct 7, 2020
    Copy the full SHA
    6750a77 View commit details
  2. Copy the full SHA
    96d0c99 View commit details
  3. Copy the full SHA
    7fdf891 View commit details
  4. Copy the full SHA
    3f06008 View commit details
  5. Copy the full SHA
    f7ea60d View commit details
  6. merge

    eggrobin committed Oct 7, 2020
    Copy the full SHA
    a6589c9 View commit details
  7. Revert "it still mostly works in double precision"

    This reverts commit 3f06008.
    eggrobin committed Oct 7, 2020
    Copy the full SHA
    991e4db View commit details
  8. Copy the full SHA
    d211fe4 View commit details
  9. Copy the full SHA
    7e46120 View commit details
  10. revert spurious changes

    eggrobin committed Oct 7, 2020
    Copy the full SHA
    9c83b5a View commit details
  11. remove traces

    eggrobin committed Oct 7, 2020
    Copy the full SHA
    a0ee0a6 View commit details

Commits on Oct 8, 2020

  1. A test

    eggrobin committed Oct 8, 2020
    Copy the full SHA
    b66a0be View commit details
  2. Copy the full SHA
    250a9cc View commit details
  3. cite and document

    eggrobin committed Oct 8, 2020
    Copy the full SHA
    5a94014 View commit details
  4. the tests are very broken

    eggrobin committed Oct 8, 2020
    Copy the full SHA
    bca7fc2 View commit details
  5. 7

    eggrobin committed Oct 8, 2020
    Copy the full SHA
    e430c33 View commit details
  6. iwyu

    eggrobin committed Oct 8, 2020
    Copy the full SHA
    b68cc00 View commit details
  7. after pleroy’s review

    eggrobin committed Oct 8, 2020
    Copy the full SHA
    5aa021a View commit details
  8. encoding

    eggrobin committed Oct 8, 2020
    Copy the full SHA
    3597ca0 View commit details

Commits on Oct 9, 2020

  1. Merge pull request #2756 from eggrobin/clenshaw-curtis

    Implement an automatic Clenshaw-Curtis quadrature and use it in the Poisson series inner product and norm
    eggrobin authored Oct 9, 2020
    Copy the full SHA
    46eb777 View commit details
35 changes: 35 additions & 0 deletions documentation/bibliography.bib
Original file line number Diff line number Diff line change
@@ -448,6 +448,41 @@ @article{Fukushima2011c
volume = {236},
}

@article{Gentleman1972a,
author = {Gentleman, W Morven},
date = {1972-05},
doi = {10.1145/355602.355603},
journaltitle = {Communications of the ACM},
number = {5},
pages = {353--355},
title = {Algorithm 424: Clenshaw-Curtis quadrature [D1]},
volume = {15},
}

@article{Gentleman1972b,
author = {Gentleman, W Morven},
date = {1972-05},
doi = {10.1145/355602.361310},
issn = {0001-0782},
journaltitle = {Communications of the ACM},
number = {5},
pages = {337--342},
title = {Implementing Clenshaw-Curtis quadrature, I methodology and experience},
volume = {15},
}

@article{Gentleman1972c,
author = {Gentleman, W Morven},
date = {1972-05},
doi = {10.1145/355602.361311},
issn = {0001-0782},
journaltitle = {Communications of the ACM},
number = {5},
pages = {343--346},
title = {Implementing Clenshaw-Curtis quadrature, II computing the cosine transformation},
volume = {15},
}

@article{HairerMcLachlanRazakarivony2008,
author = {Hairer, E. and McLachlan, R. I. and Razakarivony, A.},
language = {English},
Binary file modified documentation/bibliography.pdf
Binary file not shown.
6 changes: 5 additions & 1 deletion geometry/hilbert.hpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#pragma once
#pragma once

#include <type_traits>

@@ -48,6 +48,8 @@ struct Hilbert<T, T, std::enable_if_t<is_quantity_v<T>>>
using NormType = T;
static NormType Norm(T const& t);

// TODO(egg): Hilbert::Norm².

using NormalizedType = double;
};

@@ -78,6 +80,8 @@ struct Hilbert<T, T,
using NormType = decltype(std::declval<T>().Norm());
static NormType Norm(T const& t);

// TODO(egg): Hilbert::Norm².

using NormalizedType = Quotient<T, NormType>;
};

4 changes: 2 additions & 2 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -211,6 +211,7 @@ 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);
@@ -387,7 +388,6 @@ TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {

// TODO(phl): This test produces NaNs on casanova. Revisit once we have better
// integration.
#if 0
TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionNoSecular) {
std::mt19937_64 random(42);
std::uniform_real_distribution<> frequency_distribution(2000.0, 3000.0);
@@ -452,7 +452,6 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionNoSecular) {
AllOf(Gt(5.9e-12), Lt(1.9e-6))));
}
}
#endif

TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionSecular) {
std::mt19937_64 random(42);
@@ -518,6 +517,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionSecular) {
AllOf(Gt(1.9e-16), Lt(3.5e-12))));
}
}
#endif

} // namespace frequency_analysis
} // namespace numerics
27 changes: 16 additions & 11 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -345,8 +345,11 @@ PoissonSeries<Value, degree_, Evaluator>::Norm(
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) const {
auto const integrand = PointwiseInnerProduct(*this, *this) * weight;
return Sqrt(integrand.Integrate(t_min, t_max) / (t_max - t_min));
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) /
(t_max - t_min));
}

template<typename Value, int degree_,
@@ -637,15 +640,17 @@ std::ostream& operator<<(
template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
typename Hilbert<LValue, RValue>::InnerProductType
InnerProduct(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
auto const integrand = PointwiseInnerProduct(left, right) * weight;
auto const primitive = integrand.Primitive();
return (primitive(t_max) - primitive(t_min)) / (t_max - t_min);
typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
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) /
(t_max - t_min);
}

template<typename Value, int degree_,
2 changes: 1 addition & 1 deletion numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -312,7 +312,7 @@ TEST_F(PoissonSeriesTest, InnerProduct) {
apodization::Hann<HornerEvaluator>(t_min, t_max),
t_min,
t_max),
AlmostEquals(-381.25522770148542400, 71));
AlmostEquals(-381.25522770148542400, 7));
}

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


#pragma once

#include "quantities/named_quantities.hpp"

#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>
Primitive<std::invoke_result_t<Function, Argument>, Argument> GaussLegendre(
Function const& function,
Function const& f,
Argument const& lower_bound,
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.
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);

// |points| must be of the form 2ᵖ + 1 for some p ∈ ℕ. Returns the
// Clenshaw-Curtis quadrature of f with the given number of points.
template<int points, typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> ClenshawCurtis(
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound);

template<typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> Midpoint(
Function const& function,
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
int intervals);

} // namespace internal_quadrature

using internal_quadrature::AutomaticClenshawCurtis;
using internal_quadrature::GaussLegendre;
using internal_quadrature::Midpoint;

122 changes: 114 additions & 8 deletions numerics/quadrature_body.hpp
Original file line number Diff line number Diff line change
@@ -1,20 +1,30 @@


#pragma once

#include "numerics/quadrature.hpp"

#include <memory>
#include <vector>

#include "numerics/fast_fourier_transform.hpp"
#include "numerics/gauss_legendre_weights.mathematica.h"
#include "numerics/legendre_roots.mathematica.h"
#include "numerics/quadrature.hpp"
#include "quantities/elementary_functions.hpp"

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

using base::FloorLog2;
using quantities::Angle;
using quantities::Cos;
using quantities::Difference;
using quantities::si::Radian;

template<int points, typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> Gauss(
Function const& function,
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
double const* const nodes,
@@ -24,35 +34,131 @@ Primitive<std::invoke_result_t<Function, Argument>, Argument> Gauss(
for (int i = 0; i < points; ++i) {
Argument const scaled_node = lower_bound + half_width * (nodes[i] + 1);
// TODO(phl): Consider compensated summation.
result += weights[i] * function(scaled_node);
result += weights[i] * f(scaled_node);
}
return result * half_width;
}

template<int points, typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> GaussLegendre(
Function const& function,
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound) {
static_assert(points < LegendreRoots.size,
"No table for Gauss-Legendre with the chosen number of points");
return Gauss<points>(function,
return Gauss<points>(f,
lower_bound,
upper_bound,
LegendreRoots[points],
GaussLegendreWeights[points]);
}

template<int points, typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument>
AutomaticClenshawCurtisImplementation(
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
double const relative_tolerance,
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) {
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);
}
}
return estimate;
}

template<int initial_points, 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 const relative_tolerance) {
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);
}

template<int points, typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> ClenshawCurtis(
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound) {
// We follow the notation from [Gen72b] and [Gen72c].
using Value = std::invoke_result_t<Function, Argument>;

constexpr int N = points - 1;
constexpr int log2_N = FloorLog2(N);
static_assert(N == 1 << log2_N);

Difference<Argument> const half_width = (upper_bound - lower_bound) / 2;

constexpr Angle N⁻¹π = π * Radian / N;

// We use a discrete Fourier transform rather than a cosine transform, see
// [Gen72c], equation (3).
// TODO(egg): Consider a discrete cosine transform, ideally incrementally
// computed to improve the performance of the automatic quadrature.

// If we identify [lower_bound, upper_bound] with [-1, 1],
// f_cos_N⁻¹π[s] is f(cos πs/N).
std::vector<Value> f_cos_N⁻¹π;
f_cos_N⁻¹π.resize(2 * N);
for (int s = 0; s <= N; ++s) {
// The N + 1 evaluations.
f_cos_N⁻¹π[s] = f(lower_bound + half_width * (1 + Cos(N⁻¹π * s)));
}
for (int s = N + 1; s <= 2 * N - 1; ++s) {
f_cos_N⁻¹π[s] = f_cos_N⁻¹π[2 * N - s]; // [Gen72c] (5).
}

auto const fft = std::make_unique<FastFourierTransform<Value, Angle, 2 * N>>(
f_cos_N⁻¹π, N⁻¹π);
auto const& a = *fft;

// [Gen72b] equation (7), factoring out the division by N.
Value Σʺ{};
for (std::int64_t n = 0; n <= N; n += 2) {
// The notation g is from [OLBC10], 3.5.17.
int const gₙ = n == 0 || n == N ? 1 : 2;
Σʺ += a[n].real_part() * gₙ / (1 - n * n);
}
Σʺ /= N;

return Σʺ * half_width;
}

template<typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> Midpoint(
Function const& function,
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
int const intervals) {
Difference<Argument> const h = (upper_bound - lower_bound) / intervals;
Primitive<std::invoke_result_t<Function, Argument>, Argument> result{};
for (int i = 0; i < intervals; ++i) {
result += function(lower_bound + (i + 0.5) * h) * h;
result += f(lower_bound + (i + 0.5) * h) * h;
}
return result;
}
34 changes: 32 additions & 2 deletions numerics/quadrature_test.cpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@

#include "numerics/quadrature.hpp"

#include <limits>

#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "quantities/elementary_functions.hpp"
@@ -20,11 +22,16 @@ using testing_utilities::AlmostEquals;
using testing_utilities::IsNear;
using testing_utilities::RelativeErrorFrom;
using testing_utilities::operator""_⑴;
using ::testing::Eq;

class QuadratureTest : public ::testing::Test {};

TEST_F(QuadratureTest, Sin) {
auto const f = [](Angle const x) { return Sin(x); };
int evaluations = 0;
auto const f = [&evaluations](Angle const x) {
++evaluations;
return Sin(x);
};
auto const ʃf = (Cos(2.0 * Radian) - Cos(5.0 * Radian)) * Radian;
EXPECT_THAT(GaussLegendre<1>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(10_⑴)));
@@ -56,9 +63,17 @@ TEST_F(QuadratureTest, Sin) {
AlmostEquals(ʃf, 1, 2));
EXPECT_THAT(GaussLegendre<15>(f, -2.0 * Radian, 5.0 * Radian),
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(evaluations, Eq(134));
}

TEST_F(QuadratureTest, Sin10) {
TEST_F(QuadratureTest, Sin2) {
auto const f = [](Angle const x) { return Sin(2 * x); };
auto const ʃf = (Cos(4 * Radian) - Cos(10 * Radian)) / 2 * Radian;
EXPECT_THAT(GaussLegendre<1>(f, -2.0 * Radian, 5.0 * Radian),
@@ -95,6 +110,21 @@ TEST_F(QuadratureTest, Sin10) {
AlmostEquals(ʃf, 152, 159));
}

TEST_F(QuadratureTest, Sin10) {
int evaluations = 0;
auto const f = [&evaluations](Angle const x) {
++evaluations;
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));
}

} // namespace quadrature
} // namespace numerics
} // namespace principia