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

Commits on Oct 23, 2020

  1. Add support for limiting the degree of polynomials produced by the Po…

    …issonSeriesBasisGenerator.
    pleroy committed Oct 23, 2020
    Copy the full SHA
    bd331ec View commit details
  2. Copy the full SHA
    a83f695 View commit details
  3. Lint.

    pleroy committed Oct 23, 2020
    Copy the full SHA
    9e7abae View commit details
  4. Merge pull request #2768 from pleroy/Base

    Add support for limiting the degree of polynomials produced by the basis generator
    pleroy authored Oct 23, 2020
    Copy the full SHA
    f9c97d3 View commit details
Showing with 176 additions and 48 deletions.
  1. +13 −10 numerics/frequency_analysis.hpp
  2. +30 −22 numerics/frequency_analysis_body.hpp
  3. +8 −5 numerics/poisson_series_basis.hpp
  4. +9 −9 numerics/poisson_series_basis_body.hpp
  5. +116 −2 numerics/poisson_series_basis_test.cpp
23 changes: 13 additions & 10 deletions numerics/frequency_analysis.hpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@

#pragma once

#include <algorithm>
#include <type_traits>

#include "geometry/interval.hpp"
@@ -25,26 +26,26 @@ using quantities::Time;
// obtained using FFT). The |Function| must have a member |FourierTransform|
// that returns its spectrum. See [Cha95].
template<typename Function,
int wdegree_,
int wdegree,
template<typename, typename, int> class Evaluator>
AngularFrequency PreciseMode(
Interval<AngularFrequency> const& fft_mode,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight);
PoissonSeries<double, wdegree, Evaluator> const& weight);

// In the projection functions the |Function| must have an |InnerProduct| with
// |PoissonSeries| or |PiecewisePoissonSeries|.

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

