earthsci-ast 0.6.0

//! earthsci-ast CLI Tool
//!
//! Command-line interface for working with ESM files.
//!
//! The `esm` binary target is the thin `esm.rs` next to this file; this is its
//! whole implementation, compiled for native targets only.

use clap::{Parser, Subcommand, ValueEnum};
use earthsci_ast::extension::analysis::{
    collect_unit_types, collect_variables, contains_common_subexpressions,
    contains_expensive_operations, contains_redundant_operations, count_expression_nodes,
    count_numerical_values, count_operations, expression_depth, expressions_numerically_equal,
    find_longest_dependency_chain, find_strongly_connected_components,
};
use earthsci_ast::{
    EquationDimensionVerdict, ExpressionGraphOptions, StructuralErrorCode, Unit, UnitError,
    component_exists, component_graph, equation_dimension_verdict, expression_graph,
    expression_graph_with_options, model_unit_env, parse_unit, stoichiometric_matrix, to_json,
    to_json_compact, validate, validate_text,
};
use std::collections::{BTreeMap, HashMap};
use std::fs;
use std::path::PathBuf;

#[derive(Parser)]
#[command(name = "esm")]
#[command(about = "A CLI tool for EarthSciML Abstract Syntax Tree Format")]
#[command(version)]
struct Cli {
    #[command(subcommand)]
    command: Commands,
}

#[derive(Subcommand)]
enum Commands {
    /// Validate an ESM file
    Validate {
        /// Path to the ESM file to validate
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Show detailed validation output
        #[arg(long)]
        verbose: bool,
    },
    /// Pretty-print mathematical expressions in an ESM file
    #[command(visible_alias = "display")]
    Pretty {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Output format (unicode, latex, ascii)
        #[arg(short, long, default_value = "unicode")]
        format: DisplayFormat,
    },
    /// Extract a single component from an ESM file
    Extract {
        /// Input ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Component ID to extract
        #[arg(long)]
        component: String,
        /// Output file (or stdout if not specified)
        #[arg(short, long)]
        output: Option<PathBuf>,
    },
    /// Compare two ESM files semantically
    Diff {
        /// First ESM file
        #[arg(value_name = "FILE1")]
        file1: PathBuf,
        /// Second ESM file
        #[arg(value_name = "FILE2")]
        file2: PathBuf,
    },
    /// Generate stoichiometric matrix for a reaction system
    Stoich {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Reaction system ID
        #[arg(long)]
        system: String,
    },
    /// Generate system or expression graphs
    Graph {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Graph level (component or expression)
        #[arg(long, default_value = "component")]
        level: GraphLevel,
        /// Output format (dot, mermaid, json)
        #[arg(short, long, default_value = "dot")]
        format: GraphFormat,
        /// System ID (for expression-level graphs)
        #[arg(long)]
        system: Option<String>,
    },
    /// Convert an ESM file to a different format
    Convert {
        /// Input ESM file
        #[arg(value_name = "INPUT")]
        input: PathBuf,
        /// Output file (or stdout if not specified)
        #[arg(short, long, value_name = "OUTPUT")]
        output: Option<PathBuf>,
        /// Target format (json, compact-json)
        #[arg(long, default_value = "json")]
        to: ConvertFormat,
    },
    /// Analyze system properties and metrics
    Analyze {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Analysis type (all, structure, complexity, coupling)
        #[arg(long, default_value = "all")]
        analysis_type: AnalysisType,
    },
    /// Run basic simulations
    Simulate {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Simulation duration
        #[arg(long, default_value = "10.0")]
        time: f64,
        /// Output file for results
        #[arg(short, long)]
        output: Option<PathBuf>,
        /// Observed field to write alongside the state; repeat for several.
        /// Output is the state PLUS the fields named here, never every
        /// observed by default (streaming-output-sinks RFC decision 8).
        ///
        /// On a document with nothing to integrate this selects the COLUMNS of
        /// the single evaluation; omitted, every field the build materialized
        /// is written.
        #[arg(long = "observed", value_name = "NAME")]
        observed: Vec<String>,
        /// Build this model when the document holds several.
        #[arg(long = "model", value_name = "NAME")]
        model: Option<String>,
        /// Shape of --output: `flat` is the flat state rows keyed by cell key;
        /// `grid` is the derived output plan — dimension-labeled, row-major
        /// [time, …spatial] arrays with CF coordinates and attributes; `csv`
        /// is one row per index tuple over the shape the written fields share
        /// — a RELATIONAL document's rows, for a row-by-row comparator.
        #[arg(long, default_value = "flat")]
        format: SimulateFormat,
        /// WHICH strategy builds the right-hand side (API_SPEC §5.8's closed
        /// vocabulary): `native` (the default — the tape, for every document,
        /// and strict: a rule it cannot express is a build error naming the
        /// rule), `interpreter` (the reference: every fast tier off, no
        /// performance promise), or `xla` (the emitter over the tape, run
        /// through XLA; needs a build with the `xla` feature and a usable
        /// `xla_extension`, and refuses with `compiler_unavailable`
        /// otherwise). `mtk` and `sympy` belong to the Julia and Python
        /// bindings and are always `compiler_unavailable` here.
        #[arg(long, value_name = "NAME")]
        compiler: Option<String>,
    },
    /// Show information about an ESM file
    Info {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
    },

    // ANALYSIS COMMANDS
    /// Dimensional analysis report
    Units {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Check dimensional consistency
        #[arg(long)]
        check: bool,
    },
    /// Coupling dependency analysis
    CouplingAnalysis {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Analysis depth (shallow, deep)
        #[arg(long, default_value = "shallow")]
        depth: CouplingDepth,
    },
    /// Memory/time analysis
    PerformanceProfile {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Profiling type (memory, time, both)
        #[arg(long, default_value = "both")]
        profile_type: ProfileType,
    },
    /// Model difference report
    Compare {
        /// First ESM file
        #[arg(value_name = "FILE1")]
        file1: PathBuf,
        /// Second ESM file
        #[arg(value_name = "FILE2")]
        file2: PathBuf,
        /// Comparison type (semantic, structural, numerical)
        #[arg(long, default_value = "semantic")]
        comparison_type: ComparisonType,
    },
    /// Expression optimization suggestions
    Optimize {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Optimization type (expression, structure, performance)
        #[arg(long, default_value = "expression")]
        opt_type: OptimizeTarget,
    },

    // DEVELOPMENT COMMANDS
    /// Create new ESM project template
    Init {
        /// Project name
        #[arg(value_name = "NAME")]
        name: String,
        /// Template type (minimal, atmospheric, ecosystem, coupling)
        #[arg(long, default_value = "minimal")]
        template: InitTemplate,
    },
    /// Batch validation for test suites
    ValidateFixtures {
        /// Directory containing ESM files
        #[arg(value_name = "DIRECTORY")]
        directory: PathBuf,
        /// Recursive directory traversal
        #[arg(short, long)]
        recursive: bool,
    },
    /// Performance benchmarking
    Benchmark {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Benchmark type (parse, validate, simulate)
        #[arg(long, default_value = "all")]
        bench_type: BenchType,
        /// Number of iterations
        #[arg(long, default_value = "100")]
        iterations: usize,
    },
    /// JSON schema compliance
    SchemaCheck {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Schema version to check against
        #[arg(long)]
        schema_version: Option<String>,
    },
    /// Test load/save fidelity
    RoundTrip {
        /// Path to the ESM file
        #[arg(value_name = "FILE")]
        file: PathBuf,
        /// Number of round trips
        #[arg(long, default_value = "1")]
        rounds: usize,
    },
    /// Run the inline §6.6 tests of every named file and report pass/fail
    #[command(visible_alias = "run-tests")]
    Test {
        /// ESM files, and/or directories searched recursively for `.esm` files
        #[arg(value_name = "PATH", required = true)]
        paths: Vec<PathBuf>,
        /// Only run the tests of this model
        #[arg(long, value_name = "NAME")]
        model: Option<String>,
        /// Only run tests whose id contains this substring
        #[arg(long, value_name = "SUBSTRING")]
        filter: Option<String>,
        /// Solver family for the inline-test integrations
        #[arg(long, default_value = "bdf")]
        solver: SolverAlg,
        /// Relative SOLVER tolerance — how accurately each test is integrated,
        /// NOT the §6.6.4 tolerance its assertions are judged against.
        /// Unset resolves per esm-spec §2.2.2: the document's `solver.reltol`
        /// if it declares one, else the runner default.
        #[arg(long)]
        reltol: Option<f64>,
        /// Absolute SOLVER tolerance — how accurately each test is integrated,
        /// NOT the §6.6.4 tolerance its assertions are judged against.
        /// Unset resolves per esm-spec §2.2.2, as `--reltol`.
        #[arg(long)]
        abstol: Option<f64>,
        /// Report every assertion, not just the summary table
        #[arg(short, long)]
        verbose: bool,
        /// WHICH strategy builds each test's right-hand side (API_SPEC §5.8's
        /// closed vocabulary). Unset is `native`, which is STRICT: a document
        /// whose rules the tape cannot lower fails with the refusal rather
        /// than running on a slower path. Pass `interpreter` to run the
        /// reference evaluator, which refuses nothing it can evaluate.
        #[arg(long, value_name = "NAME")]
        compiler: Option<String>,
    },
    /// Run cross-language conformance tests and write results.json to OUT_DIR
    ConformanceTest {
        /// Directory to write results.json into
        #[arg(value_name = "OUT_DIR")]
        out_dir: PathBuf,
        /// Shared corpus manifest (scripts/conformance_corpus.py). Defaults to
        /// $ESM_CONFORMANCE_MANIFEST, then <OUT_DIR>/../corpus_manifest.json.
        #[arg(value_name = "MANIFEST")]
        manifest: Option<PathBuf>,
    },
}

// === Enumerated flag values ================================================
//
// Typed `ValueEnum` flags: clap itself rejects an unknown value (with its
// standard "invalid value … possible values …" usage error, exit 2), so the
// handlers match exhaustively instead of each carrying a reject-unknown arm.

/// `graph --level`.
#[derive(Clone, Copy, ValueEnum)]
enum GraphLevel {
    Component,
    Expression,
}

/// `graph --format` (§4.8.3 rendering formats).
#[derive(Clone, Copy, ValueEnum)]
enum GraphFormat {
    Dot,
    Mermaid,
    Json,
}

/// `simulate --format`.
#[derive(Clone, Copy, PartialEq, Eq, ValueEnum)]
enum SimulateFormat {
    /// The flat solver state: one row per cell key, no axis metadata.
    Flat,
    /// The `derive_output_plan` result: one entry per emergent grid, each with
    /// its dimensions, CF dimension coordinates and row-major variable arrays.
    Grid,
    /// One row per index tuple: the index columns `i1…in`, then one column per
    /// field, over the shape every written field shares. This is the shape of a
    /// RELATIONAL document — rows computed from rows — and the shape a
    /// row-by-row comparator reads. Only for a document evaluated once; a
    /// trajectory is not a row set.
    Csv,
}

/// `convert --to`.
#[derive(Clone, Copy, ValueEnum)]
enum ConvertFormat {
    Json,
    CompactJson,
}

/// `analyze --analysis-type`.
#[derive(Clone, Copy, ValueEnum)]
enum AnalysisType {
    All,
    Structure,
    Complexity,
    Coupling,
}

/// `coupling-analysis --depth`.
#[derive(Clone, Copy, ValueEnum)]
enum CouplingDepth {
    Shallow,
    Deep,
}

/// `performance-profile --profile-type`.
#[derive(Clone, Copy, ValueEnum)]
enum ProfileType {
    Memory,
    Time,
    Both,
}

/// `compare --comparison-type`.
#[derive(Clone, Copy, ValueEnum)]
enum ComparisonType {
    Semantic,
    Structural,
    Numerical,
}

/// `optimize --opt-type`.
#[derive(Clone, Copy, ValueEnum)]
enum OptimizeTarget {
    Expression,
    Structure,
    Performance,
}

/// `init --template` (see the TEMPLATE_* consts).
#[derive(Clone, Copy, ValueEnum)]
enum InitTemplate {
    Minimal,
    Atmospheric,
    Ecosystem,
    Coupling,
}

/// `benchmark --bench-type`.
#[derive(Clone, Copy, ValueEnum)]
enum BenchType {
    All,
    Parse,
    Validate,
    Simulate,
}

/// `test --solver`. [`earthsci_ast::Alg`] is a library type and cannot derive
/// `clap::ValueEnum` here, so the CLI spelling of the solver family lives in
/// this shim; `Alg::from_name` remains the mapping for string hosts.
#[derive(Clone, Copy, ValueEnum)]
enum SolverAlg {
    Bdf,
    Sdirk,
    Erk,
}

impl From<SolverAlg> for earthsci_ast::Alg {
    fn from(value: SolverAlg) -> Self {
        match value {
            SolverAlg::Bdf => earthsci_ast::Alg::Bdf,
            SolverAlg::Sdirk => earthsci_ast::Alg::Sdirk,
            SolverAlg::Erk => earthsci_ast::Alg::Erk,
        }
    }
}

/// The kebab-case CLI spelling of a `ValueEnum` value, for the commands that
/// echo a flag back in their output.
fn arg_name<T: ValueEnum>(value: &T) -> String {
    value
        .to_possible_value()
        .expect("no skipped variants")
        .get_name()
        .to_string()
}

// === Exit discipline =======================================================

/// A deliberate CLI failure (exit 1) raised after a handler has shaped its own
/// diagnostics. [`main`] prints `message` — if any — to stderr verbatim, with
/// no `Error:` prefix and no Debug quoting; `None` means the handler already
/// printed everything (typically a report on stdout). Ordinary propagated
/// errors keep the `Error: {:?}` rendering a `Result` main gives them.
#[derive(Debug)]
struct CliFailure {
    message: Option<String>,
}

impl std::fmt::Display for CliFailure {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.write_str(self.message.as_deref().unwrap_or("command failed"))
    }
}

impl std::error::Error for CliFailure {}

/// Exit-1 failure carrying a message for stderr.
fn fail(message: impl Into<String>) -> Box<dyn std::error::Error> {
    Box::new(CliFailure {
        message: Some(message.into()),
    })
}

/// Exit-1 failure whose diagnostics the handler already printed.
fn fail_silent() -> Box<dyn std::error::Error> {
    Box::new(CliFailure { message: None })
}

fn analyze_model_units(esm_file: &earthsci_ast::EsmFile) {
    if let Some(ref models) = esm_file.models {
        for (model_id, model) in models {
            println!("Model: {model_id}");

            // Analyze variable units
            let mut unit_count = 0;
            for (var_name, var) in &model.variables {
                if let Some(ref units) = var.units {
                    println!("  Variable {var_name}: {var_name} [{units}]");
                    unit_count += 1;
                } else {
                    println!("  Variable {var_name}: {var_name} [no units specified]");
                }
            }

            println!(
                "  {} out of {} variables have units specified",
                unit_count,
                model.variables.len()
            );

            for (i, equation, verdict) in model_equation_verdicts(esm_file, model) {
                let label = match equation_lhs_label(&equation.lhs) {
                    Some(lhs) => format!("Equation {} ({lhs})", i + 1),
                    None => format!("Equation {}", i + 1),
                };
                match verdict {
                    EquationDimensionVerdict::Consistent(dim) => {
                        println!("  {label}: consistent [{dim}]");
                    }
                    EquationDimensionVerdict::Mismatch(messages) => {
                        println!("  {label}: MISMATCH — {}", messages.join("; "));
                    }
                    EquationDimensionVerdict::NotChecked(reason) => {
                        println!("  {label}: NOT CHECKED — {reason}");
                    }
                }
            }
        }
    }
}

/// The dimensional analyser's verdict on each of `model`'s equations, judged
/// against the same unit environment `validate` uses.
fn model_equation_verdicts<'m>(
    esm_file: &earthsci_ast::EsmFile,
    model: &'m earthsci_ast::Model,
) -> Vec<(usize, &'m earthsci_ast::Equation, EquationDimensionVerdict)> {
    let (env, _) = model_unit_env(esm_file, model);
    model
        .equations
        .iter()
        .enumerate()
        .map(|(i, eq)| (i, eq, equation_dimension_verdict(eq, &env)))
        .collect()
}

/// A short name for an equation's left-hand side: `x` or `D(x)`.
fn equation_lhs_label(lhs: &earthsci_ast::Expr) -> Option<String> {
    match lhs {
        earthsci_ast::Expr::Variable(name) => Some(name.clone()),
        earthsci_ast::Expr::Operator(op) if op.op == "D" => match op.args.first() {
            Some(earthsci_ast::Expr::Variable(name)) => Some(format!("D({name})")),
            _ => None,
        },
        _ => None,
    }
}

fn analyze_reaction_system_units(
    esm_file: &earthsci_ast::EsmFile,
    validation: &earthsci_ast::ValidationResult,
) {
    if let Some(ref reaction_systems) = esm_file.reaction_systems {
        for (rs_id, rs) in reaction_systems {
            println!("Reaction System: {rs_id}");

            // Analyze species units
            let mut species_with_units = 0;
            for (species_name, species) in &rs.species {
                if let Some(ref units) = species.units {
                    println!("  Species {species_name}: [{units}]");
                    species_with_units += 1;
                } else {
                    println!("  Species {species_name}: [no units]");
                }
            }

            println!(
                "  {} out of {} species have units specified",
                species_with_units,
                rs.species.len()
            );

            for (i, verdict) in reaction_rate_verdicts(rs_id, rs, validation) {
                match verdict {
                    RateVerdict::Consistent => println!(
                        "  Reaction {} rate: no stoichiometric mismatch found",
                        i + 1
                    ),
                    RateVerdict::Mismatch(message) => {
                        println!("  Reaction {} rate: MISMATCH — {message}", i + 1);
                    }
                    RateVerdict::NotChecked(reason) => {
                        println!("  Reaction {} rate: NOT CHECKED — {reason}", i + 1);
                    }
                }
            }
        }
    }
}

/// The dimensional analyser's verdict on one reaction's rate expression.
enum RateVerdict {
    /// The rate's dimension resolved and no stoichiometric mismatch was found.
    Consistent,
    /// A provable mismatch: either the stoichiometric one `validate` reports, or
    /// an inconsistency between the rate expression's own operands.
    Mismatch(String),
    /// The rate's dimension could not be determined, so nothing was compared.
    NotChecked(String),
}

/// Judge every reaction's rate in `rs`, paired with its 0-based index.
///
/// The stoichiometric rate check lives in `validate`; its findings are keyed by
/// the reaction's JSON Pointer, so they are read back out of `validation`
/// rather than recomputed. A rate whose own operands are provably incompatible
/// is a MISMATCH too — the analyser did reach a verdict — and only an
/// undeterminable dimension is "not checked".
fn reaction_rate_verdicts(
    rs_id: &str,
    rs: &earthsci_ast::ReactionSystem,
    validation: &earthsci_ast::ValidationResult,
) -> Vec<(usize, RateVerdict)> {
    let env: HashMap<String, Unit> = rs
        .species
        .iter()
        .map(|(name, s)| (name, &s.units))
        .chain(rs.parameters.iter().map(|(name, p)| (name, &p.units)))
        .filter_map(|(name, units)| {
            let unit = parse_unit(units.as_deref()?).ok()?;
            Some((name.clone(), unit))
        })
        .collect();
    rs.reactions
        .iter()
        .enumerate()
        .map(|(i, reaction)| {
            let path = format!("/reaction_systems/{rs_id}/reactions/{i}");
            let mismatches: Vec<&str> = validation
                .structural_errors
                .iter()
                .filter(|e| {
                    e.path == path && matches!(e.code, StructuralErrorCode::UnitInconsistency)
                })
                .map(|e| e.message.as_str())
                .collect();
            if !mismatches.is_empty() {
                return (i, RateVerdict::Mismatch(mismatches.join("; ")));
            }
            match Unit::propagate(&reaction.rate, &env) {
                Ok(_) => (i, RateVerdict::Consistent),
                Err(UnitError::DimensionMismatch(message)) => (i, RateVerdict::Mismatch(message)),
                Err(e) => (i, RateVerdict::NotChecked(e.to_string())),
            }
        })
        .collect()
}

