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

Commits on Sep 26, 2020

  1. New PoissonSeries constructor to avoid sorting/normalization.

    # Conflicts:
    #	numerics/poisson_series_body.hpp
    pleroy committed Sep 26, 2020
    Copy the full SHA
    f630004 View commit details
  2. Copy the full SHA
    4a6d13b View commit details
  3. Copy the full SHA
    d708eda View commit details
  4. Specify the number of points to use for Gauss integration. This is a …

    …temporary measure until we have a smarter integrator.
    pleroy committed Sep 26, 2020
    Copy the full SHA
    a848ef1 View commit details
  5. A test that doesn't work.

    pleroy committed Sep 26, 2020
    Copy the full SHA
    1d90dd7 View commit details
  6. Copy the full SHA
    42b48ad View commit details

Commits on Sep 27, 2020

  1. Fix the tests.

    pleroy committed Sep 27, 2020
    Copy the full SHA
    08541c4 View commit details
  2. Obsolete comment.

    pleroy committed Sep 27, 2020
    Copy the full SHA
    2e8485f View commit details
  3. Lint.

    pleroy committed Sep 27, 2020
    Copy the full SHA
    31222a0 View commit details
  4. Merge pull request #2736 from pleroy/PoissonOptimization

    Improve the Poisson series implementation
    pleroy authored Sep 27, 2020
    Copy the full SHA
    ddb36de View commit details
Showing with 110 additions and 49 deletions.
  1. +13 −3 numerics/frequency_analysis_test.cpp
  2. +29 −9 numerics/poisson_series.hpp
  3. +60 −29 numerics/poisson_series_body.hpp
  4. +8 −8 numerics/poisson_series_test.cpp
16 changes: 13 additions & 3 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -86,7 +86,17 @@ typename Hilbert<std::invoke_result_t<LFunction, Instant>,
DotImplementation::operator()(LFunction const& left,
RFunction const& right,
Weight const& weight) const {
return Dot(left, right, weight, t_min_, t_max_);
// We need to use a large number of points otherwise the test
// PiecewisePoissonSeriesProjection doesn't yield a reasonable solution. This
// doesn't happen with real life functions, which are much smoother than the
// test functions.
return Dot<std::invoke_result_t<LFunction, Instant>,
std::invoke_result_t<RFunction, Instant>,
LFunction::degree,
RFunction::degree,
Weight::degree,
HornerEvaluator,
/*points=*/30>(left, right, weight, t_min_, t_max_);
}

class FrequencyAnalysisTest : public ::testing::Test {
@@ -328,7 +338,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesVectorProjection) {
}

TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {
AngularFrequency const ω = 666.543 * π * Radian / Second;
AngularFrequency const ω = 6.66543 * π * 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);
@@ -370,7 +380,7 @@ TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {
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))));
AllOf(Gt(6.2e-9), Lt(3.2e-4))));
}
}

