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

Commits on Oct 7, 2020

  1. Copy the full SHA
    4510ad3 View commit details

Commits on Oct 8, 2020

  1. Merge pull request #2755 from pleroy/Serialization

    Serialization of Poisson series
    pleroy authored Oct 8, 2020
    Copy the full SHA
    7240b46 View commit details
Showing with 74 additions and 43 deletions.
  1. +6 −0 numerics/poisson_series.hpp
  2. +35 −0 numerics/poisson_series_body.hpp
  3. +22 −43 numerics/poisson_series_test.cpp
  4. +11 −0 serialization/numerics.proto
6 changes: 6 additions & 0 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -9,6 +9,7 @@
#include <utility>
#include <vector>

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

using base::not_null;
using geometry::Complexification;
using geometry::Hilbert;
using geometry::Instant;
@@ -122,6 +124,10 @@ class PoissonSeries {
template<typename V, int d, template<typename, typename, int> class E>
PoissonSeries& operator-=(PoissonSeries<V, d, E> const& right);

void WriteToMessage(not_null<serialization::PoissonSeries*> message) const;
static PoissonSeries ReadFromMessage(
serialization::PoissonSeries const& message);

private:
// Similar to the public constructor, but passing by copy allows moves, which
// is useful for internal algorithms.
35 changes: 35 additions & 0 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -229,6 +229,41 @@ PoissonSeries<Value, degree_, Evaluator>::Integrate(Instant const& t1,
return result;
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
void PoissonSeries<Value, degree_, Evaluator>::WriteToMessage(
not_null<serialization::PoissonSeries*> const message) const {
aperiodic_.WriteToMessage(message->mutable_aperiodic());
for (auto const& [ω, polynomials] : periodic_) {
auto* const polynomials_and_angular_frequency = message->add_periodic();
ω.WriteToMessage(
polynomials_and_angular_frequency->mutable_angular_frequency());
polynomials.sin.WriteToMessage(
polynomials_and_angular_frequency->mutable_sin());
polynomials.cos.WriteToMessage(
polynomials_and_angular_frequency->mutable_cos());
}
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, degree_, Evaluator>
PoissonSeries<Value, degree_, Evaluator>::ReadFromMessage(
serialization::PoissonSeries const& message) {
auto const aperiodic = Polynomial::ReadFromMessage(message.aperiodic());
PolynomialsByAngularFrequency periodic;
for (auto const& polynomial_and_angular_frequency : message.periodic()) {
auto const ω = AngularFrequency::ReadFromMessage(
polynomial_and_angular_frequency.angular_frequency());
auto const sin =
Polynomial::ReadFromMessage(polynomial_and_angular_frequency.sin());
auto const cos =
Polynomial::ReadFromMessage(polynomial_and_angular_frequency.cos());
periodic.emplace_back(ω, Polynomials{/*sin=*/sin, /*cos=*/cos});
}
return PoissonSeries(aperiodic, periodic);
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, degree_, Evaluator>::PoissonSeries(
65 changes: 22 additions & 43 deletions numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -16,8 +16,10 @@
#include "quantities/elementary_functions.hpp"
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"
#include "serialization/numerics.pb.h"
#include "testing_utilities/almost_equals.hpp"
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/matchers.hpp"
#include "testing_utilities/numerics_matchers.hpp"
#include "testing_utilities/vanishes_before.hpp"

@@ -42,6 +44,7 @@ using quantities::si::Metre;
using quantities::si::Radian;
using quantities::si::Second;
using testing_utilities::AlmostEquals;
using testing_utilities::EqualsProto;
using testing_utilities::IsNear;
using testing_utilities::VanishesBefore;
using testing_utilities::RelativeErrorFrom;
@@ -312,53 +315,29 @@ TEST_F(PoissonSeriesTest, InnerProduct) {
AlmostEquals(-381.25522770148542400, 71));
}

TEST_F(PoissonSeriesTest, DotConditioning) {
using Degree7 = PoissonSeries<double, 7, HornerEvaluator>;
using Degree2 = PoissonSeries<double, 2, HornerEvaluator>;

Instant const t_min = t0_;
Instant const t_max = t0_ + 4800 * Second;
Instant const t_mid = t0_ + 2400 * Second;

Degree7 s7(Degree7::Polynomial({-8.752190840128326e8,
-6265.007683216121 / Second,
-0.3289504016189549 / Pow<2>(Second),
+1.973941298531457e-6 / Pow<3>(Second),
+4.757589125055360e-11 / Pow<4>(Second),
-1.668764721235290e-16 / Pow<5>(Second),
-2.971788090416876e-21 / Pow<6>(Second),
+5.893739999926978e-27 / Pow<7>(Second)},
t_mid),
{{}});

Degree2 s2(Degree2::Polynomial({0, 0 / Second, 1 / Second / Second}, t_min),
{{}});

Instant const origin = s7.origin();
auto const integrand =
PointwiseInnerProduct(s7, s2.AtOrigin(origin)) *
apodization::Hann<HornerEvaluator>(t_min, t0_ + 7891200 * Second)
.AtOrigin(origin);
auto const primitive = integrand.Primitive();

// Exact value is -7.15802e13
LOG(ERROR) << primitive(t_max) << " " << primitive(t_min) << " "
<< primitive(t_max) - primitive(t_min);

auto const integral = integrand.Integrate(t_min, t_max);
LOG(ERROR) << integral;

for (int n = 1; n < 100'000; n *= 10) {
auto const better_integral =
quadrature::Midpoint(integrand, t_min, t_max, n);
LOG(ERROR) << better_integral;
}
}

TEST_F(PoissonSeriesTest, Output) {
LOG(ERROR) << *pa_;
}

TEST_F(PoissonSeriesTest, Serialization) {
serialization::PoissonSeries message;
pa_->WriteToMessage(&message);
EXPECT_TRUE(message.has_aperiodic());
EXPECT_EQ(2, message.periodic_size());

auto const poisson_series_read = Degree1::ReadFromMessage(message);
EXPECT_THAT((*pa_)(t0_ + 1 * Second),
AlmostEquals(poisson_series_read(t0_ + 1 * Second), 0));
EXPECT_THAT((*pa_)(t0_ + 2 * Second),
AlmostEquals(poisson_series_read(t0_ + 2 * Second), 0));
EXPECT_THAT((*pa_)(t0_ + 3 * Second),
AlmostEquals(poisson_series_read(t0_ + 3 * Second), 0));

serialization::PoissonSeries message2;
poisson_series_read.WriteToMessage(&message2);
EXPECT_THAT(message2, EqualsProto(message));
}

class PiecewisePoissonSeriesTest : public ::testing::Test {
protected:
using Degree0 = PiecewisePoissonSeries<double, 0, HornerEvaluator>;
11 changes: 11 additions & 0 deletions serialization/numerics.proto
Original file line number Diff line number Diff line change
@@ -42,6 +42,17 @@ message DoublePrecision {
required Error error = 2;
}

message PoissonSeries {
message PolynomialsAndAngularFrequency {
required Quantity angular_frequency = 1;
required Polynomial sin = 2;
required Polynomial cos = 3;
}
// The origin is not serialized since it's part of the polynomials.
required Polynomial aperiodic = 1;
repeated PolynomialsAndAngularFrequency periodic = 2;
}

// Polynomials follow.

message Polynomial {