fn analyze_expression_optimization(esm_file: &earthsci_ast::EsmFile) {
    println!("\n=== EXPRESSION OPTIMIZATION ANALYSIS ===");

    let mut suggestions: Vec<String> = Vec::new();
    let mut complex_expressions = 0;
    let mut redundant_operations = 0;

    // Analyze model expressions
    if let Some(ref models) = esm_file.models {
        for (model_id, model) in models {
            println!("Model: {model_id}");

            for (i, equation) in model.equations.iter().enumerate() {
                // Check expression complexity
                let lhs_depth = expression_depth(&equation.lhs);
                let rhs_depth = expression_depth(&equation.rhs);
                let max_depth = lhs_depth.max(rhs_depth);

                if max_depth > 5 {
                    complex_expressions += 1;
                    suggestions.push(format!(
                        "  Eq {}: High complexity (depth {}), consider factoring",
                        i + 1,
                        max_depth
                    ));
                }

                // Check for redundant operations
                if contains_redundant_operations(&equation.rhs) {
                    redundant_operations += 1;
                    suggestions.push(format!("  Eq {}: Contains redundant operations", i + 1));
                }

                // Check for common subexpressions
                if contains_common_subexpressions(&equation.lhs, &equation.rhs) {
                    suggestions.push(format!("  Eq {}: Common subexpressions found", i + 1));
                }
            }
        }
    }

    // Analyze reaction expressions
    if let Some(ref reaction_systems) = esm_file.reaction_systems {
        for (rs_id, rs) in reaction_systems {
            println!("Reaction System: {rs_id}");

            for (i, reaction) in rs.reactions.iter().enumerate() {
                let rate_depth = expression_depth(&reaction.rate);

                if rate_depth > 4 {
                    complex_expressions += 1;
                    suggestions.push(format!(
                        "  Reaction {}: Complex rate expression (depth {})",
                        i + 1,
                        rate_depth
                    ));
                }

                if contains_redundant_operations(&reaction.rate) {
                    redundant_operations += 1;
                    suggestions.push(format!(
                        "  Reaction {}: Rate contains redundant operations",
                        i + 1
                    ));
                }
            }
        }
    }

    // Summary
    println!("\n=== OPTIMIZATION SUGGESTIONS ===");
    if suggestions.is_empty() {
        println!("✓ No obvious expression optimizations needed");
    } else {
        println!(
            "Found {complex_expressions} complex expressions, {redundant_operations} redundant operations"
        );
        for suggestion in suggestions {
            println!("{suggestion}");
        }
    }

    // General recommendations
    println!("\n=== GENERAL RECOMMENDATIONS ===");
    println!("- Factor out common subexpressions into intermediate variables");
    println!("- Use pre-computed constants where possible");
    println!("- Consider using operator precedence to reduce parentheses");
}

fn analyze_structure_optimization(esm_file: &earthsci_ast::EsmFile) {
    println!("\n=== STRUCTURE OPTIMIZATION ANALYSIS ===");

    // Component count analysis
    let model_count = esm_file.models.as_ref().map(|m| m.len()).unwrap_or(0);
    let rs_count = esm_file
        .reaction_systems
        .as_ref()
        .map(|rs| rs.len())
        .unwrap_or(0);
    let coupling_count = esm_file.coupling.as_ref().map(|c| c.len()).unwrap_or(0);

    println!(
        "Structure: {model_count} models, {rs_count} reaction systems, {coupling_count} coupling rules"
    );

    let mut suggestions: Vec<String> = Vec::new();

    // Analyze coupling complexity
    if coupling_count > 10 {
        suggestions
            .push("High coupling complexity - consider grouping related couplings".to_string());
    }

    // Analyze model sizes
    if let Some(ref models) = esm_file.models {
        for (model_id, model) in models {
            let var_count = model.variables.len();
            let eq_count = model.equations.len();

            if var_count > 50 {
                suggestions.push(format!(
                    "Model '{model_id}' has {var_count} variables - consider splitting into smaller components"
                ));
            }

            if eq_count > 30 {
                suggestions.push(format!(
                    "Model '{model_id}' has {eq_count} equations - consider modular decomposition"
                ));
            }

            // Check for unused variables
            let mut used_vars = std::collections::HashSet::new();
            for equation in &model.equations {
                collect_variables(&equation.lhs, &mut used_vars);
                collect_variables(&equation.rhs, &mut used_vars);
            }

            let unused_count = model
                .variables
                .keys()
                .filter(|&v| !used_vars.contains(v))
                .count();
            if unused_count > 0 {
                suggestions.push(format!(
                    "Model '{model_id}' has {unused_count} potentially unused variables"
                ));
            }
        }
    }

    // Analyze reaction system sizes
    if let Some(ref reaction_systems) = esm_file.reaction_systems {
        for (rs_id, rs) in reaction_systems {
            if rs.species.len() > 100 {
                suggestions.push(format!(
                    "Reaction system '{}' has {} species - consider chemical family grouping",
                    rs_id,
                    rs.species.len()
                ));
            }

            if rs.reactions.len() > 200 {
                suggestions.push(format!(
                    "Reaction system '{}' has {} reactions - consider mechanism reduction",
                    rs_id,
                    rs.reactions.len()
                ));
            }
        }
    }

    // Output recommendations
    if suggestions.is_empty() {
        println!("✓ Structure appears well-organized");
    } else {
        println!("\n=== OPTIMIZATION SUGGESTIONS ===");
        for suggestion in suggestions {
            println!("  - {suggestion}");
        }
    }

    println!("\n=== GENERAL RECOMMENDATIONS ===");
    println!("- Keep models focused on single domains/processes");
    println!("- Use hierarchical naming for related variables");
    println!("- Group related coupling operations together");
}

fn analyze_performance_optimization(esm_file: &earthsci_ast::EsmFile) {
    println!("\n=== PERFORMANCE OPTIMIZATION ANALYSIS ===");

    // Computational complexity estimation
    let mut total_ops = 0;
    let mut expensive_ops = 0;

    if let Some(ref models) = esm_file.models {
        for (model_id, model) in models {
            println!("Model: {model_id}");
            let mut model_ops = 0;

            for equation in &model.equations {
                let lhs_ops = count_operations(&equation.lhs);
                let rhs_ops = count_operations(&equation.rhs);
                model_ops += lhs_ops + rhs_ops;

                // Check for expensive operations
                if contains_expensive_operations(&equation.rhs) {
                    expensive_ops += 1;
                }
            }

            total_ops += model_ops;
            println!("  Operations per timestep: ~{model_ops}");
        }
    }

    if let Some(ref reaction_systems) = esm_file.reaction_systems {
        for (rs_id, rs) in reaction_systems {
            println!("Reaction System: {rs_id}");
            let mut rs_ops = 0;

            for reaction in &rs.reactions {
                let rate_ops = count_operations(&reaction.rate);
                rs_ops += rate_ops;

                if contains_expensive_operations(&reaction.rate) {
                    expensive_ops += 1;
                }
            }

            // Add stoichiometric operations
            rs_ops += rs.reactions.len() * rs.species.len(); // Approximate stoichiometric matrix operations

            total_ops += rs_ops;
            println!("  Operations per timestep: ~{rs_ops}");
        }
    }

    println!("\nTotal operations per timestep: ~{total_ops}");
    println!("Expensive operations found: {expensive_ops}");

    // Performance suggestions
    println!("\n=== OPTIMIZATION SUGGESTIONS ===");
    if total_ops > 10000 {
        println!("  - High computational load - consider numerical method optimization");
    }
    if expensive_ops > 5 {
        println!(
            "  - {expensive_ops} expensive operations found - consider lookup tables or approximations"
        );
    }

    // Memory usage estimation
    let var_count = esm_file
        .models
        .as_ref()
        .map(|models| models.values().map(|m| m.variables.len()).sum::<usize>())
        .unwrap_or(0);

    let species_count = esm_file
        .reaction_systems
        .as_ref()
        .map(|rs_map| rs_map.values().map(|rs| rs.species.len()).sum::<usize>())
        .unwrap_or(0);

    let total_state_vars = var_count + species_count;
    let estimated_bytes = total_state_vars * 8; // 8 bytes per double
    println!(
        "State variables: {total_state_vars} (estimated memory: {estimated_bytes} bytes / {:.1} KiB)",
        estimated_bytes as f64 / 1024.0
    );

    println!("\n=== GENERAL RECOMMENDATIONS ===");
    println!("- Use fast numerical integrators for stiff systems");
    println!("- Cache frequently computed expressions");
    println!("- Consider sparse matrix representations for large systems");
    println!("- Profile actual runtime to identify bottlenecks");
}

fn perform_structural_comparison(
    esm_file1: &earthsci_ast::EsmFile,
    esm_file2: &earthsci_ast::EsmFile,
) {
    println!("\n=== STRUCTURAL COMPARISON ===");

    let mut differences = Vec::new();

    // Compare basic structure
    if esm_file1.esm != esm_file2.esm {
        differences.push(format!(
            "ESM version: {} vs {}",
            esm_file1.esm, esm_file2.esm
        ));
    }

    // Compare component counts
    let models1_count = esm_file1.models.as_ref().map(|m| m.len()).unwrap_or(0);
    let models2_count = esm_file2.models.as_ref().map(|m| m.len()).unwrap_or(0);
    if models1_count != models2_count {
        differences.push(format!("Model count: {models1_count} vs {models2_count}"));
    }

    let rs1_count = esm_file1
        .reaction_systems
        .as_ref()
        .map(|rs| rs.len())
        .unwrap_or(0);
    let rs2_count = esm_file2
        .reaction_systems
        .as_ref()
        .map(|rs| rs.len())
        .unwrap_or(0);
    if rs1_count != rs2_count {
        differences.push(format!("Reaction system count: {rs1_count} vs {rs2_count}"));
    }

    let coupling1_count = esm_file1.coupling.as_ref().map(|c| c.len()).unwrap_or(0);
    let coupling2_count = esm_file2.coupling.as_ref().map(|c| c.len()).unwrap_or(0);
    if coupling1_count != coupling2_count {
        differences.push(format!(
            "Coupling rule count: {coupling1_count} vs {coupling2_count}"
        ));
    }

    // Compare model structures
    match (&esm_file1.models, &esm_file2.models) {
        (Some(models1), Some(models2)) => {
            let keys1: std::collections::HashSet<_> = models1.keys().collect();
            let keys2: std::collections::HashSet<_> = models2.keys().collect();

            let only_in_1: Vec<_> = keys1.difference(&keys2).collect();
            let only_in_2: Vec<_> = keys2.difference(&keys1).collect();

            if !only_in_1.is_empty() {
                differences.push(format!("Models only in file1: {only_in_1:?}"));
            }
            if !only_in_2.is_empty() {
                differences.push(format!("Models only in file2: {only_in_2:?}"));
            }

            // Compare common models
            for model_id in keys1.intersection(&keys2) {
                let model1 = &models1[*model_id];
                let model2 = &models2[*model_id];

                if model1.variables.len() != model2.variables.len() {
                    differences.push(format!(
                        "Model {}: variable count {} vs {}",
                        model_id,
                        model1.variables.len(),
                        model2.variables.len()
                    ));
                }

                if model1.equations.len() != model2.equations.len() {
                    differences.push(format!(
                        "Model {}: equation count {} vs {}",
                        model_id,
                        model1.equations.len(),
                        model2.equations.len()
                    ));
                }

                // Compare variable names
                let vars1: std::collections::HashSet<_> = model1.variables.keys().collect();
                let vars2: std::collections::HashSet<_> = model2.variables.keys().collect();
                if vars1 != vars2 {
                    differences.push(format!("Model {model_id}: different variable sets"));
                }
            }
        }
        (Some(_), None) => differences.push("File1 has models, file2 does not".to_string()),
        (None, Some(_)) => differences.push("File2 has models, file1 does not".to_string()),
        (None, None) => {}
    }

    // Compare reaction system structures
    match (&esm_file1.reaction_systems, &esm_file2.reaction_systems) {
        (Some(rs1), Some(rs2)) => {
            let keys1: std::collections::HashSet<_> = rs1.keys().collect();
            let keys2: std::collections::HashSet<_> = rs2.keys().collect();

            let only_in_1: Vec<_> = keys1.difference(&keys2).collect();
            let only_in_2: Vec<_> = keys2.difference(&keys1).collect();

            if !only_in_1.is_empty() {
                differences.push(format!("Reaction systems only in file1: {only_in_1:?}"));
            }
            if !only_in_2.is_empty() {
                differences.push(format!("Reaction systems only in file2: {only_in_2:?}"));
            }

            // Compare common reaction systems
            for rs_id in keys1.intersection(&keys2) {
                let system1 = &rs1[*rs_id];
                let system2 = &rs2[*rs_id];

                if system1.species.len() != system2.species.len() {
                    differences.push(format!(
                        "Reaction system {}: species count {} vs {}",
                        rs_id,
                        system1.species.len(),
                        system2.species.len()
                    ));
                }

                if system1.reactions.len() != system2.reactions.len() {
                    differences.push(format!(
                        "Reaction system {}: reaction count {} vs {}",
                        rs_id,
                        system1.reactions.len(),
                        system2.reactions.len()
                    ));
                }
            }
        }
        (Some(_), None) => {
            differences.push("File1 has reaction systems, file2 does not".to_string())
        }
        (None, Some(_)) => {
            differences.push("File2 has reaction systems, file1 does not".to_string())
        }
        (None, None) => {}
    }

    // Output results
    if differences.is_empty() {
        println!("✓ Files have identical structure");
    } else {
        println!("✗ Structural differences found:");
        for diff in differences {
            println!("  - {diff}");
        }
    }
}

/// Numeric difference between two `default` values, or `None` when they agree
/// to `tolerance`. Inline ARRAY data (esm-spec §6.3) compares elementwise after
/// a row-major flatten; a shape change, a scalar-vs-array change, or a ragged
/// array is itself a reported difference rather than a silent pass.
fn inline_value_difference(
    a: &earthsci_ast::InlineValue,
    b: &earthsci_ast::InlineValue,
    tolerance: f64,
) -> Option<String> {
    match (a.as_scalar(), b.as_scalar()) {
        (Some(x), Some(y)) => ((x - y).abs() > tolerance).then(|| format!("{x} vs {y}")),
        (None, None) => {
            let (sa, va) = match a.to_dense() {
                Ok(v) => v,
                Err(e) => return Some(format!("unreadable array data ({e})")),
            };
            let (sb, vb) = match b.to_dense() {
                Ok(v) => v,
                Err(e) => return Some(format!("unreadable array data ({e})")),
            };
            if sa != sb {
                return Some(format!("array shape {sa:?} vs {sb:?}"));
            }
            va.iter()
                .zip(vb.iter())
                .position(|(x, y)| (x - y).abs() > tolerance)
                .map(|i| format!("array element {i}: {} vs {}", va[i], vb[i]))
        }
        _ => Some("scalar vs array data".to_string()),
    }
}

fn perform_numerical_comparison(
    esm_file1: &earthsci_ast::EsmFile,
    esm_file2: &earthsci_ast::EsmFile,
) {
    println!("\n=== NUMERICAL COMPARISON ===");

    let tolerance = 1e-10;
    let mut differences = Vec::new();

    // Compare numerical values in models
    if let (Some(models1), Some(models2)) = (&esm_file1.models, &esm_file2.models) {
        for model_id in models1.keys() {
            if let Some(model2) = models2.get(model_id) {
                let model1 = &models1[model_id];

                // Compare variable default values
                for (var_name, var1) in &model1.variables {
                    if let Some(var2) = model2.variables.get(var_name)
                        && let (Some(default1), Some(default2)) =
                            (var1.default.as_ref(), var2.default.as_ref())
                        && let Some(detail) = inline_value_difference(default1, default2, tolerance)
                    {
                        differences.push(format!(
                            "Model {model_id}, variable {var_name}: default {detail}"
                        ));
                    }
                }

                // Compare equation expressions numerically
                if model1.equations.len() == model2.equations.len() {
                    for (i, (eq1, eq2)) in model1
                        .equations
                        .iter()
                        .zip(model2.equations.iter())
                        .enumerate()
                    {
                        if !expressions_numerically_equal(&eq1.lhs, &eq2.lhs, tolerance) {
                            differences.push(format!(
                                "Model {}, equation {} LHS: numerical difference",
                                model_id,
                                i + 1
                            ));
                        }
                        if !expressions_numerically_equal(&eq1.rhs, &eq2.rhs, tolerance) {
                            differences.push(format!(
                                "Model {}, equation {} RHS: numerical difference",
                                model_id,
                                i + 1
                            ));
                        }
                    }
                }
            }
        }
    }

    // Compare reaction system coefficients
    if let (Some(rs1), Some(rs2)) = (&esm_file1.reaction_systems, &esm_file2.reaction_systems) {
        for rs_id in rs1.keys() {
            if let Some(system2) = rs2.get(rs_id) {
                let system1 = &rs1[rs_id];

                if system1.reactions.len() == system2.reactions.len() {
                    for (i, (reaction1, reaction2)) in system1
                        .reactions
                        .iter()
                        .zip(system2.reactions.iter())
                        .enumerate()
                    {
                        // Compare substrate coefficients (substrates/products are now
                        // Option<Vec<_>>; treat None as empty when zipping).
                        let subs1: &[_] = reaction1.substrates.as_deref().unwrap_or(&[]);
                        let subs2: &[_] = reaction2.substrates.as_deref().unwrap_or(&[]);
                        if subs1.len() == subs2.len() {
                            for (sub1, sub2) in subs1.iter().zip(subs2.iter()) {
                                if sub1.coefficient != sub2.coefficient {
                                    differences.push(format!("Reaction system {}, reaction {}: substrate coefficient {} vs {}",
                                        rs_id, i + 1, sub1.coefficient, sub2.coefficient));
                                }
                            }
                        }

                        // Compare product coefficients
                        let prods1: &[_] = reaction1.products.as_deref().unwrap_or(&[]);
                        let prods2: &[_] = reaction2.products.as_deref().unwrap_or(&[]);
                        if prods1.len() == prods2.len() {
                            for (prod1, prod2) in prods1.iter().zip(prods2.iter()) {
                                if prod1.coefficient != prod2.coefficient {
                                    differences.push(format!("Reaction system {}, reaction {}: product coefficient {} vs {}",
                                        rs_id, i + 1, prod1.coefficient, prod2.coefficient));
                                }
                            }
                        }

                        // Compare rate expressions
                        if !expressions_numerically_equal(
                            &reaction1.rate,
                            &reaction2.rate,
                            tolerance,
                        ) {
                            differences.push(format!("Reaction system {}, reaction {}: rate expression numerical difference",
                                rs_id, i + 1));
                        }
                    }
                }
            }
        }
    }

    // Output results
    if differences.is_empty() {
        println!("✓ Files are numerically equivalent (tolerance: {tolerance})");
    } else {
        println!("✗ Numerical differences found (tolerance: {tolerance}):");
        for diff in differences {
            println!("  - {diff}");
        }
    }

    // Summary statistics
    let mut total_numbers1 = 0;
    let mut total_numbers2 = 0;
    count_numerical_values(esm_file1, &mut total_numbers1);
    count_numerical_values(esm_file2, &mut total_numbers2);

    println!("\nNumerical values compared: {total_numbers1} vs {total_numbers2}");
}

