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

Commits on Sep 2, 2020

  1. Copy the full SHA
    abc2fed View commit details
  2. Use the new class.

    pleroy committed Sep 2, 2020
    Copy the full SHA
    1118f79 View commit details
  3. Copy the full SHA
    4c75b37 View commit details
  4. Lint.

    pleroy committed Sep 2, 2020
    Copy the full SHA
    4eb72d7 View commit details

Commits on Sep 5, 2020

  1. Merge pull request #2704 from pleroy/BasisTest

    Move Poisson basis construction to a separate file
    pleroy authored Sep 5, 2020
    Copy the full SHA
    2516413 View commit details
158 changes: 8 additions & 150 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -9,11 +9,10 @@

#include "base/tags.hpp"
#include "geometry/hilbert.hpp"
#include "numerics/poisson_series_basis.hpp"
#include "numerics/root_finders.hpp"
#include "numerics/unbounded_arrays.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/si.hpp"
#include "quantities/traits.hpp"

namespace principia {
namespace numerics {
@@ -23,153 +22,9 @@ namespace internal_frequency_analysis {
using base::uninitialized;
using geometry::Hilbert;
using quantities::Inverse;
using quantities::is_quantity_v;
using quantities::Sqrt;
using quantities::Square;
using quantities::SquareRoot;
namespace si = quantities::si;

// A helper struct for generating the Poisson series tⁿ sin ω t and tⁿ cos ω t.
// d is either 0 (for a 1-dimensional value type) or 0, 1, 2 (for a
// 3-dimensional value type).
template<typename Series, int n, int d>
struct SeriesGenerator {
// The series tⁿ.
static Series Aperiodic(Instant const& origin);
// The series tⁿ sin ω t.
static Series Sin(AngularFrequency const& ω, Instant const& origin);
// The series tⁿ cos ω t.
static Series Cos(AngularFrequency const& ω, Instant const& origin);

private:
// The polynomial tⁿ.
static typename Series::Polynomial Unit(Instant const& origin);
};

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

template<typename Series, int dimension, std::size_t... indices>
struct BasisGenerator<Series, dimension, std::index_sequence<indices...>> {
// Basis of aperiodic terms.
static std::array<Series, dimension * (Series::degree + 1)> Basis(
Instant const& origin);

// Basis of periodic terms.
static std::array<Series, 2 * dimension * (Series::degree + 1)> Basis(
AngularFrequency const& ω,
Instant const& origin);
};


template<typename Series, int n, int d>
Series SeriesGenerator<Series, n, d>::Aperiodic(Instant const& origin) {
return Series(Unit(origin), {});
}

template<typename Series, int n, int d>
Series SeriesGenerator<Series, n, d>::Sin(AngularFrequency const& ω,
Instant const& origin) {
typename Series::Polynomial::Coefficients const zeros;
typename Series::Polynomial const zero{zeros, origin};
return Series(zero,
{{ω,
{/*sin=*/Unit(origin),
/*cos=*/zero}}});
}

template<typename Series, int n, int d>
Series SeriesGenerator<Series, n, d>::Cos(AngularFrequency const& ω,
Instant const& origin) {
typename Series::Polynomial::Coefficients const zeros;
typename Series::Polynomial const zero{zeros, origin};
return Series(zero,
{{ω,
{/*sin=*/zero,
/*cos=*/Unit(origin)}}});
}

template<typename Series, int n, int d>
typename Series::Polynomial SeriesGenerator<Series, n, d>::Unit(
Instant const& origin) {
typename Series::Polynomial::Coefficients coefficients;
using Coefficient =
std::tuple_element_t<n, typename Series::Polynomial::Coefficients>;
Coefficient& coefficient = std::get<n>(coefficients);
if constexpr (is_quantity_v<Coefficient>) {
coefficient = si::Unit<Coefficient>;
} else {
using Scalar = typename Hilbert<Coefficient>::NormType;
if constexpr (d == 0) {
coefficient = Coefficient({si::Unit<Scalar>, Scalar{}, Scalar{}});
} else if constexpr (d == 1) {
coefficient = Coefficient({Scalar{}, si::Unit<Scalar>, Scalar{}});
} else if constexpr (d == 2) {
coefficient = Coefficient({Scalar{}, Scalar{}, si::Unit<Scalar>});
}
}
return typename Series::Polynomial(coefficients, origin);
}

// 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>
std::array<Series, dimension * (Series::degree + 1)>
BasisGenerator<Series, dimension, 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)>
BasisGenerator<Series, dimension, 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 = {
SeriesGenerator<Series, indices / dimension, indices % dimension>::Sin(
ω, origin)...,
SeriesGenerator<Series, indices / dimension, indices % dimension>::Cos(
ω, origin)...};

// Order all_series by repeatedly swapping its elements.
if constexpr (all_series.size() > 2) {
// The index of this array is the current index of a series in all_series.
// The value is the index of the final resting place of that series in
// all_series. The elements at indices 0 and
// 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)
? 2 * i
: 2 * (i - dimension * (Series::degree + 1)) + 1;
}
for (int i = 1; i < all_series.size() - 1;) {
// Swap the series currently at index i to its final resting place.
// Iterate until the series at index i is at its final resting place
// (i.e., after we have executed an entire cycle of the permutation).
// Then move to the next series.
if (i == permutation[i]) {
++i;
} else {
int const j = permutation[i];
std::swap(all_series[i], all_series[j]);
std::swap(permutation[i], permutation[j]);
}
}
}
return all_series;
}


