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

Commits on Mar 7, 2020

  1. A new test.

    pleroy committed Mar 7, 2020
    Copy the full SHA
    61d14b9 View commit details
  2. Mathematica support.

    pleroy committed Mar 7, 2020
    Copy the full SHA
    7b8b581 View commit details
  3. Drop units in Mathematica.

    pleroy committed Mar 7, 2020
    Copy the full SHA
    0ceec82 View commit details
  4. A complete test.

    pleroy committed Mar 7, 2020
    Copy the full SHA
    7736f6f View commit details
  5. Lint and tolerance.

    pleroy committed Mar 7, 2020
    Copy the full SHA
    0a218ef View commit details
  6. Merge pull request #2484 from pleroy/2400BaselineTest

    A test of the downsampling/serialization/compression of trajectories
    pleroy authored Mar 7, 2020
    Copy the full SHA
    3256aae View commit details
Showing with 132 additions and 0 deletions.
  1. +17 −0 mathematica/mathematica.hpp
  2. +41 −0 mathematica/mathematica_body.hpp
  3. +74 −0 physics/ephemeris_test.cpp
17 changes: 17 additions & 0 deletions mathematica/mathematica.hpp
Original file line number Diff line number Diff line change
@@ -10,6 +10,8 @@
#include "geometry/point.hpp"
#include "geometry/r3_element.hpp"
#include "numerics/fixed_arrays.hpp"
#include "physics/degrees_of_freedom.hpp"
#include "physics/discrete_trajectory.hpp"
#include "quantities/quantities.hpp"

