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: 6b44f7fa5d55
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 10ff8c1bc35b
Choose a head ref
  • 10 commits
  • 11 files changed
  • 2 contributors

Commits on Jul 27, 2020

  1. Too many templates.

    pleroy committed Jul 27, 2020
    Copy the full SHA
    f5accb5 View commit details

Commits on Aug 5, 2020

  1. Copy the full SHA
    39972d4 View commit details

Commits on Aug 6, 2020

  1. Copy the full SHA
    433296e View commit details
  2. Copy the full SHA
    9ca311d View commit details

Commits on Aug 7, 2020

  1. A better test.

    pleroy committed Aug 7, 2020
    Copy the full SHA
    2c69b6a View commit details
  2. Bibliography and comments.

    pleroy committed Aug 7, 2020
    Copy the full SHA
    e5258e2 View commit details
  3. Evaluate -> ().

    pleroy committed Aug 7, 2020
    Copy the full SHA
    a3bbf9c View commit details
  4. After egg's review.

    pleroy committed Aug 7, 2020
    Copy the full SHA
    6067c17 View commit details
  5. Lint.

    pleroy committed Aug 7, 2020
    Copy the full SHA
    b39a7bd View commit details
  6. Merge pull request #2666 from pleroy/PreciseMode

    Computation of the precise mode of a quasi-periodic function using dot product
    pleroy authored Aug 7, 2020
    Copy the full SHA
    10ff8c1 View commit details
10 changes: 10 additions & 0 deletions documentation/bibliography.bib
Original file line number Diff line number Diff line change
@@ -221,6 +221,16 @@ @article{Celledoni2008
volume = {30},
}

@article{Chapront1995,
author = {Chapront, J.},
date = {1995-01},
journaltitle = {Astronomy \& Astrophysics},
keywords = {ephemerides,planet and satellites: general,methods: numerical},
pages = {181-192},
title = {Representation of planetary ephemerides by frequency analysis. Applications to the five outer planets},
volume = {109},
}

@article{Chin1997,
author = {Chin, S. A.},
url = {http://www.sciencedirect.com/science/article/pii/S0375960197000030},
Binary file modified documentation/bibliography.pdf
Binary file not shown.
80 changes: 40 additions & 40 deletions numerics/apodization_test.cpp
Original file line number Diff line number Diff line change
@@ -34,82 +34,82 @@ class ApodizationTest : public ::testing::Test {

TEST_F(ApodizationTest, Dirichlet) {
auto a = apodization::Dirichlet<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(1, 0));
EXPECT_THAT(a(t1_), AlmostEquals(1, 0));
EXPECT_THAT(a(t0_), AlmostEquals(1, 0));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(1, 0));
}

TEST_F(ApodizationTest, Sine) {
auto a = apodization::Sine<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(0, 0));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(Sqrt(3) / 2, 0));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), VanishesBefore(1, 1));
EXPECT_THAT(a(t1_), AlmostEquals(0, 0));
EXPECT_THAT(a(t0_), AlmostEquals(Sqrt(3) / 2, 0));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), VanishesBefore(1, 1));
}

TEST_F(ApodizationTest, Hann) {
auto a = apodization::Hann<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(0, 0));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.75, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), VanishesBefore(1, 0));
EXPECT_THAT(a(t1_), AlmostEquals(0, 0));
EXPECT_THAT(a(t0_), AlmostEquals(0.75, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, Hamming) {
auto a = apodization::Hamming<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(2.0 / 23.0, 0));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(71.0 / 92.0, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(2.0 / 23.0, 0));
EXPECT_THAT(a(t1_), AlmostEquals(2.0 / 23.0, 0));
EXPECT_THAT(a(t0_), AlmostEquals(71.0 / 92.0, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(2.0 / 23.0, 0));
}

TEST_F(ApodizationTest, Blackman) {
auto a = apodization::Blackman<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.63, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 1));
EXPECT_THAT(a.Evaluate(t2_), VanishesBefore(1, 0));
EXPECT_THAT(a(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a(t0_), AlmostEquals(0.63, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 1));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, ExactBlackman) {
auto a = apodization::ExactBlackman<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(8.0 / 1163.0, 37));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(11843.0 / 18608.0, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(8.0 / 1163.0, 37));
EXPECT_THAT(a(t1_), AlmostEquals(8.0 / 1163.0, 37));
EXPECT_THAT(a(t0_), AlmostEquals(11843.0 / 18608.0, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(8.0 / 1163.0, 37));
}

