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

Commits on Nov 29, 2020

  1. Accuracy analysis.

    pleroy committed Nov 29, 2020
    Copy the full SHA
    a6e1963 View commit details
  2. Unverified

    This commit is not signed, but one or more authors requires that any commit attributed to them is signed.
    Copy the full SHA
    74abed9 View commit details
  3. Copy the full SHA
    9f07519 View commit details

Commits on Nov 30, 2020

  1. Copy the full SHA
    187fc9a View commit details

Commits on Dec 1, 2020

  1. Copy the full SHA
    72ec6af View commit details
  2. max_ω.

    pleroy committed Dec 1, 2020
    Copy the full SHA
    f420005 View commit details
  3. Factor out heuristics.

    pleroy committed Dec 1, 2020
    Copy the full SHA
    20f9f3a View commit details
  4. Copy the full SHA
    d22cb39 View commit details

Commits on Dec 2, 2020

  1. Cleanup.

    pleroy committed Dec 2, 2020
    Copy the full SHA
    39c5b45 View commit details
  2. After egg's review.

    pleroy committed Dec 2, 2020
    Copy the full SHA
    da45afc View commit details
  3. Merge pull request #2804 from pleroy/PwInnerProduct3

    Use Clenshaw-Curtis in the piecewise Poisson series inner product
    pleroy authored Dec 2, 2020
    Copy the full SHA
    a1621c9 View commit details
49 changes: 13 additions & 36 deletions numerics/piecewise_poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -85,6 +85,7 @@ class PiecewisePoissonSeries {

Instant t_min() const;
Instant t_max() const;
AngularFrequency max_ω() const;

// t must be in the interval [t_min, t_max].
Value operator()(Instant const& t) const;
@@ -373,14 +374,7 @@ template<
int aperiodic_ldegree, int periodic_ldegree,
int aperiodic_rdegree, int periodic_rdegree,
int aperiodic_wdegree, int periodic_wdegree,
template<typename, typename, int> class Evaluator,
int points = (std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree) +
std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree) +
std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree)) /
2>
template<typename, typename, int> class Evaluator>
typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
PoissonSeries<LValue,
aperiodic_ldegree, periodic_ldegree,
@@ -390,21 +384,15 @@ typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
Evaluator> const& right,
PoissonSeries<double,
aperiodic_wdegree, periodic_wdegree,
Evaluator> const& weight);
Evaluator> const& weight,
std::optional<int> max_points = std::nullopt);

template<
typename LValue, typename RValue,
int aperiodic_ldegree, int periodic_ldegree,
int aperiodic_rdegree, int periodic_rdegree,
int aperiodic_wdegree, int periodic_wdegree,
template<typename, typename, int> class Evaluator,
int points = (std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree) +
std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree) +
std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree)) /
2>
template<typename, typename, int> class Evaluator>
typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
PoissonSeries<LValue,
aperiodic_ldegree, periodic_ldegree,
@@ -416,21 +404,15 @@ typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
aperiodic_wdegree, periodic_wdegree,
Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max);
Instant const& t_max,
std::optional<int> max_points = std::nullopt);

template<
typename LValue, typename RValue,
int aperiodic_ldegree, int periodic_ldegree,
int aperiodic_rdegree, int periodic_rdegree,
int aperiodic_wdegree, int periodic_wdegree,
template<typename, typename, int> class Evaluator,
int points = (std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree) +
std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree) +
std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree)) /
2>
template<typename, typename, int> class Evaluator>
typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
PiecewisePoissonSeries<LValue,
aperiodic_ldegree, periodic_ldegree,
@@ -440,21 +422,15 @@ typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
Evaluator> const& right,
PoissonSeries<double,
aperiodic_wdegree, periodic_wdegree,
Evaluator> const& weight);
Evaluator> const& weight,
std::optional<int> max_points = std::nullopt);

