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: 22a3d91b6db3
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 3a9bf8153a16
Choose a head ref

Commits on Aug 29, 2020

  1. Copy the full SHA
    3d3cd61 View commit details
  2. Copy the full SHA
    4629d27 View commit details
  3. Funky product of polynomials.

    pleroy committed Aug 29, 2020
    Copy the full SHA
    79623a7 View commit details
  4. Copy the full SHA
    92379dd View commit details
  5. Copy the full SHA
    96ae92f View commit details
  6. Remove the operators.

    pleroy committed Aug 29, 2020
    Copy the full SHA
    a3c8b66 View commit details
  7. Copy the full SHA
    3d0eb82 View commit details
  8. Back to Kudryavtsev.

    pleroy committed Aug 29, 2020
    Copy the full SHA
    b757175 View commit details

Commits on Aug 30, 2020

  1. Merge.

    pleroy committed Aug 30, 2020
    Copy the full SHA
    d6b04ac View commit details
  2. Copy the full SHA
    cfc67aa View commit details
  3. Copy the full SHA
    427aa33 View commit details
  4. Copy the full SHA
    0e313d9 View commit details
  5. Copy the full SHA
    ed8de51 View commit details
  6. A test for is_quantity_v.

    pleroy committed Aug 30, 2020
    Copy the full SHA
    6a85c5e View commit details
  7. Copy the full SHA
    4c58dbf View commit details
  8. Copy the full SHA
    9c24584 View commit details
  9. Merge.

    pleroy committed Aug 30, 2020
    Copy the full SHA
    6334575 View commit details
  10. Copy the full SHA
    92ad860 View commit details

Commits on Aug 31, 2020

  1. Fix a few types.

    pleroy committed Aug 31, 2020
    Copy the full SHA
    2bba773 View commit details
  2. Copy the full SHA
    27c6bea View commit details

Commits on Sep 1, 2020

  1. Hilbert space dimension.

    pleroy committed Sep 1, 2020
    Copy the full SHA
    7c32959 View commit details
  2. Copy the full SHA
    dc99206 View commit details
  3. Polynomial output.

    pleroy committed Sep 1, 2020
    Copy the full SHA
    d40b620 View commit details

Commits on Sep 2, 2020

  1. Merge.

    pleroy committed Sep 2, 2020
    Copy the full SHA
    ce538c9 View commit details
  2. Test passing.

    pleroy committed Sep 2, 2020
    Copy the full SHA
    d576c50 View commit details
  3. Copy the full SHA
    d20b096 View commit details
  4. Copy the full SHA
    b8acc2b View commit details
  5. Fix types.

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

Commits on Sep 5, 2020

  1. Copy the full SHA
    12b005c View commit details
  2. After egg's review.

    pleroy committed Sep 5, 2020
    Copy the full SHA
    982da0b View commit details
  3. Merge pull request #2705 from pleroy/VectorKudryavtsev

    Кудрявцев projection for vector-valued functions
    pleroy authored Sep 5, 2020
    Copy the full SHA
    3a9bf81 View commit details
Showing with 139 additions and 38 deletions.
  1. +4 −0 geometry/hilbert.hpp
  2. +15 −9 numerics/frequency_analysis_body.hpp
  3. +98 −11 numerics/frequency_analysis_test.cpp
  4. +7 −7 numerics/poisson_series.hpp
  5. +13 −11 numerics/poisson_series_body.hpp
  6. +2 −0 quantities/named_quantities.hpp
4 changes: 4 additions & 0 deletions geometry/hilbert.hpp
Original file line number Diff line number Diff line change
@@ -46,6 +46,8 @@ struct Hilbert<T, T, std::enable_if_t<is_quantity_v<T>>>

using NormType = T;
static NormType Norm(T const& t);

using NormalizedType = double;
};

template<typename T1, typename T2>
@@ -74,6 +76,8 @@ struct Hilbert<T, T,