@@ -66,14 +67,16 @@ Projection(Function const& function,
// to |IncrementalProjection|.
// If the calculator cannot find a suitable frequency, or if it wants to stop
// the algorithm, it does so by returning std::nullopt.
template<int degree_,
template<int aperiodic_degree, int periodic_degree = aperiodic_degree,
typename Function,
typename AngularFrequencyCalculator, int wdegree_,
typename AngularFrequencyCalculator, int wdegree,
template<typename, typename, int> class Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>, degree_, Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>,
std::max(aperiodic_degree, periodic_degree),
Evaluator>
IncrementalProjection(Function const& function,
AngularFrequencyCalculator const& calculator,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
PoissonSeries<double, wdegree, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max);

52 changes: 30 additions & 22 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -29,12 +29,12 @@ using quantities::Square;
using quantities::SquareRoot;

template<typename Function,
int wdegree_,
int wdegree,
template<typename, typename, int> class Evaluator>
AngularFrequency PreciseMode(
Interval<AngularFrequency> const& fft_mode,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight) {
PoissonSeries<double, wdegree, Evaluator> const& weight) {
auto const weighted_function = weight * function;
auto const weighted_function_spectrum = weighted_function.FourierTransform();

@@ -49,14 +49,14 @@ AngularFrequency PreciseMode(
std::greater<>());
}

template<int degree_,
template<int degree,
typename Function,
int wdegree_,
int wdegree,
template<typename, typename, int> class Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>, degree_, Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>, degree, Evaluator>
Projection(Function const& function,
AngularFrequency const& ω,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
PoissonSeries<double, wdegree, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
std::optional<AngularFrequency> optional_ω = ω;
@@ -68,26 +68,30 @@ Projection(Function const& function,
return result;
};

return IncrementalProjection<degree_>(function,
return IncrementalProjection<degree>(function,
angular_frequency_calculator,
weight,
t_min, t_max);
}

template<int degree_,
template<int aperiodic_degree, int periodic_degree,
typename Function,
typename AngularFrequencyCalculator, int wdegree_,
typename AngularFrequencyCalculator, int wdegree,
template<typename, typename, int> class Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>, degree_, Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>,
std::max(aperiodic_degree, periodic_degree),
Evaluator>
IncrementalProjection(Function const& function,
AngularFrequencyCalculator const& calculator,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
PoissonSeries<double, wdegree, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
constexpr int degree = std::max(aperiodic_degree, periodic_degree);

using Value = std::invoke_result_t<Function, Instant>;
using Norm = typename Hilbert<Value>::NormType;
using Normalized = typename Hilbert<Value>::NormalizedType;
using Series = PoissonSeries<Value, degree_, Evaluator>;
using Series = PoissonSeries<Value, degree, Evaluator>;

// This code follows [Kud07], section 2. Our indices start at 0, unlike those
// of Кудрявцев which start at 1.
@@ -102,28 +106,30 @@ IncrementalProjection(Function const& function,
int basis_size;
if (ω.value() == AngularFrequency{}) {
auto const ω_basis =
PoissonSeriesBasisGenerator<Series, Hilbert<Value>::dimension>::Basis(
t0);
PoissonSeriesBasisGenerator<Series,
Hilbert<Value>::dimension,
aperiodic_degree>::Basis(t0);
basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
} else {
auto const ω_basis =
PoissonSeriesBasisGenerator<Series, Hilbert<Value>::dimension>::Basis(
ω.value(), t0);
PoissonSeriesBasisGenerator<Series,
Hilbert<Value>::dimension,
periodic_degree>::Basis(ω.value(), t0);
basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
}

// This is logically Q in the QR decomposition of basis.
std::vector<PoissonSeries<Normalized, degree_, Evaluator>> q;
std::vector<PoissonSeries<Normalized, degree, Evaluator>> q;

auto const a₀ = basis[0];
auto const r₀₀ = a₀.Norm(weight, t_min, t_max);
q.push_back(a₀ / r₀₀);

auto const A₀ = InnerProduct(function, q[0], weight, t_min, t_max);

PoissonSeries<Value, degree_, Evaluator> F = A₀ * q[0];
PoissonSeries<Value, degree, Evaluator> F = A₀ * q[0];
auto f = function - F;

int m_begin = 1;
@@ -153,14 +159,16 @@ IncrementalProjection(Function const& function,
int ω_basis_size;
if (ω.value() == AngularFrequency{}) {
auto const ω_basis =
PoissonSeriesBasisGenerator<Series, Hilbert<Value>::dimension>::Basis(
t0);
PoissonSeriesBasisGenerator<Series,
Hilbert<Value>::dimension,
aperiodic_degree>::Basis(t0);
ω_basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
} else {
auto const ω_basis =
PoissonSeriesBasisGenerator<Series, Hilbert<Value>::dimension>::Basis(
ω.value(), t0);
PoissonSeriesBasisGenerator<Series,
Hilbert<Value>::dimension,
periodic_degree>::Basis(ω.value(), t0);
ω_basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
}
13 changes: 8 additions & 5 deletions numerics/poisson_series_basis.hpp
Original file line number Diff line number Diff line change
@@ -15,25 +15,28 @@ using geometry::Instant;
using quantities::AngularFrequency;

// A helper struct for generating the Кудрявцев basis, i.e., functions of the
// form tⁿ sin ω t and tⁿ cos ω t properly ordered.
// form tⁿ sin ω t and tⁿ cos ω t properly ordered. |dimension| is the number
// of multivector dimensions to produce. |degree| is the maximum degree of tⁿ.
template<typename Series,
int dimension,
typename = std::make_index_sequence<dimension * (Series::degree + 1)>>
int degree,
typename = std::make_index_sequence<dimension * (degree + 1)>>
struct PoissonSeriesBasisGenerator;

// In this template, the indices encode the degree and the dimension of the
// basis term so that, in the terminology of SeriesGenerator, n (the degree) is
// indices / dimension and d (the dimension index) is indices % dimension.
template<typename Series, int dimension, std::size_t... indices>
template<typename Series, int dimension, int degree, std::size_t... indices>
struct PoissonSeriesBasisGenerator<Series,
dimension,
degree,
std::index_sequence<indices...>> {
// Basis of aperiodic terms.
static std::array<Series, dimension * (Series::degree + 1)> Basis(
static std::array<Series, dimension * (degree + 1)> Basis(
Instant const& origin);

// Basis of periodic terms.
static std::array<Series, 2 * dimension * (Series::degree + 1)> Basis(
static std::array<Series, 2 * dimension * (degree + 1)> Basis(
AngularFrequency const& ω,
Instant const& origin);
};
18 changes: 9 additions & 9 deletions numerics/poisson_series_basis_body.hpp
Original file line number Diff line number Diff line change
@@ -82,25 +82,25 @@ typename Series::Polynomial SeriesGenerator<Series, n, d>::Unit(
return typename Series::Polynomial(coefficients, origin);
}

template<typename Series, int dimension, std::size_t... indices>
std::array<Series, dimension*(Series::degree + 1)> PoissonSeriesBasisGenerator<
Series, dimension,
template<typename Series, int dimension, int degree, std::size_t... indices>
std::array<Series, dimension * (degree + 1)> PoissonSeriesBasisGenerator<
Series, dimension, degree,
std::index_sequence<indices...>>::Basis(Instant const& origin) {
return {SeriesGenerator<Series, indices / dimension, indices % dimension>::
Aperiodic(origin)...};
}

template<typename Series, int dimension, std::size_t... indices>
std::array<Series, 2 * dimension*(Series::degree + 1)>
template<typename Series, int dimension, int degree, std::size_t... indices>
std::array<Series, 2 * dimension * (degree + 1)>
PoissonSeriesBasisGenerator<
Series, dimension,
Series, dimension, degree,
std::index_sequence<indices...>>::Basis(AngularFrequency const& ω,
Instant const& origin) {
// This has the elements {Sin(ωt), t Sin(ωt), t² Sin(ωt), ..., Cos(ωt), ...}
// in the scalar case and {x Sin(ωt), y Sin(ωt), z Sin(ωt), x t Sin(ωt), ...}
// in the vector case. This is not the order we want (we want lower-degree
// polynomials first) so we'll need to reorder the terms.
std::array<Series, 2 * dimension * (Series::degree + 1)> all_series = {
std::array<Series, 2 * dimension * (degree + 1)> all_series = {
SeriesGenerator<Series, indices / dimension, indices % dimension>::Sin(
ω, origin)...,
SeriesGenerator<Series, indices / dimension, indices % dimension>::Cos(
@@ -114,9 +114,9 @@ PoissonSeriesBasisGenerator<
// 2 * dimension * (Series::degree + 1) are unused.
std::array<int, all_series.size()> permutation;
for (int i = 1; i < all_series.size() - 1; ++i) {
permutation[i] = i < dimension * (Series::degree + 1)
permutation[i] = i < dimension * (degree + 1)
? 2 * i
: 2 * (i - dimension * (Series::degree + 1)) + 1;
: 2 * (i - dimension * (degree + 1)) + 1;
}
for (int i = 1; i < all_series.size() - 1;) {
// Swap the series currently at index i to its final resting place.
Loading