38 changes: 29 additions & 9 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -69,7 +69,7 @@ template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
class PoissonSeries {
public:
static const int degree = degree_;
static constexpr int degree = degree_;
using Polynomial =
numerics::PolynomialInMonomialBasis<Value, Instant, degree_, Evaluator>;

@@ -85,6 +85,8 @@ class PoissonSeries {
using PolynomialsByAngularFrequency =
std::vector<std::pair<AngularFrequency, Polynomials>>;

// The |periodic| vector may contain frequencies that are negative or zero, as
// well as repeated frequencies. It is not expected to be ordered.
PoissonSeries(Polynomial const& aperiodic,
PolynomialsByAngularFrequency const& periodic);

@@ -108,12 +110,21 @@ class PoissonSeries {
PoissonSeries& operator-=(PoissonSeries<V, d, E> const& right);

private:
// Similar to the public constructor, but passing by copy allows moves, which
// is useful for internal algorithms.
struct PrivateConstructor {};

PoissonSeries(PrivateConstructor,
Polynomial aperiodic,
PolynomialsByAngularFrequency periodic);

// Similar to the previous constructor, except that the |periodic| vector is
// used verbatim, without sorting or normalization, which is useful for
// internal algorithms which produce positive, ordered frequencies.
struct TrustedPrivateConstructor {};
PoissonSeries(TrustedPrivateConstructor,
Polynomial aperiodic,
PolynomialsByAngularFrequency periodic);

Instant origin_; // Common to all polynomials.
Polynomial aperiodic_;
// The frequencies in this vector are positive, distinct and in increasing
@@ -248,9 +259,12 @@ std::ostream& operator<<(
// Technically the weight function must be nonnegative for this to be an inner
// product. Not sure how this works with the flat-top windows, which can be
// negative. Note that the result is normalized by dividing by (t_max - t_min).
// NOTE(phl): |points| is currently unused but ensures that the template has the
// same profile as the one for |PiecewisePoissonSeries|.
template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
template<typename, typename, int> class Evaluator,
int points = 0>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
@@ -266,6 +280,7 @@ template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
class PiecewisePoissonSeries {
public:
static constexpr int degree = degree_;
using Series = PoissonSeries<Value, degree_, Evaluator>;

PiecewisePoissonSeries(Interval<Instant> const& interval,
@@ -340,15 +355,15 @@ class PiecewisePoissonSeries {
friend operator*(PiecewisePoissonSeries<L, l, E> const& left,
PoissonSeries<R, r, E> const& right);
template<typename L, typename R, int l, int r, int w,
template<typename, typename, int> class E>
template<typename, typename, int> class E, int p>
typename Hilbert<L, R>::InnerProductType
friend Dot(PoissonSeries<L, l, E> const& left,
PiecewisePoissonSeries<R, r, E> const& right,
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>
template<typename, typename, int> class E, int p>
typename Hilbert<L, R>::InnerProductType
friend Dot(PiecewisePoissonSeries<L, l, E> const& left,
PoissonSeries<R, r, E> const& right,
@@ -439,15 +454,17 @@ operator*(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
template<typename, typename, int> class Evaluator,
int points = (ldegree_ + rdegree_ + wdegree_) / 2>
typename Hilbert<LValue, RValue>::InnerProductType
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>
template<typename, typename, int> class Evaluator,
int points = (ldegree_ + rdegree_ + wdegree_) / 2>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<RValue, rdegree_, Evaluator> const& right,
@@ -457,15 +474,17 @@ Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
template<typename, typename, int> class Evaluator,
int points = (ldegree_ + rdegree_ + wdegree_) / 2>
typename Hilbert<LValue, RValue>::InnerProductType
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>
template<typename, typename, int> class Evaluator,
int points = (ldegree_ + rdegree_ + wdegree_) / 2>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
@@ -475,6 +494,7 @@ Dot(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,

} // namespace internal_poisson_series

using internal_poisson_series::Dot;
using internal_poisson_series::PiecewisePoissonSeries;
using internal_poisson_series::PoissonSeries;

89 changes: 60 additions & 29 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/quadrature.hpp"
#include "numerics/ulp_distance.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/named_quantities.hpp"
@@ -161,14 +162,19 @@ PoissonSeries<Value, degree_, Evaluator>::AtOrigin(
PolynomialsByAngularFrequency periodic;
periodic.reserve(periodic_.size());
for (auto const& [ω, polynomials] : periodic_) {
Polynomial const sin = polynomials.sin.AtOrigin(origin);
Polynomial const cos = polynomials.cos.AtOrigin(origin);
periodic.emplace_back(
ω,
Polynomials{/*sin=*/sin * Cos(ω * shift) - cos * Sin(ω * shift),
/*cos=*/sin * Sin(ω * shift) + cos * Cos(ω * shift)});
double const cos_ω_shift = Cos(ω * shift);
double const sin_ω_shift = Sin(ω * shift);
Polynomial const sin_at_origin = polynomials.sin.AtOrigin(origin);
Polynomial const cos_at_origin = polynomials.cos.AtOrigin(origin);
periodic.emplace_back(ω,
Polynomials{/*sin=*/sin_at_origin * cos_ω_shift -
cos_at_origin * sin_ω_shift,
/*cos=*/sin_at_origin * sin_ω_shift +
cos_at_origin * cos_ω_shift});
}
return {PrivateConstructor{}, std::move(aperiodic), std::move(periodic)};
return {TrustedPrivateConstructor{},
std::move(aperiodic),
std::move(periodic)};
}

template<typename Value, int degree_,
@@ -285,6 +291,16 @@ PoissonSeries<Value, degree_, Evaluator>::PoissonSeries(
}
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, degree_, Evaluator>::PoissonSeries(
TrustedPrivateConstructor,
Polynomial aperiodic,
PolynomialsByAngularFrequency periodic)
: origin_(aperiodic.origin()),
aperiodic_(std::move(aperiodic)),
periodic_(std::move(periodic)) {}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<typename V, int d, template<typename, typename, int> class E>
@@ -326,7 +342,7 @@ operator-(PoissonSeries<Value, rdegree_, Evaluator> const& right) {
typename Result::Polynomials{/*sin=*/-polynomials.sin,
/*cos=*/-polynomials.cos});
}
return {typename Result::PrivateConstructor{},
return {typename Result::TrustedPrivateConstructor{},
std::move(aperiodic),
std::move(periodic)};
}
@@ -379,7 +395,9 @@ operator+(PoissonSeries<Value, ldegree_, Evaluator> const& left,
++it_right;
}
}
return {typename Result::PrivateConstructor{},
// Because we have done a merge on the periodic vectors, we can use the
// trusted constructor.
return {typename Result::TrustedPrivateConstructor{},
std::move(aperiodic),
std::move(periodic)};
}
@@ -432,7 +450,9 @@ operator-(PoissonSeries<Value, ldegree_, Evaluator> const& left,
++it_right;
}
}
return {typename Result::PrivateConstructor{},
// Because we have done a merge on the periodic vectors, we can use the
// trusted constructor.
return {typename Result::TrustedPrivateConstructor{},
std::move(aperiodic),
std::move(periodic)};
}
@@ -472,7 +492,7 @@ operator*(PoissonSeries<Value, degree_, Evaluator> const& left,
typename Result::Polynomials{/*sin=*/polynomials.sin * right,
/*cos=*/polynomials.cos * right});
}
return {typename Result::PrivateConstructor{},
return {typename Result::TrustedPrivateConstructor{},
std::move(aperiodic),
std::move(periodic)};
}
@@ -492,7 +512,7 @@ operator/(PoissonSeries<Value, degree_, Evaluator> const& left,
typename Result::Polynomials{/*sin=*/polynomials.sin / right,
/*cos=*/polynomials.cos / right});
}
return {typename Result::PrivateConstructor{},
return {typename Result::TrustedPrivateConstructor{},
std::move(aperiodic),
std::move(periodic)};
}
@@ -568,7 +588,8 @@ std::ostream& operator<<(

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
template<typename, typename, int> class Evaluator,
int points>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
@@ -832,17 +853,20 @@ operator*(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
template<typename, typename, int> class Evaluator,
int points>
typename Hilbert<LValue, RValue>::InnerProductType
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());
return Dot<LValue, RValue, ldegree_, rdegree_, wdegree_, Evaluator, points>(
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>
template<typename, typename, int> class Evaluator,
int points>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<RValue, rdegree_, Evaluator> const& right,
@@ -853,28 +877,33 @@ Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
Primitive<typename Hilbert<LValue, RValue>::InnerProductType, Time>;
Result result{};
for (int i = 0; i < right.series_.size(); ++i) {
Instant const origin = right.series_[i].origin();
auto const integrand = PointwiseInnerProduct(
left.AtOrigin(origin) * weight.AtOrigin(origin), right.series_[i]);
auto const primitive = integrand.Primitive();
result += primitive(right.bounds_[i + 1]) - primitive(right.bounds_[i]);
auto integrand = [i, &left, &right, &weight](Instant const& t) {
return Hilbert<LValue, RValue>::InnerProduct(left(t) * weight(t),
right.series_[i](t));
};
auto const integral = quadrature::GaussLegendre<points>(
integrand, right.bounds_[i], right.bounds_[i + 1]);
result += integral;
}
return result / (t_max - t_min);
}

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
template<typename, typename, int> class Evaluator,
int points>
typename Hilbert<LValue, RValue>::InnerProductType
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());
return Dot<LValue, RValue, ldegree_, rdegree_, wdegree_, Evaluator, points>(
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>
template<typename, typename, int> class Evaluator,
int points>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
@@ -885,11 +914,13 @@ Dot(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
Primitive<typename Hilbert<LValue, RValue>::InnerProductType, Time>;
Result result{};
for (int i = 0; i < left.series_.size(); ++i) {
Instant const origin = left.series_[i].origin();
auto const integrand = PointwiseInnerProduct(
left.series_[i], right.AtOrigin(origin) * weight.AtOrigin(origin));
auto const primitive = integrand.Primitive();
result += primitive(left.bounds_[i + 1]) - primitive(left.bounds_[i]);
auto integrand = [i, &left, &right, &weight](Instant const& t) {
return Hilbert<LValue, RValue>::InnerProduct(left.series_[i](t),
right(t) * weight(t));
};
auto const integral = quadrature::GaussLegendre<points>(
integrand, left.bounds_[i], left.bounds_[i + 1]);
result += integral;
}
return result / (t_max - t_min);
}
16 changes: 8 additions & 8 deletions numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -451,22 +451,22 @@ TEST_F(PiecewisePoissonSeriesTest, ActionMultiorigin) {
}

TEST_F(PiecewisePoissonSeriesTest, Dot) {
double const d1 = Dot(
double const d1 = Dot<double, double, 0, 0, 0, HornerEvaluator, 8>(
pp_, p_, apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second));
double const d2 = Dot(
double const d2 = Dot<double, double, 0, 0, 0, HornerEvaluator, 8>(
p_, pp_, apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second));
EXPECT_THAT(d1, AlmostEquals((3 * π - 26) / (8 * π), 1));
EXPECT_THAT(d2, AlmostEquals((3 * π - 26) / (8 * π), 1));
EXPECT_THAT(d1, AlmostEquals((3 * π - 26) / (8 * π), 0));
EXPECT_THAT(d2, AlmostEquals((3 * π - 26) / (8 * π), 0));
}

TEST_F(PiecewisePoissonSeriesTest, DotMultiorigin) {
auto const p = p_.AtOrigin(t0_ + 2 * Second);
double const d1 = Dot(
double const d1 = Dot<double, double, 0, 0, 0, HornerEvaluator, 8>(
pp_, p, apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second));
double const d2 = Dot(
double const d2 = Dot<double, double, 0, 0, 0, HornerEvaluator, 8>(
p, pp_, apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second));
EXPECT_THAT(d1, AlmostEquals((3 * π - 26) / (8 * π), 1));
EXPECT_THAT(d2, AlmostEquals((3 * π - 26) / (8 * π), 1));
EXPECT_THAT(d1, AlmostEquals((3 * π - 26) / (8 * π), 0));
EXPECT_THAT(d2, AlmostEquals((3 * π - 26) / (8 * π), 0));
}

} // namespace numerics