Made the tolerance trait separable from the integrator itself

2026-09-30 13:48:20 -04:00
1 parent 812b5adaf9
commit da479a378e
@@ -0,0 +1,174 @@
+# Feature: Runtime Integrator Selection
+
+**Status**: proposed (2026-09-30)
+
+## Overview
+
+Which integrator runs a problem is fixed at compile time, because `Integrator<D>`
+keeps its structural facts in associated constants:
+
+```rust
+pub trait Integrator<const D: usize> {
+    const ORDER: usize;
+    const STAGES: usize;
+    const ADAPTIVE: bool;
+    const DENSE: bool;
+    const MAIN_STAGES: usize = Self::STAGES;
+    const EXTRA_STAGES: usize = 0;
+    const TOTAL_DENSE_STAGES: usize = Self::MAIN_STAGES + Self::EXTRA_STAGES;
+}
+```
+
+`Problem` is generic over `S`, and `Problem::solve` reads `S::ADAPTIVE` to choose
+its stepping branch and `S::EXTRA_STAGES` / `S::TOTAL_DENSE_STAGES` to decide
+whether lazy stages still need computing. A value cannot answer those questions —
+only a type can.
+
+So the obvious workaround does not work. An enum wrapping the four integrators
+cannot implement `Integrator<D>`, because it would have to give one
+`EXTRA_STAGES` for all of them, and Vern7 has six where Dormand-Prince has none.
+
+`Box<dyn Integrator<D>>` is not available either, and will not become available:
+`step<T: Scalar, P>` and `interpolate<T: Scalar>` are generic methods, so the
+trait is not object-safe. Erasing the scalar would defeat the genericity that
+makes autodiff work. The enum is the practical route, not `dyn`.
+
+## Why This Feature Matters
+
+This surfaced building `astrokit` on top of this crate. Three things a dependent
+wants and cannot have:
+
+1. **Accuracy presets that include a method.** `Accuracy::High` should mean "tight
+   tolerances *and* a high-order method", since the two belong together — a 7th
+   order method at 1e-6 is waste, and a 3rd order method at 1e-12 is hopeless. As
+   it stands a preset can only set tolerances, and the method has to be named
+   separately as a type parameter.
+
+2. **An integrator named in a configuration file.** `astrokit` is built for a
+   study where every run parameter must be declared in a committed file and every
+   propagation must record enough to recreate itself. `integrator: "Vern7"`
+   deserialized from that file cannot select anything, because the choice is a
+   type. The integrator ends up recorded in the provenance but not settable from
+   it, which is backwards.
+
+3. **A heterogeneous collection.** Benchmarking or comparing methods over the same
+   problem means either a macro over method types or duplicated code per method,
+   where a `Vec<Method>` would do.
+
+None of these is fatal. `astrokit` has settled on Vern7 as a type-level default,
+which is the right method for its accuracy range anyway. But the constraint is
+structural rather than deliberate, and it is worth fixing before more dependents
+grow around it.
+
+## Dependencies
+
+None, though it interacts with [Tolerances Trait](40-tolerances-trait.md): an
+enum wrapper needs to forward tolerance setters too, which is much easier once
+those are trait methods.
+
+## Implementation Approach
+
+### Step 1: constants become methods
+
+The structural facts move from associated constants to methods taking `&self`:
+
+```rust
+pub trait Integrator<const D: usize> {
+    fn order(&self) -> usize;
+    fn stages(&self) -> usize;
+    fn adaptive(&self) -> bool;
+    fn dense(&self) -> bool;
+
+    fn main_stages(&self) -> usize { self.stages() }
+    fn extra_stages(&self) -> usize { 0 }
+    fn total_dense_stages(&self) -> usize { self.main_stages() + self.extra_stages() }
+
+    // step, interpolate, compute_extra_stages unchanged
+}
+```
+
+**This costs nothing measurable.** Every use of these constants today is already a
+runtime comparison, not a compile-time size:
+
+- `if S::ADAPTIVE` in `Problem::solve` becomes `if self.integrator.adaptive()`,
+  which monomorphizes to the same constant and folds the same way for a concrete
+  integrator.
+- `borrowed.len() < 2 + S::TOTAL_DENSE_STAGES` compares against a `Vec` length.
+- `&dense[2..2 + S::MAIN_STAGES]` slices a `Vec`.
+
+None of them appears in type position, so nothing needs the values to be consts.
+For a concrete `S` the method calls inline to the same constants they are now.
+
+This is a breaking change to anyone implementing `Integrator` by hand — but a
+mechanical one, and this crate has no external implementors.
+
+### Step 2: an enum in this crate
+
+With step 1 done, an enum can dispatch:
+
+```rust
+/// An integrator chosen at runtime.
+#[non_exhaustive]
+#[derive(Debug, Clone, Copy)]
+pub enum Method<const D: usize> {
+    BS3(BS3<D>),
+    DormandPrince45(DormandPrince45<D>),
+    Vern7(Vern7<D>),
+    Rosenbrock23(Rosenbrock23<D>),
+}
+
+impl<const D: usize> Integrator<D> for Method<D> {
+    fn extra_stages(&self) -> usize {
+        match self {
+            Method::Vern7(m) => m.extra_stages(),
+            Method::BS3(m) => m.extra_stages(),
+            ..
+        }
+    }
+    // and so on, including step and interpolate
+}
+```
+
+Each method body is a `match` forwarding to the wrapped integrator. Verbose, and
+a good candidate for a macro over the variant list.
+
+Worth deriving `Serialize`/`Deserialize` on a companion `MethodName` enum so a
+method can be named in a configuration file, and adding
+`MethodName::instantiate(&self) -> Method<D>`. That is what makes point 2 above
+work.
+
+### Step 3: presets
+
+Once a method is a value, an accuracy preset can carry one, pairing a tolerance
+with a method whose order suits it. That belongs in whichever crate owns the
+notion of a preset — probably this one, given
+[Default Algorithm Selection](13-default-algorithm-selection.md) is already on
+the roadmap and needs the same machinery.
+
+### The cost of not doing step 1
+
+If the constants must stay, the fallback is a parallel object-ish layer: a second
+trait with `&self` methods, blanket-implemented for everything that implements
+`Integrator`, and `Problem` generic over that instead. It works, it duplicates the
+trait surface, and every implementor carries both. Step 1 is cheaper and cleaner.
+
+## Testing Requirements
+
+- Every existing test passes after step 1, with method calls substituted for
+  constant reads.
+- A benchmark confirms no regression from the constants becoming methods — the
+  claim is that they inline, and that should be measured rather than assumed, on
+  the existing `simple_1d` and `orbit` benches.
+- `Method` gives bit-identical results to the integrator it wraps, for each
+  variant, on a problem with dense output interpolation (which exercises the
+  lazy-stage path where `extra_stages` matters).
+- A method name round-trips through serde and instantiates the right variant.
+
+## Success Criteria
+
+- `Problem` can be constructed with an integrator chosen at runtime.
+- No measurable performance change on the existing criterion benches.
+- A dependent crate can offer an accuracy preset that selects both a tolerance and
+  a method, and can read a method name out of a configuration file.
+
+**Effort**: Medium — step 1 is small and mechanical, step 2 is repetitive