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

Commits on Aug 10, 2020

  1. Compiles, doesn't link.

    pleroy committed Aug 10, 2020
    Copy the full SHA
    69284de View commit details
  2. Comments and renaming.

    pleroy committed Aug 10, 2020
    Copy the full SHA
    65679f6 View commit details
  3. Copy the full SHA
    c26cba0 View commit details
  4. Cleanup.

    pleroy committed Aug 10, 2020
    Copy the full SHA
    49b1c50 View commit details
  5. Copy the full SHA
    992ec06 View commit details
  6. More conversions.

    pleroy committed Aug 10, 2020
    Copy the full SHA
    f729a0e View commit details
  7. Copy the full SHA
    f37a2e8 View commit details

Commits on Aug 11, 2020

  1. Merge pull request #2670 from pleroy/DotProduct

    Change the dot product passed to frequency analysis to be a templated functor
    pleroy authored Aug 11, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    d102b2a View commit details
Showing with 131 additions and 70 deletions.
  1. +27 −20 numerics/frequency_analysis.hpp
  2. +16 −17 numerics/frequency_analysis_body.hpp
  3. +56 −24 numerics/frequency_analysis_test.cpp
  4. +4 −2 numerics/poisson_series.hpp
  5. +28 −7 numerics/poisson_series_body.hpp
47 changes: 27 additions & 20 deletions numerics/frequency_analysis.hpp
Original file line number Diff line number Diff line change
@@ -20,43 +20,50 @@ using quantities::Primitive;
using quantities::Product;
using quantities::Time;

// A function that implements the dot product between a |Function| and a Poisson
// series with an apodization function |weight| function that is itself a
// Poisson series. |Function| must be a functor taking an Instant.
template<typename Function,
typename RValue, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
using DotProduct =
std::function<Product<std::invoke_result_t<Function, Instant>, RValue>(
Function const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight)>;
// In this file Dot is a templated functor that implements the dot product
// between two functors and a weight. Its declaration must look like:
//
// class Dot {
// public:
// ...
// template<typename LFunction, typename RFunction, typename Weight>
// Product<std::invoke_result_t<LFunction, Instant>,
// std::invoke_result_t<RFunction, Instant>>
// operator()(LFunction const& left,
// RFunction const& right,
// Weight const& weight) const;
// ...
//};
//
// Where the implementation of Dot may assume that Weight is a Poisson series
// returning a double.

// Computes the precise mode of a quasi-periodic |function|, assuming that the
// mode is over the interval |fft_mode| (so named because it has presumably been
// obtained using FFT). See [Cha95].
template<typename Function,
int wdegree_,
typename Dot,
template<typename, typename, int> class Evaluator>
AngularFrequency PreciseMode(
Interval<AngularFrequency> const& fft_mode,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
DotProduct<Function, double, 0, wdegree_, Evaluator> const& dot);
Dot const& dot);

// Computes the Кудрявцев projection of |function| on a basis with angular
// frequency ω and maximum degree |degree_|. See [Kud07].
// TODO(phl): We really need multiple angular frequencies.
template<typename Function,
int degree_, int wdegree_,
template<int degree_,
typename Function,
int wdegree_,
typename Dot,
template<typename, typename, int> class Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>, degree_, Evaluator>
Projection(
AngularFrequency const& ω,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
DotProduct<Function, std::invoke_result_t<Function, Instant>,
degree_, wdegree_, Evaluator> const& dot);
Projection(AngularFrequency const& ω,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Dot const& dot);

} // namespace internal_frequency_analysis

