diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 67427aa..6c25047 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -233,11 +233,15 @@ jobs:
       # The point is not to publish, it is that a docstring which no longer
       # parses -- a bad cross-reference, a broken code block -- fails here
       # rather than silently degrading. Output is not committed (#30).
+      #
+      # The module list is derived, not written here: `pdoc pyimr` alone documents
+      # only the package page, because `__init__` sets `__all__`, and the explicit
+      # list this replaced had drifted to seven modules out of nineteen.
       - name: Build API documentation
         if: matrix.python-version == '3.12'
         run: |
           python -m pip install "pdoc>=15"
-          python -m pdoc -o site pyimr pyimr.sensitivity pyimr.inference pyimr.data pyimr.design pyimr.pymc_op pyimr.assimilation pyimr.optimize
+          python tools/api_docs.py -o site
 
       - name: Build distributions
         if: matrix.python-version == '3.12'
diff --git a/README.md b/README.md
index 24656a2..69b71ce 100644
--- a/README.md
+++ b/README.md
@@ -54,8 +54,10 @@ constitutive selector, no shared bag of parameters.
 Materials compose: pick an elastic law and a viscous law and combine them, or
 reach for a closed-form memory model when it applies. Neo-Hookean, Mooney-Rivlin,
 Yeoh, Fung, Gent, Arruda-Boyce and Ogden on the elastic side; Carreau-Yasuda,
-Cross, Powell-Eyring, Herschel-Bulkley and Bingham on the viscous; Zener,
-Oldroyd-B, Giesekus and linear PTT for memory. See
+Cross, Powell-Eyring, Herschel-Bulkley and Bingham on the viscous; the Zener
+family (linear, quadratic, cubic, two-mode and Carreau-thinning), Oldroyd-B and
+linear Maxwell in closed form; Giesekus and linear PTT as distributed memory.
+`RelaxingMaterial` puts a Maxwell arm on any elastic law. See
 **[docs/materials.md](docs/materials.md)**.
 
 ## Beyond a forward solve
@@ -76,6 +78,19 @@ redundancy and Occam penalties.
 tolerance meeting an accuracy target on *your* problem, and raises rather than
 guessing when the target is out of reach.
 
+**Experiment design.** `pyimr.design` scores a design that has not been run;
+`pyimr.measure` optimises over a batch of them and returns a certificate of
+global optimality with the answer; `pyimr.gain` scores any question, a material
+parameter or a model label, in nats, so one batch can serve several;
+`pyimr.discriminate` screens which rivals are still worth an experiment;
+`pyimr.pareto` traces the front when the criteria disagree.
+
+**Knowing why a solve failed.** `pyimr.diagnose` separates a step budget, a
+material-domain violation and an ill-conditioned trajectory, which raise the
+same error. `pyimr.store` caches solves and failures on disk so a re-run of an
+unchanged study costs nothing. `pyimr.assimilation` runs 4D-Var and ensemble
+state estimation on the exact tangents.
+
 See **[docs/usage.md](docs/usage.md)**.
 
 ## Validation
@@ -102,7 +117,7 @@ visible.
 
 | | |
 |---|---|
-| [Usage](docs/usage.md) | solving, sensitivities, inference, model selection, resolution |
+| [Usage](docs/usage.md) | solving, sensitivities, inference, model selection, experiment design, diagnostics |
 | [Materials](docs/materials.md) | every constitutive law, and what each one requires |
 | [Accuracy](docs/accuracy.md) | what error each tolerance and discretization actually buys |
 | [Discretization](docs/discretization.md) | stress quadrature and the two thermal backends |
@@ -110,7 +125,8 @@ visible.
 | [Upstream](docs/upstream.md) | defects found in IMRv2, and what PyIMR does instead |
 | [Boundaries](docs/boundaries.md) | where PyIMR stops, and where it diverges deliberately |
 
-API reference: `pip install 'PyIMR[docs]'`, then `python -m pdoc pyimr`.
+API reference: `pip install 'PyIMR[docs]'`, then `python tools/api_docs.py`
+(every public module; `python -m pdoc pyimr` alone stops at the package page).
 
 ## Citation
 
