earthsci_ast/

problem.rs

//! The EsmProblem / `solve` surface — one noun and one verb
//! (`esm-libraries-spec.md` §2.5, `API_SPEC.md` §5.8).
//!
//! ```text
//! let prob = esm_problem(&file, (0.0, 10.0), ProblemOptions::default())?;  // build once
//! let sol  = solve(&prob, &SolveOptions::default())?;                      // run per knob-set
//! ```
//!
//! [`ProblemOptions`] is taken **by value**: it owns the boxed data providers
//! the build hands to the runtime, and there is nothing left in it worth
//! reusing afterwards. [`SolveOptions`] is taken by reference, because a caller
//! sweeping knobs reuses it.
//!
//! ## What replaced what
//!
//! `simulate` is **gone**, in all of its forms — the one-shot free function,
//! `simulate_with_inspection`, `simulate_with_providers_inspect`, and
//! `Prepared`/`prepare`. It conflated two operations whose costs differ by
//! orders of magnitude and whose inputs differ in kind, which is exactly why
//! this crate had grown a second, `prepare`-shaped entry point beside it.
//!
//! [`esm_problem`] absorbs the whole deterministic-per-document pipeline: the
//! projection-pushdown rewrite, value invention, the gated fetch of provider
//! data, CONST-provider materialization, and the compile of the right-hand
//! side. [`solve`] varies only per-run knobs — `alg`, `abstol`, `reltol`,
//! `saveat`, `callback`, `maxiters`.
//!
//! ## The solver is optional
//!
//! Per §2.5.9, **constructing a EsmProblem does not require the solver.**
//! [`esm_problem`], [`remake`], [`callbacks`] and [`observed_field`] are all
//! available with the `solve` Cargo feature switched off; only [`solve`],
//! [`init`], [`step`] and [`solve_to_completion`] need it. `diffsol` is behind
//! that feature and is not an unconditional dependency.

use std::collections::{BTreeMap, HashMap};
use std::path::{Path, PathBuf};
use std::rc::Rc;

use ndarray::ArrayD;
use serde_json::Value as JsonValue;

use crate::classification::DeclaredRank;
use crate::flatten::FlattenedSystem;
use crate::precision::{self, Precision};
#[cfg_attr(not(feature = "solve"), allow(unused_imports))]
use crate::simulate::Solution;
use crate::simulate::{Flow, Progress, ProgressFn, SimulateError, SolveOptions};
use crate::types::EsmFile;

use crate::simulate_array::{ArrayCompiled, BuildInspection};

// =============================================================================
// Callbacks
// =============================================================================

/// One callback: an observer invoked once before the first step and then after
/// every accepted step, with the interval covered so far AND the integrator's
/// current state vector ([`Progress::u`]).
///
/// Returning [`Flow::Cancel`] stops the integration; the solution comes back
/// with [`crate::ReturnCode::Terminated`] and the trajectory computed so far
/// rather than as an error, because a caller who stops a run deliberately still
/// wants what it produced.
#[cfg(target_arch = "wasm32")]
pub type CallbackFn = std::rc::Rc<dyn for<'p> Fn(&Progress<'p>) -> Flow>;
/// One callback. See the `wasm32` twin for the documentation.
#[cfg(not(target_arch = "wasm32"))]
pub type CallbackFn = std::sync::Arc<dyn for<'p> Fn(&Progress<'p>) -> Flow + Send + Sync>;

/// The callbacks declared on a [`EsmProblem`].
///
/// Callbacks live on the EsmProblem, not on a run: a callback that refreshes
/// provider buffers or writes an output stream belongs to the *document*, not
/// to a particular run's tolerances (§2.5.4).
///
/// Entries are named so that [`CallbackSet::names`] can tell a caller what a
/// EsmProblem already carries before they decide whether to replace it.
#[derive(Clone, Default)]
pub struct CallbackSet {
    entries: Vec<(String, CallbackFn)>,
}

impl std::fmt::Debug for CallbackSet {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("CallbackSet")
            .field("names", &self.names())
            .finish()
    }
}

impl CallbackSet {
    /// An empty set.
    pub fn new() -> Self {
        Self::default()
    }

    /// A one-entry set.
    pub fn of(name: impl Into<String>, f: CallbackFn) -> Self {
        let mut s = Self::new();
        s.push(name, f);
        s
    }

    /// Append one named callback.
    pub fn push(&mut self, name: impl Into<String>, f: CallbackFn) -> &mut Self {
        self.entries.push((name.into(), f));
        self
    }

    /// How many callbacks the set holds.
    pub fn len(&self) -> usize {
        self.entries.len()
    }

    /// Whether the set is empty.
    pub fn is_empty(&self) -> bool {
        self.entries.is_empty()
    }

    /// The names, in invocation order.
    pub fn names(&self) -> Vec<&str> {
        self.entries.iter().map(|(n, _)| n.as_str()).collect()
    }

    /// Invoke every callback in order. The first [`Flow::Cancel`] wins, and the
    /// remaining callbacks are still called — an output callback must not be
    /// skipped for the final step just because an earlier watchdog asked to
    /// stop on it.
    pub fn invoke(&self, p: &Progress<'_>) -> Flow {
        let mut flow = Flow::Continue;
        for (_, f) in &self.entries {
            if matches!(f(p), Flow::Cancel) {
                flow = Flow::Cancel;
            }
        }
        flow
    }
}

/// Concatenate two callback sets into a new one.
///
/// This is the explicit composition §2.5.4 requires: a `callback` argument to
/// [`solve`] REPLACES the EsmProblem's set entirely — it does not append, merge or
/// wrap — so a caller who wants to *extend* rather than replace reads the
/// existing set back with [`callbacks`] and composes:
///
/// ```text
/// let cb = compose(callbacks(&prob), &my_extra);
/// solve(&prob, &SolveOptions { callback: Some(cb), ..Default::default() })
/// ```
pub fn compose(a: &CallbackSet, b: &CallbackSet) -> CallbackSet {
    let mut out = a.clone();
    out.entries.extend(b.entries.iter().cloned());
    out
}

// =============================================================================
// Construction inputs and options
// =============================================================================

/// What [`esm_problem`] can be built from: a path on disk, a raw JSON
/// document, a typed document, or an already-flattened system.
pub enum ProblemInput<'a> {
    /// A path to an `.esm` file.
    Path(&'a Path),
    /// A raw (un-typed) JSON document — the form the build pipeline rewrites.
    Json(&'a JsonValue),
    /// A typed document.
    File(&'a EsmFile),
    /// An already-flattened system (scalar ODE path only).
    Flattened(&'a FlattenedSystem),
}

impl<'a> From<&'a Path> for ProblemInput<'a> {
    fn from(p: &'a Path) -> Self {
        ProblemInput::Path(p)
    }
}
impl<'a> From<&'a JsonValue> for ProblemInput<'a> {
    fn from(v: &'a JsonValue) -> Self {
        ProblemInput::Json(v)
    }
}
impl<'a> From<&'a EsmFile> for ProblemInput<'a> {
    fn from(f: &'a EsmFile) -> Self {
        ProblemInput::File(f)
    }
}
impl<'a> From<&'a FlattenedSystem> for ProblemInput<'a> {
    fn from(f: &'a FlattenedSystem) -> Self {
        ProblemInput::Flattened(f)
    }
}

/// Whether [`esm_problem`] compiles a right-hand side.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum Rhs {
    /// Compile when the selected document declares differential equations,
    /// skip when it does not. A document with no ODEs — a dispatched static
    /// evaluation — still gets a EsmProblem, whose build-time products
    /// [`observed_field`] reads; [`solve`] on it raises
    /// [`SimulateError::NotDynamic`].
    #[default]
    Auto,
    /// Always compile; a document with no differential equations is an error.
    Always,
    /// Never compile. Build-time products only.
    Never,
}

/// The former spelling of [`Rhs`], kept for one minor (API_SPEC §8 item 22).
///
/// `Compile` and [`Compiler`] differ by one letter and answer different
/// questions — `Rhs` decides WHETHER a right-hand side is built at all, and
/// `Compiler` decides WHICH strategy builds it — so the two sat on the same
/// `ProblemOptions` as `compile` and `compiler`. The type is `extension` tier,
/// so the rename is a minor.
#[deprecated(
    since = "0.2.0",
    note = "renamed to `Rhs`; the field is now `ProblemOptions::rhs`"
)]
pub type Compile = Rhs;

/// Which strategy builds a Problem's right-hand side (API_SPEC §5.8,
/// esm-libraries-spec §2.5.10).
///
/// The vocabulary is CLOSED and the same in every binding. A value this
/// binding does not provide is REFUSED — never answered by quietly selecting
/// another one, because the point of naming a compiler is that the caller
/// knows which one ran.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum Compiler {
    /// The universally fast option with no heavy external dependency, and the
    /// default for that reason: the array runtime's TAPE, for every document
    /// whatever its shape — never a routing decision taken on document
    /// content.
    ///
    /// **Strict.** A rule the tape cannot lower is a construction error naming
    /// the rule and the deepest decline reason
    /// ([`crate::compile_error::CompileError::CompilerRefusedRule`]), not a
    /// quiet demotion to the per-cell oracle.
    ///
    /// It covers the whole evaluable core. The only documents it may refuse
    /// are those carrying a construct outside its scope, each refused by name:
    /// events, implicit equations and SDE noise (`unsupported_construct`) and
    /// a spatial operator with no discretization rule (`unlowered_operator`,
    /// refused under every compiler). Every compiler-driven conformance tier
    /// requires `native` for every fixture, so any other refusal is a defect
    /// in this crate.
    #[default]
    Native,
    /// The reference, and only the reference: the array runtime with the tape
    /// AND the whole-array overlay off, so every rule is evaluated by the
    /// per-cell oracle. Complete over the evaluable core, and carrying no
    /// performance promise of any kind. A caller selects it to CHECK another
    /// compiler, not to run a model.
    Interpreter,
    /// StableHLO through XLA: the specialty compiler that needs a heavy
    /// external dependency to exist at all.
    ///
    /// The tape is lowered to an `rhs(u, p, t) -> du` XLA computation
    /// (`crate::simulate_array::tape::xla_emit`) and run through PJRT
    /// (`crate::xla_runtime`), so the Problem's right-hand side — and the
    /// finite-difference Jacobian an implicit solve differences out of it —
    /// ARE the compiled executable. The passes that are not the right-hand
    /// side (the build-time materialization of constants and static
    /// observeds, the per-segment seed, the inspection snapshot, and the
    /// observeds reported at output times) are served from the SAME tape the
    /// emitter was built from, because the emitted program's only output is
    /// `du`. That is one evaluator, not two.
    ///
    /// **Strict twice over, and never a fallback.** A rule the tape cannot
    /// lower and an instruction the emitter cannot lower are both
    /// [`crate::compile_error::CompileError::CompilerRefusedRule`] at
    /// construction, naming the rule and the reason; neither demotes to the
    /// tape interpreter or to the per-cell oracle.
    ///
    /// It needs the `xla` Cargo feature AND a usable XLA runtime in this
    /// process (an unpacked `xla_extension` release — see
    /// `scripts/fetch-xla-extension.sh`). Without either, `esm_problem`
    /// answers [`crate::SimulateError::CompilerUnavailable`], never a quiet
    /// substitution.
    Xla,
    /// A ModelingToolkit system. Julia only; refused here with
    /// `compiler_unavailable`.
    Mtk,
    /// A lambdified SymPy scalar right-hand side. Python only; refused here
    /// with `compiler_unavailable`.
    Sympy,
}

impl Compiler {
    /// The closed vocabulary's spelling of this value — the one a diagnostic,
    /// a CLI argument and a conformance manifest all use.
    pub fn as_str(self) -> &'static str {
        match self {
            Compiler::Native => "native",
            Compiler::Interpreter => "interpreter",
            Compiler::Xla => "xla",
            Compiler::Mtk => "mtk",
            Compiler::Sympy => "sympy",
        }
    }

    /// Parse a vocabulary member. `None` for anything outside it — the caller
    /// raises `compiler_unknown`, which is a different failure from
    /// `compiler_unavailable`.
    pub fn parse(name: &str) -> Option<Compiler> {
        match name {
            "native" => Some(Compiler::Native),
            "interpreter" => Some(Compiler::Interpreter),
            "xla" => Some(Compiler::Xla),
            "mtk" => Some(Compiler::Mtk),
            "sympy" => Some(Compiler::Sympy),
            _ => None,
        }
    }

    /// [`Self::parse`] with the `compiler_unknown` diagnostic (esm-spec
    /// §9.6.6) already written, for the entry points that take the compiler as
    /// TEXT — the CLI flag, the wasm handle, a conformance adapter.
    ///
    /// One wording and one place the closed vocabulary is listed: a value
    /// outside it is a different failure from one inside it this binding
    /// cannot provide, and a caller who conflates them cannot tell a typo from
    /// a missing runtime.
    pub fn parse_named(name: &str) -> Result<Compiler, String> {
        Compiler::parse(name).ok_or_else(|| {
            format!(
                "{}: '{name}' is not a compiler. The vocabulary is closed: {}",
                crate::diagnostic::codes::COMPILER_UNKNOWN,
                Compiler::vocabulary()
                    .iter()
                    .map(|c| c.as_str())
                    .collect::<Vec<_>>()
                    .join(", ")
            )
        })
    }

    /// Whether this compiler is STRICT: a rule it cannot express is a
    /// construction error naming the rule, never a quiet demotion to a slower
    /// evaluator (API_SPEC §5.8, "The default is `native`, and it is strict").
    ///
    /// Both of the values this binding implements over the tape are strict —
    /// `native`, which is the tape, and `xla`, which is the emitter over it —
    /// and they are strict for the same reason: a compiler either runs the
    /// whole document or names what it refused. `interpreter` is not strict
    /// because it declines nothing it can evaluate; it IS the fallback the
    /// other two refuse to take.
    pub fn is_strict(self) -> bool {
        matches!(self, Compiler::Native | Compiler::Xla)
    }

    /// Every member, in the order API_SPEC §5.8 tabulates them.
    pub fn vocabulary() -> &'static [Compiler] {
        &[
            Compiler::Native,
            Compiler::Interpreter,
            Compiler::Xla,
            Compiler::Mtk,
            Compiler::Sympy,
        ]
    }
}

impl std::fmt::Display for Compiler {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.write_str(self.as_str())
    }
}

/// Where ONE rule of a Problem landed, and at what cadence (API_SPEC §5.8,
/// "Every Problem reports what ran").
#[derive(Clone, Debug)]
pub struct CompilerRuleReport {
    /// The rule — an equation or an observed — component-qualified, the same
    /// spelling [`observed_field`] resolves and a refusal names.
    pub rule: String,
    /// `"observed"` or `"state derivative"` for a rule of the compiled model.
    /// The evaluations construction performs OUTSIDE that rule set are
    /// reported too, under their own kinds: `"build-time observed"` and
    /// `"value invention"` (the build pipeline), and `"initial condition"` and
    /// `"initial-condition scope"` (a field `ic` and the state-free observeds
    /// it reads).
    pub kind: &'static str,
    /// The cadence tier the rule runs at: `"const"` (once per solve, at
    /// setup), `"discrete"` (once per refresh segment) or `"continuous"`
    /// (every right-hand-side call). The tier is the cost story: a decline at
    /// `const` costs once, the same decline at `continuous` costs per step.
    pub cadence: &'static str,
    /// Where the rule landed: `"taped"`, `"xla"` (lowered from the tape into
    /// the compiled XLA program, under [`Compiler::Xla`]), `"fallback"` (the
    /// per-cell oracle beneath the tape, reachable only outside a strict
    /// compiler), or `"oracle"` (every rule under [`Compiler::Interpreter`],
    /// by design). An evaluation outside the compiled rule set is
    /// `"vectorized"` when the whole-array overlay served it and `"oracle"`
    /// when it was walked per cell, which a strict compiler refuses. A
    /// `"value invention"` row is `"relational"`: the producer's member set,
    /// computed once at setup by the relational engine from build-time factors
    /// (CONFORMANCE_SPEC §5.7.6 keeps that engine off the hot path), under
    /// every compiler alike. It is neither the oracle nor a fallback, so a
    /// build that succeeds under a strict compiler still has
    /// [`CompilerReport::n_oracle`] `== 0`.
    pub tier: &'static str,
    /// For a `"fallback"`, the DEEPEST decline reason reached while trying to
    /// lower the rule; for an `"oracle"` evaluation outside the compiled rule
    /// set, why it was walked per cell. `None` otherwise.
    pub reason: Option<String>,
}

/// What built a Problem's right-hand side, per rule (API_SPEC §5.8).
///
/// The names and their meanings are stable surface; the record's exact shape
/// is per-binding and is deliberately not pinned across the five.
///
/// **Fusion is reported per PROGRAM, not per rule.** The tape's Step-4 fusion
/// pass merges instructions across rule boundaries and the fused program keeps
/// no rule attribution, so a `tier` of `"taped"` covers both a rule that runs
/// as discrete instructions and one absorbed into a fused group;
/// [`Self::fused_groups`] and [`Self::fused_instructions`] report the fusion
/// at the level the tape actually tracks it.
#[derive(Clone, Debug)]
pub struct CompilerReport {
    compiler: Compiler,
    rules: Vec<CompilerRuleReport>,
    fused_groups: usize,
    fused_instructions: usize,
}

impl CompilerReport {
    /// The report of a Problem with no right-hand side: no rules, because no
    /// compiler was asked to build anything.
    fn empty(compiler: Compiler) -> CompilerReport {
        CompilerReport {
            compiler,
            rules: Vec::new(),
            fused_groups: 0,
            fused_instructions: 0,
        }
    }

    /// The compiler that built the Problem.
    pub fn compiler(&self) -> Compiler {
        self.compiler
    }

    /// Every rule, in program order: the evaluations construction performs
    /// outside the compiled model (the build pipeline's, then a state-free
    /// document's observeds, then the field initial conditions), then the
    /// compiled model's observed rules in dependency order, then its state
    /// derivatives. A state-free document's observeds are here too: it has no
    /// right-hand side, but its observed graph is still evaluated by the
    /// compiler the Problem names.
    pub fn rules(&self) -> &[CompilerRuleReport] {
        &self.rules
    }

    /// How many rules landed on the tape.
    pub fn n_taped(&self) -> usize {
        self.rules.iter().filter(|r| r.tier == "taped").count()
    }

    /// How many rules were lowered out of the tape into the compiled XLA
    /// program. Non-zero only under [`Compiler::Xla`], where it is every rule
    /// — the build refuses the document otherwise.
    pub fn n_xla(&self) -> usize {
        self.rules.iter().filter(|r| r.tier == "xla").count()
    }

    /// How many rules landed on the per-cell oracle — as a fallback beneath
    /// the tape, or, under [`Compiler::Interpreter`], by design.
    ///
    /// Counts the two oracle tiers by NAME rather than everything that is not
    /// `"taped"`: an `"xla"` rule is compiled, and counting it here would
    /// report the fastest tier in the crate as the slowest.
    pub fn n_oracle(&self) -> usize {
        self.rules
            .iter()
            .filter(|r| matches!(r.tier, "oracle" | "fallback"))
            .count()
    }

    /// How many evaluations outside the compiled rule set the whole-array
    /// overlay served (tier `"vectorized"`).
    pub fn n_vectorized(&self) -> usize {
        self.rules.iter().filter(|r| r.tier == "vectorized").count()
    }

    /// Fused instruction groups in the compiled program. `0` under
    /// [`Compiler::Interpreter`], which compiles no program.
    pub fn fused_groups(&self) -> usize {
        self.fused_groups
    }

    /// Instructions absorbed into those groups.
    pub fn fused_instructions(&self) -> usize {
        self.fused_instructions
    }
}

impl std::fmt::Display for CompilerReport {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(
            f,
            "compiler {}: {} rules, {} taped, {} on the per-cell oracle",
            self.compiler,
            self.rules.len(),
            self.n_taped(),
            self.n_oracle()
        )?;
        if self.n_xla() > 0 {
            write!(f, ", {} lowered into the XLA program", self.n_xla())?;
        }
        if self.n_vectorized() > 0 {
            write!(f, ", {} on the whole-array overlay", self.n_vectorized())?;
        }
        let n_relational = self.rules.iter().filter(|r| r.tier == "relational").count();
        if n_relational > 0 {
            write!(f, ", {n_relational} by the relational engine at setup")?;
        }
        if self.fused_groups > 0 {
            write!(
                f,
                "; {} fused groups absorbing {} instructions",
                self.fused_groups, self.fused_instructions
            )?;
        }
        Ok(())
    }
}

/// Build-time bindings for [`esm_problem`] — everything that fixes a
/// *document* rather than a *run*.
///
/// The per-run knobs (`alg`, `abstol`, `reltol`, `saveat`, `callback`,
/// `maxiters`) live on [`SolveOptions`] instead. That split is the whole point
/// of §2.5: construction is deterministic per document and expensive, `solve`
/// is cheap and varies.
#[derive(Default)]
pub struct ProblemOptions {
    /// Parameter bindings (canonical SciML `p`). Exact or bare names.
    pub p: HashMap<String, f64>,
    /// Initial state bindings (canonical SciML `u0`). Exact or bare names.
    pub u0: HashMap<String, f64>,
    /// Select one model when the document holds several.
    pub model_name: Option<String>,
    /// Metaparameter bindings closed at load (esm-spec §9.7.6 site 4).
    pub metaparameters: BTreeMap<String, i64>,
    /// Base path anchoring relative `{ref}`s.
    pub base_path: Option<PathBuf>,
    /// Callbacks declared on the EsmProblem (§2.5.4).
    pub callbacks: CallbackSet,
    /// Whether to compile a right-hand side. See [`Rhs`].
    pub rhs: Rhs,
    /// WHICH strategy builds the right-hand side (§2.5.10). `None` is
    /// [`Compiler::Native`], and `native` is STRICT: a rule it cannot express
    /// is a construction error, not a quiet demotion.
    ///
    /// `Option` rather than a bare `Compiler` so that "the caller said
    /// nothing" and "the caller asked for native" stay distinguishable in a
    /// `ProblemOptions` a test prints; they build identically.
    pub compiler: Option<Compiler>,