33 changes: 16 additions & 17 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -124,12 +124,13 @@ BasisGenerator<Series, std::index_sequence<indices...>>::Basis(

template<typename Function,
int wdegree_,
typename Dot,
template<typename, typename, int> class Evaluator>
AngularFrequency PreciseMode(
Interval<AngularFrequency> const& fft_mode,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
DotProduct<Function, double, 0, wdegree_, Evaluator> const& dot) {
Dot const& dot) {
using Value = std::invoke_result_t<Function, Instant>;
using Degree0 = PoissonSeries<double, 0, Evaluator>;

@@ -154,16 +155,16 @@ AngularFrequency PreciseMode(
std::greater<Square<Value>>());
}

template<typename Function,
int degree_, int wdegree_,
template<int degree_,
typename Function,
int wdegree_,
typename Dot,
template<typename, typename, int> class Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>, degree_, Evaluator>
Projection(
AngularFrequency const& ω,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
DotProduct<Function, std::invoke_result_t<Function, Instant>,
degree_, wdegree_, Evaluator> const& dot) {
Projection(AngularFrequency const& ω,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Dot const& dot) {
using Value = std::invoke_result_t<Function, Instant>;
using Series = PoissonSeries<Value, degree_, Evaluator>;

@@ -174,22 +175,21 @@ Projection(
// This code follows [Kud07], section 2. Our indices start at 0, unlike those
// of Кудрявцев which start at 1.
FixedLowerTriangularMatrix<Inverse<Value>, basis_size> α;
std::vector<Function> f;

// Only indices 0 to m - 1 are used in this array. At the beginning of
// iteration m it contains Aⱼ⁽ᵐ⁻¹⁾.
std::array<double, basis_size> A;

auto const F₀ = dot(function, basis[0], weight);
// TODO(phl): This does not work if basis does not have the same type as
// Function, i.e., if the degrees don't match.
auto const Q₀₀ = dot(basis[0], basis[0], weight);
α[0][0] = 1 / Sqrt(Q₀₀);
A[0] = F₀ / Q₀₀;
f.emplace_back(function - A[0] * basis[0]);

// At the beginning of iteration m this contains fₘ₋₁.
auto f = function - A[0] * basis[0];
for (int m = 1; m < basis_size; ++m) {
// Contains Fₘ.
auto const F = dot(f[m - 1], basis[m], weight);
auto const F = dot(f, basis[m], weight);

// Only indices 0 to m are used in this array. It contains Qₘⱼ.
std::array<Square<Value>, basis_size> Q;
@@ -231,12 +231,11 @@ Projection(
}

{
PoissonSeries<double, degree_, Evaluator> Σ_αₘᵢ_eᵢ =
α[m][0] * basis[0];
PoissonSeries<double, degree_, Evaluator> Σ_αₘᵢ_eᵢ = α[m][0] * basis[0];
for (int i = 1; i <= m; ++i) {
Σ_αₘᵢ_eᵢ += α[m][i] * basis[i];
}
f.emplace_back(f[m - 1] - α[m][m] * F * Σ_αₘᵢ_eᵢ);
f -= α[m][m] * F * Σ_αₘᵢ_eᵢ;
}
}

80 changes: 56 additions & 24 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -27,6 +27,7 @@ using geometry::Instant;
using quantities::AngularFrequency;
using quantities::Length;
using quantities::Pow;
using quantities::Product;
using quantities::Sin;
using quantities::Square;
using quantities::Time;
@@ -37,8 +38,41 @@ using testing_utilities::AlmostEquals;
using testing_utilities::IsNear;
using testing_utilities::RelativeErrorFrom;
using testing_utilities::operator""_⑴;
using ::testing::AllOf;
using ::testing::Gt;
using ::testing::Lt;
namespace si = quantities::si;

class DotImplementation {
public:
DotImplementation(Instant const& t_min, Instant const& t_max);

template<typename LFunction, typename RFunction, typename Weight>
Product<std::invoke_result_t<LFunction, Instant>,
std::invoke_result_t<RFunction, Instant>>
operator()(LFunction const& left,
RFunction const& right,
Weight const& weight) const;

private:
Instant const t_min_;
Instant const t_max_;
};

DotImplementation::DotImplementation(Instant const& t_min,
Instant const& t_max)
: t_min_(t_min),
t_max_(t_max) {}

template<typename LFunction, typename RFunction, typename Weight>
Product<std::invoke_result_t<LFunction, Instant>,
std::invoke_result_t<RFunction, Instant>>
DotImplementation::operator()(LFunction const& left,
RFunction const& right,
Weight const& weight) const {
return Dot(left, right, weight, t_min_, t_max_);
}

class FrequencyAnalysisTest : public ::testing::Test {};

TEST_F(FrequencyAnalysisTest, PreciseMode) {
@@ -86,16 +120,7 @@ TEST_F(FrequencyAnalysisTest, PreciseMode) {
EXPECT_THAT(mode.midpoint(),
RelativeErrorFrom(ω, IsNear(8.1e-4_⑴)));

using Double = PoissonSeries<double, 0, HornerEvaluator>;

std::function<Length(Series const& left,
Double const& right,
Double const& weight)> const dot =
[t_min, t_max](Series const& left,
Double const& right,
Double const& weight) {
return Dot(left, right, weight, t_min, t_max);
};
DotImplementation dot(t_min, t_max);

// The precise analysis is only limited by our ability to pinpoint the
// maximum.
@@ -131,24 +156,31 @@ TEST_F(FrequencyAnalysisTest, Projection) {

Instant const t_min = t0;
Instant const t_max = t0 + 100 * Radian / ω;
DotImplementation dot(t_min, t_max);

using Double = PoissonSeries<double, 0, HornerEvaluator>;

std::function<Square<Length>(Series const& left,
Series const& right,
Double const& weight)> const dot =
[t_min, t_max](Series const& left,
Series const& right,
Double const& weight) {
return Dot(left, right, weight, t_min, t_max);
};
// Projection on a 4-th degree basis accurately reconstructs the function.
auto const projection4 = Projection<4>(
ω, series, apodization::Hann<HornerEvaluator>(t_min, t_max), dot);
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(projection4(t0 + i * Radian / ω),
AlmostEquals(series(t0 + i * Radian / ω), 0, 2688));
}

auto const projection = Projection(
ω, series, apodization::Dirichlet<HornerEvaluator>(t_min, t_max), dot);
// Projection on a 5-th degree basis is also accurate.
auto const projection5 = Projection<5>(
ω, series, apodization::Hann<HornerEvaluator>(t_min, t_max), dot);
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(projection5(t0 + i * Radian / ω),
AlmostEquals(series(t0 + i * Radian / ω), 0, 8000));
}

// Projection on a 3-rd degree basis introduces significant errors.
auto const projection3 = Projection<3>(
ω, series, apodization::Hann<HornerEvaluator>(t_min, t_max), dot);
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(projection(t0 + i * Radian / ω),
AlmostEquals(series(t0 + i * Radian / ω), 0, 2432));
EXPECT_THAT(projection3(t0 + i * Radian / ω),
RelativeErrorFrom(series(t0 + i * Radian / ω),
AllOf(Gt(3.6e-13), Lt(9.0e-6))));
}
}

6 changes: 4 additions & 2 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -56,8 +56,10 @@ class PoissonSeries {
PoissonSeries<quantities::Primitive<Value, Time>, degree_ + 1, Evaluator>
Primitive() const;

PoissonSeries& operator+=(PoissonSeries const& right);
PoissonSeries& operator-=(PoissonSeries const& right);
template<typename V, int d, template<typename, typename, int> class E>
PoissonSeries& operator+=(PoissonSeries<V, d, E> const& right);
template<typename V, int d, template<typename, typename, int> class E>
PoissonSeries& operator-=(PoissonSeries<V, d, E> const& right);

private:
Instant origin_; // Common to all polynomials.
35 changes: 28 additions & 7 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -133,18 +133,20 @@ PoissonSeries<Value, degree_, Evaluator>::Primitive() const {

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<typename V, int d, template<typename, typename, int> class E>
PoissonSeries<Value, degree_, Evaluator>&
PoissonSeries<Value, degree_, Evaluator>::operator+=(
PoissonSeries const& right) {
PoissonSeries<V, d, E> const& right) {
*this = *this + right;
return *this;
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
template<typename V, int d, template<typename, typename, int> class E>
PoissonSeries<Value, degree_, Evaluator>&
PoissonSeries<Value, degree_, Evaluator>::operator-=(
PoissonSeries const& right) {
PoissonSeries<V, d, E> const& right) {
*this = *this - right;
return *this;
}
@@ -192,10 +194,22 @@ operator+(PoissonSeries<Value, ldegree_, Evaluator> const& left,
? Infinity<AngularFrequency>
: it_right->first;
if (ωl < ωr) {
periodic.insert(periodic.cend(), *it_left);
auto const& polynomials_left = it_left->second;
periodic.emplace_hint(
periodic.cend(),
ωl,
typename Result::Polynomials{
/*sin=*/typename Result::Polynomial(polynomials_left.sin),
/*cos=*/typename Result::Polynomial(polynomials_left.cos)});
++it_left;
} else if (ωr < ωl) {
periodic.insert(periodic.cend(), *it_right);
auto const& polynomials_right = it_right->second;
periodic.emplace_hint(
periodic.cend(),
ωr,
typename Result::Polynomials{
/*sin=*/typename Result::Polynomial(polynomials_right.sin),
/*cos=*/typename Result::Polynomial(polynomials_right.cos)});
++it_right;
} else {
DCHECK_EQ(ωl, ωr);
@@ -233,15 +247,22 @@ operator-(PoissonSeries<Value, ldegree_, Evaluator> const& left,
? Infinity<AngularFrequency>
: it_right->first;
if (ωl < ωr) {
periodic.insert(periodic.cend(), *it_left);
auto const& polynomials_left = it_left->second;
periodic.emplace_hint(
periodic.cend(),
ωl,
typename Result::Polynomials{
/*sin=*/typename Result::Polynomial(polynomials_left.sin),
/*cos=*/typename Result::Polynomial(polynomials_left.cos)});
++it_left;
} else if (ωr < ωl) {
auto const& polynomials_right = it_right->second;
periodic.emplace_hint(
periodic.cend(),
ωr,
typename Result::Polynomials{/*sin=*/-polynomials_right.sin,
/*cos=*/-polynomials_right.cos});
typename Result::Polynomials{
/*sin=*/typename Result::Polynomial(-polynomials_right.sin),
/*cos=*/typename Result::Polynomial(-polynomials_right.cos)});
++it_right;
} else {
DCHECK_EQ(ωl, ωr);