Merge pull request #4407 from pleroy/DownsamplingBenchmark

A benchmark for downsampling

Author: pleroy

benchmarks/discrete_trajectory_benchmark.cpp

@@ -3,38 +3,72 @@
 #include <optional>
 #include <vector>
 
+#include "astronomy/time_scales.hpp"
 #include "base/status_utilities.hpp"  // 🧙 For CHECK_OK.
 #include "benchmark/benchmark.h"
 #include "geometry/barycentre_calculator.hpp"
 #include "geometry/frame.hpp"
 #include "geometry/instant.hpp"
 #include "geometry/space.hpp"
 #include "glog/logging.h"
+#include "integrators/methods.hpp"
+#include "integrators/symmetric_linear_multistep_integrator.hpp"
 #include "ksp_plugin/frames.hpp"
 #include "physics/degrees_of_freedom.hpp"
 #include "physics/discrete_trajectory.hpp"
+#include "physics/discrete_trajectory_segment.hpp"
 #include "physics/discrete_trajectory_types.hpp"
+#include "physics/ephemeris.hpp"
+#include "physics/kepler_orbit.hpp"
+#include "physics/massless_body.hpp"
+#include "physics/solar_system.hpp"
 #include "quantities/quantities.hpp"
 #include "quantities/si.hpp"
 #include "testing_utilities/discrete_trajectory_factories.hpp"