    /// Opt in to the automatic projection-pushdown desugar. Turning this on
    /// also turns the build pipeline on.
    #[cfg(not(target_arch = "wasm32"))]
    pub pushdown_rewrite: bool,
    /// Run-time data providers, keyed by the forcing VARIABLE name they feed.
    ///
    /// CONST providers are materialized ONCE, here, at construction — the
    /// "gated fetch" §2.5.2 puts in the build half. DISCRETE providers
    /// contribute the refresh anchors [`solve`] segments the integration on.
    #[cfg(not(target_arch = "wasm32"))]
    pub providers: HashMap<String, Box<dyn crate::provider::CadenceProvider>>,
    /// Build-time array providers consumed by the pushdown / gated-fetch
    /// pipeline, keyed by provider key. Supplying any turns the build pipeline
    /// on.
    #[cfg(not(target_arch = "wasm32"))]
    pub build_providers: Vec<(String, Box<dyn crate::prepare::PrepareProvider>)>,
    /// Constant arrays injected into the build-time evaluation scope.
    #[cfg(not(target_arch = "wasm32"))]
    pub const_arrays: HashMap<String, ArrayD<f64>>,
    /// Force the deterministic build pipeline on (it is on automatically when
    /// `pushdown_rewrite`, `build_providers` or `const_arrays` is set).
    #[cfg(not(target_arch = "wasm32"))]
    pub build_pipeline: bool,
    /// Split a gated provider's pre-sliced fetch into requests of at most this
    /// many native indices along the gated axis.
    #[cfg(not(target_arch = "wasm32"))]
    pub gated_fetch_batch: Option<usize>,
    /// Print per-phase build progress lines on stdout.
    #[cfg(not(target_arch = "wasm32"))]
    pub verbose: bool,
    /// Build-time progress observer. **Extension seam, not stable API** — a
    /// binding MAY expose its own build observability (§2.5.2).
    #[cfg(not(target_arch = "wasm32"))]
    pub progress: Option<crate::prepare::PrepareProgressFn>,

    /// Collect the array runtime's named build-time products so
    /// [`observed_field`] can read them back.
    ///
    /// This is the **construction-time build-observability seam** that replaced
    /// the old `simulate_with_inspection` / `observed_field(prep, inspection,
    /// name)` threading: the caller asks for observability when it builds,
    /// and reads it back off the EsmProblem, instead of carrying a
    /// `BuildInspection` through the run.
    pub inspect: bool,

    /// Time at which CONST providers are sampled. Defaults to `tspan.0`.
    pub sample_time: Option<f64>,

    /// Leave every shaped parameter that has a value and is not fed from
    /// outside the model (no `update`) as a run-time input of
    /// [`crate::rhs_program`]: an [`crate::SlotKind::ArrayParameter`] forcing
    /// slot, seeded with that value, instead of build-time data (API_SPEC
    /// §5.8). The value is the document's (inline array data, or one scalar
    /// broadcast over its grid) or the caller's [`Self::const_arrays`] entry,
    /// which takes precedence. The Problem itself
    /// reads it from its forcing buffer, where construction serves that
    /// value, so `solve` integrates the same numbers either way. Off by
    /// default.
    pub runtime_array_parameters: bool,
}

impl std::fmt::Debug for ProblemOptions {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        let mut d = f.debug_struct("ProblemOptions");
        d.field("p", &self.p)
            .field("u0", &self.u0)
            .field("model_name", &self.model_name)
            .field("metaparameters", &self.metaparameters)
            .field("base_path", &self.base_path)
            .field("callbacks", &self.callbacks)
            .field("rhs", &self.rhs)
            .field("compiler", &self.compiler)
            .field("inspect", &self.inspect)
            .field("sample_time", &self.sample_time)
            .field("runtime_array_parameters", &self.runtime_array_parameters);
        #[cfg(not(target_arch = "wasm32"))]
        d.field("pushdown_rewrite", &self.pushdown_rewrite)
            .field("providers", &self.providers.keys().collect::<Vec<_>>())
            .field(
                "build_providers",
                &self
                    .build_providers
                    .iter()
                    .map(|(k, _)| k)
                    .collect::<Vec<_>>(),
            )
            .field(
                "const_arrays",
                &self.const_arrays.keys().collect::<Vec<_>>(),
            )
            .field("build_pipeline", &self.build_pipeline)
            .field("gated_fetch_batch", &self.gated_fetch_batch)
            .field("verbose", &self.verbose)
            // Hand-written because the observer is a trait object: it cannot
            // derive `Debug`, and a `ProblemOptions` that no longer prints
            // would be a regression for every `{:?}` on a build error path.
            .field(
                "progress",
                &self
                    .progress
                    .as_ref()
                    .map(|_| "<observer>")
                    .unwrap_or("None"),
            );
        d.finish()
    }
}

// =============================================================================
// The EsmProblem
// =============================================================================

/// The compiled right-hand side a [`EsmProblem`] integrates, if any.
pub(crate) enum Backend {
    /// The array / spatial runtime — the tape under `native` / `xla`, its
    /// per-cell oracle under `interpreter` (API_SPEC §5.8), for every document.
    Array(Rc<ArrayCompiled>),
    /// No right-hand side: the document declares no differential equations, or
    /// the caller asked for [`Rhs::Never`]. Carries the reason.
    Static(String),
}

/// The build-time products construction materialized.
#[derive(Default)]
pub(crate) struct BuildProducts {
    /// Observed name → build-time field (dependency-complete).
    pub fields: HashMap<String, ArrayD<f64>>,
    /// Value-invention producer id (`from_faq`) → sorted-distinct member ids.
    pub members: HashMap<String, Vec<i64>>,
    /// Value-invention producer id → derived index-set extent.
    pub extents: HashMap<String, i64>,
    /// Provider keys that were deferred + fetched pre-sliced (sorted).
    pub gated_provider_keys: Vec<String>,
    /// Parameters baked into the build. Substituting one of these needs a
    /// rebuild, so [`remake`] refuses rather than lying.
    pub baked_parameters: Vec<String>,
    /// Why a state-free document's observed graph could not be evaluated at
    /// construction — a fail-closed evaluation fault such as
    /// `table_lookup_out_of_bounds` — so [`observed_field`] can name the cause
    /// instead of only the name it cannot answer.
    pub static_eval_error: Option<String>,
}

/// The raw document a Problem carries. A document read from a file and built
/// without rewriting keeps the file's text and parses it on the first
/// [`EsmProblem::document`] call, so the build hands its one parsed copy to
/// the loader instead of cloning it.
pub(crate) struct RawDocument {
    text: Option<String>,
    value: std::cell::OnceCell<JsonValue>,
}

impl RawDocument {
    fn parsed(value: JsonValue) -> Self {
        RawDocument {
            text: None,
            value: std::cell::OnceCell::from(value),
        }
    }

    fn from_text(text: String) -> Self {
        RawDocument {
            text: Some(text),
            value: std::cell::OnceCell::new(),
        }
    }

    fn get(&self) -> &JsonValue {
        self.value.get_or_init(|| {
            // The text parsed once already, when the Problem was built.
            self.text
                .as_deref()
                .and_then(|t| serde_json::from_str(t).ok())
                .unwrap_or(JsonValue::Null)
        })
    }
}

/// A simulation problem: a document, an interval, and the bindings that fix
/// the document rather than the run.
///
/// Built once by [`esm_problem`]; run by [`solve`] (or stepped by [`init`]).
/// Re-parameterized without rebuilding by [`remake`].
pub struct EsmProblem {
    /// The (possibly rewritten) raw document.
    pub(crate) doc: Rc<RawDocument>,
    /// The name of the model this EsmProblem was built from, when one was
    /// selected.
    pub(crate) model_name: Option<String>,
    /// The working precision the document declared (`domain.element_type`,
    /// esm-spec §11.3), captured at construction. Every run entry re-arms it
    /// (`precision::enter`) so a problem carries its own precision rather than
    /// depending on what the calling thread last evaluated.
    pub(crate) precision: precision::Env,
    /// The integration interval.
    pub(crate) tspan: (f64, f64),
    /// Parameter bindings.
    pub(crate) p: HashMap<String, f64>,
    /// Initial-state bindings.
    pub(crate) u0: HashMap<String, f64>,
    /// The compiled right-hand side. `Rc` so [`remake`] shares it rather than
    /// recompiling.
    pub(crate) backend: Rc<Backend>,
    /// Build-time products, shared with every [`remake`] descendant.
    pub(crate) build: Rc<BuildProducts>,
    /// Build observability, filled at construction when
    /// [`ProblemOptions::inspect`] is set and topped up at the first solve.
    pub(crate) inspection: std::cell::RefCell<BuildInspection>,
    /// Whether the caller asked for build observability.
    pub(crate) inspect: bool,
    /// The compiler that built this Problem's right-hand side (§2.5.10).
    pub(crate) compiler: Compiler,
    /// Where every rule landed, computed at CONSTRUCTION — which is also where
    /// a strict `native` refuses, so the record and the refusal are read off
    /// the same tape build.
    pub(crate) compiler_report: CompilerReport,
    /// Callbacks declared on the EsmProblem (§2.5.4).
    pub(crate) callbacks: CallbackSet,
    /// The document's §2.2 `solver` block, captured at construction.
    ///
    /// Its OWN field rather than a read-back out of `doc`, because `doc` is
    /// `JsonValue::Null` whenever the problem was built from a typed
    /// [`EsmFile`] and the build pipeline was not requested — which made the
    /// §2.2.2 chain dead on exactly the paths a library caller uses. `None` when
    /// the document declares no block, or when the problem was built from a bare
    /// [`FlattenedSystem`], which carries no document at all.
    pub(crate) solver: Option<crate::types::Solver>,
    /// Bound run-time providers, already CONST-materialized.
    #[cfg(not(target_arch = "wasm32"))]
    pub(crate) refresh: Option<std::cell::RefCell<crate::provider::RefreshExecutor>>,
    /// Forcing variables fed by a DISCRETE provider.
    #[cfg(not(target_arch = "wasm32"))]
    pub(crate) discrete_forcing: std::collections::HashSet<String>,
    /// Refresh anchors strictly inside `tspan`.
    #[cfg(not(target_arch = "wasm32"))]
    pub(crate) refresh_boundaries: Vec<f64>,
    /// The typed input the right-hand side was compiled from, kept so
    /// [`crate::range_rhs::rhs_program`] can flatten it once and select terms
    /// by provenance (API_SPEC §5.8). `None` for a problem with no array
    /// backend.
    pub(crate) program_source: Option<Rc<ProgramSource>>,
    /// [`ProblemOptions::runtime_array_parameters`], for `rhs_program`.
    pub(crate) runtime_array_parameters: bool,
}

/// What [`EsmProblem::program_source`] holds: the prepared typed document, or
/// the caller's flattened system.
pub(crate) enum ProgramSource {
    File(Box<EsmFile>),
    Flattened(Box<FlattenedSystem>),
}

impl std::fmt::Debug for EsmProblem {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("EsmProblem")
            .field("model_name", &self.model_name)
            .field("tspan", &self.tspan)
            .field("p", &self.p.len())
            .field("u0", &self.u0.len())
            .field("backend", &self.backend_kind())
            .field("callbacks", &self.callbacks.names())
            .finish()
    }
}

impl EsmProblem {
    /// The document's §2.2 `solver` block, or `None` when it declares none.
    ///
    /// The middle level of the esm-spec §2.2.2 tolerance chain, captured at
    /// construction so it holds however the problem was built.
    pub fn solver(&self) -> Option<&crate::types::Solver> {
        self.solver.as_ref()
    }

    /// The integration interval.
    pub fn tspan(&self) -> (f64, f64) {
        self.tspan
    }

    /// The parameter bindings this EsmProblem carries.
    pub fn p(&self) -> &HashMap<String, f64> {
        &self.p
    }

    /// The initial-state bindings this EsmProblem carries.
    pub fn u0(&self) -> &HashMap<String, f64> {
        &self.u0
    }

    /// The model this EsmProblem was built from, when one was selected.
    pub fn model_name(&self) -> Option<&str> {
        self.model_name.as_deref()
    }

    /// The (possibly rewritten) raw document.
    pub fn document(&self) -> &JsonValue {
        self.doc.get()
    }

    /// The compiler that built this Problem's right-hand side (API_SPEC §5.8).
    ///
    /// The vocabulary value the build ACTUALLY used — which, because a refused
    /// compiler is never answered by substituting another, is always the one
    /// the caller asked for.
    pub fn compiler(&self) -> Compiler {
        self.compiler
    }

    /// Where every rule of this Problem landed (API_SPEC §5.8).
    ///
    /// Informative even under a strict [`Compiler::Native`], where every rule
    /// is taped by construction: the interesting question there is the cadence
    /// tier each rule runs at, which is what makes a cost readable.
    pub fn compiler_report(&self) -> &CompilerReport {
        &self.compiler_report
    }

    /// `"array"` or `"static"`.
    ///
    /// Kept for compatibility beside [`Self::compiler`], which is the stable
    /// answer to "what built this". Under both `native` and `interpreter` a
    /// dynamic Problem is `"array"`, so this does not distinguish the two;
    /// `compiler()` does.
    pub fn backend_kind(&self) -> &'static str {
        match &*self.backend {
            Backend::Array(_) => "array",
            Backend::Static(_) => "static",
        }
    }

    /// The working precision this problem evaluates in — `domain.element_type`
    /// (esm-spec §11.3), [`Precision::Float64`] unless the document said
    /// `"Float32"`.
    pub fn precision(&self) -> Precision {
        self.precision.document
    }

    /// The compiled array model behind this Problem's right-hand side, or
    /// `None` for a static Problem.
    ///
    /// Read-only access for measurement tooling (the scaling conformance
    /// tier's adapter and the wasm suite), which drive the same compiled
    /// artifact through the `debug_*` entry points of [`ArrayCompiled`]
    /// rather than through a solve.
    #[doc(hidden)]
    pub fn debug_array_compiled(&self) -> Option<Rc<ArrayCompiled>> {
        match &*self.backend {
            Backend::Array(c) => Some(Rc::clone(c)),
            Backend::Static(_) => None,
        }
    }

    /// Whether this EsmProblem has a right-hand side to integrate.
    pub fn is_dynamic(&self) -> bool {
        !matches!(&*self.backend, Backend::Static(_))
    }

    /// The state-variable names of the compiled right-hand side, in the
    /// flattened state-vector order. Empty for a static EsmProblem.
    pub fn state_variable_names(&self) -> Vec<String> {
        match &*self.backend {
            Backend::Array(c) => c.qualified_state_names().to_vec(),
            Backend::Static(_) => Vec::new(),
        }
    }

    /// The parameter names of the compiled right-hand side. Empty for a static
    /// EsmProblem.
    pub fn parameter_names(&self) -> Vec<String> {
        match &*self.backend {
            Backend::Array(c) => qualify_array_names(c, c.parameter_names()),
            Backend::Static(_) => Vec::new(),
        }
    }

    /// The value-invention members this EsmProblem's build materialized.
    pub fn members(&self) -> &HashMap<String, Vec<i64>> {
        &self.build.members
    }

    /// The derived index-set extents this EsmProblem's build materialized.
    pub fn extents(&self) -> &HashMap<String, i64> {
        &self.build.extents
    }

    /// Provider keys that were deferred and fetched pre-sliced.
    pub fn gated_provider_keys(&self) -> &[String] {
        &self.build.gated_provider_keys
    }

    /// Take the build-observability record this EsmProblem collected.
    ///
    /// Only populated when [`ProblemOptions::inspect`] was set. This is the
    /// construction-time seam that replaced threading a `&mut BuildInspection`
    /// through the run.
    pub fn take_inspection(&self) -> BuildInspection {
        std::mem::take(&mut self.inspection.borrow_mut())
    }

    /// Re-arm the build-observability record to the state CONSTRUCTION leaves
    /// it in (empty).
    ///
    /// [`esm_problem`] initialises `inspection` empty and only [`solve`] fills
    /// it — and only on the array backend, which OVERWRITES rather than merges;
    /// [`Self::take_inspection`] DRAINS it. So a caller that solves one problem
    /// more than once (the inline-test runner, which memoises the build across
    /// consecutive tests that share it) must re-arm between solves, or the
    /// second read would see what the first left behind instead of what a
    /// rebuild would have produced. Idempotent, and a no-op on a problem that
    /// has never been read.
    pub(crate) fn reset_inspection(&self) {
        *self.inspection.borrow_mut() = BuildInspection::default();
    }

    /// Every build-time field this EsmProblem's construction materialized, keyed
    /// by observed name.
    pub fn observed_fields(&self) -> &HashMap<String, ArrayD<f64>> {
        &self.build.fields
    }

    /// The names of the OBSERVED VARIABLES [`observed_trajectory`] can report a
    /// trajectory for.
    ///
    /// Disjoint from [`Self::observed_field_names`] and answering a different
    /// question: those are constants the BUILD materialized, these vary along a
    /// solution. Empty on the static backend, which has no run to vary along.
    pub fn observed_variable_names(&self) -> Vec<String> {
        match &*self.backend {
            Backend::Array(c) => qualify_array_names(c, &c.observed_variable_names()),
            Backend::Static(_) => Vec::new(),
        }
    }

    /// The names of every build-time field [`observed_field`] can return,
    /// component-qualified — the spelling that resolves whatever the problem's
    /// component count, and the one to report to an author whose bare name was
    /// refused.
    pub fn observed_field_names(&self) -> Vec<String> {
        let model = self.model_name.as_deref().unwrap_or("");
        let mut names: Vec<String> = self
            .build
            .fields
            .keys()
            .chain(self.inspection.borrow().setup_arrays.keys())
            .map(|k| qualify(model, k))
            .collect();
        names.sort();
        names.dedup();
        names
    }
}

/// The callbacks declared on `prob` (§2.5.4).
///
/// Stable API in every simulation-capable binding, and for a specific reason:
/// a `callback` argument to [`solve`] REPLACES this set entirely, so without a
/// way to read it back, a EsmProblem-level callback would be impossible to extend.
/// See [`compose`].
pub fn callbacks(prob: &EsmProblem) -> &CallbackSet {
    &prob.callbacks
}

/// `key` qualified by the component path `model`, idempotently.
///
/// The two field maps an [`EsmProblem`] carries are keyed differently. The
/// build pipeline keys a field by the SELECTED MODEL's own variable name
/// (`E_PM25`, or `North.u` for a mounted subsystem), because it evaluates one
/// model's observed graph rather than the flattened document's; the state-free
/// static evaluation keys it by the FLATTENED name (`ISRM.E_PM25`), because
/// it evaluates the flattened system. Qualifying the first spelling and leaving
/// the second alone puts both into the one namespace Julia and Python key
/// their build-time fields by, so `observed_field` answers the same spellings
/// in all three bindings.
///
/// A key EQUAL to `model` is a variable named like its model, not a component
/// path, so it is qualified too: `fuel_mce` in model `fuel_mce` is
/// `fuel_mce.fuel_mce`, the spelling the flattened name has (API_SPEC §5.8)
/// and the one Julia and Python report.
pub(crate) fn qualify(model: &str, key: &str) -> String {
    let already = key
        .strip_prefix(model)
        .is_some_and(|rest| rest.starts_with('.'));
    if model.is_empty() || already {
        key.to_string()
    } else {
        format!("{model}.{key}")
    }
}

/// `names` under the array runtime's single-model namespace, when that build
/// left them bare.
///
/// Presentation only: the compiled artifact keeps its own spelling, and the
/// override-key canonicalization (`canonicalize_override_keys`) already
/// resolves a `<namespace>.`-prefixed key against a bare slot, so a caller may
/// hand back exactly the names it was given.
///
/// This exists because the array runtime names its slots differently on its two
/// build routes: the single-model build leaves them bare, the flattened one
/// qualifies them. Without it the SAME document would report `M.y` or `y`
/// depending on a routing decision the caller cannot see, which is the failure
/// §2.5.10 exists to prevent.
fn qualify_array_names(
    compiled: &crate::simulate_array::ArrayCompiled,
    names: &[String],
) -> Vec<String> {
    let Some(ns) = compiled.namespace() else {
        return names.to_vec();
    };
    names.iter().map(|n| qualify(ns, n)).collect()
}

/// The component that owns a qualified name — everything before the final
/// segment, or `""` when the name is unqualified.
fn component_of(qualified: &str) -> &str {
    qualified.rsplit_once('.').map_or("", |(head, _)| head)
}

/// The distinct components owning `prob`'s build-time fields (API_SPEC §5.8).
///
/// A bare name is only resolvable when this holds exactly one: with two
/// mounted components a bare `u` designates `North.u` and `South.u` equally,
/// and answering with either is worse than refusing.
fn field_components(prob: &EsmProblem) -> std::collections::BTreeSet<String> {
    let model = prob.model_name.as_deref().unwrap_or("");
    let mut out = std::collections::BTreeSet::new();
    for k in prob.build.fields.keys() {
        out.insert(component_of(&qualify(model, k)).to_string());
    }
    for k in prob.inspection.borrow().setup_arrays.keys() {
        out.insert(component_of(&qualify(model, k)).to_string());
    }
    out
}

/// Outcome of applying §5.8's precedence to one namespace of stored keys.
enum NameResolution {
    /// The stored key `name` designates.
    Key(String),
    /// `name` is a bare name whose tail matches these keys (sorted), but the
    /// gate refused to bind one: a multi-component problem, or a second tail
    /// match within ONE component — where two keys cannot share a tail, so a
    /// second match means the component set was miscounted; refuse rather
    /// than pick.
    Ambiguous(Vec<String>),
    /// No rule matched.
    Miss,
}

/// The §5.8 name-resolution rule, shared by [`observed_field`] and the
/// observed-trajectory resolver, in precedence order:
///
/// 1. **Exact hit** — `name` is a stored key.
/// 2. **Component-qualified hit** — `name` is a stored key's [`qualify`]d
///    spelling (`ISRM.E_PM25` for the pipeline's `E_PM25`). Done as a strip
///    rather than a scan: `name` is the qualified spelling of key `k` exactly
///    when it is `model` + `.` + `k`. With an empty `model` the inner strip
///    fails and this is a no-op.
/// 3. **Bare name** — `name` carries no `.` and the problem has exactly ONE
///    component; it then resolves to the unique key with that tail.
///
/// A bare name whose tail matches but which rule 3's gate refuses comes back
/// as [`NameResolution::Ambiguous`] carrying every candidate, so a caller can
/// put the remedy in its diagnostic — the author has to qualify the name, and
/// cannot know which spellings exist without being told.
fn resolve_observed_key<'a>(
    keys: impl Iterator<Item = &'a String> + Clone,
    model: &str,
    single_component: bool,
    name: &str,
) -> NameResolution {
    if keys.clone().any(|k| k == name) {
        return NameResolution::Key(name.to_string()); // rule 1
    }
    if let Some(rest) = name.strip_prefix(model).and_then(|r| r.strip_prefix('.'))
        && keys.clone().any(|k| k == rest)
    {
        return NameResolution::Key(rest.to_string()); // rule 2
    }
    if name.contains('.') {
        return NameResolution::Miss; // rule 3 is for bare names only
    }
    let mut matches: Vec<&String> = keys
        .filter(|k| k.rsplit('.').next() == Some(name))
        .collect();
    matches.sort();
    if single_component && matches.len() == 1 {
        return NameResolution::Key(matches[0].clone());
    }
    if matches.is_empty() {
        NameResolution::Miss
    } else {
        NameResolution::Ambiguous(matches.into_iter().cloned().collect())
    }
}