template<typename Function,
int wdegree_,
@@ -253,12 +108,14 @@ IncrementalProjection(Function const& function,
int basis_size;
if (ω.value() == AngularFrequency{}) {
auto const ω_basis =
BasisGenerator<Series, Hilbert<Value>::dimension>::Basis(t0);
PoissonSeriesBasisGenerator<Series, Hilbert<Value>::dimension>::Basis(
t0);
basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
} else {
auto const ω_basis =
BasisGenerator<Series, Hilbert<Value>::dimension>::Basis(ω.value(), t0);
PoissonSeriesBasisGenerator<Series, Hilbert<Value>::dimension>::Basis(
ω.value(), t0);
basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
}
@@ -371,12 +228,13 @@ IncrementalProjection(Function const& function,
int ω_basis_size;
if (ω.value() == AngularFrequency{}) {
auto const ω_basis =
BasisGenerator<Series, Hilbert<Value>::dimension>::Basis(t0);
PoissonSeriesBasisGenerator<Series, Hilbert<Value>::dimension>::Basis(
t0);
ω_basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
} else {
auto const ω_basis =
BasisGenerator<Series, Hilbert<Value>::dimension>::Basis(
PoissonSeriesBasisGenerator<Series, Hilbert<Value>::dimension>::Basis(
ω.value(), t0);
ω_basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
3 changes: 3 additions & 0 deletions numerics/numerics.vcxproj
Original file line number Diff line number Diff line change
@@ -47,6 +47,8 @@
<ClInclude Include="pid.hpp" />
<ClInclude Include="pid_body.hpp" />
<ClInclude Include="poisson_series.hpp" />
<ClInclude Include="poisson_series_basis.hpp" />
<ClInclude Include="poisson_series_basis_body.hpp" />
<ClInclude Include="poisson_series_body.hpp" />
<ClInclude Include="polynomial.hpp" />
<ClInclude Include="polynomial_body.hpp" />
@@ -84,6 +86,7 @@
<ClCompile Include="max_abs_normalized_associated_legendre_functions_test.cc" />
<ClCompile Include="newhall_test.cpp" />
<ClCompile Include="pid_test.cpp" />
<ClCompile Include="poisson_series_basis_test.cpp" />
<ClCompile Include="poisson_series_test.cpp" />
<ClCompile Include="polynomial_evaluators_test.cpp" />
<ClCompile Include="polynomial_test.cpp" />
9 changes: 9 additions & 0 deletions numerics/numerics.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -155,6 +155,12 @@
<ClInclude Include="unbounded_arrays_body.hpp">
<Filter>Source Files</Filter>
</ClInclude>
<ClInclude Include="poisson_series_basis.hpp">
<Filter>Header Files</Filter>
</ClInclude>
<ClInclude Include="poisson_series_basis_body.hpp">
<Filter>Source Files</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<ClCompile Include="чебышёв_series_test.cpp">
@@ -238,6 +244,9 @@
<ClCompile Include="unbounded_arrays_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
<ClCompile Include="poisson_series_basis_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<Text Include="xgscd.proto.txt">
48 changes: 48 additions & 0 deletions numerics/poisson_series_basis.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
#pragma once

#include <array>
#include <utility>

#include "geometry/named_quantities.hpp"
#include "numerics/poisson_series.hpp"
#include "quantities/named_quantities.hpp"

namespace principia {
namespace numerics {
namespace internal_poisson_series_basis {

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.
template<typename Series,
int dimension,
typename = std::make_index_sequence<dimension * (Series::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>
struct PoissonSeriesBasisGenerator<Series,
dimension,
std::index_sequence<indices...>> {
// Basis of aperiodic terms.
static std::array<Series, dimension * (Series::degree + 1)> Basis(
Instant const& origin);

// Basis of periodic terms.
static std::array<Series, 2 * dimension * (Series::degree + 1)> Basis(
AngularFrequency const& ω,
Instant const& origin);
};

} // namespace internal_poisson_series_basis

using internal_poisson_series_basis::PoissonSeriesBasisGenerator;

} // namespace numerics
} // namespace principia

#include "numerics/poisson_series_basis_body.hpp"
Loading