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

Commits on May 19, 2021

  1. New polynomials in EllipticK.

    pleroy committed May 19, 2021
    Copy the full SHA
    197d2ae View commit details
  2. Copy the full SHA
    50bdd83 View commit details
  3. Copy the full SHA
    1bb171d View commit details

Commits on May 21, 2021

  1. Mixed affine operators.

    pleroy committed May 21, 2021
    Copy the full SHA
    bbe6df3 View commit details

Commits on May 25, 2021

  1. Copy the full SHA
    04fda73 View commit details
  2. Copy the full SHA
    141d576 View commit details
  3. Copy the full SHA
    c2b82d6 View commit details
  4. Copy the full SHA
    7515b19 View commit details
  5. Value, not Argument.

    pleroy committed May 25, 2021
    Copy the full SHA
    3710124 View commit details
  6. Copy the full SHA
    6aad468 View commit details

Commits on May 26, 2021

  1. Merge.

    pleroy committed May 26, 2021
    Copy the full SHA
    db84fbe View commit details
  2. Botched merge.

    pleroy committed May 26, 2021
    Copy the full SHA
    f597fa0 View commit details
  3. Affine-valued trajectories.

    pleroy committed May 26, 2021
    Copy the full SHA
    4ba7122 View commit details
  4. Compatibility code.

    pleroy committed May 26, 2021
    Copy the full SHA
    9cdf4ec View commit details
  5. Merge pull request #2999 from pleroy/UsePolynomial

    Use affine-valued polynomials for continuous trajectories
    pleroy authored May 26, 2021
    Copy the full SHA
    10d34b9 View commit details
Showing with 120 additions and 43 deletions.
  1. +1 −2 numerics/newhall_body.hpp
  2. +5 −6 physics/continuous_trajectory.hpp
  3. +44 −21 physics/continuous_trajectory_body.hpp
  4. +70 −14 physics/continuous_trajectory_test.cpp