/// The build-time field named `name`, from `prob`'s construction.
///
/// Resolution is the cross-binding rule of API_SPEC §5.8
/// ([`resolve_observed_key`]), in precedence order:
///
/// 1. **Exact hit** — `name` is a stored key.
/// 2. **Component-qualified hit** — `name` is a stored key's [`qualify`]d
///    spelling (`ISRM.E_PM25` for the pipeline's `E_PM25`).
/// 3. **Bare name** — `name` carries no `.`, and the problem has exactly ONE
///    component; it then resolves to the unique field with that tail.
///
/// A bare name against a MULTI-component problem is refused rather than bound
/// to an arbitrary candidate. Reads the build pipeline's fields, the state-free
/// static evaluation's fields, and — when [`ProblemOptions::inspect`] was set —
/// the array runtime's materialized setup arrays. One arity, `(prob, name)`,
/// in every binding.
pub fn observed_field(prob: &EsmProblem, name: &str) -> Result<ArrayD<f64>, SimulateError> {
    // Re-arm the document's working precision for the duration of this call
    // (`domain.element_type`, esm-spec §11.3); a no-op for a Float64 document.
    let _precision_guard = prob.precision.enter();
    // A name an `operator_compose` renaming match DELETED addresses a field that
    // MOVED (issue #230). The merge only ever removes a spelling, so a hit here
    // can never shadow a live field.
    let name = problem_merged_renames(prob)
        .get(name)
        .map(String::as_str)
        .unwrap_or(name);
    let model = prob.model_name.as_deref().unwrap_or("");
    let components = field_components(prob);
    let single = components.len() == 1;

    if let NameResolution::Key(k) =
        resolve_observed_key(prob.build.fields.keys(), model, single, name)
    {
        return Ok(prob.build.fields[&k].clone());
    }
    let inspection = prob.inspection.borrow();
    if let NameResolution::Key(k) =
        resolve_observed_key(inspection.setup_arrays.keys(), model, single, name)
    {
        return Ok(inspection.setup_arrays[&k].clone());
    }
    // A bare name that WOULD have resolved but for the component gate gets the
    // remedy in the diagnostic: the author has to qualify it, and cannot know
    // which spellings exist without being told. (`Ambiguous` under a SINGLE
    // component — the miscounted-set refusal — keeps the generic error below,
    // as it always has.)
    if !single
        && let NameResolution::Ambiguous(matches) = resolve_observed_key(
            prob.build
                .fields
                .keys()
                .chain(inspection.setup_arrays.keys()),
            model,
            single,
            name,
        )
    {
        let mut cands: Vec<String> = matches.iter().map(|k| qualify(model, k)).collect();
        cands.sort();
        cands.dedup();
        return Err(SimulateError::Compile(
            crate::compile_error::CompileError::build_err(format!(
                "observed_field: '{name}' is a bare name and this EsmProblem has {} \
                     components ({}); qualify it as one of: {}",
                components.len(),
                components
                    .iter()
                    .filter(|c| !c.is_empty())
                    .cloned()
                    .collect::<Vec<_>>()
                    .join(", "),
                cands.join(", ")
            )),
        ));
    }
    let cause = prob
        .build
        .static_eval_error
        .as_deref()
        .map(|e| format!("; evaluating the state-free observed graph failed: {e}"))
        .unwrap_or_default();
    Err(SimulateError::Compile(
        crate::compile_error::CompileError::build_err(format!(
            "observed_field: '{name}' is not a build-time-evaluable observed of this \
             EsmProblem{cause}"
        )),
    ))
}

/// The trajectory of the observed variable `name` over `sol`'s output grid.
///
/// The solution-aware half of [`observed_field`], and the answer to a question
/// neither argument can answer alone. An observed is a pure function of
/// `(state, params, t)`: `prob` holds the function — the topo-sorted observed
/// graph and the parameter bindings — and `sol` holds the arguments. That is
/// why this takes both, and why [`observed_field`] can only ever report fields
/// the BUILD materialized. A [`Solution`] carries state rows only, so an
/// observed of a model that does integrate has no other way out.
///
/// Named `observed_trajectory` rather than being a second arity of
/// `observed_field` because the return RANK differs: a field is the shape the
/// document declares, a trajectory is that shape with a time axis added. Two
/// ranks under one name is a contract a caller cannot read off the call. (A
/// binding that can overload — Julia, Python — may still spell it as an extra
/// arity of `observed_field`; API_SPEC §5.8 records the transliteration.)
///
/// Name resolution is §5.8's rule, the same one [`observed_field`] applies:
/// exact hit, then the component-qualified spelling, then a bare name only when
/// the problem has exactly one component.
///
/// A static problem has no trajectory at all and wants [`observed_field`].
#[cfg(feature = "solve")]
pub fn observed_trajectory(
    prob: &EsmProblem,
    sol: &Solution,
    name: &str,
) -> Result<Vec<f64>, SimulateError> {
    // Re-arm the document's working precision for the duration of this call
    // (`domain.element_type`, esm-spec §11.3); a no-op for a Float64 document.
    let _precision_guard = prob.precision.enter();
    let one = [name.to_string()];
    // The bulk form omits what it cannot resolve; the singular one must not,
    // so an empty result becomes the diagnostic the resolver would have given.
    match observed_trajectories(prob, sol, &one)?.pop() {
        Some((_, values)) => Ok(values),
        None => {
            let declared: Vec<String> = match &*prob.backend {
                // The name resolved to nothing, and the diagnostic has to be
                // built against the runtime's own declared names.
                Backend::Array(c) => c.observed_variable_names(),
                Backend::Static(_) => Vec::new(),
            };
            let model = prob.model_name.as_deref().unwrap_or("");
            Err(
                resolve_observed_name(&declared, model, components_of(&declared, model) == 1, name)
                    .expect_err("resolution failed, or the bulk form would have answered"),
            )
        }
    }
}

/// [`observed_trajectory`] for several names, in ONE pass over the output grid.
///
/// The graph is walked once per output time however many names are asked for,
/// so a caller wanting five observeds should ask once rather than five times.
///
/// **Tolerant where the singular form is strict**, and returns `(name, values)`
/// pairs so a caller can tell which is which. A name that is not an observed
/// variable — most often a STATE, which the caller already has in `sol` — is
/// omitted rather than failing the whole call. That is what a host asking "give
/// me whichever of these are observed" needs: a test harness reading a model's
/// authored assertions knows the variable names but not which kind each is, and
/// one state in the list must not cost it the other four answers.
///
/// The returned names are the spellings that were ASKED FOR, not the resolved
/// ones, so a caller can key its own lookup by them.
#[cfg(feature = "solve")]
pub fn observed_trajectories(
    prob: &EsmProblem,
    sol: &Solution,
    names: &[String],
) -> Result<Vec<(String, Vec<f64>)>, SimulateError> {
    // Re-arm the document's working precision for the duration of this call
    // (`domain.element_type`, esm-spec §11.3); a no-op for a Float64 document.
    let _precision_guard = prob.precision.enter();
    match &*prob.backend {
        // The array runtime has ALREADY materialized its observeds at every
        // saved time — `append_observed_trajectories` appends each 0-D one to
        // the solution as a row beside the states, in the same cell-key
        // spelling — so the answer is a lookup in `sol`, not a second
        // evaluation.
        Backend::Array(c) => array_observed_trajectories(prob, c, sol, names),
        Backend::Static(reason) => Err(SimulateError::NotDynamic {
            details: format!("{reason}; a static document's results are read with observed_field"),
        }),
    }
}

/// [`observed_trajectories`] on the array runtime: read the rows that runtime
/// already appended to `sol`.
///
/// **Not a second evaluation, on purpose.** The array runtime materializes its
/// observeds at every output node during the solve (under `native`, from the
/// tape) and appends each 0-D one to the solution. Re-deriving them here would
/// be exactly the off-tape output-time evaluation `esm-libraries-spec.md`
/// §2.5.10 puts under the refusal — and it would also be a second answer to a
/// question that already has one.
///
/// **Rank.** An ARRAY-valued observed has no scalar trajectory, so it is not
/// answered here. It is reachable — `SolveOptions::output_observed` appends it
/// as one row per cell, named `base[i,j]` — and the diagnostic says so, because
/// a caller who asked for a name that exists deserves the remedy rather than
/// "not an observed variable".
///
/// Tolerant: a name that resolves to no observed is omitted rather than failing
/// the call.
#[cfg(feature = "solve")]
fn array_observed_trajectories(
    prob: &EsmProblem,
    compiled: &crate::simulate_array::ArrayCompiled,
    sol: &Solution,
    names: &[String],
) -> Result<Vec<(String, Vec<f64>)>, SimulateError> {
    let declared = compiled.observed_variable_names();
    // The COMPILED model's namespace, not `prob.model_name`: a problem built
    // from a typed single-model document carries no `model_name` (the caller
    // named none), while the artifact knows the `models` key it consumed. With
    // the empty string here, §5.8's rule-2 strip never fires and a caller's
    // component-qualified spelling resolves to nothing.
    let model = compiled
        .namespace()
        .or(prob.model_name.as_deref())
        .unwrap_or("");
    let single = components_of(&declared, model) == 1;
    let mut out = Vec::new();
    for name in names {
        let Ok(resolved) = resolve_observed_name(&declared, model, single, name) else {
            continue;
        };
        // The solution's rows carry the qualified spelling (see
        // `qualify_array_names`); the compiled model's rule names are bare.
        let row_name = match compiled.namespace() {
            Some(ns) => qualify(ns, &resolved),
            None => resolved.clone(),
        };
        match sol.state_variable_names.iter().position(|r| *r == row_name) {
            Some(i) => out.push((name.clone(), sol.state[i].clone())),
            None => {
                return Err(SimulateError::Compile(
                    crate::compile_error::CompileError::build_err(format!(
                        "observed_trajectory: '{name}' is an observed of this EsmProblem but \
                         carries no scalar trajectory — it is array-valued, and the array \
                         runtime appends one row per cell only for an observed the solve was \
                         asked for. Name it in `SolveOptions::output_observed` and read the \
                         `{resolved}[i,j]` rows of the solution."
                    )),
                ));
            }
        }
    }
    Ok(out)
}

/// §5.8's precedence ([`resolve_observed_key`]), against a flat list of
/// declared observed names.
#[cfg(feature = "solve")]
fn resolve_observed_name(
    declared: &[String],
    model: &str,
    single_component: bool,
    name: &str,
) -> Result<String, SimulateError> {
    match resolve_observed_key(declared.iter(), model, single_component, name) {
        NameResolution::Key(k) => Ok(k),
        // A bare name that WOULD have resolved but for the component count gets
        // the remedy in the diagnostic — the author cannot qualify it without
        // being told which spellings exist.
        NameResolution::Ambiguous(matches) => Err(SimulateError::Compile(
            crate::compile_error::CompileError::build_err(format!(
                "observed_trajectory: '{name}' is a bare name and this EsmProblem has {}                          components; qualify it as one of: {}",
                components_of(declared, model),
                matches
                    .iter()
                    .map(|k| qualify(model, k))
                    .collect::<Vec<_>>()
                    .join(", ")
            )),
        )),
        NameResolution::Miss => Err(SimulateError::Compile(
            crate::compile_error::CompileError::build_err(format!(
                "observed_trajectory: '{name}' is not an observed variable of this EsmProblem"
            )),
        )),
    }
}

#[cfg(feature = "solve")]
fn components_of(declared: &[String], model: &str) -> usize {
    declared
        .iter()
        .map(|k| component_of(&qualify(model, k)).to_string())
        .collect::<std::collections::BTreeSet<_>>()
        .len()
}

// =============================================================================
// remake
// =============================================================================

/// The substitutions [`remake`] applies. Every field is optional; an omitted
/// field is inherited unchanged.
#[derive(Debug, Clone, Default)]
pub struct Remake {
    /// Replacement parameter bindings, merged over the EsmProblem's.
    pub p: HashMap<String, f64>,
    /// Replacement initial-state bindings, merged over the EsmProblem's.
    pub u0: HashMap<String, f64>,
    /// A different integration interval.
    pub tspan: Option<(f64, f64)>,
    /// Replacement callbacks. `None` inherits the EsmProblem's set.
    pub callbacks: Option<CallbackSet>,
}

/// A NEW EsmProblem with `changes` applied and everything else shared (§2.5.5).
///
/// Does not mutate `prob`, and does not redo the parts of construction the
/// substitution cannot have invalidated: the compiled right-hand side, the
/// build-time fields and the materialized provider data are shared by `Rc`, not
/// rebuilt. A changed parameter value does not re-fetch provider data or
/// recompile.
///
/// **Refusal is deliberate.** A substitution the EsmProblem cannot honour without
/// a rebuild raises [`SimulateError::UnsubstitutableBinding`], naming the
/// binding and the class that makes it un-substitutable, rather than silently
/// rebuilding or silently ignoring it. Two classes refuse: a parameter that was
/// baked into the build (it is a load-time constant of the compiled RHS, not a
/// solver input), and a name that is not a parameter of the compiled system at
/// all.
pub fn remake(prob: &EsmProblem, changes: &Remake) -> Result<EsmProblem, SimulateError> {
    // Re-arm the document's working precision for the duration of this call
    // (`domain.element_type`, esm-spec §11.3); a no-op for a Float64 document.
    let _precision_guard = prob.precision.enter();
    let known: std::collections::HashSet<String> = prob.parameter_names().into_iter().collect();
    let known_bare: HashMap<String, usize> = {
        let mut counts: HashMap<String, usize> = HashMap::new();
        for n in &known {
            *counts
                .entry(n.rsplit('.').next().unwrap_or(n).to_string())
                .or_insert(0) += 1;
        }
        counts
    };

    for name in changes.p.keys() {
        if prob.build.baked_parameters.iter().any(|b| b == name) {
            return Err(SimulateError::UnsubstitutableBinding {
                name: name.clone(),
                class: "baked into the build as a load-time constant — build a new EsmProblem"
                    .to_string(),
            });
        }
        if prob.is_dynamic()
            && !known.contains(name)
            && known_bare.get(name).copied().unwrap_or(0) == 0
        {
            return Err(SimulateError::UnsubstitutableBinding {
                name: name.clone(),
                class: "not a parameter of the compiled system".to_string(),
            });
        }
    }

    let mut p = prob.p.clone();
    p.extend(changes.p.iter().map(|(k, v)| (k.clone(), *v)));
    let mut u0 = prob.u0.clone();
    u0.extend(changes.u0.iter().map(|(k, v)| (k.clone(), *v)));
    let tspan = changes.tspan.unwrap_or(prob.tspan);

    #[cfg(not(target_arch = "wasm32"))]
    let boundaries: Vec<f64> = prob
        .refresh_boundaries
        .iter()
        .copied()
        .filter(|&b| b > tspan.0 && b < tspan.1)
        .collect();

    Ok(EsmProblem {
        doc: Rc::clone(&prob.doc),
        solver: prob.solver.clone(),
        model_name: prob.model_name.clone(),
        precision: prob.precision.clone(),
        tspan,
        p,
        u0,
        backend: Rc::clone(&prob.backend),
        build: Rc::clone(&prob.build),
        inspection: std::cell::RefCell::new(prob.inspection.borrow().clone()),
        inspect: prob.inspect,
        // The backend is SHARED (`Rc::clone` above), so the compiler and its
        // per-rule record carry over verbatim: a remake re-binds values, it
        // never rebuilds a right-hand side, and reporting a compiler the
        // derivative did not run would be a lie about which one produced its
        // numbers.
        compiler: prob.compiler,
        compiler_report: prob.compiler_report.clone(),
        callbacks: changes
            .callbacks
            .clone()
            .unwrap_or_else(|| prob.callbacks.clone()),
        // The materialized provider data is shared by moving the executor's
        // ownership question out of `remake`'s way: a remade EsmProblem reads the
        // SAME forcing buffer, because the provider fetch is exactly the work
        // §2.5.5 forbids redoing. A EsmProblem with providers therefore cannot be
        // remade into two live copies; the derivative borrows nothing and
        // carries no executor of its own.
        #[cfg(not(target_arch = "wasm32"))]
        refresh: None,
        #[cfg(not(target_arch = "wasm32"))]
        discrete_forcing: prob.discrete_forcing.clone(),
        #[cfg(not(target_arch = "wasm32"))]
        refresh_boundaries: boundaries,
        program_source: prob.program_source.clone(),
        runtime_array_parameters: prob.runtime_array_parameters,
    })
}

// =============================================================================
// Ensembles
// =============================================================================

/// A EsmProblem plus a per-trajectory rewrite, and the family it stands for
/// (§2.5.8) — the canonical form for parameter sweeps, Monte Carlo over
/// declared distributions, and perturbed initial conditions.
pub struct EnsembleProblem<'a> {
    prob: &'a EsmProblem,
    trajectories: usize,
    #[allow(clippy::type_complexity)]
    rewrite: Box<dyn Fn(&EsmProblem, usize) -> Result<Remake, SimulateError> + 'a>,
}

impl std::fmt::Debug for EnsembleProblem<'_> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("EnsembleProblem")
            .field("trajectories", &self.trajectories)
            .field("prob", self.prob)
            .finish()
    }
}

impl<'a> EnsembleProblem<'a> {
    /// Wrap `prob` with a rewrite applied once per trajectory index.
    ///
    /// The rewrite returns a [`Remake`], not a EsmProblem, so every trajectory
    /// goes through the same refusal rules and shares the same compiled
    /// right-hand side.
    pub fn new(
        prob: &'a EsmProblem,
        trajectories: usize,
        rewrite: impl Fn(&EsmProblem, usize) -> Result<Remake, SimulateError> + 'a,
    ) -> Self {
        Self {
            prob,
            trajectories,
            rewrite: Box::new(rewrite),
        }
    }

    /// The base EsmProblem.
    pub fn problem(&self) -> &EsmProblem {
        self.prob
    }

    /// How many trajectories the family holds.
    pub fn trajectories(&self) -> usize {
        self.trajectories
    }

    /// The EsmProblem for trajectory `i`.
    pub fn trajectory(&self, i: usize) -> Result<EsmProblem, SimulateError> {
        let changes = (self.rewrite)(self.prob, i)?;
        remake(self.prob, &changes)
    }
}

/// Solve every trajectory of `ens`, in index order (§2.5.8).
#[cfg(feature = "solve")]
pub fn solve_ensemble(
    ens: &EnsembleProblem<'_>,
    opts: &SolveOptions,
) -> Result<Vec<Solution>, SimulateError> {
    let mut out = Vec::with_capacity(ens.trajectories());
    for i in 0..ens.trajectories() {
        let prob = ens.trajectory(i)?;
        out.push(solve(&prob, opts)?);
    }
    Ok(out)
}

// =============================================================================
// Construction
// =============================================================================