fn perform_deep_coupling_analysis(esm_file: &earthsci_ast::EsmFile) {
    println!("\n=== DEEP COUPLING DEPENDENCY ANALYSIS ===");

    let graph = component_graph(esm_file);
    println!(
        "Analyzing {} components with {} coupling relationships",
        graph.nodes.len(),
        graph.edges.len()
    );

    // Build dependency maps
    let mut dependencies = std::collections::HashMap::new();
    let mut dependents = std::collections::HashMap::new();

    for edge in &graph.edges {
        dependencies
            .entry(edge.to.clone())
            .or_insert_with(Vec::new)
            .push(edge.from.clone());
        dependents
            .entry(edge.from.clone())
            .or_insert_with(Vec::new)
            .push(edge.to.clone());
    }

    // Analyze coupling patterns
    analyze_coupling_patterns(&graph);

    // Find strongly connected components; a component is a genuine cycle iff
    // it has more than one node or a self-edge.
    let sccs = find_strongly_connected_components(&graph);
    let has_self_edge = |id: &str| graph.edges.iter().any(|e| e.from == id && e.to == id);
    let cycles: Vec<&Vec<String>> = sccs
        .iter()
        .filter(|c| c.len() > 1 || (c.len() == 1 && has_self_edge(&c[0])))
        .collect();
    if cycles.is_empty() {
        println!("\n✓ No circular dependencies found");
    } else {
        println!("\n=== CIRCULAR DEPENDENCIES FOUND ===");
        for (i, component) in cycles.iter().enumerate() {
            println!("Circular dependency group {}: {:?}", i + 1, component);
        }
    }

    // Analyze coupling strength
    analyze_coupling_strength(&graph, esm_file);

    // Detect coupling anti-patterns
    detect_coupling_antipatterns(&graph, esm_file);

    // Generate coupling recommendations
    println!("\n=== COUPLING RECOMMENDATIONS ===");
    if graph.edges.len() > graph.nodes.len() * 2 {
        println!(
            "- High coupling density ({} edges, {} nodes) - consider reducing interdependencies",
            graph.edges.len(),
            graph.nodes.len()
        );
    }

    // Check for hub components (components with many connections)
    let mut connection_counts = std::collections::HashMap::new();
    for edge in &graph.edges {
        *connection_counts.entry(&edge.from).or_insert(0) += 1;
        *connection_counts.entry(&edge.to).or_insert(0) += 1;
    }

    let max_connections = connection_counts.values().max().unwrap_or(&0);
    if *max_connections > 5 {
        println!(
            "- Hub components detected (max {max_connections} connections) - consider architectural refactoring"
        );
        for (component, count) in &connection_counts {
            if *count > 5 {
                println!("  High connectivity: {component} ({count} connections)");
            }
        }
    }

    // Analyze coupling by type
    let mut coupling_types = std::collections::HashMap::new();
    for edge in &graph.edges {
        *coupling_types.entry(&edge.coupling_type).or_insert(0) += 1;
    }

    println!("\n=== COUPLING TYPE ANALYSIS ===");
    for (coupling_type, count) in coupling_types {
        println!("  {coupling_type}: {count} occurrences");
    }

    // Calculate coupling metrics. The density denominator is n*(n-1), the
    // number of possible directed edges; for a graph with fewer than two
    // components it is zero, so dividing would print `NaN`. Report N/A instead
    // (and skip the qualitative verdict, which is meaningless without a pair).
    // Guarding `n < 2` also keeps the `n - 1` usize subtraction underflow-safe.
    let n_nodes = graph.nodes.len();
    if n_nodes < 2 {
        println!("\nCoupling density: N/A (fewer than two components)");
    } else {
        let coupling_density = graph.edges.len() as f64 / (n_nodes as f64 * (n_nodes - 1) as f64);
        println!("\nCoupling density: {coupling_density:.3} (0=no coupling, 1=fully coupled)");

        if coupling_density > 0.5 {
            println!("⚠ Very high coupling density - system may be difficult to maintain");
        } else if coupling_density > 0.2 {
            println!("- Moderate coupling density - monitor for growth");
        } else {
            println!("✓ Low coupling density - good architectural separation");
        }
    }
}

fn analyze_coupling_patterns(graph: &earthsci_ast::ComponentGraph) {
    println!("\n=== COUPLING PATTERNS ===");

    // Pattern 1: Fan-out (one component affects many others)
    let mut out_degree = std::collections::HashMap::new();
    let mut in_degree = std::collections::HashMap::new();

    for edge in &graph.edges {
        *out_degree.entry(&edge.from).or_insert(0) += 1;
        *in_degree.entry(&edge.to).or_insert(0) += 1;
    }

    let mut fan_out_components = Vec::new();
    let mut fan_in_components = Vec::new();

    for (component, degree) in &out_degree {
        if *degree > 3 {
            fan_out_components.push((component, degree));
        }
    }

    for (component, degree) in &in_degree {
        if *degree > 3 {
            fan_in_components.push((component, degree));
        }
    }

    if !fan_out_components.is_empty() {
        println!("Fan-out patterns (components affecting many others):");
        for (component, degree) in fan_out_components {
            println!("  {component} -> {degree} components");
        }
    }

    if !fan_in_components.is_empty() {
        println!("Fan-in patterns (components affected by many others):");
        for (component, degree) in fan_in_components {
            println!("  {component} <- {degree} components");
        }
    }

    // Pattern 2: Chains (A->B->C->D)
    let chain_length = find_longest_dependency_chain(graph);
    if chain_length > 4 {
        println!("Long dependency chain detected: {chain_length} components");
        println!("  Consider reducing chain length for better parallelization");
    }
}

fn analyze_coupling_strength(
    _graph: &earthsci_ast::ComponentGraph,
    esm_file: &earthsci_ast::EsmFile,
) {
    println!("\n=== COUPLING STRENGTH ANALYSIS ===");

    // Analyze based on coupling rules
    if let Some(ref coupling) = esm_file.coupling {
        let mut strong_couplings = 0;
        let mut weak_couplings = 0;

        for rule in coupling {
            match rule {
                earthsci_ast::CouplingEntry::Couple { systems, .. } => {
                    // Bidirectional coupling is stronger
                    strong_couplings += 1;
                    if systems.len() >= 2 {
                        println!(
                            "  Strong coupling: {} <-> {} (bidirectional)",
                            systems[0], systems[1]
                        );
                    }
                }
                earthsci_ast::CouplingEntry::OperatorCompose { systems, .. } => {
                    // Composition is medium strength
                    if systems.len() >= 2 {
                        println!(
                            "  Medium coupling: {} -> {} (composition)",
                            systems[0], systems[1]
                        );
                    }
                }
                earthsci_ast::CouplingEntry::VariableMap { from, to, .. } => {
                    // Variable mapping is weaker
                    weak_couplings += 1;
                    println!("  Weak coupling: {from} -> {to} (variable mapping)");
                }
                earthsci_ast::CouplingEntry::OperatorApply { operator, .. } => {
                    println!("  Medium coupling: {operator} (operator application)");
                }
                earthsci_ast::CouplingEntry::Callback { callback_id, .. } => {
                    println!("  Dynamic coupling: {callback_id} (callback)");
                }
                earthsci_ast::CouplingEntry::Event { name, affects, .. } => {
                    let event_name = name.as_deref().unwrap_or("unnamed");
                    let affected_count = affects.as_ref().map(|a| a.len()).unwrap_or(0);
                    println!("  Event coupling: {event_name} affects {affected_count} components");
                }
                earthsci_ast::CouplingEntry::CouplingImport { reference, .. } => {
                    println!("  Library import: {reference} (coupling_import)");
                }
            }
        }

        println!("\nCoupling strength distribution:");
        println!("  Strong (bidirectional): {strong_couplings}");
        println!(
            "  Medium (composition/operators): {}",
            coupling.len() - strong_couplings - weak_couplings
        );
        println!("  Weak (variable mapping): {weak_couplings}");

        if strong_couplings > weak_couplings {
            println!("⚠ Many strong couplings - consider loosening dependencies");
        }
    }
}

fn detect_coupling_antipatterns(
    graph: &earthsci_ast::ComponentGraph,
    _esm_file: &earthsci_ast::EsmFile,
) {
    println!("\n=== COUPLING ANTI-PATTERNS ===");

    let mut antipatterns_found = Vec::new();

    // Anti-pattern 1: God component (connected to everything)
    let mut connection_counts = std::collections::HashMap::new();
    for edge in &graph.edges {
        *connection_counts.entry(&edge.from).or_insert(0) += 1;
        *connection_counts.entry(&edge.to).or_insert(0) += 1;
    }

    for (component, count) in &connection_counts {
        if *count > graph.nodes.len() / 2 {
            antipatterns_found.push(format!(
                "God component: {} connected to {}/{} components",
                component,
                count,
                graph.nodes.len()
            ));
        }
    }

    // Anti-pattern 2: Inappropriate intimacy (too many coupling types between same components)
    let mut component_pairs = std::collections::HashMap::new();
    for edge in &graph.edges {
        let key = if edge.from < edge.to {
            format!("{}->{}", edge.from, edge.to)
        } else {
            format!("{}->{}", edge.to, edge.from)
        };
        *component_pairs.entry(key).or_insert(0) += 1;
    }

    for (pair, count) in &component_pairs {
        if *count > 2 {
            antipatterns_found.push(format!(
                "Inappropriate intimacy: {pair} ({count} coupling types)"
            ));
        }
    }

    // Anti-pattern 3: Circular dependencies (already detected in SCC analysis)

    if antipatterns_found.is_empty() {
        println!("✓ No major coupling anti-patterns detected");
    } else {
        for antipattern in antipatterns_found {
            println!("  ⚠ {antipattern}");
        }
    }
}

/// `esm init` project templates. Every template MUST load cleanly with
/// [`earthsci_ast::load_string`] — the `init_templates_are_loadable` test and a
/// pre-write self-check in the `init` command both enforce this, so the CLI
/// can never scaffold a project its own `validate` rejects. `{name}` is
/// replaced with the project name.
const TEMPLATE_MINIMAL: &str = r#"{
  "esm": "1.0.0",
  "metadata": {
    "name": "{name}",
    "description": "A minimal ESM model",
    "authors": ["Generated by earthsci-ast CLI"]
  },
  "models": {
    "{name}": {
      "variables": {},
      "equations": []
    }
  }
}"#;

const TEMPLATE_ATMOSPHERIC: &str = r#"{
  "esm": "1.0.0",
  "metadata": {
    "name": "{name}",
    "description": "An atmospheric chemistry model",
    "authors": ["Generated by earthsci-ast CLI"]
  },
  "reaction_systems": {
    "atmospheric_chemistry": {
      "species": {
        "O2": {"units": "molec/cm^3", "default": 1.0e18},
        "O3": {"units": "molec/cm^3", "default": 1.0e12}
      },
      "parameters": {
        "k": {"units": "1/s", "default": 1.0e-6, "description": "Rate constant"}
      },
      "reactions": [
        {
          "id": "o2_to_o3",
          "name": "O2_to_O3",
          "substrates": [{"species": "O2", "stoichiometry": 1}],
          "products": [{"species": "O3", "stoichiometry": 1}],
          "rate": "k"
        }
      ]
    }
  }
}"#;

const TEMPLATE_ECOSYSTEM: &str = r#"
    {
      "esm": "1.0.0",
      "metadata": {
        "name": "{name}",
        "description": "An ecosystem model",
        "authors": [
          "Generated by earthsci-ast CLI"
        ]
      },
      "models": {
        "ecosystem": {
          "variables": {
            "biomass": {
              "type": "unknown",
              "units": "kg/m^2",
              "default": 1.0
            },
            "r": {
              "type": "parameter",
              "units": "1/s",
              "default": 0.1
            }
          },
          "equations": [
            {
              "lhs": {
                "op": "D",
                "args": [
                  "biomass"
                ],
                "wrt": "t"
              },
              "rhs": {
                "op": "*",
                "args": [
                  "r",
                  "biomass"
                ]
              }
            }
          ]
        }
      }
    }
    "#;

const TEMPLATE_COUPLING: &str = r#"
    {
      "esm": "1.0.0",
      "metadata": {
        "name": "{name}",
        "description": "A coupled system model",
        "authors": [
          "Generated by earthsci-ast CLI"
        ]
      },
      "models": {
        "atmosphere": {
          "variables": {
            "temperature": {
              "type": "unknown",
              "units": "K",
              "default": 288.0
            },
            "relaxation": {
              "type": "parameter",
              "units": "1/s",
              "default": 0.01
            }
          },
          "equations": [
            {
              "lhs": {
                "op": "D",
                "args": [
                  "temperature"
                ],
                "wrt": "t"
              },
              "rhs": {
                "op": "*",
                "args": [
                  {
                    "op": "-",
                    "args": [
                      "relaxation"
                    ]
                  },
                  "temperature"
                ]
              }
            }
          ]
        },
        "ocean": {
          "variables": {
            "sea_temp": {
              "type": "unknown",
              "units": "K",
              "default": 285.0
            },
            "heat_flux": {
              "type": "parameter",
              "units": "K",
              "default": 0.0
            }
          },
          "equations": [
            {
              "lhs": {
                "op": "D",
                "args": [
                  "sea_temp"
                ],
                "wrt": "t"
              },
              "rhs": {
                "op": "*",
                "args": [
                  0.5,
                  "heat_flux"
                ]
              }
            }
          ]
        }
      },
      "coupling": [
        {
          "type": "variable_map",
          "from": "atmosphere.temperature",
          "to": "ocean.heat_flux",
          "transform": {
            "op": "*",
            "args": [
              "atmosphere.temperature",
              0.1
            ]
          },
          "description": "Ocean heat flux driven by atmospheric temperature"
        }
      ]
    }
    "#;

/// Print the semantic comparison of two loaded files (canonical-JSON
/// equality plus a coarse structural breakdown on divergence). Returns true
/// when the files are semantically identical; the caller decides whether a
/// divergence is a failing exit. Shared by `diff` and `compare --comparison-type
/// semantic`.
fn semantic_compare(
    esm_file1: &earthsci_ast::EsmFile,
    esm_file2: &earthsci_ast::EsmFile,
) -> Result<bool, Box<dyn std::error::Error>> {
    let json1 = serde_json::to_value(esm_file1)?;
    let json2 = serde_json::to_value(esm_file2)?;

    if json1 == json2 {
        println!("✓ Files are semantically identical");
        return Ok(true);
    }
    println!("✗ Files differ semantically");

    // Basic structural comparison
    if esm_file1.esm != esm_file2.esm {
        println!("  ESM version: {} vs {}", esm_file1.esm, esm_file2.esm);
    }

    let models1_count = esm_file1.models.as_ref().map(|m| m.len()).unwrap_or(0);
    let models2_count = esm_file2.models.as_ref().map(|m| m.len()).unwrap_or(0);
    if models1_count != models2_count {
        println!("  Model count: {models1_count} vs {models2_count}");
    }

    let rs1_count = esm_file1
        .reaction_systems
        .as_ref()
        .map(|rs| rs.len())
        .unwrap_or(0);
    let rs2_count = esm_file2
        .reaction_systems
        .as_ref()
        .map(|rs| rs.len())
        .unwrap_or(0);
    if rs1_count != rs2_count {
        println!("  Reaction system count: {rs1_count} vs {rs2_count}");
    }

    Ok(false)
}

/// Read a UTF-8 input file, mapping the failure to a message that names the
/// offending path (a bare `io::Error` would surface as a Debug dump with no
/// filename).
fn read_input(path: &std::path::Path) -> Result<String, Box<dyn std::error::Error>> {
    fs::read_to_string(path).map_err(|e| format!("cannot read {}: {e}", path.display()).into())
}

/// The base directory a document loaded FROM `path` must resolve against: its
/// own directory (esm-spec §4.7 for a relative `ref`, §8.2.1 for a relative
/// `data_sources` location), never the process working directory.
fn base_of(path: &std::path::Path) -> earthsci_ast::LoadOptions {
    earthsci_ast::LoadOptions {
        base_path: Some(
            path.parent()
                .filter(|p| !p.as_os_str().is_empty())
                .unwrap_or_else(|| std::path::Path::new("."))
                .to_path_buf(),
        ),
        ..Default::default()
    }
}

/// `load_string`, but anchored on the file the text came FROM.
///
/// Every subcommand here reads a document with [`read_input`] and then loaded
/// it with a bare `load_string`, which defaults the base to the process working
/// directory. That is downstream finding F7 -- recorded against `round-trip`,
/// but in fact true of ~17 subcommands, `convert` and `simulate` among them --
/// and esm-spec §8.2.1 turns it from cosmetic into dangerous. A CWD-anchored
/// relative `ref` fails loudly with "file not found"; a CWD-anchored relative
/// `url_template` can RESOLVE, succeed, and read a different file, which no
/// assertion about errors can see. Measured before this fix: the same document
/// converted from three working directories emitted three different
/// `url_template` values.
///
/// `validate` and `test` were always right because they use `load_path`. This
/// is the same anchoring for every other command that has a path in hand.
fn load_at(
    path: &std::path::Path,
    content: &str,
) -> Result<earthsci_ast::EsmFile, earthsci_ast::EsmError> {
    earthsci_ast::load_string_with_options(content, &base_of(path))
}

/// Collect every `.esm` file under `dir` (recursing when asked), in sorted
/// order per directory.
fn collect_esm_files(
    dir: &std::path::Path,
    recursive: bool,
    out: &mut Vec<PathBuf>,
) -> Result<(), Box<dyn std::error::Error>> {
    let mut entries: Vec<_> = fs::read_dir(dir)
        .map_err(|e| format!("cannot read directory {}: {e}", dir.display()))?
        .collect::<Result<_, _>>()?;
    entries.sort_by_key(|e| e.path());
    for entry in entries {
        let path = entry.path();
        if path.is_dir() {
            if recursive {
                collect_esm_files(&path, recursive, out)?;
            }
        } else if path.extension().and_then(|e| e.to_str()) == Some("esm") {
            out.push(path);
        }
    }
    Ok(())
}

/// Time `iterations` repeated parses of `content` — the loop shared by
/// `benchmark --bench-type parse` and both `performance-profile` timing
/// reports, which print the same measurement in different shapes.
fn time_parse(
    path: &std::path::Path,
    content: &str,
    iterations: usize,
) -> Result<std::time::Duration, Box<dyn std::error::Error>> {
    let start = std::time::Instant::now();
    for _ in 0..iterations {
        // Anchored on the document, like every other command: a benchmark that
        // resolved differently from `validate` would be timing a different
        // amount of work (esm-spec §4.7, §8.2.1).
        let _ = load_at(path, content)?;
    }
    Ok(start.elapsed())
}

/// Time `iterations` repeated validations of an already-loaded file.
fn time_validate(esm_file: &earthsci_ast::EsmFile, iterations: usize) -> std::time::Duration {
    let start = std::time::Instant::now();
    for _ in 0..iterations {
        let _ = validate(esm_file);
    }
    start.elapsed()
}

fn bench_parse(
    path: &std::path::Path,
    content: &str,
    iterations: usize,
) -> Result<(), Box<dyn std::error::Error>> {
    let duration = time_parse(path, content, iterations)?;
    println!(
        "Parse: {} iterations in {:?} ({:?}/iter)",
        iterations,
        duration,
        duration / (iterations as u32)
    );
    Ok(())
}

fn bench_validate(
    path: &std::path::Path,
    content: &str,
    iterations: usize,
) -> Result<(), Box<dyn std::error::Error>> {
    let esm_file = load_at(path, content)?;
    let duration = time_validate(&esm_file, iterations);
    println!(
        "Validate: {} iterations in {:?} ({:?}/iter)",
        iterations,
        duration,
        duration / (iterations as u32)
    );
    Ok(())
}

/// Time repeated 1-time-unit simulations from default initial conditions.
fn bench_simulate(
    path: &std::path::Path,
    content: &str,
    iterations: usize,
) -> Result<(), Box<dyn std::error::Error>> {
    let esm_file = load_at(path, content)?;
    let opts = earthsci_ast::SolveOptions::default();
    // Build the EsmProblem ONCE, outside the loop, and time only the solve. That
    // split is the point of the EsmProblem/`solve` surface: construction is
    // deterministic per document, `solve` is what a benchmark should measure.
    let prob = earthsci_ast::esm_problem(
        &esm_file,
        (0.0, 1.0),
        earthsci_ast::ProblemOptions {
            rhs: earthsci_ast::Rhs::Always,
            ..Default::default()
        },
    )
    .map_err(|e| format!("problem build failed: {e}"))?;
    let start = std::time::Instant::now();
    for _ in 0..iterations {
        earthsci_ast::solve(&prob, &opts).map_err(|e| format!("solve failed: {e}"))?;
    }
    let duration = start.elapsed();
    println!(
        "Simulate: {} iterations in {:?} ({:?}/iter)",
        iterations,
        duration,
        duration / (iterations as u32)
    );
    Ok(())
}