TEST_F(ApodizationTest, Nuttall) {
auto a = apodization::Nuttall<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.514746, 2));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), VanishesBefore(1, 0));
EXPECT_THAT(a(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a(t0_), AlmostEquals(0.514746, 2));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, BlackmanNuttall) {
auto a = apodization::BlackmanNuttall<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(0.0003628, 703));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.5292298, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(0.0003628, 703));
EXPECT_THAT(a(t1_), AlmostEquals(0.0003628, 703));
EXPECT_THAT(a(t0_), AlmostEquals(0.5292298, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(0.0003628, 703));
}

TEST_F(ApodizationTest, BlackmanHarris) {
auto a = apodization::BlackmanHarris<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(0.00006, 151));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.520575, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(0.00006, 151));
EXPECT_THAT(a(t1_), AlmostEquals(0.00006, 151));
EXPECT_THAT(a(t0_), AlmostEquals(0.520575, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(0.00006, 151));
}

TEST_F(ApodizationTest, ISO18431_2) {
auto a = apodization::ISO18431_2<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(-0.00195313 / 4.63867187, 440));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.9194336 / 4.63867187, 6));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 1));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(-0.00195313 / 4.63867187, 440));
EXPECT_THAT(a(t1_), AlmostEquals(-0.00195313 / 4.63867187, 440));
EXPECT_THAT(a(t0_), AlmostEquals(0.9194336 / 4.63867187, 6));
EXPECT_THAT(a(mid_), AlmostEquals(1, 1));
EXPECT_THAT(a(t2_), AlmostEquals(-0.00195313 / 4.63867187, 440));
}

} // namespace numerics
50 changes: 50 additions & 0 deletions numerics/frequency_analysis.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@

#pragma once

#include <type_traits>

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

