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

Commits on Nov 8, 2020

  1. orthogonal decomposition

    eggrobin committed Nov 8, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    bcbb78b View commit details
  2. space

    eggrobin committed Nov 8, 2020
    Copy the full SHA
    581b537 View commit details

Commits on Nov 15, 2020

  1. Copy the full SHA
    18ba8fb View commit details

Commits on Nov 22, 2020

  1. merge

    eggrobin committed Nov 22, 2020
    Copy the full SHA
    0aa6403 View commit details
  2. make it build

    eggrobin committed Nov 22, 2020
    Copy the full SHA
    fe35575 View commit details
  3. lint

    eggrobin committed Nov 22, 2020
    Copy the full SHA
    386311f View commit details

Commits on Nov 25, 2020

  1. after pleroy’s review

    eggrobin committed Nov 25, 2020
    Copy the full SHA
    d2d4a04 View commit details
  2. Copy the full SHA
    35ab830 View commit details

Commits on Nov 26, 2020

  1. Merge pull request #2783 from eggrobin/orthogonal-decomposition

    Orthogonal decomposition
    eggrobin authored Nov 26, 2020
    Copy the full SHA
    d797277 View commit details
Showing with 163 additions and 5 deletions.
  1. +32 −2 numerics/frequency_analysis_body.hpp
  2. +36 −0 numerics/poisson_series_basis.hpp
  3. +60 −3 numerics/poisson_series_basis_body.hpp
  4. +35 −0 numerics/poisson_series_basis_test.cpp