/// Expression-complexity metrics over every model equation, observed-variable
/// expression, and reaction rate: AST node count, maximum nesting depth, and
/// operator frequency (sorted by operator name for stable output).
fn print_complexity_analysis(esm_file: &earthsci_ast::EsmFile) {
    use std::collections::BTreeMap;

    fn walk(
        expr: &earthsci_ast::Expr,
        depth: usize,
        nodes: &mut usize,
        max_depth: &mut usize,
        ops: &mut BTreeMap<String, usize>,
    ) {
        *nodes += 1;
        *max_depth = (*max_depth).max(depth);
        if let earthsci_ast::Expr::Operator(node) = expr {
            *ops.entry(node.op.clone()).or_insert(0) += 1;
            node.for_each_child(&mut |child| walk(child, depth + 1, nodes, max_depth, ops));
        }
    }

    println!("\n=== COMPLEXITY ANALYSIS ===");
    let mut total_nodes = 0usize;
    let mut total_max_depth = 0usize;
    let mut ops: BTreeMap<String, usize> = BTreeMap::new();

    if let Some(ref models) = esm_file.models {
        let mut ids: Vec<&String> = models.keys().collect();
        ids.sort();
        for model_id in ids {
            let model = &models[model_id];
            let mut nodes = 0usize;
            let mut max_depth = 0usize;
            for eq in &model.equations {
                walk(&eq.lhs, 1, &mut nodes, &mut max_depth, &mut ops);
                walk(&eq.rhs, 1, &mut nodes, &mut max_depth, &mut ops);
            }
            // The Expressions a VARIABLE still carries in esm 1.0.0 are the
            // parameter-update ones (§5.4); an unknown's defining expression is
            // one of the equations already walked above.
            let mut var_names: Vec<&String> = model.variables.keys().collect();
            var_names.sort();
            for name in var_names {
                model.variables[name].for_each_expression(&mut |expr| {
                    walk(expr, 1, &mut nodes, &mut max_depth, &mut ops);
                });
            }
            println!("  Model {model_id}: {nodes} expression nodes, max depth {max_depth}");
            total_nodes += nodes;
            total_max_depth = total_max_depth.max(max_depth);
        }
    }
    if let Some(ref reaction_systems) = esm_file.reaction_systems {
        let mut ids: Vec<&String> = reaction_systems.keys().collect();
        ids.sort();
        for rs_id in ids {
            let rs = &reaction_systems[rs_id];
            let mut nodes = 0usize;
            let mut max_depth = 0usize;
            for reaction in &rs.reactions {
                walk(&reaction.rate, 1, &mut nodes, &mut max_depth, &mut ops);
            }
            println!(
                "  Reaction System {rs_id}: {nodes} rate expression nodes, max depth {max_depth}"
            );
            total_nodes += nodes;
            total_max_depth = total_max_depth.max(max_depth);
        }
    }
    println!("Total: {total_nodes} expression nodes, max depth {total_max_depth}");
    if !ops.is_empty() {
        let mut by_count: Vec<(&String, &usize)> = ops.iter().collect();
        by_count.sort_by(|a, b| b.1.cmp(a.1).then(a.0.cmp(b.0)));
        println!("Operator frequency:");
        for (op, count) in by_count.iter().take(10) {
            println!("  {op}: {count}");
        }
    }
}

// === CLI helpers ===========================================================

/// The three expression rendering formats accepted by `pretty --format`.
#[derive(Clone, Copy, ValueEnum)]
enum DisplayFormat {
    Unicode,
    Latex,
    Ascii,
}

// === Per-subcommand handlers ===============================================

fn run_validate(file: PathBuf, verbose: bool) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;

    // Perform complete validation (schema + structural). Anchor relative
    // §4.7/§9.7 refs at the input file's own directory when validating a
    // real path so a document with relative refs is not spuriously
    // rejected (stdin has no directory → process CWD).
    let base_dir = file.parent();
    let validation_result = validate_text(&content, base_dir);

    if validation_result.is_valid {
        println!("✓ Validation passed");
        if verbose && !validation_result.unit_warnings.is_empty() {
            println!("Warnings:");
            for warning in &validation_result.unit_warnings {
                // `message` is the exact string this line printed when
                // `unit_warnings` was a `Vec<String>`; the record's `path` /
                // `code` are for machine consumers.
                println!("  ⚠ {}", warning.message);
            }
        }
    } else {
        println!("✗ Validation failed");

        if !validation_result.schema_errors.is_empty() {
            println!("Schema errors:");
            for error in &validation_result.schema_errors {
                println!("  ✗ {}: {}", error.path, error.message);
            }
        }

        if !validation_result.structural_errors.is_empty() {
            println!("Structural errors:");
            for error in &validation_result.structural_errors {
                if verbose {
                    println!("  ✗ {} ({}): {}", error.path, error.code, error.message);
                    if !error.details.is_null() {
                        println!("    Details: {}", error.details);
                    }
                } else {
                    println!("  ✗ {}: {}", error.path, error.message);
                }
            }
        }
        return Err(fail_silent());
    }
    Ok(())
}

fn run_pretty(file: PathBuf, format: DisplayFormat) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    let render = |expr: &earthsci_ast::Expr| match format {
        DisplayFormat::Unicode => earthsci_ast::to_unicode(expr),
        DisplayFormat::Latex => earthsci_ast::to_latex(expr),
        DisplayFormat::Ascii => earthsci_ast::to_ascii(expr),
    };

    println!("Displaying expressions in {}:", file.display());

    // Print model equations
    if let Some(ref models) = esm_file.models {
        for (model_id, model) in models {
            println!(
                "\nModel: {} ({})",
                model_id,
                model.name.as_deref().unwrap_or("unnamed")
            );
            for (i, equation) in model.equations.iter().enumerate() {
                let lhs = render(&equation.lhs);
                let rhs = render(&equation.rhs);
                println!("  Eq {}: {} = {}", i + 1, lhs, rhs);
            }
        }
    }

    // Print reaction rates
    if let Some(ref reaction_systems) = esm_file.reaction_systems {
        for (rs_id, rs) in reaction_systems {
            println!("\nReaction System: {rs_id}");
            for (i, reaction) in rs.reactions.iter().enumerate() {
                let rate = render(&reaction.rate);
                println!("  Reaction {}: rate = {}", i + 1, rate);
            }
        }
    }
    Ok(())
}

fn run_extract(
    file: PathBuf,
    component: String,
    output: Option<PathBuf>,
) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    // Check if component exists
    if !component_exists(&esm_file, &component) {
        return Err(fail(format!(
            "Component '{component}' not found in the ESM file"
        )));
    }

    // Create a new ESM file with just the requested component. The extracted
    // component is a NEW document. Any template call sites it had were
    // already expanded at load, and metaparameters already folded, so it
    // declares neither.
    let mut extracted_esm = earthsci_ast::EsmFile {
        esm: esm_file.esm.clone(),
        metadata: esm_file.metadata.clone(),
        // Document-scoped registry (v0.8.0): preserve it on the
        // extracted component so its `{ "from": <set> }` ranges resolve.
        index_sets: esm_file.index_sets.clone(),
        domain: esm_file.domain.clone(),
        function_tables: esm_file.function_tables.clone(),
        ..Default::default()
    };

    // Extract the specific component
    if let Some(ref models) = esm_file.models
        && let Some(model) = models.get(&component)
    {
        let mut extracted_models = indexmap::IndexMap::new();
        extracted_models.insert(component.clone(), model.clone());
        extracted_esm.models = Some(extracted_models);
    }

    if let Some(ref reaction_systems) = esm_file.reaction_systems
        && let Some(rs) = reaction_systems.get(&component)
    {
        let mut extracted_rs = indexmap::IndexMap::new();
        extracted_rs.insert(component.clone(), rs.clone());
        extracted_esm.reaction_systems = Some(extracted_rs);
    }

    if let Some(ref data_sources) = esm_file.data_sources
        && let Some(dl) = data_sources.get(&component)
    {
        let mut extracted_dl = indexmap::IndexMap::new();
        extracted_dl.insert(component.clone(), dl.clone());
        extracted_esm.data_sources = Some(extracted_dl);
    }

    if let Some(ref operators) = esm_file.operators
        && let Some(op) = operators.get(&component)
    {
        let mut extracted_op = indexmap::IndexMap::new();
        extracted_op.insert(component.clone(), op.clone());
        extracted_esm.operators = Some(extracted_op);
    }

    let output_content = to_json(&extracted_esm)?;

    if let Some(output_path) = output {
        fs::write(&output_path, output_content)?;
        println!(
            "Extracted component '{}' to {}",
            component,
            output_path.display()
        );
    } else {
        print!("{output_content}");
    }
    Ok(())
}

fn run_diff(file1: PathBuf, file2: PathBuf) -> Result<(), Box<dyn std::error::Error>> {
    let content1 = read_input(&file1)?;
    let content2 = read_input(&file2)?;
    let esm_file1 = load_at(&file1, &content1)?;
    let esm_file2 = load_at(&file2, &content2)?;

    println!("Comparing {} and {}:", file1.display(), file2.display());

    if !semantic_compare(&esm_file1, &esm_file2)? {
        return Err(fail_silent());
    }
    Ok(())
}

fn run_stoich(file: PathBuf, system: String) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    if let Some(ref reaction_systems) = esm_file.reaction_systems {
        if let Some(rs) = reaction_systems.get(&system) {
            println!("Stoichiometric matrix for reaction system '{system}':");

            // The matrix is the LIBRARY's `stoichiometric_matrix`, not a second
            // computation. This used to rebuild it by hand from a SORTED
            // species order, which was harmless while the library sorted too and
            // wrong the moment it stopped: species order is declaration order in
            // every binding now (API_SPEC.md §5.10), so a hand-rolled sort here
            // would print the rows of one matrix under the column labels of
            // another.
            //
            // The library returns species-by-reaction; this view is its
            // TRANSPOSE — one line per reaction, one column per species — which
            // is a presentation choice, not a second answer.
            let matrix = stoichiometric_matrix(rs);
            let species: Vec<&String> = rs.species.keys().collect();

            print!("{:>15}", "");
            for name in &species {
                print!("{name:>10}");
            }
            println!();

            for reaction_idx in 0..rs.reactions.len() {
                print!("Reaction {:>3}:", reaction_idx + 1);
                for row in &matrix {
                    print!("{:>10.1}", row[reaction_idx]);
                }
                println!();
            }
        } else {
            return Err(fail(format!("Reaction system '{system}' not found")));
        }
    } else {
        return Err(fail("No reaction systems found in the file"));
    }
    Ok(())
}

/// Render a graph in one of the three §4.8.3 formats.
///
/// The renderers are the LIBRARY's — `earthsci_ast::to_dot` / `to_mermaid` /
/// `to_json_graph`, generic over the `Graph` trait that both `ComponentGraph`
/// and `ExpressionGraph` implement. This subcommand used to roll its own, which
/// made `esm graph` a SECOND, unpinned renderer of the same graphs, and the two
/// had drifted: the CLI wrote `[label="x", shape=y]` where the library writes
/// `[label="x" shape=y]`, it interpolated ids and labels without `escape_dot` /
/// `mermaid_label` (so a component name containing a quote produced malformed
/// DOT), and its `--format json` emitted a different document entirely — node
/// key `type` instead of `component_type`, edge keys `from`/`to`/`type` instead
/// of `source`/`target`/`data`, no `metadata`, and no `adjacency` map at all,
/// even though esm-libraries-spec §5.4.5 calls that output "JSON adjacency
/// list".
///
/// Phase 6 made the library renderers canonical and pinned them with
/// `tests/conformance/graph/cases.json`. A second path to the same renderings is
/// the defect §8 item 8 removed from Go; this removes it from the CLI.
fn render_graph<G: earthsci_ast::Graph>(graph: &G, format: GraphFormat) {
    match format {
        GraphFormat::Dot => println!("{}", earthsci_ast::to_dot(graph)),
        GraphFormat::Mermaid => println!("{}", earthsci_ast::to_mermaid(graph)),
        GraphFormat::Json => println!("{}", earthsci_ast::to_json_graph(graph)),
    }
}

fn run_graph(
    file: PathBuf,
    level: GraphLevel,
    format: GraphFormat,
    system: Option<String>,
) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    match level {
        GraphLevel::Component => render_graph(&component_graph(&esm_file), format),
        GraphLevel::Expression => {
            // Both forms esm-libraries-spec §5.4.5 documents. `--level=expression`
            // used to print "not yet implemented" and exit 1 — writing it meant
            // writing a second expression renderer, which routing through the
            // library makes unnecessary.
            match system {
                // A single named system: the ExpressionGraphInput impls for
                // Model and ReactionSystem.
                Some(name) => {
                    if let Some(model) = esm_file.models.as_ref().and_then(|m| m.get(&name)) {
                        render_graph(&expression_graph(model), format)
                    } else if let Some(rs) = esm_file
                        .reaction_systems
                        .as_ref()
                        .and_then(|r| r.get(&name))
                    {
                        render_graph(&expression_graph(rs), format)
                    } else {
                        return Err(fail(format!(
                            "Unsupported system: {name}. Name a model or reaction system \
                             declared in this file."
                        )));
                    }
                }
                // "all systems merged": the whole-document graph with
                // `variable_map` coupling folded into cross-system edges.
                None => render_graph(
                    &expression_graph_with_options(
                        &esm_file,
                        &ExpressionGraphOptions {
                            merge_coupled: true,
                        },
                    ),
                    format,
                ),
            }
        }
    }
    Ok(())
}

fn run_convert(
    input: PathBuf,
    output: Option<PathBuf>,
    to: ConvertFormat,
) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&input)?;
    let esm_file = load_at(&input, &content)?;

    let output_content = match to {
        ConvertFormat::Json => to_json(&esm_file)?,
        ConvertFormat::CompactJson => to_json_compact(&esm_file)?,
    };

    if let Some(output_path) = output {
        fs::write(&output_path, output_content)?;
        println!("Converted {} to {}", input.display(), output_path.display());
    } else {
        print!("{output_content}");
    }
    Ok(())
}

fn run_analyze(
    file: PathBuf,
    analysis_type: AnalysisType,
) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    println!("System Analysis for: {}", file.display());
    println!("Analysis Type: {}", arg_name(&analysis_type));

    if matches!(analysis_type, AnalysisType::All | AnalysisType::Structure) {
        {
            println!("\n=== STRUCTURAL ANALYSIS ===");

            // Component count analysis
            let model_count = esm_file.models.as_ref().map(|m| m.len()).unwrap_or(0);
            let rs_count = esm_file
                .reaction_systems
                .as_ref()
                .map(|rs| rs.len())
                .unwrap_or(0);
            let dl_count = esm_file
                .data_sources
                .as_ref()
                .map(|dl| dl.len())
                .unwrap_or(0);
            let op_count = esm_file.operators.as_ref().map(|op| op.len()).unwrap_or(0);

            println!(
                "Components: {model_count} models, {rs_count} reaction systems, {dl_count} data loaders, {op_count} operators"
            );

            // Equation and variable analysis
            let mut total_vars = 0;
            let mut total_eqs = 0;
            if let Some(ref models) = esm_file.models {
                for (model_id, model) in models {
                    total_vars += model.variables.len();
                    total_eqs += model.equations.len();
                    println!(
                        "  Model {}: {} variables, {} equations",
                        model_id,
                        model.variables.len(),
                        model.equations.len()
                    );
                }
            }

            // Species and reaction analysis
            let mut total_species = 0;
            let mut total_reactions = 0;
            if let Some(ref reaction_systems) = esm_file.reaction_systems {
                for (rs_id, rs) in reaction_systems {
                    total_species += rs.species.len();
                    total_reactions += rs.reactions.len();
                    println!(
                        "  Reaction System {}: {} species, {} reactions",
                        rs_id,
                        rs.species.len(),
                        rs.reactions.len()
                    );
                }
            }

            println!(
                "Total: {total_vars} variables, {total_eqs} equations, {total_species} species, {total_reactions} reactions"
            );
        }
    }
    if matches!(analysis_type, AnalysisType::All | AnalysisType::Complexity) {
        print_complexity_analysis(&esm_file);
    }
    if matches!(analysis_type, AnalysisType::All | AnalysisType::Coupling) {
        {
            println!("\n=== COUPLING ANALYSIS ===");
            if let Some(ref coupling) = esm_file.coupling {
                println!("Coupling rules: {}", coupling.len());
                for (i, rule) in coupling.iter().enumerate() {
                    match rule {
                        earthsci_ast::CouplingEntry::OperatorCompose { systems, .. } => {
                            if systems.len() >= 2 {
                                println!(
                                    "  Rule {}: {} -> {} (OperatorCompose)",
                                    i + 1,
                                    systems[0],
                                    systems[1]
                                );
                            }
                        }
                        earthsci_ast::CouplingEntry::Couple { systems, .. } => {
                            if systems.len() >= 2 {
                                println!(
                                    "  Rule {}: {} <-> {} (Couple)",
                                    i + 1,
                                    systems[0],
                                    systems[1]
                                );
                            }
                        }
                        earthsci_ast::CouplingEntry::VariableMap { from, to, .. } => {
                            println!("  Rule {}: {} -> {} (VariableMap)", i + 1, from, to);
                        }
                        earthsci_ast::CouplingEntry::OperatorApply { operator, .. } => {
                            println!("  Rule {}: {} (OperatorApply)", i + 1, operator);
                        }
                        earthsci_ast::CouplingEntry::Callback { callback_id, .. } => {
                            println!("  Rule {}: {} (Callback)", i + 1, callback_id);
                        }
                        earthsci_ast::CouplingEntry::Event { name, affects, .. } => {
                            let event_name = name.as_deref().unwrap_or("unnamed_event");
                            let systems: Vec<String> = affects
                                .as_ref()
                                .map(|affects| {
                                    affects
                                        .iter()
                                        .filter_map(|affect| {
                                            if affect.lhs.contains('.') {
                                                Some(
                                                    affect
                                                        .lhs
                                                        .split('.')
                                                        .next()
                                                        .unwrap_or("")
                                                        .to_string(),
                                                )
                                            } else {
                                                None
                                            }
                                        })
                                        .collect()
                                })
                                .unwrap_or_default();
                            println!("  Rule {}: {} -> {:?} (Event)", i + 1, event_name, systems);
                        }
                        earthsci_ast::CouplingEntry::CouplingImport { reference, .. } => {
                            println!("  Rule {}: {} (CouplingImport)", i + 1, reference);
                        }
                    }
                }
            } else {
                println!("No coupling rules defined");
            }
        }
    }
    Ok(())
}