namespace principia {
namespace numerics {
namespace frequency_analysis {
namespace internal_frequency_analysis {

using geometry::Instant;
using geometry::Interval;
using quantities::AngularFrequency;
using quantities::Primitive;
using quantities::Product;
using quantities::Time;

// Computes the precise mode of a quasi-periodic function, assuming that the
// mode is over the interval fft_mode (so named because it has presumably been
// obtained using FFT). The function weight is the apodization function. The
// function dot is the dot product between function and Poisson series. See
// [Cha95].
template<typename Function,
typename RValue, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
AngularFrequency PreciseMode(
Interval<AngularFrequency> const& fft_mode,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
std::function<Primitive<Product<std::invoke_result_t<Function, Instant>,
RValue>,
Time>(
Function const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight)> const& dot);

} // namespace internal_frequency_analysis

using internal_frequency_analysis::PreciseMode;

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

#include "numerics/frequency_analysis_body.hpp"
58 changes: 58 additions & 0 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@

#pragma once

#include "numerics/frequency_analysis.hpp"

#include <functional>

#include "numerics/root_finders.hpp"

namespace principia {
namespace numerics {
namespace frequency_analysis {
namespace internal_frequency_analysis {

using quantities::Square;

template<typename Function,
typename RValue, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
AngularFrequency PreciseMode(
Interval<AngularFrequency> const& fft_mode,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
std::function<Primitive<Product<std::invoke_result_t<Function, Instant>,
RValue>,
Time>(
Function const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight)> const& dot) {
using DotResult =
Primitive<Product<std::invoke_result_t<Function, Instant>, RValue>, Time>;
using Degree0 = PoissonSeries<double, 0, Evaluator>;

auto amplitude = [&dot, &function, &weight](AngularFrequency const& ω) {
Instant const& t0 = weight.origin();
Degree0 const sin(typename Degree0::Polynomial({0}, t0),
{{ω,
{/*sin=*/typename Degree0::Polynomial({1}, t0),
/*cos=*/typename Degree0::Polynomial({0}, t0)}}});
Degree0 const cos(typename Degree0::Polynomial({0}, t0),
{{ω,
{/*sin=*/typename Degree0::Polynomial({0}, t0),
/*cos=*/typename Degree0::Polynomial({1}, t0)}}});
auto const sin_amplitude = dot(function, sin, weight);
auto const cos_amplitude = dot(function, cos, weight);
return sin_amplitude * sin_amplitude + cos_amplitude * cos_amplitude;
};

return GoldenSectionSearch(amplitude,
fft_mode.min,
fft_mode.max,
std::greater<Square<DotResult>>());
}

} // namespace internal_frequency_analysis
} // namespace frequency_analysis
} // namespace numerics
} // namespace principia
103 changes: 103 additions & 0 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,103 @@

#include "numerics/frequency_analysis.hpp"

#include <random>
#include <vector>

#include "geometry/named_quantities.hpp"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "numerics/apodization.hpp"
#include "numerics/fast_fourier_transform.hpp"
#include "numerics/polynomial_evaluators.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/approximate_quantity.hpp"
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/numerics_matchers.hpp"

namespace principia {
namespace numerics {
namespace frequency_analysis {

using geometry::Instant;
using quantities::AngularFrequency;
using quantities::Length;
using quantities::Sin;
using quantities::Time;
using quantities::si::Metre;
using quantities::si::Radian;
using quantities::si::Second;
using testing_utilities::IsNear;
using testing_utilities::RelativeErrorFrom;
using testing_utilities::operator""_⑴;

class FrequencyAnalysisTest : public ::testing::Test {};

TEST_F(FrequencyAnalysisTest, PreciseMode) {
using FFT = FastFourierTransform<Length, 1 << 16>;
Instant const t0;
AngularFrequency const ω = 666.543 * π / FFT::size * Radian / Second;
Time const Δt = 1 * Second;
std::mt19937_64 random(42);
std::uniform_real_distribution<> amplitude_noise(-0.1, 0.1);
std::uniform_real_distribution<> frequency_noise(-100.0, 100.0);

using Degree0 = PoissonSeries<double, 0, HornerEvaluator>;
Degree0::PolynomialsByAngularFrequency polynomials;

// Main harmonic.
polynomials.emplace(
ω,
Degree0::Polynomials{/*sin=*/Degree0::Polynomial({1}, t0),
/*cos=*/Degree0::Polynomial({0}, t0)});

// Noise with lower amplitude and higher frequency.
for (int i = 0; i < 10; ++i) {
auto const sin_amplitude = amplitude_noise(random);
auto const cos_amplitude = amplitude_noise(random);
polynomials.emplace(
ω * frequency_noise(random),
Degree0::Polynomials{/*sin=*/Degree0::Polynomial({sin_amplitude}, t0),
/*cos=*/Degree0::Polynomial({cos_amplitude}, t0)});
}
Degree0 const sin(Degree0::Polynomial({amplitude_noise(random)}, t0),
polynomials);

Instant const t_min = t0;
Instant const t_max = t0 + (FFT::size - 1) * Δt;
std::vector<Length> signal;
for (int n = 0; n < FFT::size; ++n) {
signal.push_back(sin(t0 + n * Δt) * Metre);
}

// Won't fit on the stack.
auto transform = std::make_unique<FFT>(signal, Δt);

// The FFT gives us an accuracy which is of the order of the number of points.
auto const mode = transform->Mode();
EXPECT_THAT(mode.midpoint(),
RelativeErrorFrom(ω, IsNear(8.1e-4_⑴)));

std::function<Time(Degree0 const& left,
Degree0 const& right,
Degree0 const& weight)> const dot =
[t_min, t_max](Degree0 const& left,
Degree0 const& right,
Degree0 const& weight) {
return Dot(left, right, weight, t_min, t_max);
};

// The precise analysis is only limited by our ability to pinpoint the
// maximum.
auto const precise_mode = PreciseMode(
mode, sin, apodization::Hann<HornerEvaluator>(t_min, t_max), dot);
EXPECT_THAT(precise_mode,
RelativeErrorFrom(ω, IsNear(6.4e-11_⑴)));
}

} // namespace frequency_analysis
} // namespace numerics
} // namespace principia
3 changes: 3 additions & 0 deletions numerics/numerics.vcxproj
Original file line number Diff line number Diff line change
@@ -33,6 +33,8 @@
<ClInclude Include="fixed_arrays_body.hpp" />
<ClInclude Include="double_precision.hpp" />
<ClInclude Include="double_precision_body.hpp" />
<ClInclude Include="frequency_analysis.hpp" />
<ClInclude Include="frequency_analysis_body.hpp" />
<ClInclude Include="hermite3.hpp" />
<ClInclude Include="hermite3_body.hpp" />
<ClInclude Include="legendre.hpp" />
@@ -74,6 +76,7 @@
<ClCompile Include="finite_difference_test.cpp" />
<ClCompile Include="fit_hermite_spline_test.cpp" />
<ClCompile Include="fixed_arrays_test.cpp" />
<ClCompile Include="frequency_analysis_test.cpp" />
<ClCompile Include="hermite3_test.cpp" />
<ClCompile Include="legendre_test.cpp" />
<ClCompile Include="max_abs_normalized_associated_legendre_functions_test.cc" />
9 changes: 9 additions & 0 deletions numerics/numerics.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -143,6 +143,12 @@
<ClInclude Include="fast_fourier_transform_body.hpp">
<Filter>Source Files</Filter>
</ClInclude>
<ClInclude Include="frequency_analysis.hpp">
<Filter>Header Files</Filter>
</ClInclude>
<ClInclude Include="frequency_analysis_body.hpp">
<Filter>Source Files</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<ClCompile Include="чебышёв_series_test.cpp">
@@ -220,6 +226,9 @@
<ClCompile Include="fast_fourier_transform_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
<ClCompile Include="frequency_analysis_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<Text Include="xgscd.proto.txt">
Loading