From 1c74b342e3095dfd251570188bfae8d622cc11b1 Mon Sep 17 00:00:00 2001 From: David Legg Date: Mon, 5 Feb 2024 17:54:22 -0800 Subject: [PATCH] Improve performance and stability for Polynomial root-finding Improve performance of polynomial root-finding by solving linear polynomials directly, instead of handing off to the Laguerre solver. For nonlinear polynomials, improve the stability of root-finding for polynomials with very large or very small coefficients by first conditioning the problem, normalizing the constant coefficient to 1. --- .../modeling/polynomial/Polynomial.java | 27 ++++++- .../modeling/polynomial/ComparisonsTest.java | 74 +++++++++++++++++++ 2 files changed, 97 insertions(+), 4 deletions(-) diff --git a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/Polynomial.java b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/Polynomial.java index ea8db7c3ff..bc7be13927 100644 --- a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/Polynomial.java +++ b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/Polynomial.java @@ -256,19 +256,38 @@ public Expiring max(Polynomial other) { * Finds all roots of this function in the future */ private Stream findFutureRoots() { + // TODO: In some sense, isn't having an infinite coefficient the same as a vertical line, + // hence the same as having a root at x = 0? + // Unless the value itself is non-finite, that is... // If this polynomial can never have a root, fail immediately if (this.isNonFinite() || this.isConstant()) { return Stream.empty(); } - // Defining epsilon keeps the Laguerre solver fast and stable for poorly-behaved polynomials. - final double epsilon = 2 * Arrays.stream(coefficients).map(Math::ulp).max().orElseThrow(); + if (coefficients[0] == 0.0) { + return Stream.of(ZERO); + } + + // If the polynomial is linear, solve it analytically for performance + if (this.degree() <= 1) { + double t = -getCoefficient(0) / getCoefficient(1); + if (t >= -ABSOLUTE_ACCURACY_FOR_DURATIONS / 2 && t <= MAX_SECONDS_FOR_DURATION) { + return Stream.of(Duration.roundNearest(t, SECOND)); + } else { + return Stream.empty(); + } + } + + // Condition the problem by dividing through by the first coefficient: + double[] conditionedCoefficients = Arrays.stream(coefficients).map(c -> c / coefficients[0]).toArray(); + // Defining epsilon keeps the Laguerre solver faster and more stable for poorly-behaved polynomials. + final double epsilon = 2 * Arrays.stream(conditionedCoefficients).map(Math::ulp).max().orElseThrow(); final Complex[] solutions = new LaguerreSolver(0, ABSOLUTE_ACCURACY_FOR_DURATIONS, epsilon) - .solveAllComplex(coefficients, 0); + .solveAllComplex(conditionedCoefficients, 0); return Arrays.stream(solutions) .filter(solution -> Math.abs(solution.getImaginary()) < epsilon) .map(Complex::getReal) - .filter(t -> t >= 0 && t <= MAX_SECONDS_FOR_DURATION) + .filter(t -> t >= -ABSOLUTE_ACCURACY_FOR_DURATIONS / 2 && t <= MAX_SECONDS_FOR_DURATION) .sorted() .map(t -> Duration.roundNearest(t, SECOND)); } diff --git a/contrib/src/test/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/ComparisonsTest.java b/contrib/src/test/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/ComparisonsTest.java index f707482e72..a7cad57998 100644 --- a/contrib/src/test/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/ComparisonsTest.java +++ b/contrib/src/test/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/ComparisonsTest.java @@ -299,6 +299,80 @@ void comparing_converging_nonlinear_terms_with_fine_precision() { check_extrema(false, true); } + // Unrepresentable convergence: + // These tests reflect polynomials that in theory converge, but do so in timespans + // that are too large to represent. Thus, they should be treated as non-converging. + + @Test + void comparing_linear_terms_with_convergence_unrepresentable_by_double() { + setup(() -> { + set(p, polynomial(Double.MAX_VALUE)); + set(q, polynomial(0, 0.1)); + }); + + check_comparison(p_lt_q, false, false); + check_comparison(p_lte_q, false, false); + check_comparison(p_gt_q, true, false); + check_comparison(p_gte_q, true, false); + check_extrema(true, false); + } + + @Test + void comparing_linear_terms_with_convergence_unrepresentable_by_duration() { + setup(() -> { + set(p, polynomial(Duration.MAX_VALUE.ratioOver(SECOND))); + set(q, polynomial(0, 0.1)); + }); + + check_comparison(p_lt_q, false, false); + check_comparison(p_lte_q, false, false); + check_comparison(p_gt_q, true, false); + check_comparison(p_gte_q, true, false); + check_extrema(true, false); + } + + @Test + void comparing_nonlinear_terms_with_convergence_unrepresentable_by_double() { + setup(() -> { + set(p, polynomial(Double.MAX_VALUE)); + set(q, polynomial(0, 0, 0.1)); + }); + + check_comparison(p_lt_q, false, false); + check_comparison(p_lte_q, false, false); + check_comparison(p_gt_q, true, false); + check_comparison(p_gte_q, true, false); + check_extrema(true, false); + } + + @Test + void comparing_nonlinear_terms_with_convergence_unrepresentable_by_duration() { + setup(() -> { + set(p, polynomial(Duration.MAX_VALUE.ratioOver(SECOND) * Duration.MAX_VALUE.ratioOver(SECOND))); + set(q, polynomial(0, 0, 0.1)); + }); + + check_comparison(p_lt_q, false, false); + check_comparison(p_lte_q, false, false); + check_comparison(p_gt_q, true, false); + check_comparison(p_gte_q, true, false); + check_extrema(true, false); + } + + @Test + void comparing_pathological_nonlinear_terms_with_convergence_unrepresentable_by_duration() { + setup(() -> { + set(p, polynomial(Duration.MAX_VALUE.ratioOver(SECOND) * Duration.MAX_VALUE.ratioOver(SECOND))); + set(q, polynomial(0, Duration.MIN_VALUE.ratioOver(SECOND), 1.0 + Math.ulp(1.0))); + }); + + check_comparison(p_lt_q, false, false); + check_comparison(p_lte_q, false, false); + check_comparison(p_gt_q, true, false); + check_comparison(p_gte_q, true, false); + check_extrema(true, false); + } + private void check_comparison(Resource> result, boolean expectedValue, boolean expectCrossover) { reset(); var resultDynamics = result.getDynamics().getOrThrow();