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

Commits on Sep 27, 2020

  1. Copy the full SHA
    c2dad76 View commit details
  2. merge

    eggrobin committed Sep 27, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    59cf89e View commit details
  3. Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    8afc111 View commit details

Commits on Sep 28, 2020

  1. Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    1dd7bab View commit details
  2. iwyu

    eggrobin committed Sep 28, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    bd0b372 View commit details
  3. After pleroy’s review

    eggrobin committed Sep 28, 2020
    Copy the full SHA
    a0c771e View commit details
  4. convention

    eggrobin committed Sep 28, 2020
    Copy the full SHA
    3637337 View commit details
  5. Merge pull request #2738 from eggrobin/fourier

    Fourier
    eggrobin authored Sep 28, 2020
    Copy the full SHA
    1832ab1 View commit details
Showing with 118 additions and 12 deletions.
  1. +7 −6 geometry/cartesian_product_body.hpp
  2. +1 −1 geometry/interval_body.hpp
  3. +19 −0 numerics/poisson_series.hpp
  4. +26 −0 numerics/poisson_series_body.hpp
  5. +65 −5 numerics/poisson_series_test.cpp
13 changes: 7 additions & 6 deletions geometry/cartesian_product_body.hpp
Original file line number Diff line number Diff line change
@@ -144,7 +144,9 @@ class CartesianProductVectorSpace<Scalar, Tuple,
Tuple const& left,
Scalar const& right);

FORCE_INLINE(static constexpr) Apply<Quotient, Tuple> Divide(
// SFINAEd out unless the division is defined.
template<typename TupleQuotient = Apply<Quotient, Tuple>>
FORCE_INLINE(static constexpr) TupleQuotient Divide(
Tuple const& left,
Scalar const& right);
};
@@ -168,11 +170,10 @@ constexpr auto CartesianProductVectorSpace<
}