#[allow(clippy::too_many_arguments)]
fn run_simulate(
    file: PathBuf,
    time: f64,
    output: Option<PathBuf>,
    observed: Vec<String>,
    model: Option<String>,
    format: SimulateFormat,
    compiler: Option<String>,
) -> Result<(), Box<dyn std::error::Error>> {
    // A value outside the closed vocabulary is `compiler_unknown` (esm-spec
    // §9.6.6) and is refused HERE, before anything is loaded: it is the one of
    // the three failures that says nothing about the document.
    let compiler = match compiler.as_deref() {
        None => None,
        Some(name) => Some(earthsci_ast::Compiler::parse_named(name).map_err(fail)?),
    };
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    println!("Running simulation for: {}", file.display());
    println!("Simulation time: 0 → {time}");

    // RFC decision 8: the runner must EVALUATE and append a requested observed
    // before an output plan can find a slot for it, so the request travels with
    // the solve rather than being applied to the finished solution.
    let opts = earthsci_ast::SolveOptions {
        output_observed: observed.clone(),
        ..Default::default()
    };
    // The document's own `data_sources`, bound to real readers before anything
    // is evaluated. A document that declares data and cannot read it must stop
    // here: running on with every data-fed parameter at its `default` produces
    // a trajectory that looks like an answer and is not one.
    //
    // A FACTORY, because the build below may happen twice (`ProblemOptions`
    // consumes its providers by value). `close_declared_extents` also re-loads
    // the document with every metaparameter its sources measure closed, so an
    // index set sized by one is sized by the data rather than by the
    // declaration's placeholder default.
    let (esm_file, providers) = close_declared_extents(esm_file, |mp| {
        earthsci_ast::load_string_with_options(
            &content,
            &earthsci_ast::LoadOptions {
                base_path: file.parent().map(std::path::Path::to_path_buf),
                metaparameters: mp.clone(),
            },
        )
        .map_err(|e| format!("re-loading with the discovered source extents: {e}"))
    })
    .map_err(fail)?;
    let ingesting = providers.is_some();
    let build = |pipeline: bool| -> Result<earthsci_ast::EsmProblem, Box<dyn std::error::Error>> {
        let mut popts = earthsci_ast::ProblemOptions {
            // `Auto`, not `Always`. A document with no differential equations
            // has nothing to integrate, and compiling a right-hand side for it
            // anyway is what made this command die INSIDE THE SOLVER
            // ("Exceeded maximum number of nonlinear solver failures at
            // time = 0") on a document that evaluates perfectly well. `Auto`
            // gives it the static backend and the single evaluation below.
            rhs: earthsci_ast::Rhs::Auto,
            compiler,
            model_name: model.clone(),
            build_pipeline: pipeline,
            ..Default::default()
        };
        if let Some(make) = &providers {
            popts.build_providers = make()?;
            popts.build_pipeline = true;
        }
        earthsci_ast::esm_problem(&esm_file, (0.0, time), popts)
            .map_err(|e| fail(format!("problem build failed: {e}")))
    };

    let prob = build(false)?;
    // Filled only on the static path: `--format csv` writes from the FIELDS,
    // not from flattened cell keys.
    let mut evaluated: Vec<(String, ndarray::ArrayD<f64>)> = Vec::new();
    let sol = if prob.is_dynamic() {
        // API_SPEC §5.8, "Every Problem reports what ran": one line, so a
        // number in the output below is tied to the way it was produced.
        println!("{}", prob.compiler_report());
        let sol = earthsci_ast::solve(&prob, &opts).map_err(|e| format!("solve failed: {e}"))?;
        println!(
            "✓ Simulation complete: {} output points, alg {}, retcode {}",
            sol.time.len(),
            sol.metadata.alg,
            sol.retcode
        );
        sol
    } else {
        // A RELATIONAL document: it computes rows from rows and integrates
        // nothing. Its answers are the fields the build materialized, so it is
        // EVALUATED ONCE rather than solved.
        //
        // The first build already evaluated its observed graph on the named
        // compiler — the tape under `native` — and those fields are the
        // answer. Only when it produced none (a document the array runtime
        // cannot build on its own: value invention, a derived index set) is
        // the document rebuilt with the build pipeline on, and that rebuild is
        // under the named compiler like any other build: the pipeline reports
        // each observed it evaluates, and a strict compiler refuses one it
        // would have to walk per cell. A command that wrote an empty file
        // would be the silent-empty twin of the silent zero. The pipeline has
        // already run whenever the document ingests data. The report printed
        // is the one for the build whose fields are written.
        //
        // Some fields cannot mean "some are missing": the first build's
        // state-free evaluation is all or nothing (one stateless pass over
        // EVERY observed rule of the flattened system, or no fields at all),
        // and the pipeline evaluates a subset of those rules — one model's
        // observed definitions, value-invention assignments excluded.
        let prob = if ingesting || !prob.observed_fields().is_empty() {
            prob
        } else {
            build(true)?
        };
        println!("{}", prob.compiler_report());
        evaluated = static_fields(&prob, &observed)?;
        println!(
            "✓ Static evaluation complete: {} field(s). The document declares no \
             differential equations, so it was evaluated once rather than integrated.",
            evaluated.len()
        );
        static_solution(&evaluated)
    };
    if !sol.retcode.is_success() {
        println!(
            "⚠ The run did NOT reach t = {time}: it stopped at t = {} with retcode {}",
            sol.time.last().copied().unwrap_or(0.0),
            sol.retcode
        );
    }
    if let Some(&t_final) = sol.time.last() {
        println!("Final state at t = {t_final}:");
        for (name, series) in sol.state_variable_names.iter().zip(sol.state.iter()) {
            if let Some(v) = series.last() {
                println!("  {name} = {v}");
            }
        }
    }

    if let Some(output_path) = output {
        let out = match format {
            SimulateFormat::Flat => serde_json::json!({
                "time": sol.time,
                "state": sol.state,
                "state_variable_names": sol.state_variable_names,
                "alg": sol.metadata.alg,
                "retcode": sol.retcode.name(),
            }),
            SimulateFormat::Grid => gridded_results(&esm_file, &sol, &observed)?,
            SimulateFormat::Csv => {
                if evaluated.is_empty() {
                    return Err(fail(
                        "--format csv writes ONE row set, and this document integrates: \
                         its result is a trajectory. Use --format flat or --format grid.",
                    ));
                }
                write_csv_rows(&evaluated, &output_path)?;
                println!("Results written to: {}", output_path.display());
                return Ok(());
            }
        };
        // Streamed, and pretty so a 600k-cell field is one value per line and
        // `diff` against a reference dataset reports the cells that moved.
        let mut w = std::io::BufWriter::new(fs::File::create(&output_path)?);
        serde_json::to_writer_pretty(&mut w, &out)?;
        std::io::Write::write_all(&mut w, b"\n")?;
        std::io::Write::flush(&mut w)?;
        println!("Results written to: {}", output_path.display());
    }
    Ok(())
}

// =============================================================================
// Single evaluation of a document with nothing to integrate
// =============================================================================
//
// A MOVES-shaped calculator is a RELATIONAL document: every equation computes
// rows from rows and none of them is a `D(x)/dt`. Such a document evaluates
// correctly — `esm test` asserts on it — but until now there was no route from
// its computed rows to a file: `simulate --output` was the only subcommand that
// writes values, and it compiled a right-hand side for a system with nothing to
// integrate and then failed inside the solver. `analyze` and `info` report
// structure, not values.
//
// These three functions are that route. They do not add a mode to the solver;
// they add a SINK to the evaluation the build pipeline already performs.

/// The materialized static fields a run emits: one named column per field,
/// each an N-dimensional array.
type StaticFields = Vec<(String, ndarray::ArrayD<f64>)>;

/// The build-time fields to write, in the order they will be columns.
///
/// `observed` names them explicitly (the `--observed` flag); empty means every
/// field the build materialized. A named field that the build did not produce
/// is an ERROR that says so — never an omitted column, which a comparator would
/// read as a schema change rather than as a failure.
fn static_fields(
    prob: &earthsci_ast::EsmProblem,
    observed: &[String],
) -> Result<StaticFields, Box<dyn std::error::Error>> {
    let names: Vec<String> = if observed.is_empty() {
        prob.observed_field_names()
    } else {
        observed.to_vec()
    };
    if names.is_empty() {
        return Err(fail(
            "the document declares no differential equations and its build \
             materialized no fields, so there is nothing to write. Check that the \
             values you want are OBSERVED variables (an equation whose LHS is the \
             bare name), and that any `data_sources` they read were ingested.",
        ));
    }
    let mut out = Vec::with_capacity(names.len());
    for name in names {
        let arr = earthsci_ast::observed_field(prob, &name)
            .map_err(|e| fail(format!("evaluating '{name}': {e}")))?;
        out.push((name, arr));
    }
    Ok(out)
}

/// The single-record [`earthsci_ast::Solution`] the `flat` and `grid` writers
/// render, built from `fields`.
///
/// Cell keys are `name[i,j,…]`, 1-based and COLUMN-MAJOR — the identical
/// spelling the array runtime's own slot table gives a state variable, so
/// `derive_output_plan` grids these fields exactly as it grids a trajectory's.
fn static_solution(fields: &[(String, ndarray::ArrayD<f64>)]) -> earthsci_ast::Solution {
    let mut state_variable_names = Vec::new();
    let mut state = Vec::new();
    for (name, arr) in fields {
        let shape = arr.shape();
        if shape.is_empty() {
            state_variable_names.push(name.clone());
            state.push(vec![arr[ndarray::IxDyn(&[])]]);
            continue;
        }
        for flat in 0..arr.len() {
            let mut rem = flat;
            let mut idx = vec![0usize; shape.len()];
            for (d, dim) in shape.iter().enumerate() {
                idx[d] = rem % dim;
                rem /= dim;
            }
            let label = idx
                .iter()
                .map(|i| (i + 1).to_string())
                .collect::<Vec<_>>()
                .join(",");
            state_variable_names.push(format!("{name}[{label}]"));
            state.push(vec![arr[ndarray::IxDyn(&idx)]]);
        }
    }
    earthsci_ast::Solution {
        time: vec![0.0],
        state,
        state_variable_names,
        retcode: earthsci_ast::ReturnCode::Success,
        metadata: earthsci_ast::SolutionMetadata {
            alg: "static".to_string(),
            ..Default::default()
        },
    }
}

/// Write `fields` as one row per index tuple: `i1…in`, then one column per
/// field, in natural (row-major) reading order.
///
/// Every field must share one shape, because a row is one index tuple across
/// all of them — that is what makes the table a row set rather than several
/// unrelated arrays side by side. A mismatch names the shapes rather than
/// truncating or padding, and a rank-0 field is refused for the same reason: it
/// is one value, not a column. Select a consistent set with `--observed`.
fn write_csv_rows(
    fields: &[(String, ndarray::ArrayD<f64>)],
    path: &std::path::Path,
) -> Result<(), Box<dyn std::error::Error>> {
    use std::io::Write;

    let shape = fields[0].1.shape().to_vec();
    if shape.is_empty() {
        return Err(fail(format!(
            "--format csv writes a ROW SET, and '{}' is a single value, not a column. \
             Name the row-shaped fields with --observed.",
            fields[0].0
        )));
    }
    if let Some((name, arr)) = fields.iter().find(|(_, a)| a.shape() != shape.as_slice()) {
        return Err(fail(format!(
            "--format csv writes ONE row set, so every written field must share a \
             shape: '{}' is {:?} but '{}' is {:?}. Name a consistent set with \
             --observed.",
            fields[0].0,
            shape,
            name,
            arr.shape()
        )));
    }

    let mut w = std::io::BufWriter::new(fs::File::create(path)?);
    let mut header: Vec<String> = (1..=shape.len()).map(|d| format!("i{d}")).collect();
    header.extend(fields.iter().map(|(n, _)| n.clone()));
    writeln!(w, "{}", header.join(","))?;

    let total: usize = shape.iter().product();
    let mut idx = vec![0usize; shape.len()];
    for _ in 0..total {
        let mut row: Vec<String> = idx.iter().map(|i| (i + 1).to_string()).collect();
        row.extend(
            fields
                .iter()
                .map(|(_, a)| format_csv_float(a[ndarray::IxDyn(&idx)])),
        );
        writeln!(w, "{}", row.join(","))?;
        // Row-major odometer: last axis fastest, so the file reads in the same
        // order as the nested arrays `--format grid` writes.
        for d in (0..shape.len()).rev() {
            idx[d] += 1;
            if idx[d] < shape[d] {
                break;
            }
            idx[d] = 0;
        }
    }
    w.flush()?;
    Ok(())
}

/// A float rendered for CSV without losing a bit: `{}` on an f64 is Rust's
/// shortest round-tripping form, which is what a comparator reading with
/// `float()` needs. Non-finite values are written as-is (`NaN`, `inf`) rather
/// than blanked, so a wrong answer stays visible instead of reading as missing.
fn format_csv_float(v: f64) -> String {
    format!("{v}")
}

/// The `simulate --format grid` document: [`earthsci_ast::derive_output_plan`]
/// applied to the solution's own flat slot names, then filled.
///
/// The plan is the whole of the derivation — which variables are output, how
/// the flat column-major state inverts into row-major `[time, …spatial]`
/// arrays, what the axes are really called, and which CF coordinates and
/// attributes describe them — so this function only renders it. The rendering
/// is field-for-field `earthsciio::format::OutputSchema` (see
/// `tests/data_output_zarr_roundtrip.rs`, which builds exactly this from a
/// plan), so pointing the CLI at a real Zarr sink later is a swap of the
/// serializer, not of the derivation. It is deliberately NOT done here:
/// EarthSciIO is behind the opt-in `esio` feature, which raises the crate's
/// MSRV above its declared 1.89, and `esm` ships on the default features.
fn gridded_results(
    esm_file: &earthsci_ast::EsmFile,
    sol: &earthsci_ast::Solution,
    observed: &[String],
) -> Result<serde_json::Value, Box<dyn std::error::Error>> {
    // A `--observed` name an `operator_compose` renaming match DELETED addresses
    // a variable that MOVED (issue #230), and `derive_output_plan` sees only the
    // request and the slot names — so it would report the dead spelling as
    // unknown. The solution carries the merge's map for exactly this.
    let renames = &sol.metadata.merged_variable_renames;
    let resolved: Vec<String> = observed
        .iter()
        .map(|n| {
            renames
                .get(n.as_str())
                .cloned()
                .unwrap_or_else(|| n.clone())
        })
        .collect();
    let plan = earthsci_ast::derive_output_plan(esm_file, &sol.state_variable_names, &resolved)
        .map_err(|e| fail(format!("deriving the output plan: {e}")))?;

    let n_records = sol.time.len();
    let n_slots = sol.state.len();
    let mut grids = Vec::with_capacity(plan.grids.len());
    for grid in &plan.grids {
        // One buffer per variable, `[time, …spatial]` row-major, filled a whole
        // record at a time so the flat state is transposed once.
        let mut data: Vec<Vec<f64>> = grid
            .vars
            .iter()
            .map(|v| vec![0.0f64; n_records * v.values_per_record()])
            .collect();
        let mut flat = vec![0.0f64; n_slots];
        for r in 0..n_records {
            for (i, slot) in flat.iter_mut().enumerate() {
                *slot = sol.state[i][r];
            }
            for (v, buf) in grid.vars.iter().zip(data.iter_mut()) {
                v.gridding
                    .scatter_record(&flat, buf, r)
                    .map_err(|e| fail(format!("gridding '{}': {e}", v.name)))?;
            }
        }

        let mut dims = vec![serde_json::json!([grid.time_dim, n_records])];
        dims.extend(grid.dims.iter().map(|(n, l)| serde_json::json!([n, l])));

        let mut coords = vec![serde_json::json!({
            "name": grid.time_dim,
            "values": sol.time,
            "attrs": {},
        })];
        coords.extend(
            grid.coords.iter().map(
                |c| serde_json::json!({ "name": c.name, "values": c.values, "attrs": c.attrs }),
            ),
        );

        let vars: Vec<serde_json::Value> = grid
            .vars
            .iter()
            .zip(data)
            .map(|(v, buf)| {
                serde_json::json!({
                    "name": v.name,
                    "dims": v.dims,
                    "attrs": v.attrs,
                    "dtype": v.dtype,
                    "data": buf,
                })
            })
            .collect();

        grids.push(serde_json::json!({
            "dims": dims,
            "time_dim": grid.time_dim,
            "coords": coords,
            "vars": vars,
        }));
    }

    println!(
        "Output plan: {} grid(s), {} variable(s)",
        plan.grids.len(),
        plan.n_vars()
    );
    for grid in &plan.grids {
        let axes: Vec<String> = grid.dims.iter().map(|(n, l)| format!("{n}={l}")).collect();
        println!(
            "  [{}] {}",
            if axes.is_empty() {
                "scalar".to_string()
            } else {
                axes.join(", ")
            },
            grid.vars
                .iter()
                .map(|v| v.name.as_str())
                .collect::<Vec<_>>()
                .join(", ")
        );
    }

    Ok(serde_json::json!({
        "alg": sol.metadata.alg,
        "retcode": sol.retcode.name(),
        "grids": grids,
    }))
}

fn run_info(file: PathBuf) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    println!("ESM File Information for: {}", file.display());
    println!("ESM Version: {}", esm_file.esm);

    if let Some(ref name) = esm_file.metadata.name {
        println!("Model Name: {name}");
    }
    if let Some(ref description) = esm_file.metadata.description {
        println!("Description: {description}");
    }
    if let Some(ref authors) = esm_file.metadata.authors {
        println!("Authors: {}", authors.join(", "));
    }

    if let Some(ref models) = esm_file.models {
        println!("Models: {} component(s)", models.len());
        for (id, model) in models {
            println!(
                "  - {} ({} vars, {} eqs)",
                id,
                model.variables.len(),
                model.equations.len()
            );
        }
    }

    if let Some(ref reaction_systems) = esm_file.reaction_systems {
        println!("Reaction Systems: {} component(s)", reaction_systems.len());
        for (id, rs) in reaction_systems {
            println!(
                "  - {} ({} species, {} reactions)",
                id,
                rs.species.len(),
                rs.reactions.len()
            );
        }
    }

    if let Some(ref data_sources) = esm_file.data_sources {
        println!("Data Loaders: {} component(s)", data_sources.len());
    }

    if let Some(ref operators) = esm_file.operators {
        println!("Operators: {} component(s)", operators.len());
    }

    if let Some(ref coupling) = esm_file.coupling {
        println!("Coupling Rules: {} rule(s)", coupling.len());
    }
    Ok(())
}

fn run_units(file: PathBuf, check: bool) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    println!("Units Analysis for: {}", file.display());

    if check {
        println!("\n=== DIMENSIONAL CONSISTENCY CHECK ===");
        // A provable mismatch or an unresolvable unit string is a structural
        // error, not a warning (esm-spec §4.8.4).
        let validation_result = validate(&esm_file);
        let unit_errors: Vec<_> = validation_result
            .structural_errors
            .iter()
            .filter(|e| {
                matches!(
                    e.code,
                    StructuralErrorCode::UnitInconsistency | StructuralErrorCode::UnitParseError
                )
            })
            .collect();

        let mut error_lines: Vec<String> = unit_errors
            .iter()
            .map(|e| format!("{}: {}", e.path, e.message))
            .collect();

        let (mut consistent, mut mismatched, mut not_checked) = (0, 0, 0);
        for model in esm_file.models.iter().flat_map(|m| m.values()) {
            for (_, _, verdict) in model_equation_verdicts(&esm_file, model) {
                match verdict {
                    EquationDimensionVerdict::Consistent(_) => consistent += 1,
                    EquationDimensionVerdict::Mismatch(_) => mismatched += 1,
                    EquationDimensionVerdict::NotChecked(_) => not_checked += 1,
                }
            }
        }

        // Reaction rates are judged here too. `validate` routes only the
        // STOICHIOMETRIC rate check to a structural error, so a rate whose own
        // operands are provably incompatible reaches no error list — and
        // reporting it in the reaction report while `--check` calls the document
        // consistent would be the same wrong-way-round confidence #303 is about.
        let (mut rates_ok, mut rates_bad, mut rates_unchecked) = (0, 0, 0);
        for (rs_id, rs) in esm_file.reaction_systems.iter().flat_map(|m| m.iter()) {
            for (i, verdict) in reaction_rate_verdicts(rs_id, rs, &validation_result) {
                match verdict {
                    RateVerdict::Consistent => rates_ok += 1,
                    RateVerdict::Mismatch(message) => {
                        rates_bad += 1;
                        let line = format!("/reaction_systems/{rs_id}/reactions/{i}: {message}");
                        if !error_lines.contains(&line) {
                            error_lines.push(line);
                        }
                    }
                    RateVerdict::NotChecked(_) => rates_unchecked += 1,
                }
            }
        }

        if !error_lines.is_empty() {
            println!("Unit errors:");
            for line in &error_lines {
                println!("  ✗ {line}");
            }
        }
        if !validation_result.unit_warnings.is_empty() {
            println!("Unit warnings:");
            for warning in &validation_result.unit_warnings {
                println!("  ⚠ {}", warning.message);
            }
        }

        println!(
            "Equations: {consistent} consistent, {mismatched} mismatched, {not_checked} not checked"
        );
        if rates_ok + rates_bad + rates_unchecked > 0 {
            println!(
                "Reaction rates: {rates_ok} consistent, {rates_bad} mismatched, {rates_unchecked} not checked"
            );
        }

        // The exit status follows the verdicts this command just printed, not
        // only `validate`'s error list, so the two can never disagree.
        if !error_lines.is_empty() || mismatched > 0 || rates_bad > 0 {
            return Err(fail_silent());
        }
        if not_checked == 0 && rates_unchecked == 0 {
            println!("✓ All units are dimensionally consistent");
        } else {
            println!(
                "✓ No dimensional mismatch found; run without --check to see why equations were not checked"
            );
        }
    } else {
        println!("\n=== COMPREHENSIVE UNIT ANALYSIS ===");

        // Analyze units in all components
        let validation_result = validate(&esm_file);
        analyze_model_units(&esm_file);
        analyze_reaction_system_units(&esm_file, &validation_result);

        // Unit system summary
        println!("\n=== UNIT SYSTEM SUMMARY ===");
        let unit_summary = collect_unit_types(&esm_file);

        if unit_summary.is_empty() {
            println!("No units specified in the model");
        } else {
            println!("Unit types found:");
            for unit_type in &unit_summary {
                println!("  - {unit_type}");
            }

            // Unit system recommendations
            println!("\n=== RECOMMENDATIONS ===");
            if unit_summary.contains(&"m".to_string()) && unit_summary.contains(&"cm".to_string()) {
                println!("⚠ Mixed length units detected (m, cm) - consider standardizing");
            }
            if unit_summary.contains(&"s".to_string()) && unit_summary.contains(&"min".to_string())
            {
                println!("⚠ Mixed time units detected (s, min) - consider standardizing");
            }

            // Check for SI compliance
            let si_base_units = ["m", "kg", "s", "A", "K", "mol", "cd"];
            let using_si = unit_summary
                .iter()
                .any(|u| si_base_units.contains(&u.as_str()));
            if using_si {
                println!("✓ SI base units detected - good for scientific consistency");
            } else {
                println!("ℹ Consider using SI base units for scientific consistency");
            }
        }
    }
    Ok(())
}

