Merge remote-tracking branch 'la-vache/master' into container-barycentre

geometry/barycentre_calculator.hpp

@@ -46,6 +46,8 @@ template<affine Point, homogeneous_field Weight, std::size_t size>
 Point Barycentre(Point const (&points)[size], Weight const (&weights)[size]);
 template<real_affine_space Point, std::size_t size>
 Point Barycentre(Point const (&points)[size], double const (&weights)[size]);
+template<real_affine_space Point, std::size_t size>
+Point Barycentre(Point const (&points)[size]);
 
 }  // namespace internal
 

geometry/barycentre_calculator_body.hpp

@@ -71,6 +71,25 @@ Point Barycentre(Point const (&points)[size], double const (&weights)[size]) {
   return Barycentre<Point, double>(points, weights);
 }
 
+template<real_affine_space Point, std::size_t size>
+Point Barycentre(Point const (&points)[size]) {
+  static_assert(size != 0);
+  Difference<Point> total{};
+  static const Point origin{};
+  for (int i = 0; i < size; ++i) {
+    if constexpr (additive_group<Point>) {
+      total += points[i];
+    } else {
+      total += points[i] - origin;
+    }
+  }
+  if constexpr (additive_group<Point>) {
+    return total / size;
+  } else {
+    return origin + total / size;
+  }
+}
+
 }  // namespace internal
 }  // namespace _barycentre_calculator
 }  // namespace geometry

geometry/point_test.cpp

@@ -171,7 +171,7 @@ TEST_F(PointTest, Barycentres) {
   Instant const t1 = mjd0 + 1 * Day;
   Instant const t2 = mjd0 - 3 * Day;
   Instant const b1 = Barycentre({t1, t2}, {3 * Litre, 1 * Litre});
-  Instant const b2 = Barycentre({t2, t1}, {1, 1});
+  Instant const b2 = Barycentre({t2, t1});
   EXPECT_THAT(b1, AlmostEquals(mjd0, 1));
   EXPECT_THAT(b2, Eq(mjd0 - 1 * Day));
 }

ksp_plugin_test/flight_plan_optimizer_test.cpp

@@ -533,7 +533,7 @@ class MetricTest
     EXPECT_OK(ephemeris_->Prolong(desired_final_time));
     auto const earth_dof = solar_system_1950_.degrees_of_freedom("Earth");
     auto const mars_dof = solar_system_1950_.degrees_of_freedom("Mars");
