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

Commits on Aug 18, 2020

  1. Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    8d1d596 View commit details
  2. Stray / and tolerances.

    pleroy committed Aug 18, 2020
    Copy the full SHA
    34e915f View commit details
  3. Copy the full SHA
    e93288e View commit details
  4. Merge pull request #2679 from pleroy/IllConditioned

    An attempt at addressing ill-conditioned basis elements.
    pleroy authored Aug 18, 2020
    Copy the full SHA
    62b9d20 View commit details
Showing with 14 additions and 13 deletions.
  1. +12 −11 numerics/frequency_analysis_body.hpp
  2. +2 −2 numerics/frequency_analysis_test.cpp
23 changes: 12 additions & 11 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -249,17 +249,20 @@ IncrementalProjection(Function const& function,
for (int s = 0; s < m; ++s) {
Σ_Bₛ⁽ᵐ⁾² += B[s] * B[s];
}
if (Q[m] <= Σ_Bₛ⁽ᵐ⁾²) {
if (Q[m] <= Σ_Bₛ⁽ᵐ⁾² ||
(Q[m] - Σ_Bₛ⁽ᵐ⁾²) / std::max(Q[m], Σ_Bₛ⁽ᵐ⁾²) < 0x1.0p-24) {
// We arrive here when the norm of Σₛ Bₛ⁽ᵐ⁾bₛ + eₘ is small (see
// [SN97] for the notation) and, due to rounding errors, the computed
// value of that norm ends up negative or zero. It makes no sense to
// have complex numbers (or infinities) here because our function is
// real and bounded. Geometrically, we are in a situation where eₘ
// is very close to the space spanned by the (bₛ), that is, by the
// (eₛ) for i < m. The fact that the basis elements and no longer
// independent when the degree increases is duely noted by [CV84].
// Given that eₘ effectively doesn't have benefit for the projection,
// we just drop it and continue with the algorithm.
// value of the square of that norm ends up negative, zero, or very
// small. It makes no sense to have complex numbers (or infinities)
// here because our function is real and bounded. But even if the
// norm could be computed but was very small, we would end up with an
// ill-conditioned solution. Geometrically, we are in a situation
// where eₘ is very close to the space spanned by the (bₛ), that is,
// by the (eₛ) for i < m. The fact that the basis elements and no
// longer independent when the degree increases is duely noted by
// [CV84]. Given that eₘ effectively doesn't have benefit for the
// projection, we just drop it and continue with the algorithm.
LOG(ERROR) << "Q[m]: " << Q[m] << " Σ_Bₛ⁽ᵐ⁾² " << Σ_Bₛ⁽ᵐ⁾²
<< " difference: " << Q[m] - Σ_Bₛ⁽ᵐ⁾²;
LOG(ERROR) << "Dropping " << basis[m];
@@ -273,8 +276,6 @@ IncrementalProjection(Function const& function,
--m;
continue;
} else {
// TODO(phl): If we have a huge cancellation here, we should probably
// drop basis[m] too.
α[m][m] = 1 / Sqrt(Q[m] - Σ_Bₛ⁽ᵐ⁾²);
}
}
4 changes: 2 additions & 2 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -285,7 +285,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjection) {
? AllOf(Gt(6.7e-2 * Metre), Lt(7.9 * Metre))
: ω_index == 2
? AllOf(Gt(1.1e-4 * Metre), Lt(9.7e-1 * Metre))
: AllOf(Gt(1.7e-9 * Metre), Lt(1.9e-4 * Metre)))
: AllOf(Gt(4.2e-10 * Metre), Lt(1.7e-5 * Metre)))
<< ω_index;
}
if (ω_index == ωs.size()) {
@@ -306,7 +306,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjection) {
EXPECT_THAT(
projection4(t_min + i * (t_max - t_min) / 100),
RelativeErrorFrom(series.value()(t_min + i * (t_max - t_min) / 100),
AllOf(Gt(2.4e-9), Lt(3.7e-3))));
AllOf(Gt(1.3e-10), Lt(5.4e-4))));
}
}