From 70d1aff928c0d2bf11bb63517e1e95c19591553f Mon Sep 17 00:00:00 2001 From: David Legg Date: Mon, 5 Feb 2024 17:32:09 -0800 Subject: [PATCH 1/6] Make Expiry comparable --- .../java/gov/nasa/jpl/aerie/contrib/streamline/core/Expiry.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Expiry.java b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Expiry.java index f8242572a9..f7e245fa16 100644 --- a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Expiry.java +++ b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Expiry.java @@ -8,7 +8,7 @@ /** * The time at which a value expires. */ -public record Expiry(Optional value) { +public record Expiry(Optional value) implements Comparable { public static Expiry NEVER = expiry(Optional.empty()); public static Expiry at(Duration t) { From 06e840c0e2be0206e4c890d2f54a8a7c385339c0 Mon Sep 17 00:00:00 2001 From: David Legg Date: Mon, 5 Feb 2024 17:33:43 -0800 Subject: [PATCH 2/6] Account for expiry in dynamicsChange condition There's an edge case that was missed by dynamicsChange conditions: when a cell is set to a dynamics with the same data, but a later expiry. This was discovered while debugging clamped integrals whose integrands changed from one limited-out value to another. The solver would emit the same effective rate for the integral (0), but with an expiry that changed from zero time to some positive time. By only reacting to the first change, the zero expiry time was propagated forward, but the corrected non-zero time wasn't. --- .../gov/nasa/jpl/aerie/contrib/streamline/core/Resources.java | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Resources.java b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Resources.java index 3ed1821ec3..998851f285 100644 --- a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Resources.java +++ b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Resources.java @@ -99,9 +99,11 @@ public static > Condition dynamicsChange(Resource re final Duration startTime = currentTime(); Condition result = (positive, atEarliest, atLatest) -> { var currentDynamics = resource.getDynamics(); + var elapsedTime = currentTime().minus(startTime); boolean haveChanged = startingDynamics.match( start -> currentDynamics.match( - current -> !current.data().equals(start.data().step(currentTime().minus(startTime))), + current -> !current.data().equals(start.data().step(elapsedTime)) || + !current.expiry().equals(start.expiry().minus(elapsedTime)), ignored -> true), startException -> currentDynamics.match( ignored -> true, From c783a868ab57e138a9fc244e8b64e10ac05f04ad Mon Sep 17 00:00:00 2001 From: David Legg Date: Mon, 5 Feb 2024 17:39:13 -0800 Subject: [PATCH 3/6] Take expiry into account for approximation interval When solving for an interval over which to take an approximation, we take the expiry of the dynamics to be approximated into account. In the case when no solution is found, we should still take expiry into account. This is especially important when that expiry is much shorter than the typical approximation interval. In this case, we can often choose the full expiry as the approximation interval without reaching our maximum error tolerance, resulting in a NoBracketingException from the solver. In this case, we should choose the expiry, not the maximum approximation interval, to be consistent with the time range we searched. --- .../modeling/black_box/IntervalFunctions.java | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/black_box/IntervalFunctions.java b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/black_box/IntervalFunctions.java index 3024d654d9..6e91504987 100644 --- a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/black_box/IntervalFunctions.java +++ b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/black_box/IntervalFunctions.java @@ -59,23 +59,26 @@ public static > Function, Duration> byBound exp.data(), t, maximumError)); + var effectiveMinSamplePeriod = Duration.min(e, minimumSamplePeriod); + var effectiveMaxSamplePeriod = Duration.min(e, maximumSamplePeriod); + try { double intervalSize = solver.solve( 100, errorFn, - Duration.min(e, minimumSamplePeriod).ratioOver(SECOND), - Duration.min(e, maximumSamplePeriod).ratioOver(SECOND)); + effectiveMinSamplePeriod.ratioOver(SECOND), + effectiveMaxSamplePeriod.ratioOver(SECOND)); return Duration.roundNearest(intervalSize, SECOND); } catch (NoBracketingException x) { if (errorFn.value(minimumSamplePeriod.ratioOver(SECOND)) > 0) { // maximum error > estimated error, best case - return maximumSamplePeriod; + return effectiveMaxSamplePeriod; } else { // maximum error < estimated error, worst case - return minimumSamplePeriod; + return effectiveMinSamplePeriod; } } catch (TooManyEvaluationsException | NumberIsTooLargeException x) { - return minimumSamplePeriod; + return effectiveMinSamplePeriod; } }; } From 8f9ebbda9243c5b55d8f5e06c8af8cc76823af84 Mon Sep 17 00:00:00 2001 From: David Legg Date: Mon, 5 Feb 2024 17:53:56 -0800 Subject: [PATCH 4/6] Prevent overflow in Polynomial root-finding --- .../contrib/streamline/modeling/polynomial/Polynomial.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 7e00dd83a9..ea8db7c3ff 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 @@ -195,7 +195,7 @@ private Expiry findExpiryNearRoot(Predicate expires) { // Do a binary search to find the exact transition time while (end.longerThan(start.plus(EPSILON))) { - Duration midpoint = start.plus(end).dividedBy(2); + Duration midpoint = start.plus(end.minus(start).dividedBy(2)); if (expires.test(midpoint)) { end = midpoint; } else { From 762a534a51a5d5f90425c0bccf625f804acb88b2 Mon Sep 17 00:00:00 2001 From: David Legg Date: Mon, 5 Feb 2024 17:54:22 -0800 Subject: [PATCH 5/6] 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(); From fac658253faed0f3feca2df2f65c89bd5d42b326 Mon Sep 17 00:00:00 2001 From: David Legg Date: Mon, 5 Feb 2024 18:09:01 -0800 Subject: [PATCH 6/6] Pass expiry through LBCS and clampedIntegrate Removes the use of eraseExpiry in LinearBoundaryConsistencySolver. This was a patch put in place to support feedback loops involving the solver, but it erases too much information. Leaving the expiry in place and erasing it only when actually building a feedback loop lets downstream components like approximations take that expiry into account. Also reconfigures the use of eraseExpiry in clampedIntegrate to pass expiry information through. --- .../polynomial/LinearBoundaryConsistencySolver.java | 3 +-- .../modeling/polynomial/PolynomialResources.java | 9 +++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/LinearBoundaryConsistencySolver.java b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/LinearBoundaryConsistencySolver.java index 69cbdb5c1d..11515d3be9 100644 --- a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/LinearBoundaryConsistencySolver.java +++ b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/LinearBoundaryConsistencySolver.java @@ -27,7 +27,6 @@ import static gov.nasa.jpl.aerie.contrib.streamline.core.Expiring.neverExpiring; import static gov.nasa.jpl.aerie.contrib.streamline.core.Reactions.whenever; import static gov.nasa.jpl.aerie.contrib.streamline.core.monads.ExpiringMonad.bind; -import static gov.nasa.jpl.aerie.contrib.streamline.core.Resources.eraseExpiry; import static gov.nasa.jpl.aerie.contrib.streamline.debugging.Context.contextualized; import static gov.nasa.jpl.aerie.contrib.streamline.debugging.Dependencies.addDependency; import static gov.nasa.jpl.aerie.contrib.streamline.debugging.Naming.getName; @@ -314,7 +313,7 @@ private Stream directionalConstraints(Variable constraine // Expiry for driven terms is captured by re-solving rather than expiring the solution. // If solver has a feedback loop from last iteration (which is common) // feeding that expiry in here can loop the solver forever. - var result = eraseExpiry(drivenTerm); + var result = drivenTerm; for (var drivingVariable : drivingVariables) { var scale = controlledTerm.get(drivingVariable); var domain = domains.get(drivingVariable); diff --git a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/PolynomialResources.java b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/PolynomialResources.java index fc12ae44e2..40c26b3ecb 100644 --- a/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/PolynomialResources.java +++ b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/PolynomialResources.java @@ -364,6 +364,7 @@ public static ClampedIntegrateResult clampedIntegrate( Resource integrand, Resource lowerBound, Resource upperBound, double startingValue) { LinearBoundaryConsistencySolver rateSolver = new LinearBoundaryConsistencySolver("clampedIntegrate rate solver"); var integral = resource(polynomial(startingValue)); + var neverExpiringIntegral = eraseExpiry(integral); // Solve for the rate as a function of value var overflowRate = rateSolver.variable("overflowRate", Domain::lowerBound); @@ -377,11 +378,11 @@ public static ClampedIntegrateResult clampedIntegrate( // Set up rate clamping conditions var integrandUB = choose( - greaterThanOrEquals(integral, upperBound), + greaterThanOrEquals(neverExpiringIntegral, upperBound), differentiate(upperBound), constant(Double.POSITIVE_INFINITY)); var integrandLB = choose( - lessThanOrEquals(integral, lowerBound), + lessThanOrEquals(neverExpiringIntegral, lowerBound), differentiate(lowerBound), constant(Double.NEGATIVE_INFINITY)); @@ -390,10 +391,10 @@ public static ClampedIntegrateResult clampedIntegrate( // Use a simple feedback loop on volumes to do the integration and clamping. // Clamping here takes care of discrete under-/overflows and overshooting bounds due to discrete time steps. - var clampedCell = clamp(integral, lowerBound, upperBound); + var clampedCell = clamp(neverExpiringIntegral, lowerBound, upperBound); var correctedCell = map(clampedCell, rate.resource(), (v, r) -> r.integral(v.extract())); // Use the corrected integral values to set volumes, but erase expiry information in the process to avoid loops - forward(eraseExpiry(correctedCell), integral); + forward(correctedCell, integral); name(integral, "Clamped Integral (%s)", integrand); name(overflowRate.resource(), "Overflow of %s", integral);