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

Commits on Aug 11, 2020

  1. Copy the full SHA
    c6efa74 View commit details
  2. Copy the full SHA
    29822c6 View commit details
  3. Merge pull request #2672 from pleroy/PiecewiseTest

    A test for the projection of piecewise Poisson series
    pleroy authored Aug 11, 2020
    Copy the full SHA
    24ea8de View commit details
Showing with 156 additions and 27 deletions.
  1. +79 −22 numerics/frequency_analysis_test.cpp
  2. +31 −2 numerics/poisson_series.hpp
  3. +46 −3 numerics/poisson_series_body.hpp
101 changes: 79 additions & 22 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -68,8 +68,8 @@ template<typename LFunction, typename RFunction, typename Weight>
Product<std::invoke_result_t<LFunction, Instant>,
std::invoke_result_t<RFunction, Instant>>
DotImplementation::operator()(LFunction const& left,
RFunction const& right,
Weight const& weight) const {
RFunction const& right,
Weight const& weight) const {
return Dot(left, right, weight, t_min_, t_max_);
}

@@ -81,8 +81,8 @@ TEST_F(FrequencyAnalysisTest, PreciseMode) {
AngularFrequency const ω = 666.543 * π / FFT::size * Radian / Second;
Time const Δt = 1 * Second;
std::mt19937_64 random(42);
std::uniform_real_distribution<> amplitude_noise(-0.1, 0.1);
std::uniform_real_distribution<> frequency_noise(-100.0, 100.0);
std::uniform_real_distribution<> amplitude_distribution(-0.1, 0.1);
std::uniform_real_distribution<> frequency_distribution(-100.0, 100.0);

using Series = PoissonSeries<Length, 0, HornerEvaluator>;
Series::PolynomialsByAngularFrequency polynomials;
@@ -95,15 +95,16 @@ TEST_F(FrequencyAnalysisTest, PreciseMode) {

// Noise with lower amplitude and higher frequency.
for (int i = 0; i < 10; ++i) {
auto const sin_amplitude = amplitude_noise(random) * Metre;
auto const cos_amplitude = amplitude_noise(random) * Metre;
auto const sin_amplitude = amplitude_distribution(random) * Metre;
auto const cos_amplitude = amplitude_distribution(random) * Metre;
polynomials.emplace(
ω * frequency_noise(random),
ω * frequency_distribution(random),
Series::Polynomials{/*sin=*/Series::Polynomial({sin_amplitude}, t0),
/*cos=*/Series::Polynomial({cos_amplitude}, t0)});
}
Series const sin(Series::Polynomial({amplitude_noise(random) * Metre}, t0),
polynomials);
Series const sin(
Series::Polynomial({amplitude_distribution(random) * Metre}, t0),
polynomials);

Instant const t_min = t0;
Instant const t_max = t0 + (FFT::size - 1) * Δt;
@@ -117,33 +118,31 @@ TEST_F(FrequencyAnalysisTest, PreciseMode) {

// The FFT gives us an accuracy which is of the order of the number of points.
auto const mode = transform->Mode();
EXPECT_THAT(mode.midpoint(),
RelativeErrorFrom(ω, IsNear(8.1e-4_⑴)));
EXPECT_THAT(mode.midpoint(), RelativeErrorFrom(ω, IsNear(8.1e-4_⑴)));

DotImplementation dot(t_min, t_max);

// The precise analysis is only limited by our ability to pinpoint the
// maximum.
auto const precise_mode = PreciseMode(
mode, sin, apodization::Hann<HornerEvaluator>(t_min, t_max), dot);
EXPECT_THAT(precise_mode,
RelativeErrorFrom(ω, IsNear(6.4e-11_⑴)));
EXPECT_THAT(precise_mode, RelativeErrorFrom(ω, IsNear(6.4e-11_⑴)));
}

TEST_F(FrequencyAnalysisTest, Projection) {
TEST_F(FrequencyAnalysisTest, PoissonSeriesProjection) {
Instant const t0;
AngularFrequency const ω = 666.543 * π * Radian / Second;
std::mt19937_64 random(42);
std::uniform_real_distribution<> noise(-10.0, 10.0);
std::uniform_real_distribution<> amplitude_distribution(-10.0, 10.0);

using Series = PoissonSeries<Length, 4, HornerEvaluator>;

auto random_polynomial = [&noise, &random, &t0]() {
auto const c0 = noise(random) * Metre;
auto const c1 = noise(random) * Metre / Second;
auto const c2 = noise(random) * Metre / Pow<2>(Second);
auto const c3 = noise(random) * Metre / Pow<3>(Second);
auto const c4 = noise(random) * Metre / Pow<4>(Second);
auto random_polynomial = [&amplitude_distribution, &random, &t0]() {
auto const c0 = amplitude_distribution(random) * Metre;
auto const c1 = amplitude_distribution(random) * Metre / Second;
auto const c2 = amplitude_distribution(random) * Metre / Pow<2>(Second);
auto const c3 = amplitude_distribution(random) * Metre / Pow<3>(Second);
auto const c4 = amplitude_distribution(random) * Metre / Pow<4>(Second);

return Series::Polynomial({c0, c1, c2, c3, c4}, t0);
};
@@ -156,7 +155,7 @@ TEST_F(FrequencyAnalysisTest, Projection) {

Instant const t_min = t0;
Instant const t_max = t0 + 100 * Radian / ω;
DotImplementation dot(t_min, t_max);
DotImplementation const dot(t_min, t_max);

// Projection on a 4-th degree basis accurately reconstructs the function.
auto const projection4 = Projection<4>(
@@ -184,6 +183,64 @@ TEST_F(FrequencyAnalysisTest, Projection) {
}
}

TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {
Instant const t0;
AngularFrequency const ω = 666.543 * π * Radian / Second;
std::mt19937_64 random(42);
std::uniform_real_distribution<> amplitude_distribution(-10.0, 10.0);
std::uniform_real_distribution<> perturbation_distribution(-1e-6, 1e-6);

using Series = PoissonSeries<Length, 4, HornerEvaluator>;
using PiecewiseSeries = PiecewisePoissonSeries<Length, 4, HornerEvaluator>;

auto random_polynomial =
[&random, &t0](std::uniform_real_distribution<> const distribution) {
auto const c0 = distribution(random) * Metre;
auto const c1 = distribution(random) * Metre / Second;
auto const c2 = distribution(random) * Metre / Pow<2>(Second);
auto const c3 = distribution(random) * Metre / Pow<3>(Second);
auto const c4 = distribution(random) * Metre / Pow<4>(Second);

return Series::Polynomial({c0, c1, c2, c3, c4}, t0);
};

auto const sin = random_polynomial(amplitude_distribution);
auto const cos = random_polynomial(amplitude_distribution);
Series const series(
Series::Polynomial(Series::Polynomial::Coefficients{}, t0),
{{ω, Series::Polynomials{sin, cos}}});

// Build a series that is based on |series| with different perturbations over
// different intervals.
PiecewiseSeries piecewise_series({t0, t0 + 1 * Second}, series);
for (int i = 1; i < 3; ++i) {
auto const perturbation_sin = random_polynomial(perturbation_distribution);
auto const perturbation_cos = random_polynomial(perturbation_distribution);
Series const perturbation_series(
Series::Polynomial(Series::Polynomial::Coefficients{}, t0),
{{ω, Series::Polynomials{perturbation_sin, perturbation_cos}}});
piecewise_series.Append({t0 + i * Second, t0 + (i + 1) * Second},
series + perturbation_series);
}

Instant const t_min = piecewise_series.t_min();
Instant const t_max = piecewise_series.t_max();
DotImplementation const dot(t_min, t_max);

// Projection on a 4-th degree basis. The errors are of the order of the
// perturbation.
auto const projection4 =
Projection<4>(ω,
piecewise_series,
apodization::Hann<HornerEvaluator>(t_min, t_max),
dot);
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(projection4(t0 + i * Radian / ω),
RelativeErrorFrom(series(t0 + i * Radian / ω),
AllOf(Gt(2.1e-7), Lt(8.8e-4))));
}
}

} // namespace frequency_analysis
} // namespace numerics
} // namespace principia
33 changes: 31 additions & 2 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -209,6 +209,11 @@ class PiecewisePoissonSeries {
// t must be in the interval [t_min, t_max[.
Value operator()(Instant const& t) 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>
PiecewisePoissonSeries& operator-=(PoissonSeries<V, d, E> const& right);

private:
PiecewisePoissonSeries(std::vector<Instant> const& bounds,
std::vector<Series> const& series);
@@ -265,13 +270,17 @@ class PiecewisePoissonSeries {
Product<L, R>
friend Dot(PoissonSeries<L, l, E> const& left,
PiecewisePoissonSeries<R, r, E> const& right,
PoissonSeries<double, w, E> const& weight);
PoissonSeries<double, w, E> const& weight,
Instant const& t_min,
Instant const& t_max);
template<typename L, typename R, int l, int r, int w,
template<typename, typename, int> class E>
Product<L, R>
friend Dot(PiecewisePoissonSeries<L, l, E> const& left,
PoissonSeries<R, r, E> const& right,
PoissonSeries<double, w, E> const& weight);
PoissonSeries<double, w, E> const& weight,
Instant const& t_min,
Instant const& t_max);
};

// Some of the vector space operations for piecewise Poisson series. Note that
@@ -356,6 +365,16 @@ Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight);

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
Product<LValue, RValue>
Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max);

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
@@ -364,6 +383,16 @@ Dot(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight);

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
Product<LValue, RValue>
Dot(PiecewisePoissonSeries<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);

} // namespace internal_poisson_series