fn run_coupling_analysis(
    file: PathBuf,
    depth: CouplingDepth,
) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    println!("Coupling Dependency Analysis for: {}", file.display());
    println!("Depth: {}", arg_name(&depth));

    let graph = component_graph(&esm_file);
    println!(
        "Found {} components with {} coupling relationships",
        graph.nodes.len(),
        graph.edges.len()
    );

    for edge in &graph.edges {
        println!("  {} --[{}]--> {}", edge.from, edge.coupling_type, edge.to);
    }

    match depth {
        CouplingDepth::Shallow => {
            println!("Shallow analysis complete");
        }
        CouplingDepth::Deep => {
            perform_deep_coupling_analysis(&esm_file);
        }
    }
    Ok(())
}

fn run_performance_profile(
    file: PathBuf,
    profile_type: ProfileType,
) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    println!("Performance Profile for: {}", file.display());
    println!("Profile type: {}", arg_name(&profile_type));

    match profile_type {
        ProfileType::Memory => {
            println!("\n=== MEMORY ANALYSIS ===");

            // File size analysis
            let file_size = content.len();
            println!(
                "File size: {} bytes ({:.2} KB)",
                file_size,
                file_size as f64 / 1024.0
            );

            // Component memory estimation
            let _model_count = esm_file.models.as_ref().map(|m| m.len()).unwrap_or(0);
            let _rs_count = esm_file
                .reaction_systems
                .as_ref()
                .map(|rs| rs.len())
                .unwrap_or(0);
            let _dl_count = esm_file
                .data_sources
                .as_ref()
                .map(|dl| dl.len())
                .unwrap_or(0);
            let _op_count = esm_file.operators.as_ref().map(|op| op.len()).unwrap_or(0);

            // Estimate memory usage
            let mut estimated_memory = file_size; // Base JSON representation

            if let Some(ref models) = esm_file.models {
                for (model_id, model) in models {
                    let vars_mem = model.variables.len() * 64; // ~64 bytes per variable
                    let eqs_mem = model.equations.len() * 128; // ~128 bytes per equation
                    estimated_memory += vars_mem + eqs_mem;
                    println!(
                        "  Model {}: ~{} bytes ({} vars, {} eqs)",
                        model_id,
                        vars_mem + eqs_mem,
                        model.variables.len(),
                        model.equations.len()
                    );
                }
            }

            if let Some(ref reaction_systems) = esm_file.reaction_systems {
                for (rs_id, rs) in reaction_systems {
                    let species_mem = rs.species.len() * 32; // ~32 bytes per species
                    let reactions_mem = rs.reactions.len() * 256; // ~256 bytes per reaction
                    estimated_memory += species_mem + reactions_mem;
                    println!(
                        "  Reaction System {}: ~{} bytes ({} species, {} reactions)",
                        rs_id,
                        species_mem + reactions_mem,
                        rs.species.len(),
                        rs.reactions.len()
                    );
                }
            }

            println!(
                "Total estimated memory: {} bytes ({:.2} KB)",
                estimated_memory,
                estimated_memory as f64 / 1024.0
            );

            // Memory optimization suggestions
            if file_size > 10000 {
                // > 10KB
                println!("\n=== MEMORY OPTIMIZATION SUGGESTIONS ===");
                println!("- Consider using compact JSON format for smaller file sizes");
                println!("- Use references instead of duplicating common expressions");
            }
        }
        ProfileType::Time => {
            println!("\n=== TIME ANALYSIS ===");

            // Parsing benchmark
            let iterations = 100;
            let parse_duration = time_parse(&file, &content, iterations)?;
            println!(
                "Parse time: {} iterations in {:?} ({:?}/iter)",
                iterations,
                parse_duration,
                parse_duration / iterations as u32
            );

            // Validation benchmark
            let validate_duration = time_validate(&esm_file, iterations);
            println!(
                "Validate time: {} iterations in {:?} ({:?}/iter)",
                iterations,
                validate_duration,
                validate_duration / iterations as u32
            );

            // Expression complexity analysis
            if let Some(ref models) = esm_file.models {
                println!("\n=== EXPRESSION COMPLEXITY ===");
                for (model_id, model) in models {
                    let mut max_depth = 0;
                    let mut total_nodes = 0;

                    for equation in &model.equations {
                        let lhs_depth = expression_depth(&equation.lhs);
                        let rhs_depth = expression_depth(&equation.rhs);
                        max_depth = max_depth.max(lhs_depth).max(rhs_depth);

                        total_nodes += count_expression_nodes(&equation.lhs);
                        total_nodes += count_expression_nodes(&equation.rhs);
                    }

                    println!(
                        "  Model {model_id}: max depth = {max_depth}, total nodes = {total_nodes}"
                    );
                }
            }

            // Performance suggestions
            println!("\n=== PERFORMANCE SUGGESTIONS ===");
            if parse_duration.as_millis() > 100 {
                println!("- Parse time is high; consider SIMD JSON parsing for large files");
            }
            if validate_duration.as_millis() > 50 {
                println!("- Validation time is high; consider caching validation results");
            }
        }
        ProfileType::Both => {
            // Run both memory and time analysis
            println!("=== COMBINED PERFORMANCE ANALYSIS ===");
            // Reuse the above logic for both types

            // Memory analysis (simplified)
            let file_size = content.len();
            println!(
                "File size: {} bytes ({:.2} KB)",
                file_size,
                file_size as f64 / 1024.0
            );

            // Time analysis (simplified)
            let iterations = 50;
            let parse_duration = time_parse(&file, &content, iterations)?;
            println!("Parse time: {:?}/iter", parse_duration / iterations as u32);

            // Overall assessment
            println!("\n=== OVERALL ASSESSMENT ===");
            if file_size < 1000 {
                println!("✓ Small file size - good performance expected");
            } else if file_size < 10000 {
                println!("- Medium file size - acceptable performance");
            } else {
                println!("⚠ Large file size - consider optimization");
            }
        }
    }
    Ok(())
}

fn run_compare(
    file1: PathBuf,
    file2: PathBuf,
    comparison_type: ComparisonType,
) -> Result<(), Box<dyn std::error::Error>> {
    let content1 = read_input(&file1)?;
    let content2 = read_input(&file2)?;
    let esm_file1 = load_at(&file1, &content1)?;
    let esm_file2 = load_at(&file2, &content2)?;

    println!("Comparing {} and {}:", file1.display(), file2.display());
    println!("Comparison type: {}", arg_name(&comparison_type));

    match comparison_type {
        ComparisonType::Semantic => {
            // Report-only: unlike `diff`, `compare` does not set a
            // failing exit code on divergence.
            let _ = semantic_compare(&esm_file1, &esm_file2)?;
        }
        ComparisonType::Structural => {
            perform_structural_comparison(&esm_file1, &esm_file2);
        }
        ComparisonType::Numerical => {
            perform_numerical_comparison(&esm_file1, &esm_file2);
        }
    }
    Ok(())
}

fn run_optimize(file: PathBuf, opt_type: OptimizeTarget) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;
    let esm_file = load_at(&file, &content)?;

    println!("Optimization Analysis for: {}", file.display());
    println!("Optimization type: {}", arg_name(&opt_type));

    match opt_type {
        OptimizeTarget::Expression => {
            analyze_expression_optimization(&esm_file);
        }
        OptimizeTarget::Structure => {
            analyze_structure_optimization(&esm_file);
        }
        OptimizeTarget::Performance => {
            analyze_performance_optimization(&esm_file);
        }
    }
    Ok(())
}

fn run_init(name: String, template: InitTemplate) -> Result<(), Box<dyn std::error::Error>> {
    println!("Creating new ESM project: {name}");
    println!("Template: {}", arg_name(&template));

    let project_dir = PathBuf::from(&name);

    if project_dir.exists() {
        return Err(fail(format!("Directory '{name}' already exists")));
    }

    fs::create_dir(&project_dir)?;

    // Create basic ESM template (see the TEMPLATE_* consts; each is
    // load_string()-checked by the `init_templates_are_loadable` unit test).
    let template_content = match template {
        InitTemplate::Minimal => TEMPLATE_MINIMAL,
        InitTemplate::Atmospheric => TEMPLATE_ATMOSPHERIC,
        InitTemplate::Ecosystem => TEMPLATE_ECOSYSTEM,
        InitTemplate::Coupling => TEMPLATE_COUPLING,
    };

    let final_content = template_content.replace("{name}", &name);

    let esm_file = project_dir.join(format!("{name}.esm"));

    // Defensive self-check: never write a file our own loader rejects. Anchored
    // on the path we are about to WRITE, so this check resolves exactly as a
    // later `validate` of that same file will (esm-spec §4.7, §8.2.1).
    if let Err(e) = load_at(&esm_file, &final_content) {
        return Err(fail(format!(
            "internal error: generated template does not load: {e}"
        )));
    }

    fs::write(&esm_file, final_content)?;

    println!(
        "Created project '{name}' with template '{}'",
        arg_name(&template)
    );
    println!("Main file: {}", esm_file.display());
    Ok(())
}

fn run_validate_fixtures(
    directory: PathBuf,
    recursive: bool,
) -> Result<(), Box<dyn std::error::Error>> {
    println!("Validating fixtures in: {}", directory.display());
    println!("Recursive: {recursive}");

    if !directory.is_dir() {
        return Err(fail(format!(
            "Directory does not exist: {}",
            directory.display()
        )));
    }

    let mut files = Vec::new();
    collect_esm_files(&directory, recursive, &mut files)?;
    files.sort();

    if files.is_empty() {
        println!("No .esm files found");
    }

    let mut passed = 0usize;
    let mut failed = 0usize;
    for path in &files {
        let content = match fs::read_to_string(path) {
            Ok(c) => c,
            Err(e) => {
                println!("✗ {} (cannot read: {e})", path.display());
                failed += 1;
                continue;
            }
        };
        let result = validate_text(&content, path.parent());
        if result.is_valid {
            println!("✓ {}", path.display());
            passed += 1;
        } else {
            let n = result.schema_errors.len() + result.structural_errors.len();
            println!("✗ {} ({n} error(s))", path.display());
            failed += 1;
        }
    }

    println!("\n{passed} passed, {failed} failed ({} total)", files.len());
    if failed > 0 {
        return Err(fail_silent());
    }
    Ok(())
}

fn run_benchmark(
    file: PathBuf,
    bench_type: BenchType,
    iterations: usize,
) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;

    println!("Benchmarking: {}", file.display());
    println!("Type: {}, Iterations: {iterations}", arg_name(&bench_type));

    match bench_type {
        BenchType::Parse => bench_parse(&file, &content, iterations)?,
        BenchType::Validate => bench_validate(&file, &content, iterations)?,
        BenchType::Simulate => bench_simulate(&file, &content, iterations)?,
        BenchType::All => {
            bench_parse(&file, &content, iterations)?;
            bench_validate(&file, &content, iterations)?;
            // A file may be valid but not simulatable (e.g. no ODE
            // system); report that rather than failing the whole run.
            if let Err(e) = bench_simulate(&file, &content, iterations) {
                println!("Simulate: skipped ({e})");
            }
        }
    }
    Ok(())
}

fn run_schema_check(
    file: PathBuf,
    schema_version: Option<String>,
) -> Result<(), Box<dyn std::error::Error>> {
    let content = read_input(&file)?;

    println!("Schema check for: {}", file.display());
    if let Some(ref version) = schema_version {
        println!("Target schema version: {version}");
    }

    // Load performs schema validation. Anchored on the document, like every
    // other command (esm-spec §4.7, §8.2.1) -- a schema check that resolved
    // against the CWD could pass or fail on where it was run from.
    match load_at(&file, &content) {
        Ok(_) => println!("✓ Schema validation passed"),
        Err(e) => {
            println!("✗ Schema validation failed: {e}");
            return Err(fail_silent());
        }
    }
    Ok(())
}

fn run_round_trip(file: PathBuf, rounds: usize) -> Result<(), Box<dyn std::error::Error>> {
    println!(
        "Round-trip test for: {} ({} rounds)",
        file.display(),
        rounds
    );

    let original_content = read_input(&file)?;
    let mut content = original_content.clone();

    // esm-spec §4.7 / §8.2.1 both anchor on the directory of the REFERENCING
    // FILE. This command used to `load_string` with no base, so every relative
    // `ref` resolved against the process working directory instead — a document
    // that `validate` and `test` loaded fine failed to round-trip unless you
    // happened to `cd` into its directory first (downstream finding F7). With
    // §8.2.1 in place that would have been worse than cosmetic: a relative
    // `url_template` would have meant one thing under `validate` and another
    // here, which is the bug class §8.2.1 exists to close. Every round uses the
    // file's own directory, including the rounds that load EMITTED text — that
    // text stands in for the same file, so it must resolve against the same
    // base.
    let opts = earthsci_ast::LoadOptions {
        base_path: Some(
            file.parent()
                .filter(|p| !p.as_os_str().is_empty())
                .unwrap_or_else(|| std::path::Path::new("."))
                .to_path_buf(),
        ),
        ..Default::default()
    };

    for round in 1..=rounds {
        // Load
        let esm_file = earthsci_ast::load_string_with_options(&content, &opts)
            .map_err(|e| fail(format!("Round {round}: Load failed: {e}")))?;

        // Save
        content =
            to_json(&esm_file).map_err(|e| fail(format!("Round {round}: Save failed: {e}")))?;

        println!("Round {round}: OK");
    }

    // Check if final content matches original (semantically)
    let original_esm = earthsci_ast::load_string_with_options(&original_content, &opts)?;
    let final_esm = earthsci_ast::load_string_with_options(&content, &opts)?;

    let original_json = serde_json::to_value(&original_esm)?;
    let final_json = serde_json::to_value(&final_esm)?;

    if original_json == final_json {
        println!("✓ Round-trip fidelity maintained");
    } else {
        println!("✗ Round-trip fidelity lost");
        return Err(fail_silent());
    }
    Ok(())
}

// === Inline-test runner (esm-spec §6.6) ====================================

/// The SOLVER tolerances the inline-test integrations run at — deliberately
/// NOT the library defaults [`earthsci_ast::DEFAULT_RELTOL`] /
/// [`earthsci_ast::DEFAULT_ABSTOL`], and deliberately not the §6.6.4 assertion
/// tolerances either. A default is what a document gets when its author
/// expressed no opinion about accuracy; a test asserting a number to sixteen
/// digits has very much expressed one, so integrating at the library default
/// would make the test assert something about the library's default rather
/// than about the model.
///
/// `TEST_ABSTOL` matches the Python binding's `1e-14` rather than the Julia
/// runner's `1e-12`, for the reason Python records: the shared fixture corpus
/// this command is pointed at integrates concentrations of order `1e-6` while
/// asserting a RELATIVE `1e-6`, an absolute gate near `4e-13`. A solver
/// `abstol` of `1e-12` is looser than the assertion it has to clear, so the
/// run could not be accurate enough to be worth asserting on. That is a
/// property of the fixtures, not a divergence between the bindings.
const TEST_RELTOL: f64 = 1e-10;
const TEST_ABSTOL: f64 = 1e-14;

/// PASS / FAIL / ERROR — the tri-state verdict the Julia runner reports.
///
/// [`earthsci_ast::AssertionResult`] carries a two-state `passed: bool` and
/// stays that way: it is `Serialize`d verbatim by `examples/pde_conformance.rs`
/// as a payload an external runner consumes. The third state is RECOVERED from
/// it instead. `run_inline_tests` sets `actual: Some(_)` exactly on the path that
/// got as far as comparing a number against the resolved tolerance, and
/// `actual: None` on every path that failed before then — discretization
/// injection, the problem build, the solve, the solver retcode, and assertion
/// evaluation. So a non-passing result with no `actual` never reached the
/// comparison and is an ERROR; one with an `actual` reached it and lost, and
/// is a FAIL.
#[derive(Clone, Copy, PartialEq, Eq)]
enum Verdict {
    Pass,
    Fail,
    Error,
}

impl Verdict {
    /// The verdict of one finished assertion, per the `actual`-is-`None`
    /// invariant documented on [`Verdict`].
    fn of(passed: bool, actual: Option<f64>) -> Self {
        match (passed, actual) {
            (true, _) => Verdict::Pass,
            (false, None) => Verdict::Error,
            (false, Some(_)) => Verdict::Fail,
        }
    }

    fn label(self) -> &'static str {
        match self {
            Verdict::Pass => "PASS",
            Verdict::Fail => "FAIL",
            Verdict::Error => "ERROR",
        }
    }

    /// Column of the summary table's `pass | fail | err` triple.
    fn column(self) -> usize {
        match self {
            Verdict::Pass => 0,
            Verdict::Fail => 1,
            Verdict::Error => 2,
        }
    }
}

// =============================================================================
// data_sources ingest
// =============================================================================
//
// A document whose parameters declare `update: {kind: "data", …}` reads its
// numbers off disk (esm-spec §8.9). Nothing in the library opens a file on its
// own: the build-time provider contract `earthsci_ast::PrepareProvider` is
// supplied by the CALLER, and until this module existed the `esm` binary
// supplied none — so `esm test` and `esm simulate` evaluated every data-fed
// parameter at its `default` and reported the result as if it were the
// document's answer. A declared source that is never read is a WRONG ANSWER,
// not a missing feature, so every path below either ingests or says why it
// cannot, naming the source.
//
// The reader itself is EarthSciIO, behind the crate's opt-in `esio` feature
// (it pulls netcdf/hdf5 C libraries no other consumer should have to build).
// A binary built without it still detects a consumed source and refuses,
// which is the whole point: the two builds differ in what they can read, never
// in whether they tell you.

/// The `data_sources` entries some parameter actually READS, mapped to the
/// flattened names of the parameters that read them.
///
/// A source no parameter consumes is inert — a document may carry sources only
/// some of its models use — so it raises nothing. Subsystems are walked because
/// a mounted component's parameters consume the document's sources exactly as a
/// top-level model's do.
fn consumed_data_sources(doc: &serde_json::Value) -> BTreeMap<String, Vec<String>> {
    fn visit(node: &serde_json::Value, prefix: &str, out: &mut BTreeMap<String, Vec<String>>) {
        if let Some(vars) = node.get("variables").and_then(|v| v.as_object()) {
            for (vname, var) in vars {
                let Some(update) = var.get("update") else {
                    continue;
                };
                if update.get("kind").and_then(|k| k.as_str()) != Some("data") {
                    continue;
                }
                let Some(src) = update.get("source").and_then(|s| s.as_str()) else {
                    continue;
                };
                out.entry(src.to_string())
                    .or_default()
                    .push(format!("{prefix}.{vname}"));
            }
        }
        if let Some(subs) = node.get("subsystems").and_then(|v| v.as_object()) {
            for (sname, sub) in subs {
                visit(sub, &format!("{prefix}.{sname}"), out);
            }
        }
    }

    let declared: Vec<&String> = match doc.get("data_sources").and_then(|v| v.as_object()) {
        Some(m) => m.keys().collect(),
        None => return BTreeMap::new(),
    };
    let mut out = BTreeMap::new();
    if let Some(models) = doc.get("models").and_then(|v| v.as_object()) {
        for (mname, model) in models {
            visit(model, mname, &mut out);
        }
    }
    // A binding naming a source the document does not declare is a validation
    // error, not an ingest one — leave it to `validate` and ignore it here.
    out.retain(|src, _| declared.contains(&src));
    out
}