template<typename Scalar, typename Tuple, std::size_t... indices>
constexpr auto CartesianProductVectorSpace<
Scalar, Tuple,
std::index_sequence<indices...>>::Divide(Tuple const& left,
Scalar const& right)
-> Apply<Quotient, Tuple> {
template<typename TupleQuotient>
constexpr TupleQuotient
CartesianProductVectorSpace<Scalar, Tuple, std::index_sequence<indices...>>::
Divide(Tuple const& left, Scalar const& right) {
return {std::get<indices>(left) / right...};
}

2 changes: 1 addition & 1 deletion geometry/interval_body.hpp
Original file line number Diff line number Diff line change
@@ -17,7 +17,7 @@ Difference<T> Interval<T>::measure() const {

template<typename T>
T Interval<T>::midpoint() const {
return max >= min ? min + measure() / 2 : NaN<T>;
return max >= min ? min + measure() / 2 : min + NaN<Difference<T>>;
}

template<typename T>
19 changes: 19 additions & 0 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -2,12 +2,14 @@
#pragma once

#include <algorithm>
#include <functional>
#include <map>
#include <optional>
#include <string>
#include <utility>
#include <vector>

#include "geometry/complexification.hpp"
#include "geometry/hilbert.hpp"
#include "geometry/interval.hpp"
#include "geometry/named_quantities.hpp"
@@ -51,6 +53,7 @@ FORWARD_DECLARE_FUNCTION_FROM(
namespace numerics {
namespace internal_poisson_series {

using geometry::Complexification;
using geometry::Hilbert;
using geometry::Instant;
using geometry::Interval;
@@ -283,6 +286,8 @@ class PiecewisePoissonSeries {
public:
static constexpr int degree = degree_;
using Series = PoissonSeries<Value, degree_, Evaluator>;
using Spectrum = std::function<Complexification<Primitive<Value, Time>>(
AngularFrequency const&)>;

PiecewisePoissonSeries(Interval<Instant> const& interval,
Series const& series);
@@ -299,6 +304,20 @@ class PiecewisePoissonSeries {
// t must be in the interval [t_min, t_max].
Value operator()(Instant const& t) const;

// Returns the Fourier transform of this piecewise Poisson series.
// The function is taken to be 0 outside [t_min, t_max].
// The convention used is ∫ f(t) exp(-iωt) dt, corresponding to Mathematica’s
// FourierParameters -> {1, -1} for FourierTransform (the “pure mathematics;
// systems engineering”).
// When evaluated at a given frequency, the Fourier transform is computed by
// Gauss-Legendre quadrature on each subinterval, where the number of points
// is chosen assuming that the periods of periodic terms are all large
// compared to the subintervals.
// If apodization is desired, |*this| should be multiplied by an apodization
// function, and |FourierTransform| should be called on the product.
// |*this| must outlive the resulting function.
Spectrum FourierTransform() const;

template<typename V, int d, template<typename, typename, int> class E>
PiecewisePoissonSeries& operator+=(PoissonSeries<V, d, E> const& right);
template<typename V, int d, template<typename, typename, int> class E>
26 changes: 26 additions & 0 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -653,6 +653,32 @@ Value PiecewisePoissonSeries<Value, degree_, Evaluator>::operator()(
return series_[it - bounds_.cbegin() - 1](t);
}

template<typename Value,
int degree_,
template<typename, typename, int> class Evaluator>
auto PiecewisePoissonSeries<Value, 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.
return [this](AngularFrequency const& ω) {
Interval<Instant> const time_domain{t_min(), t_max()};
Instant const t0 = time_domain.midpoint();
Primitive<Complexification<Value>, Instant> integral;
for (int k = 0; k < series_.size(); ++k) {
integral += quadrature::GaussLegendre<std::max(1, (degree_ + 1) / 2)>(
[&f = series_[k], t0, ω](
Instant const& t) -> Complexification<Value> {
return f(t) * Complexification<double>{Cos(ω * (t - t0)),
-Sin(ω * (t - t0))};
},
bounds_[k],
bounds_[k + 1]);
}
return integral;
};
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<typename V, int d, template<typename, typename, int> class E>
70 changes: 65 additions & 5 deletions numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@

#include "numerics/poisson_series.hpp"

#include <functional>
#include <limits>
#include <memory>

@@ -11,10 +12,13 @@
#include "numerics/apodization.hpp"
#include "numerics/polynomial_evaluators.hpp"
#include "numerics/quadrature.hpp"
#include "numerics/root_finders.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/almost_equals.hpp"
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/numerics_matchers.hpp"
#include "testing_utilities/vanishes_before.hpp"

namespace principia {
@@ -38,7 +42,10 @@ using quantities::si::Metre;
using quantities::si::Radian;
using quantities::si::Second;
using testing_utilities::AlmostEquals;
using testing_utilities::IsNear;
using testing_utilities::VanishesBefore;
using testing_utilities::RelativeErrorFrom;
using testing_utilities::operator""_⑴;

class PoissonSeriesTest : public ::testing::Test {
protected:
@@ -47,6 +54,7 @@ class PoissonSeriesTest : public ::testing::Test {
Handedness::Right,
serialization::Frame::TEST>;

using Degree0 = PoissonSeries<double, 0, HornerEvaluator>;
using Degree1 = PoissonSeries<double, 1, HornerEvaluator>;

PoissonSeriesTest()
@@ -199,11 +207,63 @@ TEST_F(PoissonSeriesTest, PointwiseInnerProduct) {
AlmostEquals(5 * Metre * Metre, 0));
}

TEST_F(PoissonSeriesTest, Fourier) {
Degree0::Polynomial constant({1.0}, t0_);
AngularFrequency const ω = 4 * Radian / Second;
PoissonSeries<Displacement<World>, 0, HornerEvaluator> signal(
constant * Displacement<World>{},
{{ω,
{/*sin=*/constant *
Displacement<World>({2 * Metre, -3 * Metre, 5 * Metre}),
/*cos=*/constant *
Displacement<World>({-7 * Metre, 11 * Metre, -13 * Metre})}}});
// Slice our signal into segments short enough that one-point Gauss-Legendre
// (also known as midpoint) does the job.
constexpr int segments = 100;
PiecewisePoissonSeries<Displacement<World>, 0, HornerEvaluator> f(
{t0_, t0_ + π * Second / segments}, signal);
for (int k = 1; k < segments; ++k) {
f.Append({t0_ + k * π * Second / segments,
t0_ + (k + 1) * π * Second / segments},
signal);
}
auto const f_fourier_transform = f.FourierTransform();
auto const f_power_spectrum =
[&f_fourier_transform](AngularFrequency const& ω) {
return f_fourier_transform(ω).Norm²();
};
EXPECT_THAT(Brent(f_power_spectrum, 0.9 * ω, 1.1 * ω, std::greater<>{}),
RelativeErrorFrom(ω, IsNear(0.03_⑴)));
// A peak arising from the finite interval.
EXPECT_THAT(Brent(f_power_spectrum,
0 * Radian / Second,
2 * Radian / Second,
std::greater<>{}),
IsNear(1.209_⑴ * Radian / Second));

auto const fw = f * apodization::Hann<HornerEvaluator>(f.t_min(), f.t_max());

auto const fw_fourier_transform = fw.FourierTransform();
auto const fw_power_spectrum =
[&fw_fourier_transform](AngularFrequency const& ω) {
return fw_fourier_transform(ω).Norm²();
};
EXPECT_THAT(Brent(fw_power_spectrum, 0.9 * ω, 1.1 * ω, std::greater<>{}),
RelativeErrorFrom(ω, IsNear(0.005_⑴)));
// Hann filters out the interval; this search for a second maximum converges
// to its boundary.
EXPECT_THAT(Brent(fw_power_spectrum,
0 * Radian / Second,
2 * Radian / Second,
std::greater<>{}),
IsNear(1.9999999_⑴ * Radian / Second));
}

TEST_F(PoissonSeriesTest, Primitive) {
auto const actual_primitive = pb_->Primitive();

// The primitive was computed using Mathematica.
auto const expected_primitive = [=](Time const& t){
auto const expected_primitive = [=](Time const& t) {
auto const a0 = 3;
auto const a1 = 4 / Second;
auto const b0 = 9;
@@ -400,11 +460,11 @@ TEST_F(PiecewisePoissonSeriesTest, Action) {
auto const p1 = p_ * pp_;
auto const p2 = pp_ * p_;
EXPECT_THAT(p1(t0_ + 0.5 * Second),
AlmostEquals((7 - 4* Sqrt(2))/4, 0, 4));
AlmostEquals((7 - 4 * Sqrt(2)) / 4, 0, 4));
EXPECT_THAT(p1(t0_ + 1.5 * Second),
AlmostEquals((-3 - 3 * Sqrt(2)) / 4, 1));
EXPECT_THAT(p2(t0_ + 0.5 * Second),
AlmostEquals((7 - 4* Sqrt(2))/4, 0, 4));
AlmostEquals((7 - 4 * Sqrt(2)) / 4, 0, 4));
EXPECT_THAT(p2(t0_ + 1.5 * Second),
AlmostEquals((-3 - 3 * Sqrt(2)) / 4, 1));
}
@@ -440,11 +500,11 @@ TEST_F(PiecewisePoissonSeriesTest, ActionMultiorigin) {
auto const p1 = p * pp_;
auto const p2 = pp_ * p;
EXPECT_THAT(p1(t0_ + 0.5 * Second),
AlmostEquals((7 - 4* Sqrt(2))/4, 0, 4));
AlmostEquals((7 - 4 * Sqrt(2)) / 4, 0, 4));
EXPECT_THAT(p1(t0_ + 1.5 * Second),
AlmostEquals((-3 - 3 * Sqrt(2)) / 4, 1));
EXPECT_THAT(p2(t0_ + 0.5 * Second),
AlmostEquals((7 - 4* Sqrt(2))/4, 0, 4));
AlmostEquals((7 - 4 * Sqrt(2)) / 4, 0, 4));
EXPECT_THAT(p2(t0_ + 1.5 * Second),
AlmostEquals((-3 - 3 * Sqrt(2)) / 4, 1));
}