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

Commits on Aug 19, 2020

  1. Generate the secular basis.

    pleroy committed Aug 19, 2020
    Copy the full SHA
    dd6fcc5 View commit details
  2. A test of the secular case.

    pleroy committed Aug 19, 2020

    Unverified

    This commit is not signed, but one or more authors requires that any commit attributed to them is signed.
    Copy the full SHA
    66661f8 View commit details
  3. Merge pull request #2682 from pleroy/Secular

    Support for secular terms in the Кудрявцев projection
    pleroy authored Aug 19, 2020
    Copy the full SHA
    73de885 View commit details
Showing with 105 additions and 11 deletions.
  1. +38 −8 numerics/frequency_analysis_body.hpp
  2. +67 −3 numerics/frequency_analysis_test.cpp
46 changes: 38 additions & 8 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -28,6 +28,8 @@ namespace si = quantities::si;
// A helper struct for generating the Poisson series tⁿ sin ω t and tⁿ cos ω t.
template<typename Series, int n>
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.
@@ -46,13 +48,21 @@ struct BasisGenerator;

template<typename Series, std::size_t... indices>
struct BasisGenerator<Series, std::index_sequence<indices...>> {
// TODO(phl): Will need to properly handle ω = 0.
// Basis of aperiodic terms.
static std::array<Series, Series::degree + 1> Basis(Instant const& origin);

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


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

template<typename Series, int n>
Series SeriesGenerator<Series, n>::Sin(AngularFrequency const& ω,
Instant const& origin) {
@@ -84,6 +94,13 @@ typename Series::Polynomial SeriesGenerator<Series, n>::Unit(
return typename Series::Polynomial(coefficients, origin);
}

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

template<typename Series, std::size_t... indices>
std::array<Series, 2 * Series::degree + 2>
BasisGenerator<Series, std::index_sequence<indices...>>::Basis(
@@ -203,10 +220,16 @@ IncrementalProjection(Function const& function,

std::vector<Series> basis;

// TODO(phl): Add support for secular terms.
auto const ω_basis = BasisGenerator<Series>::Basis(ω.value(), t0);
int basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
int basis_size;
if (ω.value() == AngularFrequency{}) {
auto const ω_basis = BasisGenerator<Series>::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>::Basis(ω.value(), t0);
basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
}

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

@@ -313,9 +336,16 @@ IncrementalProjection(Function const& function,
return result;
}

auto const ω_basis = BasisGenerator<Series>::Basis(ω.value(), t0);
constexpr int ω_basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
int ω_basis_size;
if (ω.value() == AngularFrequency{}) {
auto const ω_basis = BasisGenerator<Series>::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>::Basis(ω.value(), t0);
ω_basis_size = std::tuple_size_v<decltype(ω_basis)>;
std::move(ω_basis.begin(), ω_basis.end(), std::back_inserter(basis));
}
α.Extend(ω_basis_size, uninitialized);
A.Extend(ω_basis_size, uninitialized);
m_begin = basis_size;
70 changes: 67 additions & 3 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -195,7 +195,6 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesProjection) {
}
}

#if 0
TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {
AngularFrequency const ω = 666.543 * π * Radian / Second;
std::mt19937_64 random(42);
@@ -242,9 +241,8 @@ TEST_F(FrequencyAnalysisTest, PiecewisePoissonSeriesProjection) {
AllOf(Gt(2.1e-7), Lt(8.8e-4))));
}
}
#endif

TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjection) {
TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionNoSecular) {
std::mt19937_64 random(42);
std::uniform_real_distribution<> frequency_distribution(2000.0, 3000.0);

@@ -310,6 +308,72 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjection) {
}
}

TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionSecular) {
std::mt19937_64 random(42);
std::uniform_real_distribution<> frequency_distribution(2000.0, 3000.0);
std::uniform_real_distribution<> secular_distribution(-30.0, 30.0);

std::vector<AngularFrequency> ωs = {AngularFrequency{}};
Series4 series(
random_polynomial4_(t0_, random, secular_distribution), {});
for (int i = 3; i >= 1; --i) {
std::uniform_real_distribution<> amplitude_distribution(-(1 << i),
(1 << i));
ωs.push_back(frequency_distribution(random) * Radian / Second);
auto const sin = random_polynomial4_(t0_, random, amplitude_distribution);
auto const cos = random_polynomial4_(t0_, random, amplitude_distribution);
series += Series4(
Series4::Polynomial(Series4::Polynomial::Coefficients{}, t0_),
{{ωs.back(), Series4::Polynomials{sin, cos}}});
}

Instant const t_min = t0_;
Instant const t_max =
t0_ + 200 * Radian / *std::max_element(ωs.cbegin(), ωs.cend());
DotImplementation const dot(t_min, t_max);

// A perfect calculator for the frequencies of the series.
int ω_index = 0;
auto angular_frequency_calculator =
[&series, t_min, t_max, &ω_index, &ωs](
auto const& residual) -> std::optional<AngularFrequency> {
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(
Abs(residual(t_min + i * (t_max - t_min) / 100)),
ω_index == 0
? AllOf(Gt(12.4 * Metre), Lt(19.5 * Metre))
: ω_index == 1
? AllOf(Gt(8.4e-3 * Metre), Lt(3.7 * Metre))
: ω_index == 2
? AllOf(Gt(3.3e-2 * Metre), Lt(3.6 * Metre))
: ω_index == 3
? AllOf(Gt(7.5e-3 * Metre), Lt(5.4 * Metre))
: AllOf(Gt(2.2e-13 * Metre),
Lt(7.4e-10 * Metre)))
<< ω_index;
}
if (ω_index == ωs.size()) {
return std::nullopt;
} else {
return ωs[ω_index++];
}
};

// Projection on a 4-th degree basis reconstructs the function with a decent
// accuracy.
auto const projection4 =
IncrementalProjection<4>(series,
angular_frequency_calculator,
apodization::Hann<HornerEvaluator>(t_min, t_max),
dot);
for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(
projection4(t_min + i * (t_max - t_min) / 100),
RelativeErrorFrom(series(t_min + i * (t_max - t_min) / 100),
AllOf(Gt(1.2e-14), Lt(4.0e-11))));
}
}

} // namespace frequency_analysis
} // namespace numerics
} // namespace principia