/// Where EarthSciIO's content-addressed cache lives for a CLI run.
///
/// Overridable so a batch runner can point it at fast local scratch; the
/// default is a stable per-user directory rather than a fresh temporary one, so
/// re-running `esm test` over the same sources re-reads bytes it already has.
#[cfg(feature = "esio")]
fn cache_root() -> PathBuf {
    match std::env::var_os("ESM_CACHE_DIR") {
        Some(dir) => PathBuf::from(dir),
        None => std::env::temp_dir().join("earthsci-esm-cache"),
    }
}

/// The sources being ingested and the parameters that read them, rendered for a
/// diagnostic.
///
/// Every ingest failure carries this, whatever went wrong underneath: the
/// reader's own errors name the CONSUMER (`provider Ingest.annual: …`), which
/// does not tell an author which `data_sources` entry to go and look at.
fn consumed_summary(consumed: &BTreeMap<String, Vec<String>>) -> String {
    consumed
        .iter()
        .map(|(src, readers)| format!("'{src}' (read by {})", readers.join(", ")))
        .collect::<Vec<_>>()
        .join("; ")
}

/// The message a build with no reader gives for a document it cannot ingest.
#[cfg(not(feature = "esio"))]
fn no_reader_diagnostic(consumed: &BTreeMap<String, Vec<String>>) -> String {
    format!(
        "this `esm` binary cannot ingest data_sources: it was built WITHOUT the \
         `esio` feature, so it links no data reader, but the document declares \
         {} that {} read: {}. Rebuild the binary with \
         `cargo build --release --features esio` (see pkg/earthsci-ast-rs/Cargo.toml). \
         Running anyway would evaluate every data-fed parameter at its `default` \
         and report the result as the document's answer.",
        if consumed.len() == 1 {
            "a source"
        } else {
            "sources"
        },
        if consumed.len() == 1 {
            "it"
        } else {
            "its parameters"
        },
        consumed_summary(consumed)
    )
}

/// Close every metaparameter a `data_sources` `extent` binds, by SAMPLING the
/// sources that measure themselves, and re-load the document with the sizes the
/// data gave.
///
/// esm-spec §8.9.4 puts this sampling BEFORE the document's metaparameters are
/// closed (§9.7.6 site 4), and it has to be there. An index set declared
/// `size: "N_REC"` is folded to an integer AT LOAD, so a runner that loads
/// first has already frozen the row count at the declaration's default — which,
/// for a source whose count is not knowable until the table is read, is zero.
/// Every shape downstream is then zero-length, and the two ways that surfaces
/// are both bad: the array runtime panics building a rank-1 array of length 0,
/// and an aggregate over the empty range returns a sum over no rows — 0, from a
/// document that validates, which is this toolchain's characteristic failure.
/// The alternative, a hardcoded row count in the document, is exactly what
/// `extent` exists to avoid: the next snapshot has a different number of rows.
///
/// **Best effort by construction.** Anything that goes wrong in the sampling —
/// a source that cannot be constructed, a file that is not there, a build with
/// no reader linked — is left to the RUN to report, because the run reports it
/// against the assertion or the command that needed the data, with the source
/// named. Reporting it here as well would only make it arrive twice, and the
/// first time in a worse place. Only a re-load that fails with the discovered
/// extents in hand is an error of this pass's own.
fn close_declared_extents(
    esm_file: earthsci_ast::EsmFile,
    reload: impl FnOnce(&BTreeMap<String, i64>) -> Result<earthsci_ast::EsmFile, String>,
) -> Result<
    (
        earthsci_ast::EsmFile,
        Option<Box<earthsci_ast::BuildProviderFactory<'static>>>,
    ),
    String,
> {
    let factory = |file: &earthsci_ast::EsmFile| {
        serde_json::to_value(file)
            .ok()
            .and_then(|doc| data_source_providers(&doc))
    };
    let Some(make) = factory(&esm_file) else {
        return Ok((esm_file, None));
    };
    let Ok(provs) = make() else {
        return Ok((esm_file, Some(make)));
    };

    let mut extents: BTreeMap<String, i64> = BTreeMap::new();
    for (key, mut prov) in provs {
        let Some(mp) = earthsci_ast::PrepareProvider::extent_metaparameter(&*prov) else {
            continue;
        };
        let Ok(arr) = prov.sample() else {
            return Ok((esm_file, Some(make)));
        };
        // The LEADING axis is the record count: a record-table provider serves
        // one column, rank-1 over its rows.
        let n = arr.shape().first().copied().unwrap_or(0) as i64;
        if let Some(prev) = extents.insert(mp.clone(), n)
            && prev != n
        {
            return Err(format!(
                "the data sources binding metaparameter '{mp}' disagree on its extent: \
                 one measures {prev} records and '{key}' measures {n}. Every source \
                 declaring the same `extent` metaparameter must measure the same number \
                 of records — that agreement is what keeps a table's columns aligned."
            ));
        }
    }
    if extents.is_empty() {
        return Ok((esm_file, Some(make)));
    }
    let reloaded = reload(&extents)?;
    // Rebuilt from the RELOADED document, so the providers the run ingests
    // through and the shapes it fills come from one and the same document.
    let make = factory(&reloaded);
    Ok((reloaded, make))
}

/// The build-time provider factory for `doc`'s consumed `data_sources`, or
/// `None` when the document consumes none (the overwhelmingly common case,
/// which must stay free of every cost above).
///
/// A factory rather than a set because `ProblemOptions` takes providers by
/// value and `esm test` builds one problem per inline test.
#[cfg(feature = "esio")]
fn data_source_providers(
    doc: &serde_json::Value,
) -> Option<Box<earthsci_ast::BuildProviderFactory<'static>>> {
    let consumed = consumed_data_sources(doc);
    if consumed.is_empty() {
        return None;
    }
    let doc = doc.clone();
    let cache_root = cache_root();
    let summary = consumed_summary(&consumed);
    let names: Vec<String> = consumed.into_keys().collect();
    Some(Box::new(move || {
        // NAMED, not an unrestricted sweep: `providers_from_document` silently
        // SKIPS a source with no `metadata.esio_format` when it is sweeping,
        // and a silently skipped source is the defect this whole module exists
        // to close. Naming them makes a source that cannot be constructed an
        // error that names it.
        let selected: Vec<&str> = names.iter().map(String::as_str).collect();
        earthsci_ast::esio_provider::providers_from_document(
            &doc,
            &cache_root,
            Some(&selected),
            &HashMap::new(),
        )
        .map(|provs| {
            provs
                .into_iter()
                .map(|(k, p)| (k, Box::new(p) as Box<dyn earthsci_ast::PrepareProvider>))
                .collect()
        })
        .map_err(|e| format!("data_sources ingest failed for {summary}: {}", e.0))
    }))
}

#[cfg(not(feature = "esio"))]
fn data_source_providers(
    doc: &serde_json::Value,
) -> Option<Box<earthsci_ast::BuildProviderFactory<'static>>> {
    let consumed = consumed_data_sources(doc);
    if consumed.is_empty() {
        return None;
    }
    let msg = no_reader_diagnostic(&consumed);
    Some(Box::new(move || Err(msg.clone())))
}

/// One assertion's outcome, tagged with the file it came from.
///
/// [`earthsci_ast::AssertionResult`] now carries a `file` of its own, but it is
/// empty on every row this command produces: the engine is handed an
/// already-loaded document, which has no path. This command also has to report
/// a file that fails to LOAD, and a `Verdict` the library's two-state `passed`
/// does not carry — so the runner keeps its own row type rather than the
/// library's.
struct TestRow {
    file: PathBuf,
    container: String,
    test_id: String,
    assertion_idx: usize,
    variable: String,
    time: f64,
    verdict: Verdict,
    message: String,
}

/// Render a float the way Julia's `print` does, so the two runners' summaries
/// read the same: Rust's `{}` prints `100` where Julia prints `100.0`.
fn julia_float(v: f64) -> String {
    let rendered = format!("{v}");
    if v.is_finite() && !rendered.contains(['.', 'e', 'E']) {
        format!("{rendered}.0")
    } else {
        rendered
    }
}

/// Every `.esm` file named by `paths`: a path that is a directory is walked
/// RECURSIVELY (matching the Julia runner's `discover_esm_files`, which has no
/// non-recursive mode), any other path is taken as named — including one that
/// does not exist, which then surfaces as a load ERROR naming it rather than
/// vanishing from the run. Globally sorted and deduplicated, so overlapping
/// arguments run each file exactly once.
fn discover_test_inputs(paths: &[PathBuf]) -> Result<Vec<PathBuf>, Box<dyn std::error::Error>> {
    let mut out = Vec::new();
    for path in paths {
        if path.is_dir() {
            collect_esm_files(path, true, &mut out)?;
        } else {
            out.push(path.clone());
        }
    }
    out.sort();
    out.dedup();
    Ok(out)
}

/// `path` relative to the working directory when it lies under it, else `path`
/// verbatim — the Julia summary's `rel`, which relativizes against its root
/// and otherwise prints the path as given.
fn relative_to_cwd(path: &std::path::Path) -> String {
    std::env::current_dir()
        .ok()
        .and_then(|cwd| {
            path.strip_prefix(&cwd)
                .ok()
                .map(std::path::Path::to_path_buf)
        })
        .unwrap_or_else(|| path.to_path_buf())
        .display()
        .to_string()
}

/// The top-level `models.<k>` MOUNT EDGES a document declares in its SOURCE, as
/// `"<k> ← <ref>"` lines in document order.
///
/// Read from the raw file because a LOADED document no longer shows them: the
/// mount splices the referenced leaf's model in under the same key (§4.7 /
/// §9.7.10), leaving nothing to distinguish it from a model the document wrote
/// itself. `esm test` reports them so the §6.6 rule — a mount does not carry the
/// mounted component's inline tests — is VISIBLE rather than silent: the reader
/// is told which components were not asserted on here, and where their own
/// assertions live.
///
/// Best-effort: an unreadable or unparseable file yields nothing, because the
/// run itself already reports that failure as a load ERROR row.
fn mounted_components(path: &std::path::Path) -> Vec<String> {
    let Ok(text) = fs::read_to_string(path) else {
        return Vec::new();
    };
    let Ok(raw) = serde_json::from_str::<serde_json::Value>(&text) else {
        return Vec::new();
    };
    let Some(models) = raw.get("models").and_then(|v| v.as_object()) else {
        return Vec::new();
    };
    models
        .iter()
        // The mount-edge shape `inline_toplevel_model_refs` recognises: a `ref`
        // and no inline `variables`.
        .filter(|(_, m)| m.get("ref").is_some() && m.get("variables").is_none())
        .map(|(k, m)| {
            let target = m.get("ref").and_then(|v| v.as_str()).unwrap_or("<ref>");
            format!("{k} ← {target}")
        })
        .collect()
}

#[allow(clippy::too_many_arguments)]
fn run_test(
    paths: Vec<PathBuf>,
    model: Option<String>,
    filter: Option<String>,
    solver: SolverAlg,
    reltol: Option<f64>,
    abstol: Option<f64>,
    verbose: bool,
    compiler: Option<String>,
) -> Result<(), Box<dyn std::error::Error>> {
    // Refused before anything is discovered: a value outside the closed
    // vocabulary is `compiler_unknown` and says nothing about any document.
    let compiler = match compiler.as_deref() {
        None => None,
        Some(name) => Some(earthsci_ast::Compiler::parse_named(name).map_err(fail)?),
    };
    let files = discover_test_inputs(&paths)?;
    if files.is_empty() {
        // A mis-rooted invocation must not silently pass, but nor is finding
        // nothing a test failure: warn loudly and exit 0, as the Julia runner
        // does.
        let named: Vec<String> = paths.iter().map(|p| p.display().to_string()).collect();
        println!("⚠ No .esm files discovered under: {}", named.join(", "));
        return Ok(());
    }

    // Tolerances are NOT pinned here: they resolve per FILE, below, because
    // level 2 of the esm-spec §2.2.2 chain is the document's own `solver`
    // block and this command runs many documents.
    let opts = earthsci_ast::SolveOptions {
        alg: solver.into(),
        ..Default::default()
    };

    let mut rows: Vec<TestRow> = Vec::new();
    for path in &files {
        let first_row = rows.len();
        match earthsci_ast::load_path(path) {
            // A file that will not parse is one ERROR row, not a crash that
            // abandons the remaining files.
            Err(e) => rows.push(TestRow {
                file: path.clone(),
                container: "<parse>".to_string(),
                test_id: "<load>".to_string(),
                assertion_idx: 0,
                variable: String::new(),
                time: f64::NAN,
                verdict: Verdict::Error,
                message: format!("Parse failed: {e}"),
            }),
            Ok(esm_file) => {
                // The document's own `data_sources`, bound to real readers, and
                // every metaparameter one of them measures closed to the value
                // the DATA gave before any shape is taken from it.
                //
                // The providers are derived from the LOADED file rather than
                // from the file on disk so a §4.7 subsystem `ref` has already
                // been resolved and a mounted component's data bindings are
                // visible here.
                let ready = close_declared_extents(esm_file, |mp| {
                    earthsci_ast::load_path_with_options(path, mp)
                        .map_err(|e| format!("re-loading with the discovered source extents: {e}"))
                });
                let (esm_file, providers) = match ready {
                    Ok(pair) => pair,
                    Err(message) => {
                        rows.push(TestRow {
                            file: path.clone(),
                            container: "<data_sources>".to_string(),
                            test_id: "<extent>".to_string(),
                            assertion_idx: 0,
                            variable: String::new(),
                            time: f64::NAN,
                            verdict: Verdict::Error,
                            message,
                        });
                        continue;
                    }
                };
                // `path.parent()` anchors §6.6.5 `from_file` references at the
                // document's own directory — the pinned cross-binding
                // convention, and the reason this loads with `load_path`
                // rather than `load_string` (only `load_path` resolves §4.7
                // subsystem `ref`s relative to the file).
                // `--filter` is applied INSIDE the runner, before anything is
                // built or solved, so narrowing to one test costs one test
                // rather than the whole document (it used to select rows out of
                // an already-evaluated `Vec`). The surviving rows are the same
                // either way: a result's `test_id` is its test's `id`.
                // esm-spec §2.2.2, resolved per file: an explicit `--reltol` /
                // `--abstol` wins, else this document's `solver` block, else
                // the runner defaults. The runner defaults sit at the BOTTOM of
                // the chain (they are binding defaults, not a caller's
                // opinion), which is what lets a stiff document ask for its own
                // integration accuracy without every invocation naming it.
                let file_opts = earthsci_ast::SolveOptions {
                    reltol: Some(
                        reltol
                            .or(esm_file.solver.as_ref().and_then(|s| s.reltol))
                            .unwrap_or(TEST_RELTOL),
                    ),
                    abstol: Some(
                        abstol
                            .or(esm_file.solver.as_ref().and_then(|s| s.abstol))
                            .unwrap_or(TEST_ABSTOL),
                    ),
                    ..opts.clone()
                };
                let results = earthsci_ast::run_inline_tests_with_options(
                    &esm_file,
                    &earthsci_ast::InlineTestOptions {
                        model_name: model.clone(),
                        solve: file_opts,
                        base_dir: path.parent().map(std::path::Path::to_path_buf),
                        test_filter: filter.clone(),
                        compiler,
                        ..Default::default()
                    },
                    providers.as_deref(),
                );
                for r in results {
                    rows.push(TestRow {
                        file: path.clone(),
                        container: r.model,
                        test_id: r.test_id,
                        assertion_idx: r.assertion_idx,
                        variable: r.variable,
                        time: r.time,
                        verdict: Verdict::of(r.passed, r.actual),
                        message: r.message,
                    });
                }
            }
        }
        if verbose {
            println!("{}", relative_to_cwd(path));
            for row in &rows[first_row..] {
                println!(
                    "  {:<5} {}/{}[{}] ({}@t={})",
                    row.verdict.label(),
                    row.container,
                    row.test_id,
                    row.assertion_idx,
                    row.variable,
                    julia_float(row.time),
                );
                if !row.message.is_empty() {
                    println!("        {}", row.message);
                }
            }
        }
    }

    print_test_summary(&files, &rows);
    if rows.iter().any(|r| r.verdict != Verdict::Pass) {
        return Err(fail_silent());
    }
    Ok(())
}

/// The per-file `pass | fail | err` table, its `TOTAL` row and the `Failures:`
/// block — a port of the Julia runner's `_print_summary`, so the two runners'
/// output can be read side by side.
fn print_test_summary(files: &[PathBuf], rows: &[TestRow]) {
    println!();
    println!("================ ESM inline-test summary ================");
    println!("Files discovered: {}", files.len());
    println!("Assertions:       {}", rows.len());

    // esm-spec §6.6: a mount does not carry the mounted component's inline
    // tests. Naming the mount edges keeps that VISIBLE — the reader sees which
    // components this run did not assert on, and where their assertions do run —
    // rather than losing them silently. Printed before the `rows.is_empty()`
    // exit, so a document that is nothing but mounts and coupling still says so.
    let mounts: Vec<(String, String)> = files
        .iter()
        .flat_map(|path| {
            let file = relative_to_cwd(path);
            mounted_components(path)
                .into_iter()
                .map(move |edge| (file.clone(), edge))
        })
        .collect();
    if !mounts.is_empty() {
        println!(
            "Mounted:          {} (esm-spec §6.6 — a mounted component's inline tests are not \
             run here; they run when its own file is a test target)",
            mounts.len()
        );
        for (file, edge) in &mounts {
            println!("  - {file} :: {edge}");
        }
    }

    if rows.is_empty() {
        println!("(no inline tests found)");
        println!("=========================================================");
        return;
    }

    let mut by_file: BTreeMap<String, [usize; 3]> = BTreeMap::new();
    let mut totals = [0usize; 3];
    for row in rows {
        by_file.entry(relative_to_cwd(&row.file)).or_default()[row.verdict.column()] += 1;
        totals[row.verdict.column()] += 1;
    }

    let pad = by_file
        .keys()
        .map(|name| name.chars().count())
        .max()
        .unwrap_or(20)
        .max(20);
    let rule = "-".repeat(pad + 25);
    println!(
        "  {:<pad$}  {:>5}  {:>5}  {:>5}",
        "file", "pass", "fail", "err"
    );
    println!("  {rule}");
    for (name, counts) in &by_file {
        println!(
            "  {:<pad$}  {:>5}  {:>5}  {:>5}",
            name, counts[0], counts[1], counts[2]
        );
    }
    println!("  {rule}");
    println!(
        "  {:<pad$}  {:>5}  {:>5}  {:>5}",
        "TOTAL", totals[0], totals[1], totals[2]
    );

    if totals[1] + totals[2] > 0 {
        println!();
        println!("Failures:");
        for row in rows.iter().filter(|r| r.verdict != Verdict::Pass) {
            println!(
                "  - {} :: {}/{}[{}] ({}@t={}) — {}",
                relative_to_cwd(&row.file),
                row.container,
                row.test_id,
                row.assertion_idx,
                row.variable,
                julia_float(row.time),
                row.verdict.label(),
            );
            if !row.message.is_empty() {
                println!("      {}", row.message);
            }
        }
    }
    println!("=========================================================");
}