+#include "testing_utilities/solar_system_factory.hpp"
 
 namespace principia {
 namespace physics {
 
+using namespace principia::astronomy::_time_scales;
 using namespace principia::geometry::_barycentre_calculator;
 using namespace principia::geometry::_frame;
 using namespace principia::geometry::_instant;
 using namespace principia::geometry::_space;
+using namespace principia::integrators::_methods;
+using namespace principia::integrators::_symmetric_linear_multistep_integrator;
 using namespace principia::ksp_plugin::_frames;
 using namespace principia::physics::_degrees_of_freedom;
 using namespace principia::physics::_discrete_trajectory;
+using namespace principia::physics::_discrete_trajectory_segment;
 using namespace principia::physics::_discrete_trajectory_types;
+using namespace principia::physics::_ephemeris;
+using namespace principia::physics::_kepler_orbit;
+using namespace principia::physics::_massless_body;
+using namespace principia::physics::_solar_system;
 using namespace principia::quantities::_quantities;
 using namespace principia::quantities::_si;
 using namespace principia::testing_utilities::_discrete_trajectory_factories;
+using namespace principia::testing_utilities::_solar_system_factory;
 
 namespace {
 
+Ephemeris<Barycentric>::FixedStepParameters EphemerisParameters() {
+  return Ephemeris<Barycentric>::FixedStepParameters(
+      SymmetricLinearMultistepIntegrator<
+          QuinlanTremaine1990Order12,
+          Ephemeris<Barycentric>::NewtonianMotionEquation>(),
+      /*step=*/10 * Minute);
+}
+
+Ephemeris<Barycentric>::FixedStepParameters HistoryParameters() {
+  return Ephemeris<Barycentric>::FixedStepParameters(
+      SymmetricLinearMultistepIntegrator<
+          Quinlan1999Order8A,
+          Ephemeris<Barycentric>::NewtonianMotionEquation>(),
+      /*step=*/10 * Second);
+}
+
 // Constructs a trajectory by assigning the points in `timeline` to segments
 // defined by `splits`, which must be doubles in [0, 1].
 DiscreteTrajectory<World> MakeTrajectory(Timeline<World> const& timeline,
@@ -396,6 +430,69 @@ void BM_DiscreteTrajectoryEvaluateDegreesOfFreedomInterpolated(
   }
 }
 
+void BM_DiscreteTrajectoryDownsampling(benchmark::State& state) {
+  static DiscreteTrajectory<Barycentric> const* const goes_8_trajectory = []() {
+    SolarSystem<Barycentric> const solar_system(
+        SOLUTION_DIR / "astronomy" / "sol_gravity_model.proto.txt",
+        SOLUTION_DIR / "astronomy" /
+            "sol_initial_state_jd_2451545_000000000.proto.txt",
+        /*ignore_frame=*/true);
+    auto const ephemeris = solar_system.MakeEphemeris(
+        /*accuracy_parameters=*/{/*fitting_tolerance=*/1 * Milli(Metre),
+                                 /*geopotential_tolerance=*/0x1p-24},
+        EphemerisParameters());
+    auto const earth = solar_system.massive_body(
+        *ephemeris, SolarSystemFactory::name(SolarSystemFactory::Earth));
+
+    // Two-line elements for GOES-8:
+    // 1 23051U 94022A   00004.06628221 -.00000243  00000-0  00000-0 0  9630
+    // 2 23051   0.4232  97.7420 0004776 192.8349 121.5613  1.00264613 28364
+    constexpr Instant goes_8_epoch = "JD2451548.56628221"_UT1;
+    KeplerianElements<Barycentric> goes_8_elements;
+    goes_8_elements.inclination = 0.4232 * Degree;
+    goes_8_elements.longitude_of_ascending_node = 97.7420 * Degree;
+    goes_8_elements.eccentricity = 0.0004776;
+    goes_8_elements.argument_of_periapsis = 192.8349 * Degree;
+    goes_8_elements.mean_anomaly = 121.5613 * Degree;
+    goes_8_elements.mean_motion = 1.00264613 * (2 * π * Radian / Day);
+    CHECK_OK(ephemeris->Prolong(goes_8_epoch));
+    KeplerOrbit<Barycentric> const goes_8_orbit(
+        *earth, MasslessBody{}, goes_8_elements, goes_8_epoch);
+
+    // First build a realistic trajectory without downsampling.
+    auto* const goes_8_trajectory  = new DiscreteTrajectory<Barycentric>;
+    CHECK_OK(goes_8_trajectory->Append(
+        goes_8_epoch,
+        ephemeris->trajectory(earth)->EvaluateDegreesOfFreedom(goes_8_epoch) +
+            goes_8_orbit.StateVectors(goes_8_epoch)));
+    auto goes_8_instance =
+        ephemeris->NewInstance({goes_8_trajectory},
+                               Ephemeris<Barycentric>::NoIntrinsicAccelerations,
+                               HistoryParameters());
+    CHECK_OK(ephemeris->FlowWithFixedStep(goes_8_epoch + 100 * Day,
+                                          *goes_8_instance));
+
+    return goes_8_trajectory;
+  }();
+
+  // Now append the same points to a trajectory with downsampling.
+  DiscreteTrajectory<Barycentric> downsampled_trajectory;
+  for (auto _ : state) {
+    downsampled_trajectory.clear();
+    downsampled_trajectory.segments().begin()->SetDownsampling(
+        DiscreteTrajectorySegment<Barycentric>::DownsamplingParameters{
+            .max_dense_intervals = 10'000, .tolerance = 10 * Metre});
+    for (auto const& [t, degrees_of_freedom] : *goes_8_trajectory) {
+      CHECK_OK(downsampled_trajectory.Append(t, degrees_of_freedom));
+    }
+  }
+  state.SetLabel((std::stringstream()
+                  << goes_8_trajectory->size()
+                  << " points before downsampling, "
+                  << downsampled_trajectory.size() << " after")
+                     .str());
+}
+
 BENCHMARK(BM_DiscreteTrajectoryFront);
 BENCHMARK(BM_DiscreteTrajectoryFrontEmpty);
 BENCHMARK(BM_DiscreteTrajectoryBack);
@@ -419,6 +516,7 @@ BENCHMARK(BM_DiscreteTrajectoryFind)->Range(8, 1024);
 BENCHMARK(BM_DiscreteTrajectoryLowerBound)->Range(8, 1024);
 BENCHMARK(BM_DiscreteTrajectoryEvaluateDegreesOfFreedomExact);
 BENCHMARK(BM_DiscreteTrajectoryEvaluateDegreesOfFreedomInterpolated);
+BENCHMARK(BM_DiscreteTrajectoryDownsampling);
 
 }  // namespace physics
 }  // namespace principia

benchmarks/planetarium_benchmark.cpp

@@ -302,8 +302,7 @@ class Satellites {
   DownsamplingParameters() {
     return DiscreteTrajectorySegment<Barycentric>::DownsamplingParameters{
         .max_dense_intervals = 10'000,
-        .tolerance = 10 * Metre,
-    };
+        .tolerance = 10 * Metre};
   }
 
   not_null<std::unique_ptr<SolarSystem<Barycentric>>> const solar_system_;