Add streamline resource framework

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/CellRefV2.java

@@ -0,0 +1,133 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.ErrorCatchingMonad;
+import gov.nasa.jpl.aerie.merlin.framework.CellRef;
+import gov.nasa.jpl.aerie.merlin.protocol.model.CellType;
+import gov.nasa.jpl.aerie.merlin.protocol.model.EffectTrait;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+import java.util.List;
+import java.util.function.BinaryOperator;
+import java.util.stream.Collectors;
+import java.util.stream.Stream;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching.failure;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Expiring.expiring;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Labelled.labelled;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.ZERO;
+
+public final class CellRefV2 {
+  private CellRefV2() {}
+
+  /**
+   * Allocate a new resource with an explicitly given effect type and effect trait.
+   */
+  public static <D extends Dynamics<?, D>, E extends DynamicsEffect<D>> CellRef<Labelled<E>, Cell<D>> allocate(ErrorCatching<Expiring<D>> initialDynamics, EffectTrait<Labelled<E>> effectTrait) {
+    return CellRef.allocate(new Cell<>(initialDynamics), new CellType<>() {
+      @Override
+      public EffectTrait<Labelled<E>> getEffectType() {
+        return effectTrait;
+      }
+
+      @Override
+      public Cell<D> duplicate(Cell<D> cell) {
+        return new Cell<>(cell.initialDynamics, cell.dynamics, cell.elapsedTime);
+      }
+
+      @Override
+      public void apply(Cell<D> cell, Labelled<E> effect) {
+        cell.initialDynamics = effect.data().apply(cell.dynamics).match(
+            ErrorCatching::success,
+            error -> failure(new RuntimeException(
+                "Applying '%s' failed.".formatted(effect.name()), error)));
+        cell.dynamics = cell.initialDynamics;
+        cell.elapsedTime = ZERO;
+      }
+
+      @Override
+      public void step(Cell<D> cell, Duration duration) {
+        // Avoid accumulated round-off error in imperfect stepping
+        // by always stepping up from the initial dynamics
+        cell.elapsedTime = cell.elapsedTime.plus(duration);
+        cell.dynamics = ErrorCatchingMonad.map(cell.initialDynamics, d ->
+            expiring(d.data().step(cell.elapsedTime), d.expiry().minus(cell.elapsedTime)));
+      }
+    });
+  }
+
+  public static <D extends Dynamics<?, D>> EffectTrait<Labelled<DynamicsEffect<D>>> noncommutingEffects() {
+    return resolvingConcurrencyBy((left, right) -> x -> {
+          throw new UnsupportedOperationException(
+              "Concurrent effects are not supported on this resource.");
+        });
+  }
+
+  public static <D extends Dynamics<?, D>> EffectTrait<Labelled<DynamicsEffect<D>>> commutingEffects() {
+    return resolvingConcurrencyBy((left, right) -> x -> right.apply(left.apply(x)));
+  }
+
+  public static <D extends Dynamics<?, D>> EffectTrait<Labelled<DynamicsEffect<D>>> autoEffects() {
+    return resolvingConcurrencyBy((left, right) -> x -> {
+      final var lrx = left.apply(right.apply(x));
+      final var rlx = right.apply(left.apply(x));
+      if (lrx.equals(rlx)) {
+        return lrx;
+      } else {
+        throw new UnsupportedOperationException(
+            "Detected non-commuting concurrent effects!");
+      }
+    });
+  }
+
+  public static <D extends Dynamics<?, D>> EffectTrait<Labelled<DynamicsEffect<D>>> resolvingConcurrencyBy(BinaryOperator<DynamicsEffect<D>> combineConcurrent) {
+    return new EffectTrait<>() {
+      @Override
+      public Labelled<DynamicsEffect<D>> empty() {
+        return labelled("No-op", x -> x);
+      }
+
+      @Override
+      public Labelled<DynamicsEffect<D>> sequentially(final Labelled<DynamicsEffect<D>> prefix, final Labelled<DynamicsEffect<D>> suffix) {
+        return new Labelled<>(
+            x -> suffix.data().apply(prefix.data().apply(x)),
+            "(%s) then (%s)".formatted(prefix.name(), suffix.name()));
+      }
+
+      @Override
+      public Labelled<DynamicsEffect<D>> concurrently(final Labelled<DynamicsEffect<D>> left, final Labelled<DynamicsEffect<D>> right) {
+        try {
+          final DynamicsEffect<D> combined = combineConcurrent.apply(left.data(), right.data());
+          return new Labelled<>(
+              x -> {
+                try {
+                  return combined.apply(x);
+                } catch (Exception e) {
+                  return failure(e);
+                }
+              },
+              "(%s) and (%s)".formatted(left.name(), right.name()));
+        } catch (Throwable e) {
+          return new Labelled<>(
+              $ -> failure(e),
+              "Failed to combine concurrent effects: (%s) and (%s)".formatted(left.name(), right.name()));
+        }
+      }
+    };
+  }
+
+  public static class Cell<D> {
+    public ErrorCatching<Expiring<D>> initialDynamics;
+    public ErrorCatching<Expiring<D>> dynamics;
+    public Duration elapsedTime;
+
+    public Cell(ErrorCatching<Expiring<D>> dynamics) {
+      this(dynamics, dynamics, ZERO);
+    }
+
+    public Cell(ErrorCatching<Expiring<D>> initialDynamics, ErrorCatching<Expiring<D>> dynamics, Duration elapsedTime) {
+      this.initialDynamics = initialDynamics;
+      this.dynamics = dynamics;
+      this.elapsedTime = elapsedTime;
+    }
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/CellResource.java

@@ -0,0 +1,131 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.ErrorCatchingMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.debugging.Context;
+import gov.nasa.jpl.aerie.merlin.framework.CellRef;
+import gov.nasa.jpl.aerie.contrib.streamline.core.CellRefV2.Cell;
+import gov.nasa.jpl.aerie.merlin.framework.Scoped;
+import gov.nasa.jpl.aerie.merlin.protocol.model.EffectTrait;
+
+import java.util.LinkedList;
+import java.util.List;
+import java.util.function.Supplier;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellRefV2.allocate;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellRefV2.autoEffects;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Labelled.labelled;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad.unit;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Unit.UNIT;
+import static java.util.stream.Collectors.joining;
+
+/**
+ * Resource which is backed directly by a cell.
+ * Effects can be applied to this resource.
+ * Effect names are augmented with this resource's name(s).
+ */
+public interface CellResource<D extends Dynamics<?, D>> extends Resource<D> {
+  void emit(Labelled<DynamicsEffect<D>> effect);
+  default void emit(String effectName, DynamicsEffect<D> effect) {
+    emit(labelled(effectName, effect));
+  }
+  default void emit(DynamicsEffect<D> effect) {
+    emit("anonymous effect", effect);
+  }
+
+  static <D extends Dynamics<?, D>> CellResource<D> cellResource(D initial) {
+    return cellResource(unit(initial));
+  }
+
+  static <D extends Dynamics<?, D>> CellResource<D> cellResource(D initial, EffectTrait<Labelled<DynamicsEffect<D>>> effectTrait) {
+    return cellResource(unit(initial), effectTrait);
+  }
+
+  static <D extends Dynamics<?, D>> CellResource<D> cellResource(ErrorCatching<Expiring<D>> initial) {
+    return cellResource(initial, autoEffects());
+  }
+
+  static <D extends Dynamics<?, D>> CellResource<D> cellResource(ErrorCatching<Expiring<D>> initial, EffectTrait<Labelled<DynamicsEffect<D>>> effectTrait) {
+    return new CellResource<>() {
+      // Use autoEffects for a generic CellResource, on the theory that most resources
+      // have relatively few effects, and even fewer concurrent effects, so this is performant enough.
+      // If that doesn't hold, a more specialized solution can be constructed directly.
+      private final CellRef<Labelled<DynamicsEffect<D>>, Cell<D>> cell = allocate(initial, effectTrait);
+      private final List<String> names = new LinkedList<>();
+
+      @Override
+      public void emit(final Labelled<DynamicsEffect<D>> effect) {
+        cell.emit(labelled(augmentEffectName(effect.name()), effect.data()));
+      }
+
+      @Override
+      public ErrorCatching<Expiring<D>> getDynamics() {
+        return cell.get().dynamics;
+      }
+
+      @Override
+      public void registerName(final String name) {
+        names.add(name);
+      }
+
+      private String augmentEffectName(String effectName) {
+        var resourceName = switch (names.size()) {
+          case 0 -> "anonymous resource";
+          case 1 -> names.get(0);
+          default -> names.get(0) + " (aka. %s)".formatted(String.join(", ", names.subList(1, names.size())));
+        };
+        return effectName + " on " + resourceName + Context.get().stream().map(c -> " during " + c).collect(joining());
+      }
+    };
+  }
+
+  static <D extends Dynamics<?, D>> CellResource<D> staticallyCreated(Supplier<CellResource<D>> constructor) {
+    return new CellResource<>() {
+      private CellResource<D> delegate = constructor.get();
+
+      @Override
+      public void emit(final Labelled<DynamicsEffect<D>> effect) {
+        actOnCell(() -> delegate.emit(effect));
+      }
+
+      @Override
+      public ErrorCatching<Expiring<D>> getDynamics() {
+        // Keep the field access using () -> ... form, don't simplify to delegate::getDynamics
+        // Simplifying will access delegate before calling actOnCell, failing if we need to re-allocate delegate.
+        return actOnCell(() -> delegate.getDynamics());
+      }
+
+      @Override
+      public void registerName(final String name) {
+        delegate.registerName(name);
+      }
+
+      private void actOnCell(Runnable action) {
+        actOnCell(() -> {
+          action.run();
+          return UNIT;
+        });
+      }
+
+      private <R> R actOnCell(Supplier<R> action) {
+        try {
+          return action.get();
+        } catch (Scoped.EmptyDynamicCellException | IllegalArgumentException e) {
+          // If we're running unit tests, several simulations can happen without reloading the Resources class.
+          // In that case, we'll have discarded the clock resource we were using, and get the above exception.
+          // REVIEW: Is there a cleaner way to make sure this resource gets (re-)initialized?
+          delegate = constructor.get();
+          return action.get();
+        }
+      }
+    };
+  }
+
+  static <D extends Dynamics<?, D>> void set(CellResource<D> resource, D newDynamics) {
+    resource.emit("Set " + newDynamics, DynamicsMonad.effect(x -> newDynamics));
+  }
+
+  static <D extends Dynamics<?, D>> void set(CellResource<D> resource, Expiring<D> newDynamics) {
+    resource.emit("Set " + newDynamics, ErrorCatchingMonad.<Expiring<D>, Expiring<D>>lift($ -> newDynamics)::apply);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Dynamics.java

@@ -0,0 +1,9 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+public interface Dynamics<V, D extends Dynamics<V, D>> {
+    V extract();
+
+    D step(Duration t);
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/DynamicsEffect.java

@@ -0,0 +1,5 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+public interface DynamicsEffect<D extends Dynamics<?, D>> {
+    ErrorCatching<Expiring<D>> apply(ErrorCatching<Expiring<D>> dynamics);
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/ErrorCatching.java

@@ -0,0 +1,43 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.ErrorCatchingMonad;
+
+import java.util.function.Function;
+
+public sealed interface ErrorCatching<T> {
+  <R> R match(Function<T, R> onSuccess, Function<Throwable, R> onError);
+
+  static <T> ErrorCatching<T> success(T result) {
+    return new Success<>(result);
+  }
+
+  static <T> ErrorCatching<T> failure(Throwable exception) {
+    return new Failure<>(exception);
+  }
+
+  default <R> ErrorCatching<R> map(Function<T, R> f) {
+    return ErrorCatchingMonad.map(this, f);
+  }
+
+  default T getOrThrow() {
+    return match(
+        Function.identity(),
+        e -> {
+          throw new RuntimeException(e);
+        });
+  }
+
+  record Success<T>(T result) implements ErrorCatching<T> {
+    @Override
+    public <R> R match(final Function<T, R> onSuccess, final Function<Throwable, R> onError) {
+      return onSuccess.apply(result);
+    }
+  }
+
+  record Failure<T>(Throwable exception) implements ErrorCatching<T> {
+    @Override
+    public <R> R match(final Function<T, R> onSuccess, final Function<Throwable, R> onError) {
+      return onError.apply(exception);
+    }
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Expiring.java

@@ -0,0 +1,19 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Expiry.NEVER;
+
+public record Expiring<D>(D data, Expiry expiry) {
+  public static <D> Expiring<D> expiring(D data, Expiry expiry) {
+    return new Expiring<>(data, expiry);
+  }
+
+  public static <D> Expiring<D> neverExpiring(D data) {
+    return expiring(data, NEVER);
+  }
+
+  public static <D> Expiring<D> expiring(D data, Duration expiry) {
+    return expiring(data, Expiry.at(expiry));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Expiry.java

@@ -0,0 +1,68 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+import java.util.Optional;
+import java.util.stream.Stream;
+
+public record Expiry(Optional<Duration> value) {
+  public static Expiry NEVER = expiry(Optional.empty());
+
+  public static Expiry at(Duration t) {
+    return expiry(Optional.of(t));
+  }
+
+  public static Expiry expiry(Optional<Duration> value) {
+    return new Expiry(value);
+  }
+
+  public Expiry or(Expiry other) {
+    return expiry(
+        Stream.concat(value().stream(), other.value().stream()).reduce(Duration::min));
+  }
+
+  public Expiry minus(Duration t) {
+    return expiry(value().map(v -> v.minus(t)));
+  }
+
+  public boolean isNever() {
+    return value().isEmpty();
+  }
+
+  public int compareTo(Expiry other) {
+    if (this.isNever()) {
+      if (other.isNever()) {
+        return 0;
+      } else {
+        return 1;
+      }
+    } else {
+      if (other.isNever()) {
+        return -1;
+      } else {
+        return this.value().get().compareTo(other.value().get());
+      }
+    }
+  }
+
+  public boolean shorterThan(Expiry other) {
+    return this.compareTo(other) < 0;
+  }
+
+  public boolean noShorterThan(Expiry other) {
+    return this.compareTo(other) >= 0;
+  }
+
+  public boolean longerThan(Expiry other) {
+    return this.compareTo(other) > 0;
+  }
+
+  public boolean noLongerThan(Expiry other) {
+    return this.compareTo(other) <= 0;
+  }
+
+  @Override
+  public String toString() {
+    return value.map(Duration::toString).orElse("NEVER");
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Labelled.java

@@ -0,0 +1,25 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+import java.util.Objects;
+
+/**
+ * Attaches name and context to a datum.
+ */
+public record Labelled<V>(V data, String name) {
+  public static <V> Labelled<V> labelled(String name, V data) {
+    return new Labelled<>(data, name);
+  }
+
+  @Override
+  public boolean equals(final Object o) {
+    if (this == o) return true;
+    if (o == null || getClass() != o.getClass()) return false;
+    Labelled<?> labelled = (Labelled<?>) o;
+    return Objects.equals(data, labelled.data);
+  }
+
+  @Override
+  public int hashCode() {
+    return Objects.hash(data);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Reactions.java

@@ -0,0 +1,49 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+import gov.nasa.jpl.aerie.merlin.framework.Condition;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+import java.util.function.Consumer;
+import java.util.function.Supplier;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Resources.dynamicsChange;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.DiscreteResources.when;
+import static gov.nasa.jpl.aerie.merlin.framework.ModelActions.delay;
+import static gov.nasa.jpl.aerie.merlin.framework.ModelActions.replaying;
+import static gov.nasa.jpl.aerie.merlin.framework.ModelActions.spawn;
+import static gov.nasa.jpl.aerie.merlin.framework.ModelActions.waitUntil;
+
+public final class Reactions {
+  private Reactions() {}
+
+  public static void whenever(Resource<Discrete<Boolean>> conditionResource, Runnable action) {
+    whenever(when(conditionResource), action);
+  }
+
+  public static void whenever(Condition condition, Runnable action) {
+    whenever(() -> condition, action);
+  }
+
+  public static void whenever(Supplier<Condition> trigger, Runnable action) {
+    final Condition condition = trigger.get();
+    // Use replaying tasks to avoid threading overhead.
+    spawn(replaying(() -> {
+      waitUntil(condition);
+      action.run();
+      // Trampoline off this task to avoid replaying.
+      whenever(trigger, action);
+    }));
+  }
+
+  // Special case for dynamicsChange condition, since it's non-obvious that this needs to be run in lambda form
+  public static <D extends Dynamics<?, D>> void wheneverDynamicsChange(Resource<D> resource, Consumer<ErrorCatching<Expiring<D>>> reaction) {
+    whenever(() -> dynamicsChange(resource), () -> reaction.accept(resource.getDynamics()));
+  }
+
+  public static void every(Duration period, Runnable action) {
+    delay(period);
+    action.run();
+    spawn(replaying(() -> every(period, action)));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Resource.java

@@ -0,0 +1,8 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+public interface Resource<D> {
+  ErrorCatching<Expiring<D>> getDynamics();
+
+  // By default, resources don't track their names
+  default void registerName(String name) {}
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/Resources.java

@@ -0,0 +1,162 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core;
+
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.clocks.Clock;
+import gov.nasa.jpl.aerie.merlin.framework.Condition;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Unit;
+
+import java.util.List;
+import java.util.Optional;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellResource.cellResource;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellResource.staticallyCreated;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Reactions.wheneverDynamicsChange;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.clocks.Clock.clock;
+import static gov.nasa.jpl.aerie.merlin.framework.ModelActions.*;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.ZERO;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete.discrete;
+
+/**
+ * Utility methods for {@link Resource}s.
+ */
+public final class Resources {
+  private Resources() {}
+
+  /**
+   * Ensure that Resources are initialized.
+   *
+   * <p>
+   *   This method needs to be called during simulation initialization.
+   *   This method is idempotent; calling it multiple times is the same as calling it once.
+   * </p>
+   */
+  public static void init() {
+    currentTime();
+  }
+
+  private static final Resource<Clock> CLOCK = staticallyCreated(() -> cellResource(clock(ZERO)));
+  public static Duration currentTime() {
+    return currentValue(CLOCK);
+  }
+
+  public static <D> D currentData(Resource<D> resource) {
+    return resource.getDynamics().getOrThrow().data();
+  }
+
+  public static <V, D extends Dynamics<V, D>> V currentValue(Resource<D> resource) {
+    return currentData(resource).extract();
+  }
+
+  public static <D extends Dynamics<?, D>> Condition dynamicsChange(Resource<D> resource) {
+    final var startingDynamics = resource.getDynamics();
+    final Duration startTime = currentTime();
+    return (positive, atEarliest, atLatest) -> {
+      var currentDynamics = resource.getDynamics();
+      boolean haveChanged = startingDynamics.match(
+          start -> currentDynamics.match(
+              current -> !current.data().equals(start.data().step(currentTime().minus(startTime))),
+              ignored -> true),
+          startException -> currentDynamics.match(
+              ignored -> true,
+              // Use semantic comparison for exceptions, since derivation can generate the exception each invocation.
+              currentException -> !(startException.getClass().equals(currentException.getClass())
+                                    && startException.getMessage().equals(currentException.getMessage()))));
+
+      return positive == haveChanged
+          ? Optional.of(atEarliest)
+          : positive
+            ? currentDynamics.match(
+                expiring -> expiring.expiry().value().filter(atLatest::noShorterThan),
+                exception -> Optional.empty())
+            : Optional.empty();
+    };
+  }
+
+  public static <D extends Dynamics<?, D>> Condition dynamicsChange(List<Resource<D>> resources) {
+    assert resources.size() > 0;
+    var result = dynamicsChange(resources.get(0));
+    for (Resource<D> r : resources) {
+      result = result.or(dynamicsChange(r));
+    }
+    return result;
+  }
+
+  /**
+   * Cache this resource in a resource.
+   *
+   * <p>
+   *   Updates the resource when resource changes dynamics.
+   *   This can be used to isolate a resource from effects
+   *   which don't change the dynamics, so Aerie samples that
+   *   resource only when strictly necessary.
+   * </p>
+   * <p>
+   *   This introduces a small delay in deriving values.
+   *   Specifically, the cached version of a resource changes two
+   *   simulation engine cycles after its uncached version.
+   *   It will show up as the same instant in the results,
+   *   but beware that it could be momentarily out-of-sync
+   *   with its sources during simulation.
+   * </p>
+   */
+  public static <D extends Dynamics<?, D>> Resource<D> cache(Resource<D> resource) {
+    var cell = cellResource(resource.getDynamics());
+    wheneverDynamicsChange(resource, newDynamics -> cell.emit($ -> newDynamics));
+    return cell;
+  }
+
+  /**
+   * Signal discrete changes in this resource's dynamics.
+   *
+   * <p>
+   *   For Aerie's resource sampling to work correctly,
+   *   there must be an effect every time a resource changes dynamics.
+   *   For most derived resources, this happens automatically.
+   *   For some derivations, though, continuous changes in the source state
+   *   can cause discrete changes in the result.
+   *   For example, imagine a continuous numeric resource R,
+   *   and a derived resource S := "R > 0".
+   *   If R changes continuously from positive to negative,
+   *   then S changes discretely from true to false, *without* an effect.
+   *   If used directly, Aerie would not re-sample S at this time.
+   *   This method emits a trivial effect when this happens so that S
+   *   *would* be resampled correctly.
+   * </p>
+   * <p>
+   *   Unlike {@link Resources#cache}, this method does *not* introduce
+   *   a delay between the source and derived resources.
+   *   Signalling resources use a resource "in parallel" rather than "in series"
+   *   with the derivation process, thereby avoiding the delay.
+   *   Like regular derived resources, signalling resources calculate their value
+   *   through the derivation every time they are sampled.
+   * </p>
+   */
+  // REVIEW: Suggestion from Jonathan Castello to remove this method
+  // in favor of allowing resources to report expiry information directly.
+  // This would be cleaner and potentially more performant.
+  public static <D extends Dynamics<?, D>> Resource<D> signalling(Resource<D> resource) {
+    var cell = cellResource(discrete(Unit.UNIT));
+    wheneverDynamicsChange(resource, ignored -> cell.emit($ -> $));
+    return () -> {
+      cell.getDynamics();
+      return resource.getDynamics();
+    };
+  }
+
+  public static <D extends Dynamics<?, D>> Resource<D> shift(Resource<D> resource, Duration interval, D initialDynamics) {
+    var cell = cellResource(initialDynamics);
+    delayedSet(cell, resource.getDynamics(), interval);
+    wheneverDynamicsChange(resource, newDynamics ->
+        delayedSet(cell, newDynamics, interval));
+    return cell;
+  }
+
+  private static <D extends Dynamics<?, D>> void delayedSet(
+      CellResource<D> cell, ErrorCatching<Expiring<D>> newDynamics, Duration interval)
+  {
+    spawn(replaying(() -> {
+      delay(interval);
+      cell.emit($ -> newDynamics);
+    }));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/monads/DynamicsMonad.java

@@ -0,0 +1,45 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Dynamics;
+import gov.nasa.jpl.aerie.contrib.streamline.core.DynamicsEffect;
+import gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+
+import java.util.function.Function;
+
+public final class DynamicsMonad {
+  private DynamicsMonad() {}
+
+  public static <A> ErrorCatching<Expiring<A>> unit(A a) {
+    return ExpiringMonadTransformer.unit(ErrorCatchingMonad::unit, a);
+  }
+
+  public static <A, B> ErrorCatching<Expiring<B>> bind(ErrorCatching<Expiring<A>> a, Function<A, ErrorCatching<Expiring<B>>> f) {
+    return ExpiringMonadTransformer.<A, ErrorCatching<Expiring<A>>, B, ErrorCatching<Expiring<B>>>bind(
+        ErrorCatchingMonad::unit,
+        ErrorCatchingMonad::bind,
+        ErrorCatchingMonad::bind,
+        a,
+        f);
+  }
+
+  // Convenient methods defined in terms of bind and unit:
+
+  public static <A, B> ErrorCatching<Expiring<B>> map(ErrorCatching<Expiring<A>> a, Function<A, B> f) {
+    return bind(a, f.andThen(DynamicsMonad::unit));
+  }
+
+  public static <A, B> Function<ErrorCatching<Expiring<A>>, ErrorCatching<Expiring<B>>> lift(Function<A, B> f) {
+    return a -> map(a, f);
+  }
+
+  // Not fully monadic since we intentionally ignore expiry information, but useful nonetheless.
+
+  public static <A extends Dynamics<?, A>> DynamicsEffect<A> effect(Function<A, A> f) {
+    return bindEffect(f.andThen(DynamicsMonad::unit));
+  }
+
+  public static <A extends Dynamics<?, A>> DynamicsEffect<A> bindEffect(Function<A, ErrorCatching<Expiring<A>>> f) {
+    return ea -> ErrorCatchingMonad.bind(ea, a -> f.apply(a.data()));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/monads/ErrorCatchingMonad.java

@@ -0,0 +1,31 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching;
+
+import java.util.function.Function;
+
+public final class ErrorCatchingMonad {
+  private ErrorCatchingMonad() {}
+
+  public static <A> ErrorCatching<A> unit(A a) {
+    return ErrorCatchingMonadTransformer.unit(IdentityMonad::unit, a);
+  }
+
+  public static <A, B> ErrorCatching<B> bind(ErrorCatching<A> a, Function<A, ErrorCatching<B>> f) {
+    return ErrorCatchingMonadTransformer.<A, ErrorCatching<A>, B, ErrorCatching<B>>bind(
+        IdentityMonad::unit,
+        IdentityMonad::bind,
+        a,
+        f);
+  }
+
+  // Convenient methods defined in terms of bind and unit:
+
+  public static <A, B> ErrorCatching<B> map(ErrorCatching<A> a, Function<A, B> f) {
+    return bind(a, f.andThen(ErrorCatchingMonad::unit));
+  }
+
+  public static <A, B> Function<ErrorCatching<A>, ErrorCatching<B>> lift(Function<A, B> f) {
+    return a -> map(a, f);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/monads/ErrorCatchingMonadTransformer.java

@@ -0,0 +1,41 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching;
+
+import java.util.function.BiFunction;
+import java.util.function.Function;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching.*;
+
+/**
+ * Monad transformer for {@link ErrorCatching}, M A -> M {@link ErrorCatching}<A>
+ *
+ * <p>
+ *   Bind for this operation does two jobs:
+ *   First, if a computation fails by throwing an exception, that exception is caught so as not to crash the simulation.
+ *   Second, if a prior computation failed, the exception is passed through and further computations are skipped.
+ * </p>
+ */
+public final class ErrorCatchingMonadTransformer {
+  private ErrorCatchingMonadTransformer() {}
+
+  public static <A, MEA> MEA unit(Function<ErrorCatching<A>, MEA> mUnit, A a) {
+    return mUnit.apply(success(a));
+  }
+
+  public static <A, MEA, B, MEB> MEB bind(
+      Function<ErrorCatching<B>, MEB> mUnit,
+      BiFunction<MEA, Function<ErrorCatching<A>, MEB>, MEB> mBind,
+      MEA mea,
+      Function<A, MEB> f) {
+    return mBind.apply(mea, eb -> eb.match(
+        a -> {
+          try {
+            return f.apply(a);
+          } catch (Throwable e) {
+            return mUnit.apply(failure(e));
+          }
+        },
+        e -> mUnit.apply(failure(e))));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/monads/ErrorCatchingToResourceMonad.java

@@ -0,0 +1,21 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+
+import java.util.function.Function;
+
+public final class ErrorCatchingToResourceMonad {
+  private ErrorCatchingToResourceMonad() {}
+
+  public static <A> Resource<A> unit(ErrorCatching<Expiring<A>> a) {
+    return () -> a;
+  }
+
+  public static <A, B> Resource<B> bind(
+      Resource<A> a,
+      Function<ErrorCatching<Expiring<A>>, Resource<B>> f) {
+    return () -> f.apply(a.getDynamics()).getDynamics();
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/monads/ExpiringMonad.java

@@ -0,0 +1,37 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+
+import java.util.function.Function;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Expiring.expiring;
+
+/**
+ * The {@link Expiring} monad, which demands derived values expire no later than their sources.
+ */
+public final class ExpiringMonad {
+  private ExpiringMonad() {}
+
+  public static <A> Expiring<A> unit(A data) {
+    return ExpiringMonadTransformer.unit(IdentityMonad::unit, data);
+  }
+
+  public static <A, B> Expiring<B> bind(Expiring<A> a, Function<A, Expiring<B>> f) {
+    return ExpiringMonadTransformer.<A, Expiring<A>, B, Expiring<B>>bind(
+        IdentityMonad::unit,
+        IdentityMonad::bind,
+        IdentityMonad::bind,
+        a,
+        f);
+  }
+
+  // Convenient methods defined in terms of bind and unit:
+
+  public static <A, B> Expiring<B> map(Expiring<A> a, Function<A, B> f) {
+    return bind(a, f.andThen(ExpiringMonad::unit));
+  }
+
+  public static <A, B> Function<Expiring<A>, Expiring<B>> lift(Function<A, B> f) {
+    return a -> map(a, f);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/monads/ExpiringMonadTransformer.java

@@ -0,0 +1,35 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+
+import java.util.function.BiFunction;
+import java.util.function.Function;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Expiring.*;
+
+/**
+ * Monad transformer for {@link Expiring}, M A -> M {@link Expiring}<A>
+ *
+ * <p>
+ *   Bind for this monad ensures that results expire no later than the values they're derived from.
+ * </p>
+ */
+public final class ExpiringMonadTransformer {
+  private ExpiringMonadTransformer() {}
+
+  public static <A, MEA> MEA unit(Function<Expiring<A>, MEA> mUnit, A a) {
+    return mUnit.apply(neverExpiring(a));
+  }
+
+  public static <A, MEA, B, MEB> MEB bind(
+      Function<Expiring<B>, MEB> mUnit,
+      BiFunction<MEA, Function<Expiring<A>, MEB>, MEB> mBind1,
+      BiFunction<MEB, Function<Expiring<B>, MEB>, MEB> mBind2,
+      MEA mea,
+      Function<A, MEB> f) {
+    // TODO: for performance, this could take mMap instead of mUnit and mBind2
+    return mBind1.apply(mea, ea ->
+        mBind2.apply(f.apply(ea.data()), eb ->
+            mUnit.apply(expiring(eb.data(), ea.expiry().or(eb.expiry())))));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/monads/ExpiringToResourceMonad.java

@@ -0,0 +1,28 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+
+import java.util.function.Function;
+
+public final class ExpiringToResourceMonad {
+  private ExpiringToResourceMonad() {}
+
+  public static <A> Resource<A> unit(Expiring<A> a) {
+    return ErrorCatchingMonadTransformer.unit(ErrorCatchingToResourceMonad::unit, a);
+  }
+
+  public static <A, B> Resource<B> bind(Resource<A> a, Function<Expiring<A>, Resource<B>> f) {
+    return ErrorCatchingMonadTransformer.<Expiring<A>, Resource<A>, Expiring<B>, Resource<B>>bind(
+        ErrorCatchingToResourceMonad::unit,
+        ErrorCatchingToResourceMonad::bind,
+        a,
+        f);
+  }
+
+  // Convenient methods defined in terms of bind and unit:
+
+  public static <A, B> Resource<B> map(Resource<A> a, Function<Expiring<A>, Expiring<B>> f) {
+    return bind(a, f.andThen(ExpiringToResourceMonad::unit));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/monads/IdentityMonad.java

@@ -0,0 +1,18 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core.monads;
+
+import java.util.function.Function;
+
+/**
+ * The trivial monad A -> A
+ */
+public final class IdentityMonad {
+  private IdentityMonad() {}
+
+  public static <A> A unit(A a) {
+                                return a;
+                                         }
+
+  public static <A, B> B bind(A a, Function<A, B> f) {
+                                                     return f.apply(a);
+    }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/core/monads/ResourceMonad.java

@@ -0,0 +1,59 @@
+package gov.nasa.jpl.aerie.contrib.streamline.core.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+import org.apache.commons.lang3.function.TriFunction;
+
+import java.util.function.BiFunction;
+import java.util.function.Function;
+
+/**
+ * Monad A -> Resource<A>.
+ * This is the primary monad for model authors,
+ * handling both expiry and "stitching together" of derived resources.
+ */
+public final class ResourceMonad {
+  private ResourceMonad() {}
+
+  public static <A> Resource<A> unit(A a) {
+    return ExpiringMonadTransformer.unit(ExpiringToResourceMonad::unit, a);
+  }
+
+  public static <A, B> Resource<B> bind(Resource<A> a, Function<A, Resource<B>> f) {
+    return ExpiringMonadTransformer.<A, Resource<A>, B, Resource<B>>bind(
+        ExpiringToResourceMonad::unit,
+        ExpiringToResourceMonad::bind,
+        ExpiringToResourceMonad::bind,
+        a,
+        f);
+  }
+
+  // Convenient methods defined in terms of bind and unit:
+
+  public static <A, B> Resource<B> map(Resource<A> a, Function<A, B> f) {
+    return bind(a, f.andThen(ResourceMonad::unit));
+  }
+
+  public static <A, B, C> Resource<C> map(Resource<A> a, Resource<B> b, BiFunction<A, B, C> f) {
+    return bind(a, a$ -> map(b, b$ -> f.apply(a$, b$)));
+  }
+
+  public static <A, B, C, D> Resource<D> map(Resource<A> a, Resource<B> b, Resource<C> c, TriFunction<A, B, C, D> f) {
+    return bind(a, a$ -> map(b, c, (b$, c$) -> f.apply(a$, b$, c$)));
+  }
+
+  public static <A, B, C> Resource<C> bind(Resource<A> a, Resource<B> b, BiFunction<A, B, Resource<C>> f) {
+    return bind(a, a$ -> bind(b, b$ -> f.apply(a$, b$)));
+  }
+
+  public static <A, B, C, D> Resource<D> bind(Resource<A> a, Resource<B> b, Resource<C> c, TriFunction<A, B, C, Resource<D>> f) {
+    return bind(a, a$ -> bind(b, c, (b$, c$) -> f.apply(a$, b$, c$)));
+  }
+
+  public static <A, B> Function<Resource<A>, Resource<B>> lift(Function<A, B> f) {
+    return a -> map(a, f);
+  }
+
+  public static <A, B, C> BiFunction<Resource<A>, Resource<B>, Resource<C>> lift(BiFunction<A, B, C> f) {
+    return (a, b) -> map(a, b, f);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/debugging/Context.java

@@ -0,0 +1,147 @@
+package gov.nasa.jpl.aerie.contrib.streamline.debugging;
+
+import gov.nasa.jpl.aerie.merlin.protocol.types.Unit;
+
+import java.util.ArrayDeque;
+import java.util.Deque;
+import java.util.List;
+import java.util.function.Supplier;
+
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Unit.UNIT;
+
+/**
+ * Thread-local scope-bound description of the current context.
+ */
+public final class Context {
+  private Context() {}
+
+  private static final ThreadLocal<Deque<String>> contexts = ThreadLocal.withInitial(ArrayDeque::new);
+
+  /**
+   * @see Context#inContext(String, Supplier)
+   */
+  public static void inContext(String contextName, Runnable action) {
+    inContext(contextName, asSupplier(action));
+  }
+
+  /**
+   * Run action in a globally-visible context.
+   * Contexts stack, and contexts are removed when control leaves action for any reason.
+   */
+  public static <R> R inContext(String contextName, Supplier<R> action) {
+    // Using a thread-local context stack maintains isolation for threaded tasks.
+    try {
+      contexts.get().push(contextName);
+      return action.get();
+      // TODO: Should we add a catch clause here that would add context to the error?
+    } finally {
+      // Doing the tear-down in a finally block maintains isolation for replaying tasks.
+      contexts.get().pop();
+    }
+  }
+
+  /**
+   * @see Context#inContext(List, Supplier)
+   */
+  public static void inContext(List<String> contextStack, Runnable action) {
+    inContext(contextStack, asSupplier(action));
+  }
+
+  /**
+   * Run action in a context stack like that returned by {@link Context#get}.
+   *
+   * <p>
+   *   This can be used to "copy" a context into another task, e.g.
+   *   <pre>
+   *     var context = Context.get();
+   *     spawn(() -> inContext(context, () -> { ... });
+   *   </pre>
+   * </p>
+   *
+   * @see Context#contextualized
+   */
+  public static <R> R inContext(List<String> contextStack, Supplier<R> action) {
+    if (contextStack.isEmpty()) {
+      return action.get();
+    } else {
+      int n = contextStack.size() - 1;
+      return inContext(contextStack.get(n), () ->
+          inContext(contextStack.subList(0, n), action));
+    }
+  }
+
+  /**
+   * @see Context#contextualized(Supplier)
+   */
+  public static Runnable contextualized(Runnable action) {
+    return contextualized(asSupplier(action))::get;
+  }
+
+  /**
+   * Adds the current context into action.
+   *
+   * <p>
+   *   This can be used to contextualize sub-tasks with their parents context:
+   *   <pre>
+   *     inContext("parent", () -> {
+   *       // Capture parent context while calling spawn:
+   *       spawn(contextualized(() -> {
+   *         // Runs child task in context "parent"
+   *       }));
+   *     });
+   *   </pre>
+   * </p>
+   *
+   * @see Context#contextualized(String, Runnable)
+   * @see Context#inContext(List, Runnable)
+   * @see Context#inContext(String, Runnable)
+   */
+  public static <R> Supplier<R> contextualized(Supplier<R> action) {
+    final var context = get();
+    return () -> inContext(context, action);
+  }
+
+  /**
+   * @see Context#contextualized(String, Supplier)
+   */
+  public static Runnable contextualized(String childContext, Runnable action) {
+    return contextualized(childContext, asSupplier(action))::get;
+  }
+
+  /**
+   * Adds the current context into action, as well as an additional child context.
+   *
+   * <p>
+   *   This can be used to contextualize sub-tasks with their parents context:
+   *   <pre>
+   *     inContext("parent", () -> {
+   *       // Capture parent context while calling spawn:
+   *       spawn(contextualized("child", () -> {
+   *         // Runs child task in context ("child", "parent")
+   *       }));
+   *     });
+   *   </pre>
+   * </p>
+   *
+   * @see Context#contextualized(Runnable)
+   * @see Context#inContext(List, Runnable)
+   * @see Context#inContext(String, Runnable)
+   */
+  public static <R> Supplier<R> contextualized(String childContext, Supplier<R> action) {
+    return contextualized(() -> inContext(childContext, action));
+  }
+
+  /**
+   * Returns the list of contexts, from innermost context out.
+   */
+  public static List<String> get() {
+    return contexts.get().stream().toList();
+  }
+
+  private static Supplier<Unit> asSupplier(Runnable action) {
+    return () -> {
+      action.run();
+      return UNIT;
+    };
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/debugging/Tracing.java

@@ -0,0 +1,117 @@
+package gov.nasa.jpl.aerie.contrib.streamline.debugging;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.CellResource;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Dynamics;
+import gov.nasa.jpl.aerie.contrib.streamline.core.DynamicsEffect;
+import gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Labelled;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+import gov.nasa.jpl.aerie.merlin.framework.Condition;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Unit;
+
+import java.util.Stack;
+import java.util.function.Consumer;
+import java.util.function.Supplier;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching.failure;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Resources.currentTime;
+
+/**
+ * Functions for debugging resources by tracing their calculation.
+ */
+public final class Tracing {
+  private Tracing() {}
+
+  private final static Stack<String> activeTracePoints = new Stack<>();
+
+  public static <D> Resource<D> trace(String name, Resource<D> resource) {
+    return traceFull(name, $ -> {}, resource);
+  }
+
+  public static <D> Resource<D> trace(String name, Consumer<D> assertion, Resource<D> resource) {
+    return traceExpiring(name, $ -> assertion.accept($.data()), resource);
+  }
+
+  public static <D> Resource<D> traceExpiring(String name, Consumer<Expiring<D>> assertion, Resource<D> resource) {
+    return traceFull(name, $ -> $.match(d -> {
+      assertion.accept(d);
+      return Unit.UNIT;
+    }, e -> {
+      throw new AssertionError("%s failed while computing".formatted(formatStack()), e);
+    }), resource);
+  }
+
+  public static <D> Resource<D> traceFull(String name, Consumer<ErrorCatching<Expiring<D>>> assertion, Resource<D> resource) {
+    return () -> traceAction(name, () -> {
+      var result = resource.getDynamics();
+      try {
+        assertion.accept(result);
+      } catch (Exception e) {
+        result = failure(e);
+      }
+      return result;
+    });
+  }
+
+  public static <D extends Dynamics<?, D>> CellResource<D> trace(String name, CellResource<D> resource) {
+    return traceFull(name, $ -> {}, resource);
+  }
+
+  public static <D extends Dynamics<?, D>> CellResource<D> trace(String name, Consumer<D> assertion, CellResource<D> resource) {
+    return traceExpiring(name, $ -> assertion.accept($.data()), resource);
+  }
+
+  public static <D extends Dynamics<?, D>> CellResource<D> traceExpiring(String name, Consumer<Expiring<D>> assertion, CellResource<D> resource) {
+    return traceFull(name, $ -> $.match(d -> {
+      assertion.accept(d);
+      return Unit.UNIT;
+    }, e -> {
+      throw new AssertionError("%s failed while computing".formatted(formatStack()), e);
+    }), resource);
+  }
+
+  public static <D extends Dynamics<?, D>> CellResource<D> traceFull(String name, Consumer<ErrorCatching<Expiring<D>>> assertion, CellResource<D> resource) {
+    return new CellResource<>() {
+      private final Resource<D> tracedResource = traceFull(name, assertion, (Resource<D>)resource);
+
+      @Override
+      public void emit(final Labelled<DynamicsEffect<D>> effect) {
+        resource.emit(effect);
+      }
+
+      @Override
+      public ErrorCatching<Expiring<D>> getDynamics() {
+        return tracedResource.getDynamics();
+      }
+
+      @Override
+      public void registerName(final String name) {
+        resource.registerName(name);
+      }
+    };
+  }
+
+  public static Condition trace(String name, Condition condition) {
+    return (positive, atEarliest, atLatest) ->
+        traceAction(name + " evaluate (%s, %s, %s)".formatted(positive, atEarliest, atLatest), () -> condition.nextSatisfied(positive, atEarliest, atLatest));
+  }
+
+  public static Supplier<Condition> trace(String name, Supplier<Condition> condition) {
+    // Trace calling the supplier separately from tracing the condition itself.
+    return () -> traceAction(name + " (generation)", () -> trace(name, condition.get()));
+  }
+
+  private static <T> T traceAction(String name, Supplier<T> action) {
+    activeTracePoints.push(name);
+    System.out.printf("TRACE: %s - %s start...%n", currentTime(), formatStack());
+    T result = action.get();
+    System.out.printf("TRACE: %s - %s: %s%n", currentTime(), formatStack(), result);
+    activeTracePoints.pop();
+    return result;
+  }
+
+  private static String formatStack() {
+    return String.join("->", activeTracePoints);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/Demo.java

@@ -0,0 +1,105 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling;
+
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.DiscreteResources;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.Polynomial;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAware;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+import gov.nasa.jpl.aerie.contrib.streamline.core.CellResource;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete.discrete;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.DiscreteEffects.set;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.DiscreteEffects.toggle;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.DiscreteEffects.using;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads.DiscreteResourceMonad.map;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.Polynomial.polynomial;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.PolynomialEffects.consume;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.PolynomialResources.asPolynomial;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.PolynomialResources.clamp;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.PolynomialResources.constant;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.PolynomialResources.integrate;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.PolynomialResources.lessThan;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.PolynomialResources.lessThan$;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.PolynomialResources.unitAware;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.Quantities.quantity;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.StandardUnits.*;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAwareOperations.simplify;
+import static gov.nasa.jpl.aerie.merlin.framework.ModelActions.*;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellResource.cellResource;
+
+public final class Demo {
+
+  // Unit-naive version of a model, to demonstrate some core concepts:
+
+  // Consumable, continuous:
+  CellResource<Polynomial> fuel_kg = cellResource(polynomial(20.0));
+  // Non-consumable, discrete:
+  CellResource<Discrete<Double>> power_w = cellResource(discrete(120.0));
+  // Atomic non-consumable:
+  CellResource<Discrete<Integer>> rwaControl = cellResource(discrete(1));
+  // Settable / enum state:
+  CellResource<Discrete<OnOff>> enumSwitch = cellResource(discrete(OnOff.ON));
+  // Toggle / flag:
+  CellResource<Discrete<Boolean>> boolSwitch = cellResource(discrete(true));
+
+  // Derived states:
+  Resource<Discrete<OnOff>> derivedEnumSwitch = map(boolSwitch, b -> b ? OnOff.ON : OnOff.OFF);
+  Resource<Polynomial> batterySOC_J = integrate(asPolynomial(power_w), 100);
+  Resource<Discrete<Double>> clampedPower_w = map(power_w, p -> p < 0 ? 0 : p);
+  Resource<Polynomial> clampedBatterySOC_J = clamp(batterySOC_J, constant(0), constant(100));
+  Resource<Discrete<Boolean>> lowPower = lessThan(batterySOC_J, 20);
+  Resource<Discrete<Boolean>> badness = map(
+      lowPower, enumSwitch,
+      (lowPower$, switch$) ->
+          lowPower$ && switch$ == OnOff.OFF);
+
+  {
+    using(power_w, 10, () -> {
+      using(rwaControl, () -> {
+        // Consume 5.4 kg of fuel over the next minute, linearly
+        consume(fuel_kg, 5.4, Duration.MINUTE);
+        // Separately, we could be doing things during that minute.
+        delay(Duration.MINUTE);
+      });
+      set(enumSwitch, OnOff.OFF);
+      toggle(boolSwitch);
+    });
+
+    set(boolSwitch, false);
+  }
+
+
+  // The exact same model again, but this time made unit-aware throughout.
+  // States without units have been re-used instead of being re-defined
+
+  // Consumable, continuous:
+  // CellResource<Polynomial> fuel_kg = cellResource(polynomial(20.0));
+  UnitAware<CellResource<Polynomial>> fuel = unitAware(
+      cellResource(polynomial(20.0)), KILOGRAM);
+  // Non-consumable, discrete:
+  UnitAware<CellResource<Discrete<Double>>> power = DiscreteResources.unitAware(
+      cellResource(discrete(120.0)), WATT);
+
+  UnitAware<Resource<Polynomial>> batterySOC = integrate(asPolynomial(simplify(power)), quantity(100, JOULE));
+  UnitAware<Resource<Discrete<Double>>> clampedPower = DiscreteResources.unitAware(map(power.value(WATT), p -> p < 0 ? 0 : p), WATT);
+  UnitAware<Resource<Discrete<Double>>> clampedPower_v2 = /* map(power, p -> lessThan(p, quantity(0, WATT)) ? quantity(0, WATT) : p) */
+      null;
+  UnitAware<Resource<Polynomial>> clampedBatterySOC = clamp(batterySOC, constant(quantity(0, JOULE)), constant(quantity(100, JOULE)));
+  Resource<Discrete<Boolean>> lowPower$ = lessThan$(batterySOC, quantity(20, JOULE));
+
+  {
+    using(power, quantity(10, WATT), () -> {
+      using(rwaControl, () -> {
+        // Consume 5.4 kg of fuel over the next minute, linearly
+        consume(fuel, quantity(5.4, KILOGRAM), Duration.MINUTE);
+        // Separately, we could be doing things during that minute.
+        delay(Duration.MINUTE);
+      });
+      set(enumSwitch, OnOff.OFF);
+      toggle(boolSwitch);
+    });
+  }
+
+  public enum OnOff { ON, OFF }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/Registrar.java

@@ -0,0 +1,121 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling;
+
+import gov.nasa.jpl.aerie.contrib.serialization.mappers.IntegerValueMapper;
+import gov.nasa.jpl.aerie.contrib.serialization.mappers.NullableValueMapper;
+import gov.nasa.jpl.aerie.contrib.serialization.mappers.StringValueMapper;
+import gov.nasa.jpl.aerie.contrib.streamline.core.CellResource;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Dynamics;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resources;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.linear.Linear;
+import gov.nasa.jpl.aerie.merlin.framework.ValueMapper;
+import gov.nasa.jpl.aerie.merlin.protocol.types.RealDynamics;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Unit;
+import org.apache.commons.lang3.exception.ExceptionUtils;
+
+import java.util.Collection;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.Map;
+import java.util.Set;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellResource.cellResource;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Reactions.wheneverDynamicsChange;
+import static gov.nasa.jpl.aerie.contrib.streamline.debugging.Tracing.trace;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads.DiscreteDynamicsMonad.effect;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads.DiscreteResourceMonad.map;
+import static java.util.stream.Collectors.joining;
+
+public class Registrar {
+  private final gov.nasa.jpl.aerie.merlin.framework.Registrar baseRegistrar;
+  private boolean trace = false;
+  private final CellResource<Discrete<Map<Throwable, Set<String>>>> errors;
+
+  public Registrar(final gov.nasa.jpl.aerie.merlin.framework.Registrar baseRegistrar) {
+    Resources.init();
+    this.baseRegistrar = baseRegistrar;
+    errors = cellResource(Discrete.discrete(Map.of()));
+    var errorString = map(errors, errors$ -> errors$.entrySet().stream().map(entry -> formatError(entry.getKey(), entry.getValue())).collect(joining("\n\n")));
+    discrete("errors", errorString, new StringValueMapper());
+    discrete("numberOfErrors", map(errors, Map::size), new IntegerValueMapper());
+  }
+
+  private static String formatError(Throwable e, Collection<String> affectedResources) {
+    return "Error affecting %s:%n%s".formatted(
+        String.join(", ", affectedResources),
+        formatException(e));
+  }
+
+  private static String formatException(Throwable e) {
+    return ExceptionUtils.stream(e)
+        .map(ExceptionUtils::getMessage)
+        .collect(joining("\nCaused by: "));
+  }
+
+  public void setTrace() {
+    trace = true;
+  }
+
+  public void clearTrace() {
+    trace = false;
+  }
+
+  public <Value> void discrete(final String name, final Resource<Discrete<Value>> resource, final ValueMapper<Value> mapper) {
+    resource.registerName(name);
+    var registeredResource = trace ? trace(name, resource) : resource;
+    baseRegistrar.discrete(
+        name,
+        () -> registeredResource.getDynamics().match(v -> v.data().extract(), e -> null),
+        new NullableValueMapper<>(mapper));
+    logErrors(name, registeredResource);
+  }
+
+  public void real(final String name, final Resource<Linear> resource) {
+    resource.registerName(name);
+    var registeredResource = trace ? trace(name, resource) : resource;
+    baseRegistrar.real(name, () -> registeredResource.getDynamics().match(
+        v -> RealDynamics.linear(v.data().extract(), v.data().rate()),
+        e -> RealDynamics.constant(0)));
+    logErrors(name, registeredResource);
+  }
+
+  private <D extends Dynamics<?, D>> void logErrors(String name, Resource<D> resource) {
+     wheneverDynamicsChange(resource, ec -> ec.match($ -> null, e -> logError(name, e)));
+  }
+
+  // TODO: Consider pulling in a Guava MultiMap instead of doing this by hand below
+  private Unit logError(String resourceName, Throwable e) {
+    errors.emit(effect(s -> {
+      var s$ = new HashMap<>(s);
+      s$.compute(e, (e$, affectedResources) -> {
+        if (affectedResources == null) {
+          return Set.of(resourceName);
+        } else {
+          var affectedResources$ = new HashSet<>(affectedResources);
+          affectedResources$.add(resourceName);
+          return affectedResources$;
+        }
+      });
+      return s$;
+    }));
+    return Unit.UNIT;
+  }
+
+  private Unit removeError(String resourceName, Throwable e) {
+    errors.emit(effect(s -> {
+      var s$ = new HashMap<>(s);
+      s$.compute(e, (e$, affectedResources) -> {
+        if (affectedResources == null) {
+          return null;
+        } else {
+          var affectedResources$ = new HashSet<>(affectedResources);
+          affectedResources$.remove(resourceName);
+          return affectedResources$.isEmpty() ? null : affectedResources$;
+        }
+      });
+      return s$;
+    }));
+    return Unit.UNIT;
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/clocks/Clock.java

@@ -0,0 +1,15 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.clocks;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Dynamics;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+public record Clock(Duration extract) implements Dynamics<Duration, Clock> {
+  @Override
+  public Clock step(Duration t) {
+    return clock(extract().plus(t));
+  }
+
+  public static Clock clock(Duration startingTime) {
+    return new Clock(startingTime);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/clocks/ClockResources.java

@@ -0,0 +1,49 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.clocks;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.ExpiringToResourceMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads.DiscreteResourceMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellRefV2.allocate;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellResource.cellResource;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching.success;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Expiring.*;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Resources.signalling;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.monads.ResourceMonad.bind;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete.discrete;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.DiscreteResources.not;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.EPSILON;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.ZERO;
+
+public final class ClockResources {
+  private ClockResources() {}
+
+  public static Resource<Clock> clock() {
+    return cellResource(Clock.clock(ZERO));
+  }
+
+  public static Resource<Discrete<Boolean>> lessThan(Resource<Clock> clock, Duration threshold) {
+    return bind(clock, (Clock c) -> {
+      final Duration crossoverTime = threshold.minus(c.extract());
+      return ExpiringToResourceMonad.unit(
+          crossoverTime.isPositive()
+              ? expiring(discrete(true), crossoverTime)
+              : neverExpiring(discrete(false)));
+    });
+  }
+
+  public static Resource<Discrete<Boolean>> lessThanOrEquals(Resource<Clock> clock, Duration threshold) {
+    // Since Duration is an integral type, implement strictness through EPSILON stepping
+    return lessThan(clock, threshold.plus(EPSILON));
+  }
+
+  public static Resource<Discrete<Boolean>> greaterThan(Resource<Clock> clock, Duration threshold) {
+    return not(lessThanOrEquals(clock, threshold));
+  }
+
+  public static Resource<Discrete<Boolean>> greaterThanOrEquals(Resource<Clock> clock, Duration threshold) {
+    return not(lessThan(clock, threshold));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/discrete/Discrete.java

@@ -0,0 +1,15 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Dynamics;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+public record Discrete<V>(V extract) implements Dynamics<V, Discrete<V>> {
+  @Override
+  public Discrete<V> step(Duration t) {
+    return this;
+  }
+
+  public static <V> Discrete<V> discrete(V value) {
+    return new Discrete<>(value);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/discrete/DiscreteEffects.java

@@ -0,0 +1,100 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.CellResource;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAware;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads.DiscreteDynamicsMonad.effect;
+
+public final class DiscreteEffects {
+  private DiscreteEffects() {}
+
+  // More convenient overload of "set" when using discrete dynamics
+
+  public static <A> void set(CellResource<Discrete<A>> resource, A newValue) {
+    resource.emit("Set " + newValue, effect(x -> newValue));
+  }
+
+  // Flag/Switch style operations
+
+  public static void set(CellResource<Discrete<Boolean>> resource) {
+    set(resource, true);
+  }
+
+  public static void unset(CellResource<Discrete<Boolean>> resource) {
+    set(resource, false);
+  }
+
+  public static void toggle(CellResource<Discrete<Boolean>> resource) {
+    resource.emit(effect(x -> !x));
+  }
+
+  // Counter style operations
+
+  public static void increment(CellResource<Discrete<Integer>> resource) {
+    increment(resource, 1);
+  }
+
+  public static void increment(CellResource<Discrete<Integer>> resource, int amount) {
+    resource.emit(effect(x -> x + amount));
+  }
+
+  public static void decrement(CellResource<Discrete<Integer>> resource) {
+    decrement(resource, 1);
+  }
+
+  public static void decrement(CellResource<Discrete<Integer>> resource, int amount) {
+    resource.emit(effect(x -> x - amount));
+  }
+
+  // Consumable style operations
+
+  public static void consume(CellResource<Discrete<Double>> resource, double amount) {
+    resource.emit(effect(x -> x - amount));
+  }
+
+  public static void restore(CellResource<Discrete<Double>> resource, double amount) {
+    resource.emit(effect(x -> x + amount));
+  }
+
+  // Non-consumable style operations
+
+  public static void using(CellResource<Discrete<Double>> resource, double amount, Runnable action) {
+    consume(resource, amount);
+    action.run();
+    restore(resource, amount);
+  }
+
+  // Atomic style operations
+
+  public static void using(CellResource<Discrete<Integer>> resource, Runnable action) {
+    decrement(resource);
+    action.run();
+    increment(resource);
+  }
+
+  // Unit-aware effects:
+
+  // More convenient overload of "set" when using discrete dynamics
+
+  public static <A> void set(UnitAware<CellResource<Discrete<A>>> resource, UnitAware<A> newValue) {
+    set(resource.value(), newValue.value(resource.unit()));
+  }
+
+  // Consumable style operations
+
+  public static void consume(UnitAware<CellResource<Discrete<Double>>> resource, UnitAware<Double> amount) {
+    consume(resource.value(), amount.value(resource.unit()));
+  }
+
+  public static void restore(UnitAware<CellResource<Discrete<Double>>> resource, UnitAware<Double> amount) {
+    restore(resource.value(), amount.value(resource.unit()));
+  }
+
+  // Non-consumable style operations
+
+  public static void using(UnitAware<CellResource<Discrete<Double>>> resource, UnitAware<Double> amount, Runnable action) {
+    consume(resource, amount);
+    action.run();
+    restore(resource, amount);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/discrete/DiscreteResources.java

@@ -0,0 +1,87 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads.DiscreteMonad;
+import gov.nasa.jpl.aerie.merlin.framework.Condition;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+import gov.nasa.jpl.aerie.contrib.streamline.core.CellResource;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.Unit;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAware;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAwareResources;
+
+import java.util.Arrays;
+import java.util.Optional;
+import java.util.function.BiPredicate;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellResource.cellResource;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Reactions.whenever;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete.discrete;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads.DiscreteResourceMonad.*;
+
+public final class DiscreteResources {
+  private DiscreteResources() {}
+
+  public static Condition when(Resource<Discrete<Boolean>> resource) {
+    return (positive, atEarliest, atLatest) ->
+        resource.getDynamics().match(
+            dynamics -> Optional.of(atEarliest).filter($ -> dynamics.data().extract() == positive),
+            error -> Optional.empty());
+  }
+
+  public static <V> Resource<Discrete<V>> cache(Resource<Discrete<V>> resource, BiPredicate<V, V> updatePredicate) {
+    final var cell = cellResource(resource.getDynamics());
+    BiPredicate<ErrorCatching<Expiring<Discrete<V>>>, ErrorCatching<Expiring<Discrete<V>>>> liftedUpdatePredicate = (eCurrent, eNew) ->
+        eCurrent.match(
+            current -> eNew.match(
+                value -> updatePredicate.test(current.data().extract(), value.data().extract()),
+                newException -> true),
+            currentException -> eNew.match(
+                value -> true,
+                newException -> !currentException.equals(newException)));
+    whenever(() -> {
+      var currentDynamics = resource.getDynamics();
+      return when(() -> DynamicsMonad.unit(discrete(liftedUpdatePredicate.test(currentDynamics, resource.getDynamics()))));
+    }, () -> {
+      final var newDynamics = resource.getDynamics();
+      cell.emit($ -> newDynamics);
+    });
+    return cell;
+  }
+
+  public static UnitAware<Resource<Discrete<Double>>> unitAware(Resource<Discrete<Double>> resource, Unit unit) {
+    return UnitAwareResources.unitAware(resource, unit, DiscreteResources::discreteScaling);
+  }
+
+  public static UnitAware<CellResource<Discrete<Double>>> unitAware(CellResource<Discrete<Double>> resource, Unit unit) {
+    return UnitAwareResources.unitAware(resource, unit, DiscreteResources::discreteScaling);
+  }
+
+  private static Discrete<Double> discreteScaling(Discrete<Double> d, Double scale) {
+    return DiscreteMonad.map(d, $ -> $ * scale);
+  }
+
+  // Boolean logic
+
+  @SafeVarargs
+  public static Resource<Discrete<Boolean>> and(Resource<Discrete<Boolean>>... operands) {
+    return Arrays.stream(operands).reduce(unit(true), lift(Boolean::logicalAnd)::apply);
+  }
+
+  @SafeVarargs
+  public static Resource<Discrete<Boolean>> or(Resource<Discrete<Boolean>>... operands) {
+    return Arrays.stream(operands).reduce(unit(false), lift(Boolean::logicalOr)::apply);
+  }
+
+  public static Resource<Discrete<Boolean>> not(Resource<Discrete<Boolean>> operand) {
+    return map(operand, $ -> !$);
+  }
+
+  public static Resource<Discrete<Boolean>> assertThat(String description, Resource<Discrete<Boolean>> assertion) {
+    return map(assertion, a -> {
+      if (a) return true;
+      throw new AssertionError(description);
+    });
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/discrete/monads/DiscreteDynamicsMonad.java

@@ -0,0 +1,37 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.DynamicsEffect;
+import gov.nasa.jpl.aerie.contrib.streamline.core.ErrorCatching;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+
+import java.util.function.Function;
+
+public final class DiscreteDynamicsMonad {
+  private DiscreteDynamicsMonad() {}
+
+  public static <A> ErrorCatching<Expiring<Discrete<A>>> unit(A a) {
+    return DiscreteMonadTransformer.unit(DynamicsMonad::unit, a);
+  }
+
+  public static <A, B> ErrorCatching<Expiring<Discrete<B>>> bind(ErrorCatching<Expiring<Discrete<A>>> a, Function<A, ErrorCatching<Expiring<Discrete<B>>>> f) {
+    return DiscreteMonadTransformer.bind(DynamicsMonad::bind, a, f);
+  }
+
+  // Convenient methods defined in terms of bind and unit:
+
+  public static <A, B> ErrorCatching<Expiring<Discrete<B>>> map(ErrorCatching<Expiring<Discrete<A>>> a, Function<A, B> f) {
+    return bind(a, f.andThen(DiscreteDynamicsMonad::unit));
+  }
+
+  public static <A, B> Function<ErrorCatching<Expiring<Discrete<A>>>, ErrorCatching<Expiring<Discrete<B>>>> lift(Function<A, B> f) {
+    return a -> map(a, f);
+  }
+
+  // Not monadic, strictly speaking, but useful nonetheless.
+
+  public static <A> DynamicsEffect<Discrete<A>> effect(Function<A, A> f) {
+      return DynamicsMonad.effect(DiscreteMonad.lift(f));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/discrete/monads/DiscreteExpiringMonad.java

@@ -0,0 +1,25 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.ExpiringMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+
+import java.util.function.Function;
+
+public final class DiscreteExpiringMonad {
+  private DiscreteExpiringMonad() {}
+
+  public static <A> Expiring<Discrete<A>> unit(A a) {
+    return DiscreteMonadTransformer.unit(ExpiringMonad::unit, a);
+  }
+
+  public static <A, B> Expiring<Discrete<B>> bind(Expiring<Discrete<A>> a, Function<A, Expiring<Discrete<B>>> f) {
+    return DiscreteMonadTransformer.bind(ExpiringMonad::bind, a, f);
+  }
+
+  // Convenient methods defined in terms of bind and unit:
+
+  public static <A, B> Expiring<Discrete<B>> map(Expiring<Discrete<A>> a, Function<A, B> f) {
+    return bind(a, f.andThen(DiscreteExpiringMonad::unit));
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/discrete/monads/DiscreteMonad.java

@@ -0,0 +1,31 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.IdentityMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+
+import java.util.function.Function;
+
+/**
+ * {@link Discrete} monad
+ */
+public final class DiscreteMonad {
+  private DiscreteMonad() {}
+
+  public static <A> Discrete<A> unit(A a) {
+    return DiscreteMonadTransformer.unit(IdentityMonad::unit, a);
+  }
+
+  public static <A, B> Discrete<B> bind(Discrete<A> a, Function<A, Discrete<B>> f) {
+    return DiscreteMonadTransformer.bind(IdentityMonad::bind, a, f);
+  }
+
+  // Convenient methods defined in terms of bind and unit:
+
+  public static <A, B> Discrete<B> map(Discrete<A> a, Function<A, B> f) {
+    return bind(a, f.andThen(DiscreteMonad::unit));
+  }
+
+  public static <A, B> Function<Discrete<A>, Discrete<B>> lift(Function<A, B> f) {
+    return a -> map(a, f);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/discrete/monads/DiscreteMonadTransformer.java

@@ -0,0 +1,25 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+
+import java.util.function.BiFunction;
+import java.util.function.Function;
+
+/**
+ * Monad transformer for {@link Discrete}, M A -> M ({@link Discrete}<A>)
+ */
+public final class DiscreteMonadTransformer {
+  private DiscreteMonadTransformer() {}
+
+  public static <A, MDA> MDA unit(Function<Discrete<A>, MDA> mUnit, A a) {
+    return mUnit.apply(Discrete.discrete(a));
+  }
+
+  public static <A, MDA, MDB> MDB bind(BiFunction<MDA, Function<Discrete<A>, MDB>, MDB> mBind, MDA m, Function<A, MDB> f) {
+    return mBind.apply(m, d -> f.apply(d.extract()));
+  }
+
+  public static <A, MA, MDA> MDA lift(BiFunction<MA, Function<A, Discrete<A>>, MDA> mBind, MA m) {
+    return mBind.apply(m, Discrete::discrete);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/discrete/monads/DiscreteResourceMonad.java

@@ -0,0 +1,48 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resources;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.ResourceMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAware;
+import org.apache.commons.lang3.function.TriFunction;
+
+import java.util.function.BiFunction;
+import java.util.function.Function;
+
+public final class DiscreteResourceMonad {
+  private DiscreteResourceMonad() {}
+
+  public static <A> Resource<Discrete<A>> unit(A a) {
+    return DiscreteMonadTransformer.unit(ResourceMonad::unit, a);
+  }
+
+  public static <A, B> Resource<Discrete<B>> bind(Resource<Discrete<A>> a, Function<A, Resource<Discrete<B>>> f) {
+    return DiscreteMonadTransformer.bind(ResourceMonad::bind, a, f);
+  }
+
+  // Convenient methods defined in terms of bind and unit:
+
+  public static <A, B> Resource<Discrete<B>> map(Resource<Discrete<A>> a, Function<A, B> f) {
+    return bind(a, f.andThen(DiscreteResourceMonad::unit));
+  }
+
+  // Map functions with higher arities are defined for discrete resources,
+  // because deriving discrete values with many sources is fairly common.
+
+  public static <A, B, C> Resource<Discrete<C>> map(Resource<Discrete<A>> a, Resource<Discrete<B>> b, BiFunction<A, B, C> f) {
+    return bind(a, a$ -> map(b, b$ -> f.apply(a$, b$)));
+  }
+
+  public static <A, B, C, D> Resource<Discrete<D>> map(Resource<Discrete<A>> a, Resource<Discrete<B>> b, Resource<Discrete<C>> c, TriFunction<A, B, C, D> f) {
+    return bind(a, a$ -> map(b, c, (b$, c$) -> f.apply(a$, b$, c$)));
+  }
+
+  public static <A, B> Function<Resource<Discrete<A>>, Resource<Discrete<B>>> lift(Function<A, B> f) {
+    return a -> map(a, f);
+  }
+
+  public static <A, B, C> BiFunction<Resource<Discrete<A>>, Resource<Discrete<B>>, Resource<Discrete<C>>> lift(BiFunction<A, B, C> f) {
+    return (a, b) -> map(a, b, f);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/linear/Linear.java

@@ -0,0 +1,20 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.linear;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Dynamics;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+import java.util.Objects;
+
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.SECOND;
+
+// TODO: Implement better support for going to/from Linear
+public record Linear(Double extract, Double rate) implements Dynamics<Double, Linear> {
+  @Override
+  public Linear step(Duration t) {
+    return linear(extract() + t.ratioOver(SECOND) * rate(), rate());
+  }
+
+  public static Linear linear(double value, double rate) {
+    return new Linear(value, rate);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/Polynomial.java

@@ -0,0 +1,259 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Dynamics;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Expiring;
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.ExpiringMonad;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+import org.apache.commons.math3.analysis.solvers.LaguerreSolver;
+import org.apache.commons.math3.complex.Complex;
+
+import java.util.Arrays;
+import java.util.function.DoublePredicate;
+import java.util.function.Predicate;
+import java.util.stream.IntStream;
+import java.util.stream.Stream;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Expiring.expiring;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Expiry.NEVER;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Expiry.expiry;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.EPSILON;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.SECOND;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete.discrete;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.ZERO;
+import static org.apache.commons.math3.analysis.polynomials.PolynomialsUtils.shift;
+
+public record Polynomial(double[] coefficients) implements Dynamics<Double, Polynomial> {
+  /**
+   * Maximum imaginary component allowed in a root to be considered "real" when performing root-finding.
+   * Should be a very small number to avoid spurious roots.
+   */
+  private static final double ROOT_FINDING_IMAGINARY_COMPONENT_TOLERANCE = 1e-12;
+
+  /**
+   * Maximum number of time steps to search in either direction around near-roots
+   * to find the corresponding discretized transition point.
+   */
+  private static final int MAX_RANGE_FOR_ROOT_SEARCH = 2;
+
+  public static Polynomial polynomial(double... coefficients) {
+    int n = coefficients.length;
+    if (n == 0) {
+      return new Polynomial(new double[] { 0.0 });
+    }
+    while (n > 1 && coefficients[n - 1] == 0) --n;
+    return new Polynomial(Arrays.copyOf(coefficients, n));
+  }
+
+  @Override
+  public Double extract() {
+    return coefficients()[0];
+  }
+
+  @Override
+  public Polynomial step(Duration t) {
+    return t.isEqualTo(ZERO) ? this : polynomial(shift(coefficients(), t.ratioOver(SECOND)));
+  }
+
+  public int degree() {
+    return coefficients().length - 1;
+  }
+
+  public boolean isConstant() {
+    return degree() == 0;
+  }
+
+  public boolean isNonFinite() {
+    return Arrays.stream(coefficients()).anyMatch(c -> !Double.isFinite(c));
+  }
+
+  public Polynomial add(Polynomial other) {
+    final double[] coefficients = coefficients();
+    final double[] otherCoefficients = other.coefficients();
+    final int minLength = Math.min(coefficients.length, otherCoefficients.length);
+    final int maxLength = Math.max(coefficients.length, otherCoefficients.length);
+    final double[] newCoefficients = new double[maxLength];
+    for (int i = 0; i < minLength; ++i) {
+      newCoefficients[i] = coefficients[i] + otherCoefficients[i];
+    }
+    if (coefficients.length > minLength)
+      System.arraycopy(coefficients, minLength, newCoefficients, minLength, coefficients.length - minLength);
+    if (otherCoefficients.length > minLength)
+      System.arraycopy(
+          otherCoefficients, minLength, newCoefficients, minLength, otherCoefficients.length - minLength);
+    return polynomial(newCoefficients);
+  }
+
+  public Polynomial subtract(Polynomial other) {
+    return add(other.multiply(polynomial(-1)));
+  }
+
+  public Polynomial multiply(Polynomial other) {
+    final double[] coefficients = coefficients();
+    final double[] otherCoefficients = other.coefficients();
+    // Length = degree + 1, so
+    // new length = 1 + new degree
+    //   = 1 + (degree + other.degree)
+    //   = 1 + (length - 1 + other.length - 1)
+    //   = length + other.length - 1
+    final double[] newCoefficients = new double[coefficients.length + otherCoefficients.length - 1];
+    for (int exponent = 0; exponent < newCoefficients.length; ++exponent) {
+      newCoefficients[exponent] = 0.0;
+      // 0 <= k < length and 0 <= exponent - k < other.length
+      // implies k >= 0, k > exponent - other.length,
+      // k < length, and k <= exponent
+      for (int k = Math.max(0, exponent - otherCoefficients.length + 1);
+           k < Math.min(coefficients.length, exponent + 1);
+           ++k) {
+        newCoefficients[exponent] += coefficients[k] * otherCoefficients[exponent - k];
+      }
+    }
+    return polynomial(newCoefficients);
+  }
+
+  public Polynomial divide(double scalar) {
+    final double[] coefficients = coefficients();
+    final double[] newCoefficients = new double[coefficients.length];
+    for (int i = 0; i < coefficients.length; ++i) {
+      newCoefficients[i] = coefficients[i] / scalar;
+    }
+    return polynomial(newCoefficients);
+  }
+
+  public Polynomial integral(double startingValue) {
+    final double[] coefficients = coefficients();
+    final double[] newCoefficients = new double[coefficients.length + 1];
+    newCoefficients[0] = startingValue;
+    for (int i = 0; i < coefficients.length; ++i) {
+      newCoefficients[i + 1] = coefficients[i] / (i + 1);
+    }
+    return polynomial(newCoefficients);
+  }
+
+  public Polynomial derivative() {
+    final double[] coefficients = coefficients();
+    final double[] newCoefficients = new double[coefficients.length - 1];
+    for (int i = 1; i < coefficients.length; ++i) {
+      newCoefficients[i - 1] = coefficients[i] * i;
+    }
+    return polynomial(newCoefficients);
+  }
+
+  public double evaluate(Duration t) {
+    return evaluate(t.ratioOver(SECOND));
+  }
+
+  public double evaluate(double x) {
+    // Horner's method of polynomial evaluation:
+    // Transforms a_0 + a_1 x + a_2 x^2 + ... + a_n x^n
+    // into a_0 + x (a_1 + x (a_2 + ... x ( a_n ) ... ))
+    // Which can be done with one addition and one multiplication per coefficient,
+    // as opposed to the traditional method, which takes one addition and multiple multiplications.
+    final double[] coefficients = coefficients();
+    double accumulator = coefficients[coefficients.length - 1];
+    for (int i = coefficients.length - 2; i >= 0; --i) {
+      accumulator *= x;
+      accumulator += coefficients[i];
+    }
+    return accumulator;
+  }
+
+  private Expiring<Discrete<Boolean>> compare(DoublePredicate predicate, double threshold) {
+    return find(t -> predicate.test(evaluate(t)), threshold);
+  }
+
+  private Expiring<Discrete<Boolean>> find(Predicate<Duration> timePredicate, double target) {
+    final boolean currentValue = timePredicate.test(ZERO);
+    final var expiry = this.isConstant() || this.isNonFinite() ? NEVER : expiry(findFuturePreImage(target)
+        .flatMap(t -> IntStream.rangeClosed(-MAX_RANGE_FOR_ROOT_SEARCH, MAX_RANGE_FOR_ROOT_SEARCH)
+            .mapToObj(i -> t.plus(EPSILON.times(i))))
+        .filter(t -> (timePredicate.test(t) ^ currentValue) && t.isPositive())
+        .findFirst());
+    return expiring(discrete(currentValue), expiry);
+  }
+
+  public Expiring<Discrete<Boolean>> greaterThan(double threshold) {
+    return compare(x -> x > threshold, threshold);
+  }
+
+  public Expiring<Discrete<Boolean>> greaterThanOrEquals(double threshold) {
+    return compare(x -> x >= threshold, threshold);
+  }
+
+  public Expiring<Discrete<Boolean>> lessThan(double threshold) {
+    return compare(x -> x < threshold, threshold);
+  }
+
+  public Expiring<Discrete<Boolean>> lessThanOrEquals(double threshold) {
+    return compare(x -> x <= threshold, threshold);
+  }
+
+  private boolean dominates$(Polynomial other) {
+    for (int i = 0; i <= Math.max(this.degree(), other.degree()); ++i) {
+      if (this.getCoefficient(i) > other.getCoefficient(i)) return true;
+      if (this.getCoefficient(i) < other.getCoefficient(i)) return false;
+    }
+    // Equal, so either answer is correct
+    return true;
+  }
+
+  private Expiring<Discrete<Boolean>> dominates(Polynomial other) {
+    return this.subtract(other).find(t -> this.step(t).dominates$(other.step(t)), 0);
+  }
+
+  public Expiring<Polynomial> min(Polynomial other) {
+    return ExpiringMonad.map(this.dominates(other), d -> d.extract() ? other : this);
+  }
+
+  public Expiring<Polynomial> max(Polynomial other) {
+    return ExpiringMonad.map(this.dominates(other), d -> d.extract() ? this : other);
+  }
+
+  /**
+   * Finds all occasions in the future when this function will reach the target value.
+   */
+  private Stream<Duration> findFuturePreImage(double target) {
+    // add a check for an infinite target (i.e. unbounded above/below) or a poorly behaved polynomial
+    if (!Double.isFinite(target) || this.isNonFinite()) {
+      return Stream.empty();
+    }
+
+    final double[] shiftedCoefficients = add(polynomial(-target)).coefficients();
+    final Complex[] solutions = new LaguerreSolver().solveAllComplex(shiftedCoefficients, 0);
+    return Arrays.stream(solutions)
+                 .filter(solution -> Math.abs(solution.getImaginary()) < ROOT_FINDING_IMAGINARY_COMPONENT_TOLERANCE)
+                 .map(Complex::getReal)
+                 .filter(t -> t >= 0)
+                 .sorted()
+                 .map(t -> Duration.roundNearest(t, SECOND));
+  }
+
+  /**
+   * Get the nth coefficient.
+   * @param n the n.
+   * @return the nth coefficient
+   */
+  public double getCoefficient(int n) {
+    return n >= coefficients().length ? 0.0 : coefficients()[n];
+  }
+
+  @Override
+  public boolean equals(final Object o) {
+    if (this == o) return true;
+    if (o == null || getClass() != o.getClass()) return false;
+    Polynomial that = (Polynomial) o;
+    return Arrays.equals(coefficients, that.coefficients);
+  }
+
+  @Override
+  public int hashCode() {
+    return Arrays.hashCode(coefficients);
+  }
+
+  @Override
+  public String toString() {
+    return "Polynomial{" +
+           "coefficients=" + Arrays.toString(coefficients) +
+           '}';
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/PolynomialEffects.java

@@ -0,0 +1,104 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resources;
+import gov.nasa.jpl.aerie.contrib.streamline.core.CellResource;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAware;
+
+import static gov.nasa.jpl.aerie.merlin.framework.ModelActions.delay;
+import static gov.nasa.jpl.aerie.merlin.framework.ModelActions.replaying;
+import static gov.nasa.jpl.aerie.merlin.framework.ModelActions.spawn;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad.effect;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.SECOND;
+
+public final class PolynomialEffects {
+  private PolynomialEffects() {}
+
+  // Consumable style operations
+
+  public static void consuming(CellResource<Polynomial> resource, Polynomial profile, Runnable action) {
+    resource.emit(effect($ -> $.subtract(profile)));
+    final Duration start = Resources.currentTime();
+    action.run();
+    final Duration elapsedTime = Resources.currentTime().minus(start);
+    // Nullify ongoing effects by adding a profile with the same behavior,
+    // but with an initial value of 0
+    final Polynomial steppedProfile = profile.step(elapsedTime);
+    final Polynomial counteractingProfile = steppedProfile.subtract(Polynomial.polynomial(steppedProfile.extract()));
+    resource.emit(effect($ -> $.add(counteractingProfile)));
+  }
+
+  public static void consume(CellResource<Polynomial> resource, double amount) {
+    resource.emit(effect($ -> $.subtract(Polynomial.polynomial(amount))));
+  }
+
+  public static void consume(CellResource<Polynomial> resource, double amount, Duration time) {
+    final Polynomial profile = Polynomial.polynomial(0, amount / time.ratioOver(SECOND));
+    spawn(replaying(() -> consuming(resource, profile, () -> delay(time))));
+  }
+
+  public static void restoring(CellResource<Polynomial> resource, Polynomial profile, Runnable action) {
+    consuming(resource, profile.multiply(Polynomial.polynomial(-1)), action);
+  }
+
+  public static void restore(CellResource<Polynomial> resource, double amount) {
+    consume(resource, -amount);
+  }
+
+  public static void restore(CellResource<Polynomial> resource, double amount, Duration time) {
+    consume(resource, -amount, time);
+  }
+
+  // Non-consumable style operations
+
+  public static void using(CellResource<Polynomial> resource, Polynomial profile, Runnable action) {
+    resource.emit(effect($ -> $.subtract(profile)));
+    final Duration start = Resources.currentTime();
+    action.run();
+    final Duration elapsedTime = Resources.currentTime().minus(start);
+    // Reset by adding a counteracting profile
+    final Polynomial counteractingProfile = profile.step(elapsedTime);
+    resource.emit(effect($ -> $.add(counteractingProfile)));
+  }
+
+  public static void using(CellResource<Polynomial> resource, double amount, Runnable action) {
+    using(resource, Polynomial.polynomial(amount), action);
+  }
+
+  // Consumable style operations
+
+  public static void consuming(UnitAware<CellResource<Polynomial>> resource, UnitAware<Polynomial> profile, Runnable action) {
+    consuming(resource.value(), profile.value(resource.unit()), action);
+  }
+
+  public static void consume(UnitAware<CellResource<Polynomial>> resource, UnitAware<Double> amount) {
+    consume(resource.value(), amount.value(resource.unit()));
+  }
+
+  public static void consume(UnitAware<CellResource<Polynomial>> resource, UnitAware<Double> amount, Duration time) {
+    consume(resource.value(), amount.value(resource.unit()), time);
+  }
+
+  public static void restoring(UnitAware<CellResource<Polynomial>> resource, UnitAware<Polynomial> profile, Runnable action) {
+    restoring(resource.value(), profile.value(resource.unit()), action);
+  }
+
+  public static void restore(UnitAware<CellResource<Polynomial>> resource, UnitAware<Double> amount) {
+    restore(resource.value(), amount.value(resource.unit()));
+  }
+
+  public static void restore(UnitAware<CellResource<Polynomial>> resource, UnitAware<Double> amount, Duration time) {
+    restore(resource.value(), amount.value(resource.unit()), time);
+  }
+
+  // Non-consumable style operations
+
+  // Ugly $ suffix to avoid overload conflict because of erasure.
+  public static void using$(UnitAware<CellResource<Polynomial>> resource, UnitAware<Polynomial> profile, Runnable action) {
+    using(resource.value(), profile.value(resource.unit()), action);
+  }
+
+  public static void using(UnitAware<CellResource<Polynomial>> resource, UnitAware<Double> amount, Runnable action) {
+    using(resource.value(), amount.value(resource.unit()), action);
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/PolynomialResources.java

@@ -0,0 +1,232 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial;
+
+import gov.nasa.jpl.aerie.contrib.streamline.core.CellResource;
+import gov.nasa.jpl.aerie.contrib.streamline.core.Resource;
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.ExpiringToResourceMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.core.monads.ResourceMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.Discrete;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.discrete.monads.DiscreteResourceMonad;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.Unit;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAware;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAwareOperations;
+import gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAwareResources;
+import gov.nasa.jpl.aerie.merlin.protocol.types.Duration;
+
+import java.util.Arrays;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.core.CellResource.cellResource;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Reactions.whenever;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Reactions.wheneverDynamicsChange;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Resources.currentValue;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.Resources.shift;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad.bindEffect;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad.effect;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad.map;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.monads.DynamicsMonad.unit;
+import static gov.nasa.jpl.aerie.contrib.streamline.core.monads.ResourceMonad.bind;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.polynomial.Polynomial.polynomial;
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAwareResources.extend;
+import static gov.nasa.jpl.aerie.merlin.protocol.types.Duration.SECOND;
+
+public final class PolynomialResources {
+  private PolynomialResources() {}
+
+  public static Resource<Polynomial> constant(double value) {
+    var dynamics = unit(polynomial(value));
+    return () -> dynamics;
+  }
+
+  public static UnitAware<Resource<Polynomial>> constant(UnitAware<Double> quantity) {
+    return unitAware(constant(quantity.value()), quantity.unit());
+  }
+
+  public static Resource<Polynomial> asPolynomial(Resource<Discrete<Double>> discrete) {
+    return ResourceMonad.map(discrete, d -> polynomial(d.extract()));
+  }
+
+  public static UnitAware<Resource<Polynomial>> asPolynomial(UnitAware<Resource<Discrete<Double>>> discrete) {
+    return unitAware(asPolynomial(discrete.value()), discrete.unit());
+  }
+
+  @SafeVarargs
+  public static Resource<Polynomial> add(Resource<Polynomial>... summands) {
+    return Arrays.stream(summands)
+        .reduce(constant(0), (p, q) -> ResourceMonad.map(p, q, Polynomial::add));
+  }
+
+  public static Resource<Polynomial> subtract(Resource<Polynomial> p, Resource<Polynomial> q) {
+    return ResourceMonad.map(p, q, Polynomial::subtract);
+  }
+
+  public static Resource<Polynomial> negate(Resource<Polynomial> p) {
+    return multiply(constant(-1), p);
+  }
+
+  @SafeVarargs
+  public static Resource<Polynomial> multiply(Resource<Polynomial>... factors) {
+    return Arrays.stream(factors)
+                 .reduce(constant(1), (p, q) -> ResourceMonad.map(p, q, Polynomial::multiply));
+  }
+
+  public static Resource<Polynomial> divide(Resource<Polynomial> p, Resource<Discrete<Double>> q) {
+    return ResourceMonad.map(p, q, (p$, q$) -> p$.divide(q$.extract()));
+  }
+
+  public static Resource<Polynomial> integrate(Resource<Polynomial> integrand, double startingValue) {
+    var cell = cellResource(map(integrand.getDynamics(), (Polynomial $) -> $.integral(startingValue)));
+    // Use integrand's expiry but not integral's, since we're refreshing the integral
+    wheneverDynamicsChange(integrand, integrandDynamics ->
+        cell.emit(bindEffect(integral -> DynamicsMonad.map(integrandDynamics, integrand$ ->
+            integrand$.integral(integral.extract())))));
+    return cell;
+  }
+
+  public static Resource<Polynomial> differentiate(Resource<Polynomial> p) {
+    return ResourceMonad.map(p, Polynomial::derivative);
+  }
+
+  public static Resource<Polynomial> movingAverage(Resource<Polynomial> p, Duration interval) {
+    var pIntegral = integrate(p, 0);
+    var shiftedIntegral = shift(pIntegral, interval, polynomial(0));
+    return divide(subtract(pIntegral, shiftedIntegral), DiscreteResourceMonad.unit(interval.ratioOver(SECOND)));
+  }
+
+  public static Resource<Discrete<Boolean>> greaterThan(Resource<Polynomial> p, double threshold) {
+    return bind(p, p$ -> ExpiringToResourceMonad.unit(p$.greaterThan(threshold)));
+  }
+
+  public static Resource<Discrete<Boolean>> greaterThanOrEquals(Resource<Polynomial> p, double threshold) {
+    return bind(p, p$ -> ExpiringToResourceMonad.unit(p$.greaterThanOrEquals(threshold)));
+  }
+
+  public static Resource<Discrete<Boolean>> lessThan(Resource<Polynomial> p, double threshold) {
+    return bind(p, p$ -> ExpiringToResourceMonad.unit(p$.lessThan(threshold)));
+  }
+
+  public static Resource<Discrete<Boolean>> lessThanOrEquals(Resource<Polynomial> p, double threshold) {
+    return bind(p, p$ -> ExpiringToResourceMonad.unit(p$.lessThanOrEquals(threshold)));
+  }
+
+  public static Resource<Discrete<Boolean>> greaterThan(Resource<Polynomial> p, Resource<Polynomial> q) {
+    return greaterThan(subtract(p, q), 0);
+  }
+
+  public static Resource<Discrete<Boolean>> greaterThanOrEquals(Resource<Polynomial> p, Resource<Polynomial> q) {
+    return greaterThanOrEquals(subtract(p, q), 0);
+  }
+
+  public static Resource<Discrete<Boolean>> lessThan(Resource<Polynomial> p, Resource<Polynomial> q) {
+    return lessThan(subtract(p, q), 0);
+  }
+
+  public static Resource<Discrete<Boolean>> lessThanOrEquals(Resource<Polynomial> p, Resource<Polynomial> q) {
+    return lessThanOrEquals(subtract(p, q), 0);
+  }
+
+  public static Resource<Polynomial> min(Resource<Polynomial> p, Resource<Polynomial> q) {
+    return ResourceMonad.bind(p, q, (p$, q$) -> ExpiringToResourceMonad.unit(p$.min(q$)));
+  }
+
+  public static Resource<Polynomial> max(Resource<Polynomial> p, Resource<Polynomial> q) {
+    return ResourceMonad.bind(p, q, (p$, q$) -> ExpiringToResourceMonad.unit(p$.max(q$)));
+  }
+
+  public static Resource<Polynomial> abs(Resource<Polynomial> p) {
+    return max(p, negate(p));
+  }
+
+  public static Resource<Polynomial> clamp(Resource<Polynomial> p, Resource<Polynomial> lowerBound, Resource<Polynomial> upperBound) {
+    return ResourceMonad.bind(
+        lessThan(upperBound, lowerBound),
+        impossible -> {
+          if (impossible.extract()) {
+            throw new IllegalStateException(
+                "Inverted bounds for clamp: maximum %f < minimum %f"
+                    .formatted(currentValue(upperBound), currentValue(lowerBound)));
+          }
+          return max(lowerBound, min(upperBound, p));
+        });
+  }
+
+  private static Polynomial scalePolynomial(Polynomial p, double s) {
+    return p.multiply(polynomial(s));
+  }
+
+  public static UnitAware<Resource<Polynomial>> unitAware(Resource<Polynomial> p, Unit unit) {
+    return UnitAwareResources.unitAware(p, unit, PolynomialResources::scalePolynomial);
+  }
+
+  public static UnitAware<CellResource<Polynomial>> unitAware(CellResource<Polynomial> p, Unit unit) {
+    return UnitAwareResources.unitAware(p, unit, PolynomialResources::scalePolynomial);
+  }
+
+  public static UnitAware<Resource<Polynomial>> add(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> q) {
+    return UnitAwareOperations.add(extend(PolynomialResources::scalePolynomial), p, q, PolynomialResources::add);
+  }
+
+  public static UnitAware<Resource<Polynomial>> subtract(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> q) {
+    return UnitAwareOperations.subtract(extend(PolynomialResources::scalePolynomial), p, q, PolynomialResources::subtract);
+  }
+
+  public static UnitAware<Resource<Polynomial>> multiply(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> q) {
+    return UnitAwareOperations.multiply(extend(PolynomialResources::scalePolynomial), p, q, PolynomialResources::multiply);
+  }
+
+  public static UnitAware<Resource<Polynomial>> divide(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Discrete<Double>>> q) {
+    return UnitAwareOperations.divide(extend(PolynomialResources::scalePolynomial), p, q, PolynomialResources::divide);
+  }
+
+  public static UnitAware<Resource<Polynomial>> integrate(UnitAware<Resource<Polynomial>> p, UnitAware<Double> startingValue) {
+    return UnitAwareOperations.integrate(extend(PolynomialResources::scalePolynomial), p, startingValue, PolynomialResources::integrate);
+  }
+
+  public static UnitAware<Resource<Polynomial>> differentiate(UnitAware<Resource<Polynomial>> p) {
+    return UnitAwareOperations.differentiate(extend(PolynomialResources::scalePolynomial), p, PolynomialResources::differentiate);
+  }
+
+  // Ugly $ suffix is to avoid ambiguous overloading after erasure.
+  public static Resource<Discrete<Boolean>> greaterThan$(UnitAware<Resource<Polynomial>> p, UnitAware<Double> threshold) {
+    return greaterThan(p.value(), threshold.value(p.unit()));
+  }
+
+  public static Resource<Discrete<Boolean>> greaterThanOrEquals$(UnitAware<Resource<Polynomial>> p, UnitAware<Double> threshold) {
+    return greaterThanOrEquals(p.value(), threshold.value(p.unit()));
+  }
+
+  public static Resource<Discrete<Boolean>> lessThan$(UnitAware<Resource<Polynomial>> p, UnitAware<Double> threshold) {
+    return lessThan(p.value(), threshold.value(p.unit()));
+  }
+
+  public static Resource<Discrete<Boolean>> lessThanOrEquals$(UnitAware<Resource<Polynomial>> p, UnitAware<Double> threshold) {
+    return lessThanOrEquals(p.value(), threshold.value(p.unit()));
+  }
+
+  public static Resource<Discrete<Boolean>> greaterThan(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> q) {
+    return greaterThan(subtract(p, q).value(), 0);
+  }
+
+  public static Resource<Discrete<Boolean>> greaterThanOrEquals(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> q) {
+    return greaterThanOrEquals(subtract(p, q).value(), 0);
+  }
+
+  public static Resource<Discrete<Boolean>> lessThan(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> q) {
+    return lessThan(subtract(p, q).value(), 0);
+  }
+
+  public static Resource<Discrete<Boolean>> lessThanOrEquals(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> q) {
+    return lessThanOrEquals(subtract(p, q).value(), 0);
+  }
+
+  public static UnitAware<Resource<Polynomial>> min(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> q) {
+    return unitAware(min(p.value(), q.value(p.unit())), p.unit());
+  }
+
+  public static UnitAware<Resource<Polynomial>> max(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> q) {
+    return unitAware(max(p.value(), q.value(p.unit())), p.unit());
+  }
+
+  public static UnitAware<Resource<Polynomial>> clamp(UnitAware<Resource<Polynomial>> p, UnitAware<Resource<Polynomial>> lowerBound, UnitAware<Resource<Polynomial>> upperBound) {
+    return unitAware(clamp(p.value(), lowerBound.value(p.unit()), upperBound.value(p.unit())), p.unit());
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/unit_aware/Dimension.java

@@ -0,0 +1,175 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware;
+
+import java.util.HashMap;
+import java.util.Map;
+import java.util.Objects;
+
+import static java.util.Collections.reverseOrder;
+import static java.util.Map.Entry.comparingByValue;
+import static java.util.stream.Collectors.joining;
+
+/**
+ * A kind of quantity, which can be measured.
+ * For example, length, time, energy, or data rate.
+ *
+ * <p>
+ * Quantities with the same dimension but different units, like meters and miles,
+ * can be compared, added, and subtracted.
+ * Quantities with different dimensions, like meters and seconds,
+ * cannot be added or subtracted, and are never equal.
+ * </p>
+ *
+ * <p>
+ * Base dimensions are declared using {@link Dimension#createBase}. Base dimensions are definitionally all distinct
+ * from each other, and are distinct from all combinations of other base dimensions.
+ * </p>
+ *
+ * <p>
+ * Dimensions can be composed by multiplication, division, and exponentiation by a constant
+ * to derive new dimensions, and composite units correlate to the composite dimension.
+ * For example, Energy is the dimension defined as Mass * Length^2 / Time^2, and
+ * Newton is a unit of Energy defined as Kilogram * Meter^2 / Second^2.
+ * Internally, all dimensions are a map from base dimensions to their power. For example, Mass is stored (loosely) as
+ * {@code {"Mass": 1}} and Energy is {@code {"Mass": 1, "Length": 2, "Time": -2}}.
+ * </p>
+ */
+public sealed interface Dimension {
+  Dimension SCALAR = new DerivedDimension(Map.of());
+
+  Map<BaseDimension, Rational> basePowers();
+  boolean isBase();
+
+
+  default Dimension multiply(Dimension other) {
+    var resultBasePowers = new HashMap<>(basePowers());
+    for (var dimensionPower : other.basePowers().entrySet()) {
+      var power = dimensionPower.getValue();
+      resultBasePowers.compute(
+          dimensionPower.getKey(),
+          (k, p) -> p == null ? power : power.add(p));
+    }
+    return create(resultBasePowers);
+  }
+
+  default Dimension divide(Dimension other) {
+    var resultBasePowers = new HashMap<>(basePowers());
+    for (var dimensionPower : other.basePowers().entrySet()) {
+      var power = dimensionPower.getValue().negate();
+      resultBasePowers.compute(
+          dimensionPower.getKey(),
+          (k, p) -> p == null ? power : power.add(p));
+    }
+    return create(resultBasePowers);
+  }
+
+  default Dimension power(Rational power) {
+    var resultBasePowers = new HashMap<BaseDimension, Rational>();
+    for (var dimensionPower : basePowers().entrySet()) {
+      resultBasePowers.put(dimensionPower.getKey(), dimensionPower.getValue().multiply(power));
+    }
+    return create(resultBasePowers);
+  }
+
+  private static Dimension create(Map<BaseDimension, Rational> basePowers) {
+    var normalizedBasePowers = new HashMap<BaseDimension, Rational>();
+    for (var entry : basePowers.entrySet()) {
+      if (!entry.getValue().equals(Rational.ZERO)) {
+        normalizedBasePowers.put(entry.getKey(), entry.getValue());
+      }
+    }
+
+    if (normalizedBasePowers.isEmpty()) {
+      return Dimension.SCALAR;
+    } else if (normalizedBasePowers.size() == 1) {
+      final var solePower = normalizedBasePowers.entrySet().stream().findAny().get();
+      if (solePower.getValue().equals(Rational.ONE)) {
+        // This actually *is* the base dimension, so return that instead
+        // Normalizing like this lets us bootstrap using reference equality on base dimensions.
+        return solePower.getKey();
+      }
+    }
+
+    // Otherwise, this is some composite dimension, build it anew.
+    return new DerivedDimension(normalizedBasePowers);
+  }
+
+  static Dimension createBase(String name) {
+    return new BaseDimension(name);
+  }
+
+  final class BaseDimension implements Dimension {
+    public final String name;
+
+    private BaseDimension(final String name) {
+      this.name = name;
+    }
+
+    @Override
+    public Map<BaseDimension, Rational> basePowers() {
+      return Map.of(this, Rational.ONE);
+    }
+
+    @Override
+    public boolean isBase() {
+      return true;
+    }
+
+    // Reference equality is sufficient here. Do *not* override equals/hashCode.
+
+    @Override
+    public String toString() {
+      return name;
+    }
+  }
+
+  final class DerivedDimension implements Dimension {
+    private final Map<BaseDimension, Rational> basePowers;
+
+    private DerivedDimension(final Map<BaseDimension, Rational> basePowers) {
+      this.basePowers = basePowers;
+    }
+
+    @Override
+    public Map<BaseDimension, Rational> basePowers() {
+      return basePowers;
+    }
+
+    @Override
+    public boolean isBase() {
+      return false;
+    }
+
+    // Use semantic equality defined by the base powers map
+
+    @Override
+    public boolean equals(final Object o) {
+      if (this == o) return true;
+      if (o == null || getClass() != o.getClass()) return false;
+      DerivedDimension that = (DerivedDimension) o;
+      return Objects.equals(basePowers, that.basePowers);
+    }
+
+    @Override
+    public int hashCode() {
+      return Objects.hash(basePowers);
+    }
+
+    @Override
+    public String toString() {
+      return basePowers.entrySet().stream()
+          .sorted(reverseOrder(comparingByValue()))
+          .map(basePower -> formatBasePower(basePower.getKey(), basePower.getValue()))
+          .collect(joining(" "));
+    }
+
+    private static String formatBasePower(BaseDimension d, Rational p) {
+      if (p.equals(Rational.ONE)) {
+        return "[%s]".formatted(d.name);
+      } else if (p.denominator() == 1) {
+        return "[%s]^%s".formatted(d.name, p.numerator());
+      } else {
+        return "[%s]^(%d/%d)".formatted(d.name, p.numerator(), p.denominator());
+      }
+    }
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/unit_aware/Quantities.java

@@ -0,0 +1,29 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware;
+
+import static gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware.UnitAware.unitAware;
+
+public final class Quantities {
+  public static UnitAware<Double> quantity(double amount, Unit unit) {
+    return unitAware(amount, unit, Quantities::scaling);
+  }
+
+  public static UnitAware<Double> add(UnitAware<Double> p, UnitAware<Double> q) {
+    return UnitAwareOperations.add(Quantities::scaling, p, q, (x, y) -> x + y);
+  }
+
+  public static UnitAware<Double> subtract(UnitAware<Double> p, UnitAware<Double> q) {
+    return UnitAwareOperations.subtract(Quantities::scaling, p, q, (x, y) -> x - y);
+  }
+
+  public static UnitAware<Double> multiply(UnitAware<Double> p, UnitAware<Double> q) {
+    return UnitAwareOperations.multiply(Quantities::scaling, p, q, (x, y) -> x * y);
+  }
+
+  public static UnitAware<Double> divide(UnitAware<Double> p, UnitAware<Double> q) {
+    return UnitAwareOperations.divide(Quantities::scaling, p, q, (x, y) -> x / y);
+  }
+
+  private static double scaling(double x, double y) {
+    return x * y;
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/unit_aware/Rational.java

@@ -0,0 +1,66 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware;
+
+import static java.lang.Integer.signum;
+import static org.apache.commons.math3.util.ArithmeticUtils.gcd;
+
+public record Rational(int numerator, int denominator) implements Comparable<Rational> {
+  public static final Rational ZERO = new Rational(0, 1);
+  public static final Rational ONE = new Rational(1, 1);
+
+  public Rational(final int numerator, final int denominator) {
+    if (denominator == 0) {
+      throw new ArithmeticException("Cannot create a Rational with 0 denominator.");
+    }
+
+    // Normalize by dividing by the GCD and forcing the denominator to be positive.
+    final int gcd = gcd(numerator, denominator);
+    final int s = signum(denominator);
+    this.numerator = s * numerator / gcd;
+    this.denominator = s * denominator / gcd;
+  }
+
+  public static Rational rational(final int numerator, final int denominator) {
+    return new Rational(numerator, denominator);
+  }
+
+  public static Rational rational(final int value) {
+    return new Rational(value, 1);
+  }
+
+  public Rational add(Rational other) {
+    return new Rational(
+        numerator * other.denominator + denominator * other.numerator,
+        denominator * other.denominator);
+  }
+
+  public Rational negate() {
+    return new Rational(-numerator, denominator);
+  }
+
+  public Rational subtract(Rational other) {
+    return this.add(other.negate());
+  }
+
+  public Rational multiply(Rational other) {
+    return new Rational(
+        numerator * other.numerator,
+        denominator * other.denominator);
+  }
+
+  public Rational invert() {
+    return new Rational(denominator, numerator);
+  }
+
+  public Rational divide(Rational other) {
+    return this.multiply(other.invert());
+  }
+
+  @Override
+  public int compareTo(final Rational o) {
+    return Integer.compare(numerator * o.denominator, denominator * o.numerator);
+  }
+
+  public double doubleValue() {
+    return ((double) numerator) / denominator;
+  }
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/unit_aware/StandardDimensions.java

@@ -0,0 +1,18 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware;
+
+public final class StandardDimensions {
+  private StandardDimensions() {}
+
+  // SI base dimensions:
+  public static final Dimension TIME = Dimension.createBase("Time");
+  public static final Dimension LENGTH = Dimension.createBase("Length");
+  public static final Dimension MASS = Dimension.createBase("Mass");
+  public static final Dimension CURRENT = Dimension.createBase("Current");
+  public static final Dimension TEMPERATURE = Dimension.createBase("Temperature");
+  public static final Dimension LUMINOUS_INTENSITY = Dimension.createBase("Luminous Intensity");
+  public static final Dimension AMOUNT = Dimension.createBase("Amount of Substance");
+
+  // Additional base dimensions we've found useful in practice
+  public static final Dimension INFORMATION = Dimension.createBase("Information");
+  public static final Dimension ANGLE = Dimension.createBase("Angle");
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/unit_aware/StandardUnits.java

@@ -0,0 +1,28 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware;
+
+public final class StandardUnits {
+  private StandardUnits() {}
+
+  // Base units, should correspond 1-1 with base Dimensions
+  public static final Unit SECOND = Unit.createBase("s", "second", StandardDimensions.TIME);
+  public static final Unit METER = Unit.createBase("m", "meter", StandardDimensions.LENGTH);
+  public static final Unit KILOGRAM = Unit.createBase("kg", "kilogram", StandardDimensions.MASS);
+  public static final Unit AMPERE = Unit.createBase("A", "ampere", StandardDimensions.CURRENT);
+  public static final Unit KELVIN = Unit.createBase("K", "Kelvin", StandardDimensions.TEMPERATURE);
+  public static final Unit CANDELA = Unit.createBase("cd", "candela", StandardDimensions.LUMINOUS_INTENSITY);
+  public static final Unit MOLE = Unit.createBase("mol", "mole", StandardDimensions.AMOUNT);
+  public static final Unit BIT = Unit.createBase("b", "bit", StandardDimensions.INFORMATION);
+  public static final Unit RADIAN = Unit.createBase("rad", "radian", StandardDimensions.ANGLE);
+
+  // REVIEW: What derived units should be included here?
+  // Including a few arbitrarily to show a few different styles of derivation
+  public static final Unit MILLISECOND = Unit.derived("ms", "millisecond", 1e-3, SECOND);
+  public static final Unit MINUTE = Unit.derived("min", "minute", 60, SECOND);
+  public static final Unit HOUR = Unit.derived("hr", "hour", 60, MINUTE);
+  public static final Unit BYTE = Unit.derived("B", "byte", 8, BIT);
+  public static final Unit NEWTON = Unit.derived("N", "newton", KILOGRAM.multiply(METER).divide(SECOND.power(2)));
+  public static final Unit MEGABIT_PER_SECOND = Unit.derived("Mbps", "megabit per second", 1e6, BIT.divide(SECOND));
+
+  public static final Unit JOULE = Unit.derived("J", "joule", NEWTON.multiply(METER));
+  public static final Unit WATT = Unit.derived("W", "watt", JOULE.divide(SECOND));
+}

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/unit_aware/Unit.java

@@ -0,0 +1,101 @@
+package gov.nasa.jpl.aerie.contrib.streamline.modeling.unit_aware;
+
+import java.util.Objects;
+
+public final class Unit {
+  public final Dimension dimension;
+  public final double multiplier;
+  public final String longName;
+  public final String shortName;
+
+  private Unit(final Dimension dimension, final double multiplier) {
+    this(dimension, multiplier, null, null);
+  }
+
+  private Unit(final Dimension dimension, final double multiplier, final String longName, final String shortName) {
+    this.dimension = dimension;
+    this.multiplier = multiplier;
+    this.longName = longName;
+    this.shortName = shortName;
+  }
+
+  public static Unit createBase(final String shortName, final String longName, final Dimension dimension) {
+    // TODO: Track base units, to detect and prevent collisions
+    assert(dimension.isBase());
+    return new Unit(dimension, 1, longName, shortName);
+  }
+
+  /**
+   * Create a "local" unit. Local units denote a locally-relevant concept that doesn't need to be integrated across models into the broader unit system.
+   *
+   * <p>
+   * Local units are given their own unique base dimension with the same name, so are distinct from all other base dimensions.
+   * </p>
+   *
+   * <p>
+   *     For example, to record that instrument A takes 6 observations per hour, we could declare a local unit in the instrument A model for observations and use it to derive "observations / hour", like so:
+   *     <br/>
+   *     <code>
+   *         Unit Observations = Unit.createLocalUnit("observations");<br/>
+   *         UnitAware&lt;Double&gt; observationRate = quantity(6, Observations.divide(HOUR));
+   *     </code>
+   *     <br/>
+   *     If instrument B also declared a local unit called "observations", this would be incommensurate with instrument A's "observations" unit.
+   * </p>
+   */
+  public static Unit createLocalUnit(final String name) {
+    return createBase(name, name, Dimension.createBase(name));
+  }
+
+  // Used for naming compound units
+  public static Unit derived(final String shortName, final String longName, final Unit definingUnit) {
+    return derived(shortName, longName, 1, definingUnit);
+  }
+
+  public static Unit derived(final String shortName, final String longName, final UnitAware<Double> definingQuantity) {
+    return derived(shortName, longName, definingQuantity.value(), definingQuantity.unit());
+  }
+
+  public static Unit derived(final String shortName, final String longName, final double multiplier, final Unit baseUnit) {
+    return new Unit(baseUnit.dimension, baseUnit.multiplier * multiplier, longName, shortName);
+  }
+
+  public Unit multiply(Unit other) {
+    return new Unit(dimension.multiply(other.dimension), multiplier * other.multiplier);
+  }
+
+  public Unit divide(Unit other) {
+    return new Unit(dimension.divide(other.dimension), multiplier / other.multiplier);
+  }
+
+  public Unit power(int power) {
+    return power(Rational.rational(power));
+  }
+
+  public Unit power(Rational power) {
+    return new Unit(dimension.power(power), Math.pow(multiplier, power.doubleValue()));
+  }
+
+  @Override
+  public boolean equals(final Object o) {
+    if (this == o) return true;
+    if (o == null || getClass() != o.getClass()) return false;
+    Unit unit = (Unit) o;
+    return Double.compare(unit.multiplier, multiplier) == 0 && Objects.equals(dimension, unit.dimension);
+  }
+
+  @Override
+  public int hashCode() {
+    return Objects.hash(dimension, multiplier);
+  }
+
+  @Override
+  public String toString() {
+     if (longName != null) {
+       return longName;
+     } else {
+       // TODO: Better name derivation, by tracking named base units
+       return "Unit{" + multiplier + " in " + dimension + "}";
+     }
+  }
+}