/// Build a [`EsmProblem`] (§2.5.2).
///
/// Absorbs the whole deterministic-per-document pipeline — the pushdown
/// rewrite, value invention, the gated fetch of provider data, CONST-provider
/// materialization, and the compile of the right-hand side — and takes the
/// bindings that fix a *document*: parameters, initial state, data providers,
/// metaparameters, and the model to build when the file holds several.
///
/// **Does not require the solver.** With the `solve` Cargo feature off this
/// still builds; only [`solve`] and the stepping entry points need `diffsol`.
pub fn esm_problem<'a>(
    input: impl Into<ProblemInput<'a>>,
    tspan: (f64, f64),
    opts: ProblemOptions,
) -> Result<EsmProblem, SimulateError> {
    let input = input.into();
    #[allow(unused_mut)]
    let mut opts = opts;

    // ---- (1) Resolve the input to a raw document and/or a typed one. -------
    let mut owned_json: Option<JsonValue> = None;
    let mut owned_file: Option<EsmFile> = None;
    let mut flat_only: Option<&FlattenedSystem> = None;
    // `owned_json` is exactly what the parser made of the file's text, so it
    // nests within the parser's recursion limit (see stage (3)).
    let mut json_from_text = false;
    // Whether the build pipeline rewrote the document. Only a rewritten one
    // may fail the typed parse and still build (see stage (3)).
    #[cfg_attr(target_arch = "wasm32", allow(unused_mut))]
    let mut rewritten_by_pipeline = false;

    // The file's text, while `owned_json` is still exactly its parse.
    let mut source_text: Option<String> = None;
    match input {
        ProblemInput::Path(path) => {
            let text = std::fs::read_to_string(path).map_err(|e| {
                SimulateError::Compile(crate::compile_error::CompileError::build_err(format!(
                    "reading {}: {e}",
                    path.display()
                )))
            })?;
            let raw: JsonValue = serde_json::from_str(&text).map_err(|e| {
                SimulateError::Compile(crate::compile_error::CompileError::build_err(format!(
                    "parsing {}: {e}",
                    path.display()
                )))
            })?;
            owned_json = Some(raw);
            source_text = Some(text);
            json_from_text = true;
            // A relative `{ref}` or template import resolves against the
            // referencing file's directory (esm-spec §4.7), not the process's
            // working directory; a caller's explicit `base_path` still wins.
            if opts.base_path.is_none() {
                opts.base_path = Some(match path.parent() {
                    Some(dir) if !dir.as_os_str().is_empty() => dir.to_path_buf(),
                    _ => PathBuf::from("."),
                });
            }
        }
        ProblemInput::Json(v) => owned_json = Some(v.clone()),
        ProblemInput::File(f) => owned_file = Some(f.clone()),
        ProblemInput::Flattened(f) => flat_only = Some(f),
    }

    // ---- (1b) The working precision, and the guard that arms it. ----------
    // `domain.element_type` (esm-spec §11.3). This has to be read from the RAW
    // input and armed BEFORE the build pipeline, not after the typed parse:
    // value invention, the pushdown rewrite and the build-time field
    // materialization all evaluate expressions, and a constant folded in
    // binary64 at build time would disagree with the same subexpression
    // evaluated in binary32 at run time.
    let prec = element_type_of(owned_json.as_ref(), owned_file.as_ref(), flat_only)?;
    // Per-variable overrides (`ModelVariable.element_type`, esm-spec §11.3.1)
    // are read from the same raw input and armed alongside the document's
    // precision, because the ingress sites below — and the provider fetches in
    // `prepare` — fill NAMED variables and must round each to ITS declared
    // precision, not to the document's. Empty for every document that declares
    // none, which is the inert path.
    let var_precisions = variable_element_types_of(
        owned_json.as_ref(),
        owned_file.as_ref(),
        flat_only,
        prec,
        opts.model_name.as_deref(),
    )?;
    let _precision_guard = precision::enter_env(prec, std::sync::Arc::new(var_precisions));

    // Host-supplied numbers are the one class of value the evaluator does not
    // itself produce, so they are rounded ONCE here rather than on each of the
    // O(N) reads of them. Everything computed downstream is binary32 already,
    // every operation having rounded. Each is rounded to the precision of the
    // variable it binds: an override that spared a key column on the provider
    // path but not on the host-supplied one would be no exemption at all.
    if prec.is_f32() || precision::has_variable_overrides() {
        for (k, v) in opts.p.iter_mut() {
            *v = precision::of_variable(k).round(*v);
        }
        for (k, v) in opts.u0.iter_mut() {
            *v = precision::of_variable(k).round(*v);
        }
        // `const_arrays` is a host-side build-pipeline input and does not
        // exist on wasm32, where `ProblemOptions` carries no build half.
        #[cfg(not(target_arch = "wasm32"))]
        for (k, a) in opts.const_arrays.iter_mut() {
            let kp = precision::of_variable(k);
            if kp.is_f32() {
                a.mapv_inplace(|x| kp.round(x));
            }
        }
    }

    // ---- (1c) Static precision inference, on the TYPED document. ---------
    // Propagate every `ModelVariable.element_type` over the equations and mark
    // the precision boundaries (esm-spec §11.3.1, `crate::precision_infer`).
    //
    // It runs on the TYPED form and BEFORE the build pipeline, which is the
    // only point that satisfies both constraints. Typed, because that is what
    // has `expression_templates` expanded and `$ref` imports resolved, so no
    // template body escapes the pass. Before the pipeline, because the pipeline
    // MATERIALIZES a relational document's whole observed graph — for a
    // state-free document that is the entire evaluation — and a marker
    // inserted after it would arrive too late to decide anything.
    //
    // Skipped entirely, not merely a no-op walk, for a document that declares
    // no per-variable element type.
    let mut inferred = false;
    if precision::has_variable_overrides()
        && let Some(f) = owned_file.as_mut()
        && let Some(models) = f.models.as_mut()
    {
        crate::precision_infer::annotate_models(models, prec).map_err(SimulateError::Compile)?;
        inferred = true;
    }

    // ---- (1d) Lower `table_lookup` (esm-spec §9.5.3). ---------------------
    // `table_lookup` is sugar over `interp.linear` / `interp.bilinear` /
    // `index`, and NOTHING downstream evaluates it — the array runtime's
    // stage-(0) gate refuses it as `unevaluable_operator`. Lowering happens
    // here rather than in `parse::load_value` because §9.5.4 requires the
    // AUTHORED form to round-trip and this binding serializes the typed
    // document it loaded; a load-time rewrite would emit the lowered `fn`
    // tree (issue #188).
    //
    // Done for a TYPED input before the block below re-serializes it, so the
    // build pipeline's own observed-graph evaluation sees the lowered form
    // too; a raw-JSON input gets the same treatment at stage (3c), after its
    // typed parse. The pass is idempotent, so a document that goes through
    // both is lowered once.
    if let Some(f) = owned_file.as_mut() {
        crate::lower_table_lookup::lower_table_lookups(f).map_err(SimulateError::Compile)?;
    }
    // A caller-flattened system has no document, but `flatten` carries
    // `function_tables` so that this carrier can be lowered too. The lowered
    // system is a copy the build owns; the caller's keeps the authored form.
    let lowered_flat;
    if let Some(flat) = flat_only
        && let Some(lowered) = crate::lower_table_lookup::lowered_flattened_copy(flat)
            .map_err(SimulateError::Compile)?
    {
        lowered_flat = lowered;
        flat_only = Some(&lowered_flat);
    }

    // ---- (2) The deterministic build pipeline. ----------------------------
    // `mut` on wasm32 only in the sense that the pipeline that writes these is
    // native-only; the bindings themselves exist on both targets.
    let mut build = BuildProducts::default();
    #[cfg_attr(target_arch = "wasm32", allow(unused_mut))]
    let mut model_name = opts.model_name.clone();
    let compiler = opts.compiler.unwrap_or_default();
    // A compiler this build or process cannot provide is refused before
    // anything is evaluated for it: that failure says nothing about the
    // document, so no document-level refusal may pre-empt it.
    if let Some(details) = unavailable_reason(compiler) {
        return Err(SimulateError::CompilerUnavailable {
            compiler: compiler.as_str(),
            details,
        });
    }
    // What construction evaluates OUTSIDE the compiled model's rule set — the
    // build pipeline, a state-free document's observed graph — reported ahead
    // of that rule set (API_SPEC §5.8, "Every Problem reports what ran").
    #[cfg_attr(target_arch = "wasm32", allow(unused_mut))]
    let mut route_rows: Vec<CompilerRuleReport> = Vec::new();

    // What the caller supplied for the document's shaped parameters, captured
    // before the build pipeline consumes `const_arrays` and `build_providers`:
    // the compiled model binds the arrays too (stage 3d), and the missing-data
    // gate asks which names any channel serves.
    let supplied = SuppliedData::capture(&mut opts, owned_json.as_ref(), owned_file.as_ref());

    // A TYPED document is a legitimate input to the build pipeline, and it used
    // to be the one input shape that silently was not: the pipeline reads raw
    // JSON, so `ProblemInput::File` skipped it entirely — `build_providers`,
    // `const_arrays` and `pushdown_rewrite` were accepted and then DROPPED,
    // with no error, and the caller got a document evaluated against every
    // data-fed parameter's `default`. Re-serializing here routes it through the
    // same path `Json` takes; the typed form is dropped so step (3) re-parses
    // the PREPARED document rather than the authored one, exactly as the JSON
    // input does.
    //
    // `ProblemInput::Flattened` still cannot: a flattened system has no
    // document to rewrite. It is reachable only from inside this crate.
    #[cfg(not(target_arch = "wasm32"))]
    if owned_json.is_none()
        && let Some(file) = owned_file.as_ref().filter(|_| wants_build_pipeline(&opts))
    {
        owned_json = Some(serde_json::to_value(file).map_err(|e| {
            SimulateError::Compile(crate::compile_error::CompileError::build_err(format!(
                "re-serializing the typed document for the build pipeline: {e}"
            )))
        })?);
        owned_file = None;
    }

    #[cfg(not(target_arch = "wasm32"))]
    if let Some(raw) = owned_json.as_mut().filter(|_| wants_build_pipeline(&opts)) {
        {
            let prepared = crate::prepare::run_build_pipeline(raw, &mut opts)?;
            // The pipeline evaluates the observed graph through the reference
            // evaluator whatever the compiler, and under a strict one it
            // stopped at the first observed it had to walk per cell.
            if let Some(r) = &prepared.refused {
                // The walk stopped at its first cell (`StopAtFirstCell`).
                let reason = match &r.route {
                    crate::prepare::Route::PerCell(why) => {
                        format!("{why}; {}", crate::simulate_array::ONE_CELL_NOTE)
                    }
                    _ => String::new(),
                };
                return Err(SimulateError::Compile(
                    crate::compile_error::CompileError::CompilerRefusedRule {
                        compiler: compiler.as_str(),
                        kind: r.kind,
                        rule: qualify(&prepared.model_name, &r.name),
                        tier: "const",
                        reason,
                    },
                ));
            }
            // The pipeline keeps the whole-array overlay under every compiler,
            // `interpreter` included, so its rows say which observeds the
            // overlay served.
            route_rows.extend(prepared.rules.iter().map(|r| {
                let rule = qualify(&prepared.model_name, &r.name);
                match &r.route {
                    crate::prepare::Route::Relational => CompilerRuleReport {
                        rule,
                        kind: r.kind,
                        cadence: "const",
                        tier: "relational",
                        reason: None,
                    },
                    crate::prepare::Route::Overlay => {
                        outside_rule_row(compiler, rule, r.kind, None, false)
                    }
                    crate::prepare::Route::Tape => CompilerRuleReport {
                        rule,
                        kind: r.kind,
                        cadence: "const",
                        tier: "taped",
                        reason: None,
                    },
                    crate::prepare::Route::PerCell(why) => {
                        outside_rule_row(compiler, rule, r.kind, Some(why.clone()), false)
                    }
                }
            }));
            *raw = prepared.doc;
            json_from_text = false;
            source_text = None;
            rewritten_by_pipeline = true;
            model_name = Some(prepared.model_name);
            build.fields = prepared.fields;
            build.members = prepared.members;
            build.extents = prepared.extents;
            build.gated_provider_keys = prepared.gated_provider_keys;
            build.baked_parameters = opts.p.keys().cloned().collect();
            build.baked_parameters.sort();
        }
    }

    // ---- (3) Typed parse. -------------------------------------------------
    // The document's own solver block, read here when the raw document is
    // handed to the loader below rather than kept.
    let mut doc_solver_of_text: Option<Option<crate::Solver>> = None;
    if owned_file.is_none() && owned_json.is_some() {
        {
            // The prepared document goes to the loader as the value it already
            // is. Writing it out and parsing the text back gives the same
            // document -- this crate parses with `float_roundtrip`, so every
            // number reads back as the value written -- except that the parser
            // refuses nesting past its recursion limit, so only a document
            // nested that deep still takes the round trip, to get the same
            // refusal. A document parsed from a file's text and not rewritten
            // since is within that limit already.
            // The loader-API metaparameters close the AUTHORED document
            // (esm-spec §9.7.6 site 4); a pipeline-rewritten one has had them
            // closed already, and may no longer declare the names.
            let load_opts = crate::parse::LoadOptions {
                base_path: opts.base_path.clone(),
                metaparameters: if rewritten_by_pipeline {
                    BTreeMap::new()
                } else {
                    opts.metaparameters.clone()
                },
            };
            let loaded = if json_from_text && source_text.is_some() {
                // The Problem keeps the file's text, so the loader may consume
                // the parse rather than a copy of it.
                let raw = owned_json.take().expect("checked above");
                doc_solver_of_text = Some(document_solver(&raw));
                crate::parse::load_document_owned(raw, &load_opts)
            } else if let Some(raw) = owned_json.as_ref()
                && (json_from_text || nesting_within(raw, NESTING_WITHOUT_ROUND_TRIP))
            {
                crate::parse::load_document_with_options(raw, &load_opts)
            } else {
                let raw = owned_json.as_ref().expect("checked above");
                let text = serde_json::to_string(raw).map_err(|e| {
                    SimulateError::Compile(crate::compile_error::CompileError::build_err(format!(
                        "re-serializing the prepared document: {e}"
                    )))
                })?;
                crate::parse::load_string_with_options(&text, &load_opts)
            };
            match loaded {
                Ok(f) => owned_file = Some(f),
                Err(e) => {
                    // The AUTHORED document failing to load — a schema
                    // violation, an unresolved import, a removed op, a version
                    // the library does not read — is fatal under every `Rhs`
                    // mode: esm-libraries-spec §2.1a says a library MUST NOT
                    // silently accept an invalid file, and a Problem built
                    // from one would answer for a document that does not load.
                    if !rewritten_by_pipeline {
                        return Err(SimulateError::Compile(
                            crate::compile_error::CompileError::build_err(format!(
                                "loading the document: {e}"
                            )),
                        ));
                    }
                    // A document the build pipeline rewrote may no longer be a
                    // *typed* ESM document (the pushdown desugar emits engine
                    // constructs). Its build-time products are still valid, so
                    // this is only fatal when the caller wanted a solve.
                    if opts.rhs == Rhs::Always {
                        return Err(SimulateError::Compile(
                            crate::compile_error::CompileError::build_err(format!(
                                "typed parse of the prepared document: {e}"
                            )),
                        ));
                    }
                }
            }
        }
    }

    // ---- (3b) Static precision inference, for a RAW-JSON input. -----------
    // Stage (1c) already ran when the caller handed in a typed document. A raw
    // JSON input has no typed form until stage (3) — this is its first chance,
    // and the pass has to happen before stage (4) lowers the equations.
    if !inferred
        && precision::has_variable_overrides()
        && let Some(f) = owned_file.as_mut()
        && let Some(models) = f.models.as_mut()
    {
        crate::precision_infer::annotate_models(models, prec).map_err(SimulateError::Compile)?;
    }

    // ---- (3c) Lower `table_lookup`, for a RAW-JSON input. ------------------
    // The counterpart of stage (1d) for a document that had no typed form
    // until stage (3). Idempotent, so a typed input already lowered above
    // walks nothing here.
    if let Some(f) = owned_file.as_mut() {
        crate::lower_table_lookup::lower_table_lookups(f).map_err(SimulateError::Compile)?;
    }

    // ---- (3d) Caller arrays for shaped parameters. ------------------------
    // CONFORMANCE_SPEC §5.32.5: a caller's `const_arrays` entry that gives a
    // shaped parameter its value is the same binding as an inline-array
    // override, so it lands where that override does — on the parameter's
    // `default`, which the array compile lowers. The build pipeline above only
    // evaluates the observed graph; without this the compiled right-hand side
    // never saw the array and ran on the declared default instead.
    #[cfg(not(target_arch = "wasm32"))]
    if let Some(f) = owned_file.as_mut() {
        supplied.bind_arrays(f, &opts.p, model_name.as_deref());
    }
    // A scalar parameter value invention reads (a skolem bin width) is
    // build-time data: the invented set's size depends on it, so the caller's
    // `p` entry lands on its `default`, which the compile's value invention
    // reads, and the parameter is baked into the build.
    if let Some(f) = owned_file.as_mut() {
        for name in bind_value_invention_params(f, &opts.p) {
            if !build.baked_parameters.contains(&name) {
                build.baked_parameters.push(name);
            }
        }
        // A scalar parameter a join key column reads decides the join's match
        // set, which the build constructs once (CONFORMANCE_SPEC §5.5.8): it
        // is baked the same way, under every compiler.
        for name in bind_join_key_params(f, &opts.p) {
            if !build.baked_parameters.contains(&name) {
                build.baked_parameters.push(name);
            }
        }
        build.baked_parameters.sort();
    }

    // ---- (3d) A `p` pin BINDS a data-fed parameter (esm-spec §9.6.6). -----
    // The caller passed a value for a parameter the document says is read from
    // a file, so for THIS build it is not read from a file — that is the
    // documented escape hatch for running a data-fed document offline, and for
    // its own inline tests (`Test.parameter_overrides`, esm-spec §6.6, is this
    // argument by another name). Stripping the `update` here, before the
    // backend is built, is what makes the pin REACH the right-hand side: a
    // parameter that still carries a `data` update is routed to the external
    // forcing channel, where a `p` binding never lands and which the tape
    // cannot lower at all, so the pin used to change nothing and `native`
    // refused the document anyway.
    #[cfg(not(target_arch = "wasm32"))]
    if !opts.p.is_empty()
        && let Some(file) = owned_file.as_mut()
    {
        crate::data_fed::pin_data_fed_parameters(file, opts.p.keys());
    }

    // ---- (4) Compile the right-hand side. ---------------------------------
    let backend = compile_backend(
        owned_file.as_ref(),
        flat_only,
        model_name.as_deref(),
        opts.rhs,
        compiler,
        opts.runtime_array_parameters,
    )?;

    // ---- (4a) The missing-data gate (esm-spec §10.10). ---------------------
    // A shaped parameter with neither a default nor a supplied value is a
    // construction error, raised here naming the parameter rather than as an
    // unbound read at the first evaluation. An array model is checked on its
    // compiled tables (after step 5 has bound the providers); a static one on
    // its flattened parameters, before its observed graph is evaluated.
    if let Backend::Static(_) = &backend
        && let Some(f) = owned_file.as_ref()
    {
        supplied.refuse_missing_static(f, &opts.p, model_name.as_deref())?;
    }

    // ---- (4b) State-free static evaluation. -------------------------------
    // A document that declares no differential equations has nothing to
    // integrate — `solve` refuses it with `NotDynamic` — but its whole content
    // is its observed graph, which is exactly what `observed_field` reads. This
    // evaluates that graph so the read works with NO options set, which is what
    // a stable-API function has to do (API_SPEC §5.8): before this, reaching a
    // state-free document's own results through the Rust binding meant knowing
    // to pass `build_pipeline: true` AND to hand in raw JSON rather than a
    // typed document, neither of which the surface contract mentions.
    //
    // Skipped when the build pipeline already produced fields, so the ISRM /
    // pushdown path keeps its own (model-local) keys untouched.
    //
    // Evaluated by the compiler the Problem names, and gated by it: the tape
    // under `native` / `xla`, the per-cell oracle under `interpreter`.
    if build.fields.is_empty()
        && let Backend::Static(_) = &backend
    {
        let t0 = opts.sample_time.unwrap_or(tspan.0);
        let (fields, rows, eval_error) = static_observed_fields(
            owned_file.as_ref(),
            flat_only,
            &opts.p,
            t0,
            model_name.as_deref(),
            compiler,
        )?;
        build.fields = fields;
        build.static_eval_error = eval_error;
        route_rows.extend(rows);
    }

    // ---- (5) Bind and CONST-materialize the run-time providers. -----------
    #[cfg(not(target_arch = "wasm32"))]
    let (refresh, discrete_forcing, refresh_boundaries) =
        bind_providers(&backend, &mut opts, tspan)?;
    if let Backend::Array(c) = &backend {
        supplied.refuse_missing_array(c, owned_file.as_ref(), &opts.p)?;
    }

    // ---- (5c) Refuse a data-fed parameter nothing bound (§9.6.6). ---------
    // After the providers are bound and BEFORE the compiler's gate below builds
    // the tape, so the refusal is the same under `native` and `interpreter`:
    // this asks whether the DOCUMENT's inputs are bound, not what a compiler
    // can lower. Without it `native` reported the same document as
    // `compiler_refused_rule` ("wholesale: unresolved symbol"), naming the
    // tape's limits in place of the defect, while `interpreter` built and then
    // failed inside the run with an uncoded treewalk fault.
    #[cfg(not(target_arch = "wasm32"))]
    if let Backend::Array(compiled) = &backend {
        crate::data_fed::refuse_unbound(
            compiled,
            &opts.const_arrays,
            &discrete_forcing,
            opts.p.keys(),
        )?;
    }

    // esm-spec §2.2: the document's own solver hints. Read from whichever
    // carrier survived to here — the raw JSON when there is one, else the typed
    // file — so the §2.2.2 chain holds however the problem was built.
    let doc_solver = match (doc_solver_of_text, owned_json.as_ref(), owned_file.as_ref()) {
        (Some(solver), _, _) => solver,
        (None, Some(raw), _) => document_solver(raw),
        (None, None, Some(file)) => crate::solver::normalize_empty(file.solver.clone()),
        (None, None, None) => None,
    };

    // ---- (5b) The compiler's own gate, and its report. --------------------
    // §2.5.10 puts the refusal at CONSTRUCTION, and puts every evaluation the
    // compiler performs for the Problem under it — not the right-hand side
    // alone. Both need the tape, so the tape is built HERE rather than inside
    // `solve`, and with the DISCRETE forcing set step (5) just resolved, which
    // is what decides a rule's cadence tier.
    let mut compiler_report = build_compiler_report(
        &backend,
        model_name.as_deref(),
        compiler,
        #[cfg(not(target_arch = "wasm32"))]
        &discrete_forcing,
        #[cfg(target_arch = "wasm32")]
        &std::collections::HashSet::new(),
        &opts.p,
    )?;
    compiler_report.rules.splice(0..0, route_rows);

    // ---- (5b') Every state slot has a starting value (esm-spec §11.4). ------
    // After (5b), which resolved the field `ic`s this reads.
    #[cfg(feature = "solve")]
    if let Backend::Array(c) = &backend {
        refuse_missing_initial_values(c, &opts.u0)?;
    }

    // ---- (5c) `xla`: emit and compile, here, once. ------------------------
    // AFTER the gate above, so a rule that never reached the tape is refused
    // by name with its cadence tier and component, rather than by the
    // emitter's own unqualified refusal of the declined tape. At construction rather than at the first
    // step, because §2.5.10 puts the refusal at construction and because XLA
    // compilation is the expensive step: a Problem pays it once, and every
    // segment and every Jacobian call reuses the one executable.
    install_xla_program(&backend, model_name.as_deref(), compiler)?;

    let program_source = match (&backend, owned_file, flat_only) {
        (Backend::Array(_), Some(file), _) => Some(Rc::new(ProgramSource::File(Box::new(file)))),
        (Backend::Array(_), None, Some(flat)) => {
            Some(Rc::new(ProgramSource::Flattened(Box::new(flat.clone()))))
        }
        _ => None,
    };
    let prob = EsmProblem {
        doc: Rc::new(match (owned_json, source_text) {
            (None, Some(text)) => RawDocument::from_text(text),
            (raw, _) => RawDocument::parsed(raw.unwrap_or(JsonValue::Null)),
        }),
        solver: doc_solver,
        model_name,
        precision: precision::Env::capture(),
        tspan,
        p: std::mem::take(&mut opts.p),
        u0: std::mem::take(&mut opts.u0),
        backend: Rc::new(backend),
        build: Rc::new(build),
        inspection: std::cell::RefCell::new(BuildInspection::default()),
        inspect: opts.inspect,
        compiler,
        compiler_report,
        callbacks: std::mem::take(&mut opts.callbacks),
        #[cfg(not(target_arch = "wasm32"))]
        refresh,
        #[cfg(not(target_arch = "wasm32"))]
        discrete_forcing,
        #[cfg(not(target_arch = "wasm32"))]
        refresh_boundaries,
        program_source,
        runtime_array_parameters: opts.runtime_array_parameters,
    };
    Ok(prob)
}

/// Read `domain.element_type` off whichever form of the input the caller gave
/// (esm-spec §11.3).
///
/// The raw-JSON branch reads the field textually because it runs *before* the
/// typed parse — the build pipeline needs the precision armed already.
///
/// # Errors
///
/// [`SimulateError::Compile`] wrapping
/// [`crate::compile_error::CompileError::UnsupportedElementType`] for a
/// spelling that is neither `"Float64"` nor `"Float32"`. An unrecognised
/// element type is not silently binary64.
fn element_type_of(
    raw: Option<&JsonValue>,
    file: Option<&EsmFile>,
    flat: Option<&FlattenedSystem>,
) -> Result<Precision, SimulateError> {
    let named = raw
        .and_then(|v| v.get("domain"))
        .and_then(|d| d.get("element_type"))
        .and_then(JsonValue::as_str)
        .map(str::to_string)
        .or_else(|| {
            file.and_then(|f| f.domain.as_ref())
                .and_then(|d| d.element_type.clone())
        })
        .or_else(|| {
            flat.and_then(|f| f.domain.as_ref())
                .and_then(|d| d.element_type.clone())
        });
    Precision::from_element_type(named.as_deref()).map_err(SimulateError::Compile)
}

