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: f9c97d3ac233
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 6c66f2f79d6d
Choose a head ref
  • 11 commits
  • 5 files changed
  • 1 contributor

Commits on Oct 24, 2020

  1. Copy the full SHA
    f2af8b2 View commit details
  2. More sampling for the FFT.

    pleroy committed Oct 24, 2020
    Copy the full SHA
    6cca88f View commit details
  3. Copy the full SHA
    7df1c63 View commit details
  4. Cleanup fast/slow norm.

    pleroy committed Oct 24, 2020
    Copy the full SHA
    69581f8 View commit details
  5. Copy the full SHA
    8536c5b View commit details
  6. Handle no slow terms.

    pleroy committed Oct 24, 2020
    Copy the full SHA
    c03c8ab View commit details
  7. Different degrees in test.

    # Conflicts:
    #	physics/analytical_series_test.cpp
    pleroy committed Oct 24, 2020
    Copy the full SHA
    96702ba View commit details
  8. Comments and cleanups.

    pleroy committed Oct 24, 2020
    Copy the full SHA
    fdedfe9 View commit details
  9. Adjust test tolerances.

    pleroy committed Oct 24, 2020
    Copy the full SHA
    646cf96 View commit details
  10. Test cleanup.

    pleroy committed Oct 24, 2020
    Copy the full SHA
    660d1c5 View commit details

Commits on Oct 25, 2020

  1. Merge pull request #2770 from pleroy/FastSlow

    A hybrid algorithm for InnerProduct and Norm of Poisson series
    pleroy authored Oct 25, 2020
    Copy the full SHA
    6c66f2f View commit details
Showing with 202 additions and 89 deletions.
  1. +2 −2 numerics/frequency_analysis_test.cpp
  2. +8 −0 numerics/poisson_series.hpp
  3. +111 −24 numerics/poisson_series_body.hpp
  4. +1 −1 numerics/poisson_series_test.cpp
  5. +80 −62 physics/analytical_series_test.cpp