/// The single exit point: every handler returns `Result`, and failure exits 1.
/// A [`CliFailure`] prints its own message (if any) verbatim to stderr; any
/// other error is rendered `Error: {:?}` — the same report `std`'s
/// `Termination` gives a `Result`-returning main. (Clap usage errors exit 2
/// from `Cli::parse`, as ever.)
pub fn main() -> std::process::ExitCode {
    let cli = Cli::parse();

    let result = match cli.command {
        Commands::Validate { file, verbose } => run_validate(file, verbose),
        Commands::Pretty { file, format } => run_pretty(file, format),
        Commands::Extract {
            file,
            component,
            output,
        } => run_extract(file, component, output),
        Commands::Diff { file1, file2 } => run_diff(file1, file2),
        Commands::Stoich { file, system } => run_stoich(file, system),
        Commands::Graph {
            file,
            level,
            format,
            system,
        } => run_graph(file, level, format, system),
        Commands::Convert { input, output, to } => run_convert(input, output, to),
        Commands::Analyze {
            file,
            analysis_type,
        } => run_analyze(file, analysis_type),
        Commands::Simulate {
            file,
            time,
            output,
            observed,
            model,
            format,
            compiler,
        } => run_simulate(file, time, output, observed, model, format, compiler),
        Commands::Info { file } => run_info(file),
        Commands::Units { file, check } => run_units(file, check),
        Commands::CouplingAnalysis { file, depth } => run_coupling_analysis(file, depth),
        Commands::PerformanceProfile { file, profile_type } => {
            run_performance_profile(file, profile_type)
        }
        Commands::Compare {
            file1,
            file2,
            comparison_type,
        } => run_compare(file1, file2, comparison_type),
        Commands::Optimize { file, opt_type } => run_optimize(file, opt_type),
        Commands::Init { name, template } => run_init(name, template),
        Commands::ValidateFixtures {
            directory,
            recursive,
        } => run_validate_fixtures(directory, recursive),
        Commands::Benchmark {
            file,
            bench_type,
            iterations,
        } => run_benchmark(file, bench_type, iterations),
        Commands::SchemaCheck {
            file,
            schema_version,
        } => run_schema_check(file, schema_version),
        Commands::RoundTrip { file, rounds } => run_round_trip(file, rounds),
        Commands::Test {
            paths,
            model,
            filter,
            solver,
            reltol,
            abstol,
            verbose,
            compiler,
        } => run_test(
            paths, model, filter, solver, reltol, abstol, verbose, compiler,
        ),
        Commands::ConformanceTest { out_dir, manifest } => {
            run_conformance_test(&out_dir, manifest.as_deref())
        }
    };

    match result {
        Ok(()) => std::process::ExitCode::SUCCESS,
        Err(e) => {
            match e.downcast::<CliFailure>() {
                Ok(failure) => {
                    if let Some(message) = failure.message {
                        eprintln!("{message}");
                    }
                }
                Err(other) => eprintln!("Error: {other:?}"),
            }
            std::process::ExitCode::FAILURE
        }
    }
}

// === Cross-language conformance PRODUCER ====================================
//
// Reads the shared CORPUS MANIFEST (`scripts/conformance_corpus.py`) and emits a
// record for every entry in it. The producer does NOT enumerate the corpus
// itself: each producer used to walk `tests/valid` / `tests/invalid`
// NON-recursively and all four skipped the same 69 fixtures — the entire
// `faq` and `template_imports` corpora — plus `lib/**`, which nothing swept
// at all (audit 2026-07-14, F5; CONFORMANCE_SPEC §2.2.1).
//
// Every validation entry runs the full **load → resolve → validate** pipeline.
// `load_path` resolves §4.7 subsystem refs against the file's own directory; the
// producer used to call `validate_text(&content)` on the file's TEXT, which
// takes no base path, so no `{ref}` could ever resolve here.
//
// See `scripts/run-python-conformance.py` for the emitted wire shape; every
// producer emits the same one.

#[derive(serde::Deserialize)]
struct ConformanceManifest {
    validation_files: Vec<ManifestValidationEntry>,
    display_cases: Vec<ManifestDisplayCase>,
    substitution_cases: Vec<ManifestSubstitutionCase>,
}

#[derive(serde::Deserialize)]
struct ManifestValidationEntry {
    id: String,
    path: String,
}

#[derive(serde::Deserialize)]
struct ManifestDisplayCase {
    id: String,
    // The manifest's `kind` ("formula" | "expression") is deliberately not read
    // here: `Expr` deserializes a bare string as `Expr::Variable`, and the
    // renderers chemical-subscript a variable name exactly as they do a formula,
    // so both kinds take the same path. Bindings whose renderer dispatches on a
    // static type (Julia, Go) do need the tag. Serde ignores the extra key.
    input: serde_json::Value,
}

#[derive(serde::Deserialize)]
struct ManifestSubstitutionCase {
    id: String,
    input: serde_json::Value,
    bindings: std::collections::BTreeMap<String, serde_json::Value>,
}

/// Walk up from the current directory to the repo root, so the manifest's
/// repo-relative paths resolve no matter where the binary is invoked from.
fn find_repo_root() -> Option<std::path::PathBuf> {
    let mut dir = std::env::current_dir().ok()?;
    loop {
        if dir.join("tests").join("valid").is_dir() {
            return Some(dir);
        }
        if !dir.pop() {
            return None;
        }
    }
}

/// The manifest path is passed by the harness; there is no fallback sweep. A
/// producer that invents its own corpus when the manifest is missing is a
/// producer that can silently under-report coverage. Fail instead.
fn locate_manifest(
    out_dir: &std::path::Path,
    explicit: Option<&std::path::Path>,
) -> Result<std::path::PathBuf, Box<dyn std::error::Error>> {
    let candidate = match explicit {
        Some(p) => p.to_path_buf(),
        None => match std::env::var("ESM_CONFORMANCE_MANIFEST") {
            Ok(v) if !v.is_empty() => std::path::PathBuf::from(v),
            _ => out_dir
                .parent()
                .unwrap_or(std::path::Path::new("."))
                .join("corpus_manifest.json"),
        },
    };
    if !candidate.is_file() {
        return Err(format!(
            "Corpus manifest not found: {}\nGenerate it with: python3 scripts/conformance_corpus.py --output <path>",
            candidate.display()
        )
        .into());
    }
    Ok(candidate)
}

fn schema_error_json(e: &earthsci_ast::SchemaError) -> serde_json::Value {
    serde_json::json!({
        "path": e.path,
        "message": e.message,
        "keyword": e.keyword,
        "code": e.keyword,
        "details": {},
    })
}

fn structural_error_json(e: &earthsci_ast::StructuralError) -> serde_json::Value {
    // `code` is an enum whose Display is the snake_case spec spelling
    // (`undefined_variable`, …) — the same spelling `expected_errors.json` pins.
    let code = e.code.to_string();
    serde_json::json!({
        "path": e.path,
        "message": e.message,
        "code": code,
        "keyword": code,
        "details": e.details,
    })
}

/// Recover subsystem-`ref` resolution errors (esm-spec §4.7) as structured
/// `(code, path)` findings when the load aborted before the typed validator
/// could run.
///
/// The §4.7 ref mechanism resolves against the referencing FILE'S OWN directory,
/// and a bad ref fails inside `load` (`ref_loading`) as an opaque string error —
/// so the producer, which owns the file path, reconstructs the finding here.
/// Scans every `models.<M>.subsystems.<N>` and
/// `reaction_systems.<R>.subsystems.<N>` edge for a `{ "ref": <str> }` and
/// classifies it: a ref whose file cannot be found (or is remote) is
/// `unresolved_subsystem_ref`; a ref resolving to a file that does NOT hold
/// exactly one top-level system is `ambiguous_subsystem_ref` (matching
/// `ref_loading::extract_single_system`). A ref that resolves cleanly
/// contributes nothing — the load failed for some other reason. The pointer is
/// the offending subsystem slot, `/<container>/<system>/subsystems/<name>`.
fn collect_subsystem_ref_errors(
    content: &str,
    file_path: &std::path::Path,
) -> Vec<earthsci_ast::StructuralError> {
    use earthsci_ast::{StructuralError, StructuralErrorCode};
    let mut out = Vec::new();
    let Ok(json) = serde_json::from_str::<serde_json::Value>(content) else {
        return out;
    };
    let base = file_path
        .parent()
        .unwrap_or_else(|| std::path::Path::new("."));

    // A top-level `reaction_systems.<k>` `{ref}` is a mount form this binding
    // does not implement: `mount_form_unsupported` at the entry itself
    // (`ref_loading::refuse_toplevel_reaction_system_refs`).
    if let Some(rs) = json.get("reaction_systems").and_then(|v| v.as_object()) {
        for (name, entry) in rs {
            if let Some(ref_str) = entry.get("ref").and_then(|v| v.as_str())
                && entry.get("species").is_none()
            {
                out.push(StructuralError {
                    path: format!("/reaction_systems/{name}"),
                    code: StructuralErrorCode::MountFormUnsupported,
                    message: format!(
                        "a top-level `reaction_systems.<k>` `{{ref}}` mount (ref '{ref_str}') is \
                         not supported by this binding"
                    ),
                    details: serde_json::json!({ "ref": ref_str }),
                });
            }
        }
    }

    for container in ["models", "reaction_systems"] {
        let Some(systems) = json.get(container).and_then(|v| v.as_object()) else {
            continue;
        };
        for (sys_name, sys_val) in systems {
            let Some(subs) = sys_val.get("subsystems").and_then(|v| v.as_object()) else {
                continue;
            };
            for (sub_name, sub_val) in subs {
                let Some(ref_str) = sub_val.get("ref").and_then(|v| v.as_str()) else {
                    continue;
                };
                let path = format!("/{container}/{sys_name}/subsystems/{sub_name}");
                let details = serde_json::json!({
                    "ref": ref_str,
                    "subsystem": sub_name,
                    "parent_model": sys_name,
                });

                // Remote refs are not resolvable by the local loader.
                if ref_str.starts_with("http://") || ref_str.starts_with("https://") {
                    out.push(StructuralError {
                        path,
                        code: StructuralErrorCode::UnresolvedSubsystemRef,
                        message: format!(
                            "Subsystem reference '{ref_str}' is a remote URL and could not be resolved"
                        ),
                        details,
                    });
                    continue;
                }

                let Ok(canonical) = base.join(ref_str).canonicalize() else {
                    out.push(StructuralError {
                        path,
                        code: StructuralErrorCode::UnresolvedSubsystemRef,
                        message: format!(
                            "Subsystem reference '{ref_str}' could not be resolved — file does not exist"
                        ),
                        details,
                    });
                    continue;
                };

                let doc = fs::read_to_string(&canonical)
                    .ok()
                    .and_then(|s| serde_json::from_str::<serde_json::Value>(&s).ok());
                let Some(doc) = doc else {
                    out.push(StructuralError {
                        path,
                        code: StructuralErrorCode::UnresolvedSubsystemRef,
                        message: format!(
                            "Subsystem reference '{ref_str}' could not be read or parsed"
                        ),
                        details,
                    });
                    continue;
                };

                // "Exactly one top-level system" counts models, reaction systems
                // and data loaders TOGETHER (mirrors `extract_single_system`).
                let count = ["models", "reaction_systems", "data_sources"]
                    .iter()
                    .filter_map(|k| doc.get(*k).and_then(|v| v.as_object()))
                    .map(|m| m.len())
                    .sum::<usize>();
                if count != 1 {
                    out.push(StructuralError {
                        path,
                        code: StructuralErrorCode::AmbiguousSubsystemRef,
                        message: format!(
                            "Subsystem reference '{ref_str}' resolves to a file containing \
                             multiple top-level systems; exactly one is required"
                        ),
                        details,
                    });
                }
            }
        }
    }
    out
}

fn run_conformance_test(
    out_dir: &std::path::Path,
    manifest_path: Option<&std::path::Path>,
) -> Result<(), Box<dyn std::error::Error>> {
    use serde_json::json;

    fs::create_dir_all(out_dir)?;

    let manifest_file = locate_manifest(out_dir, manifest_path)?;
    let manifest: ConformanceManifest = serde_json::from_str(&fs::read_to_string(&manifest_file)?)?;
    let repo_root = find_repo_root()
        .ok_or("Could not locate repo root by walking up from the current directory")?;

    let validation = conformance_validation_sweep(&manifest, &repo_root);
    let display = conformance_display_results(&manifest);
    let substitution = conformance_substitution_results(&manifest)?;

    let results = json!({
        "language": "rust",
        "validation_results": validation,
        "display_results": display,
        "substitution_results": substitution,
        "errors": Vec::<String>::new(),
    });

    let results_file = out_dir.join("results.json");
    fs::write(&results_file, serde_json::to_string_pretty(&results)?)?;
    println!(
        "Rust conformance results written to: {}",
        results_file.display()
    );
    println!(
        "✓ Validation sweep completed ({} files); display {} cases; substitution {} cases",
        manifest.validation_files.len(),
        manifest.display_cases.len(),
        manifest.substitution_cases.len()
    );

    Ok(())
}

/// Validation phase: load -> resolve -> validate, one record per manifest file.
fn conformance_validation_sweep(
    manifest: &ConformanceManifest,
    repo_root: &std::path::Path,
) -> std::collections::BTreeMap<String, serde_json::Value> {
    let mut validation: std::collections::BTreeMap<String, serde_json::Value> =
        std::collections::BTreeMap::new();
    for entry in &manifest.validation_files {
        let path = repo_root.join(&entry.path);
        validation.insert(
            entry.id.clone(),
            serde_json::Value::Object(conformance_validation_record(&path)),
        );
    }
    validation
}

/// One validation file's conformance record: schema/structural findings, the
/// answering phase, and the accept/reject verdict.
fn conformance_validation_record(
    path: &std::path::Path,
) -> serde_json::Map<String, serde_json::Value> {
    use serde_json::json;

    let mut record = serde_json::Map::new();
    record.insert("schema_errors".into(), json!([]));
    record.insert("structural_errors".into(), json!([]));

    let content = fs::read_to_string(path).unwrap_or_default();

    // load_path anchors §4.7 subsystem refs and §9.7 template imports at the
    // file's own directory: this IS the resolve phase.
    let loaded = earthsci_ast::load_path(path);
    let resolve_ok = loaded.is_ok();
    record.insert("resolve_ok".into(), json!(resolve_ok));

    match &loaded {
        Ok(esm_file) => {
            // Schema validation judges the document AS WRITTEN, so it runs
            // on the raw text; structural validation judges the RESOLVED
            // form, so it runs on the loaded file (which has its §4.7 refs
            // spliced in). `validate_text` is the only exported entry
            // that reports schema errors, so its structural half is
            // discarded — it saw the unresolved document.
            let schema = validate_text(&content, path.parent());
            let structural = validate(esm_file);
            record.insert(
                "schema_errors".into(),
                json!(
                    schema
                        .schema_errors
                        .iter()
                        .map(schema_error_json)
                        .collect::<Vec<_>>()
                ),
            );
            record.insert(
                "structural_errors".into(),
                json!(
                    structural
                        .structural_errors
                        .iter()
                        .map(structural_error_json)
                        .collect::<Vec<_>>()
                ),
            );
            record.insert(
                "is_valid".into(),
                json!(schema.schema_errors.is_empty() && structural.structural_errors.is_empty()),
            );
            record.insert("phase".into(), json!("validate"));
        }
        Err(e) => {
            record.insert("error".into(), json!(e.to_string()));
            record.insert("error_type".into(), json!("LoadError"));
            record.insert("phase".into(), json!("load"));
            record.insert("is_valid".into(), json!(false));
            // Raw document: the load-phase rejection still yields whatever
            // structured findings the binding is able to enumerate.
            // `validate_text` now itself recovers the typed structural
            // `(code, path)` records on a load rejection (best-effort raw parse
            // + typed `validate()`), so a document rejected at load no longer
            // records `is_valid:false` with an EMPTY `structural_errors`
            // (CONFORMANCE_SPEC §7.1.2). We must therefore NOT re-run
            // `validate()` here, or every such record would be duplicated.
            let raw = validate_text(&content, path.parent());
            let mut structural_errors: Vec<earthsci_ast::StructuralError> =
                raw.structural_errors.clone();
            // A subsystem `ref` that could not be resolved (missing file) or is
            // ambiguous (resolves to a file with != 1 top-level system) aborts
            // the load before any typed pass AND before `validate_text` can
            // deserialize the document, so those findings are recovered here
            // directly from the raw document against the file's own directory.
            structural_errors.extend(collect_subsystem_ref_errors(&content, path));

            record.insert(
                "schema_errors".into(),
                json!(
                    raw.schema_errors
                        .iter()
                        .map(schema_error_json)
                        .collect::<Vec<_>>()
                ),
            );
            record.insert(
                "structural_errors".into(),
                json!(
                    structural_errors
                        .iter()
                        .map(structural_error_json)
                        .collect::<Vec<_>>()
                ),
            );
        }
    }

    // The verdict is "did this binding accept the document", regardless of
    // WHICH phase answered. A rejection at resolve is still a rejection.
    let is_valid = record
        .get("is_valid")
        .and_then(|v| v.as_bool())
        .unwrap_or(false)
        && resolve_ok;
    record.insert("is_valid".into(), json!(is_valid));
    record.insert(
        "outcome".into(),
        json!(if is_valid { "valid" } else { "invalid" }),
    );
    record
}

/// Display phase: render every manifest case in all three formats.
///
/// The producer used to emit a literal `"display_results": {}` here, which the
/// comparator scored as a category with zero divergences — a permanent,
/// structurally unfailable 1.00 that diluted the one real category (audit C2).
fn conformance_display_results(
    manifest: &ConformanceManifest,
) -> std::collections::BTreeMap<String, serde_json::Value> {
    use serde_json::{Value, json};

    let mut display: std::collections::BTreeMap<String, Value> = std::collections::BTreeMap::new();
    for case in &manifest.display_cases {
        let parsed: Result<earthsci_ast::Expr, _> = serde_json::from_value(case.input.clone());
        match parsed {
            Ok(expr) => {
                display.insert(
                    case.id.clone(),
                    json!({
                        "unicode": earthsci_ast::to_unicode(&expr),
                        "latex": earthsci_ast::to_latex(&expr),
                        "ascii": earthsci_ast::to_ascii(&expr),
                    }),
                );
            }
            Err(e) => {
                display.insert(
                    case.id.clone(),
                    json!({
                        "unicode": Value::Null,
                        "latex": Value::Null,
                        "ascii": Value::Null,
                        "errors": { "parse": e.to_string() },
                    }),
                );
            }
        }
    }
    display
}

/// Substitution phase.
fn conformance_substitution_results(
    manifest: &ConformanceManifest,
) -> Result<std::collections::BTreeMap<String, serde_json::Value>, Box<dyn std::error::Error>> {
    use serde_json::{Value, json};

    let mut substitution: std::collections::BTreeMap<String, Value> =
        std::collections::BTreeMap::new();
    for case in &manifest.substitution_cases {
        let expr: Result<earthsci_ast::Expr, _> = serde_json::from_value(case.input.clone());
        let mut bindings: HashMap<String, earthsci_ast::Expr> = HashMap::new();
        let mut binding_err: Option<String> = None;
        for (name, value) in &case.bindings {
            match serde_json::from_value::<earthsci_ast::Expr>(value.clone()) {
                Ok(e) => {
                    bindings.insert(name.clone(), e);
                }
                Err(e) => binding_err = Some(e.to_string()),
            }
        }
        match (expr, binding_err) {
            (Ok(e), None) => {
                let result = earthsci_ast::substitute(&e, &bindings);
                substitution.insert(
                    case.id.clone(),
                    json!({ "result": serde_json::to_value(&result)? }),
                );
            }
            (Err(e), _) => {
                substitution.insert(
                    case.id.clone(),
                    json!({ "result": Value::Null, "error": e.to_string() }),
                );
            }
            (_, Some(e)) => {
                substitution.insert(
                    case.id.clone(),
                    json!({ "result": Value::Null, "error": e }),
                );
            }
        }
    }
    Ok(substitution)
}

#[cfg(test)]
mod tests {
    use super::*;

    /// Every `esm init` template must load cleanly — the CLI must never
    /// scaffold a project its own `validate` rejects.
    #[test]
    fn init_templates_are_loadable() {
        for (name, template) in [
            ("minimal", TEMPLATE_MINIMAL),
            ("atmospheric", TEMPLATE_ATMOSPHERIC),
            ("ecosystem", TEMPLATE_ECOSYSTEM),
            ("coupling", TEMPLATE_COUPLING),
        ] {
            let content = template.replace("{name}", "demo_project");
            let esm_file = earthsci_ast::load_string(&content)
                .unwrap_or_else(|e| panic!("template '{name}' does not load: {e}"));
            let result = validate(&esm_file);
            assert!(
                result.is_valid,
                "template '{name}' fails validation: schema={:?} structural={:?}",
                result.schema_errors, result.structural_errors
            );
        }
    }
}