34 changes: 32 additions & 2 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -112,20 +112,36 @@ IncrementalProjection(Function const& function,
CHECK(ω.has_value());

std::vector<Series> basis;
// The Poisson series basis[k] belongs to the subspace basis_subspaces[k];
// this remains true after orthonormalization, i.e., q[k] belongs to the
// subspace basis_subspaces[k] below.
std::vector<PoissonSeriesSubspace> basis_subspaces;

int basis_size;
// TODO(phl): This is replicated below.
if (ω.value() == AngularFrequency{}) {
auto const ω_basis =
PoissonSeriesBasisGenerator<Series,
aperiodic_degree>::Basis(t0);
auto const ω_basis_subspaces =
PoissonSeriesBasisGenerator<Series, aperiodic_degree>::Subspaces(t0);
basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
std::move(ω_basis_subspaces.begin(),
ω_basis_subspaces.end(),
std::back_inserter(basis_subspaces));
} else {
auto const ω_basis =
PoissonSeriesBasisGenerator<Series,
periodic_degree>::Basis(ω.value(), t0);
auto const ω_basis_subspaces =
PoissonSeriesBasisGenerator<Series, periodic_degree>::Subspaces(
ω.value(), t0);
basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
std::move(ω_basis_subspaces.begin(),
ω_basis_subspaces.end(),
std::back_inserter(basis_subspaces));
}

// This is logically Q in the QR decomposition of basis.
@@ -147,8 +163,11 @@ IncrementalProjection(Function const& function,
for (int m = m_begin; m < basis_size; ++m) {
auto aₘ⁽ᵏ⁾ = basis[m];
for (int k = 0; k < m; ++k) {
auto const rₖₘ = InnerProduct(q[k], aₘ⁽ᵏ⁾, weight, t_min, t_max);
aₘ⁽ᵏ⁾ -= rₖₘ * q[k];
if (!PoissonSeriesSubspace::orthogonal(basis_subspaces[k],
basis_subspaces[m])) {
auto const rₖₘ = InnerProduct(q[k], aₘ⁽ᵏ⁾, weight, t_min, t_max);
aₘ⁽ᵏ⁾ -= rₖₘ * q[k];
}
}

auto const rₘₘ = aₘ⁽ᵏ⁾.Norm(weight, t_min, t_max);
@@ -171,14 +190,25 @@ IncrementalProjection(Function const& function,
auto const ω_basis =
PoissonSeriesBasisGenerator<Series,
aperiodic_degree>::Basis(t0);
auto const ω_basis_subspaces =
PoissonSeriesBasisGenerator<Series, aperiodic_degree>::Subspaces(t0);
ω_basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
std::move(ω_basis_subspaces.begin(),
ω_basis_subspaces.end(),
std::back_inserter(basis_subspaces));
} else {
auto const ω_basis =
PoissonSeriesBasisGenerator<Series,
periodic_degree>::Basis(ω.value(), t0);
auto const ω_basis_subspaces =
PoissonSeriesBasisGenerator<Series, periodic_degree>::Subspaces(
ω.value(), t0);
ω_basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
std::move(ω_basis_subspaces.begin(),
ω_basis_subspaces.end(),
std::back_inserter(basis_subspaces));
}
m_begin = basis_size;
basis_size += ω_basis_size;
36 changes: 36 additions & 0 deletions numerics/poisson_series_basis.hpp
Original file line number Diff line number Diff line change
@@ -17,6 +17,35 @@ using geometry::Hilbert;
using geometry::Instant;
using quantities::AngularFrequency;

// A |PoissonSeriesSubspace| represents a linear subspace of the space of
// Poisson series.
// The type |PoissonSeriesSubspace| defines an orthogonal decomposition of the
// space of Poisson series, i.e., The space of Poisson series is the orthogonal
// sum of all values of |PoissonSeriesSubspace|.
class PoissonSeriesSubspace {
public:
// Whether the subspaces |v| and |w| are orthogonal.
// TODO(egg): When we take parity into account, orthogonality will be defined
// with respect to an inner product over an interval centred on the origin of
// the Poisson series, with an even apodization.
static bool orthogonal(PoissonSeriesSubspace v, PoissonSeriesSubspace w);

private:
enum class Coordinate { X = 0, Y = 1, Z = 2 };
enum class Parity { Even = 0, Odd = 1 };

PoissonSeriesSubspace(Coordinate coordinate, Parity parity);

Coordinate coordinate_;
// The parity of the Poisson series is defined with respect to its origin.
Parity parity_;

template<typename Series, int degree, int dimension, typename>
friend struct AperiodicSeriesGenerator;
template<typename Series, int degree, int dimension, typename>
friend struct PeriodicSeriesGenerator;
};

// A generator for the Кудрявцев basis, i.e., functions of the
// form tⁿ sin ω t and tⁿ cos ω t properly ordered. |degree| is the maximum
// degree of tⁿ.
@@ -29,16 +58,23 @@ class PoissonSeriesBasisGenerator {
// Basis of aperiodic terms.
static std::array<Series, dimension * (degree + 1)> Basis(
Instant const& origin);
// The subspaces to which the above terms belong.
static std::array<PoissonSeriesSubspace, dimension * (degree + 1)> Subspaces(
Instant const& origin);

// Basis of periodic terms.
static std::array<Series, 2 * dimension * (degree + 1)> Basis(
AngularFrequency const& ω,
Instant const& origin);
// The subspaces to which the above terms belong.
static std::array<PoissonSeriesSubspace, 2 * dimension * (degree + 1)>
Subspaces(AngularFrequency const& ω, Instant const& origin);
};

} // namespace internal_poisson_series_basis

using internal_poisson_series_basis::PoissonSeriesBasisGenerator;
using internal_poisson_series_basis::PoissonSeriesSubspace;

} // namespace numerics
} // namespace principia
63 changes: 60 additions & 3 deletions numerics/poisson_series_basis_body.hpp
Original file line number Diff line number Diff line change
@@ -136,6 +136,14 @@ Polynomials PolynomialGenerator<Polynomial, dimension>::UnitPolynomials(
}
}

inline bool PoissonSeriesSubspace::orthogonal(PoissonSeriesSubspace const v,
PoissonSeriesSubspace const w) {
return v.coordinate_ != w.coordinate_;
}

inline PoissonSeriesSubspace::PoissonSeriesSubspace(Coordinate coordinate,
Parity parity)
: coordinate_(coordinate), parity_(parity) {}

// A helper to generate aperiodic series, i.e., series of the form tⁿ along some
// dimension. Note how a single index_sequence is generated that covers the
@@ -151,6 +159,8 @@ struct AperiodicSeriesGenerator<Series,
std::index_sequence<indices...>> {
static std::array<Series, sizeof...(indices)> BasisElements(
Instant const& origin);
static std::array<PoissonSeriesSubspace, sizeof...(indices)> Subspaces(
Instant const& origin);
};

