Record accepted-mesh provenance and enforce replay compatibility

CHANGELOG.md

@@ -14,6 +14,12 @@ versioning while the project remains pre-1.0.
   The capability matrix separates method, loop, sparse, and hardware evidence
   and documents matching `just qualify` commands (#206).
 
+- Live adaptive schedules carry lightweight method/layout/kernel provenance and
+  descriptive discovery controls, namespace, precision and caller regime labels.
+  Replay rejects known structural mismatches while retaining portable dynamic
+  parameter/state reuse and compatible manual meshes. A structured analytic
+  example demonstrates accuracy-driven refresh outside differentiation (#205).
+
 ### Fixed
 
 - Namespace-aware advection and diffusion builders honor a floating reference

docs/guides/backends.md

@@ -65,8 +65,11 @@ The resulting derivative is conditional on the recorded mesh. This is also the
 branchwise meaning of an eager live-adaptive gradient: neither mode differentiates
 the discrete accept/reject decision. Refresh the schedule as optimization
 parameters move, and always refresh it after a material model or tolerance
-change. A schedule is validated against the output grid, while matching the
-method and other configuration is the caller's responsibility.
+change. A recorded schedule checks producing method, state layout, relevant kernel
+configuration, and optional structural model identity before replay. Backend,
+dtype and fitted values can change; manually constructed meshes remain
+unverified unless provenance is required. See the
+[replay contract and refresh example](replay-provenance.md).
 
 Core schedule replay retains the eager Python loop by default. Select
 `RunConfig(replay_loop="auto")` to use a namespace's registered scan driver,

docs/guides/replay-provenance.md

@@ -0,0 +1,65 @@
+# Accepted meshes: compatibility and accuracy
+
+A frozen accepted mesh is static controller output. Replay invokes its numerical
+steps and bypasses adaptive error estimates and acceptance decisions. Its AD
+results are derivatives **conditional on that mesh**. A compatible mesh may
+still become inaccurate when model parameters or initial conditions move.
+
+The core contract distinguishes these questions before testing accuracy:
+
+| Metadata | Replay treatment |
+| --- | --- |
+| Finite/increasing output grid and positive interval steps summing to each interval | Always validated by `AdaptiveStepSchedule` |
+| Producing resolved method, state shape/axis order, relevant operator axis, TR-BDF2 gamma, nonlinear stage solver/configuration | Known mismatches reject before stepping |
+| Caller-supplied structural model identity | Compared when recording/replay use an identity; it must describe equations/geometry/layout, not fitted parameter values |
+| Array namespace and dtype | Descriptive; portable meshes may change NumPy/JAX/Torch or precision |
+| Discovery tolerances, controller limits, and reference parameter regime | Descriptive; changing fitted values or requested accuracy requires an external accuracy decision, not an automatic incompatibility |
+| Callback bodies and arbitrary custom solver internals | No automatic semantic identity guarantee; use the structural identity and independently qualify changed implementations |
+
+`AdaptiveScheduleProvenance` is attached to schedules created by live adaptive
+runs. `RunConfig.replay_model_id` and `replay_reference_regime` are optional
+static labels supplied by the caller. Record concrete nominal values or a
+stable parameter-set identifier in the reference label at an eager boundary;
+replay does not inspect or convert dynamic parameter arrays.
+
+Manually constructed meshes and legacy schedules with `provenance=None` remain
+accepted by default after mesh/output-grid validation. Their method/model
+compatibility is unverified and remains the caller's responsibility. Use
+`RunConfig(replay_require_provenance=True)` to reject unverified schedules.
+Adding a provenance record is not an accuracy certificate, and copying labels
+from another model does not make that model compatible.
+
+Controller settings are recorded to explain discovery. Replay ignores error
+control even if the supplied configuration requests tighter tolerances. A
+provider may impose stronger configuration/context checks on its own artifact;
+these core rules do not remove those provider checks.
+
+A practical reuse policy is:
+
+1. Discover a live mesh at a documented nominal parameter regime.
+2. Reuse it while changing dynamic values, retaining structural model identity.
+3. Compare a replay trajectory to a fresh accurate solve or to a refined mesh
+   at representative new parameter values. Check the scientific objective and
+   needed derivatives as well as state error.
+4. If the independent discrepancy exceeds the intended accuracy, rediscover
+   the mesh at the new regime, recheck agreement, and freeze that mesh for the
+   next differentiation phase.
+
+The structured example accompanying this guide executes that policy and shows
+an accuracy failure despite compatible provenance. Discrete refresh decisions
+remain outside the differentiated numerical solve.
+
+Run the [structured refresh example](https://github.com/ACCIDDA/op_engine/blob/main/examples/replay_mesh_refresh.py):
+
+```shell
+uv run python examples/replay_mesh_refresh.py
+```
+
+It transfers mass between two generic pools with `A'=-k*A^2`, `B'=+k*A^2`,
+a shared rate and group-local multipliers. The independent solution is
+`A(t)=A0/(1+k*A0*t)`, with conserved `A+B`. A mesh discovered at rate 0.05
+remains accurate nearby; rate 3.0 is structurally compatible but exceeds the
+example's `1e-4` state-error policy on the old mesh. Fresh discovery restores
+that accuracy. JSON reports both discrepancies and accepted-step counts.
+This policy checks this fixture's state error; different scientific objectives
+and derivative accuracy need their own independent checks.

examples/replay_mesh_refresh.py

@@ -0,0 +1,122 @@
+"""Check and refresh a compatible mesh after a nonlinear parameter change.
+
+Run ``uv run python examples/replay_mesh_refresh.py``. A shared transfer rate
+and two group-local multipliers drive a generic closed two-pool model. The
+analytic solution independently checks replay accuracy before any mesh reuse
+is trusted. Refresh is an eager decision outside a differentiated solve.
+"""
+
+from __future__ import annotations
+
+import json
+from dataclasses import replace
+from typing import cast
+
+import numpy as np
+
+from op_engine import Array, CoreSolver, ModelCore, array_namespace
+from op_engine.core_solver import AdaptiveConfig, AdaptiveStepSchedule, RunConfig
+from op_engine.model_core import ModelCoreOptions
+
+INITIAL = np.asarray([[0.8, 0.6], [0.2, 0.4]])
+LOCAL = np.asarray([1.0, 1.3])
+CONFIG = RunConfig(
+    method="heun",
+    adaptive=True,
+    adaptive_cfg=AdaptiveConfig(rtol=1e-4, atol=1e-7, dt_init=0.5),
+    replay_model_id="quadratic-transfer:pool,group:v1",
+    replay_require_provenance=True,
+)
+
+
+def trajectory(
+    rate: float, schedule: AdaptiveStepSchedule | None = None
+) -> tuple[np.ndarray, AdaptiveStepSchedule]:
+    """Discover or replay this structured model at a shared transfer rate.
+
+    Returns:
+        Final state and its accepted mesh.
+
+    Raises:
+        RuntimeError: If live discovery did not return a mesh.
+    """
+    core = ModelCore(
+        2,
+        2,
+        np.asarray([0.0, 1.0]),
+        options=ModelCoreOptions(
+            axis_names=("pool", "group"),
+        ),
+    )
+    core.set_initial_state(INITIAL)
+    solver = CoreSolver(core)
+
+    def rhs(_time: object, state: Array) -> Array:
+        xp = array_namespace(state)
+        flow = rate * xp.asarray(LOCAL, dtype=state.dtype) * state[0] ** 2
+        return cast("Array", xp.stack((-flow, flow)))
+
+    config = replace(
+        CONFIG, replay_reference_regime=f"shared_rate={rate}; local=[1,1.3]; A0=[.8,.6]"
+    )
+    if schedule is None:
+        solver.run(rhs, config=config)
+    else:
+        solver.replay_adaptive_schedule(rhs, schedule, config=config)
+    recorded = solver.last_adaptive_schedule
+    if recorded is None:
+        msg = "Adaptive solve did not return an accepted mesh"
+        raise RuntimeError(msg)
+    return np.asarray(core.get_current_state()), recorded
+
+
+def exact(rate: float) -> np.ndarray:
+    """Evaluate the independent closed-form two-pool solution at time one.
+
+    Returns:
+        Final state from A'=-k*A^2 and the conserved total A+B.
+    """
+    first = INITIAL[0] / (1 + rate * LOCAL * INITIAL[0])
+    return np.stack((first, INITIAL.sum(axis=0) - first))
+
+
+def assess_and_refresh() -> dict[str, float | int | bool | str]:
+    """Apply an independent state-error policy at a changed parameter regime.
+
+    Returns:
+        Accuracy and mesh-count evidence for the reusable and refreshed cases.
+
+    Raises:
+        RuntimeError: If the selected mesh fails the independent accuracy check.
+    """
+    _, nominal = trajectory(0.05)
+    near, _ = trajectory(0.055, nominal)
+    changed_rate = 3.0
+    changed, _ = trajectory(changed_rate, nominal)
+    reference = exact(changed_rate)
+    maximum_error = float(np.max(np.abs(changed - reference)))
+    accuracy_limit = 1e-4
+    selected = nominal
+    result = changed
+    if maximum_error > accuracy_limit:
+        result, selected = trajectory(changed_rate)
+    refreshed_error = float(np.max(np.abs(result - reference)))
+    if refreshed_error > accuracy_limit:
+        msg = "The selected mesh still fails the independent accuracy check"
+        raise RuntimeError(msg)
+    return {
+        "model_id": CONFIG.replay_model_id or "",
+        "nominal_rate": 0.05,
+        "changed_rate": changed_rate,
+        "accuracy_limit": accuracy_limit,
+        "nearby_error": float(np.max(np.abs(near - exact(0.055)))),
+        "old_mesh_error": maximum_error,
+        "selected_mesh_error": refreshed_error,
+        "nominal_steps": sum(map(len, nominal.step_sizes)),
+        "selected_steps": sum(map(len, selected.step_sizes)),
+        "refreshed": selected is not nominal,
+    }
+
+
+if __name__ == "__main__":
+    print(json.dumps(assess_and_refresh(), indent=2))  # noqa: T201

justfile

@@ -102,6 +102,7 @@ build-test-core:
 	uv pip install --python "${CLEANROOM}/venv/bin/python" -r "${CLEANROOM}/dev-requirements.txt"
 	cp pyproject.toml "${CLEANROOM}/pyproject.toml"
 	cp -R tests "${CLEANROOM}/tests"
+	cp -R examples "${CLEANROOM}/examples"
 	cd "${CLEANROOM}"
 	"${CLEANROOM}/venv/bin/pytest" --import-mode=importlib tests --quiet --exitfirst
 

mkdocs.yml

@@ -57,6 +57,8 @@ nav:
       - Biogeochemical Network: guides/biogeochemical-network.md
       - Trait-Structured Chemostat: guides/chemostat-pde.md
       - Backend Qualification: guides/backend-qualification.md
+
+      - Replay Provenance: guides/replay-provenance.md
       - Backend Strategies: guides/backends.md
   - API Reference:
       - Adaptive Tau: api-reference/adaptive_tau.md

src/op_engine/__init__.py

@@ -9,6 +9,7 @@
 from ._typing import Scalar as Scalar
 from .adaptive_tau import AdaptiveTauLeapingConfig, AdaptiveTauLeapingSolver
 from .core_solver import (
+    AdaptiveScheduleProvenance,
     CoreSolver,
     NonlinearIntegrationConvergenceError,
     NonlinearIntegrationDiagnostics,
@@ -88,6 +89,7 @@
 )
 
 __all__ = [
+    "AdaptiveScheduleProvenance",
     "AdaptiveTauLeapingConfig",
     "AdaptiveTauLeapingSolver",
     "Array",