diff --git a/docs/accuracy.md b/docs/accuracy.md
index b775142..5e9d806 100644
--- a/docs/accuracy.md
+++ b/docs/accuracy.md
@@ -47,9 +47,10 @@ augmented state/tangent system is integrated. On the coupled fd case at
 
 | `rtol` / `atol` | relative error |
 |---|---|
-| `1e-9` / `1e-11` (default) | 8.43e-05 |
+| `1e-9` / `1e-11` (the suite's setting for this check) | 8.43e-05 |
 | `1e-12` / `1e-14` | 1.53e-06 |
 
+The package default is `1e-8` / `1e-10`, one order looser than the first row.
 So three orders of tolerance buys a factor of about 55, at a large cost in
 runtime -- the coupled tangent solve is already the slowest operation in the
 package. Tightening further stops helping, because the centered-difference
diff --git a/docs/materials.md b/docs/materials.md
index 671af2c..9582b79 100644
--- a/docs/materials.md
+++ b/docs/materials.md
@@ -85,8 +85,8 @@ and the equivalent composable material agree to solver tolerance.
 
 ## Closed-form memory
 
-The finite-dimensional hot paths are `Zener`, `QuadraticZener`, `OldroydB` and
-`LinearMaxwell`:
+The finite-dimensional hot paths are `Zener`, `QuadraticZener`, `CubicZener`,
+`QuadraticKelvinVoigt`, `OldroydB` and `LinearMaxwell`:
 
 ```python
 material = Zener(
@@ -97,6 +97,43 @@ material = Zener(
 )
 ```
 
+`QuadraticKelvinVoigt` is `NeoHookeanKelvinVoigt` with a stiffening term
+quadratic in `I1 - 3`, which is Yeoh truncated after `c2`; `QuadraticZener`
+adds the Maxwell arm to it, and `Zener` is `QuadraticZener` at zero stiffening.
+`CubicZener` carries the next term of the same expansion, with `cubic = 0`
+reducing to `QuadraticZener` exactly. On measured collapses `I1 - 3` reaches
+24 to 119, where the cubic term is not a small correction. `LinearMaxwell` is
+Zener without the parallel spring: no modulus, no retardation, an elastic
+target of zero.
+
+### Relaxation for any elastic law
+
+Each Zener above carries its own hand-derived stress integral, so relaxation
+was available to the neo-Hookean family and to nothing else. Gent, Fung,
+Arruda-Boyce, Mooney-Rivlin, Yeoh and Ogden existed only as
+`InstantaneousMaterial`, elastic with no memory, which is the wrong comparison
+to make against a relaxing model. `RelaxingMaterial` is the Zener construction
+freed from one potential: any `ElasticModel` as the equilibrium target, with a
+Maxwell arm on top.
+
+```python
+from pyimr import RelaxingMaterial, Yeoh
+
+material = RelaxingMaterial(
+    elastic=Yeoh(c1_pa=1250.0, c2_pa=100.0, c3_pa=10.0),
+    viscosity_pa_s=0.1,
+    relaxation_time_s=40e-6,
+    retardation_time_s=8e-6,
+)
+```
+
+The equilibrium target is taken by the same quadrature `InstantaneousMaterial`
+uses rather than by a closed form derived per law, so one class covers every
+elastic model in the package and any added later. That costs the quadrature at
+each step (`quadrature_points`, default 32) and buys a controlled comparison.
+With `elastic=NeoHookean(...)` it agrees with `Zener`; the closed form is the
+one to use where it applies.
+
 ### Beyond one relaxation time
 
 A single relaxation time is a strong assumption for a crosslinked biopolymer, and the
@@ -143,7 +180,7 @@ material = CarreauZener(
 
 Both are comparison candidates, in `EXTENDED_MODELS` rather than `STANDARD_MODELS`: at six
 free parameters the grid quadrature in `pyimr.selection` costs `count**6`, so they are
-scored by `candidate_log_evidence` instead. See [selection](../README.md#model-selection).
+scored by `candidate_log_evidence` instead. See [model selection](usage.md#model-selection).
 
 ## Distributed nonlinear memory
 
diff --git a/docs/upstream.md b/docs/upstream.md
index 621f595..2ec6ffc 100644
--- a/docs/upstream.md
+++ b/docs/upstream.md
@@ -1,7 +1,7 @@
 # Reference implementation
 
 Defects found in IMRv2 at `dea31cd`, all reproduced with MATLAB R2025a via
-`tools/gen_imrv2_cases.m` and `tools/probe_viscosity.m`. The full list is below;
+`tools/gen_imrv2_cases.m` and `tools/probe_viscosity.m`. The eight are listed below;
 the original scoping notes are in git history, in a `PLAN.md` retired in #218.
 
 - **Giesekus and linear PTT cannot be run.** `f_call_params.m` dispatches
@@ -34,7 +34,8 @@ the original scoping notes are in git history, in a `PLAN.md` retired in #218.
   as `A -> 0`; `f_radial_eq.m` takes `(-b + sqrt(d))/(2a)`, which is the
   `-1/nog` branch -- a 32.5% density deficit at ambient pressure, a 48-60%
   enthalpy error, and a negative `c^2`. That negative `c^2` is why `radial = 6`
-  returns complex radii: it is the only branch that evaluates the sound speed
+  (Gilmore/Mie-Gruneisen) returns complex radii, reaching `max|imag(R/R0)| =
+  4.069` without raising: it is the only branch that evaluates the sound speed
   from the EoS. The branch also omits the stress term from `Pb`, which
   `radial = 3` and `4` both include.
 
@@ -48,22 +49,16 @@ the original scoping notes are in git history, in a `PLAN.md` retired in #218.
 - **`f_init_stress.m` uses an undefined `z1`** in the `De == 0 || De == Inf`
   branch. Unreachable for the memory models that call it, so latent rather
   than active.
+- **`calc_omega_N` treats the gas pressure at `Rmax` as the equilibrium
+  value.** This one is in IMR-vanilla rather than IMRv2. It inflates the
+  linearised stiffness by `alpha**(-3*kappa)` and overpredicts the natural
+  frequency by 42x on the reference case, which is why PyIMR's
+  `data.natural_frequency` is a reimplementation rather than a port.
 
 These are the reason several PyIMR models are validated by reduction limit
 rather than against a pinned upstream trajectory: for those models, no working
 upstream implementation exists to pin against.
 
-[Back to the README](../README.md)
-
-- **`calc_omega_N` (IMR-vanilla) treats the gas pressure at `Rmax` as the
-  equilibrium value.** That inflates the linearised stiffness by
-  `alpha**(-3*kappa)` and overpredicts the natural frequency by 42x on the
-  reference case, which is why PyIMR's `data.natural_frequency` is a
-  reimplementation rather than a port.
-- **`radial = 6` (Gilmore/Mie-Gruneisen) returns complex radii.** Upstream
-  reaches `max|imag(R/R0)| = 4.069` without raising, from a wrong root of the
-  Mie-Gruneisen density quadratic.
-
 ## Which branches replicate upstream, and which correct it
 
 Moved here from the package docstring, where four of its claims had gone stale
@@ -114,3 +109,5 @@ is internally inconsistent, and the reduction limit to `LinearMaxwell` converges
 only with the `LAM` factor restored. Three Zener reference trajectories were
 regenerated from PyIMR as a result, and pin regressions rather than
 cross-checking upstream (#174, IMRv2#18).
+
+[Back to the README](../README.md)
diff --git a/docs/usage.md b/docs/usage.md
index 493ef9b..bb50874 100644
--- a/docs/usage.md
+++ b/docs/usage.md
@@ -70,7 +70,7 @@ config = SimulationConfig(
 ## Sensitivities
 
 The tangent-linear solver differentiates the production RHS rather than a
-reduced surrogate. It covers every operator but `gilmore/mie-gruneisen`, every typed material, thermal
+reduced surrogate. It covers every operator, every typed material, thermal
 and mass-transfer states, distributed nonlinear memory, forcing, geometry,
 initial conditions, and continuous physical parameters.
 
@@ -139,7 +139,14 @@ against experimental radius data; keep `1e-10, 1e-12` for sensitivities and
 
 Tolerance does not bound how long one solve can take. `max_steps` (default
 `1_000_000`) turns a trajectory that will not finish into a `SimulationError` at
-a point of your choosing, which is what makes a grid sweep affordable.
+a point of your choosing, which is what makes a grid sweep affordable. Three
+guards refuse a trajectory the model has left rather than returning it:
+`max_radius_ratio` (default 50) raises when the bubble runs away,
+`min_radius_ratio` (off by default) when it collapses past where the material
+is trustworthy, and `max_wall_mach` (off by default) when a rebound approaches
+the sonic line. The last two are different guards for a reason;
+[boundaries.md](boundaries.md) has the measurement. When a solve fails and the
+message does not say why, `pyimr.diagnose` does (below).
 
 `Nt` and tolerance requirements depend on record length, material stiffness and
 which observable is fitted, so a setting adequate for one collapse can be badly
@@ -223,8 +230,9 @@ just as confident. Worked studies are in `examples/`.
 ### Models the grid cannot reach
 
 `solve_grid` is a Cartesian product at one count on every axis, so it costs
-`count**dimension` and runs out at four or five parameters. At `count = 12`, six
-axes is 5,971,968 solves against 168,072 for the whole of `STANDARD_MODELS`.
+`count**dimension` and runs out at four or five parameters. At `count = 12`, one
+six-axis candidate is 2,985,984 solves against 139,272 for the whole of
+`STANDARD_MODELS`.
 
 Candidates past that limit live in `EXTENDED_MODELS` and are scored by expansion
 about a fit instead of by quadrature over a grid:
@@ -252,10 +260,6 @@ coordinates the prior is uniform on, which is both where the Occam factor has to
 be measured and what lets it handle candidates whose axes are not material fields
 -- `qSLS2` has `tau_ratio`, a ratio of two of them, and `oldroydb` likewise.
 
-You must supply the fitted point. Nothing in the package yet locates one for a
-candidate, which is the current gap between having these models and being able to
-rank them; it is tracked in the issue list.
-
 ### The forward operator is a model choice too
 
 `DYNAMICS_MODELS` names the fourteen operators, as `(dynamics, liquid_eos)` pairs. Six
@@ -279,8 +283,8 @@ for dynamics, liquid_eos in DYNAMICS_MODELS:
 
 The parameter space is identical across the set, so the Occam terms cancel and the difference
 in log evidence is a Bayes factor between operators. Every candidate in this package assumes
-`dynamics="keller-miksis"`; on the records analysed in the companion analysis repository, two other operators beat it. See
-the issue list for what that comparison does and does not establish --- in
+`dynamics="keller-miksis"`; on the records analysed in the companion analysis repository, two other operators beat it.
+[#294](https://github.com/sbryngelson/PyIMR/issues/294) records what that comparison does and does not establish --- in
 particular, it must be run in identified coordinates, or the ranking follows the prior box
 rather than the data, and the operator is the most absorbed of the model axes, so most of an
 operator change can be mimicked by refitting the material.
@@ -301,7 +305,139 @@ form to round-off wherever every direction is already sharper than the prior, so
 it changes an answer only where the plain form was not entitled to one.
 `candidate_log_evidence` turns it on, because it compares across dimensions.
 
-## Designing when the criteria disagree
+## Designing experiments
+
+Before a record exists, the question is which experiment to run. Four modules
+answer it. Every number they return is in nats, so the answers can be compared
+and added.
+
+### Scoring one design
+
+`pyimr.design` scores a design that has not been run. A `DesignInference` is a
+`PreparedInference` whose observations are placeholders: the same
+configuration, time grid, noise level and parameter box, with no data.
+`expected_information_gain` averages the Laplace/Fisher gain over prior draws
+and returns its Monte Carlo error bar alongside:
+
+```python
+from pyimr.design import design_inference, expected_information_gain
+
+design = design_inference(config, times, standard_deviation_m=2e-6, parameters=parameters)
+score = expected_information_gain(design, draws=128, workers=4)
+score.expected_information_gain, score.standard_error   # nats
+```
+
+`optimize_design` in `pyimr.optimize` searches a continuous design space for
+the largest gain. You supply `build_inference(design) -> DesignInference`, so a
+design is whatever you say it is -- pulse amplitude, window, radius, or a
+mixture -- and the error bar is passed to the surrogate as observation noise so
+that a point which scored well by luck is not chased:
+
+```python
+from pyimr.optimize import optimize_design
+
+search = optimize_design(build_inference, [(0.0, 8e4)], draws=64, evaluations=24)
+search.best_point, search.best_value
+```
+
+### A batch with a certificate
+
+A single best design is the answer to the wrong question: an experimenter runs
+`n_1` bubbles at one setting and `n_2` at another, and a search over one point
+is not convex, so nothing in its result says whether it stopped at the optimum.
+`pyimr.measure` optimises over a probability measure on a candidate set instead.
+The averaged information matrix is linear in the weights, the criterion is
+concave, and the Kiefer-Wolfowitz equivalence theorem turns the first-order
+condition into a proof of global optimality: `gap` is the largest directional
+derivative toward any candidate, and it is zero at the optimum.
+
+```python
+import numpy as np
+from pyimr.design import design_information
+from pyimr.measure import apportion, optimal_measure
+
+matrices = []
+for R0 in (100e-6, 200e-6, 400e-6):
+    design = design_inference(config_at(R0), times, 2e-6, parameters)
+    information, requested, failed = design_information(design, draws=64)
+    matrices.append(information.mean(axis=0))   # one prior-averaged J^T J per candidate
+
+measure = optimal_measure(np.array(matrices))
+measure.weights, measure.support, measure.gap    # weights [0, 0.54, 0.46], gap ~1e-10
+
+batch = apportion(measure.weights, 12, np.array(matrices))
+batch.counts, batch.efficiency                   # [0, 6, 6], D-efficiency 0.998
+```
+
+`apportion` is the efficient rounding of Pukelsheim and Rieder. Passing
+`matrices` makes it measure the D-efficiency of the integer batch against the
+measure rather than assume the rounding was harmless, which for a dozen runs it
+need not be.
+
+A plain measure is not always the right object. An information criterion is
+happiest concentrating, and a batch on one setting cannot detect that every
+model is wrong, so `constrained_measure` forces weight onto at least `settings`
+distinct candidates. Counting runs asserts that every candidate costs the same,
+which fails once a trace of `2N` frames competes with a trace of `N`, so
+`budgeted_measure` normalises by cost. `identification_front` returns integer
+batches that trade parameter precision against separating two models. Each
+reports whether its certificate still holds under the constraint.
+
+### What would this experiment change your mind about
+
+`log det` scores the material parameters. A Schur complement scores the model
+label. Neither says which question a batch should serve, because they do not
+share a scale. `pyimr.gain` puts every question in nats: the expected
+information gain about any subset of coordinates, with the rest treated as
+nuisance, is `(1/2) log det` of a Schur complement of `I + M` in
+prior-standardised coordinates. It is concave, so `optimal_measure` still
+certifies it, and it is finite when `M` is singular, so a design that cannot
+determine every parameter scores low rather than `-inf`.
+
+```python
+from pyimr.gain import Question, expected_gain, gain_criterion
+
+questions = [Question("modulus", (0,)), Question("viscosity", (1,), weight=0.5)]
+expected_gain(matrices[0], questions).per_question           # nats, unweighted
+measure = optimal_measure(np.array(matrices), criterion=gain_criterion(questions))
+```
+
+Nats say how much a batch teaches, not whether it settles anything.
+`runs_to_settle` inverts the Gaussian-linear Bayes factor distribution to say
+how many runs clear a threshold in favour of the truth at a stated confidence;
+`runs_to_precision` asks the same of one parameter; `lack_of_fit_degrees`
+reports whether the batch leaves any degrees of freedom to discover that every
+model in the catalogue is wrong.
+
+### Which rivals still matter
+
+`pyimr.discriminate` scores model discrimination as an integral rather than a
+minimum. T-optimality takes `min` over the rival's parameters, and that inner
+problem is multimodal here: a local method landing in the wrong basin returns a
+wrong answer with nothing to signal it. `expected_log_bayes_factor` replaces it
+with the expected log Bayes factor between prior-predictive banks, and every
+result carries an effective sample size so the silent collapse of a
+prior-sample evidence onto one draw is visible. `laplace_log_evidence` is the
+fallback for when it collapses, and is what `candidate_log_evidence` uses.
+
+`screen_models` decides which rivals are worth designing for at all. A rival
+the records have already decided against by more than `decisive` nats needs no
+experiment, and a design that spends runs separating it is spending them on a
+settled matter:
+
+```python
+from pyimr.discriminate import screen_models
+
+screen = screen_models(evidences)   # one log evidence per model, on the data in hand
+screen.live, screen.decided, screen.weights
+```
+
+The survivors' weights feed `gain_criterion` as question weights, and the
+close pairs are what `measure.augmented_information` and `separability` are for:
+the model label becomes one more Jacobian column, and what matters is how much
+of the difference between two models survives refitting the material.
+
+### When the criteria disagree
 
 `optimize_design` maximises one number. That is the right question only while
 the design criteria happen to agree, and on the qSLS study they do not.
@@ -348,6 +484,76 @@ qSLS-versus-SLS separation by about an order of magnitude — the honest figure
 the present design is 0.668 noise units, meaning SLS imitates qSLS to well within
 the noise.
 
+## Diagnosing a failed solve
+
+`SimulationError: maximum number of solver steps was reached` is true and
+nearly useless: it is the same message whether the trajectory needs a bigger
+budget, violates a material's domain, or is so sensitive to its inputs that no
+tolerance will be met, and those want opposite responses. `pyimr.diagnose` runs
+the ladder that distinguishes them:
+
+```python
+from pyimr.diagnose import diagnose
+
+report = diagnose(config, times)
+report.outcome    # ok, runaway, ill-conditioned, domain, budget, unresolved
+report.summary
+```
+
+It is several solves, meant for a point that already misbehaved, not for a
+sweep. A solve that succeeds is checked too. `ill-conditioned` means a `1e-9`
+change in `R0` is amplified until the trajectory has shed digits, and the
+result is one draw from a sensitive system rather than an answer. `budget`
+means a larger `max_steps` finishes. `domain` means the material refused a state
+the collapse reached. `runaway` means no tolerance or solver will help.
+`unresolved` means a loose tolerance completes and the requested one does not,
+which is worth a look by hand.
+
+## Caching solves
+
+A parameter study is re-run many times while the analysis around it changes.
+`ResultStore` keys on the content of `(times, config)`, so the second run is
+free, and it caches failures too, which is most of the value: a point that
+exhausts its step budget spends the whole budget before saying so, every time
+it is asked.
+
+```python
+from pyimr.store import ResultStore
+
+store = ResultStore("~/.cache/pyimr")
+result = store.simulate(times, config)   # same signature and result as pyimr.simulate
+```
+
+The key cannot see a change to PyIMR itself; pass `version=` to invalidate a
+directory by hand after one.
+
+## State estimation
+
+`pyimr.assimilation` estimates the initial state of a prepared problem from a
+window of observations. `four_dvar` minimises the strong-constraint 4D-Var
+cost with exact gradients from the tangent-linear operator, which is why true
+4D-Var is reachable here rather than an ensemble approximation to it.
+`ensemble_smoother` answers the same question from ensemble statistics,
+`ienks` re-linearises about the current estimate each iteration and closes on
+the minimum `four_dvar` finds, and `enkf_analysis`, `ensemble_update` and
+`kalman_analysis` are the one-step analyses those are built from. The smoothers
+take the `PreparedProblem`, the times, the observations and a linear
+observation operator; `noise` is a scalar or vector of standard deviations, or
+a full covariance. The docstrings give the shapes.
+
+## Running in parallel
+
+Every `workers=` argument in the package goes through `pyimr.parallel`. It
+exists because XLA sizes its CPU thread pool from the process's affinity mask
+and OpenBLAS does the same at import: sixteen spawn workers on a 128-core host
+were measured asking for roughly 6,500 threads between them. `worker_pool` gives
+each worker one core and one thread per library, `map_work` decides serial or
+pooled by timing the first item rather than counting them, and
+`default_workers` reads the affinity mask rather than `cpu_count`, so a job on a
+shared node claims only what the scheduler granted. `PYIMR_WORKERS` overrides
+the default. Set the thread-count environment variables before importing numpy
+if you build your own pool; `limit_worker_threads` does it for you.
+
 ## Trace estimators
 
 `pyimr.data` covers the step before inference: getting from a measured `R(t)`
diff --git a/pyimr/measure.py b/pyimr/measure.py
index ed7876d..63ae908 100644
--- a/pyimr/measure.py
+++ b/pyimr/measure.py
@@ -444,11 +444,19 @@ def separability(jacobian, differences, *, weights=None):
   """
   material, columns, amplitude = _augmented_parts(jacobian, differences, weights)
   count, rivals = material.shape[1], columns.shape[0]
-  try:
-    variance = np.diag(np.linalg.inv(_augmented_gram(material, columns, amplitude)))[count:].copy()
-  except np.linalg.LinAlgError:
-    variance = np.full(rivals, np.inf)
-  variance = np.where(variance > 0.0, variance, np.inf)
+  # Not `inv` with `LinAlgError` as the unidentified signal. A difference in the span of the
+  # material makes the Gram singular only to roundoff, and LAPACK raises only on an exactly
+  # zero pivot: OpenBLAS on x86 happened to produce one, Accelerate on arm64 produced 1e-14
+  # and a variance of 2.8e14 where the promise is `inf`. The null space is found at
+  # `matrix_rank`'s tolerance instead, and a coordinate with a component in it is unidentified.
+  gram = _augmented_gram(material, columns, amplitude)
+  values, vectors = np.linalg.eigh(gram)
+  tolerance = values.max() * gram.shape[0] * np.finfo(float).eps
+  kept = values > tolerance
+  null = vectors[count:, ~kept]
+  unidentified = np.linalg.norm(null, axis=1) > np.sqrt(np.finfo(float).eps)
+  pseudo = (vectors[:, kept] / values[kept]) @ vectors[:, kept].T
+  variance = np.where(unidentified, np.inf, np.diag(pseudo)[count:])
 
   basis, _ = np.linalg.qr(material)
   absorbed = np.empty(rivals)
diff --git a/pyimr/pymc_op.py b/pyimr/pymc_op.py
index 58ecc08..b677c53 100644
--- a/pyimr/pymc_op.py
+++ b/pyimr/pymc_op.py
@@ -10,10 +10,12 @@
 _MISSING = "pyimr.pymc_op requires PyMC: pip install 'PyIMR[inference]'"
 
 def _pymc():
+  # PyMC is an optional extra, so the type check has to pass in an environment without it;
+  # the ImportError below is the runtime answer, and these ignores are the static one.
   try:
-    import pymc
-    import pytensor.tensor as tensor
-    from pytensor.graph.op import Op
+    import pymc  # pyright: ignore[reportMissingImports]
+    import pytensor.tensor as tensor  # pyright: ignore[reportMissingImports]
+    from pytensor.graph.op import Op  # pyright: ignore[reportMissingImports]
   except ImportError as error:  # pragma: no cover - exercised only without pymc
     raise ImportError(_MISSING) from error
   return pymc, tensor, Op
diff --git a/tests/test_pymc_op.py b/tests/test_pymc_op.py
index 3da3821..c179bd4 100644
--- a/tests/test_pymc_op.py
+++ b/tests/test_pymc_op.py
@@ -54,7 +54,7 @@ def test_one_solve_serves_both_halves(inference, monkeypatch):
   calls = []
   real = type(inference).evaluate_with_jacobian
   monkeypatch.setattr(type(inference), "evaluate_with_jacobian", lambda self, unit: (calls.append(np.asarray(unit).copy()), real(self, unit))[1])
-  import pytensor.tensor as tensor
+  import pytensor.tensor as tensor  # pyright: ignore[reportMissingImports] - optional extra, see importorskip above
 
   operation = pymc_op.IMRLogLikelihood(inference)
   unit = tensor.as_tensor_variable(np.array([0.5, 0.5]))
diff --git a/tests/test_validation_estimators.py b/tests/test_validation_estimators.py
index 347a92d..fed5696 100644
--- a/tests/test_validation_estimators.py
+++ b/tests/test_validation_estimators.py
@@ -255,11 +255,17 @@ def test_the_correlated_gradient_is_still_the_derivative(correlated, measured):
 
 
 def test_vanishing_correlation_time_reduces_to_independent_noise(correlated):
-  """The limit that must hold exactly, not approximately: at tau -> 0 the"""
+  """At tau -> 0 the kernel is the identity and the two likelihoods are the same number.
+
+  To roundoff, not bitwise: the correlated path sums its log-determinant as
+  `n log 2pi + 2 sum log L_ii` and whitens by a triangular solve, the independent path as
+  `sum log(2 pi s^2)` and a division, and the two agree to an ulp that moved between scipy
+  1.17 and 1.18.
+  """
   times, observed, independent, _ = correlated
   tiny = prepare_inference(independent.config, FieldObservation("radius_m", times, observed, 5e-7, correlation_time_s=1e-15), independent.parameters)
   unit = np.array([0.42, 0.37])
-  assert tiny.evaluate(unit).log_likelihood == independent.evaluate(unit).log_likelihood
+  assert tiny.evaluate(unit).log_likelihood == pytest.approx(independent.evaluate(unit).log_likelihood, rel=1e-12)
 
 
 def test_correlated_noise_carries_less_information(correlated, measured):
diff --git a/tools/api_docs.py b/tools/api_docs.py
new file mode 100644
index 0000000..76a7dae
--- /dev/null
+++ b/tools/api_docs.py
@@ -0,0 +1,33 @@
+"""Build the API reference over every public module, not just the package.
+
+`pyimr/__init__.py` defines `__all__`, and pdoc takes that as the complete list of
+what to document -- so `python -m pdoc pyimr` produces one page and silently omits
+`pyimr.inference`, `pyimr.selection` and the rest. The CI step carried a hand-written
+module list instead, which was forgotten as modules were added: seven names out of
+nineteen by the time anyone looked.
+
+The list is derived here, by the rule the package already uses for "public": a
+module that does not start with an underscore and declares `__all__`.
+
+  python tools/api_docs.py            serve at http://localhost:8080
+  python tools/api_docs.py -o site    write HTML to site/
+"""
+
+import importlib
+import pathlib
+import pkgutil
+import subprocess
+import sys
+
+# So the source tree works uninstalled too; pip's editable install makes this a no-op.
+sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+import pyimr
+
+
+def public_modules():
+  names = (f"pyimr.{m.name}" for m in pkgutil.iter_modules(pyimr.__path__) if not m.name.startswith("_"))
+  return ["pyimr"] + [n for n in names if hasattr(importlib.import_module(n), "__all__")]
+
+
+if __name__ == "__main__":
+  sys.exit(subprocess.call([sys.executable, "-m", "pdoc", *sys.argv[1:], *public_modules()]))
diff --git a/tools/validate_bubtherm_adiabatic.py b/tools/validate_bubtherm_adiabatic.py
index 97d5cae..7d8a958 100644
--- a/tools/validate_bubtherm_adiabatic.py
+++ b/tools/validate_bubtherm_adiabatic.py
@@ -19,6 +19,11 @@
 initial condition and only the RHS structure is being tested.
 """
 
+import sys
+from pathlib import Path
+
+# Runnable from a checkout without `pip install -e .`, as benchmarks/run.py is.
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
 import numpy as np
 from scipy.integrate import solve_ivp
 
diff --git a/tools/validate_thermal_fd.py b/tools/validate_thermal_fd.py
index 1dd64d0..138464a 100644
--- a/tools/validate_thermal_fd.py
+++ b/tools/validate_thermal_fd.py
@@ -13,6 +13,11 @@
 behaviour explicit rather than hiding it behind a loose tolerance.
 """
 
+import sys
+from pathlib import Path
+
+# Runnable from a checkout without `pip install -e .`, as benchmarks/run.py is.
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
 import numpy as np
 
 from pyimr.thermal_fd import finite_diff_mat