namespace principia {
@@ -20,6 +22,8 @@ using geometry::Point;
using geometry::R3Element;
using geometry::Vector;
using numerics::FixedVector;
using physics::DegreesOfFreedom;
using physics::DiscreteTrajectory;
using quantities::Quantity;
using quantities::Quotient;

@@ -40,6 +44,9 @@ std::string PlottableDataset(std::vector<T> const& x, std::vector<U> const& y);
template<typename T>
std::string ToMathematica(std::vector<T> const& list);

template<typename It>
std::string ToMathematica(It begin, It end);

std::string ToMathematica(double const& real);

template<typename T, int size>
@@ -57,9 +64,17 @@ std::string ToMathematica(Vector<S, F> const& vector);
template<typename V>
std::string ToMathematica(Point<V> const& point);

template<typename F>
std::string ToMathematica(DegreesOfFreedom<F> const& degrees_of_freedom);

template<typename... Types>
std::string ToMathematica(std::tuple<Types...> const& tuple);

template<typename R,
typename = std::void_t<decltype(std::declval<R>().time)>,
typename = std::void_t<decltype(std::declval<R>().degrees_of_freedom)>>
std::string ToMathematica(R ref);

std::string ToMathematica(
astronomy::OrbitalElements::EquinoctialElements const& elements);

@@ -70,6 +85,8 @@ std::string ToMathematica(std::string const& str);
// TODO(egg): escape things properly.
std::string Escape(std::string const& str);

// TODO(phl): This doesn't work well for complex structures like orbital
// elements or trajectory iterators. Surely we can do better.
template<typename T>
struct RemoveUnit;

41 changes: 41 additions & 0 deletions mathematica/mathematica_body.hpp
Original file line number Diff line number Diff line change
@@ -80,6 +80,15 @@ std::string ToMathematica(std::vector<T> const& list) {
return Apply("List", expressions);
}

template<typename It>
std::string ToMathematica(It const begin, It const end) {
std::vector<std::string> expressions;
for (auto it = begin; it != end; ++it) {
expressions.emplace_back(ToMathematica(*it));
}
return Apply("List", expressions);
}

inline std::string ToMathematica(double const& real) {
if (std::isinf(real)) {
if (real > 0.0) {
@@ -138,6 +147,16 @@ std::string ToMathematica(Point<V> const & point) {
return ToMathematica(point - Point<V>());
}

template<typename F>
std::string ToMathematica(DegreesOfFreedom<F> const& degrees_of_freedom) {
return Apply(
"List",
std::vector<std::string>{
ToMathematica(ExpressIn(Metre, degrees_of_freedom.position())),
ToMathematica(
ExpressIn(Metre / Second, degrees_of_freedom.velocity()))});
}

template<typename... Types>
std::string ToMathematica(std::tuple<Types...> const& tuple) {
std::vector<std::string> expressions;
@@ -147,6 +166,14 @@ std::string ToMathematica(std::tuple<Types...> const& tuple) {
return Apply("List", expressions);
}

template<typename R, typename, typename>
std::string ToMathematica(R const ref) {
return Apply(
"List",
std::vector<std::string>{ToMathematica(ExpressIn(Second, ref.time)),
ToMathematica(ref.degrees_of_freedom)});
}

inline std::string ToMathematica(
astronomy::OrbitalElements::EquinoctialElements const& elements) {
return ToMathematica(std::make_tuple((elements.t - J2000) / Second,
@@ -185,6 +212,12 @@ struct RemoveUnit<Vector<T, F>> {
using Unitless = Vector<typename RemoveUnit<T>::Unitless, F>;
};

template<typename V>
struct RemoveUnit<Point<V>> {
using Unit = typename RemoveUnit<V>::Unit;
using Unitless = Point<typename RemoveUnit<V>::Unitless>;
};

template<typename T>
struct RemoveUnit<std::vector<T>> {
using Unit = typename RemoveUnit<T>::Unit;
@@ -198,6 +231,14 @@ typename RemoveUnit<T>::Unitless ExpressIn(
return value / unit;
}

template<typename V>
typename RemoveUnit<Point<V>>::Unitless ExpressIn(
typename RemoveUnit<Point<V>>::Unit const& unit,
Point<V> const& value) {
return (value - Point<V>{}) / unit +
typename RemoveUnit<Point<V>>::Unitless{};
}

template<typename T>
typename RemoveUnit<std::vector<T>>::Unitless ExpressIn(
typename RemoveUnit<std::vector<T>>::Unit const& unit,
74 changes: 74 additions & 0 deletions physics/ephemeris_test.cpp
Original file line number Diff line number Diff line change
@@ -5,19 +5,23 @@
#include <map>
#include <optional>
#include <set>
#include <string>
#include <vector>

#include "astronomy/frames.hpp"
#include "base/file.hpp"
#include "base/macros.hpp"
#include "geometry/barycentre_calculator.hpp"
#include "geometry/frame.hpp"
#include "gipfeli/gipfeli.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "integrators/embedded_explicit_generalized_runge_kutta_nyström_integrator.hpp"
#include "integrators/embedded_explicit_runge_kutta_nyström_integrator.hpp"
#include "integrators/methods.hpp"
#include "integrators/symmetric_linear_multistep_integrator.hpp"
#include "integrators/symplectic_runge_kutta_nyström_integrator.hpp"
#include "mathematica/mathematica.hpp"
#include "physics/kepler_orbit.hpp"
#include "physics/massive_body.hpp"
#include "physics/oblate_body.hpp"
@@ -44,6 +48,7 @@ namespace internal_ephemeris {

using astronomy::ICRS;
using base::not_null;
using base::OFStream;
using geometry::Barycentre;
using geometry::AngularVelocity;
using geometry::Displacement;
@@ -59,6 +64,7 @@ using integrators::methods::Fine1987RKNG34;
using integrators::methods::McLachlanAtela1992Order4Optimal;
using integrators::methods::McLachlanAtela1992Order5Optimal;
using integrators::methods::Quinlan1999Order8A;
using integrators::methods::QuinlanTremaine1990Order12;
using quantities::Abs;
using quantities::ArcTan;
using quantities::Area;
@@ -1057,6 +1063,74 @@ TEST_P(EphemerisTest, ComputeApsidesContinuousTrajectory) {
}
}

#if !defined(_DEBUG)
// This trajectory is similar to the second trajectory in the first save in
// #2400. It exhibits oscillations with a period close to 5600 s and its
// downsampling period alternates between 120 and 130 s.
TEST(EphemerisTestNoFixture, DiscreteTrajectoryCompression) {
SolarSystem<ICRS> solar_system(
SOLUTION_DIR / "astronomy" / "sol_gravity_model.proto.txt",
SOLUTION_DIR / "astronomy" /
"sol_initial_state_jd_2433282_500000000.proto.txt");

auto ephemeris = solar_system.MakeEphemeris(
/*accuracy_parameters=*/{/*fitting_tolerance=*/1 * Milli(Metre),
/*geopotential_tolerance=*/0x1p-24},
/*fixed_step_parameters=*/{
SymmetricLinearMultistepIntegrator<QuinlanTremaine1990Order12,
Position<ICRS>>(),
/*step=*/10 * Minute});

Instant const t0 = Instant() - 1.323698948906726e9 * Second;
Instant const t1 = Instant() - 1.323595217786725e9 * Second;
Position<ICRS> const q0 =
Position<ICRS>{} + Displacement<ICRS>({-7.461169719467950e10 * Metre,
1.165327644733623e11 * Metre,
5.049935298178532e10 * Metre});
Velocity<ICRS> const p0({-30169.49165384562 * Metre / Second,
-11880.03394238412 * Metre / Second,
200.2551546021678 * Metre / Second});

MasslessBody probe;
DiscreteTrajectory<ICRS> trajectory;
trajectory.SetDownsampling(/*max_dense_intervals=*/10'000,
/*tolerance=*/10 * Metre);
trajectory.Append(t0, DegreesOfFreedom<ICRS>(q0, p0));

auto instance = ephemeris->NewInstance(
{&trajectory},
Ephemeris<ICRS>::NoIntrinsicAccelerations,
Ephemeris<ICRS>::FixedStepParameters(
SymmetricLinearMultistepIntegrator<Quinlan1999Order8A,
Position<ICRS>>(),
10 * Second));
EXPECT_OK(ephemeris->FlowWithFixedStep(t1, *instance));
EXPECT_EQ(1162, trajectory.Size());

serialization::DiscreteTrajectory message;
trajectory.WriteToMessage(&message, {});
std::string uncompressed;
message.SerializePartialToString(&uncompressed);
EXPECT_EQ(178'982, uncompressed.size());

std::string compressed;
auto compressor = google::compression::NewGipfeliCompressor();
compressor->Compress(uncompressed, &compressed);

// We want a change detector, but the actual compressed size varies depending
// on the exact numerical values, and therefore on the mathematical library.
EXPECT_LE(69'883, compressed.size());
EXPECT_GE(69'910, compressed.size());

{
OFStream file(TEMP_DIR / "discrete_trajectory_compression.generated.wl");
file << mathematica::Assign(
"trajectory",
mathematica::ToMathematica(trajectory.begin(), trajectory.end()));
}
}
#endif

INSTANTIATE_TEST_CASE_P(
AllEphemerisTests,
EphemerisTest,