Skip to content

Commit

Permalink
Merge remote-tracking branch 'la-vache/master' into concept-barycentre
Browse files Browse the repository at this point in the history
  • Loading branch information
eggrobin committed Mar 28, 2024
2 parents 30be0ee + fe7f2d4 commit 43f9717
Show file tree
Hide file tree
Showing 35 changed files with 899 additions and 903 deletions.
2 changes: 1 addition & 1 deletion astronomy/orbit_analysis_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -227,7 +227,7 @@ class OrbitAnalysisTest : public ::testing::Test {
// general assembly (1976), commission 4, recommendation 5, note 2, and 21st
// general assembly (1991), resolution A4, recommendation IV, note 4. We
// can therefore work with this formula in TT.
PolynomialInMonomialBasis<Angle, Instant, 2, EstrinEvaluator> const
PolynomialInMonomialBasis<Angle, Instant, 2> const
newcomb_mean_longitude(
{279 * Degree + 41 * ArcMinute + 48.04 * ArcSecond,
129'602'768.13 * ArcSecond / (100 * JulianYear),
Expand Down
5 changes: 2 additions & 3 deletions benchmarks/newhall_benchmark.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -117,13 +117,12 @@ BENCHMARK_TEMPLATE(
BENCHMARK_TEMPLATE(
BM_NewhallApproximationDouble,
ResultMonomialDouble,
(&NewhallApproximationInMonomialBasis<double, EstrinEvaluator>))
(&NewhallApproximationInMonomialBasis<double, Estrin>))
->Arg(4)->Arg(8)->Arg(16);
BENCHMARK_TEMPLATE(
BM_NewhallApproximationDisplacement,
ResultMonomialDisplacement,
(&NewhallApproximationInMonomialBasis<Displacement<ICRS>,
EstrinEvaluator>))
(&NewhallApproximationInMonomialBasis<Displacement<ICRS>, Estrin>))
->Arg(4)->Arg(8)->Arg(16);

} // namespace numerics
Expand Down
28 changes: 10 additions & 18 deletions benchmarks/polynomial_in_monomial_basis_benchmark.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -80,11 +80,11 @@ struct RandomTupleGenerator<Tuple, size, size> {
template<typename Value, typename Argument, int degree,
template<typename, typename, int> class Evaluator>
void EvaluatePolynomialInMonomialBasis(benchmark::State& state) {
using P = PolynomialInMonomialBasis<Value, Argument, degree, Evaluator>;
using P = PolynomialInMonomialBasis<Value, Argument, degree>;
std::mt19937_64 random(42);
typename P::Coefficients coefficients;
RandomTupleGenerator<typename P::Coefficients, 0>::Fill(coefficients, random);
P const p(coefficients);
P const p(coefficients, with_evaluator<Evaluator>);

auto const min = ValueGenerator<Argument>::Get(random);
auto const max = ValueGenerator<Argument>::Get(random);
Expand Down Expand Up @@ -218,29 +218,21 @@ void BM_EvaluatePolynomialInMonomialBasisDisplacement(benchmark::State& state) {
}
}

BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisDouble,
EstrinEvaluator)
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisDouble, Estrin)
->Arg(4)->Arg(8)->Arg(12)->Arg(16)->Unit(benchmark::kMicrosecond);
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisQuantity,
EstrinEvaluator)
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisQuantity, Estrin)
->Arg(4)->Arg(8)->Arg(12)->Arg(16)->Unit(benchmark::kMicrosecond);
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisVectorDouble,
EstrinEvaluator)
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisVectorDouble, Estrin)
->Arg(4)->Arg(8)->Arg(12)->Arg(16)->Unit(benchmark::kMicrosecond);
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisDisplacement,
EstrinEvaluator)
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisDisplacement, Estrin)
->Arg(4)->Arg(8)->Arg(12)->Arg(16)->Unit(benchmark::kMicrosecond);
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisDouble,
HornerEvaluator)
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisDouble, Horner)
->Arg(4)->Arg(8)->Arg(12)->Arg(16)->Unit(benchmark::kMicrosecond);
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisQuantity,
HornerEvaluator)
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisQuantity, Horner)
->Arg(4)->Arg(8)->Arg(12)->Arg(16)->Unit(benchmark::kMicrosecond);
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisVectorDouble,
HornerEvaluator)
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisVectorDouble, Horner)
->Arg(4)->Arg(8)->Arg(12)->Arg(16)->Unit(benchmark::kMicrosecond);
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisDisplacement,
HornerEvaluator)
BENCHMARK_TEMPLATE1(BM_EvaluatePolynomialInMonomialBasisDisplacement, Horner)
->Arg(4)->Arg(8)->Arg(12)->Arg(16)->Unit(benchmark::kMicrosecond);

} // namespace numerics
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
#include "integrators/methods.hpp"
#include "integrators/ordinary_differential_equations.hpp"
#include "numerics/legendre.hpp"
#include "numerics/polynomial_evaluators.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/quantities.hpp"
Expand All @@ -39,7 +38,6 @@ using namespace principia::integrators::_integrators;
using namespace principia::integrators::_methods;
using namespace principia::integrators::_ordinary_differential_equations;
using namespace principia::numerics::_legendre;
using namespace principia::numerics::_polynomial_evaluators;
using namespace principia::quantities::_elementary_functions;
using namespace principia::quantities::_named_quantities;
using namespace principia::quantities::_quantities;
Expand Down Expand Up @@ -135,12 +133,10 @@ TEST_F(EmbeddedExplicitGeneralizedRungeKuttaNyströmIntegratorTest, Legendre) {
Variation<double> max_derivative_error{};
for (ODE::State const& state : solution) {
double const x = (state.time.value - t_initial) / (1 * Second);
double const error =
AbsoluteError(LegendrePolynomial<degree, EstrinEvaluator>()(x),
state.positions[0].value);
double const error = AbsoluteError(LegendrePolynomial<degree>()(x),
state.positions[0].value);
Variation<double> const derivative_error = AbsoluteError(
LegendrePolynomial<degree, EstrinEvaluator>().Derivative()(x) /
(1 * Second),
LegendrePolynomial<degree>().Derivative()(x) / (1 * Second),
state.velocities[0].value);
max_error = std::max(max_error, error);
max_derivative_error = std::max(max_derivative_error, derivative_error);
Expand Down
6 changes: 2 additions & 4 deletions mathematica/mathematica.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -257,17 +257,15 @@ std::string ToMathematica(DiscreteTrajectoryValueType<F> const& v,
OptionalExpressIn express_in = std::nullopt);

template<typename V, typename A, int d,
template<typename, typename, int> class E,
typename OptionalExpressIn = std::nullopt_t>
std::string ToMathematicaBody(
PolynomialInMonomialBasis<V, A, d, E> const& polynomial,
PolynomialInMonomialBasis<V, A, d> const& polynomial,
OptionalExpressIn express_in);

template<typename V, typename A, int d,
template<typename, typename, int> class E,
typename OptionalExpressIn = std::nullopt_t>
std::string ToMathematica(
PolynomialInMonomialBasis<V, A, d, E> const& polynomial,
PolynomialInMonomialBasis<V, A, d> const& polynomial,
OptionalExpressIn express_in = std::nullopt);

template<typename V, typename A, int d,
Expand Down
8 changes: 3 additions & 5 deletions mathematica/mathematica_body.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -83,13 +83,12 @@ struct TupleHelper<0, Tuple, OptionalExpressIn> : not_constructible {
};

template<typename V, typename A, int d,
template<typename, typename, int> class E,
typename OptionalExpressIn>
std::string ToMathematicaBody(
PolynomialInMonomialBasis<V, A, d, E> const& polynomial,
PolynomialInMonomialBasis<V, A, d> const& polynomial,
OptionalExpressIn express_in) {
using Coefficients =
typename PolynomialInMonomialBasis<V, A, d, E>::Coefficients;
typename PolynomialInMonomialBasis<V, A, d>::Coefficients;
std::vector<std::string> coefficients;
coefficients.reserve(std::tuple_size_v<Coefficients>);
TupleHelper<std::tuple_size_v<Coefficients>,
Expand Down Expand Up @@ -537,10 +536,9 @@ std::string ToMathematica(DiscreteTrajectoryValueType<F> const& v,
}

template<typename V, typename A, int d,
template<typename, typename, int> class E,
typename OptionalExpressIn>
std::string ToMathematica(
PolynomialInMonomialBasis<V, A, d, E> const& polynomial,
PolynomialInMonomialBasis<V, A, d> const& polynomial,
OptionalExpressIn express_in) {
return RawApply("Function", {ToMathematicaBody(polynomial, express_in)});
}
Expand Down
9 changes: 4 additions & 5 deletions mathematica/mathematica_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ TEST_F(MathematicaTest, ToMathematica) {
ToMathematica(elements, PreserveUnits));
}
{
PolynomialInMonomialBasis<Length, Time, 2, HornerEvaluator> polynomial1(
PolynomialInMonomialBasis<Length, Time, 2> polynomial1(
{2 * Metre, -3 * Metre / Second, 4 * Metre / Second / Second});
EXPECT_EQ(
absl::StrReplaceAll(
Expand All @@ -240,7 +240,7 @@ TEST_F(MathematicaTest, ToMathematica) {
{"β", ToMathematica(-3 * Metre / Second, PreserveUnits)},
{"γ", ToMathematica(4 * Metre / Second / Second, PreserveUnits)}}),
ToMathematica(polynomial1, PreserveUnits));
PolynomialInMonomialBasis<Length, Instant, 2, HornerEvaluator> polynomial2(
PolynomialInMonomialBasis<Length, Instant, 2> polynomial2(
{5 * Metre, 6 * Metre / Second, -7 * Metre / Second / Second},
Instant() + 2 * Second);
EXPECT_EQ(
Expand Down Expand Up @@ -275,7 +275,7 @@ TEST_F(MathematicaTest, ToMathematica) {
ToMathematica(series, PreserveUnits));
}
{
using Series = PoissonSeries<double, 0, 0, HornerEvaluator>;
using Series = PoissonSeries<double, 0, 0, Horner>;
Instant const t0 = Instant() + 3 * Second;
Series::AperiodicPolynomial secular({1.5}, t0);
Series::PeriodicPolynomial sin({0.5}, t0);
Expand All @@ -297,8 +297,7 @@ TEST_F(MathematicaTest, ToMathematica) {
ToMathematica(series, PreserveUnits));
}
{
using PiecewiseSeries =
PiecewisePoissonSeries<double, 0, 0, HornerEvaluator>;
using PiecewiseSeries = PiecewisePoissonSeries<double, 0, 0, Horner>;
using Series = PiecewiseSeries::Series;
Instant const t0 = Instant() + 3 * Second;
Series series(Series::AperiodicPolynomial({1.5}, t0),
Expand Down
20 changes: 10 additions & 10 deletions numerics/apodization_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,79 +34,79 @@ class ApodizationTest : public ::testing::Test {
};

TEST_F(ApodizationTest, Dirichlet) {
auto a = _apodization::Dirichlet<HornerEvaluator>(t1_, t2_);
auto a = _apodization::Dirichlet<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(1, 0));
EXPECT_THAT(a(t0_), AlmostEquals(1, 0));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(1, 0));
}

TEST_F(ApodizationTest, Sine) {
auto a = _apodization::Sine<HornerEvaluator>(t1_, t2_);
auto a = _apodization::Sine<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a(t0_), AlmostEquals(Sqrt(3) / 2, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, Hann) {
auto a = _apodization::Hann<HornerEvaluator>(t1_, t2_);
auto a = _apodization::Hann<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(0, 0));
EXPECT_THAT(a(t0_), AlmostEquals(0.75, 0));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, Hamming) {
auto a = _apodization::Hamming<HornerEvaluator>(t1_, t2_);
auto a = _apodization::Hamming<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(2.0 / 23.0, 0));
EXPECT_THAT(a(t0_), AlmostEquals(71.0 / 92.0, 0));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(2.0 / 23.0, 0));
}

TEST_F(ApodizationTest, Blackman) {
auto a = _apodization::Blackman<HornerEvaluator>(t1_, t2_);
auto a = _apodization::Blackman<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a(t0_), AlmostEquals(0.63, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 1));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, ExactBlackman) {
auto a = _apodization::ExactBlackman<HornerEvaluator>(t1_, t2_);
auto a = _apodization::ExactBlackman<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(8.0 / 1163.0, 37));
EXPECT_THAT(a(t0_), AlmostEquals(11843.0 / 18608.0, 2));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(8.0 / 1163.0, 37));
}

