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

Commits on Aug 30, 2020

  1. Copy the full SHA
    70dd4ad View commit details
  2. Incomprehensible message trying to use PointwiseInnerProduct in Dot.

    # Conflicts:
    #	numerics/frequency_analysis_test.cpp
    pleroy committed Aug 30, 2020
    Copy the full SHA
    098a059 View commit details
  3. A test for is_quantity_v.

    pleroy committed Aug 30, 2020
    Copy the full SHA
    04a6869 View commit details
  4. Copy the full SHA
    2d7ca62 View commit details
  5. Fix merge.

    pleroy committed Aug 30, 2020
    Copy the full SHA
    39065a7 View commit details
  6. Comments.

    pleroy committed Aug 30, 2020
    Copy the full SHA
    7858023 View commit details
  7. Remove debugging declaration.

    pleroy committed Aug 30, 2020
    Copy the full SHA
    5ee0d90 View commit details
  8. Merge pull request #2702 from pleroy/PSPWIP

    Add pointwise inner product for Poisson series and use it in Dot
    pleroy authored Aug 30, 2020
    Copy the full SHA
    5aa970a View commit details
22 changes: 18 additions & 4 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -7,6 +7,7 @@
#include <string>
#include <vector>

#include "geometry/hilbert.hpp"
#include "geometry/interval.hpp"
#include "geometry/named_quantities.hpp"
#include "numerics/polynomial.hpp"
@@ -36,6 +37,7 @@ FORWARD_DECLARE_FUNCTION_FROM(
namespace numerics {
namespace internal_poisson_series {

using geometry::Hilbert;
using geometry::Instant;
using geometry::Interval;
using quantities::AngularFrequency;
@@ -121,10 +123,11 @@ class PoissonSeries {
Scalar const& right);
template<typename L, typename R,
int l, int r,
template<typename, typename, int> class E>
PoissonSeries<Product<L, R>, l + r, E>
friend operator*(PoissonSeries<L, l, E> const& left,
PoissonSeries<R, r, E> const& right);
template<typename, typename, int> class E,
typename P>
auto friend Multiply(PoissonSeries<L, l, E> const& left,
PoissonSeries<R, r, E> const& right,
P const& product);
template<typename V, int d, template<typename, typename, int> class E>
friend std::ostream& operator<<(std::ostream& out,
PoissonSeries<V, d, E> const& series);
@@ -190,6 +193,17 @@ PoissonSeries<Product<LValue, RValue>, ldegree_ + rdegree_, Evaluator>
operator*(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right);

// Returns a scalar-valued Poisson series obtained by pointwise inner product of
// two vector-valued series.
template<typename LValue, typename RValue,
int ldegree_, int rdegree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<typename Hilbert<LValue, RValue>::InnerProductType,
ldegree_ + rdegree_,
Evaluator>
PointwiseInnerProduct(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right);

// Output.

template<typename Value, int degree_,
164 changes: 106 additions & 58 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -59,6 +59,80 @@ AngularFrequencyPrimitive(
}
}

// A helper for multiplication of Poisson series and pointwise inner product.
// The functor Product must take a pair of Poisson series with the types of left
// and right and return a suitable Poisson series.
template<typename LValue, typename RValue,
int ldegree_, int rdegree_,
template<typename, typename, int> class Evaluator,
typename Product>
auto Multiply(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
Product const& product) {
using Result =
PoissonSeries<typename Hilbert<LValue, RValue>::InnerProductType,
ldegree_ + rdegree_,
Evaluator>;

// Compute all the individual terms using elementary trigonometric identities
// and put them in a multimap, because the same frequency may appear multiple
// times.
std::multimap<AngularFrequency, typename Result::Polynomials> terms;
auto const aperiodic = product(left.aperiodic_, right.aperiodic_);
for (auto const& [ω, polynomials] : left.periodic_) {
terms.emplace(ω,
typename Result::Polynomials{
/*sin=*/product(polynomials.sin, right.aperiodic_),
/*cos=*/product(polynomials.cos, right.aperiodic_)});
}
for (auto const& [ω, polynomials] : right.periodic_) {
terms.emplace(ω,
typename Result::Polynomials{
/*sin=*/product(left.aperiodic_, polynomials.sin),
/*cos=*/product(left.aperiodic_, polynomials.cos)});
}
for (auto const& [ωl, polynomials_left] : left.periodic_) {
for (auto const& [ωr, polynomials_right] : right.periodic_) {
terms.emplace(
ωl - ωr,
typename Result::Polynomials{
/*sin=*/(product(-polynomials_left.cos, polynomials_right.sin) +
product(polynomials_left.sin, polynomials_right.cos)) /
2,
/*cos=*/(product(polynomials_left.sin, polynomials_right.sin) +
product(polynomials_left.cos, polynomials_right.cos)) /
2});
terms.emplace(
ωl + ωr,
typename Result::Polynomials{
/*sin=*/(product(polynomials_left.cos, polynomials_right.sin) +
product(polynomials_left.sin, polynomials_right.cos)) /
2,
/*cos=*/(product(-polynomials_left.sin, polynomials_right.sin) +
product(polynomials_left.cos, polynomials_right.cos)) /
2});
}
}

// Now group the terms together by frequency.
typename Result::PolynomialsByAngularFrequency periodic;
std::optional<AngularFrequency> previous_ω;
for (auto it = terms.cbegin(); it != terms.cend(); ++it) {
auto const& ω = it->first;
auto const& polynomials = it->second;
if (previous_ω.has_value() && previous_ω.value() == ω) {
auto& previous_polynomials = periodic.rbegin()->second;
previous_polynomials.sin += polynomials.sin;
previous_polynomials.cos += polynomials.cos;
} else {
periodic.insert(*it);
}
previous_ω = ω;
}

return Result{aperiodic, std::move(periodic)};
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, degree_, Evaluator>::PoissonSeries(
@@ -340,64 +414,36 @@ template<typename LValue, typename RValue,
PoissonSeries<Product<LValue, RValue>, ldegree_ + rdegree_, Evaluator>
operator*(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right) {
using Result =
PoissonSeries<Product<LValue, RValue>, ldegree_ + rdegree_, Evaluator>;
auto product =
[](typename PoissonSeries<LValue, ldegree_, Evaluator>::Polynomial const&
left,
typename PoissonSeries<RValue, rdegree_, Evaluator>::Polynomial const&
right) {
return left * right;
};

// Compute all the individual terms using elementary trigonometric identities
// and put them in a multimap, because the same frequency may appear multiple
// times.
std::multimap<AngularFrequency, typename Result::Polynomials> terms;
auto const aperiodic = left.aperiodic_ * right.aperiodic_;
for (auto const& [ω, polynomials] : left.periodic_) {
terms.emplace(ω,
typename Result::Polynomials{
/*sin=*/polynomials.sin * right.aperiodic_,
/*cos=*/polynomials.cos * right.aperiodic_});
}
for (auto const& [ω, polynomials] : right.periodic_) {
terms.emplace(ω,
typename Result::Polynomials{
/*sin=*/left.aperiodic_ * polynomials.sin,
/*cos=*/left.aperiodic_ * polynomials.cos});
}
for (auto const& [ωl, polynomials_left] : left.periodic_) {
for (auto const& [ωr, polynomials_right] : right.periodic_) {
terms.emplace(ωl - ωr,
typename Result::Polynomials{
/*sin=*/(-polynomials_left.cos * polynomials_right.sin +
polynomials_left.sin * polynomials_right.cos) /
2,
/*cos=*/(polynomials_left.sin * polynomials_right.sin +
polynomials_left.cos * polynomials_right.cos) /
2});
terms.emplace(ωl + ωr,
typename Result::Polynomials{
/*sin=*/(polynomials_left.cos * polynomials_right.sin +
polynomials_left.sin * polynomials_right.cos) /
2,
/*cos=*/(-polynomials_left.sin * polynomials_right.sin +
polynomials_left.cos * polynomials_right.cos) /
2});
}
}
return Multiply<LValue, RValue, ldegree_, rdegree_, Evaluator>(
left, right, product);
}

// Now group the terms together by frequency.
typename Result::PolynomialsByAngularFrequency periodic;
std::optional<AngularFrequency> previous_ω;
for (auto it = terms.cbegin(); it != terms.cend(); ++it) {
auto const& ω = it->first;
auto const& polynomials = it->second;
if (previous_ω.has_value() && previous_ω.value() == ω) {
auto& previous_polynomials = periodic.rbegin()->second;
previous_polynomials.sin += polynomials.sin;
previous_polynomials.cos += polynomials.cos;
} else {
periodic.insert(*it);
}
previous_ω = ω;
}
template<typename LValue, typename RValue,
int ldegree_, int rdegree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<typename Hilbert<LValue, RValue>::InnerProductType,
ldegree_ + rdegree_,
Evaluator>
PointwiseInnerProduct(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right) {
auto product =
[](typename PoissonSeries<LValue, ldegree_, Evaluator>::Polynomial const&
left,
typename PoissonSeries<RValue, rdegree_, Evaluator>::Polynomial const&
right) {
return PointwiseInnerProduct(left, right);
};

return {aperiodic, std::move(periodic)};
return Multiply<LValue, RValue, ldegree_, rdegree_, Evaluator>(
left, right, product);
}

template<typename Value, int degree_,
@@ -440,7 +486,7 @@ Dot(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
Instant const& t_min,
Instant const& t_max) {
auto const integrand = left * right * weight;
auto const integrand = PointwiseInnerProduct(left, right) * weight;
auto const primitive = integrand.Primitive();
return (primitive(t_max) - primitive(t_min)) / (t_max - t_min);
}
@@ -692,7 +738,8 @@ Product<LValue, RValue> Dot(
using Result = Primitive<Product<LValue, RValue>, Time>;
Result result;
for (int i = 0; i < right.series_.size(); ++i) {
auto const integrand = left * right.series_[i] * weight;
auto const integrand =
PointwiseInnerProduct(left, right.series_[i]) * weight;
auto const primitive = integrand.Primitive();
result += primitive(right.bounds_[i + 1]) - primitive(right.bounds_[i]);
}
@@ -721,7 +768,8 @@ Product<LValue, RValue> Dot(
using Result = Primitive<Product<LValue, RValue>, Time>;
Result result;
for (int i = 0; i < left.series_.size(); ++i) {
auto const integrand = left.series_[i] * right * weight;
auto const integrand =
PointwiseInnerProduct(left.series_[i], right) * weight;
auto const primitive = integrand.Primitive();
result += primitive(left.bounds_[i + 1]) - primitive(left.bounds_[i]);
}
45 changes: 45 additions & 0 deletions numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -4,6 +4,8 @@
#include <limits>
#include <memory>

#include "geometry/frame.hpp"
#include "geometry/grassmann.hpp"
#include "geometry/named_quantities.hpp"
#include "gtest/gtest.h"
#include "numerics/apodization.hpp"
@@ -17,19 +19,32 @@
namespace principia {
namespace numerics {

using geometry::Displacement;
using geometry::Frame;
using geometry::Handedness;
using geometry::Inertial;
using geometry::Instant;
using geometry::Vector;
using geometry::Velocity;
using quantities::Acceleration;
using quantities::AngularFrequency;
using quantities::Cos;
using quantities::Sin;
using quantities::Sqrt;
using quantities::Time;
using quantities::si::Metre;
using quantities::si::Radian;
using quantities::si::Second;
using testing_utilities::AlmostEquals;
using testing_utilities::VanishesBefore;

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

using Degree1 = PoissonSeries<double, 1, HornerEvaluator>;

PoissonSeriesTest()
@@ -138,6 +153,36 @@ TEST_F(PoissonSeriesTest, Algebra) {
(*pb_)(t0_ + 1 * Second), 6, 38));
}

TEST_F(PoissonSeriesTest, PointwiseInnerProduct) {
using Degree2 = PoissonSeries<Displacement<World>, 2, HornerEvaluator>;
Degree2::Polynomial::Coefficients const coefficients_a({
Displacement<World>({0 * Metre,
0 * Metre,
1 * Metre}),
Velocity<World>({0 * Metre / Second,
1 * Metre / Second,
0 * Metre / Second}),
Vector<Acceleration, World>({1 * Metre / Second / Second,
0 * Metre / Second / Second,
0 * Metre / Second / Second})});
Degree2::Polynomial::Coefficients const coefficients_b({
Displacement<World>({0 * Metre,
2 * Metre,
3 * Metre}),
Velocity<World>({-1 * Metre / Second,
1 * Metre / Second,
0 * Metre / Second}),
Vector<Acceleration, World>({1 * Metre / Second / Second,
1 * Metre / Second / Second,
-2 * Metre / Second / Second})});
Degree2 const pa(Degree2::Polynomial({coefficients_a}, t0_), {{}});
Degree2 const pb(Degree2::Polynomial({coefficients_b}, t0_), {{}});

auto const product = PointwiseInnerProduct(pa, pb);
EXPECT_THAT(product(t0_ + 1 * Second),
AlmostEquals(5 * Metre * Metre, 0));
}

TEST_F(PoissonSeriesTest, Primitive) {
auto const actual_primitive = pb_->Primitive();

1 change: 1 addition & 0 deletions quantities/quantities.vcxproj
Original file line number Diff line number Diff line change
@@ -39,6 +39,7 @@
<ClCompile Include="elementary_functions_test.cpp" />
<ClCompile Include="parser_test.cpp" />
<ClCompile Include="quantities_test.cpp" />
<ClCompile Include="traits_test.cpp" />
<ClCompile Include="wide_test.cpp" />
</ItemGroup>
</Project>
3 changes: 3 additions & 0 deletions quantities/quantities.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -109,5 +109,8 @@
<ClCompile Include="..\numerics\cbrt.cpp">
<Filter>Source Files</Filter>
</ClCompile>
<ClCompile Include="traits_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
</ItemGroup>
</Project>
2 changes: 2 additions & 0 deletions quantities/traits.hpp
Original file line number Diff line number Diff line change
@@ -23,6 +23,8 @@ template<typename T>
struct is_quantity : std::is_arithmetic<T>, not_constructible {};
template<typename D>
struct is_quantity<Quantity<D>> : std::true_type, not_constructible {};
template<typename T>
struct is_quantity<T const> : is_quantity<T> {};

template<typename T>
constexpr bool is_quantity_v = is_quantity<T>::value;
22 changes: 22 additions & 0 deletions quantities/traits_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@

#include "quantities/traits.hpp"

#include "gtest/gtest.h"
#include "quantities/named_quantities.hpp"

namespace principia {
namespace quantities {

TEST(Traits, IsQuantityV) {
static_assert(is_quantity_v<int>);
static_assert(is_quantity_v<double>);
static_assert(is_quantity_v<Area>);
static_assert(is_quantity_v<Frequency const>);
// Not sure if the following is what we want, but at least let's nail it in a
// test.
static_assert(!is_quantity_v<Entropy&>);
static_assert(!is_quantity_v<float const&>);
}

} // namespace quantities
} // namespace principia