3 changes: 1 addition & 2 deletions numerics/newhall_body.hpp
Original file line number Diff line number Diff line change
@@ -150,7 +150,6 @@ NewhallApproximationInЧебышёвBasis(int degree,
CHECK_EQ(divisions + 1, q.size());
CHECK_EQ(divisions + 1, v.size());

Vector const origin{};
Time const duration_over_two = 0.5 * (t_max - t_min);

// Tricky. The order in Newhall's matrices is such that the entries for the
@@ -159,7 +158,7 @@ NewhallApproximationInЧебышёвBasis(int degree,
for (int i = 0, j = 2 * divisions;
i < divisions + 1 && j >= 0;
++i, j -= 2) {
qv[j] = q[i] - origin;
qv[j] = q[i];
qv[j + 1] = v[i] * duration_over_two;
}

11 changes: 5 additions & 6 deletions physics/continuous_trajectory.hpp
Original file line number Diff line number Diff line change
@@ -138,14 +138,13 @@ class ContinuousTrajectory : public Trajectory<Frame> {
// never need to extract their |t_min|. Logically, the |t_min| for a
// polynomial is the |t_max| of the previous one. The first polynomial has a
// |t_min| which is |*first_time_|.
// TODO(phl): These should be polynomials returning Position<Frame>.
struct InstantPolynomialPair {
InstantPolynomialPair(
Instant t_max,
not_null<std::unique_ptr<Polynomial<Displacement<Frame>, Instant>>>
not_null<std::unique_ptr<Polynomial<Position<Frame>, Instant>>>
polynomial);
Instant t_max;
not_null<std::unique_ptr<Polynomial<Displacement<Frame>, Instant>>>
not_null<std::unique_ptr<Polynomial<Position<Frame>, Instant>>>
polynomial;
};
using InstantPolynomialPairs = std::vector<InstantPolynomialPair>;
@@ -160,10 +159,10 @@ class ContinuousTrajectory : public Trajectory<Frame> {
Instant t_max_locked() const REQUIRES_SHARED(lock_);

// Really a static method, but may be overridden for testing.
virtual not_null<std::unique_ptr<Polynomial<Displacement<Frame>, Instant>>>
virtual not_null<std::unique_ptr<Polynomial<Position<Frame>, Instant>>>
NewhallApproximationInMonomialBasis(
int degree,
std::vector<Displacement<Frame>> const& q,
std::vector<Position<Frame>> const& q,
std::vector<Velocity<Frame>> const& v,
Instant const& t_min,
Instant const& t_max,
@@ -175,7 +174,7 @@ class ContinuousTrajectory : public Trajectory<Frame> {
// instabilities.
Status ComputeBestNewhallApproximation(
Instant const& time,
std::vector<Displacement<Frame>> const& q,
std::vector<Position<Frame>> const& q,
std::vector<Velocity<Frame>> const& v) REQUIRES(lock_);

// Returns an iterator to the polynomial applicable for the given |time|, or
65 changes: 44 additions & 21 deletions physics/continuous_trajectory_body.hpp
Original file line number Diff line number Diff line change
@@ -105,16 +105,16 @@ Status ContinuousTrajectory<Frame>::Append(
if (last_points_.size() == divisions) {
// These vectors are thread-local to avoid deallocation/reallocation each
// time we go through this code path.
thread_local std::vector<Displacement<Frame>> q(divisions + 1);
thread_local std::vector<Position<Frame>> q(divisions + 1);
thread_local std::vector<Velocity<Frame>> v(divisions + 1);
q.clear();
v.clear();

for (auto const& [_, degrees_of_freedom] : last_points_) {
q.push_back(degrees_of_freedom.position() - Frame::origin);
q.push_back(degrees_of_freedom.position());
v.push_back(degrees_of_freedom.velocity());
}
q.push_back(degrees_of_freedom.position() - Frame::origin);
q.push_back(degrees_of_freedom.position());
v.push_back(degrees_of_freedom.velocity());

status = ComputeBestNewhallApproximation(time, q, v);
@@ -192,7 +192,7 @@ Position<Frame> ContinuousTrajectory<Frame>::EvaluatePosition(
auto const it = FindPolynomialForInstant(time);
CHECK(it != polynomials_.end());
auto const& polynomial = *it->polynomial;
return polynomial(time) + Frame::origin;
return polynomial(time);
}

template<typename Frame>
@@ -216,7 +216,7 @@ DegreesOfFreedom<Frame> ContinuousTrajectory<Frame>::EvaluateDegreesOfFreedom(
auto const it = FindPolynomialForInstant(time);
CHECK(it != polynomials_.end());
auto const& polynomial = *it->polynomial;
return DegreesOfFreedom<Frame>(polynomial(time) + Frame::origin,
return DegreesOfFreedom<Frame>(polynomial(time),
polynomial.EvaluateDerivative(time));
}

@@ -261,8 +261,8 @@ int ContinuousTrajectory<Frame>::PiecewisePoissonSeriesDegree(
template<typename Frame>
template<int aperiodic_degree, int periodic_degree>
PiecewisePoissonSeries<Displacement<Frame>,
aperiodic_degree, periodic_degree,
EstrinEvaluator>
aperiodic_degree, periodic_degree,
EstrinEvaluator>
ContinuousTrajectory<Frame>::ToPiecewisePoissonSeries(
Instant const& t_min,
Instant const& t_max) const {
@@ -383,6 +383,11 @@ ContinuousTrajectory<Frame>::ReadFromMessage(
message.has_is_unstable() &&
message.has_degree() &&
message.has_degree_age();
bool const is_pre_gröbner =
message.instant_polynomial_pair_size() > 0 &&
message.instant_polynomial_pair(0).polynomial().
GetExtension(serialization::PolynomialInMonomialBasis::extension).
coefficient(0).has_multivector();

not_null<std::unique_ptr<ContinuousTrajectory<Frame>>> continuous_trajectory =
std::make_unique<ContinuousTrajectory<Frame>>(
@@ -396,10 +401,10 @@ ContinuousTrajectory<Frame>::ReadFromMessage(
ЧебышёвSeries<Displacement<Frame>>::ReadFromMessage(s);
Time const step = (series.t_max() - series.t_min()) / divisions;
Instant t = series.t_min();
std::vector<Displacement<Frame>> q;
std::vector<Position<Frame>> q;
std::vector<Velocity<Frame>> v;
for (int i = 0; i <= divisions; t += step, ++i) {
q.push_back(series.Evaluate(t));
q.push_back(series.Evaluate(t) + Frame::origin);
v.push_back(series.EvaluateDerivative(t));
}
Displacement<Frame> error_estimate; // Should we do something with this?
@@ -413,10 +418,28 @@ ContinuousTrajectory<Frame>::ReadFromMessage(
}
} else {
for (auto const& pair : message.instant_polynomial_pair()) {
continuous_trajectory->polynomials_.emplace_back(
Instant::ReadFromMessage(pair.t_max()),
Polynomial<Displacement<Frame>, Instant>::template ReadFromMessage<
EstrinEvaluator>(pair.polynomial()));
if (is_pre_gröbner) {
// The easiest way to implement compatibility is to patch the serialized
// form.
serialization::Polynomial polynomial = pair.polynomial();
auto const coefficient0_multivector = polynomial.GetExtension(
serialization::PolynomialInMonomialBasis::extension).
coefficient(0).multivector();
auto* const coefficient0_point = polynomial.MutableExtension(
serialization::PolynomialInMonomialBasis::extension)->
mutable_coefficient(0)->mutable_point();
*coefficient0_point->mutable_multivector() = coefficient0_multivector;

continuous_trajectory->polynomials_.emplace_back(
Instant::ReadFromMessage(pair.t_max()),
Polynomial<Position<Frame>, Instant>::template ReadFromMessage<
EstrinEvaluator>(polynomial));
} else {
continuous_trajectory->polynomials_.emplace_back(
Instant::ReadFromMessage(pair.t_max()),
Polynomial<Position<Frame>, Instant>::template ReadFromMessage<
EstrinEvaluator>(pair.polynomial()));
}
}
}
if (message.has_first_time()) {
@@ -474,7 +497,7 @@ ContinuousTrajectory<Frame>::ContinuousTrajectory()
template<typename Frame>
ContinuousTrajectory<Frame>::InstantPolynomialPair::InstantPolynomialPair(
Instant const t_max,
not_null<std::unique_ptr<Polynomial<Displacement<Frame>, Instant>>>
not_null<std::unique_ptr<Polynomial<Position<Frame>, Instant>>>
polynomial)
: t_max(t_max),
polynomial(std::move(polynomial)) {}
@@ -554,25 +577,25 @@ Instant ContinuousTrajectory<Frame>::t_max_locked() const {
}

template<typename Frame>
not_null<std::unique_ptr<Polynomial<Displacement<Frame>, Instant>>>
not_null<std::unique_ptr<Polynomial<Position<Frame>, Instant>>>
ContinuousTrajectory<Frame>::NewhallApproximationInMonomialBasis(
int degree,
std::vector<Displacement<Frame>> const& q,
std::vector<Position<Frame>> const& q,
std::vector<Velocity<Frame>> const& v,
Instant const& t_min,
Instant const& t_max,
Displacement<Frame>& error_estimate) const {
return numerics::NewhallApproximationInMonomialBasis<
Displacement<Frame>, EstrinEvaluator>(degree,
q, v,
t_min, t_max,
error_estimate);
Position<Frame>, EstrinEvaluator>(degree,
q, v,
t_min, t_max,
error_estimate);
}

template<typename Frame>
Status ContinuousTrajectory<Frame>::ComputeBestNewhallApproximation(
Instant const& time,
std::vector<Displacement<Frame>> const& q,
std::vector<Position<Frame>> const& q,
std::vector<Velocity<Frame>> const& v) {
lock_.AssertHeld();
Length const previous_adjusted_tolerance = adjusted_tolerance_;
84 changes: 70 additions & 14 deletions physics/continuous_trajectory_test.cpp
Original file line number Diff line number Diff line change
@@ -66,10 +66,10 @@ class TestableContinuousTrajectory : public ContinuousTrajectory<Frame> {
using ContinuousTrajectory<Frame>::ContinuousTrajectory;

// Mock the Newhall factory.
not_null<std::unique_ptr<Polynomial<Displacement<Frame>, Instant>>>
not_null<std::unique_ptr<Polynomial<Position<Frame>, Instant>>>
NewhallApproximationInMonomialBasis(
int degree,
std::vector<Displacement<Frame>> const& q,
std::vector<Position<Frame>> const& q,
std::vector<Velocity<Frame>> const& v,
Instant const& t_min,
Instant const& t_max,
@@ -78,17 +78,17 @@ class TestableContinuousTrajectory : public ContinuousTrajectory<Frame> {
MOCK_CONST_METHOD7_T(
FillNewhallApproximationInMonomialBasis,
void(int degree,
std::vector<Displacement<Frame>> const& q,
std::vector<Position<Frame>> const& q,
std::vector<Velocity<Frame>> const& v,
Instant const& t_min,
Instant const& t_max,
Displacement<Frame>& error_estimate,
not_null<std::unique_ptr<Polynomial<Displacement<Frame>, Instant>>>&
not_null<std::unique_ptr<Polynomial<Position<Frame>, Instant>>>&
polynomial));

Status LockAndComputeBestNewhallApproximation(
Instant const& time,
std::vector<Displacement<Frame>> const& q,
std::vector<Position<Frame>> const& q,
std::vector<Velocity<Frame>> const& v);

// Helpers to access the internal state of the Newhall optimization.
@@ -99,19 +99,19 @@ class TestableContinuousTrajectory : public ContinuousTrajectory<Frame> {
};

template<typename Frame>
not_null<std::unique_ptr<Polynomial<Displacement<Frame>, Instant>>>
not_null<std::unique_ptr<Polynomial<Position<Frame>, Instant>>>
TestableContinuousTrajectory<Frame>::NewhallApproximationInMonomialBasis(
int degree,
std::vector<Displacement<Frame>> const& q,
std::vector<Position<Frame>> const& q,
std::vector<Velocity<Frame>> const& v,
Instant const& t_min,
Instant const& t_max,
Displacement<Frame>& error_estimate) const {
using P = PolynomialInMonomialBasis<
Displacement<Frame>, Instant, /*degree=*/1, HornerEvaluator>;
typename P::Coefficients const coefficients = {Displacement<Frame>(),
Position<Frame>, Instant, /*degree=*/1, HornerEvaluator>;
typename P::Coefficients const coefficients = {Position<Frame>(),
Velocity<Frame>()};
not_null<std::unique_ptr<Polynomial<Displacement<Frame>, Instant>>>
not_null<std::unique_ptr<Polynomial<Position<Frame>, Instant>>>
polynomial = make_not_null_unique<P>(coefficients, Instant());
FillNewhallApproximationInMonomialBasis(degree,
q, v,
@@ -125,7 +125,7 @@ template<typename Frame>
Status
TestableContinuousTrajectory<Frame>::LockAndComputeBestNewhallApproximation(
Instant const& time,
std::vector<Displacement<Frame>> const& q,
std::vector<Position<Frame>> const& q,
std::vector<Velocity<Frame>> const& v) {
absl::MutexLock l(&this->lock_);
return this->ComputeBestNewhallApproximation(time, q, v);
@@ -183,7 +183,7 @@ TEST_F(ContinuousTrajectoryTest, BestNewhallApproximation) {
Time const step = 1 * Second;
Length const tolerance = 1 * Metre;
Instant t = t0_;
std::vector<Displacement<World>> const q;
std::vector<Position<World>> const q;
std::vector<Velocity<World>> const v;

auto const trajectory = std::make_unique<TestableContinuousTrajectory<World>>(
@@ -817,8 +817,8 @@ TEST_F(ContinuousTrajectoryTest, PreCohenCompatibility) {
message2.instant_polynomial_pair(0).polynomial().GetExtension(
serialization::PolynomialInMonomialBasis::extension);
EXPECT_EQ(-2,
polynomial_in_monomial_basis.coefficient(0).multivector().vector().
x().quantity().magnitude());
polynomial_in_monomial_basis.coefficient(0).point().multivector().
vector().x().quantity().magnitude());
EXPECT_EQ(-14,
polynomial_in_monomial_basis.coefficient(1).multivector().vector().
x().quantity().magnitude());
@@ -887,6 +887,62 @@ TEST_F(ContinuousTrajectoryTest, PreGrassmannCompatibility) {
EXPECT_THAT(message2, EqualsProto(message1));
}

TEST_F(ContinuousTrajectoryTest, PreGröbnerCompatibility) {
int const number_of_steps = 30;
Time const step = 0.01 * Second;
Length const tolerance = 0.1 * Metre;

// Fill a ContinuousTrajectory and take a checkpoint.
auto position_function =
[this](Instant const t) {
return World::origin +
Displacement<World>({(t - t0_) * 3 * Metre / Second,
(t - t0_) * 5 * Metre / Second,
(t - t0_) * (-2) * Metre / Second});
};
auto velocity_function =
[](Instant const t) {
return Velocity<World>({3 * Metre / Second,
5 * Metre / Second,
-2 * Metre / Second});
};

auto const trajectory1 = std::make_unique<ContinuousTrajectory<World>>(
step, tolerance);
FillTrajectory(number_of_steps,
step,
position_function,
velocity_function,
t0_,
*trajectory1);
Instant const checkpoint_time = trajectory1->t_max();
trajectory1->checkpointer().WriteToCheckpoint(checkpoint_time);

serialization::ContinuousTrajectory message1;
trajectory1->WriteToMessage(&message1);

// Fill the pre-Gröbner fields and clear the post-Gröbner fields.
serialization::ContinuousTrajectory pre_gröbner = message1;
for (int i = 0; i < pre_gröbner.instant_polynomial_pair_size(); ++i) {
auto const coefficient0 = pre_gröbner.instant_polynomial_pair(i).
polynomial().
GetExtension(serialization::PolynomialInMonomialBasis::extension).
coefficient(0).point().multivector();
*pre_gröbner.mutable_instant_polynomial_pair(i)->mutable_polynomial()->
MutableExtension(serialization::PolynomialInMonomialBasis::extension)->
mutable_coefficient(0)->mutable_multivector() = coefficient0;
}

// Read from the pre-Gröbner message, write to a second message, and check
// that we get the same result.
auto const trajectory2 =
ContinuousTrajectory<World>::ReadFromMessage(pre_gröbner);
serialization::ContinuousTrajectory message2;
trajectory2->WriteToMessage(&message2);

EXPECT_THAT(message2, EqualsProto(message1));
}

TEST_F(ContinuousTrajectoryTest, Checkpoint) {
int const number_of_steps1 = 30;
int const number_of_steps2 = 20;