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) { 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, 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; } }; } 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/Polynomial.java b/contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/Polynomial.java index 7e00dd83a9..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 @@ -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 { @@ -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/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); 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();