4 changes: 2 additions & 2 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -488,8 +488,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.5e-17 * Metre),
Lt(7.5e-14 * Metre)))
: AllOf(Gt(1.5e-16 * Metre),
Lt(9.7e-14 * Metre)))
<< ω_index;
}
if (ω_index == ωs.size()) {
8 changes: 8 additions & 0 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -151,6 +151,14 @@ class PoissonSeries {
Polynomial aperiodic,
PolynomialsByAngularFrequency periodic);

// Splits this series into two copies, with frequencies lower and higher than
// ω_cutoff, respectively.
struct SplitPoissonSeries {
PoissonSeries slow;
PoissonSeries fast;
};
SplitPoissonSeries Split(AngularFrequency const& ω_cutoff) const;

Instant origin_; // Common to all polynomials.
Polynomial aperiodic_;
// The frequencies in this vector are positive, distinct and in increasing
135 changes: 111 additions & 24 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -10,6 +10,7 @@
#include <utility>
#include <vector>

#include "numerics/double_precision.hpp"
#include "numerics/quadrature.hpp"
#include "numerics/ulp_distance.hpp"
#include "quantities/elementary_functions.hpp"
@@ -33,8 +34,21 @@ namespace si = quantities::si;

// These parameters have been tuned for approximation of the Moon over 3 months
// with 10 periods.

// In the fast/slow algorithms, specifies the maximum number of periods over the
// time interval below which we use Clenshaw-Curtis integration.
constexpr int clenshaw_curtis_max_periods_overall = 40;

// The minimum value of the max_point parameter passed to Clenshaw-Curtis
// integration, irrespective of the frequencies of the argument function.
constexpr int clenshaw_curtis_min_points_overall = 33;

// The maximum number of points use in Clenshaw-Curtis integration for each
// period of the highest frequency of the argument function.
constexpr int clenshaw_curtis_point_per_period = 4;

// The desired relative error on Clenshaw-Curtis integration, as determined by
// two successive computations with increasing number of points.
constexpr double clenshaw_curtis_relative_error = 0x1p-32;

template<typename Value, int degree_,
@@ -218,7 +232,9 @@ PoissonSeries<Value, degree_, Evaluator>::Integrate(Instant const& t1,
{{ω,
{/*sin=*/typename FirstPart::Polynomial(polynomials.cos),
/*cos=*/typename FirstPart::Polynomial(-polynomials.sin)}}});
result += (first_part(t2) - first_part(t1)) / ω * Radian;
DoublePrecision<Value> sum;
sum += first_part(t2);
sum -= first_part(t1);

if constexpr (degree_ != 0) {
auto const sin_polynomial =
@@ -229,8 +245,9 @@ PoissonSeries<Value, degree_, Evaluator>::Integrate(Instant const& t1,
{{ω,
{/*sin=*/sin_polynomial,
/*cos=*/cos_polynomial}}});
result += second_part.Integrate(t1, t2) / ω * Radian;
sum += second_part.Integrate(t1, t2);
}
result += (sum.value + sum.error) / ω * Radian;
}
return result;
}
@@ -343,6 +360,49 @@ PoissonSeries<Value, degree_, Evaluator>::PoissonSeries(
aperiodic_(std::move(aperiodic)),
periodic_(std::move(periodic)) {}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
typename PoissonSeries<Value, degree_, Evaluator>::SplitPoissonSeries
PoissonSeries<Value, degree_, Evaluator>::Split(
AngularFrequency const& ω_cutoff) const {
// TODO(phl): Should we try to avoid a linear search and copies?
typename PoissonSeries::PolynomialsByAngularFrequency slow_periodic;
typename PoissonSeries::PolynomialsByAngularFrequency fast_periodic;
for (auto const& [ω, polynomials] : periodic_) {
if (ω <= ω_cutoff) {
slow_periodic.emplace_back(ω, polynomials);
} else {
fast_periodic.emplace_back(ω, polynomials);
}
}

// The reason for having the slow/fast split is to handle cancellations that
// occurs between the aperiodic component and the components with low
// frequencies. If there are no low frequencies, we might as well put the
// aperiodic component in the high-frequencies half.
if (slow_periodic.empty()) {
PoissonSeries slow(TrustedPrivateConstructor{},
typename PoissonSeries::Polynomial(
typename PoissonSeries::Polynomial::Coefficients{},
aperiodic_.origin()),
std::move(slow_periodic));
PoissonSeries fast(TrustedPrivateConstructor{},
aperiodic_,
std::move(fast_periodic));
return {/*slow=*/std::move(slow), /*fast=*/std::move(fast)};
} else {
PoissonSeries slow(TrustedPrivateConstructor{},
aperiodic_,
std::move(slow_periodic));
PoissonSeries fast(TrustedPrivateConstructor{},
typename PoissonSeries::Polynomial(
typename PoissonSeries::Polynomial::Coefficients{},
aperiodic_.origin()),
std::move(fast_periodic));
return {/*slow=*/std::move(slow), /*fast=*/std::move(fast)};
}
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<int higher_degree_,
@@ -374,9 +434,13 @@ PoissonSeries<Value, degree_, Evaluator>::Norm(
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) const {
AngularFrequency const ω_cutoff =
2 * π * Radian * clenshaw_curtis_max_periods_overall / (t_max - t_min);
auto const split = Split(ω_cutoff);

AngularFrequency const max_ω =
(periodic_.empty() ? AngularFrequency{}
: 2 * periodic_.back().first) +
(split.slow.periodic_.empty() ? AngularFrequency{}
: 2 * split.slow.periodic_.back().first) +
(weight.periodic_.empty() ? AngularFrequency{}
: weight.periodic_.back().first);
std::optional<int> max_points =
@@ -387,15 +451,21 @@ PoissonSeries<Value, degree_, Evaluator>::Norm(
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>::Norm²((*this)(t)) * weight(t);
auto slow_integrand = [&split, &weight](Instant const& t) {
return Hilbert<Value>::Norm²(split.slow(t)) * weight(t);
};
return Sqrt(quadrature::AutomaticClenshawCurtis(
integrand,
t_min, t_max,
/*max_relative_error=*/clenshaw_curtis_relative_error,
/*max_points=*/max_points) /
(t_max - t_min));
auto const slow_quadrature = quadrature::AutomaticClenshawCurtis(
slow_integrand,
t_min, t_max,
/*max_relative_error=*/clenshaw_curtis_relative_error,
/*max_points=*/max_points);

auto const fast_integrand =
PointwiseInnerProduct(split.fast, split.fast + 2 * split.slow) *
weight;
auto const fast_quadrature = fast_integrand.Integrate(t_min, t_max);

return Sqrt((slow_quadrature + fast_quadrature) / (t_max - t_min));
}

template<typename Value, int degree_,
@@ -694,11 +764,18 @@ typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
AngularFrequency const ω_cutoff =
2 * π * Radian * clenshaw_curtis_max_periods_overall / (t_max - t_min);
auto const left_split = left.Split(ω_cutoff);
auto const right_split = right.Split(ω_cutoff);

AngularFrequency const max_ω =
(left.periodic_.empty() ? AngularFrequency{}
: left.periodic_.back().first) +
(right.periodic_.empty() ? AngularFrequency{}
: right.periodic_.back().first) +
(left_split.slow.periodic_.empty()
? AngularFrequency{}
: left_split.slow.periodic_.back().first) +
(right_split.slow.periodic_.empty()
? AngularFrequency{}
: right_split.slow.periodic_.back().first) +
(weight.periodic_.empty() ? AngularFrequency{}
: weight.periodic_.back().first);
std::optional<int> max_points =
@@ -709,15 +786,25 @@ typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
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);
auto slow_integrand = [&left_split, &right_split, &weight](Instant const& t) {
return Hilbert<LValue, RValue>::InnerProduct(left_split.slow(t),
right_split.slow(t)) *
weight(t);
};
return quadrature::AutomaticClenshawCurtis(
integrand,
t_min, t_max,
/*max_relative_error=*/clenshaw_curtis_relative_error,
/*max_points=*/max_points) /
(t_max - t_min);
auto const slow_quadrature = quadrature::AutomaticClenshawCurtis(
slow_integrand,
t_min, t_max,
/*max_relative_error=*/clenshaw_curtis_relative_error,
/*max_points=*/max_points);

auto const fast_integrand =
(PointwiseInnerProduct(left_split.fast, right_split.slow) +
PointwiseInnerProduct(left_split.slow, right_split.fast) +
PointwiseInnerProduct(left_split.fast, right_split.fast)) *
weight;
auto const fast_quadrature = fast_integrand.Integrate(t_min, t_max);

return (slow_quadrature + fast_quadrature) / (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
@@ -310,7 +310,7 @@ TEST_F(PoissonSeriesTest, Primitive) {
EXPECT_THAT(
pb_->Integrate(t0_ + 5 * Second, t0_ + 13 * Second),
AlmostEquals(expected_primitive(13 * Second) -
expected_primitive(5 * Second), 1, 2));
expected_primitive(5 * Second), 0));
}

TEST_F(PoissonSeriesTest, InnerProduct) {
Loading