/// Collect every `ModelVariable.element_type` the document declares
/// (esm-spec §11.3.1), keyed by both the bare and the `Model.name` spelling.
///
/// Read off whichever form of the input the caller gave, and — like
/// [`element_type_of`] — off the RAW JSON textually where there is one,
/// because it has to be armed before the build pipeline runs.
///
/// Returns an EMPTY table for the overwhelmingly common document that declares
/// none, which is what keeps the whole per-variable path inert.
///
/// # Errors
///
/// [`SimulateError::Compile`] wrapping
/// [`crate::compile_error::CompileError::UnsupportedElementType`] for a
/// spelling that is neither `"Float64"` nor `"Float32"`.
fn variable_element_types_of(
    raw: Option<&JsonValue>,
    file: Option<&EsmFile>,
    flat: Option<&FlattenedSystem>,
    document: Precision,
    selected: Option<&str>,
) -> Result<precision::VarPrecisions, SimulateError> {
    if let Some(models) = raw
        .and_then(|v| v.get("models"))
        .and_then(JsonValue::as_object)
    {
        let decls: Vec<precision::ModelDeclarations<'_>> = models
            .iter()
            .map(|(mname, model)| precision::ModelDeclarations {
                model: mname,
                variables: model
                    .get("variables")
                    .and_then(JsonValue::as_object)
                    .map(|vars| {
                        vars.iter()
                            .map(|(v, d)| {
                                (
                                    v.as_str(),
                                    d.get("element_type").and_then(JsonValue::as_str),
                                )
                            })
                            .collect()
                    })
                    .unwrap_or_default(),
                subsystem_variables: precision::subsystem_declarations(model.get("subsystems")),
            })
            .collect();
        return precision::VarPrecisions::from_declarations(&decls, document, selected)
            .map_err(SimulateError::Compile);
    }
    if let Some(models) = file.and_then(|f| f.models.as_ref()) {
        let decls = crate::precision_infer::model_declarations(models);
        return precision::VarPrecisions::from_declarations(&decls, document, selected)
            .map_err(SimulateError::Compile);
    }
    // A flattened system carries namespaced names and no per-variable element
    // types of its own; nothing to collect. `document` stays the answer for
    // every name.
    let _ = (flat, document);
    Ok(precision::VarPrecisions::default())
}

/// Reject integrating a Float32 document.
///
/// The ODE/DAE solver (diffsol) is instantiated over `f64` — its step-size
/// control, error norms and Newton iteration are binary64 and have no binary32
/// form here. Integrating a Float32 document would therefore produce an answer
/// whose right-hand side was binary32 and whose time-stepping was binary64,
/// with nothing in the result to say so. Per esm-spec §11.3 that is an error
/// naming the construct, not a partial honouring of the declaration.
///
/// Algebraic, observed and relational evaluation — everything `observed_field`
/// and the inline `tests` blocks read — is unaffected and runs in binary32.
fn reject_f32_integration(prob: &EsmProblem) -> Result<(), SimulateError> {
    // Only a system that actually INTEGRATES is affected. A model whose every
    // equation is algebraic — the relational / lookup-table shape this mode
    // exists for — is evaluated entirely by the (binary32) expression
    // evaluator, and `Rhs::Always` giving it a non-static backend does not
    // make it dynamic.
    let integrates = match &*prob.backend {
        Backend::Static(_) => false,
        Backend::Array(c) => c.has_differential_equations(),
    };
    if prob.precision.is_f32() && integrates {
        return Err(SimulateError::Compile(
            crate::compile_error::CompileError::Float32Unsupported {
                construct: "time integration of a dynamic model".to_string(),
                reason: "the ODE/DAE solver is instantiated over binary64 (step-size control,                          error norms and the Newton solve), so the trajectory would be                          binary64 even though the right-hand side is binary32"
                    .to_string(),
            },
        ));
    }
    Ok(())
}

/// Arrays and objects nested this deep or less parse from text well inside
/// `serde_json`'s recursion limit (128).
const NESTING_WITHOUT_ROUND_TRIP: usize = 64;

/// Whether no array or object in `v` sits more than `limit` containers deep
/// (the root container is depth 1).
fn nesting_within(v: &JsonValue, limit: usize) -> bool {
    match v {
        JsonValue::Array(items) => limit > 0 && items.iter().all(|x| nesting_within(x, limit - 1)),
        JsonValue::Object(map) => limit > 0 && map.values().all(|x| nesting_within(x, limit - 1)),
        _ => true,
    }
}

#[cfg(not(target_arch = "wasm32"))]
fn wants_build_pipeline(opts: &ProblemOptions) -> bool {
    opts.build_pipeline
        || opts.pushdown_rewrite
        || !opts.build_providers.is_empty()
        || !opts.const_arrays.is_empty()
}

/// The compiled model of a problem with **nothing to integrate** whose
/// observeds are a scalar function of the parameters and `t` alone — the
/// graph [`evaluate_static_observeds_over`] evaluates.
///
/// `esm simulate` evaluates such a document ONCE, at `tspan.0`, because a
/// single evaluation is all that command promises. The §6.6 inline-test runner
/// cannot stop there: esm-spec §6.6.3 defines an assertion's `time` as
/// "Simulation time at which to evaluate the assertion", and a document with
/// no differential equations is still a FUNCTION OF `t` — a solar-geometry
/// component's declination, hour angle and zenith cosine are exactly that. So
/// the runner evaluates this once per asserted time (issue #406).
///
/// `None` when there is something to integrate — the trajectory, not this, is
/// the answer — and `None` for a SHAPED state-free document (an array,
/// tensor or geometry observed, [`ArrayCompiled::observeds_are_scalar`]), which
/// is answered from the fields its BUILD materialized instead: see the
/// `has_nothing_to_integrate` arm of the inline-test runner, which is why that
/// runner asks the build pipeline for them. `None` again when a DISCRETE
/// provider feeds the document, whose forcing moves with `t` only through the
/// segmented refresh a solve drives.
#[cfg(all(not(target_arch = "wasm32"), feature = "solve"))]
pub(crate) fn static_observed_graph(prob: &EsmProblem) -> Option<Rc<ArrayCompiled>> {
    if !prob.discrete_forcing.is_empty() {
        return None;
    }
    match &*prob.backend {
        Backend::Array(c)
            if !c.has_differential_equations()
                && c.state_variable_names().is_empty()
                && c.observeds_are_scalar() =>
        {
            Some(Rc::clone(c))
        }
        _ => None,
    }
}

/// Every 0-D observed of [`static_observed_graph`]'s result at each of
/// `times`, against the problem's resolved parameters — `(name, values)` in
/// the document's own observed order, one value per time, named as a
/// [`solve`] of the same problem names its rows.
///
/// Evaluated by the array runtime on the compiler the problem was built with,
/// so under `native` the values come off the tape and under `interpreter` off
/// the per-cell oracle, exactly as a solve's reported observeds would.
#[cfg(all(not(target_arch = "wasm32"), feature = "solve"))]
pub(crate) fn evaluate_static_observeds_over(
    prob: &EsmProblem,
    graph: &ArrayCompiled,
    times: &[f64],
) -> Result<Vec<(String, Vec<f64>)>, SimulateError> {
    let _precision_guard = prob.precision.enter();
    let rows = graph.evaluate_stateless_observeds(&prob.p, times)?;
    let names: Vec<String> = rows.iter().map(|(name, _)| name.clone()).collect();
    Ok(qualify_array_names(graph, &names)
        .into_iter()
        .zip(rows.into_iter().map(|(_, values)| values))
        .collect())
}

/// Every observed field of a problem with **nothing to integrate** whose
/// compiled model is state-free, evaluated at `tspan.0` on the compiler the
/// problem was built with — `(name, field)` in dependency order, named as the
/// model's rules are.
///
/// The inline-test runner's answer for a SHAPED state-free document, which
/// [`static_observed_graph`] does not serve: the same array runtime the
/// problem already built, so under `native` the fields come off the tape
/// rather than from a second, pipeline build through the reference evaluator.
/// `None` when there is state or something to integrate, when a DISCRETE
/// provider feeds the document (its forcing moves only through a solve's
/// refresh), or when the evaluation fails or produces nothing.
#[cfg(all(not(target_arch = "wasm32"), feature = "solve"))]
pub(crate) fn stateless_fields(prob: &EsmProblem) -> Option<Vec<(String, ArrayD<f64>)>> {
    if !prob.discrete_forcing.is_empty() {
        return None;
    }
    let Backend::Array(c) = &*prob.backend else {
        return None;
    };
    if c.has_differential_equations() || !c.state_variable_names().is_empty() {
        return None;
    }
    let _precision_guard = prob.precision.enter();
    c.evaluate_stateless_fields(&prob.p, prob.tspan.0)
        .ok()
        .filter(|f| !f.is_empty())
}

/// Whether this problem has **nothing to integrate**: no backend at all, or a
/// backend that declares no differential equation (esm-spec §6.3.1 derives
/// `system_kind: "nonlinear"` for such a document).
///
/// This is the derivation `Rhs::Auto` consults to pick [`Backend::Static`],
/// asked of a problem the inline-test runner deliberately built with
/// `Rhs::Always` — which still compiles, and therefore still refuses, a
/// construct no evaluator supports, but which also hands a right-hand side over
/// an empty state vector to the solver. The runner asks this BEFORE solving, so
/// a document with nothing to integrate is answered from what was evaluated or
/// built rather than reported as a nonlinear-solver failure (issue #406).
///
/// [`static_observed_graph`] serves the scalar-observed documents; this
/// predicate is what lets the runner tell "nothing to integrate, and no scalar
/// graph to evaluate here" (a shaped document) apart from "something to
/// integrate".
#[cfg(all(not(target_arch = "wasm32"), feature = "solve"))]
pub(crate) fn has_nothing_to_integrate(prob: &EsmProblem) -> bool {
    match &*prob.backend {
        Backend::Static(_) => true,
        Backend::Array(c) => !c.has_differential_equations(),
    }
}

/// Evaluate a state-free system's observed graph at `t0`, on the compiler the
/// Problem names, and report every rule it evaluated.
///
/// The build-time half of [`observed_field`] for a document with no
/// differential equations: one stateless evaluation of the array runtime
/// ([`ArrayCompiled::evaluate_stateless_observeds`]) — no solver, no build
/// pipeline, and no dependence on how the caller spelled its input.
///
/// Keys are FLATTENED names (`Sites.North.u`), matching Julia's
/// `BuildInspection.observed_exprs` and Python's `static_observed_values`.
///
/// **The rank comes from the DECLARATION**, through the same
/// [`crate::classification::DeclaredRank`] the build pipeline uses (API_SPEC
/// `observed_trajectories`: "A field is the shape the document declares").
/// This path and the pipeline serve the SAME documents — this one runs exactly
/// when the pipeline produced no fields — so deciding the rank differently
/// here would give one document two spellings, `Rel.total` or `Rel.total[1]`,
/// picked by a backend choice the author cannot see.
///
/// **Under the compiler's refusal.** esm-libraries-spec §2.5.10 puts "the
/// materialization of constants and static observeds at construction" under
/// the refusal, and this IS that materialization for a document with nothing
/// to integrate: it has no right-hand side, but its observed graph is its
/// whole answer and something evaluates it. So the runtime is built in the
/// named compiler's mode and the SAME per-rule gate a dynamic document meets
/// ([`build_compiler_report`]) runs over it: under `native` / `xla` every rule
/// must lower to the tape, which then evaluates them, and a rule that does not
/// is a `compiler_refused_rule` naming it; under `interpreter` every rule is
/// the per-cell oracle's. `xla` emits no program here — there is no
/// right-hand side to emit — so its rules are reported as `"taped"`, which is
/// what evaluates them.
///
/// **Every observed, scalar or shaped.** A shaped state-free document — a
/// relational row set, a grid-geometry table — is evaluated here as a scalar
/// one is: its observeds are rules of the same array runtime, so the compiler
/// that answers a scalar observed answers a shaped one, rather than leaving it
/// to the build pipeline's reference evaluator. A 0-D or one-cell observed takes its rank from the declaration; a
/// shaped one is the field the runtime materialized, reported only when the
/// document declares it an observed (a shaped parameter the runtime lowered to
/// a `const` rule is not one).
///
/// **Tolerant of what no compiler can build.** A system the runtime cannot
/// build or evaluate — a `v1`-unsupported feature, a parameter with no default
/// — yields NO fields rather than failing the build, under every compiler
/// alike. Construction was not asked to compile a right-hand side
/// (`Rhs::Auto` chose `Backend::Static`), so a failure to evaluate must surface
/// as `observed_field` reporting the name it cannot answer, not as a document
/// that will not build. Python makes the same call for the same reason
/// (`problem.py`, the scalar no-state branch). Only the compiler's own refusal is
/// raised: that one says the document evaluates under `interpreter` and not
/// under the compiler the caller named. An EVALUATION failure of a runtime
/// that did build (a fail-closed fault) comes back as the third element, for
/// [`observed_field`] to report.
#[allow(clippy::type_complexity)]
fn static_observed_fields(
    file: Option<&EsmFile>,
    flat_only: Option<&FlattenedSystem>,
    p: &HashMap<String, f64>,
    t0: f64,
    model_name: Option<&str>,
    compiler: Compiler,
) -> Result<
    (
        HashMap<String, ArrayD<f64>>,
        Vec<CompilerRuleReport>,
        Option<String>,
    ),
    SimulateError,
> {
    let nothing = || Ok((HashMap::new(), Vec::new(), None));
    let owned_flat;
    let flat = match (flat_only, file) {
        (Some(f), _) => f,
        // A document that does not flatten — an unresolved coupling endpoint,
        // a coupling-import error — is a build failure, raised as one
        // (esm-libraries-spec §2.5.2) rather than answered with no fields.
        (None, Some(file)) => {
            owned_flat = crate::flatten::flatten(file).map_err(|e| {
                SimulateError::Compile(crate::compile_error::CompileError::Flatten(e))
            })?;
            &owned_flat
        }
        (None, None) => return nothing(),
    };
    // The evaluable-core gate runs first (esm-spec §9.6.3 constraint 6): it
    // walks every equation, so an unknown no equation defines — which makes
    // the system not state-free below — does not let a rewrite-target operator
    // through unreported.
    if let Some(op) = crate::flatten::first_unlowered_operator(flat) {
        return Err(SimulateError::Compile(
            crate::compile_error::CompileError::UnloweredOperatorError { op },
        ));
    }
    // A system WITH state is not state-free evaluable: an observed may read
    // state, and there is none to read. Reachable when `model_name` selects an
    // ODE-free model out of a document that has ODEs elsewhere, since
    // `flatten` is document-wide.
    // An unknown defined by an element-wise equation (`index(a, i) = …`) is
    // listed as an observed as well as a state; only a state no equation
    // defines has anything to integrate.
    if flat
        .state_variables
        .keys()
        .any(|k| !flat.observed_variables.contains_key(k))
    {
        return nothing();
    }
    // An ill-formed document is not a compiler's to refuse: a name declared
    // nowhere fails the same free-variable gate the array build and the build
    // pipeline run (CONFORMANCE_SPEC §5.23), and is raised as that build
    // failure rather than surfacing as a strict compiler "declining" a rule no
    // compiler could evaluate — or as a Problem with no fields at all.
    #[cfg(not(target_arch = "wasm32"))]
    if let Some(file) = file
        && crate::simulate::whole_document_is_one_model(file)
        && let Some((_, model)) = file.models.as_ref().and_then(|m| m.iter().next())
    {
        let index_sets: HashMap<String, crate::types::IndexSet> = file
            .index_sets
            .clone()
            .unwrap_or_default()
            .into_iter()
            .collect();
        crate::simulate_array::check_free_variables(model, &index_sets, &[])
            .map_err(SimulateError::Compile)?;
    }
    // The same gates the array route runs — the evaluable-core operator gate
    // (esm-spec §9.6.3 constraint 6), events, operator arity — fire here too:
    // a state-free document is evaluated, so it is built.
    let mut compiled = ArrayCompiled::from_flattened(flat).map_err(SimulateError::Compile)?;
    compiled.runtime_mode = match compiler {
        Compiler::Interpreter => crate::simulate_array::RuntimeMode::Interpreter,
        _ => crate::simulate_array::RuntimeMode::Native,
    };
    let compiled = Rc::new(compiled);
    let report = rule_report(
        &compiled,
        model_name,
        compiler,
        &std::collections::HashSet::new(),
        false,
    )?;
    let values = match compiled.evaluate_stateless_fields(p, t0) {
        Ok(values) => values,
        // A value the spec makes invalid (a degenerate polygon operand,
        // esm-spec §8.6.1) is refused, not handed back as a Problem with no
        // fields; a value that is only not yet available (a parameter with no
        // default, a field a provider fills later) still is.
        Err(e)
            if e.to_string()
                .contains(crate::simulate_array::GEOMETRY_CLIP_CODE) =>
        {
            return Err(e);
        }
        Err(e) => return Ok((HashMap::new(), report.rules, Some(e.to_string()))),
    };
    let fields = values
        .into_iter()
        .filter_map(|(name, arr)| {
            let declared = flat.observed_variables.get(&name);
            if arr.len() == 1 {
                let shape = DeclaredRank::of_declaration(declared).scalar_shape();
                let v = arr.iter().copied().next().expect("one element");
                let arr = ArrayD::from_shape_vec(ndarray::IxDyn(shape), vec![v])
                    .expect("`scalar_shape` is `[]` or `[1]`, and both hold exactly one element");
                return Some((name, arr));
            }
            declared.map(|_| (name, arr))
        })
        .collect();
    Ok((fields, report.rules, None))
}

/// Build the per-rule record, and — under a STRICT compiler ([`Compiler::Native`]
/// or [`Compiler::Xla`]) — REFUSE the document when any rule failed to lower
/// (esm-libraries-spec §2.5.10).
///
/// This is the whole of that strictness, and it lives at construction for the
/// reason §2.5.2 gives: a document that never became a Problem has no run to
/// describe, so a build failure must be raised where the build happens.
///
/// **It gates more than the right-hand side.** §2.5.10 puts "the
/// materialization of constants and static observeds at construction, the
/// per-segment seed, the right-hand side, and the observeds reported at output
/// times" all under the refusal. In this binding every one of those is an
/// evaluation of the SAME rule set, so one gate over the rule set covers them
/// all — and the driver then SERVES those passes from the tape under both
/// strict compilers (`simulate_array::driver`) instead of from the whole-array
/// overlay, which is what makes the gate honest rather than merely necessary.
///
/// `xla` is gated HERE as well as by its emitter, and the order matters: this
/// gate names the rule, its component and its cadence tier, while an emitter
/// refusal names only the rule as the compiled model spells it, or the
/// instruction it choked on. A rule that never reached the tape is therefore
/// reported here, in full.
fn build_compiler_report(
    backend: &Backend,
    model_name: Option<&str>,
    compiler: Compiler,
    discrete_forcing: &std::collections::HashSet<String>,
    #[cfg_attr(not(feature = "solve"), allow(unused_variables))] p: &HashMap<String, f64>,
) -> Result<CompilerReport, SimulateError> {
    let compiled = match backend {
        // No compiled model. A state-free document's observed graph was
        // evaluated — and gated — by `static_observed_fields`, whose rows the
        // caller adds.
        Backend::Static(_) => return Ok(CompilerReport::empty(compiler)),
        Backend::Array(c) => c,
    };
    #[cfg_attr(not(feature = "solve"), allow(unused_mut))]
    let mut report = rule_report(compiled, model_name, compiler, discrete_forcing, true)?;
    // The field initial conditions are evaluated by the reference evaluator
    // under every compiler — no tier has a form for initial-state assembly —
    // so they are exercised HERE, where a strict compiler refuses one it would
    // have to walk per cell, and reported ahead of the model's rules.
    #[cfg(feature = "solve")]
    {
        let model = model_name.unwrap_or("");
        let mut ic_rows = Vec::new();
        for r in compiled.field_ic_records(p, compiler.is_strict()) {
            let reason = r.per_cell.then(|| {
                "evaluated by the reference evaluator at setup, and a `faq` in it is not one \
                 the whole-array overlay takes, so it was walked once per cell"
                    .to_string()
            });
            if compiler.is_strict()
                && let Some(reason) = reason
            {
                // The walk stopped at its first cell (`field_ic_records`).
                return Err(SimulateError::Compile(
                    crate::compile_error::CompileError::CompilerRefusedRule {
                        compiler: compiler.as_str(),
                        kind: r.kind,
                        rule: qualify(model, &r.name),
                        tier: "const",
                        reason: format!("{reason}; {}", crate::simulate_array::ONE_CELL_NOTE),
                    },
                ));
            }
            ic_rows.push(outside_rule_row(
                compiler,
                qualify(model, &r.name),
                r.kind,
                reason,
                compiler == Compiler::Interpreter,
            ));
        }
        report.rules.splice(0..0, ic_rows);
    }
    Ok(report)
}

/// The report row of one evaluation construction performs OUTSIDE the
/// compiled model's rule set, at setup (`"const"`). `per_cell` is why it was
/// walked per cell, if it was; `overlay_off` says the evaluation ran with the
/// whole-array overlay off, as a field `ic` does under
/// [`Compiler::Interpreter`].
///
/// The tier says what happened: the per-cell oracle, or the whole-array
/// overlay. Under the interpreter a per-cell walk is that compiler's design
/// rather than a decline, so its rows carry no reason, exactly as its rule
/// rows do.
fn outside_rule_row(
    compiler: Compiler,
    rule: String,
    kind: &'static str,
    per_cell: Option<String>,
    overlay_off: bool,
) -> CompilerRuleReport {
    let tier = if overlay_off || per_cell.is_some() {
        "oracle"
    } else {
        "vectorized"
    };
    let reason = match compiler {
        Compiler::Interpreter => None,
        _ => per_cell,
    };
    CompilerRuleReport {
        rule,
        kind,
        cadence: "const",
        tier,
        reason,
    }
}

/// The per-rule record of one compiled model, refusing under a strict compiler
/// the first rule that did not lower ([`build_compiler_report`]).
///
/// `xla_emitted` says whether the model's right-hand side is the emitted XLA
/// program. A state-free model has none, so under [`Compiler::Xla`] its rules
/// are reported as the tape's, which is what evaluates them.
fn rule_report(
    compiled: &ArrayCompiled,
    model_name: Option<&str>,
    compiler: Compiler,
    discrete_forcing: &std::collections::HashSet<String>,
    xla_emitted: bool,
) -> Result<CompilerReport, SimulateError> {
    let (records, tape_report) = compiled.tape_rule_records(discrete_forcing);
    let model = model_name.unwrap_or("");
    let mut rules = Vec::with_capacity(records.len());
    for r in records {
        let qualified = qualify(model, &r.name);
        if compiler.is_strict()
            && let Some(reason) = &r.fallback_reason
        {
            return Err(SimulateError::Compile(
                crate::compile_error::CompileError::CompilerRefusedRule {
                    compiler: compiler.as_str(),
                    kind: r.kind,
                    rule: qualified,
                    tier: r.cadence,
                    reason: reason.clone(),
                },
            ));
        }
        let tier = match (compiler, &r.fallback_reason) {
            // Under `interpreter` the tape is never installed, so every rule
            // runs on the per-cell oracle by design — not as a decline. The
            // tape build above ran only to classify each rule's cadence and is
            // discarded; it evaluates nothing.
            (Compiler::Interpreter, _) => "oracle",
            // Under `xla` the gate above has already refused anything that did
            // not lower, so every rule reaching here is in the emitted
            // program. It is reported as `"xla"` rather than as `"taped"`
            // because the tape is where it was lowered FROM, not what
            // evaluates it.
            (Compiler::Xla, _) if xla_emitted => "xla",
            (_, None) => "taped",
            (_, Some(_)) => "fallback",
        };
        rules.push(CompilerRuleReport {
            rule: qualified,
            kind: r.kind,
            cadence: r.cadence,
            tier,
            // Under `interpreter` there is no decline to report: the tape was
            // never going to run, so what it would have said about a rule is
            // not this compiler's story.
            reason: match compiler {
                Compiler::Interpreter => None,
                _ => r.fallback_reason,
            },
        });
    }
    let (fused_groups, fused_instructions) = if compiler == Compiler::Interpreter {
        (0, 0)
    } else {
        (tape_report.fuse.n_groups, tape_report.fuse.n_member_instrs)
    };
    Ok(CompilerReport {
        compiler,
        rules,
        fused_groups,
        fused_instructions,
    })
}