using internal_poisson_series::PiecewisePoissonSeries;
49 changes: 46 additions & 3 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -36,7 +36,6 @@ AngularFrequencyPrimitive(
AngularFrequency const& ω,
typename PoissonSeries<Value, degree_, Evaluator>::Polynomials const&
polynomials) {
using Argument = PoissonSeries<Value, degree_, Evaluator>;
using Result = PoissonSeries<Primitive<Value, Time>, degree_ + 1, Evaluator>;

// Integration by parts.
@@ -496,6 +495,26 @@ Value PiecewisePoissonSeries<Value, degree_, Evaluator>::operator()(
return series_[it - bounds_.cbegin() - 1](t);
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<typename V, int d, template<typename, typename, int> class E>
PiecewisePoissonSeries<Value, degree_, Evaluator>&
PiecewisePoissonSeries<Value, degree_, Evaluator>::operator+=(
PoissonSeries<V, d, E> const& right) {
*this = *this + right;
return *this;
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<typename V, int d, template<typename, typename, int> class E>
PiecewisePoissonSeries<Value, degree_, Evaluator>&
PiecewisePoissonSeries<Value, degree_, Evaluator>::operator-=(
PoissonSeries<V, d, E> const& right) {
*this = *this - right;
return *this;
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PiecewisePoissonSeries<Value, degree_, Evaluator>::PiecewisePoissonSeries(
@@ -658,14 +677,26 @@ Product<LValue, RValue> Dot(
PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight) {
return Dot(left, right, weight, right.t_min(), right.t_max());
}

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
Product<LValue, RValue> Dot(
PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
using Result = Primitive<Product<LValue, RValue>, Time>;
Result result;
for (int i = 0; i < right.series_.size(); ++i) {
auto const integrand = left * right.series_[i] * weight;
auto const primitive = integrand.Primitive();
result += primitive(right.bounds_[i + 1]) - primitive(right.bounds_[i]);
}
return result / (right.t_max() - right.t_min());
return result / (t_max - t_min);
}

template<typename LValue, typename RValue,
@@ -675,14 +706,26 @@ Product<LValue, RValue> Dot(
PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight) {
return Dot(left, right, weight, left.t_min(), left.t_max());
}

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
Product<LValue, RValue> Dot(
PiecewisePoissonSeries<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) {
using Result = Primitive<Product<LValue, RValue>, Time>;
Result result;
for (int i = 0; i < left.series_.size(); ++i) {
auto const integrand = left.series_[i] * right * weight;
auto const primitive = integrand.Primitive();
result += primitive(left.bounds_[i + 1]) - primitive(left.bounds_[i]);
}
return result / (left.t_max() - left.t_min());
return result / (t_max - t_min);
}

} // namespace internal_poisson_series