using NormType = decltype(std::declval<T>().Norm());
static NormType Norm(T const& t);

using NormalizedType = decltype(std::declval<T>() / NormType());
};

} // namespace internal_hilbert
24 changes: 15 additions & 9 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -8,6 +8,7 @@
#include <vector>

#include "base/tags.hpp"
#include "geometry/grassmann.hpp"
#include "geometry/hilbert.hpp"
#include "numerics/poisson_series_basis.hpp"
#include "numerics/root_finders.hpp"
@@ -21,6 +22,7 @@ namespace internal_frequency_analysis {

using base::uninitialized;
using geometry::Hilbert;
using geometry::Vector;
using quantities::Inverse;
using quantities::Sqrt;
using quantities::Square;
@@ -93,6 +95,9 @@ IncrementalProjection(Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Dot const& dot) {
using Value = std::invoke_result_t<Function, Instant>;
using Norm = typename Hilbert<Value>::NormType;
using Norm² = typename Hilbert<Value>::InnerProductType;
using Normalized = typename Hilbert<Value>::NormalizedType;
using Series = PoissonSeries<Value, degree_, Evaluator>;

// This code follows [Kud07], section 2. Our indices start at 0, unlike those
@@ -120,14 +125,14 @@ IncrementalProjection(Function const& function,
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
}

UnboundedLowerTriangularMatrix<Inverse<Value>> α(basis_size, uninitialized);
UnboundedLowerTriangularMatrix<Inverse<Norm>> α(basis_size, uninitialized);

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

auto const F₀ = dot(function, basis[0], weight);
auto const Q₀₀ = dot(basis[0], basis[0], weight);
Norm² const F₀ = dot(function, basis[0], weight);
Norm² const Q₀₀ = dot(basis[0], basis[0], weight);
α[0][0] = 1 / Sqrt(Q₀₀);
A[0] = F₀ / Q₀₀;

@@ -138,26 +143,26 @@ IncrementalProjection(Function const& function,
for (;;) {
for (int m = m_begin; m < basis_size; ++m) {
// Contains Fₘ.
auto const F = dot(f, basis[m], weight);
Norm² const F = dot(f, basis[m], weight);

// This vector contains Qₘⱼ.
UnboundedVector<Square<Value>> Q(m + 1, uninitialized);
UnboundedVector<Norm²> Q(m + 1, uninitialized);
for (int j = 0; j <= m; ++j) {
Q[j] = dot(basis[m], basis[j], weight);
}

// This vector contains Bⱼ⁽ᵐ⁾.
UnboundedVector<Value> B(m, uninitialized);
UnboundedVector<Norm> B(m, uninitialized);
for (int j = 0; j < m; ++j) {
Value Σ_αⱼₛ_Qₘₛ{};
Norm Σ_αⱼₛ_Qₘₛ{};
for (int s = 0; s <= j; ++s) {
Σ_αⱼₛ_Qₘₛ += α[j][s] * Q[s];
}
B[j] = -Σ_αⱼₛ_Qₘₛ;
}

{
Square<Value> Σ_Bₛ⁽ᵐ⁾²{};
Norm² Σ_Bₛ⁽ᵐ⁾²{};
for (int s = 0; s < m; ++s) {
Σ_Bₛ⁽ᵐ⁾² += B[s] * B[s];
}
@@ -207,7 +212,8 @@ IncrementalProjection(Function const& function,
}

{
PoissonSeries<double, degree_, Evaluator> Σ_αₘᵢ_eᵢ = α[m][0] * basis[0];
PoissonSeries<Normalized, degree_, Evaluator> Σ_αₘᵢ_eᵢ =
α[m][0] * basis[0];
for (int i = 1; i <= m; ++i) {
Σ_αₘᵢ_eᵢ += α[m][i] * basis[i];
}
109 changes: 98 additions & 11 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -6,6 +6,9 @@
#include <random>
#include <vector>

#include "geometry/frame.hpp"
#include "geometry/grassmann.hpp"
#include "geometry/hilbert.hpp"
#include "geometry/named_quantities.hpp"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
@@ -25,13 +28,22 @@ namespace principia {
namespace numerics {
namespace frequency_analysis {

using geometry::Handedness;
using geometry::Hilbert;
using geometry::Inertial;
using geometry::Instant;
using geometry::Frame;
using geometry::Vector;
using quantities::Abs;
using quantities::Acceleration;
using quantities::AngularFrequency;
using quantities::Jerk;
using quantities::Length;
using quantities::Pow;
using quantities::Product;
using quantities::Sin;
using quantities::Snap;
using quantities::Speed;
using quantities::Square;
using quantities::Time;
using quantities::si::Metre;
@@ -51,8 +63,8 @@ class DotImplementation {
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>>
typename Hilbert<std::invoke_result_t<LFunction, Instant>,
std::invoke_result_t<RFunction, Instant>>::InnerProductType
operator()(LFunction const& left,
RFunction const& right,
Weight const& weight) const;
@@ -68,8 +80,8 @@ DotImplementation::DotImplementation(Instant const& 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>>
typename Hilbert<std::invoke_result_t<LFunction, Instant>,
std::invoke_result_t<RFunction, Instant>>::InnerProductType
DotImplementation::operator()(LFunction const& left,
RFunction const& right,
Weight const& weight) const {
@@ -78,6 +90,11 @@ DotImplementation::operator()(LFunction const& left,

class FrequencyAnalysisTest : public ::testing::Test {
protected:
using World = Frame<serialization::Frame::TestTag,
Inertial,
Handedness::Right,
serialization::Frame::TEST>;

using Series0 = PoissonSeries<Length, 0, HornerEvaluator>;
using Series4 = PoissonSeries<Length, 4, HornerEvaluator>;

@@ -154,7 +171,7 @@ TEST_F(FrequencyAnalysisTest, PreciseMode) {
EXPECT_THAT(precise_mode, RelativeErrorFrom(ω, IsNear(6.4e-11_⑴)));
}

TEST_F(FrequencyAnalysisTest, PoissonSeriesProjection) {
TEST_F(FrequencyAnalysisTest, PoissonSeriesScalarProjection) {
AngularFrequency const ω = 666.543 * π * Radian / Second;
std::mt19937_64 random(42);
std::uniform_real_distribution<> amplitude_distribution(-10.0, 10.0);
@@ -169,23 +186,23 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesProjection) {
Instant const t_max = t0_ + 100 * Radian / ω;
DotImplementation const dot(t_min, t_max);

// Projection on a 4-th degree basis accurately reconstructs the function.
// Projection on a 4th 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));
}

// Projection on a 5-th degree basis is also accurate.
// Projection on a 5th 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.
// Projection on a 3rd 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) {
@@ -195,6 +212,76 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesProjection) {
}
}

TEST_F(FrequencyAnalysisTest, PoissonSeriesVectorProjection) {
AngularFrequency const ω = 666.543 * π * Radian / Second;
std::mt19937_64 random(42);
std::uniform_real_distribution<> amplitude_distribution(-10.0, 10.0);
using VectorSeries4 =
PoissonSeries<Vector<Length, World>, 4, HornerEvaluator>;

auto random_polynomial4 = [](Instant const& t0,
std::mt19937_64& random,
std::uniform_real_distribution<>& distribution) {
auto const c0x = distribution(random) * Metre;
auto const c1x = distribution(random) * Metre / Second;
auto const c2x = distribution(random) * Metre / Pow<2>(Second);
auto const c3x = distribution(random) * Metre / Pow<3>(Second);
auto const c4x = distribution(random) * Metre / Pow<4>(Second);
auto const c0y = distribution(random) * Metre;
auto const c1y = distribution(random) * Metre / Second;
auto const c2y = distribution(random) * Metre / Pow<2>(Second);
auto const c3y = distribution(random) * Metre / Pow<3>(Second);
auto const c4y = distribution(random) * Metre / Pow<4>(Second);
auto const c0z = distribution(random) * Metre;
auto const c1z = distribution(random) * Metre / Second;
auto const c2z = distribution(random) * Metre / Pow<2>(Second);
auto const c3z = distribution(random) * Metre / Pow<3>(Second);
auto const c4z = distribution(random) * Metre / Pow<4>(Second);
Vector<Length, World> const v0({c0x, c0y, c0z});
Vector<Speed, World> const v1({c1x, c1y, c1z});
Vector<Acceleration, World> const v2({c2x, c2y, c2z});
Vector<Jerk, World> const v3({c3x, c3y, c3z});
Vector<Snap, World> const v4({c4x, c4y, c4z});

return VectorSeries4::Polynomial({v0, v1, v2, v3, v4}, t0);
};

auto const sin = random_polynomial4(t0_, random, amplitude_distribution);
auto const cos = random_polynomial4(t0_, random, amplitude_distribution);
VectorSeries4 const series(
VectorSeries4::Polynomial(VectorSeries4::Polynomial::Coefficients{}, t0_),
{{ω, VectorSeries4::Polynomials{sin, cos}}});

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

// Projection on a 4th 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, 4016));
}

// Projection on a 5th 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, 5376));
}

// Projection on a 3rd 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(projection3(t0_ + i * Radian / ω),
RelativeErrorFrom(series(t0_ + i * Radian / ω),
AllOf(Gt(1.0e-10), Lt(2.7e-7))));
}
}

TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {
AngularFrequency const ω = 666.543 * π * Radian / Second;
std::mt19937_64 random(42);
@@ -228,7 +315,7 @@ TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {
Instant const t_max = piecewise_series.t_max();
DotImplementation const dot(t_min, t_max);

// Projection on a 4-th degree basis. The errors are of the order of the
// Projection on a 4th degree basis. The errors are of the order of the
// perturbation.
auto const projection4 =
Projection<4>(ω,
@@ -293,7 +380,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionNoSecular) {
}
};

// Projection on a 4-th degree basis reconstructs the function with a decent
// Projection on a 4th degree basis reconstructs the function with a decent
// accuracy.
auto const projection4 =
IncrementalProjection<4>(series.value(),
@@ -359,7 +446,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionSecular) {
}
};

// Projection on a 4-th degree basis reconstructs the function with a decent
// Projection on a 4th degree basis reconstructs the function with a decent
// accuracy.
auto const projection4 =
IncrementalProjection<4>(series,
14 changes: 7 additions & 7 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -220,7 +220,7 @@ std::ostream& operator<<(
template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
Product<LValue, RValue>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
@@ -310,15 +310,15 @@ class PiecewisePoissonSeries {
PoissonSeries<R, r, E> const& right);
template<typename L, typename R, int l, int r, int w,
template<typename, typename, int> class E>
Product<L, R>
typename Hilbert<L, R>::InnerProductType
friend Dot(PoissonSeries<L, l, E> const& left,
PiecewisePoissonSeries<R, r, E> const& right,
PoissonSeries<double, w, E> const& weight,
Instant const& t_min,
Instant const& t_max);
template<typename L, typename R, int l, int r, int w,
template<typename, typename, int> class E>
Product<L, R>
typename Hilbert<L, R>::InnerProductType
friend Dot(PiecewisePoissonSeries<L, l, E> const& left,
PoissonSeries<R, r, E> const& right,
PoissonSeries<double, w, E> const& weight,
@@ -403,15 +403,15 @@ operator*(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
Product<LValue, RValue>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight);

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
Product<LValue, RValue>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PiecewisePoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
@@ -421,15 +421,15 @@ Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
Product<LValue, RValue>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight);

template<typename LValue, typename RValue,
int ldegree_, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
Product<LValue, RValue>
typename Hilbert<LValue, RValue>::InnerProductType
Dot(PiecewisePoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Loading