mockingbirdnest · Nov 20, 2021 · Nov 20, 2021 · Nov 20, 2021 · Nov 20, 2021
diff --git a/benchmarks/benchmarks.vcxproj b/benchmarks/benchmarks.vcxproj
@@ -25,6 +25,7 @@
     <ClCompile Include="..\numerics\fast_sin_cos_2π.cpp" />
     <ClCompile Include="..\physics\protector.cpp" />
     <ClCompile Include="apsides.cpp" />
+    <ClCompile Include="discrete_trajectory.cpp" />
     <ClCompile Include="dynamic_frame.cpp" />
     <ClCompile Include="elliptic_integrals_benchmark.cpp" />
     <ClCompile Include="elliptic_functions_benchmark.cpp" />

diff --git a/benchmarks/benchmarks.vcxproj.filters b/benchmarks/benchmarks.vcxproj.filters
@@ -92,6 +92,9 @@
     <ClCompile Include="..\geometry\instant_output.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="discrete_trajectory.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="quantities.hpp">

diff --git a/benchmarks/discrete_trajectory.cpp b/benchmarks/discrete_trajectory.cpp
@@ -1,236 +1,283 @@
-
+
 // .\Release\x64\benchmarks.exe --benchmark_filter=DiscreteTrajectory
 
 #include "physics/discrete_trajectory.hpp"
 
+#include <vector>
+
+#include "astronomy/epoch.hpp"
 #include "base/not_null.hpp"
 #include "benchmark/benchmark.h"
+#include "geometry/barycentre_calculator.hpp"
 #include "geometry/frame.hpp"
 #include "ksp_plugin/frames.hpp"
+#include "physics/discrete_trajectory_types.hpp"
 #include "quantities/elementary_functions.hpp"
 #include "quantities/named_quantities.hpp"
 #include "quantities/quantities.hpp"
