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

Commits on Sep 6, 2020

  1. Copy the full SHA
    aea44e6 View commit details
  2. Copy the full SHA
    71a40c6 View commit details
  3. Merge pull request #2711 from pleroy/MultioriginPWPS

    Add support for multiorigin piecewise Poisson series
    pleroy authored Sep 6, 2020
    Copy the full SHA
    2c25b4b View commit details
Showing with 80 additions and 10 deletions.
  1. +1 −1 numerics/poisson_series.hpp
  2. +28 −9 numerics/poisson_series_body.hpp
  3. +51 −0 numerics/poisson_series_test.cpp
2 changes: 1 addition & 1 deletion numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -265,7 +265,7 @@ class PiecewisePoissonSeries {
Instant t_min() const;
Instant t_max() const;

// t must be in the interval [t_min, t_max[.
// 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>
37 changes: 28 additions & 9 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -529,7 +529,6 @@ void PiecewisePoissonSeries<Value, degree_, Evaluator>::Append(
Series const& series) {
CHECK_LT(Time{}, interval.measure());
CHECK_EQ(bounds_.back(), interval.min);
CHECK_EQ(series.origin(), series_.front().origin());
bounds_.push_back(interval.max);
series_.push_back(series);
}
@@ -550,6 +549,10 @@ template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
Value PiecewisePoissonSeries<Value, degree_, Evaluator>::operator()(
Instant const& t) const {
if (t == bounds_.back()) {
return series_.back()(t);
}

// If t is an element of bounds_, the returned iterator points to the next
// element. Otherwise it points to the upper bound of the interval to which
// t belongs.
@@ -651,6 +654,12 @@ PiecewisePoissonSeries<Quotient<Value, Scalar>, degree_, Evaluator> operator/(
return Result(left.bounds_, series);
}

// 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>
@@ -660,7 +669,8 @@ operator+(PoissonSeries<Value, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<Value, std::max(ldegree_, rdegree_), Evaluator>;
std::vector<typename Result::Series> series;
for (int i = 0; i < right.series_.size(); ++i) {
series.push_back(left + right.series_[i]);
Instant const origin = right.series_[i].origin();
series.push_back(left.AtOrigin(origin) + right.series_[i]);
}
return Result(right.bounds_, series);
}
@@ -674,7 +684,8 @@ operator+(PiecewisePoissonSeries<Value, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<Value, std::max(ldegree_, rdegree_), Evaluator>;
std::vector<typename Result::Series> series;
for (int i = 0; i < left.series_.size(); ++i) {
series.push_back(left.series_[i] + right);
Instant const origin = left.series_[i].origin();
series.push_back(left.series_[i] + right.AtOrigin(origin));
}
return Result(left.bounds_, series);
}
@@ -688,7 +699,8 @@ operator-(PoissonSeries<Value, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<Value, std::max(ldegree_, rdegree_), Evaluator>;
std::vector<typename Result::Series> series;
for (int i = 0; i < right.series_.size(); ++i) {
series.push_back(left - right.series_[i]);
Instant const origin = right.series_[i].origin();
series.push_back(left.AtOrigin(origin) - right.series_[i]);
}
return Result(right.bounds_, series);
}
@@ -702,7 +714,8 @@ operator-(PiecewisePoissonSeries<Value, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<Value, std::max(ldegree_, rdegree_), Evaluator>;
std::vector<typename Result::Series> series;
for (int i = 0; i < left.series_.size(); ++i) {
series.push_back(left.series_[i] - right);
Instant const origin = left.series_[i].origin();
series.push_back(left.series_[i] - right.AtOrigin(origin));
}
return Result(left.bounds_, series);
}
@@ -717,7 +730,8 @@ operator*(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
Evaluator>;
std::vector<typename Result::Series> series;
for (int i = 0; i < right.series_.size(); ++i) {
series.push_back(left * right.series_[i]);
Instant const origin = right.series_[i].origin();
series.push_back(left.AtOrigin(origin) * right.series_[i]);
}
return Result(right.bounds_, series);
}
@@ -732,7 +746,8 @@ operator*(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
Evaluator>;
std::vector<typename Result::Series> series;
for (int i = 0; i < left.series_.size(); ++i) {
series.push_back(left.series_[i] * right);
Instant const origin = left.series_[i].origin();
series.push_back(left.series_[i] * right.AtOrigin(origin));
}
return Result(left.bounds_, series);
}
@@ -760,8 +775,10 @@ 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, right.series_[i]) * weight;
PointwiseInnerProduct(left.AtOrigin(origin), right.series_[i]) *
weight.AtOrigin(origin);
auto const primitive = integrand.Primitive();
result += primitive(right.bounds_[i + 1]) - primitive(right.bounds_[i]);
}
@@ -791,8 +808,10 @@ 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) * weight;
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]);
}
51 changes: 51 additions & 0 deletions numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -290,6 +290,7 @@ TEST_F(PiecewisePoissonSeriesTest, Evaluate) {
EXPECT_THAT(pp_(t0_ + 1 * (1 + ε) * Second), VanishesBefore(1, 3));
EXPECT_THAT(pp_(t0_ + 1.5 * Second), AlmostEquals(-Sqrt(0.5), 1));
EXPECT_THAT(pp_(t0_ + 2 * (1 - ε / 2) * Second), AlmostEquals(-1, 0));
EXPECT_THAT(pp_(t0_ + 2 * Second), AlmostEquals(-1, 0));
}

TEST_F(PiecewisePoissonSeriesTest, VectorSpace) {
@@ -359,6 +360,46 @@ TEST_F(PiecewisePoissonSeriesTest, Action) {
}
}

TEST_F(PiecewisePoissonSeriesTest, ActionMultiorigin) {
auto const p = p_.AtOrigin(t0_ + 2 * Second);
{
auto const s1 = p + pp_;
auto const s2 = pp_ + p;
EXPECT_THAT(s1(t0_ + 0.5 * Second),
AlmostEquals((10 - 3 * Sqrt(2)) / 4, 1));
EXPECT_THAT(s1(t0_ + 1.5 * Second),
AlmostEquals((6 + Sqrt(2)) / 4, 0));
EXPECT_THAT(s2(t0_ + 0.5 * Second),
AlmostEquals((10 - 3 * Sqrt(2)) / 4, 1));
EXPECT_THAT(s2(t0_ + 1.5 * Second),
AlmostEquals((6 + Sqrt(2)) / 4, 0));
}
{
auto const d1 = p - pp_;
auto const d2 = pp_ - p;
EXPECT_THAT(d1(t0_ + 0.5 * Second),
AlmostEquals((2 + Sqrt(2)) / 4, 2));
EXPECT_THAT(d1(t0_ + 1.5 * Second),
AlmostEquals((6 + 5 * Sqrt(2)) / 4, 0));
EXPECT_THAT(d2(t0_ + 0.5 * Second),
AlmostEquals((-2 - Sqrt(2)) / 4, 2));
EXPECT_THAT(d2(t0_ + 1.5 * Second),
AlmostEquals((-6 - 5 * Sqrt(2)) / 4, 0));
}
{
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));
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));
EXPECT_THAT(p2(t0_ + 1.5 * Second),
AlmostEquals((-3 - 3 * Sqrt(2)) / 4, 1));
}
}

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

TEST_F(PiecewisePoissonSeriesTest, DotMultiorigin) {
auto const p = p_.AtOrigin(t0_ + 2 * Second);
double const d1 = Dot(
pp_, p, apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second));
double const d2 = Dot(
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));
}

} // namespace numerics
} // namespace principia