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

Commits on Dec 2, 2020

  1. Copy the full SHA
    408037f View commit details

Commits on Dec 3, 2020

  1. Merge pull request #2805 from pleroy/Fourier

    Cache the evaluations at the Gauss-Legendre points
    pleroy authored Dec 3, 2020
    Copy the full SHA
    eb88861 View commit details
Showing with 39 additions and 15 deletions.
  1. +3 −2 numerics/frequency_analysis_body.hpp
  2. +9 −0 numerics/piecewise_poisson_series.hpp
  3. +27 −13 numerics/piecewise_poisson_series_body.hpp
5 changes: 3 additions & 2 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -176,8 +176,9 @@ IncrementalProjection(Function const& function,

Norm const Aₘ = InnerProduct(f, q[m], weight, t_min, t_max);

f -= Aₘ * q[m];
F += Aₘ * q[m];
auto const Aₘqₘ = Aₘ * q[m];
f -= Aₘqₘ;
F += Aₘqₘ;
}

ω = calculator(f);
9 changes: 9 additions & 0 deletions numerics/piecewise_poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -2,6 +2,7 @@
#pragma once

#include <algorithm>
#include <memory>
#include <optional>
#include <string>
#include <vector>
@@ -131,6 +132,14 @@ class PiecewisePoissonSeries {
std::vector<Series> series_;
std::optional<Series> addend_;

// A cache of the evaluation of this function at the Gauss-Legendre points.
// Caching is possible because the Gauss-Legendre integration picks its points
// deterministically. Note that caching keeps working if new series get
// appended to this function: the cache will just be "too short" and new
// entries can be appended as needed. We use a shared_ptr<> to make sure that
// copies remain cheap.
mutable std::shared_ptr<std::vector<Value>> gauss_legendre_cache_;

template<typename V, int ar, int pr,
template<typename, typename, int> class E>
PiecewisePoissonSeries<V, ar, pr, E>
40 changes: 27 additions & 13 deletions numerics/piecewise_poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -12,6 +12,7 @@ namespace principia {
namespace numerics {
namespace internal_piecewise_poisson_series {

using quantities::Angle;
using quantities::Cos;
using quantities::Sin;

@@ -112,24 +113,37 @@ template<typename Value,
auto
PiecewisePoissonSeries<Value, aperiodic_degree_, periodic_degree_, Evaluator>::
FourierTransform() const -> Spectrum {
// TODO(egg): consider pre-evaluating |*this| at all points used by the
// Gaussian quadratures, removing the lifetime requirement on |*this| and
// potentially speeding up repeated evaluations of the Fourier transform.
static constexpr int gauss_legendre_points =
std::max(1, (std::max(aperiodic_degree_, periodic_degree_) + 1) / 2);

return [this](AngularFrequency const& ω) {
// Allocate a cache if there is none. Make sure that the cache has enough
// storage for all the points that we'll evaluate.
if (gauss_legendre_cache_ == nullptr) {
gauss_legendre_cache_ = std::make_shared<std::vector<Value>>();
}
gauss_legendre_cache_->reserve(series_.size() * gauss_legendre_points);

Interval<Instant> const time_domain{t_min(), t_max()};
Instant const t0 = time_domain.midpoint();
Primitive<Complexification<Value>, Instant> integral;
int cache_index = 0;
for (int k = 0; k < series_.size(); ++k) {
integral +=
quadrature::GaussLegendre<std::max(1, (aperiodic_degree_ + 1) / 2)>(
[this, &f = series_[k], t0, ω](
Instant const& t) -> Complexification<Value> {
return (f(t) + EvaluateAddend(t)) *
Complexification<double>{Cos(ω * (t - t0)),
-Sin(ω * (t - t0))};
},
bounds_[k],
bounds_[k + 1]);
integral += quadrature::GaussLegendre<gauss_legendre_points>(
[this, &cache_index, &f = series_[k], t0, ω](
Instant const& t) -> Complexification<Value> {
Angle const θ = ω * (t - t0);
auto const e⁻ⁱᶿ = Complexification<double>{Cos(θ), -Sin(θ)};

// If we reach a point that has not been cached, evaluate it now and
// cache the result.
if (gauss_legendre_cache_->size() <= cache_index) {
gauss_legendre_cache_->push_back(f(t) + EvaluateAddend(t));
}
return gauss_legendre_cache_->at(cache_index++) * e⁻ⁱᶿ;
},
bounds_[k],
bounds_[k + 1]);
}
return integral;
};