template<typename Series, int degree, int dimension, std::size_t... indices>
@@ -164,6 +174,18 @@ std::array<Series, sizeof...(indices)> AperiodicSeriesGenerator<
{}))...};
}

template<typename Series, int degree, int dimension, std::size_t... indices>
std::array<PoissonSeriesSubspace, sizeof...(indices)> AperiodicSeriesGenerator<
Series, degree, dimension,
std::index_sequence<indices...>>::Subspaces(Instant const& origin) {
return {PoissonSeriesSubspace{
static_cast<PoissonSeriesSubspace::Coordinate>(indices % dimension),
// The degree of the ith polynomial is i / dimension, so its parity is
// (i / dimension) % 2.
static_cast<PoissonSeriesSubspace::Parity>((indices / dimension) %
2)}...};
}


// A helper to generate periodic series, i.e., series of the form tⁿ sin ω t
// and tⁿ cos ω t along some dimension. Note how a single index_sequence is
@@ -174,12 +196,14 @@ template<typename Series,
struct PeriodicSeriesGenerator;

template<typename Series, int degree, int dimension, std::size_t... indices>
struct PeriodicSeriesGenerator<Series,
degree, dimension,
std::index_sequence<indices...>> {
struct PeriodicSeriesGenerator<Series, degree, dimension,
std::index_sequence<indices...>> {
static std::array<Series, sizeof...(indices)> BasisElements(
AngularFrequency const& ω,
Instant const& origin);
static std::array<PoissonSeriesSubspace, sizeof...(indices)> Subspaces(
AngularFrequency const& ω,
Instant const& origin);
};

template<typename Series, int degree, int dimension, std::size_t... indices>
@@ -197,6 +221,23 @@ std::array<Series, sizeof...(indices)> PeriodicSeriesGenerator<
origin)}})...};
}

template<typename Series, int degree, int dimension, std::size_t... indices>
std::array<PoissonSeriesSubspace, sizeof...(indices)> PeriodicSeriesGenerator<
Series, degree, dimension,
std::index_sequence<indices...>>::Subspaces(AngularFrequency const& ω,
Instant const& origin) {
return {
PoissonSeriesSubspace{
static_cast<PoissonSeriesSubspace::Coordinate>(indices % dimension),
// The parity of the trigonometric factors of the ith basis element is
// (i / dimension) % 2; the degrees its polynomial factor is
// i / (2 * dimension), so that the overall parity of that basis
// element is (i / dimension + i / (2 * dimension)) % 2, which
// simplifies to (3 * i / (2 * dimension)) % 2.
static_cast<PoissonSeriesSubspace::Parity>(
(3 * indices / (2 * dimension)) % 2)}...};
}


template<typename Series, int degree>
auto PoissonSeriesBasisGenerator<Series, degree>::Basis(Instant const& origin)
@@ -205,6 +246,13 @@ auto PoissonSeriesBasisGenerator<Series, degree>::Basis(Instant const& origin)
BasisElements(origin);
}

template<typename Series, int degree>
auto PoissonSeriesBasisGenerator<Series, degree>::Subspaces(
Instant const& origin)
-> std::array<PoissonSeriesSubspace, dimension*(degree + 1)> {
return AperiodicSeriesGenerator<Series, degree, dimension>::Subspaces(origin);
}