template<
typename LValue, typename RValue,
int aperiodic_ldegree, int periodic_ldegree,
int aperiodic_rdegree, int periodic_rdegree,
int aperiodic_wdegree, int periodic_wdegree,
template<typename, typename, int> class Evaluator,
int points = (std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree) +
std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree) +
std::max(aperiodic_ldegree,
periodic_ldegree + estimated_trigonometric_degree)) /
2>
template<typename, typename, int> class Evaluator>
typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
PiecewisePoissonSeries<LValue,
aperiodic_ldegree, periodic_ldegree,
@@ -466,7 +442,8 @@ typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
aperiodic_wdegree, periodic_wdegree,
Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max);
Instant const& t_max,
std::optional<int> max_points = std::nullopt);

} // namespace internal_piecewise_poisson_series

108 changes: 62 additions & 46 deletions numerics/piecewise_poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -15,6 +15,18 @@ namespace internal_piecewise_poisson_series {
using quantities::Cos;
using quantities::Sin;

// The minimum value of the max_point parameter passed to Clenshaw-Curtis
// integration, irrespective of the frequencies of the argument function.
constexpr int clenshaw_curtis_min_points_overall = 513;

// The maximum number of points use in Clenshaw-Curtis integration for each
// period of the highest frequency of the argument function.
constexpr int clenshaw_curtis_points_per_period = 4;

// The desired relative error on Clenshaw-Curtis integration, as determined by
// two successive computations with increasing number of points.
constexpr double clenshaw_curtis_relative_error = 0x1p-32;

template<typename Value,
int aperiodic_degree_, int periodic_degree_,
template<typename, typename, int> class Evaluator>
@@ -57,6 +69,20 @@ t_max() const {
return bounds_.back();
}

template<typename Value,
int aperiodic_degree_, int periodic_degree_,
template<typename, typename, int> class Evaluator>
AngularFrequency
PiecewisePoissonSeries<Value, aperiodic_degree_, periodic_degree_, Evaluator>::
max_ω() const {
AngularFrequency max_ω =
addend_.has_value() ? addend_->max_ω() : AngularFrequency{};
for (int i = 0; i < series_.size(); ++i) {
max_ω = std::max(max_ω, series_[i].max_ω());
}
return max_ω;
}

template<typename Value,
int aperiodic_degree_, int periodic_degree_,
template<typename, typename, int> class Evaluator>
@@ -532,8 +558,7 @@ template<typename LValue, typename RValue,
int aperiodic_ldegree, int periodic_ldegree,
int aperiodic_rdegree, int periodic_rdegree,
int aperiodic_wdegree, int periodic_wdegree,
template<typename, typename, int> class Evaluator,
int points>
template<typename, typename, int> class Evaluator>
typename Hilbert<LValue, RValue>::InnerProductType
InnerProduct(PoissonSeries<LValue,
aperiodic_ldegree, periodic_ldegree,
@@ -543,22 +568,21 @@ InnerProduct(PoissonSeries<LValue,
Evaluator> const& right,
PoissonSeries<double,
aperiodic_wdegree, periodic_wdegree,
Evaluator> const& weight) {
Evaluator> const& weight,
std::optional<int> max_points) {
return InnerProduct<LValue, RValue,
aperiodic_ldegree, periodic_ldegree,
aperiodic_rdegree, periodic_rdegree,
aperiodic_wdegree, periodic_wdegree,
Evaluator,
points>(
left, right, weight, right.t_min(), right.t_max());
Evaluator>(
left, right, weight, right.t_min(), right.t_max(), max_points);
}

template<typename LValue, typename RValue,
int aperiodic_ldegree, int periodic_ldegree,
int aperiodic_rdegree, int periodic_rdegree,
int aperiodic_wdegree, int periodic_wdegree,
template<typename, typename, int> class Evaluator,
int points>
template<typename, typename, int> class Evaluator>
typename Hilbert<LValue, RValue>::InnerProductType
InnerProduct(PoissonSeries<LValue,
aperiodic_ldegree, periodic_ldegree,
@@ -570,29 +594,16 @@ InnerProduct(PoissonSeries<LValue,
aperiodic_wdegree, periodic_wdegree,
Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
using Result =
Primitive<typename Hilbert<LValue, RValue>::InnerProductType, Time>;
Result result{};
for (int i = 0; i < right.series_.size(); ++i) {
auto integrand = [i, &left, &right, &weight](Instant const& t) {
return Hilbert<LValue, RValue>::InnerProduct(
left(t) * weight(t),
right.series_[i](t) + right.EvaluateAddend(t));
};
auto const integral = quadrature::GaussLegendre<points>(
integrand, right.bounds_[i], right.bounds_[i + 1]);
result += integral;
}
return result / (t_max - t_min);
Instant const& t_max,
std::optional<int> max_points) {
return InnerProduct(right, left, weight, t_min, t_max, max_points);
}

template<typename LValue, typename RValue,
int aperiodic_ldegree, int periodic_ldegree,
int aperiodic_rdegree, int periodic_rdegree,
int aperiodic_wdegree, int periodic_wdegree,
template<typename, typename, int> class Evaluator,
int points>
template<typename, typename, int> class Evaluator>
typename Hilbert<LValue, RValue>::InnerProductType
InnerProduct(PiecewisePoissonSeries<LValue,
aperiodic_ldegree, periodic_ldegree,
@@ -602,22 +613,21 @@ InnerProduct(PiecewisePoissonSeries<LValue,
Evaluator> const& right,
PoissonSeries<double,
aperiodic_wdegree, periodic_wdegree,
Evaluator> const& weight) {
Evaluator> const& weight,
std::optional<int> max_points) {
return InnerProduct<LValue, RValue,
aperiodic_ldegree, periodic_ldegree,
aperiodic_rdegree, periodic_rdegree,
aperiodic_wdegree, periodic_wdegree,
Evaluator,
points>(
left, right, weight, left.t_min(), left.t_max());
Evaluator>(
left, right, weight, left.t_min(), left.t_max(), max_points);
}

template<typename LValue, typename RValue,
int aperiodic_ldegree, int periodic_ldegree,
int aperiodic_rdegree, int periodic_rdegree,
int aperiodic_wdegree, int periodic_wdegree,
template<typename, typename, int> class Evaluator,
int points>
template<typename, typename, int> class Evaluator>
typename Hilbert<LValue, RValue>::InnerProductType
InnerProduct(PiecewisePoissonSeries<LValue,
aperiodic_ldegree, periodic_ldegree,
@@ -629,21 +639,27 @@ InnerProduct(PiecewisePoissonSeries<LValue,
aperiodic_wdegree, periodic_wdegree,
Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
using Result =
Primitive<typename Hilbert<LValue, RValue>::InnerProductType, Time>;
Result result{};
for (int i = 0; i < left.series_.size(); ++i) {
auto integrand = [i, &left, &right, &weight](Instant const& t) {
return Hilbert<LValue, RValue>::InnerProduct(
left.series_[i](t) + left.EvaluateAddend(t),
right(t) * weight(t));
};
auto const integral = quadrature::GaussLegendre<points>(
integrand, left.bounds_[i], left.bounds_[i + 1]);
result += integral;
}
return result / (t_max - t_min);
Instant const& t_max,
std::optional<int> max_points) {
AngularFrequency const max_ω = left.max_ω() + right.max_ω() + weight.max_ω();
std::optional<int> const max_points_heuristic =
MaxPointsHeuristicsForAutomaticClenshawCurtis(
max_ω,
t_max - t_min,
clenshaw_curtis_min_points_overall,
clenshaw_curtis_points_per_period);

auto integrand = [&left, &right, &weight](Instant const& t) {
return Hilbert<LValue, RValue>::InnerProduct(left(t), right(t)) * weight(t);
};
return quadrature::AutomaticClenshawCurtis(
integrand,
t_min,
t_max,
/*max_relative_error=*/clenshaw_curtis_relative_error,
/*max_points=*/max_points.has_value() ? max_points
: max_points_heuristic) /
(t_max - t_min);
}

} // namespace internal_piecewise_poisson_series
40 changes: 20 additions & 20 deletions numerics/piecewise_poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -203,30 +203,30 @@ TEST_F(PiecewisePoissonSeriesTest, ActionMultiorigin) {
}

TEST_F(PiecewisePoissonSeriesTest, InnerProduct) {
double const d1 =
InnerProduct<double, double, 0, 0, 0, 0, 0, 0, HornerEvaluator, 8>(
pp_, p_,
apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second));
double const d2 =
InnerProduct<double, double, 0, 0, 0, 0, 0, 0, HornerEvaluator, 8>(
p_, pp_,
apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second));
EXPECT_THAT(d1, AlmostEquals((3 * π - 26) / (8 * π), 0));
EXPECT_THAT(d2, AlmostEquals((3 * π - 26) / (8 * π), 0));
double const d1 = InnerProduct(
pp_, p_,
apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second),
/*max_points=*/1 << 20);
double const d2 = InnerProduct(
p_, pp_,
apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second),
/*max_points=*/1 << 20);
EXPECT_THAT(d1, RelativeErrorFrom((3 * π - 26) / (8 * π), IsNear(3e-11_⑴)));
EXPECT_THAT(d2, RelativeErrorFrom((3 * π - 26) / (8 * π), IsNear(3e-11_⑴)));
}