/// Why a compiler in the vocabulary cannot be provided here, and what would
/// provide it (esm-libraries-spec §2.5.10: the message names what to load or
/// build). `None` for the two this binding implements.
///
/// A refusal, never a substitution — answering `xla` by building `native`
/// would make the compiler's name describe nothing, which is the failure
/// §2.5.10 exists to prevent.
fn unavailable_reason(compiler: Compiler) -> Option<String> {
    match compiler {
        Compiler::Native | Compiler::Interpreter => None,
        Compiler::Xla => xla_unavailable_reason(),
        Compiler::Mtk => Some(
            "a ModelingToolkit system needs Julia and ModelingToolkit.jl, and this binding has \
             no ModelingToolkit runtime to load; `mtk` is provided by the Julia binding"
                .to_string(),
        ),
        Compiler::Sympy => Some(
            "a lambdified SymPy right-hand side needs Python and SymPy, and this binding has \
             neither; `sympy` is provided by the Python binding"
                .to_string(),
        ),
    }
}

/// Whether `xla` can be provided by THIS build and THIS process, and what is
/// missing when it cannot.
///
/// `xla` is the specialty compiler that needs a heavy external dependency, and
/// the dependency can be absent in two different places: the crate can be
/// built without the `xla` feature, or it can have the feature and still find
/// no usable XLA runtime to start (no unpacked `xla_extension`, no device for
/// the requested platform). Both are `compiler_unavailable` and neither says
/// anything about what the emitter can lower, so both are answered here —
/// BEFORE any document is looked at, so the answer is a property of the build
/// and not of the input.
///
/// The client is a process-wide singleton and caches its own failure, so
/// probing it per Problem costs one atomic read after the first.
#[cfg(feature = "xla")]
fn xla_unavailable_reason() -> Option<String> {
    crate::xla_runtime::client().err().map(|e| {
        format!(
            "this build has the `xla` feature, but no usable XLA runtime in this process: {e}. \
             Point XLA_EXTENSION_DIR at an unpacked xla_extension release \
             (scripts/fetch-xla-extension.sh) and rebuild against it"
        )
    })
}

#[cfg(not(feature = "xla"))]
fn xla_unavailable_reason() -> Option<String> {
    Some(
        "this binding's StableHLO emitter over the tape is behind the `xla` Cargo feature, \
         which this build does not have: rebuild with `--features xla` and XLA_EXTENSION_DIR \
         pointing at an unpacked xla_extension release (scripts/fetch-xla-extension.sh)"
            .to_string(),
    )
}

/// Emit and compile the right-hand side for [`Compiler::Xla`], once, on the
/// backend the build just produced.
///
/// A no-op for every other compiler, and for a backend with no right-hand side
/// to emit: a state-free document is not integrated, so no compiler is chosen
/// for it and there is nothing to lower.
///
/// The two failures are kept apart, because they say different things:
///
///   * the emitter REFUSED this model — a named
///     [`CompileError::CompilerRefusedRule`], the same failure a rule the tape
///     could not lower raises, so a caller (and the compiler-agreement tier)
///     reads one shape whichever half declined;
///   * the XLA runtime broke — [`SimulateError::CompilerUnavailable`], which
///     is a fact about this process and says nothing about the document.
fn install_xla_program(
    backend: &Backend,
    model_name: Option<&str>,
    compiler: Compiler,
) -> Result<(), SimulateError> {
    // `model_name` is read only on the `xla` build; without the feature the
    // whole body is the early return, and `compiler` is already known to be
    // something this build provides.
    let _ = (backend, model_name);
    if compiler != Compiler::Xla {
        return Ok(());
    }
    #[cfg(feature = "xla")]
    {
        use crate::xla_runtime::CompileRhsError;
        let Backend::Array(compiled) = backend else {
            return Ok(());
        };
        if !compiled.has_differential_equations() {
            return Ok(());
        }
        match compiled.install_xla_rhs() {
            Ok(()) => Ok(()),
            Err(CompileRhsError::Refused(e)) => {
                Err(SimulateError::Compile(
                    crate::compile_error::CompileError::CompilerRefusedRule {
                        compiler: Compiler::Xla.as_str(),
                        // The emitter refuses per INSTRUCTION, and an instruction
                        // belongs to a rule whose kind it does not carry, so the
                        // kind is the neutral one rather than a guess between
                        // "observed" and "state derivative".
                        kind: "rule",
                        rule: qualify(model_name.unwrap_or(""), &e.rule),
                        // The emitted program computes every cadence section on
                        // every call (`xla_emit`'s module docs), so a refusal
                        // anywhere in it costs on the continuous schedule.
                        tier: "continuous",
                        reason: e.reason,
                    },
                ))
            }
            Err(CompileRhsError::Runtime(details)) => Err(SimulateError::CompilerUnavailable {
                compiler: Compiler::Xla.as_str(),
                details,
            }),
        }
    }
    #[cfg(not(feature = "xla"))]
    {
        // Unreachable: `compile_backend` answered `compiler_unavailable` for
        // `xla` before any backend existed.
        Ok(())
    }
}

fn compile_backend(
    file: Option<&EsmFile>,
    flat: Option<&FlattenedSystem>,
    model_name: Option<&str>,
    mode: Rhs,
    compiler: Compiler,
    runtime_array_parameters: bool,
) -> Result<Backend, SimulateError> {
    if let Some(details) = unavailable_reason(compiler) {
        return Err(SimulateError::CompilerUnavailable {
            compiler: compiler.as_str(),
            details,
        });
    }
    // An equation that defines a parameter (esm-spec §6.3.1) is refused under
    // every compiler and every `Rhs` mode: the build used to drop it, answering
    // with the parameter's default.
    refuse_parameter_definitions(file, flat)?;
    if mode == Rhs::Never {
        return Ok(Backend::Static(
            "the caller asked for Rhs::Never".to_string(),
        ));
    }
    // The array runtime's own mode, carried on the compiled artifact so the
    // driver reads WHICH compiler it is serving rather than re-deciding per
    // call (`crate::simulate_array::RuntimeMode`).
    let runtime_mode = match compiler {
        Compiler::Interpreter => crate::simulate_array::RuntimeMode::Interpreter,
        Compiler::Xla => crate::simulate_array::RuntimeMode::Xla,
        _ => crate::simulate_array::RuntimeMode::Native,
    };
    if let Some(flat) = flat {
        // Same rule the typed-document branch applies below, on the input that
        // states it most directly: a flattened system with no state variables
        // has nothing to integrate, so `solve` reports `NotDynamic` rather than
        // failing inside the solver, and `observed_field` sees the static
        // fields.
        if mode == Rhs::Auto && flat.state_variables.is_empty() {
            return Ok(Backend::Static(
                "the flattened system declares no state variables".to_string(),
            ));
        }
        let mut compiled = crate::simulate_array::ArrayCompiled::from_flattened_with(
            flat,
            runtime_array_parameters,
        )?;
        compiled.runtime_mode = runtime_mode;
        return Ok(Backend::Array(Rc::new(compiled)));
    }
    let Some(file) = file else {
        if mode == Rhs::Always {
            return Err(SimulateError::Compile(
                crate::compile_error::CompileError::build_err("no typed document to compile"),
            ));
        }
        return Ok(Backend::Static(
            "the prepared document has no typed form to compile".to_string(),
        ));
    };

    // A document with no component at all — a template or coupling library, a
    // `data_sources` registry — has nothing to build. Julia's and Python's
    // `flatten` refuse it the same way, rather than hand back an empty Problem.
    let n_components = file.models.as_ref().map_or(0, |m| m.len())
        + file.reaction_systems.as_ref().map_or(0, |r| r.len());
    if n_components == 0 {
        return Err(SimulateError::Compile(
            crate::compile_error::CompileError::build_err(format!(
                "nothing to flatten: '{}' declares no `models` and no `reaction_systems`. A \
                 file carrying only `expression_templates` is a template LIBRARY (esm-spec \
                 §9.7); import it from a document that declares components rather than \
                 building it directly",
                file.metadata.name.as_deref().unwrap_or("")
            )),
        ));
    }

    refuse_structurally_unreachable(file)?;
    refuse_unregistered_callback_reads(file)?;
    refuse_reference_integrity_errors(file)?;

    if mode == Rhs::Auto && !has_differential_equations(file, model_name) {
        // Nothing to integrate is not nothing to run: an event, an implicit
        // equation or a noise term still changes what the document means, and
        // the static evaluation below honours none of them. The array route
        // refuses each by name, and so does this one.
        refuse_static_unsupported_construct(file, model_name)?;
        return Ok(Backend::Static(
            "the document declares no differential equations".to_string(),
        ));
    }

    // §2.5.10: "A binding MUST NOT switch strategy inside `native` on document
    // content: a scalar document and a gridded one are built by the same
    // machinery, so that what ran is a property of the name and not of the
    // input." So every document, scalar or gridded, builds the array runtime.
    let mut compiled = crate::simulate::build_array_compiled(file, runtime_array_parameters)?;
    compiled.runtime_mode = runtime_mode;
    Ok(Backend::Array(Rc::new(compiled)))
}

/// Refuse an equation whose left-hand side names a parameter
/// (`equation_defines_parameter`, esm-spec §6.3.1), in the typed document or,
/// for a flattened input, in the flattened equations. `validate` reports the
/// same finding; Julia and Python refuse it at their front doors too.
fn refuse_parameter_definitions(
    file: Option<&EsmFile>,
    flat: Option<&FlattenedSystem>,
) -> Result<(), SimulateError> {
    let refuse = |path: &str, message: &str| {
        SimulateError::Compile(crate::compile_error::CompileError::build_err(format!(
            "[{}] {path}: {message} (esm-spec §6.3.1)",
            crate::diagnostic::codes::EQUATION_DEFINES_PARAMETER
        )))
    };
    if let Some(file) = file
        && let Some(e) = crate::structural::parameter_definition_errors(file).first()
    {
        return Err(refuse(&e.path, &e.message));
    }
    if let Some(flat) = flat {
        for (k, eq) in flat.equations.iter().enumerate() {
            let mut lhs = &eq.lhs;
            let name = loop {
                match lhs {
                    crate::types::Expr::Variable(v) => break Some(v.as_str()),
                    crate::types::Expr::Operator(n) if n.op == "faq" || n.op == "aggregate" => {
                        match n.expr.as_deref() {
                            Some(x) => lhs = x,
                            None => break None,
                        }
                    }
                    crate::types::Expr::Operator(n)
                        if n.op == "index"
                            || (n.op == "D" && matches!(n.wrt.as_deref(), None | Some("t"))) =>
                    {
                        match n.args.first() {
                            Some(x) => lhs = x,
                            None => break None,
                        }
                    }
                    _ => break None,
                }
            };
            if let Some(name) = name
                && flat.parameters.contains_key(name)
            {
                return Err(refuse(
                    &format!("/equations/{k}/lhs"),
                    &format!(
                        "Equation {k} defines '{name}', which is a parameter; an equation \
                         defines unknowns only"
                    ),
                ));
            }
        }
    }
    Ok(())
}

/// Refuse two declarations the document can hold but no build may run, each a
/// structural error `validate` also reports, which built anyway would be
/// finite, plausible and wrong:
///
/// * a declaration spelled with the independent variable or `_var`
///   (`reserved_variable_name`, esm-spec §4.9.1.1): every reader of the name
///   gets the implicit symbol, never the declared quantity;
/// * an `ic` equation among a reaction system's `constraint_equations`
///   (`ic_in_reaction_system`, esm-spec §11.4.1), which no evaluator applies.
///
/// Julia and Python refuse both at their front doors too.
fn refuse_structurally_unreachable(file: &EsmFile) -> Result<(), SimulateError> {
    let refuse = |code: &str, path: &str, message: &str, section: &str| {
        SimulateError::Compile(crate::compile_error::CompileError::build_err(format!(
            "[{code}] {path}: {message} (esm-spec {section})"
        )))
    };
    if let Some(e) = crate::structural::reserved_declaration_errors(file).first() {
        return Err(refuse(&e.code.to_string(), &e.path, &e.message, "§4.9.1.1"));
    }
    if let Some(systems) = &file.reaction_systems {
        for (rs_name, rs) in systems {
            for (i, eq) in rs.constraint_equations.iter().flatten().enumerate() {
                if matches!(&eq.lhs, crate::types::Expr::Operator(n) if n.op == "ic") {
                    return Err(refuse(
                        crate::diagnostic::codes::IC_IN_REACTION_SYSTEM,
                        &format!("/reaction_systems/{rs_name}/constraint_equations/{i}"),
                        "ic equation not allowed in a reaction system; a reaction system \
                         hosts no ic equations (a species' initial value is its \
                         `species.default`, or a scoped-reference ic equation in a model)",
                        "§11.4.1",
                    ));
                }
            }
        }
    }
    Ok(())
}

/// The structural-validation codes a build refuses on (esm-libraries-spec
/// §2.5.2). They are the reference-integrity findings: a name, reference or data
/// source the document uses and does not declare. A build that went ahead would
/// read a value the document does not describe. Equation-count and unit findings
/// are not here: they stay `validate`'s to report.
const BUILD_REFUSED_VALIDATION_CODES: &[&str] = &[
    "undefined_variable",
    "undefined_parameter",
    "undefined_species",
    "undefined_system",
    "undefined_index_set",
    "unresolved_scoped_ref",
    "event_var_undeclared",
    "data_source_undefined",
    "missing_required_field",
];

/// Refuse the first [`BUILD_REFUSED_VALIDATION_CODES`] finding `validate`
/// reports for `file`, with the validator's code, pointer and message.
fn refuse_reference_integrity_errors(file: &EsmFile) -> Result<(), SimulateError> {
    let result =
        crate::validate::validate_scoped(file, crate::validate::ValidationScope::References);
    if let Some(e) = result.structural_errors.iter().find(|e| {
        // An inline test's references are the test runner's to report
        // (esm-spec §6.6); the build does not evaluate them.
        BUILD_REFUSED_VALIDATION_CODES.contains(&e.code.to_string().as_str())
            && !e.path.contains("/tests/")
    }) {
        let named = declared_name_at(&e.path)
            .map(|n| format!(" (variable '{n}')"))
            .unwrap_or_default();
        return Err(SimulateError::Compile(
            crate::compile_error::CompileError::build_err(format!(
                "[{}] {}: {}{named} (esm-libraries-spec §2.5.2)",
                e.code, e.path, e.message
            )),
        ));
    }
    Ok(())
}

/// The full name of the variable a validator pointer sits under —
/// `/models/A/subsystems/B/variables/v/…` is `A.B.v` — so a refusal names the
/// variable the way a caller addresses it (esm-spec §6.6.2).
fn declared_name_at(pointer: &str) -> Option<String> {
    let segs: Vec<&str> = pointer.trim_start_matches('/').split('/').collect();
    let mut parts = Vec::new();
    let mut i = 0;
    while i + 1 < segs.len() {
        match segs[i] {
            "models" | "reaction_systems" | "subsystems" => parts.push(segs[i + 1]),
            "variables" | "species" | "parameters" => {
                parts.push(segs[i + 1]);
                return Some(parts.join("."));
            }
            _ => return None,
        }
        i += 2;
    }
    None
}

/// Refuse, with `callback_unregistered` (esm-spec §9.6.6), an equation that
/// reads a variable a `callback` coupling injects
/// (`config.callback_variables[].name`). Only a registered callback supplies
/// that value, and nothing registers one at construction, so a build would read
/// a placeholder the document does not describe. The interpreter would read
/// an unset forcing slot at its first call instead.
fn refuse_unregistered_callback_reads(file: &EsmFile) -> Result<(), SimulateError> {
    let mut injected: BTreeMap<String, String> = BTreeMap::new();
    for entry in file.coupling.iter().flatten() {
        let crate::CouplingEntry::Callback {
            callback_id,
            config,
            ..
        } = entry
        else {
            continue;
        };
        let names = config
            .as_ref()
            .and_then(|c| c.get("callback_variables"))
            .and_then(|v| v.as_array())
            .into_iter()
            .flatten()
            .filter_map(|cv| cv.get("name").and_then(|n| n.as_str()));
        for name in names {
            injected.insert(name.to_string(), callback_id.clone());
        }
    }
    if injected.is_empty() {
        return Ok(());
    }
    let Some(models) = file.models.as_ref() else {
        return Ok(());
    };
    let mut names: Vec<&String> = models.keys().collect();
    names.sort();
    for name in names {
        callback_reads_in(&models[name], name, &injected)?;
    }
    Ok(())
}

/// [`refuse_unregistered_callback_reads`] over one model and its inline
/// subsystems. A name the model declares itself is not the injected one.
fn callback_reads_in(
    model: &crate::types::Model,
    path: &str,
    injected: &BTreeMap<String, String>,
) -> Result<(), SimulateError> {
    for eq in &model.equations {
        for side in [&eq.lhs, &eq.rhs] {
            let mut vars: Vec<String> = crate::free_variables(side).into_iter().collect();
            vars.sort();
            for v in vars {
                if model.variables.contains_key(&v) {
                    continue;
                }
                if let Some(id) = injected.get(&v) {
                    return Err(SimulateError::Compile(
                        crate::compile_error::CompileError::build_err(format!(
                            "[{}] '{path}' reads '{v}', which the `callback` coupling '{id}' \
                             supplies, but no callback is registered to supply it at \
                             construction; refusing the build rather than reading a value \
                             the document does not give (esm-spec §9.6.6)",
                            crate::diagnostic::codes::CALLBACK_UNREGISTERED
                        )),
                    ));
                }
            }
        }
    }
    if let Some(subs) = model.subsystems.as_ref() {
        let mut names: Vec<&String> = subs.keys().collect();
        names.sort();
        for name in names {
            if let Ok((sub, _)) = crate::simulate_array::parse_subsystem_model(name, &subs[name]) {
                callback_reads_in(&sub, &format!("{path}.{name}"), injected)?;
            }
        }
    }
    Ok(())
}

/// Refuse, with `unsupported_construct` (esm-spec §9.6.6), a document with
/// nothing to integrate that still carries an event, an implicit equation or a
/// Wiener-noise parameter — in the named model, or in any model when none is
/// named, subsystems included. The array route refuses the same three while
/// compiling; a state-free document never reaches it.
fn refuse_static_unsupported_construct(
    file: &EsmFile,
    model_name: Option<&str>,
) -> Result<(), SimulateError> {
    if model_name.is_none()
        && let Some((construct, event)) = crate::compile_error::first_event_in_file(file)
    {
        return Err(crate::compile_error::event_refusal(
            construct,
            crate::compile_error::ARRAY_EVALUATOR,
            event.as_deref(),
        )
        .into());
    }
    let Some(models) = file.models.as_ref() else {
        return Ok(());
    };
    for (_, model) in models
        .iter()
        .filter(|(name, _)| model_name.is_none_or(|want| want == name.as_str()))
    {
        if let Some(refusal) = static_unsupported_construct(model) {
            return Err(refusal.into());
        }
    }
    Ok(())
}

/// The first construct of `model` or its inline subsystems that
/// [`refuse_static_unsupported_construct`] refuses, as its refusal.
fn static_unsupported_construct(
    model: &crate::types::Model,
) -> Option<crate::compile_error::CompileError> {
    use crate::compile_error::ARRAY_EVALUATOR;
    if let Some((construct, event)) = crate::compile_error::first_event(model) {
        return Some(crate::compile_error::event_refusal(
            construct,
            ARRAY_EVALUATOR,
            event.as_deref(),
        ));
    }
    if let Some(eq) = crate::compile_error::first_implicit_equation(&model.equations) {
        return Some(crate::compile_error::implicit_equation_refusal(
            ARRAY_EVALUATOR,
            eq,
        ));
    }
    let class = crate::classification::Classification::of(model);
    if let Some(name) = class.brownian_parameters.first() {
        return Some(crate::compile_error::CompileError::UnsupportedConstruct {
            construct: crate::compile_error::WIENER_NOISE,
            evaluator: ARRAY_EVALUATOR,
            detail: format!("parameter '{name}'"),
        });
    }
    let subs = model.subsystems.as_ref()?;
    let mut names: Vec<&String> = subs.keys().collect();
    names.sort();
    names.into_iter().find_map(|name| {
        crate::simulate_array::parse_subsystem_model(name, &subs[name])
            .ok()
            .and_then(|(sub, _)| static_unsupported_construct(&sub))
    })
}

/// Whether the document (or the named model within it) declares at least one
/// differential equation, i.e. whether there is anything to integrate.
///
/// A document with none is a *dispatched static evaluation*: it still gets a
/// EsmProblem, and [`observed_field`] still reads its build-time products, but
/// [`solve`] on it raises [`SimulateError::NotDynamic`] rather than handing
/// back an empty trajectory.
pub(crate) fn has_differential_equations(file: &EsmFile, model_name: Option<&str>) -> bool {
    // A REACTION SYSTEM is differential, and reading only `models` said it was
    // not. Reactions lower to `D(species, t) = …` during flattening, so a
    // document whose whole content is a `reaction_systems` block has no models
    // AT ALL until that happens — and this function runs before it. Every pure
    // chemistry document in the wild is that shape (`pollu`, `superfast`,
    // `geoschem_fullchem`), and `Rhs::Auto` was calling each of them static
    // and then refusing to solve twenty-five differential equations with
    // `NotDynamic { "the document declares no differential equations" }`.
    //
    // A system with no reactions is genuinely not differential — it lowers to
    // nothing — so the check is for reactions rather than for the block.
    if let Some(systems) = file.reaction_systems.as_ref()
        && systems
            .iter()
            .filter(|(name, _)| model_name.is_none_or(|want| want == name.as_str()))
            .any(|(_, rs)| !rs.reactions.is_empty())
    {
        return true;
    }
    let Some(models) = file.models.as_ref() else {
        return false;
    };
    models
        .iter()
        .filter(|(name, _)| model_name.is_none_or(|want| want == name.as_str()))
        .any(|(_, m)| crate::simulate_array::model_tree_any(m, &model_has_derivative))
}