template<typename Series, int degree>
auto PoissonSeriesBasisGenerator<Series, degree>::Basis(
AngularFrequency const& ω,
@@ -213,6 +261,15 @@ auto PoissonSeriesBasisGenerator<Series, degree>::Basis(
BasisElements(ω, origin);
}

template<typename Series, int degree>
auto PoissonSeriesBasisGenerator<Series, degree>::Subspaces(
AngularFrequency const& ω,
Instant const& origin)
-> std::array<PoissonSeriesSubspace, 2 * dimension*(degree + 1)> {
return PeriodicSeriesGenerator<Series, degree, dimension>::Subspaces(ω,
origin);
}

} // namespace internal_poisson_series_basis
} // namespace numerics
} // namespace principia
35 changes: 35 additions & 0 deletions numerics/poisson_series_basis_test.cpp
Original file line number Diff line number Diff line change
@@ -63,6 +63,9 @@ TEST_F(PoissonSeriesBasisTest, AperiodicVector) {
auto const aperiodic = PoissonSeriesBasisGenerator<
Series3,
/*degree=*/3>::Basis(t0_);
auto const aperiodic_subspaces =
PoissonSeriesBasisGenerator<Series3,
/*degree=*/3>::Subspaces(t0_);
EXPECT_EQ(12, aperiodic.size());

Instant const t1 = t0_ + 2 * Second;
@@ -74,13 +77,27 @@ TEST_F(PoissonSeriesBasisTest, AperiodicVector) {
EXPECT_EQ(Displacement<World>({0 * Metre, 0 * Metre, 1 * Metre}),
aperiodic[2](t1));

EXPECT_FALSE(PoissonSeriesSubspace::orthogonal(aperiodic_subspaces[0],
aperiodic_subspaces[0]));
EXPECT_TRUE(PoissonSeriesSubspace::orthogonal(aperiodic_subspaces[0],
aperiodic_subspaces[1]));
EXPECT_TRUE(PoissonSeriesSubspace::orthogonal(aperiodic_subspaces[0],
aperiodic_subspaces[2]));

EXPECT_EQ(Displacement<World>({2 * Metre, 0 * Metre, 0 * Metre}),
aperiodic[3](t1));
EXPECT_EQ(Displacement<World>({0 * Metre, 2 * Metre, 0 * Metre}),
aperiodic[4](t1));
EXPECT_EQ(Displacement<World>({0 * Metre, 0 * Metre, 2 * Metre}),
aperiodic[5](t1));

EXPECT_FALSE(PoissonSeriesSubspace::orthogonal(aperiodic_subspaces[0],
aperiodic_subspaces[3]));
EXPECT_TRUE(PoissonSeriesSubspace::orthogonal(aperiodic_subspaces[0],
aperiodic_subspaces[4]));
EXPECT_TRUE(PoissonSeriesSubspace::orthogonal(aperiodic_subspaces[0],
aperiodic_subspaces[5]));

EXPECT_EQ(Displacement<World>({4 * Metre, 0 * Metre, 0 * Metre}),
aperiodic[6](t1));
EXPECT_EQ(Displacement<World>({0 * Metre, 4 * Metre, 0 * Metre}),
@@ -120,6 +137,9 @@ TEST_F(PoissonSeriesBasisTest, PeriodicVector) {
auto const periodic = PoissonSeriesBasisGenerator<
Series3,
/*degree=*/3>::Basis(ω, t0_);
auto const periodic_subspaces =
PoissonSeriesBasisGenerator<Series3,
/*degree=*/3>::Subspaces(ω, t0_);
EXPECT_EQ(24, periodic.size());

Instant const t1 = t0_ + 2 * Second;
@@ -136,6 +156,14 @@ TEST_F(PoissonSeriesBasisTest, PeriodicVector) {
periodic[2](t1),
AlmostEquals(
Displacement<World>({0 * Metre, 0 * Metre, 0.5 * Metre}), 1));

EXPECT_FALSE(PoissonSeriesSubspace::orthogonal(periodic_subspaces[0],
periodic_subspaces[0]));
EXPECT_TRUE(PoissonSeriesSubspace::orthogonal(periodic_subspaces[0],
periodic_subspaces[1]));
EXPECT_TRUE(PoissonSeriesSubspace::orthogonal(periodic_subspaces[0],
periodic_subspaces[2]));

EXPECT_THAT(
periodic[3](t1),
AlmostEquals(
@@ -149,6 +177,13 @@ TEST_F(PoissonSeriesBasisTest, PeriodicVector) {
AlmostEquals(
Displacement<World>({0 * Metre, 0 * Metre, Sqrt(3) / 2 * Metre}), 0));

EXPECT_FALSE(PoissonSeriesSubspace::orthogonal(periodic_subspaces[0],
periodic_subspaces[3]));
EXPECT_TRUE(PoissonSeriesSubspace::orthogonal(periodic_subspaces[0],
periodic_subspaces[4]));
EXPECT_TRUE(PoissonSeriesSubspace::orthogonal(periodic_subspaces[0],
periodic_subspaces[5]));

EXPECT_THAT(
periodic[6](t1),
AlmostEquals(