-    auto const midway = Barycentre({earth_dof, mars_dof}, {1, 1});
+    auto const midway = Barycentre({earth_dof, mars_dof});
     EXPECT_OK(root_.Append(epoch_, midway));
     flight_plan_ = std::make_unique<FlightPlan>(
         /*initial_mass=*/1 * Kilogram,

mathematica/mathematica_body.hpp

@@ -125,7 +125,7 @@ std::string ToMathematicaBody(
     OptionalExpressIn express_in) {
   auto const& a = polynomial.lower_bound();
   auto const& b = polynomial.upper_bound();
-  auto const midpoint = Barycentre({a, b}, {0.5, 0.5});
+  auto const midpoint = Barycentre({a, b});
   std::string const argument = RawApply(
       "Divide",
       {RawApply("Subtract", {"#", ToMathematica(midpoint, express_in)}),

numerics/apodization_body.hpp

@@ -29,7 +29,7 @@ PoissonSeries<double, 0, 0, Evaluator> Dirichlet(Instant const& t_min,
                                                  Instant const& t_max) {
   using Result = PoissonSeries<double, 0, 0, Evaluator>;
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({1}, t_mid), {});
 }
 
@@ -40,7 +40,7 @@ PoissonSeries<double, 0, 0, Evaluator> Sine(Instant const& t_min,
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
   using PeriodicPolynomial = typename Result::PeriodicPolynomial;
   AngularFrequency const ω = π * Radian / (t_max - t_min);
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({0}, t_mid),
                 {{ω,
                   {.sin = PeriodicPolynomial({0}, t_mid),
@@ -54,7 +54,7 @@ PoissonSeries<double, 0, 0, Evaluator> Hann(Instant const& t_min,
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
   using PeriodicPolynomial = typename Result::PeriodicPolynomial;
   AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({0.5}, t_mid),
                 {{ω,
                   {.sin = PeriodicPolynomial({0}, t_mid),
@@ -68,7 +68,7 @@ PoissonSeries<double, 0, 0, Evaluator> Hamming(Instant const& t_min,
   AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
   using PeriodicPolynomial = typename Result::PeriodicPolynomial;
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({25.0 / 46.0}, t_mid),
                 {{ω,
                   {.sin = PeriodicPolynomial({0}, t_mid),
@@ -82,7 +82,7 @@ PoissonSeries<double, 0, 0, Evaluator> Blackman(Instant const& t_min,
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
   using PeriodicPolynomial = typename Result::PeriodicPolynomial;
   AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({0.42}, t_mid),
                 {{ω,
                   {.sin = PeriodicPolynomial({0}, t_mid),
@@ -99,7 +99,7 @@ PoissonSeries<double, 0, 0, Evaluator> ExactBlackman(Instant const& t_min,
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
   using PeriodicPolynomial = typename Result::PeriodicPolynomial;
   AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({3969.0 / 9304.0}, t_mid),
                 {{ω,
                   {.sin = PeriodicPolynomial({0}, t_mid),
@@ -116,7 +116,7 @@ PoissonSeries<double, 0, 0, Evaluator> Nuttall(Instant const& t_min,
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
   using PeriodicPolynomial = typename Result::PeriodicPolynomial;
   AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({0.355768}, t_mid),
                 {{ω,
                   {.sin = PeriodicPolynomial({0}, t_mid),
@@ -136,7 +136,7 @@ PoissonSeries<double, 0, 0, Evaluator> BlackmanNuttall(Instant const& t_min,
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
   using PeriodicPolynomial = typename Result::PeriodicPolynomial;
   AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({0.3635819}, t_mid),
                 {{ω,
                   {.sin = PeriodicPolynomial({0}, t_mid),
@@ -156,7 +156,7 @@ PoissonSeries<double, 0, 0, Evaluator> BlackmanHarris(Instant const& t_min,
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
   using PeriodicPolynomial = typename Result::PeriodicPolynomial;
   AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({0.35875}, t_mid),
                 {{ω,
                   {.sin = PeriodicPolynomial({0}, t_mid),
@@ -176,7 +176,7 @@ PoissonSeries<double, 0, 0, Evaluator> ISO18431_2(Instant const& t_min,
   using AperiodicPolynomial = typename Result::AperiodicPolynomial;
   using PeriodicPolynomial = typename Result::PeriodicPolynomial;
   AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return Result(AperiodicPolynomial({1.0 / 4.63867187}, t_mid),
                 {{ω,
                   {.sin = PeriodicPolynomial({0}, t_mid),

numerics/approximation_body.hpp

@@ -56,7 +56,7 @@ not_null<std::unique_ptr<
     FixedVector<Value<Argument, Function>, N / 2 + 1> const& previous_aⱼ,
     Difference<Value<Argument, Function>>* const error_estimate) {
   // This implementation follows [Boy13], section 4 and appendix A.
-  auto const midpoint = Barycentre({a, b}, {0.5, 0.5});
+  auto const midpoint = Barycentre({a, b});
 
   auto чебышёв_lobato_point =
       [&a, &b, &midpoint](std::int64_t const k) -> Argument {
@@ -149,7 +149,7 @@ bool StreamingAdaptiveЧебышёвPolynomialInterpolantImplementation(
             << full_error_estimate;
     Difference<Value<Argument, Function>> upper_error_estimate;
     Difference<Value<Argument, Function>> lower_error_estimate;
-    auto const midpoint = Barycentre({lower_bound, upper_bound}, {1, 1});
+    auto const midpoint = Barycentre({lower_bound, upper_bound});
     bool const lower_interpolants_stop =
         StreamingAdaptiveЧебышёвPolynomialInterpolantImplementation<max_degree>(
             f,

numerics/nearest_neighbour_body.hpp

@@ -155,8 +155,7 @@ PrincipalComponentPartitioningTree<Value_>::BuildTree(
   auto const mid_lower = std::max_element(begin, mid_upper, projection_less);
 
   auto const anchor = Barycentre(
-      {displacements_[mid_lower->index], displacements_[mid_upper->index]},
-      {1, 1});
+      {displacements_[mid_lower->index], displacements_[mid_upper->index]});
 
   auto first_child = BuildTree(begin, mid_upper, size / 2);
   auto second_child = BuildTree(mid_upper, end, size - size / 2);

numerics/newhall_body.hpp

@@ -321,7 +321,7 @@ NewhallApproximationInMonomialBasis(std::vector<Value> const& q,
     qv[j + 1] = v[i] * duration_over_two;
   }
 
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return (origin +
           Dehomogeneize<Difference<Value>, degree, Evaluator>(
               NewhallMonomialApproximator<
@@ -330,7 +330,7 @@ NewhallApproximationInMonomialBasis(std::vector<Value> const& q,
                   Evaluator>::HomogeneousCoefficients(qv, error_estimate),
               /*scale=*/1.0 / duration_over_two,
               t_mid))
-      .WithEvaluator<Evaluator>();
+      .template WithEvaluator<Evaluator>();
 }
 
 #define PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(degree)   \

numerics/poisson_series_basis_body.hpp

@@ -173,7 +173,7 @@ std::array<Series, sizeof...(indices)> AperiodicSeriesGenerator<
     degree, dimension,
     std::index_sequence<indices...>>::BasisElements(Instant const& t_min,
                                                     Instant const& t_max) {
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   return {(Series(
       PolynomialGenerator<typename Series::AperiodicPolynomial,
                           dimension>::template UnitPolynomial<indices>(t_min,
@@ -221,7 +221,7 @@ std::array<Series, sizeof...(indices)> PeriodicSeriesGenerator<
     std::index_sequence<indices...>>::BasisElements(AngularFrequency const& ω,
                                                     Instant const& t_min,
                                                     Instant const& t_max) {
-  Instant const t_mid = Barycentre({t_min, t_max}, {1, 1});
+  Instant const t_mid = Barycentre({t_min, t_max});
   typename Series::AperiodicPolynomial const aperiodic_zero{{}, t_mid};
   return {Series(
       aperiodic_zero,

numerics/polynomial.hpp

@@ -52,9 +52,10 @@ class Polynomial {
   virtual void WriteToMessage(
       not_null<serialization::Polynomial*> message) const = 0;
 
-  // The evaluator is not part of the serialization because it's fine to read
-  // with a different evaluator than the one the polynomial was written with.
-  // TODO(phl): Revisit when we deserialize the evaluator.
+  static not_null<std::unique_ptr<Polynomial>> ReadFromMessage(
+      serialization::Polynomial const& message);
+  // Compatibility deserialization, when the evaluator is not present in
+  // |message|.
   template<template<typename, typename, int> typename Evaluator>
   static not_null<std::unique_ptr<Polynomial>> ReadFromMessage(
       serialization::Polynomial const& message);

numerics/polynomial_body.hpp

@@ -19,11 +19,16 @@ using namespace principia::numerics::_polynomial_in_monomial_basis;
     return make_not_null_unique<                                            \
         PolynomialInMonomialBasis<Value, Argument, value>>(                 \
         PolynomialInMonomialBasis<Value, Argument, value>::ReadFromMessage( \
-            message)                                                        \
-            .WithEvaluator<Evaluator>())
+            message))
+
+#define PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(value)         \
+  case value:                                                               \
+    return make_not_null_unique<                                            \
+        PolynomialInMonomialBasis<Value, Argument, value>>(                 \
+        PolynomialInMonomialBasis<Value, Argument, value>::ReadFromMessage< \
+            Evaluator>(message))
 
 template<typename Value_, typename Argument_>
-template<template<typename, typename, int> typename Evaluator>
 not_null<std::unique_ptr<Polynomial<Value_, Argument_>>>
 Polynomial<Value_, Argument_>::ReadFromMessage(
     serialization::Polynomial const& message) {
@@ -59,7 +64,45 @@ Polynomial<Value_, Argument_>::ReadFromMessage(
   }
 }
 
+template<typename Value_, typename Argument_>
+template<template<typename, typename, int> typename Evaluator>
+not_null<std::unique_ptr<Polynomial<Value_, Argument_>>>
+Polynomial<Value_, Argument_>::ReadFromMessage(
+    serialization::Polynomial const& message) {
+  // 24 is the largest exponent that we can serialize for Quantity.
+  switch (message.degree()) {
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(1);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(2);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(3);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(4);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(5);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(6);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(7);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(8);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(9);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(10);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(11);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(12);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(13);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(14);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(15);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(16);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(17);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(18);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(19);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(20);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(21);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(22);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(23);
+    PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE(24);
+    default:
+      LOG(FATAL) << "Unexpected degree " << message.degree();
+      break;
+  }
+}
+
 #undef PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_CASE
+#undef PRINCIPIA_POLYNOMIAL_DEGREE_VALUE_COMPATIBILITY_CASE
 
 }  // namespace internal
 }  // namespace _polynomial

numerics/polynomial_evaluators.hpp

@@ -1,14 +1,17 @@
 #pragma once
 
+#include "base/not_null.hpp"
 #include "quantities/concepts.hpp"
 #include "quantities/named_quantities.hpp"
 #include "quantities/tuples.hpp"
+#include "serialization/numerics.pb.h"
 
 namespace principia {
 namespace numerics {
 namespace _polynomial_evaluators {
 namespace internal {
 
+using namespace principia::base::_not_null;
 using namespace principia::quantities::_concepts;
 using namespace principia::quantities::_named_quantities;
 using namespace principia::quantities::_tuples;
@@ -30,12 +33,18 @@ struct Evaluator {
   Value Evaluate(
       Coefficients const& coefficients,
       Argument const& argument,
-      Evaluator const* evaluator);
+      not_null<Evaluator const*> evaluator);
   FORCE_INLINE(static)
   Derivative<Value, Argument> EvaluateDerivative(
       Coefficients const& coefficients,
       Argument const& argument,
-      Evaluator const* evaluator);
+      not_null<Evaluator const*> evaluator);
+
+  static void WriteToMessage(
+      not_null<serialization::PolynomialInMonomialBasis::Evaluator*> message,
+      not_null<Evaluator const*> evaluator);
+  static not_null<Evaluator const*> ReadFromMessage(
+      serialization::PolynomialInMonomialBasis::Evaluator const& message);
 };
 
 template<typename Value, typename Argument, int degree, bool allow_fma>