From 0616726db71262e44dd3f47b676d04451ce6f693 Mon Sep 17 00:00:00 2001 From: Pascal Leroy Date: Mon, 6 Feb 2023 12:04:01 +0800 Subject: [PATCH 1/2] A more efficient L-Infinity check. --- numerics/fit_hermite_spline_body.hpp | 22 +++++++++++---------- numerics/hermite3.hpp | 14 +++++++++++++- numerics/hermite3_body.hpp | 29 ++++++++++++++++++++++------ numerics/hermite3_test.cpp | 22 +++++++++++++++++++++ 4 files changed, 70 insertions(+), 17 deletions(-) diff --git a/numerics/fit_hermite_spline_body.hpp b/numerics/fit_hermite_spline_body.hpp index 00a8df0426..44c3659c00 100644 --- a/numerics/fit_hermite_spline_body.hpp +++ b/numerics/fit_hermite_spline_body.hpp @@ -25,14 +25,16 @@ absl::StatusOr> FitHermiteSpline( typename Hilbert>::NormType const& tolerance) { using Iterator = typename Samples::const_iterator; - auto interpolation_error = [get_argument, get_derivative, get_value]( - Iterator begin, Iterator last) { - return Hermite3( - {get_argument(*begin), get_argument(*last)}, - {get_value(*begin), get_value(*last)}, - {get_derivative(*begin), get_derivative(*last)}) - .LInfinityError(Range(begin, last + 1), get_argument, get_value); - }; + auto interpolation_error_is_within_tolerance = + [get_argument, get_derivative, get_value, tolerance]( + Iterator const begin, Iterator const last) { + return Hermite3( + {get_argument(*begin), get_argument(*last)}, + {get_value(*begin), get_value(*last)}, + {get_derivative(*begin), get_derivative(*last)}) + .LInfinityErrorIsWithin( + Range(begin, last + 1), get_argument, get_value, tolerance); + }; std::list tail; if (samples.size() < 3) { @@ -44,7 +46,7 @@ absl::StatusOr> FitHermiteSpline( Iterator begin = samples.begin(); Iterator const last = samples.end() - 1; while (last - begin + 1 >= 3 && - interpolation_error(begin, last) >= tolerance) { + !interpolation_error_is_within_tolerance(begin, last)) { // Look for a cubic that fits the beginning within |tolerance| and // such the cubic fitting one more sample would not fit the samples within // |tolerance|. @@ -68,7 +70,7 @@ absl::StatusOr> FitHermiteSpline( if (middle == lower) { break; } - if (interpolation_error(begin, middle) < tolerance) { + if (interpolation_error_is_within_tolerance(begin, middle)) { lower = middle; } else { upper = middle; diff --git a/numerics/hermite3.hpp b/numerics/hermite3.hpp index 37151e6a5e..6893c6d8f6 100644 --- a/numerics/hermite3.hpp +++ b/numerics/hermite3.hpp @@ -22,6 +22,7 @@ using geometry::Hilbert; // TODO(phl): Invert the two template arguments for consistency with Derivative. template class Hermite3 final { + using NormType = typename Hilbert>::NormType; public: using Derivative1 = Derivative; @@ -43,13 +44,24 @@ class Hermite3 final { // Returns the largest error (in the given |norm|) between this polynomial and // the given |samples|. template - typename Hilbert>::NormType LInfinityError( + NormType LInfinityError( Samples const& samples, std::function const& get_argument, std::function const& get_value) const; + // Returns true if the |LInfinityError| is less than |tolerance|. More + // efficient than the above function in the case where it returns false. + template + bool LInfinityErrorIsWithin( + Samples const& samples, + std::function const& + get_argument, + std::function const& + get_value, + NormType const& tolerance) const; + private: using Derivative2 = Derivative; using Derivative3 = Derivative; diff --git a/numerics/hermite3_body.hpp b/numerics/hermite3_body.hpp index 52f57a2202..a000660042 100644 --- a/numerics/hermite3_body.hpp +++ b/numerics/hermite3_body.hpp @@ -65,14 +65,13 @@ BoundedArray Hermite3::FindExtrema() const { template template -typename Hilbert>::NormType -Hermite3::LInfinityError( +auto Hermite3::LInfinityError( Samples const& samples, - std::function const& get_argument, + std::function const& + get_argument, std::function const& - get_value) const { - typename Hilbert>::NormType result{}; + get_value) const -> NormType { + NormType result{}; for (const auto& sample : samples) { result = std::max(result, Hilbert>::Norm( @@ -81,6 +80,24 @@ Hermite3::LInfinityError( return result; } +template +template +bool Hermite3::LInfinityErrorIsWithin( + Samples const& samples, + std::function const& + get_argument, + std::function const& + get_value, + NormType const& tolerance) const { + for (const auto& sample : samples) { + if (Hilbert>::Norm(Evaluate(get_argument(sample)) - + get_value(sample)) >= tolerance) { + return false; + } + } + return true; +} + } // namespace internal_hermite3 } // namespace numerics } // namespace principia diff --git a/numerics/hermite3_test.cpp b/numerics/hermite3_test.cpp index 761372ebb0..1e55312abf 100644 --- a/numerics/hermite3_test.cpp +++ b/numerics/hermite3_test.cpp @@ -110,6 +110,17 @@ TEST_F(Hermite3Test, OneDimensionalInterpolationError) { /*get_argument=*/[](auto&& pair) -> auto&& { return pair.first; }, /*get_value=*/[](auto&& pair) -> auto&& { return pair.second; }), Eq(1 / 16.0)); + + EXPECT_TRUE(not_a_quartic.LInfinityErrorIsWithin( + samples, + /*get_argument=*/[](auto&& pair) -> auto&& { return pair.first; }, + /*get_value=*/[](auto&& pair) -> auto&& { return pair.second; }, + /*tolerance=*/0.1)); + EXPECT_FALSE(not_a_quartic.LInfinityErrorIsWithin( + samples, + /*get_argument=*/[](auto&& pair) -> auto&& { return pair.first; }, + /*get_value=*/[](auto&& pair) -> auto&& { return pair.second; }, + /*tolerance=*/0.05)); } TEST_F(Hermite3Test, ThreeDimensionalInterpolationError) { @@ -138,6 +149,17 @@ TEST_F(Hermite3Test, ThreeDimensionalInterpolationError) { /*get_argument=*/[](auto&& pair) -> auto&& { return pair.first; }, /*get_value=*/[](auto&& pair) -> auto&& { return pair.second; }), IsNear(1.5_(1) * Centi(Metre))); + + EXPECT_TRUE(not_a_circle.LInfinityErrorIsWithin( + samples, + /*get_argument=*/[](auto&& pair) -> auto&& { return pair.first; }, + /*get_value=*/[](auto&& pair) -> auto&& { return pair.second; }, + /*tolerance=*/2 * Centi(Metre))); + EXPECT_FALSE(not_a_circle.LInfinityErrorIsWithin( + samples, + /*get_argument=*/[](auto&& pair) -> auto&& { return pair.first; }, + /*get_value=*/[](auto&& pair) -> auto&& { return pair.second; }, + /*tolerance=*/1 * Centi(Metre))); } } // namespace numerics From d34eb68fb153771580e969757a2c63a31606d4f9 Mon Sep 17 00:00:00 2001 From: Pascal Leroy Date: Tue, 7 Feb 2023 12:51:16 +0800 Subject: [PATCH 2/2] Lint and a better name. --- numerics/fit_hermite_spline_body.hpp | 10 +++++----- numerics/hermite3.hpp | 1 + 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/numerics/fit_hermite_spline_body.hpp b/numerics/fit_hermite_spline_body.hpp index 44c3659c00..ca30eb24b4 100644 --- a/numerics/fit_hermite_spline_body.hpp +++ b/numerics/fit_hermite_spline_body.hpp @@ -36,11 +36,11 @@ absl::StatusOr> FitHermiteSpline( Range(begin, last + 1), get_argument, get_value, tolerance); }; - std::list tail; + std::list fit; if (samples.size() < 3) { // With 0 or 1 points there is nothing to interpolate, with 2 we cannot // estimate the error. - return tail; + return fit; } Iterator begin = samples.begin(); @@ -76,7 +76,7 @@ absl::StatusOr> FitHermiteSpline( upper = middle; } } - tail.push_back(lower); + fit.push_back(lower); begin = lower; } @@ -85,9 +85,9 @@ absl::StatusOr> FitHermiteSpline( // point, except at the end where we give up because we don't have enough // points left. #if PRINCIPIA_MUST_ALWAYS_DOWNSAMPLE - CHECK_LT(tail.size(), samples.size() - 2); + CHECK_LT(fit.size(), samples.size() - 2); #endif - return tail; + return fit; } } // namespace internal_fit_hermite_spline diff --git a/numerics/hermite3.hpp b/numerics/hermite3.hpp index 6893c6d8f6..bc590973bc 100644 --- a/numerics/hermite3.hpp +++ b/numerics/hermite3.hpp @@ -23,6 +23,7 @@ using geometry::Hilbert; template class Hermite3 final { using NormType = typename Hilbert>::NormType; + public: using Derivative1 = Derivative;