TEST_F(PiecewisePoissonSeriesTest, InnerProductMultiorigin) {
auto const p = p_.AtOrigin(t0_ + 2 * Second);
double const d1 =
InnerProduct<double, double, 0, 0, 0, 0, 0, 0, HornerEvaluator, 8>(
pp_, p,
apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second));
double const d2 =
InnerProduct<double, double, 0, 0, 0, 0, 0, 0, HornerEvaluator, 8>(
p, pp_,
apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second));
EXPECT_THAT(d1, AlmostEquals((3 * π - 26) / (8 * π), 0));
EXPECT_THAT(d2, AlmostEquals((3 * π - 26) / (8 * π), 0));
double const d1 = InnerProduct(
pp_, p,
apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second),
/*max_points=*/1 << 20);
double const d2 = InnerProduct(
p, pp_,
apodization::Dirichlet<HornerEvaluator>(t0_, t0_ + 2 * Second),
/*max_points=*/1 << 20);
EXPECT_THAT(d1, RelativeErrorFrom((3 * π - 26) / (8 * π), IsNear(3e-11_⑴)));
EXPECT_THAT(d2, RelativeErrorFrom((3 * π - 26) / (8 * π), IsNear(3e-11_⑴)));
}

TEST_F(PiecewisePoissonSeriesTest, Fourier) {
10 changes: 10 additions & 0 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -100,6 +100,7 @@ class PoissonSeries {
HigherEvaluator>() const;

Instant const& origin() const;
AngularFrequency max_ω() const;

Value operator()(Instant const& t) const;

@@ -387,9 +388,18 @@ InnerProduct(PoissonSeries<LValue,
Instant const& t_min,
Instant const& t_max);

// Computes a heuristic for the maximum number of points for an oscillating
// function.
std::optional<int> MaxPointsHeuristicsForAutomaticClenshawCurtis(
AngularFrequency const& max_ω,
Time const& Δt,
int min_points_overall,
int points_per_period);

} // namespace internal_poisson_series

using internal_poisson_series::InnerProduct;
using internal_poisson_series::MaxPointsHeuristicsForAutomaticClenshawCurtis;
using internal_poisson_series::PoissonSeries;

} // namespace numerics
Loading