TEST_F(ApodizationTest, Nuttall) {
auto a = _apodization::Nuttall<HornerEvaluator>(t1_, t2_);
auto a = _apodization::Nuttall<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a(t0_), AlmostEquals(0.514746, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, BlackmanNuttall) {
auto a = _apodization::BlackmanNuttall<HornerEvaluator>(t1_, t2_);
auto a = _apodization::BlackmanNuttall<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(0.0003628, 703));
EXPECT_THAT(a(t0_), AlmostEquals(0.5292298, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(0.0003628, 703));
}

TEST_F(ApodizationTest, BlackmanHarris) {
auto a = _apodization::BlackmanHarris<HornerEvaluator>(t1_, t2_);
auto a = _apodization::BlackmanHarris<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(0.00006, 151));
EXPECT_THAT(a(t0_), AlmostEquals(0.520575, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(0.00006, 151));
}

TEST_F(ApodizationTest, ISO18431_2) {
auto a = _apodization::ISO18431_2<HornerEvaluator>(t1_, t2_);
auto a = _apodization::ISO18431_2<Horner>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(-0.00195313 / 4.63867187, 440));
EXPECT_THAT(a(t0_), AlmostEquals(0.9194336 / 4.63867187, 5));
EXPECT_THAT(a(mid_), AlmostEquals(1, 1));
Expand Down
8 changes: 4 additions & 4 deletions numerics/elliptic_functions.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,11 +34,11 @@ void JacobiSNCNDNReduced(Angle const& u,
// Maclaurin series for Fukushima b₀. These are polynomials in m that are used
// as coefficients of a polynomial in u₀². The index gives the corresponding
// power of u₀².
PolynomialInMonomialBasis<double, double, 0, HornerEvaluator>
PolynomialInMonomialBasis<double, double, 0>
fukushima_b₀_maclaurin_m_1(std::make_tuple(1.0 / 2.0));
PolynomialInMonomialBasis<double, double, 1, HornerEvaluator>
PolynomialInMonomialBasis<double, double, 1>
fukushima_b₀_maclaurin_m_2(std::make_tuple(-1.0 / 24.0, -1.0 / 6.0));
PolynomialInMonomialBasis<double, double, 2, HornerEvaluator>
PolynomialInMonomialBasis<double, double, 2>
fukushima_b₀_maclaurin_m_3(std::make_tuple(1.0 / 720.0,
11.0 / 180.0,
1.0 / 45.0));
Expand Down Expand Up @@ -80,7 +80,7 @@ void JacobiSNCNDNReduced(Angle const& u,
double const b₀1 = fukushima_b₀_maclaurin_m_1(m);
double const b₀2 = fukushima_b₀_maclaurin_m_2(m);
double const b₀3 = fukushima_b₀_maclaurin_m_3(m);
PolynomialInMonomialBasis<double, double, 3, HornerEvaluator>
PolynomialInMonomialBasis<double, double, 3>
fukushima_b₀_maclaurin_u₀²_3(std::make_tuple(0.0, b₀1, b₀2, b₀3));
double const u₀² = (u₀ * u₀) / Pow<2>(Radian);

Expand Down
Loading

0 comments on commit 43f9717

Please sign in to comment.