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

Commits on Oct 16, 2020

  1. Start adding addend_.

    pleroy committed Oct 16, 2020
    Copy the full SHA
    3f6a5d8 View commit details
  2. Copy the full SHA
    748ab04 View commit details

Commits on Oct 17, 2020

  1. Copy the full SHA
    78028a8 View commit details
  2. Copy the full SHA
    073f9d7 View commit details
  3. Merge pull request #2766 from pleroy/NoAtOrigin

    Change PiecewisePoissonSeries to avoid origin changes on additions
    pleroy authored Oct 17, 2020
    Copy the full SHA
    aa789a2 View commit details
Showing with 173 additions and 72 deletions.
  1. +22 −9 numerics/poisson_series.hpp
  2. +131 −55 numerics/poisson_series_body.hpp
  3. +19 −8 numerics/poisson_series_test.cpp
  4. +1 −0 serialization/numerics.proto
31 changes: 22 additions & 9 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -99,6 +99,13 @@ class PoissonSeries {
PoissonSeries(Polynomial const& aperiodic,
PolynomialsByAngularFrequency const& periodic);

// A Poisson series may be explicitly converted to a higher degree (possibly
// with a different evaluator).
template<int higher_degree_,
template<typename, typename, int> class HigherEvaluator>
explicit operator
PoissonSeries<Value, higher_degree_, HigherEvaluator>() const;

Instant const& origin() const;

Value operator()(Instant const& t) const;
@@ -119,10 +126,10 @@ class PoissonSeries {
Instant const& t_min,
Instant const& t_max) const;

template<typename V, int d, template<typename, typename, int> class E>
PoissonSeries& operator+=(PoissonSeries<V, d, E> const& right);
template<typename V, int d, template<typename, typename, int> class E>
PoissonSeries& operator-=(PoissonSeries<V, d, E> const& right);
template<int d>
PoissonSeries& operator+=(PoissonSeries<Value, d, Evaluator> const& right);
template<int d>
PoissonSeries& operator-=(PoissonSeries<Value, d, Evaluator> const& right);

void WriteToMessage(not_null<serialization::PoissonSeries*> message) const;
static PoissonSeries ReadFromMessage(
@@ -339,10 +346,12 @@ class PiecewisePoissonSeries {
// |*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>
PiecewisePoissonSeries& operator-=(PoissonSeries<V, d, E> const& right);
template<int d>
PiecewisePoissonSeries& operator+=(
PoissonSeries<Value, d, Evaluator> const& right);
template<int d>
PiecewisePoissonSeries& operator-=(
PoissonSeries<Value, d, Evaluator> const& right);

void WriteToMessage(
not_null<serialization::PiecewisePoissonSeries*> message) const;
@@ -351,10 +360,14 @@ class PiecewisePoissonSeries {

private:
PiecewisePoissonSeries(std::vector<Instant> const& bounds,
std::vector<Series> const& series);
std::vector<Series> const& series,
std::optional<Series> const& addend);

Value EvaluateAddend(Instant const& t) const;

std::vector<Instant> bounds_;
std::vector<Series> series_;
std::optional<Series> addend_;

template<typename V, int r, template<typename, typename, int> class E>
PiecewisePoissonSeries<V, r, E>
186 changes: 131 additions & 55 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -343,6 +343,29 @@ PoissonSeries<Value, degree_, Evaluator>::PoissonSeries(
aperiodic_(std::move(aperiodic)),
periodic_(std::move(periodic)) {}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<int higher_degree_,
template<typename, typename, int> class HigherEvaluator>
PoissonSeries<Value, degree_, Evaluator>::
operator PoissonSeries<Value, higher_degree_, HigherEvaluator>() const {
static_assert(degree_ <= higher_degree_);
using Result = PoissonSeries<Value, higher_degree_, HigherEvaluator>;
auto aperiodic = typename Result::Polynomial(aperiodic_);
typename Result::PolynomialsByAngularFrequency periodic;
periodic.reserve(periodic_.size());
for (auto const& [ω, polynomials] : periodic_) {
periodic.emplace_back(
ω,
typename Result::Polynomials{
/*sin=*/typename Result::Polynomial(polynomials.sin),
/*cos=*/typename Result::Polynomial(polynomials.cos)});
}
return Result(typename Result::TrustedPrivateConstructor{},
std::move(aperiodic),
std::move(periodic));
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<int wdegree_>
@@ -377,20 +400,22 @@ PoissonSeries<Value, degree_, Evaluator>::Norm(

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<typename V, int d, template<typename, typename, int> class E>
template<int d>
PoissonSeries<Value, degree_, Evaluator>&
PoissonSeries<Value, degree_, Evaluator>::operator+=(
PoissonSeries<V, d, E> const& right) {
PoissonSeries<Value, d, Evaluator> const& right) {
static_assert(d <= degree);
*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>
template<int d>
PoissonSeries<Value, degree_, Evaluator>&
PoissonSeries<Value, degree_, Evaluator>::operator-=(
PoissonSeries<V, d, E> const& right) {
PoissonSeries<Value, d, Evaluator> const& right) {
static_assert(d <= degree);
*this = *this - right;
return *this;
}
@@ -732,8 +757,9 @@ template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
Value PiecewisePoissonSeries<Value, degree_, Evaluator>::operator()(
Instant const& t) const {
Value const addend = EvaluateAddend(t);
if (t == bounds_.back()) {
return series_.back()(t);
return series_.back()(t) + addend;
}

// If t is an element of bounds_, the returned iterator points to the next
@@ -745,7 +771,7 @@ Value PiecewisePoissonSeries<Value, degree_, Evaluator>::operator()(
<< bounds_.front() << " .. " << bounds_.back();
CHECK(it != bounds_.cend())
<< t << " is outside of " << bounds_.front() << " .. " << bounds_.back();
return series_[it - bounds_.cbegin() - 1](t);
return series_[it - bounds_.cbegin() - 1](t) + addend;
}

template<typename Value,
@@ -762,10 +788,11 @@ auto PiecewisePoissonSeries<Value, degree_, Evaluator>::FourierTransform() const
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, ω](
[this, &f = series_[k], t0, ω](
Instant const& t) -> Complexification<Value> {
return f(t) * Complexification<double>{Cos(ω * (t - t0)),
-Sin(ω * (t - t0))};
return (f(t) + EvaluateAddend(t)) *
Complexification<double>{Cos(ω * (t - t0)),
-Sin(ω * (t - t0))};
},
bounds_[k],
bounds_[k + 1]);
@@ -776,26 +803,30 @@ auto PiecewisePoissonSeries<Value, degree_, Evaluator>::FourierTransform() const

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<typename V, int d, template<typename, typename, int> class E>
template<int d>
PiecewisePoissonSeries<Value, degree_, Evaluator>&
PiecewisePoissonSeries<Value, degree_, Evaluator>::operator+=(
PoissonSeries<V, d, E> const& right) {
for (int i = 0; i < series_.size(); ++i) {
auto& series = series_[i];
series += right.AtOrigin(series.origin());
PoissonSeries<Value, d, Evaluator> const& right) {
static_assert(d <= degree);
if (addend_.has_value()) {
addend_.value() += right;
} else {
addend_ = Series(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>
template<int d>
PiecewisePoissonSeries<Value, degree_, Evaluator>&
PiecewisePoissonSeries<Value, degree_, Evaluator>::operator-=(
PoissonSeries<V, d, E> const& right) {
for (int i = 0; i < series_.size(); ++i) {
auto& series = series_[i];
series -= right.AtOrigin(series.origin());
PoissonSeries<Value, d, Evaluator> const& right) {
static_assert(d <= degree);
if (addend_.has_value()) {
addend_.value() -= right;
} else {
addend_ = Series(-right);
}
return *this;
}
@@ -810,6 +841,9 @@ void PiecewisePoissonSeries<Value, degree_, Evaluator>::WriteToMessage(
for (Series const& series : series_) {
series.WriteToMessage(message->add_series());
}
if (addend_.has_value()) {
addend_->WriteToMessage(message->mutable_addend());
}
}

template<typename Value, int degree_,
@@ -830,16 +864,28 @@ PiecewisePoissonSeries<Value, degree_, Evaluator>::ReadFromMessage(
Instant::ReadFromMessage(message.bounds(i + 1))};
series.Append(interval, Series::ReadFromMessage(message.series(i)));
}
if (message.has_addend()) {
series.addend_ = Series::ReadFromMessage(message.addend());
}
return series;
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PiecewisePoissonSeries<Value, degree_, Evaluator>::PiecewisePoissonSeries(
std::vector<Instant> const& bounds,
std::vector<PoissonSeries<Value, degree_, Evaluator>> const& series)
std::vector<PoissonSeries<Value, degree_, Evaluator>> const& series,
std::optional<Series> const& addend)
: bounds_(bounds),
series_(series) {}
series_(series),
addend_(addend) {}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
Value PiecewisePoissonSeries<Value, degree_, Evaluator>::EvaluateAddend(
Instant const& t) const {
return addend_.has_value() ? addend_.value()(t) : Value{};
}

template<typename Value, int rdegree_,
template<typename, typename, int> class Evaluator>
@@ -858,7 +904,11 @@ PiecewisePoissonSeries<Value, rdegree_, Evaluator> operator-(
for (int i = 0; i < right.series_.size(); ++i) {
series.push_back(-right.series_[i]);
}
return Result(right.bounds_, series);
std::optional<typename Result::Series> addend;
if (right.addend_.has_value()) {
addend = -right.addend_.value();
}
return Result(right.bounds_, series, addend);
}

template<typename Scalar, typename Value, int degree_,
@@ -873,7 +923,11 @@ PiecewisePoissonSeries<Product<Scalar, Value>, degree_, Evaluator> operator*(
for (int i = 0; i < right.series_.size(); ++i) {
series.push_back(left * right.series_[i]);
}
return Result(right.bounds_, series);
std::optional<typename Result::Series> addend;
if (right.addend_.has_value()) {
addend = left * right.addend_.value();
}
return Result(right.bounds_, series, addend);
}

template<typename Scalar, typename Value, int degree_,
@@ -888,7 +942,11 @@ PiecewisePoissonSeries<Product<Value, Scalar>, degree_, Evaluator> operator*(
for (int i = 0; i < left.series_.size(); ++i) {
series.push_back(left.series_[i] * right);
}
return Result(left.bounds_, series);
std::optional<typename Result::Series> addend;
if (left.addend_.has_value()) {
addend = left.addend_.value() * right;
}
return Result(left.bounds_, series, addend);
}

template<typename Scalar, typename Value, int degree_,
@@ -903,29 +961,27 @@ PiecewisePoissonSeries<Quotient<Value, Scalar>, degree_, Evaluator> operator/(
for (int i = 0; i < left.series_.size(); ++i) {
series.push_back(left.series_[i] / right);
}
return Result(left.bounds_, series);
std::optional<typename Result::Series> addend;
if (left.addend_.has_value()) {
addend = left.addend_.value() / right;
}
return Result(left.bounds_, series, addend);
}

// In practice changing the origin of the piecewise series chunks is horribly
// ill-conditioned, so the code below changes the origin of the (single) Poisson
// series.
// TODO(phl): All these origin changes might be expensive, see if we can factor
// them.

template<typename Value, int ldegree_, int rdegree_,
template<typename, typename, int> class Evaluator>
PiecewisePoissonSeries<Value, std::max(ldegree_, rdegree_), Evaluator>
operator+(PoissonSeries<Value, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<Value, rdegree_, Evaluator> const& right) {
using Result =
PiecewisePoissonSeries<Value, std::max(ldegree_, rdegree_), Evaluator>;
std::vector<typename Result::Series> series;
series.reserve(right.series_.size());
for (int i = 0; i < right.series_.size(); ++i) {
Instant const origin = right.series_[i].origin();
series.push_back(left.AtOrigin(origin) + right.series_[i]);
std::optional<typename Result::Series> addend;
if (right.addend_.has_value()) {
addend = left + right.addend_.value();
} else {
addend = typename Result::Series(left);
}
return Result(right.bounds_, series);
return Result(right.bounds_, right.series_, addend);
}

template<typename Value, int ldegree_, int rdegree_,
@@ -935,13 +991,13 @@ operator+(PiecewisePoissonSeries<Value, ldegree_, Evaluator> const& left,
PoissonSeries<Value, rdegree_, Evaluator> const& right) {
using Result =
PiecewisePoissonSeries<Value, std::max(ldegree_, rdegree_), Evaluator>;
std::vector<typename Result::Series> series;
series.reserve(left.series_.size());
for (int i = 0; i < left.series_.size(); ++i) {
Instant const origin = left.series_[i].origin();
series.push_back(left.series_[i] + right.AtOrigin(origin));
std::optional<typename Result::Series> addend;
if (left.addend_.has_value()) {
addend = left.addend_.value() + right;
} else {
addend = typename Result::Series(right);
}
return Result(left.bounds_, series);
return Result(left.bounds_, left.series_, addend);
}

template<typename Value, int ldegree_, int rdegree_,
@@ -954,10 +1010,15 @@ operator-(PoissonSeries<Value, ldegree_, Evaluator> const& left,
std::vector<typename Result::Series> series;
series.reserve(right.series_.size());
for (int i = 0; i < right.series_.size(); ++i) {
Instant const origin = right.series_[i].origin();
series.push_back(left.AtOrigin(origin) - right.series_[i]);
series.push_back(-right.series_[i]);
}
std::optional<typename Result::Series> addend;
if (right.addend_.has_value()) {
addend = left - right.addend_.value();
} else {
addend = typename Result::Series(left);
}
return Result(right.bounds_, series);
return Result(right.bounds_, series, addend);
}

template<typename Value, int ldegree_, int rdegree_,
@@ -970,10 +1031,15 @@ operator-(PiecewisePoissonSeries<Value, ldegree_, Evaluator> const& left,
std::vector<typename Result::Series> series;
series.reserve(left.series_.size());
for (int i = 0; i < left.series_.size(); ++i) {
Instant const origin = left.series_[i].origin();
series.push_back(left.series_[i] - right.AtOrigin(origin));
series.push_back(typename Result::Series(left.series_[i]));
}
std::optional<typename Result::Series> addend;
if (left.addend_.has_value()) {
addend = left.addend_.value() - right;
} else {
addend = typename Result::Series(-right);
}
return Result(left.bounds_, series);
return Result(left.bounds_, series, addend);
}

template<typename LValue, typename RValue, int ldegree_, int rdegree_,
@@ -990,7 +1056,11 @@ operator*(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
Instant const origin = right.series_[i].origin();
series.push_back(left.AtOrigin(origin) * right.series_[i]);
}
return Result(right.bounds_, series);
std::optional<typename Result::Series> addend;
if (right.addend_.has_value()) {
addend = left * right.addend_.value();
}
return Result(right.bounds_, series, addend);
}

template<typename LValue, typename RValue, int ldegree_, int rdegree_,
@@ -1007,7 +1077,11 @@ operator*(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
Instant const origin = left.series_[i].origin();
series.push_back(left.series_[i] * right.AtOrigin(origin));
}
return Result(left.bounds_, series);
std::optional<typename Result::Series> addend;
if (left.addend_.has_value()) {
addend = left.addend_.value() * right;
}
return Result(left.bounds_, series, addend);
}

template<typename LValue, typename RValue,
@@ -1040,8 +1114,9 @@ InnerProduct(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
Result result{};
for (int i = 0; i < right.series_.size(); ++i) {
auto integrand = [i, &left, &right, &weight](Instant const& t) {
return Hilbert<LValue, RValue>::InnerProduct(left(t) * weight(t),
right.series_[i](t));
return Hilbert<LValue, RValue>::InnerProduct(
left(t) * weight(t),
right.series_[i](t) + right.EvaluateAddend(t));
};
auto const integral = quadrature::GaussLegendre<points>(
integrand, right.bounds_[i], right.bounds_[i + 1]);
@@ -1080,8 +1155,9 @@ InnerProduct(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
Result result{};
for (int i = 0; i < left.series_.size(); ++i) {
auto integrand = [i, &left, &right, &weight](Instant const& t) {
return Hilbert<LValue, RValue>::InnerProduct(left.series_[i](t),
right(t) * weight(t));
return Hilbert<LValue, RValue>::InnerProduct(
left.series_[i](t) + left.EvaluateAddend(t),
right(t) * weight(t));
};
auto const integral = quadrature::GaussLegendre<points>(
integrand, left.bounds_[i], left.bounds_[i + 1]);
27 changes: 19 additions & 8 deletions numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -120,6 +120,15 @@ TEST_F(PoissonSeriesTest, Evaluate) {
32));
}

TEST_F(PoissonSeriesTest, Conversion) {
using Degree3 = PoissonSeries<double, 3, HornerEvaluator>;
Degree3 const pa3 = Degree3(*pa_);
EXPECT_THAT(pa3(t0_ + 1 * Second),
AlmostEquals(3 + 11 * Sin(1 * Radian) + 15 * Cos(1 * Radian) +
27 * Sin(2 * Radian) + 31 * Cos(2 * Radian),
0, 1));
}

TEST_F(PoissonSeriesTest, VectorSpace) {
{
auto const identity = +*pa_;
@@ -469,23 +478,23 @@ TEST_F(PiecewisePoissonSeriesTest, Action) {
EXPECT_THAT(s1(t0_ + 0.5 * Second),
AlmostEquals((10 - 3 * Sqrt(2)) / 4, 0));
EXPECT_THAT(s1(t0_ + 1.5 * Second),
AlmostEquals((6 + Sqrt(2)) / 4, 0));
AlmostEquals((6 + Sqrt(2)) / 4, 1));
EXPECT_THAT(s2(t0_ + 0.5 * Second),
AlmostEquals((10 - 3 * Sqrt(2)) / 4, 0));
EXPECT_THAT(s2(t0_ + 1.5 * Second),
AlmostEquals((6 + Sqrt(2)) / 4, 0));
AlmostEquals((6 + Sqrt(2)) / 4, 1));
}
{
auto const d1 = p_ - pp_;
auto const d2 = pp_ - p_;
EXPECT_THAT(d1(t0_ + 0.5 * Second),
AlmostEquals((2 + Sqrt(2)) / 4, 1));
EXPECT_THAT(d1(t0_ + 1.5 * Second),
AlmostEquals((6 + 5 * Sqrt(2)) / 4, 0));
AlmostEquals((6 + 5 * Sqrt(2)) / 4, 1));
EXPECT_THAT(d2(t0_ + 0.5 * Second),
AlmostEquals((-2 - Sqrt(2)) / 4, 1));
EXPECT_THAT(d2(t0_ + 1.5 * Second),
AlmostEquals((-6 - 5 * Sqrt(2)) / 4, 0));
AlmostEquals((-6 - 5 * Sqrt(2)) / 4, 1));
}
{
auto const p1 = p_ * pp_;
@@ -507,13 +516,13 @@ TEST_F(PiecewisePoissonSeriesTest, ActionMultiorigin) {
auto const s1 = p + pp_;
auto const s2 = pp_ + p;
EXPECT_THAT(s1(t0_ + 0.5 * Second),
AlmostEquals((10 - 3 * Sqrt(2)) / 4, 1));
AlmostEquals((10 - 3 * Sqrt(2)) / 4, 0));
EXPECT_THAT(s1(t0_ + 1.5 * Second),
AlmostEquals((6 + Sqrt(2)) / 4, 0));
AlmostEquals((6 + Sqrt(2)) / 4, 1));
EXPECT_THAT(s2(t0_ + 0.5 * Second),
AlmostEquals((10 - 3 * Sqrt(2)) / 4, 1));
AlmostEquals((10 - 3 * Sqrt(2)) / 4, 0));
EXPECT_THAT(s2(t0_ + 1.5 * Second),
AlmostEquals((6 + Sqrt(2)) / 4, 0));
AlmostEquals((6 + Sqrt(2)) / 4, 1));
}
{
auto const d1 = p - pp_;
@@ -562,9 +571,11 @@ TEST_F(PiecewisePoissonSeriesTest, InnerProductMultiorigin) {

TEST_F(PiecewisePoissonSeriesTest, Serialization) {
serialization::PiecewisePoissonSeries message;
pp_ += p_;
pp_.WriteToMessage(&message);
EXPECT_EQ(3, message.bounds_size());
EXPECT_EQ(2, message.series_size());
EXPECT_TRUE(message.has_addend());

auto const piecewise_poisson_series_read = Degree0::ReadFromMessage(message);
EXPECT_THAT(
1 change: 1 addition & 0 deletions serialization/numerics.proto
Original file line number Diff line number Diff line change
@@ -45,6 +45,7 @@ message DoublePrecision {
message PiecewisePoissonSeries {
repeated Point bounds = 1;
repeated PoissonSeries series = 2;
optional PoissonSeries addend = 3;
}

message PoissonSeries {