/// Whether any equation of `model` has a TIME derivative in its LHS — the
/// esm-spec §6.3.1 derivation of `system_kind`, which is what "declares no
/// differential equations" means (API_SPEC §5.8).
///
/// A bare top-level `D` is not the only spelling. The arrayed form
/// `faq{expr: D(index(u, i))}` is how every array-op PDE is written and what
/// the pointwise lift produces, and `index(D(u), i)` also names a derivative.
/// Each of them makes the document dynamic (issue #476): the array runtime
/// then either compiles the equation or refuses it by name, where
/// [`Backend::Static`] would have `solve` answer `NotDynamic` for a document
/// that plainly integrates.
fn model_has_derivative(model: &crate::types::Model) -> bool {
    model
        .equations
        .iter()
        .any(|eq| crate::classification::has_time_derivative(&eq.lhs))
}

/// Whether `var` is a parameter with a non-empty declared `shape`.
fn is_shaped_parameter(var: &crate::types::ModelVariable) -> bool {
    var.var_type == crate::types::VariableType::Parameter
        && var.shape.as_ref().is_some_and(|s| !s.is_empty())
}

/// Whether every `update` rule of `var` takes its value from outside the model
/// — a data source (`from`) or a registered handler — which is what makes the
/// array runtime serve it from the forcing buffer.
fn is_externally_fed(var: &crate::types::ModelVariable) -> bool {
    var.update.as_ref().is_some_and(|spec| {
        spec.rules().iter().all(|rule| {
            rule.value()
                .is_some_and(|v| v.from.is_some() || v.handler.is_some())
        })
    })
}

/// Whether `name` — a compiled or flattened name — or its model-local tail is a
/// key of `keys`.
fn names_key<V>(keys: &HashMap<String, V>, name: &str) -> bool {
    keys.contains_key(name)
        || name
            .split_once('.')
            .is_some_and(|(_, tail)| keys.contains_key(tail))
}

/// The declared variable a compiled or flattened `name` refers to, for a
/// diagnostic: `Model.var` in a top-level model, or a bare name one model
/// declares.
fn declared_variable<'f>(
    file: Option<&'f EsmFile>,
    name: &str,
) -> Option<&'f crate::types::ModelVariable> {
    let models = file?.models.as_ref()?;
    if let Some((model, var)) = name.split_once('.')
        && let Some(v) = models.get(model).and_then(|m| m.variables.get(var))
    {
        return Some(v);
    }
    models.values().find_map(|m| m.variables.get(name))
}

/// The `E_TREEWALK_MISSING_DATA` refusal for shaped parameter `name`, naming
/// what the document says feeds it and how a caller supplies it.
fn missing_data_error(name: &str, var: Option<&crate::types::ModelVariable>) -> SimulateError {
    let shaped = var.is_some_and(is_shaped_parameter);
    let shape = var
        .and_then(|v| v.shape.as_ref())
        .filter(|s| !s.is_empty())
        .map(|s| format!(" (shape [{}])", s.join(", ")))
        .unwrap_or_default();
    let mut feeds: Vec<String> = Vec::new();
    for rule in var
        .and_then(|v| v.update.as_ref())
        .map(|u| u.rules())
        .unwrap_or_default()
    {
        let Some(value) = rule.value() else { continue };
        if let Some(from) = &value.from {
            feeds.push(format!(
                "the data source '{}' (file_variable '{}')",
                rule.data_source().unwrap_or("?"),
                from.file_variable
            ));
        } else if let Some(h) = &value.handler {
            feeds.push(format!(
                "the registered handler '{}' (update kind '{}'), which writes it only when it fires",
                h.handler_id,
                rule.kind()
            ));
        }
    }
    let (fed, how) = if feeds.is_empty() {
        (
            "Nothing in the document gives it a value.".to_string(),
            if shaped {
                format!(
                    "pass its array in `ProblemOptions::const_arrays` under \"{name}\", or declare \
                     a `default` on it"
                )
            } else {
                format!(
                    "pass it in `ProblemOptions::p` under \"{name}\", or declare a `default` on it"
                )
            },
        )
    } else {
        (
            format!(
                "The document feeds it from {}, and no data for it was supplied at construction.",
                feeds.join(" and ")
            ),
            format!(
                "pass a provider for it in `ProblemOptions::providers` under \"{name}\", or its \
                 array in `ProblemOptions::const_arrays` (a constant snapshot), or declare a \
                 `default` on it"
            ),
        )
    };
    SimulateError::Compile(crate::compile_error::CompileError::MissingData {
        parameter: name.to_string(),
        detail: format!(
            "no data supplied for the {kind}parameter '{name}'{shape}, which declares no \
             `default`. {fed} A parameter with neither a default nor a supplied value is an error \
             when a problem is built (esm-spec §10.10). To supply it, {how}.",
            kind = if shaped { "shaped " } else { "" }
        ),
    })
}

/// The `E_TREEWALK_MISSING_INITIAL_VALUE` refusal: a state slot no default,
/// initial condition, `ic` equation or caller `u0` gives a starting value.
#[cfg(feature = "solve")]
fn refuse_missing_initial_values(
    c: &ArrayCompiled,
    u0: &HashMap<String, f64>,
) -> Result<(), SimulateError> {
    let unset = c.unset_initial_slots(u0);
    let Some(first) = unset.first() else {
        return Ok(());
    };
    let mut shown: Vec<String> = unset.iter().take(5).map(|n| format!("'{n}'")).collect();
    if unset.len() > 5 {
        shown.push("…".to_string());
    }
    Err(SimulateError::Compile(
        crate::compile_error::CompileError::MissingInitialValue {
            slot: first.clone(),
            detail: format!(
                "no starting value for {} unknown(s) ({}): the unknown declares no `default`, \
                 and no initial condition, `ic` equation or caller `u0` sets it. An unknown with \
                 no starting value is an error when a problem is built (esm-spec §11.4). To \
                 supply it, pass it in `ProblemOptions::u0` under \"{first}\", or declare a \
                 `default` on the unknown.",
                unset.len(),
                shown.join(", ")
            ),
        },
    ))
}

/// The nested [`crate::types::InlineValue`] an array denotes, row-major — the
/// spelling a shaped parameter's inline `default` takes (esm-spec §6.3).
#[cfg(not(target_arch = "wasm32"))]
fn inline_value_of(a: &ArrayD<f64>) -> crate::types::InlineValue {
    use crate::types::InlineValue;
    fn nest(shape: &[usize], values: &[f64]) -> InlineValue {
        let Some((&axis, rest)) = shape.split_first() else {
            return InlineValue::Scalar(values[0]);
        };
        let stride = rest.iter().product::<usize>();
        InlineValue::Array(
            (0..axis)
                .map(|i| nest(rest, &values[i * stride..(i + 1) * stride]))
                .collect(),
        )
    }
    let values: Vec<f64> = a.iter().copied().collect();
    if a.ndim() == 0 {
        return InlineValue::Scalar(values[0]);
    }
    nest(a.shape(), &values)
}

/// The data a caller supplied for a document's shaped parameters, captured at
/// the top of [`esm_problem`] because the build pipeline consumes
/// `const_arrays` and `build_providers`.
#[derive(Default)]
struct SuppliedData {
    /// Each `const_arrays` entry whose key names a shaped parameter of the
    /// document, as `Model.param` or as the bare `param`.
    #[cfg(not(target_arch = "wasm32"))]
    arrays: HashMap<String, ArrayD<f64>>,
    /// The keys of the run-time providers and of the build-time providers.
    provider_keys: HashMap<String, ()>,
    /// The captured arrays for parameters the document feeds from outside
    /// (a data source or a handler) and no provider serves, keyed
    /// `Model.param`: [`Self::bind_arrays`] routes them to the forcing buffer,
    /// the channel such a parameter is read through, as a constant snapshot.
    #[cfg(not(target_arch = "wasm32"))]
    forcing_arrays: std::cell::RefCell<HashMap<String, ArrayD<f64>>>,
}

impl SuppliedData {
    /// Capture what `opts` supplies. The shaped-parameter arrays are MOVED out
    /// of `opts.const_arrays` when they are the only reason the build pipeline
    /// would run on a document of several models with none selected, which the
    /// pipeline refuses: they need nothing from it there (stage 3d binds each
    /// onto its parameter). Otherwise the pipeline keeps them — it evaluates a
    /// state-free document's observed graph with them — and they are copied
    /// here.
    fn capture(opts: &mut ProblemOptions, raw: Option<&JsonValue>, file: Option<&EsmFile>) -> Self {
        #[cfg_attr(target_arch = "wasm32", allow(unused_mut))]
        let mut out = SuppliedData::default();
        #[cfg(not(target_arch = "wasm32"))]
        {
            out.provider_keys
                .extend(opts.providers.keys().map(|k| (k.clone(), ())));
            out.provider_keys
                .extend(opts.build_providers.iter().map(|(k, _)| (k.clone(), ())));
            if !opts.const_arrays.is_empty() {
                let shaped = shaped_parameter_keys(raw, file);
                let several_models =
                    opts.model_name.is_none() && model_count(raw, file).is_some_and(|n| n > 1);
                let only_trigger = (several_models || raw.is_some_and(declares_derivative))
                    && !opts.build_pipeline
                    && !opts.pushdown_rewrite
                    && opts.build_providers.is_empty()
                    && opts.const_arrays.keys().all(|k| shaped.contains(k));
                if only_trigger {
                    out.arrays = std::mem::take(&mut opts.const_arrays);
                } else {
                    for (k, a) in &opts.const_arrays {
                        if shaped.contains(k) {
                            out.arrays.insert(k.clone(), a.clone());
                        }
                    }
                }
            }
        }
        #[cfg(target_arch = "wasm32")]
        let _ = (opts, raw, file);
        out
    }

    /// Bind each captured array to the shaped parameter it names — keyed
    /// `Model.param`, or by a bare name only one model's shaped parameter
    /// carries. A scalar override in `p` outranks it, and so does a provider
    /// for a parameter the document feeds from outside. Otherwise it becomes
    /// the parameter's inline `default`, or, for a parameter fed from outside,
    /// a constant snapshot in the forcing buffer (served by
    /// [`Self::refuse_missing_array`]).
    #[cfg(not(target_arch = "wasm32"))]
    fn bind_arrays(&self, file: &mut EsmFile, p: &HashMap<String, f64>, selected: Option<&str>) {
        if self.arrays.is_empty() {
            return;
        }
        let Some(models) = file.models.as_mut() else {
            return;
        };
        // With one model selected, a bare key names that model's parameter;
        // the other models are not built, so they cannot make it ambiguous.
        let in_scope = |mname: &str| selected.is_none_or(|s| s == mname);
        let mut owners: HashMap<String, usize> = HashMap::new();
        for (_, m) in models.iter().filter(|(n, _)| in_scope(n)) {
            for (v, var) in &m.variables {
                if is_shaped_parameter(var) {
                    *owners.entry(v.clone()).or_default() += 1;
                }
            }
        }
        for (mname, m) in models.iter_mut() {
            if !in_scope(mname) {
                continue;
            }
            for (v, var) in m.variables.iter_mut() {
                if !is_shaped_parameter(var) {
                    continue;
                }
                let qualified = format!("{mname}.{v}");
                let arr = self.arrays.get(&qualified).or_else(|| {
                    if owners.get(v) == Some(&1) {
                        self.arrays.get(v)
                    } else {
                        None
                    }
                });
                let Some(arr) = arr else { continue };
                if p.contains_key(&qualified) || p.contains_key(v) {
                    continue;
                }
                if is_externally_fed(var) {
                    if !names_key(&self.provider_keys, &qualified) {
                        self.forcing_arrays
                            .borrow_mut()
                            .insert(qualified, arr.clone());
                    }
                    continue;
                }
                var.default = Some(inline_value_of(arr));
            }
        }
    }

    /// The missing-data gate for an array model: serve a forcing name no
    /// provider supplies from the caller's array, else from its declared
    /// `default`, and refuse one with neither, or a shaped parameter the build
    /// left with no value at all.
    fn refuse_missing_array(
        &self,
        c: &ArrayCompiled,
        file: Option<&EsmFile>,
        p: &HashMap<String, f64>,
    ) -> Result<(), SimulateError> {
        // The parameter as the caller addresses it: by full name (esm-spec
        // §6.6.2), which a single-model build compiles bare under its namespace.
        let full = |name: &str| match c.namespace() {
            Some(ns) if !name.starts_with(&format!("{ns}.")) => format!("{ns}.{name}"),
            _ => name.to_string(),
        };
        let forcing = c.forcing_buffer();
        for name in c.forcing_names() {
            if forcing.borrow().contains_key(name) || names_key(&self.provider_keys, name) {
                continue;
            }
            #[cfg(not(target_arch = "wasm32"))]
            if let Some(arr) = self.forcing_array(name) {
                forcing.borrow_mut().insert(name.clone(), arr);
                continue;
            }
            // No data is supplied: a declared default is the value (esm-spec
            // §6.3), and with none the parameter is missing data (§10.10).
            let mut value = match c.forcing_default(name) {
                None => {
                    return Err(missing_data_error(
                        &full(name),
                        declared_variable(file, name),
                    ));
                }
                // A default over a shape that does not resolve yet (an
                // unmaterialized derived set) has no field to fill here.
                Some(None) => continue,
                Some(Some(value)) => value,
            };
            let prec = precision::of_variable(name);
            if prec.is_f32() {
                value.mapv_inplace(|x| prec.round(x));
            }
            forcing.borrow_mut().insert(name.clone(), value);
        }
        for name in c.unvalued_shaped_params() {
            if !names_key(p, name) {
                return Err(missing_data_error(
                    &full(name),
                    declared_variable(file, name),
                ));
            }
        }
        // A scalar parameter no default, caller `p` or `distribution` gives a
        // value (esm-spec §10.10).
        for name in c.unsupplied_params(p) {
            let var = declared_variable(file, &name);
            if c.unvalued_shaped_params().contains(&name)
                || var.is_some_and(|v| v.distribution.is_some())
            {
                continue;
            }
            return Err(missing_data_error(&full(&name), var));
        }
        Ok(())
    }

    /// The captured snapshot for forcing name `name` — a compiled name, which
    /// is `Model.param` on a coupled build and the bare `param` on a
    /// single-model one.
    #[cfg(not(target_arch = "wasm32"))]
    fn forcing_array(&self, name: &str) -> Option<ArrayD<f64>> {
        let arrays = self.forcing_arrays.borrow();
        if let Some(a) = arrays.get(name) {
            return Some(a.clone());
        }
        let mut tails = arrays
            .iter()
            .filter(|(k, _)| k.split_once('.').is_some_and(|(_, tail)| tail == name));
        match (tails.next(), tails.next()) {
            (Some((_, a)), None) => Some(a.clone()),
            _ => None,
        }
    }

    /// The missing-data gate for a document with nothing to integrate, read
    /// off its flattened parameters (coupling may bind a declared parameter to
    /// another component's variable, and then it needs no value of its own).
    /// Flattens only when some parameter declares no value.
    fn refuse_missing_static(
        &self,
        file: &EsmFile,
        p: &HashMap<String, f64>,
        selected: Option<&str>,
    ) -> Result<(), SimulateError> {
        let unvalued = |v: &crate::types::ModelVariable| {
            v.var_type == crate::types::VariableType::Parameter
                && v.default.is_none()
                && v.distribution.is_none()
        };
        let any = file
            .models
            .as_ref()
            .is_some_and(|ms| ms.values().any(|m| m.variables.values().any(unvalued)));
        if !any {
            return Ok(());
        }
        let Ok(flat) = crate::flatten::flatten(file) else {
            return Ok(());
        };
        for (name, var) in &flat.parameters {
            if !unvalued(var) || names_key(p, name) {
                continue;
            }
            // With one model selected, only its parameters are built.
            if selected.is_some_and(|m| !name.starts_with(&format!("{m}."))) {
                continue;
            }
            if is_externally_fed(var) && names_key(&self.provider_keys, name) {
                continue;
            }
            return Err(missing_data_error(name, Some(var)));
        }
        Ok(())
    }
}

/// Set each 0-D parameter a value-invention node reads (inside a `skolem`,
/// `rank` or arg-witness node, or a `distinct` / keyed `faq`) to the caller's
/// `p` entry for it — keyed `Model.param` or bare — and return the `p` keys
/// used.
fn bind_value_invention_params(file: &mut EsmFile, p: &HashMap<String, f64>) -> Vec<String> {
    use crate::types::Expr;
    use std::collections::HashSet;
    fn collect(e: &Expr, inside: bool, out: &mut HashSet<String>) {
        match e {
            Expr::Variable(v) if inside => {
                out.insert(v.clone());
            }
            Expr::Operator(n) => {
                let vi = inside
                    || matches!(n.op.as_str(), "skolem" | "rank" | "argmin" | "argmax")
                    || n.distinct == Some(true)
                    || n.key.is_some();
                n.for_each_child(&mut |c| collect(c, vi, out));
            }
            _ => {}
        }
    }
    let mut used = Vec::new();
    if p.is_empty() {
        return used;
    }
    let Some(models) = file.models.as_mut() else {
        return used;
    };
    for (mname, m) in models.iter_mut() {
        let mut read: HashSet<String> = HashSet::new();
        for eq in &m.equations {
            collect(&eq.rhs, false, &mut read);
        }
        for (v, var) in m.variables.iter_mut() {
            if var.var_type != crate::types::VariableType::Parameter
                || var.shape.as_ref().is_some_and(|s| !s.is_empty())
                || !read.contains(v)
            {
                continue;
            }
            let qualified = format!("{mname}.{v}");
            let key = if p.contains_key(&qualified) {
                qualified
            } else if p.contains_key(v) {
                v.clone()
            } else {
                continue;
            };
            var.default = Some(crate::types::InlineValue::Scalar(p[&key]));
            used.push(key);
        }
    }
    used
}

/// Bake every scalar parameter a join key column reads
/// ([`crate::join::join_key_parameters`]): a caller's `p` entry lands on its
/// `default`, which the build folds into the match set, and the parameter is
/// recorded — under both its model-qualified and its bare name, whichever a
/// later [`remake`] spells — so a substitution of it is refused. Baked whether
/// or not the caller passed a value: the declared default is a build-time
/// value just the same.
fn bind_join_key_params(file: &mut EsmFile, p: &HashMap<String, f64>) -> Vec<String> {
    let mut baked = Vec::new();
    let Some(models) = file.models.as_mut() else {
        return baked;
    };
    for (mname, m) in models.iter_mut() {
        for v in crate::join::join_key_parameters(m) {
            let qualified = format!("{mname}.{v}");
            if let Some(&x) = p.get(&qualified).or_else(|| p.get(&v))
                && let Some(var) = m.variables.get_mut(&v)
            {
                var.default = Some(crate::types::InlineValue::Scalar(x));
            }
            baked.push(qualified);
            baked.push(v);
        }
    }
    baked
}

/// Whether any equation of the raw document has a time derivative on its
/// left-hand side, or the document holds a reaction system: a document with a
/// right-hand side to compile, whose caller arrays the compile binds itself.
#[cfg(not(target_arch = "wasm32"))]
fn declares_derivative(raw: &JsonValue) -> bool {
    fn has_d(v: &JsonValue) -> bool {
        match v {
            JsonValue::Object(o) => {
                o.get("op").and_then(JsonValue::as_str) == Some("D") || o.values().any(has_d)
            }
            JsonValue::Array(a) => a.iter().any(has_d),
            _ => false,
        }
    }
    if raw
        .get("reaction_systems")
        .and_then(JsonValue::as_object)
        .is_some_and(|r| !r.is_empty())
    {
        return true;
    }
    raw.get("models")
        .and_then(JsonValue::as_object)
        .is_some_and(|ms| {
            ms.values().any(|m| {
                m.get("equations")
                    .and_then(JsonValue::as_array)
                    .is_some_and(|eqs| eqs.iter().any(|e| e.get("lhs").is_some_and(has_d)))
            })
        })
}

/// How many `models` the document declares, off the raw JSON when there is one.
#[cfg(not(target_arch = "wasm32"))]
fn model_count(raw: Option<&JsonValue>, file: Option<&EsmFile>) -> Option<usize> {
    match raw {
        Some(r) => r
            .get("models")
            .and_then(JsonValue::as_object)
            .map(|m| m.len()),
        None => file.and_then(|f| f.models.as_ref()).map(|m| m.len()),
    }
}

/// Every name a caller's `const_arrays` key may use for a shaped parameter of
/// the document: `Model.param` and the bare `param`, read off the raw JSON when
/// there is one and the typed document otherwise.
#[cfg(not(target_arch = "wasm32"))]
fn shaped_parameter_keys(
    raw: Option<&JsonValue>,
    file: Option<&EsmFile>,
) -> std::collections::HashSet<String> {
    let mut out = std::collections::HashSet::new();
    if let Some(models) = raw
        .and_then(|v| v.get("models"))
        .and_then(JsonValue::as_object)
    {
        for (m, model) in models {
            let Some(vars) = model.get("variables").and_then(JsonValue::as_object) else {
                continue;
            };
            for (v, var) in vars {
                let shaped = var
                    .get("shape")
                    .and_then(JsonValue::as_array)
                    .is_some_and(|s| !s.is_empty());
                if shaped && var.get("type").and_then(JsonValue::as_str) == Some("parameter") {
                    out.insert(format!("{m}.{v}"));
                    out.insert(v.clone());
                }
            }
        }
    } else if let Some(models) = file.and_then(|f| f.models.as_ref()) {
        for (m, model) in models {
            for (v, var) in &model.variables {
                if is_shaped_parameter(var) {
                    out.insert(format!("{m}.{v}"));
                    out.insert(v.clone());
                }
            }
        }
    }
    out
}

#[cfg(not(target_arch = "wasm32"))]
#[allow(clippy::type_complexity)]
fn bind_providers(
    backend: &Backend,
    opts: &mut ProblemOptions,
    tspan: (f64, f64),
) -> Result<
    (
        Option<std::cell::RefCell<crate::provider::RefreshExecutor>>,
        std::collections::HashSet<String>,
        Vec<f64>,
    ),
    SimulateError,