+#include "testing_utilities/discrete_trajectory_factories.hpp"
 
 namespace principia {
 namespace physics {
 
+using astronomy::InfiniteFuture;
 using base::make_not_null_unique;
 using base::not_null;
+using geometry::Barycentre;
 using geometry::Displacement;
 using geometry::Frame;
 using geometry::Handedness;
 using geometry::Inertial;
 using geometry::Instant;
 using geometry::Velocity;
 using ksp_plugin::World;
+using physics::internal_discrete_trajectory_types::Timeline;
 using quantities::AngularFrequency;
 using quantities::Cos;
 using quantities::Sin;
 using quantities::Time;
 using quantities::si::Metre;
 using quantities::si::Radian;
 using quantities::si::Second;
+using testing_utilities::NewCircularTrajectoryTimeline;
+using testing_utilities::NewMotionlessTrajectoryTimeline;
 
 namespace {
 
-// Creates a motionless trajectory with the given number of steps.
-not_null<std::unique_ptr<DiscreteTrajectory<World>>> CreateMotionlessTrajectory(
-    int const steps) {
-  auto trajectory = make_not_null_unique<DiscreteTrajectory<World>>();
-  for (Time t = 0 * Second; t < steps * Second; t += 1 * Second) {
-    trajectory->Append(Instant() + t, {World::origin, World::unmoving});
-  }
-  return trajectory;
-}
-
-// Creates a circular trajectory with the given number of steps in the XZ plane.
-not_null<std::unique_ptr<DiscreteTrajectory<World>>> CreateCircularTrajectory(
-    int const steps) {
-  auto trajectory = make_not_null_unique<DiscreteTrajectory<World>>();
-  constexpr AngularFrequency ω = 1 * Radian / Second;
-  for (Time t = 0 * Second; t < steps * Second; t += 1 * Second) {
-    double const sin_ωt = Sin(ω * t);
-    double const cos_ωt = Cos(ω * t);
-    trajectory->Append(
-        Instant() + t,
-        {World::origin +
-             Displacement<World>({cos_ωt * Metre, 0 * Metre, sin_ωt * Metre}),
-         Velocity<World>({-sin_ωt * Metre / Second,
-                          0 * Metre / Second,
-                          cos_ωt * Metre / 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,
+                                         std::vector<double> const& splits) {
+  DiscreteTrajectory<World> trajectory;
+  Instant const t_min = timeline.begin()->time;
+  Instant const t_max = timeline.rbegin()->time;
+  auto it_split = splits.begin();
+  std::optional<Instant> t_split;
+  for (auto const& [t, degrees_of_freedom] : timeline) {
+    // The computation of |t_split| is complicated enough that we want to do it
+    // at a single place.
+    if (!t_split.has_value()) {
+      if (it_split == splits.end()) {
+        t_split = InfiniteFuture;
+      } else {
+        double const split = *it_split;
+        CHECK_LE(0, split);
+        CHECK_LE(split, 1);
+        t_split = Barycentre<Instant, double>({t_min, t_max},
+                                              {1 - split, split});
+      }
+    }
+    if (t >= t_split) {
+      ++it_split;
+      t_split = std::nullopt;
+      trajectory.NewSegment();
+    }
+    trajectory.Append(t, degrees_of_freedom);
   }
   return trajectory;
 }
 
-// Forks |parent| at a position |pos| of the way through.
-// |parent| should be nonempty.
-// |pos| should be in [0, 1].
-not_null<DiscreteTrajectory<World>*> ForkAt(DiscreteTrajectory<World>& parent,
-                                            double const pos) {
-  CHECK(!parent.Empty());
-  CHECK_GE(pos, 0);
-  CHECK_LE(pos, 1);
-  Instant const desired_fork_time =
-      parent.t_min() + (parent.t_max() - parent.t_min()) * pos;
-  auto const fork_it = parent.LowerBound(desired_fork_time);
-  return parent.NewForkWithCopy(fork_it->time);
-}
-
 }  // namespace
 
 void BM_DiscreteTrajectoryFront(benchmark::State& state) {
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(4);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
+  Instant const t0;
+  auto const timeline = NewMotionlessTrajectoryTimeline(World::origin,
+                                                        /*Δt=*/1 * Second,
+                                                        /*t1=*/t0,
+                                                        /*t2=*/t0 + 4 * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
+  auto const& segment = trajectory.segments().back();
 
   for (auto _ : state) {
-    fork->front();
+    segment.front();
   }
 }
-
 void BM_DiscreteTrajectoryBack(benchmark::State& state) {
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(4);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
+  Instant const t0;
+  auto const timeline = NewMotionlessTrajectoryTimeline(World::origin,
+                                                        /*Δt=*/1 * Second,
+                                                        /*t1=*/t0,
+                                                        /*t2=*/t0 + 4 * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
+  auto const& segment = trajectory.segments().back();
 
   for (auto _ : state) {
-    fork->back();
+    segment.back();
   }
 }
 
 void BM_DiscreteTrajectoryBegin(benchmark::State& state) {
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(4);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
+  Instant const t0;
+  auto const timeline = NewMotionlessTrajectoryTimeline(World::origin,
+                                                        /*Δt=*/1 * Second,
+                                                        /*t1=*/t0,
+                                                        /*t2=*/t0 + 4 * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
+  auto const& segment = trajectory.segments().back();
 
   for (auto _ : state) {
-    fork->begin();
+    segment.begin();
   }
 }
 
 void BM_DiscreteTrajectoryEnd(benchmark::State& state) {
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(4);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
+  Instant const t0;
+  auto const timeline = NewMotionlessTrajectoryTimeline(World::origin,
+                                                        /*Δt=*/1 * Second,
+                                                        /*t1=*/t0,
+                                                        /*t2=*/t0 + 4 * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
+  auto const& segment = trajectory.segments().back();
 
   for (auto _ : state) {
-    fork->end();
+    segment.end();
   }
 }
 
 void BM_DiscreteTrajectoryTMin(benchmark::State& state) {
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(4);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
+  Instant const t0;
+  auto const timeline = NewMotionlessTrajectoryTimeline(World::origin,
+                                                        /*Δt=*/1 * Second,
+                                                        /*t1=*/t0,
+                                                        /*t2=*/t0 + 4 * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
+  auto const& segment = trajectory.segments().back();
 
   for (auto _ : state) {
-    fork->t_min();
+    segment.t_min();
   }
 }
 
 void BM_DiscreteTrajectoryTMax(benchmark::State& state) {
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(4);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
+  Instant const t0;
+  auto const timeline = NewMotionlessTrajectoryTimeline(World::origin,
+                                                        /*Δt=*/1 * Second,
+                                                        /*t1=*/t0,
+                                                        /*t2=*/t0 + 4 * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
+  auto const& segment = trajectory.segments().back();
 
   for (auto _ : state) {
-    fork->t_max();
+    segment.t_max();
   }
 }
 
 void BM_DiscreteTrajectoryCreateDestroy(benchmark::State& state) {
+  Instant const t0;
   int const steps = state.range(0);
+  auto const timeline =
+      NewMotionlessTrajectoryTimeline(World::origin,
+                                      /*Δt=*/1 * Second,
+                                      /*t1=*/t0,
+                                      /*t2=*/t0 + steps * Second);
   for (auto _ : state) {
-    CreateMotionlessTrajectory(steps);
+    MakeTrajectory(timeline, {0.5, 0.75});
   }
 }
 
 void BM_DiscreteTrajectoryIterate(benchmark::State& state) {
+  Instant const t0;
   int const steps = state.range(0);
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(steps);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
-
-  auto const begin = fork->begin();
-  auto const end = fork->end();
+  auto const timeline =
+      NewMotionlessTrajectoryTimeline(World::origin,
+                                      /*Δt=*/1 * Second,
+                                      /*t1=*/t0,
+                                      /*t2=*/t0 + steps * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
+
+  auto const begin = trajectory.begin();
+  auto const end = trajectory.end();
   for (auto _ : state) {
     for (auto it = begin; it != end; ++it) {
     }
   }
 }
 
 void BM_DiscreteTrajectoryReverseIterate(benchmark::State& state) {
+  Instant const t0;
   int const steps = state.range(0);
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(steps);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
-
-  auto const begin = fork->begin();
-  auto const end = fork->end();
+  auto const timeline =
+      NewMotionlessTrajectoryTimeline(World::origin,
+                                      /*Δt=*/1 * Second,
+                                      /*t1=*/t0,
+                                      /*t2=*/t0 + steps * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
+
+  auto const begin = trajectory.begin();
+  auto const end = trajectory.end();
   for (auto _ : state) {
     for (auto it = end; it != begin; --it) {
     }
   }
 }
 
 void BM_DiscreteTrajectoryFind(benchmark::State& state) {
+  Instant const t0;
   int const steps = state.range(0);
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(steps);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
+  auto const timeline =
+      NewMotionlessTrajectoryTimeline(World::origin,
+                                      /*Δt=*/1 * Second,
+                                      /*t1=*/t0,
+                                      /*t2=*/t0 + steps * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
 
   for (auto _ : state) {
-    // These times are in different segments of the fork.
-    fork->Find(Instant() + 1.0 / 3.0 * steps * Second);
-    fork->Find(Instant() + 2.0 / 3.0 * steps * Second);
-    fork->Find(Instant() + 5.0 / 6.0 * steps * Second);
+    // These times are in different segments of the trajectory.
+    trajectory.find(Instant() + 1.0 / 3.0 * steps * Second);
+    trajectory.find(Instant() + 2.0 / 3.0 * steps * Second);
+    trajectory.find(Instant() + 5.0 / 6.0 * steps * Second);
   }
 }
 
 void BM_DiscreteTrajectoryLowerBound(benchmark::State& state) {
+  Instant const t0;
   int const steps = state.range(0);
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateMotionlessTrajectory(steps);
-  not_null<DiscreteTrajectory<World>*> const fork =
-      ForkAt(*ForkAt(*trajectory, 0.5), 0.75);
+  auto const timeline =
+      NewMotionlessTrajectoryTimeline(World::origin,
+                                      /*Δt=*/1 * Second,
+                                      /*t1=*/t0,
+                                      /*t2=*/t0 + steps * Second);
+  auto const trajectory = MakeTrajectory(timeline, {0.5, 0.75});
 
   for (auto _ : state) {
-    // These times are in different segments of the fork.
-    fork->LowerBound(Instant() + 1.0 / 3.0 * steps * Second);
-    fork->LowerBound(Instant() + 2.0 / 3.0 * steps * Second);
-    fork->LowerBound(Instant() + 5.0 / 6.0 * steps * Second);
+    // These times are in different segments of the trajectory.
+    trajectory.lower_bound(Instant() + 1.0 / 3.0 * steps * Second);
+    trajectory.lower_bound(Instant() + 2.0 / 3.0 * steps * Second);
+    trajectory.lower_bound(Instant() + 5.0 / 6.0 * steps * Second);
   }
 }
 
 void BM_DiscreteTrajectoryEvaluateDegreesOfFreedomExact(
     benchmark::State& state) {
-  int const steps = state.range(0);
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateCircularTrajectory(4);
+  Instant const t0;
+  auto const timeline =
+      NewCircularTrajectoryTimeline<World>(/*ω=*/1 * Radian / Second,
+                                           /*r=*/1 * Metre,
+                                           /*Δt=*/1 * Second,
+                                           /*t1=*/t0,
+                                           /*t2=*/t0 + 4 * Second);
+  auto const trajectory = MakeTrajectory(timeline, {});
 
   Instant const t = Instant() + 1 * Second;
   for (auto _ : state) {
-    trajectory->EvaluateDegreesOfFreedom(t);
+    trajectory.EvaluateDegreesOfFreedom(t);
   }
 }
 
 void BM_DiscreteTrajectoryEvaluateDegreesOfFreedomInterpolated(
     benchmark::State& state) {
-  not_null<std::unique_ptr<DiscreteTrajectory<World>>> const trajectory =
-      CreateCircularTrajectory(4);
+  Instant const t0;
+  auto const timeline =
+      NewCircularTrajectoryTimeline<World>(/*ω=*/1 * Radian / Second,
+                                           /*r=*/1 * Metre,
+                                           /*Δt=*/1 * Second,
+                                           /*t1=*/t0,
+                                           /*t2=*/t0 + 4 * Second);
+  auto const trajectory = MakeTrajectory(timeline, {});
 
   Instant const t = Instant() + 1.5 * Second;
   for (auto _ : state) {
-    trajectory->EvaluateDegreesOfFreedom(t);
+    trajectory.EvaluateDegreesOfFreedom(t);
   }
 }
 

diff --git a/testing_utilities/discrete_trajectory_factories.hpp b/testing_utilities/discrete_trajectory_factories.hpp
@@ -18,6 +18,7 @@ namespace internal_discrete_trajectory_factories {
 
 using base::not_null;
 using geometry::Instant;
+using geometry::Position;
 using geometry::Velocity;
 using physics::DegreesOfFreedom;
 using physics::DiscreteTrajectory;
@@ -37,6 +38,14 @@ class DiscreteTrajectoryFactoriesFriend {
                              DiscreteTrajectorySegment<Frame>& segment);
 };
 
+// A motionless trajectory at the given position.
+template<typename Frame>
+Timeline<Frame>
+NewMotionlessTrajectoryTimeline(Position<Frame> const& position,
+                                Time const& Δt,
+                                Instant const& t1,
+                                Instant const& t2);
+
 // A linear trajectory with constant velocity, going through
 // |degrees_of_freedom.position()| at t = t0.  The first point is at time |t1|,
 // the last point at a time < |t2|.
@@ -99,6 +108,7 @@ void AppendTrajectoryTimeline(
 using internal_discrete_trajectory_factories::AppendTrajectoryTimeline;
 using internal_discrete_trajectory_factories::NewCircularTrajectoryTimeline;
 using internal_discrete_trajectory_factories::NewLinearTrajectoryTimeline;
+using internal_discrete_trajectory_factories::NewMotionlessTrajectoryTimeline;
 
 }  // namespace testing_utilities
 }  // namespace principia

diff --git a/testing_utilities/discrete_trajectory_factories_body.hpp b/testing_utilities/discrete_trajectory_factories_body.hpp
@@ -32,6 +32,15 @@ absl::Status DiscreteTrajectoryFactoriesFriend<Frame>::Append(
   return segment.Append(t, degrees_of_freedom);
 }
 
+template<typename Frame>
+Timeline<Frame> NewMotionlessTrajectoryTimeline(Position<Frame> const& position,
+                                                Time const& Δt,
+                                                Instant const& t1,
+                                                Instant const& t2) {
+  return NewLinearTrajectoryTimeline(
+      DegreesOfFreedom<Frame>(position, Velocity<Frame>()), Δt, t1, t2);
+}
+
 template<typename Frame>
 Timeline<Frame>
 NewLinearTrajectoryTimeline(DegreesOfFreedom<Frame> const& degrees_of_freedom,

diff --git a/testing_utilities/discrete_trajectory_factories_test.cpp b/testing_utilities/discrete_trajectory_factories_test.cpp
@@ -39,6 +39,30 @@ class DiscreteTrajectoryFactoriesTest : public ::testing::Test {
                       serialization::Frame::TEST>;
 };
 
+TEST_F(DiscreteTrajectoryFactoriesTest, NewMotionlessTrajectoryTimeline) {
+  auto const timeline = NewMotionlessTrajectoryTimeline<World>(
+      /*position=*/
+      World::origin + Displacement<World>({30 * Metre, 40 * Metre, 50 * Metre}),
+      /*Δt=*/0.1 * Second,
+      /*t1=*/Instant() + 4 * Second,
+      /*t2=*/Instant() + 42 * Second);
+
+  for (auto const& [time, degrees_of_freedom] : timeline) {
+    Position<World> const& position = degrees_of_freedom.position();
+    Velocity<World> const& velocity = degrees_of_freedom.velocity();
+
+    EXPECT_THAT((position - World::origin).Norm(),
+                AlmostEquals(Sqrt(5000) * Metre, 0));
+    EXPECT_THAT(velocity.Norm(),
+                AlmostEquals(0 * Metre / Second, 0));
+  }
+  EXPECT_THAT(timeline.begin()->time,
+              AlmostEquals(Instant() + 4 * Second, 0));
+  EXPECT_THAT(timeline.rbegin()->time,
+              AlmostEquals(Instant() + 41.9 * Second, 46));
+  EXPECT_EQ(380, timeline.size());
+}
+
 TEST_F(DiscreteTrajectoryFactoriesTest, NewLinearTrajectoryTimeline) {
   auto const timeline = NewLinearTrajectoryTimeline<World>(
       /*degrees_of_freedom=*/