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

Commits on Nov 20, 2021

  1. Copy the full SHA
    e884524 View commit details
  2. Add the discrete trajectory benchmark to our project and convert it t…

    …o the new DiscreteTrajectory.
    pleroy committed Nov 20, 2021
    Copy the full SHA
    2b79fdf View commit details
  3. Lint.

    pleroy committed Nov 20, 2021
    Copy the full SHA
    7be6eff View commit details
  4. Merge pull request #3206 from pleroy/Benchmark

    Benchmark for DiscreteTrajectory
    pleroy authored Nov 20, 2021
    Copy the full SHA
    3e7ad79 View commit details
1 change: 1 addition & 0 deletions benchmarks/benchmarks.vcxproj
Original file line number Diff line number Diff line change
@@ -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" />
3 changes: 3 additions & 0 deletions benchmarks/benchmarks.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -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">
265 changes: 156 additions & 109 deletions benchmarks/discrete_trajectory.cpp
Original file line number Diff line number Diff line change
@@ -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);
}
}

10 changes: 10 additions & 0 deletions testing_utilities/discrete_trajectory_factories.hpp
Original file line number Diff line number Diff line change
@@ -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
9 changes: 9 additions & 0 deletions testing_utilities/discrete_trajectory_factories_body.hpp
Original file line number Diff line number Diff line change
@@ -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,
24 changes: 24 additions & 0 deletions testing_utilities/discrete_trajectory_factories_test.cpp
Original file line number Diff line number Diff line change
@@ -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=*/