> {
    if opts.providers.is_empty() {
        return Ok((None, Default::default(), Vec::new()));
    }
    let Backend::Array(compiled) = backend else {
        // Providers feed loader fields, which only exist on the array/spatial
        // runtime; a pure-scalar or static document has nowhere to put them.
        let name = opts.providers.keys().next().cloned().unwrap_or_default();
        return Err(SimulateError::ProviderError {
            name,
            details: "providers require an array/spatial model (this document has none)"
                .to_string(),
        });
    };

    // One RefreshExecutor, classifying each provider by its `refresh_times()`.
    // CONST forcings are materialized ONCE, HERE — at construction, which is
    // where §2.5.2 puts the gated fetch — and never again.
    let providers = std::mem::take(&mut opts.providers);
    let mut exec = crate::provider::RefreshExecutor::from_providers(providers);
    let forcing = compiled.forcing_buffer();
    exec.materialize_const(&forcing)
        .map_err(|e| SimulateError::ProviderError {
            name: "<const-loader>".into(),
            details: e.to_string(),
        })?;

    let discrete: std::collections::HashSet<String> = exec
        .bindings()
        .filter(|b| b.cadence == crate::cadence::Cadence::Discrete)
        .flat_map(|b| b.variables.iter().cloned())
        .collect();

    let sample_time = opts.sample_time.unwrap_or(tspan.0);
    let mut boundaries: Vec<f64> = exec.refresh_times();
    boundaries.retain(|&b| b > sample_time && b < tspan.1);

    Ok((Some(std::cell::RefCell::new(exec)), discrete, boundaries))
}

// =============================================================================
// solve
// =============================================================================

/// Run `prob` to completion (§2.5.3).
///
/// `opts` carries only per-run knobs: `alg`, `abstol`, `reltol`, `saveat`,
/// `callback`, `maxiters`. The result carries a
/// [`retcode`](crate::ReturnCode) — a caller distinguishes "ran to `tspan.1`"
/// from "stopped early, here is why" without parsing prose.
///
/// **`opts.callback` REPLACES the EsmProblem's callback set entirely.** It does
/// not append, merge or wrap. To extend rather than replace, read the set back
/// with [`callbacks`] and [`compose`] explicitly. See §2.5.4 for why
/// replacement is the safe default.
///
/// **`opts.output_observed` adds rows.** A [`Solution`] carries state rows by
/// default; a name listed there is appended as extra rows in the flat cell-key
/// spelling the state uses, which is what [`crate::derive_output_plan`] needs
/// to write that field alongside the state (RFC decision 8).
#[cfg(feature = "solve")]
pub fn solve(prob: &EsmProblem, opts: &SolveOptions) -> Result<Solution, SimulateError> {
    // Re-arm the document's working precision for the duration of this call
    // (`domain.element_type`, esm-spec §11.3); a no-op for a Float64 document.
    let _precision_guard = prob.precision.enter();
    reject_f32_integration(prob)?;
    let effective = effective_options(prob, opts);
    match &*prob.backend {
        Backend::Static(reason) => Err(SimulateError::NotDynamic {
            details: reason.clone(),
        }),
        Backend::Array(compiled) => {
            let mut insp = BuildInspection::default();
            let sink = prob.inspect.then_some(&mut insp);

            #[cfg(not(target_arch = "wasm32"))]
            let sol = {
                match (&prob.refresh, prob.discrete_forcing.is_empty()) {
                    (Some(exec), false) => {
                        let forcing = compiled.forcing_buffer();
                        let mut exec = exec.borrow_mut();
                        let refresh_fn = |t: f64| -> Result<(), SimulateError> {
                            exec.refresh_at(t, &forcing).map(|_| ()).map_err(|e| {
                                SimulateError::ProviderError {
                                    name: "<discrete-loader>".into(),
                                    details: e.to_string(),
                                }
                            })
                        };
                        compiled.solve_with_refresh_inspect(
                            prob.tspan,
                            &prob.p,
                            &prob.u0,
                            &effective,
                            sink,
                            &prob.discrete_forcing,
                            &prob.refresh_boundaries,
                            refresh_fn,
                        )?
                    }
                    _ => compiled.solve_inspect(prob.tspan, &prob.p, &prob.u0, &effective, sink)?,
                }
            };
            #[cfg(target_arch = "wasm32")]
            let sol = compiled.solve_inspect(prob.tspan, &prob.p, &prob.u0, &effective, sink)?;

            if prob.inspect {
                *prob.inspection.borrow_mut() = insp;
            }
            // One spelling per document, whichever route it took: the
            // single-model array build names its rows bare, the flattened one
            // qualifies them. See `qualify_array_names`.
            let mut sol = sol;
            if compiled.namespace().is_some() {
                // The state rows' names are the model's (kept, already
                // qualified); only the observed rows after them are spelled here.
                let states = compiled.qualified_state_names();
                let observed = sol
                    .state_variable_names
                    .split_off(states.len().min(sol.state_variable_names.len()));
                let mut names = states.to_vec();
                names.extend(qualify_array_names(compiled, &observed));
                sol.state_variable_names = names;
            }
            // Every row now carries the namespace, so none is reached through it.
            sol.metadata.namespace = None;
            Ok(sol)
        }
    }
}

/// Run the esm-spec §2.2.2 tolerance chain against `prob`'s `solver` block and
/// return `opts` with both tolerances made CONCRETE.
///
/// Caller > the document's `solver` block > binding default, per field. The
/// chain belongs to "a document is being integrated", not to one entry point,
/// so it is applied at BOTH doors a document can be integrated through: `solve`
/// (via [`effective_options`]) and `init` (§2.5.6's stepping lifecycle). A
/// stepping caller who names no tolerance must get the document's, exactly as a
/// `solve` caller does.
///
/// **Idempotent**, which is what makes applying it at both doors safe: the
/// result has `Some` in both slots, so a second pass resolves at level 1 to the
/// same values. `init` relies on this — its `step` runs each segment through
/// `solve`, which resolves again.
///
/// Deliberately does NOT touch callbacks. Callback folding is `solve`-only for
/// the same reason: [`effective_options`] is *not* idempotent in that field
/// (a second pass would wrap the already-wrapped observer and invoke every
/// callback twice per step), so `init` must resolve tolerances without it.
fn with_resolved_tolerances(prob: &EsmProblem, opts: &SolveOptions) -> SolveOptions {
    let (abstol, reltol) =
        crate::resolve_tolerances(prob.solver.as_ref(), opts.abstol, opts.reltol);
    SolveOptions {
        abstol: Some(abstol),
        reltol: Some(reltol),
        ..opts.clone()
    }
}

/// Fold the EsmProblem's callbacks (or the run's REPLACEMENT set) and the
/// extension-seam progress observer into the one per-step hook `run_solver`
/// already drives, on top of the §2.2.2 tolerance chain.
fn effective_options(prob: &EsmProblem, opts: &SolveOptions) -> SolveOptions {
    // An OUTPUT REQUEST is keyed by name, so it addresses a state an
    // `operator_compose` renaming match may have DELETED (issue #230). Resolve
    // the request list onto the surviving spellings here, where the compiled
    // artifact's map is in hand — `derive_output_plan` sees only the request and
    // the slot names, and would report the dead spelling as unknown.
    let observed = resolve_output_request(prob, &opts.output_observed);

    // esm-spec §2.2.2, resolved once here — the single point where a `solve`
    // run's options and the document meet — so every backend below sees
    // concrete tolerances and none of them has to know about the chain.
    let resolved = with_resolved_tolerances(prob, opts);

    // §2.5.4: the run's `callback` REPLACES the EsmProblem's set. It does not
    // append, merge, or wrap.
    let set = opts
        .callback
        .clone()
        .unwrap_or_else(|| prob.callbacks.clone());
    if set.is_empty() {
        return match observed {
            Some(output_observed) => SolveOptions {
                output_observed,
                ..resolved
            },
            None => resolved,
        };
    }
    let user = opts.progress.clone();
    let observer: ProgressFn = wrap_observer(set, user);
    SolveOptions {
        progress: Some(observer),
        output_observed: observed.unwrap_or_else(|| opts.output_observed.clone()),
        ..resolved
    }
}

/// The problem's flatten-time merge map (issue #230), or an empty one for a
/// backend that carries none.
fn problem_merged_renames(prob: &EsmProblem) -> &HashMap<String, String> {
    static EMPTY: std::sync::OnceLock<HashMap<String, String>> = std::sync::OnceLock::new();
    match &*prob.backend {
        Backend::Array(c) => c.merged_renames(),
        Backend::Static(_) => EMPTY.get_or_init(HashMap::new),
    }
}

/// Rewrite an output request off names an `operator_compose` renaming match
/// DELETED. `None` when nothing changes, so the common path clones nothing.
///
/// Not gated on `solve`: its one caller, `effective_options`, is not either.
fn resolve_output_request(prob: &EsmProblem, requested: &[String]) -> Option<Vec<String>> {
    if requested.is_empty() {
        return None;
    }
    let renames = problem_merged_renames(prob);
    if renames.is_empty() || !requested.iter().any(|n| renames.contains_key(n.as_str())) {
        return None;
    }
    Some(
        requested
            .iter()
            .map(|n| {
                renames
                    .get(n.as_str())
                    .cloned()
                    .unwrap_or_else(|| n.clone())
            })
            .collect(),
    )
}

/// The document's §2.2 `solver` block, read out of a raw document.
///
/// Used at CONSTRUCTION, for the input carriers that reach `esm_problem` as raw
/// JSON (a path, a `Json` value, or a typed file the build pipeline
/// re-serialized). The typed carriers read `EsmFile::solver` directly; either
/// way the block lands on [`EsmProblem::solver`], which is what the §2.2.2
/// chain consults at solve time.
fn document_solver(doc: &JsonValue) -> Option<crate::Solver> {
    let block = doc.get("solver")?;
    // §2.2: an EMPTY block means what absence means, so it normalizes away
    // rather than becoming a `Solver` with nothing set.
    crate::solver::normalize_empty(serde_json::from_value(block.clone()).ok()?)
}

#[cfg(not(target_arch = "wasm32"))]
fn wrap_observer(set: CallbackSet, user: Option<ProgressFn>) -> ProgressFn {
    std::sync::Arc::new(move |p: &Progress<'_>| {
        let a = set.invoke(p);
        let b = user.as_ref().map(|f| f(p)).unwrap_or(Flow::Continue);
        if matches!(a, Flow::Cancel) || matches!(b, Flow::Cancel) {
            Flow::Cancel
        } else {
            Flow::Continue
        }
    })
}

#[cfg(target_arch = "wasm32")]
fn wrap_observer(set: CallbackSet, user: Option<ProgressFn>) -> ProgressFn {
    std::sync::Arc::new(move |p: &Progress<'_>| {
        let a = set.invoke(p);
        let b = user.as_ref().map(|f| f(p)).unwrap_or(Flow::Continue);
        if matches!(a, Flow::Cancel) || matches!(b, Flow::Cancel) {
            Flow::Cancel
        } else {
            Flow::Continue
        }
    })
}

// =============================================================================
// Stepping
// =============================================================================

/// A stepping integrator (§2.5.6) — the same lifecycle [`solve`] performs
/// internally, exposed for callers that need to interleave their own work with
/// the integration: the coupling driver in a host model, an interactive
/// session, a progress UI.
///
/// # What a "step" is here
///
/// [`Integrator::step`] advances to the next entry of the **output grid**, not
/// to the next internal solver step. Each advance runs the configured
/// algorithm over `[t_now, t_next]` from the current state, exactly the way the
/// array runtime's segmented provider-refresh driver already advances across a
/// refresh boundary. That means the adaptive controller restarts at every
/// interleave point, so a run stepped in N pieces is not step-for-step
/// identical to the same run solved in one call — it is the same trajectory to
/// tolerance, at some extra cost per boundary. Choose the grid accordingly:
/// coarse when you only need to interleave, fine when you need control.
///
/// (`diffsol`'s solver borrows its `OdeSolverProblem`, so a EsmProblem-owning
/// integrator that also owns a live solver would be self-referential. Restarting
/// per grid interval is the safe-Rust way to expose the lifecycle, and it is
/// the mechanism already in production here for segmented refresh.)
#[cfg(feature = "solve")]
pub struct Integrator<'a> {
    prob: &'a EsmProblem,
    opts: SolveOptions,
    grid: Vec<f64>,
    next: usize,
    t: f64,
    u: HashMap<String, f64>,
    time: Vec<f64>,
    state: Vec<Vec<f64>>,
    names: Vec<String>,
    retcode: crate::simulate::ReturnCode,
    metadata: crate::simulate::SolutionMetadata,
}

#[cfg(feature = "solve")]
impl std::fmt::Debug for Integrator<'_> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("Integrator")
            .field("t", &self.t)
            .field("remaining", &(self.grid.len().saturating_sub(self.next)))
            .field("retcode", &self.retcode)
            .finish()
    }
}

/// Whether an [`Integrator`] has more to do.
#[cfg(feature = "solve")]
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum StepStatus {
    /// The integrator advanced and has further grid points to cover.
    Advanced,
    /// The integrator reached `tspan.1`, or stopped early — read
    /// [`Integrator::retcode`].
    Done,
}

/// Build a stepping integrator for `prob` (§2.5.6).
///
/// The step grid is `opts.saveat` when the caller supplied one, else 100 evenly
/// spaced points across `tspan`.
///
/// The esm-spec §2.2.2 tolerance chain runs HERE, not only at `solve`: the
/// integrator this returns carries the document's `solver.abstol` / `reltol`
/// whenever the caller named none, so stepping a document honours its declared
/// numerics exactly as solving it does.
#[cfg(feature = "solve")]
pub fn init<'a>(
    prob: &'a EsmProblem,
    opts: &SolveOptions,
) -> Result<Integrator<'a>, SimulateError> {
    // Re-arm the document's working precision for the duration of this call
    // (`domain.element_type`, esm-spec §11.3); a no-op for a Float64 document.
    let _precision_guard = prob.precision.enter();
    reject_f32_integration(prob)?;
    if let Backend::Static(reason) = &*prob.backend {
        return Err(SimulateError::NotDynamic {
            details: reason.clone(),
        });
    }
    let (t0, t_end) = prob.tspan;
    let grid = match &opts.saveat {
        Some(g) if g.len() >= 2 => g.clone(),
        _ => {
            let n = 100usize;
            (0..=n)
                .map(|i| t0 + (t_end - t0) * (i as f64) / (n as f64))
                .collect()
        }
    };
    Ok(Integrator {
        prob,
        // esm-spec §2.2.2: caller > the document's `solver` block > binding
        // default. Resolved once, at construction, so the integrator HOLDS the
        // effective tolerances rather than leaving them to be rediscovered —
        // and so `step`'s per-segment `solve` sees them as explicit call-site
        // values (level 1) that resolve to themselves.
        opts: with_resolved_tolerances(prob, opts),
        grid,
        next: 0,
        t: t0,
        u: prob.u0.clone(),
        time: Vec::new(),
        state: Vec::new(),
        names: prob.state_variable_names(),
        retcode: crate::simulate::ReturnCode::Success,
        metadata: crate::simulate::SolutionMetadata::default(),
    })
}

#[cfg(feature = "solve")]
impl Integrator<'_> {
    /// The independent-variable value reached so far.
    pub fn t(&self) -> f64 {
        self.t
    }

    /// Why the integration stopped, once it has.
    pub fn retcode(&self) -> crate::simulate::ReturnCode {
        self.retcode
    }

    /// The current state, keyed by variable name.
    pub fn u(&self) -> &HashMap<String, f64> {
        &self.u
    }

    /// Advance to the next grid point. Returns [`StepStatus::Done`] at the end
    /// of `tspan` or when the run stopped early.
    pub fn step(&mut self) -> Result<StepStatus, SimulateError> {
        // Re-arm the document's working precision (esm-spec §11.3); a stepping
        // caller reaches the RHS without passing back through `solve`.
        let _precision_guard = self.prob.precision.enter();
        if !self.retcode.is_success() {
            return Ok(StepStatus::Done);
        }
        // Skip grid points at or behind where we already are.
        while self.next < self.grid.len() && self.grid[self.next] <= self.t {
            self.next += 1;
        }
        if self.next >= self.grid.len() {
            return Ok(StepStatus::Done);
        }
        let t_next = self.grid[self.next];

        let seg = remake(
            self.prob,
            &Remake {
                u0: self.u.clone(),
                tspan: Some((self.t, t_next)),
                ..Default::default()
            },
        )?;
        let seg_opts = SolveOptions {
            saveat: Some(vec![t_next]),
            ..self.opts.clone()
        };
        let sol = solve(&seg, &seg_opts)?;

        self.retcode = sol.retcode;
        self.metadata.alg = sol.metadata.alg.clone();
        self.metadata.n_rhs_calls += sol.metadata.n_rhs_calls;
        self.metadata.n_jacobian_calls += sol.metadata.n_jacobian_calls;
        self.metadata.n_accepted_steps += sol.metadata.n_accepted_steps;
        self.metadata.n_rejected_steps += sol.metadata.n_rejected_steps;
        self.metadata.tape_fallbacks = sol.metadata.tape_fallbacks.clone();

        if self.state.is_empty() {
            self.names = sol.state_variable_names.clone();
            self.state = vec![Vec::new(); self.names.len()];
        }
        if let Some(k) = sol.time.len().checked_sub(1) {
            self.time.push(sol.time[k]);
            for (r, row) in sol.state.iter().enumerate() {
                if let Some(v) = row.get(k) {
                    if r < self.state.len() {
                        self.state[r].push(*v);
                    }
                    if let Some(name) = self.names.get(r) {
                        self.u.insert(name.clone(), *v);
                    }
                }
            }
            self.t = sol.time[k];
        }
        self.next += 1;
        if self.next >= self.grid.len() || !self.retcode.is_success() {
            Ok(StepStatus::Done)
        } else {
            Ok(StepStatus::Advanced)
        }
    }

    /// Run to completion and take the accumulated [`Solution`]. This is the
    /// `solve!` of the SciML lifecycle; it is spelled out because `!` is not a
    /// legal suffix in Rust.
    pub fn solve_to_completion(&mut self) -> Result<Solution, SimulateError> {
        while matches!(self.step()?, StepStatus::Advanced) {}
        Ok(Solution {
            time: std::mem::take(&mut self.time),
            state: std::mem::take(&mut self.state),
            state_variable_names: self.names.clone(),
            retcode: self.retcode,
            metadata: self.metadata.clone(),
        })
    }
}

/// Advance `integrator` by one grid point (§2.5.6's `step!`).
#[cfg(feature = "solve")]
pub fn step(integrator: &mut Integrator<'_>) -> Result<StepStatus, SimulateError> {
    integrator.step()
}

/// Run `integrator` to completion (§2.5.6's `solve!`).
#[cfg(feature = "solve")]
pub fn solve_to_completion(integrator: &mut Integrator<'_>) -> Result<Solution, SimulateError> {
    integrator.solve_to_completion()
}

// ===========================================================================
// Unit tests
// ===========================================================================

/// The esm-spec §2.2.2 chain at the STEPPING door.
///
/// These live in-crate rather than in `tests/solver_block.rs` because what has
/// to be pinned is that the resolved tolerances land ON the integrator — a
/// private field. An observational test (does a looser tolerance change the
/// trajectory?) would pass either way for a well-conditioned decay and is not
/// the property; this is.
#[cfg(all(test, feature = "solve", not(target_arch = "wasm32")))]
mod stepping_tolerance_tests {
    use super::*;
    use crate::simulate::{DEFAULT_ABSTOL, DEFAULT_RELTOL};

    fn doc(solver: &str) -> serde_json::Value {
        serde_json::from_str(&format!(
            r#"{{"esm":"1.1.0","metadata":{{"name":"T"}},{solver}"models":{{"M":{{
                 "variables":{{"y":{{"type":"unknown","units":"m","default":1.0}}}},
                 "equations":[{{"lhs":{{"op":"D","args":["y"],"wrt":"t"}},
                                "rhs":{{"op":"neg","args":["y"]}}}}]}}}}}}"#
        ))
        .expect("fixture parses")
    }

    fn problem(solver: &str) -> EsmProblem {
        let json = doc(solver);
        esm_problem(
            ProblemInput::Json(&json),
            (0.0, 1.0),
            ProblemOptions::default(),
        )
        .expect("builds")
    }

    /// A stepping caller who names no tolerance gets the DOCUMENT's, not the
    /// binding default.
    ///
    /// `init` used to store the caller's `SolveOptions` verbatim. The effect was
    /// MASKED in this binding — `step` remakes each segment and re-enters
    /// `solve`, and `remake` carries `solver` across, so the chain still ran, one
    /// layer down and once per segment. That is a coincidence of the current
    /// segment-per-step implementation, not the contract: the integrator did not
    /// hold the tolerances it was running at, and a step path that does not
    /// re-enter `solve` would silently drop the document. §2.2.2 puts the
    /// resolution at the DOOR, so this pins it at the door.
    #[test]
    fn init_resolves_the_documents_tolerances() {
        let prob = problem(r#""solver":{"abstol":1e-8,"reltol":1e-6},"#);
        let integ = init(&prob, &SolveOptions::default()).expect("init");
        assert_eq!(integ.opts.abstol, Some(1e-8));
        assert_eq!(integ.opts.reltol, Some(1e-6));
    }

    /// Level 1 still wins: an explicitly-passed option beats the document.
    #[test]
    fn an_explicit_call_site_option_still_beats_the_document() {
        let prob = problem(r#""solver":{"abstol":1e-8,"reltol":1e-6},"#);
        let opts = SolveOptions {
            abstol: Some(1e-12),
            reltol: Some(1e-11),
            ..SolveOptions::default()
        };
        let integ = init(&prob, &opts).expect("init");
        assert_eq!(integ.opts.abstol, Some(1e-12));
        assert_eq!(integ.opts.reltol, Some(1e-11));
    }

    /// Per FIELD, and level 3 for whatever the document leaves unsaid.
    #[test]
    fn each_tolerance_falls_through_on_its_own() {
        let only_reltol = problem(r#""solver":{"reltol":1e-9},"#);
        let integ = init(&only_reltol, &SolveOptions::default()).expect("init");
        assert_eq!(integ.opts.reltol, Some(1e-9));
        assert_eq!(integ.opts.abstol, Some(DEFAULT_ABSTOL));

        let no_block = problem("");
        let integ = init(&no_block, &SolveOptions::default()).expect("init");
        assert_eq!(integ.opts.abstol, Some(DEFAULT_ABSTOL));
        assert_eq!(integ.opts.reltol, Some(DEFAULT_RELTOL));
    }

    /// The chain is idempotent, which is what lets `init` and the per-segment
    /// `solve` inside `step` both run it without the second pass changing the
    /// answer or displacing the document.
    #[test]
    fn resolving_twice_is_resolving_once() {
        let prob = problem(r#""solver":{"abstol":1e-8,"reltol":1e-6},"#);
        let once = with_resolved_tolerances(&prob, &SolveOptions::default());
        let twice = with_resolved_tolerances(&prob, &once);
        assert_eq!((twice.abstol, twice.reltol), (once.abstol, once.reltol));
        assert_eq!((twice.abstol, twice.reltol), (Some(1e-8), Some(1e-6)));
    }
}