#!/usr/bin/env python3
"""Audit frozen shape/field artifacts without reevaluating FEM.

Use the NumPy-enabled Python environment. The SVD is taken directly on the tall
matrix of saved block actions: squaring its condition number through a Gram
matrix would obscure precisely the weak response this experiment investigates.
Physical observations use broken scalar gradients, not RT transport or DG jump
distributions. Neither these observations nor a small restricted singular value
identify a physical gauge or establish the rank of the complete Jacobian.
"""
import argparse
import csv
import hashlib
import json
import math
import os
from pathlib import Path
import sys

# Avoid competing with a later numerical campaign through a threaded BLAS pool.
for _name in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ[_name] = "1"
import numpy as np

P = "weak_shape_"
PP = P + "physical_"
EPS = np.finfo(float).eps
RTOL = 5e-10
FIELDS = ("density", "specific_enthalpy", "gravity_potential")
INTEGRALS = ("base_integral", "pullback_integral", "transport_integral", "eulerian_integral",
             "base_squared_integral", "pullback_squared_integral", "transport_squared_integral",
             "eulerian_squared_integral", "pullback_transport_cross_integral")
SQUARES = INTEGRALS[4:8]


def digest(path):
    h = hashlib.sha256()
    with Path(path).open("rb") as stream:
        for chunk in iter(lambda: stream.read(1024 * 1024), b""):
            h.update(chunk)
    return h.hexdigest()


def rows(path):
    with Path(path).open(newline="") as stream:
        reader = csv.DictReader(stream)
        if not reader.fieldnames or len(set(reader.fieldnames)) != len(reader.fieldnames):
            raise ValueError(f"Missing or duplicate CSV columns: {path}")
        result = list(reader)
    if any(None in row or any(v is None for v in row.values()) for row in result):
        raise ValueError(f"Malformed CSV row: {path}")
    return result


def finite(value):
    number = float(value)
    if not math.isfinite(number):
        raise ValueError(f"Nonfinite available value: {value}")
    return number


def integer(value):
    number = finite(value)
    if number != int(number):
        raise ValueError(f"Nonintegral count/index: {value}")
    return int(number)


def vector(path):
    with Path(path).open() as stream:
        size = int(stream.readline().strip())
        data = np.loadtxt(stream, dtype=float, ndmin=1)
    if size <= 0 or data.ndim != 1 or data.size != size or not np.isfinite(data).all():
        raise ValueError(f"Invalid finite vector: {path}")
    return data


def norm(value):
    return float(np.linalg.norm(value))


def checked_einsum(subscripts, *operands):
    # NumPy 2.2.6 / Accelerate can flag invalid/divide/overflow for a bounded
    # 257x6 matrix product whose answer agrees with scalar fsum to 4e-15.
    # Direct contractions avoid that BLAS path; finite and reconstruction
    # checks remain mandatory, with no warning suppression or relaxed bounds.
    if any(not np.isfinite(operand).all() for operand in operands):
        raise ValueError("Nonfinite contraction input")
    result = np.einsum(subscripts, *operands, optimize=False)
    if not np.isfinite(result).all():
        raise ValueError("Nonfinite contraction result")
    return result


def near(a, b, scale=0.):
    return math.isfinite(a) and math.isfinite(b) and abs(a - b) <= (
        RTOL * max(abs(a), abs(b)) + 2048 * EPS * abs(scale) + 1e-28)


def vector_near(a, b, scale=0.):
    return a.shape == b.shape and near(norm(a - b), 0., scale + norm(a) + norm(b))


def ratio(a, b):
    return finite(a / b) if b > 0 else None


class Audit:
    def __init__(self):
        self.count = 0
        self.failures = []

    def check(self, name, condition):
        self.count += 1
        if not bool(condition):
            self.failures.append(name)

    def close(self, name, actual, expected, scale=0.):
        self.check(name, near(actual, expected, scale))


def table(path, key):
    result = {}
    for row in rows(path):
        index = key(row)
        if index in result:
            raise ValueError(f"Duplicate row {index}: {path}")
        result[index] = row
    return result


def layout(path, size):
    result = {}
    all_rows = rows(path)
    if any(r["kind"] not in ("value", "residual") for r in all_rows):
        raise ValueError("Unknown layout kind")
    for kind in ("value", "residual"):
        entries = [r for r in all_rows if r["kind"] == kind]
        offset, names, blocks = 0, set(), []
        for row in entries:
            start, count = integer(row["offset"]), integer(row["size"])
            if start != offset or count <= 0 or row["block"] in names:
                raise ValueError(f"Invalid {kind} layout partition")
            names.add(row["block"])
            blocks.append((row["block"], start, count))
            offset += count
        if offset != size:
            raise ValueError(f"Incomplete {kind} layout partition")
        result[kind] = blocks
    return result


def restricted_svd(direction, actions, blocks, audit):
    norms = np.array([norm(direction[o:o+n]) for _, o, n in blocks])
    active = [i for i, x in enumerate(norms) if x > 0.]
    for i, x in enumerate(norms):
        if x == 0.:
            audit.check(f"zero_direction_block_action:{i}", norm(actions[i]) == 0.)
    if not active:
        raise ValueError("No nonzero block directions")
    matrix = np.column_stack([actions[i] / norms[i] for i in active])
    if not np.isfinite(matrix).all():
        raise ValueError("Nonfinite scaled block action")
    u, singular, vh = np.linalg.svd(matrix, full_matrices=False)
    weights = vh[-1].copy()
    # Fix sign for reproducible reporting; repeated singular values still permit
    # different equally valid vectors, which the JSON explicitly acknowledges.
    if weights[np.argmax(np.abs(weights))] < 0:
        weights *= -1
    combination = np.zeros_like(direction)
    action = np.zeros_like(direction)
    coefficients = {}
    for weight, i in zip(weights, active):
        name, start, count = blocks[i]
        coefficient = float(weight / norms[i])
        combination[start:start+count] = coefficient * direction[start:start+count]
        action += coefficient * actions[i]
        coefficients[name] = {"orthonormal_basis_weight": float(weight),
                              "saved_block_direction_multiplier": coefficient}
    scale = float(singular[0])
    audit.close("restricted_svd:unit_direction", norm(combination), 1.)
    audit.close("restricted_svd:smallest_action_norm", norm(action), float(singular[-1]), scale)
    audit.check("restricted_svd:action_reconstruction", vector_near(action, checked_einsum("ij,j->i",matrix,weights), scale))
    audit.check("restricted_svd:factorization", vector_near(matrix, checked_einsum("ij,jk->ik",u*singular,vh), norm(matrix)))
    audit.check("restricted_svd:finite", np.isfinite(singular).all() and np.isfinite(vh).all())
    return {
        "scope": "Jacobian restricted to nonzero normalized saved block directions; not full Jacobian rank or nullspace",
        "basis": "disjoint block directions divided by their normalized coefficient two-norm",
        "active_blocks": [blocks[i][0] for i in active], "singular_values": singular.tolist(),
        "condition_ratio_in_restricted_span": ratio(float(singular[0]), float(singular[-1])),
        "smallest_direction_norm": norm(combination), "smallest_action_norm": norm(action),
        "coefficients": coefficients,
        "reconstructed_direction_float64_sha256": hashlib.sha256(combination.astype("<f8").tobytes()).hexdigest(),
        "reconstructed_action_float64_sha256": hashlib.sha256(action.astype("<f8").tobytes()).hexdigest(),
        "smallest_singular_value_gap": float(singular[-2] - singular[-1]) if len(singular) > 1 else None,
        "qualification": "basis metric is coefficient normalization; vector sign fixed, degenerate subspace basis not unique"
    }, action


def surface_fit(surface, physical, physical_direction, blocks, audit):
    shape = next(b for b in blocks if b[0] == "surface_deformation")
    _, offset, size = shape
    audit.check("surface:complete_parameter_indices", sorted(surface) == list(range(size)))
    ordered = [surface[i] for i in range(size)]
    audit.check("surface:unique_volume_dofs", len({integer(r["volume_scalar_true_dof"]) for r in ordered}) == size)
    audit.check("surface:nonnegative_volume_dofs", all(integer(r["volume_scalar_true_dof"]) >= 0 for r in ordered))
    xyz = np.array([[finite(r["unit_radial_" + axis]) for axis in "xyz"] for r in ordered])
    data = np.array([finite(r["physical_direction_parameter"]) for r in ordered])
    for i, row in enumerate(ordered):
        r0, r = finite(row["reference_radius"]), finite(row["physical_radius"])
        b = finite(row["base_parameter"])
        audit.check(f"surface:radii:{i}", r0 > 0 and r > 0)
        audit.close(f"surface:unit_radial:{i}", norm(xyz[i]), 1.)
        audit.close(f"surface:radius_identity:{i}", r, r0 + b, abs(r0) + abs(b))
        audit.close(f"surface:base:{i}", b, float(physical[offset+i]))
        audit.close(f"surface:direction:{i}", data[i], float(physical_direction[offset+i]))
    x, y, z = xyz.T
    columns = [np.ones(size), x, y, z, x*y, x*z, y*z, x*x-y*y, 3*z*z-1]
    fits = []
    for degree, count in ((0, 1), (1, 4), (2, 9)):
        matrix = np.column_stack(columns[:count])
        coef, _, rank, singular = np.linalg.lstsq(matrix, data, rcond=None)
        residual = data - checked_einsum("ij,j->i",matrix,coef)
        audit.check(f"surface:fit_finite:{degree}", np.isfinite(coef).all() and np.isfinite(singular).all())
        fits.append({"maximum_degree": degree, "rank": int(rank), "basis_coefficients": coef.tolist(),
                     "relative_nodal_l2_remainder": ratio(norm(residual), norm(data)),
                     "nodal_l2_remainder": norm(residual), "singular_values": singular.tolist()})
    return {"scope": "unweighted nodal least squares in real l=0,1,2 polynomial spans; not area-weighted harmonic coefficients or gauge identification",
            "basis": ["1", "nx", "ny", "nz", "nx*ny", "nx*nz", "ny*nz", "nx^2-ny^2", "3*nz^2-1"],
            "direction_nodal_l2": norm(data), "direction_maximum_absolute": float(np.max(np.abs(data))), "fits": fits}


def physical_analysis(folder, physical, physical_direction, blocks, audit, required):
    def read(suffix, key):
        name = PP + suffix + ".csv"
        required.add(name)
        return table(folder / name, key)
    metadata = read("metadata", lambda r: r["key"])
    meta = {k: r["value"] for k, r in metadata.items()}
    audit.check("physical:schema", meta["schema"] == "weak_shape_physical_v1")
    audit.check("physical:declared_scalar_scope", "broken" in meta["dg_derivatives"] and "not evaluated" in meta["gravity_vector"])
    elem = read("elements", lambda r: (integer(r["quadrature"]), integer(r["element"]), r["field"]))
    geometry = read("geometry", lambda r: (integer(r["quadrature"]), integer(r["element"])))
    affine = read("affine", lambda r: (integer(r["quadrature"]), integer(r["element"])))
    surface = read("surface_nodes", lambda r: integer(r["parameter"]))
    summary = read("summary", lambda r: (integer(r["quadrature"]), r["field"], r["metric"]))
    observed = {key: finite(r["value"]) for key, r in summary.items()}
    expected_summary = set()
    def compare(q, field, metric, expected, scale=0.):
        key = (q, field, metric)
        expected_summary.add(key)
        audit.close("physical:summary:" + str(key), observed[key], expected, scale)
    audit.check("physical:quadrature_coverage", {q for q, _ in geometry} == {10, 14})
    ids = {e for q, e in geometry if q == 10}
    audit.check("physical:same_elements_at_both_rules", ids and ids == {e for q, e in geometry if q == 14})
    audit.check("physical:affine_coverage", set(affine) == set(geometry))
    audit.check("physical:field_coverage", set(elem) == {(q, e, f) for q, e in geometry for f in FIELDS})
    outputs = []
    for q in (10, 14):
        geo = [geometry[(q, e)] for e in sorted(ids)]
        aff = [affine[(q, e)] for e in sorted(ids)]
        for e, g, a in zip(sorted(ids), geo, aff):
            for key in g:
                if key not in ("quadrature", "element"):
                    finite(g[key])
            volume = finite(g["physical_volume"])
            audit.check(f"physical:geometry_positive:{q}:{e}", volume > 0 and integer(g["points"]) > 0)
            for prefix in ("reference", "physical"):
                lo, hi = finite(g[prefix+"_radius_min"]), finite(g[prefix+"_radius_max"])
                audit.check(f"physical:radius_range:{q}:{e}:{prefix}", 0 <= lo <= hi)
            for key in ("base_displacement_squared_integral", "direction_displacement_squared_integral"):
                audit.check(f"physical:geometry_square:{q}:{e}:{key}", finite(g[key]) >= 0)
            for key, family in (("density_element_order", "density_family_order"),
                                ("enthalpy_element_order", "enthalpy_family_order"),
                                ("potential_element_order", "potential_family_order"),
                                ("displacement_element_order", "displacement_family_order")):
                audit.check(f"physical:order:{q}:{e}:{key}", integer(g[key]) == integer(meta[family]))
            for key in ("reference_volume", "value_error_squared_integral", "reference_value_squared_integral",
                        "gradient_error_squared_integral", "reference_gradient_squared_integral"):
                audit.check(f"physical:affine_nonnegative:{q}:{e}:{key}", finite(a[key]) >= 0)
            audit.check(f"physical:affine_volume:{q}:{e}", finite(a["reference_volume"]) > 0)
            audit.check(f"physical:affine_points:{q}:{e}", integer(a["points"]) == integer(g["points"]))
            audit.close(f"physical:affine_identity:{q}:{e}", finite(a["reference_gradient_squared_integral"]), 3*finite(a["reference_volume"]))
        totals_g = {k: math.fsum(finite(g[k]) for g in geo) for k in
                    ("physical_volume", "base_displacement_squared_integral", "direction_displacement_squared_integral")}
        for key, value in totals_g.items():
            compare(q, "geometry", key, value)
        compare(q, "geometry", "stellar_elements", len(ids))
        totals_a = {k: math.fsum(finite(a[k]) for a in aff) for k in
                    ("reference_volume", "value_error_squared_integral", "reference_value_squared_integral",
                     "gradient_error_squared_integral", "reference_gradient_squared_integral")}
        for key, value in totals_a.items():
            compare(q, "affine", key, value)
        for label, num, den in (("relative_value_l2_error", "value_error_squared_integral", "reference_value_squared_integral"),
                                ("relative_gradient_l2_error", "gradient_error_squared_integral", "reference_gradient_squared_integral")):
            compare(q, "affine", label, math.sqrt(totals_a[num] / totals_a[den]))
        field_outputs = {}
        for field in FIELDS:
            records = [elem[(q, e, field)] for e in sorted(ids)]
            for e, row, g in zip(sorted(ids), records, geo):
                values = {k: finite(row[k]) for k in INTEGRALS}
                audit.check(f"physical:gradient_kind:{q}:{e}:{field}", row["gradient_kind"] == ("continuous_H1" if field == "specific_enthalpy" else "broken_DG"))
                for key in ("attribute", "geometry", "points"):
                    audit.check(f"physical:element_binding:{q}:{e}:{field}:{key}", integer(row[key]) == integer(g[key]))
                audit.close(f"physical:volume_binding:{q}:{e}:{field}", finite(row["physical_volume"]), finite(g["physical_volume"]))
                audit.check(f"physical:squares:{q}:{e}:{field}", all(values[k] >= 0 for k in SQUARES))
                for prefix in ("base","pullback","transport","eulerian"):
                    first,second=values[prefix+"_integral"],values[prefix+"_squared_integral"]
                    product=finite(row["physical_volume"])*second
                    audit.check(f"physical:moment_cauchy:{q}:{e}:{field}:{prefix}", first*first <= product+2048*EPS*(first*first+product)+1e-28)
                # Signed element moments may themselves cancel. Cauchy bounds
                # their absolute quadrature sums using the saved positive moments.
                scale = math.sqrt(finite(row["physical_volume"])) * (math.sqrt(values["pullback_squared_integral"])+math.sqrt(values["transport_squared_integral"]))
                audit.close(f"physical:signed_identity:{q}:{e}:{field}", values["eulerian_integral"], values["pullback_integral"]-values["transport_integral"], scale)
                d2, t2, e2, cross = (values[k] for k in ("pullback_squared_integral", "transport_squared_integral", "eulerian_squared_integral", "pullback_transport_cross_integral"))
                audit.close(f"physical:squared_identity:{q}:{e}:{field}", e2, d2+t2-2*cross, d2+t2+2*abs(cross))
                audit.check(f"physical:cauchy:{q}:{e}:{field}", abs(cross) <= math.sqrt(d2)*math.sqrt(t2) + 2048*EPS*(d2+t2) + 1e-28)
            totals = {k: math.fsum(finite(r[k]) for r in records) for k in INTEGRALS}
            for key in (*SQUARES, "pullback_transport_cross_integral"):
                compare(q, field, key, totals[key])
            for prefix in ("pullback", "transport", "eulerian"):
                compare(q, field, prefix+"_l2_norm", math.sqrt(totals[prefix+"_squared_integral"]))
            d2, t2, e2, cross = (totals[k] for k in ("pullback_squared_integral", "transport_squared_integral", "eulerian_squared_integral", "pullback_transport_cross_integral"))
            compare(q, field, "squared_norm_identity_error", e2-d2-t2+2*cross, d2+t2+2*abs(cross))
            audit.close(f"physical:total_squared_identity:{q}:{field}", e2, d2+t2-2*cross, d2+t2+2*abs(cross))
            # Each element spans a radius interval. Midpoint bins are explicitly
            # element-group summaries, not integrals over exact spherical shells.
            bins = []
            outer = max(finite(g["physical_radius_max"]) for g in geo)
            for low, high in ((0., .25), (.25, .5), (.5, .75), (.75, .9), (.9, 1.0000001)):
                chosen = [i for i, g in enumerate(geo) if low <= .5*(finite(g["physical_radius_min"])+finite(g["physical_radius_max"]))/outer < high]
                bins.append({"normalized_element_midradius_interval": [low, min(high, 1.)], "elements": len(chosen),
                             "eulerian_squared_integral": math.fsum(finite(records[i]["eulerian_squared_integral"]) for i in chosen),
                             "pullback_squared_integral": math.fsum(finite(records[i]["pullback_squared_integral"]) for i in chosen)})
            audit.check(f"physical:radial_bin_coverage:{q}:{field}", sum(b["elements"] for b in bins) == len(ids))
            top = sorted(zip(sorted(ids), records), key=lambda pair: finite(pair[1]["eulerian_squared_integral"]), reverse=True)[:10]
            field_outputs[field] = {"integrals": totals, "eulerian_over_pullback_l2": ratio(math.sqrt(e2), math.sqrt(d2)),
                                    "eulerian_over_pullback_plus_transport_l2": ratio(math.sqrt(e2), math.sqrt(d2)+math.sqrt(t2)),
                                    "pullback_l2_per_unit_displacement_l2": ratio(math.sqrt(d2), math.sqrt(totals_g["direction_displacement_squared_integral"])),
                                    "transport_l2_per_unit_displacement_l2": ratio(math.sqrt(t2), math.sqrt(totals_g["direction_displacement_squared_integral"])),
                                    "eulerian_l2_per_unit_displacement_l2": ratio(math.sqrt(e2), math.sqrt(totals_g["direction_displacement_squared_integral"])),
                                    "eulerian_rms": math.sqrt(e2/totals_g["physical_volume"]),
                                    "pullback_transport_cosine": ratio(cross, math.sqrt(d2)*math.sqrt(t2)),
                                    "radial_element_groups": bins,
                                    "largest_eulerian_elements": [{"element": e, "squared_integral": finite(r["eulerian_squared_integral"]),
                                       "fraction": ratio(finite(r["eulerian_squared_integral"]), e2)} for e, r in top]}
        totals_g["direction_displacement_l2"] = math.sqrt(totals_g["direction_displacement_squared_integral"])
        totals_g["direction_displacement_rms"] = math.sqrt(totals_g["direction_displacement_squared_integral"]/totals_g["physical_volume"])
        outputs.append({"quadrature": q, "geometry": totals_g, "affine": totals_a, "fields": field_outputs})
    for metric, expected in (("accepted_trial_prepared_normalization_unchanged", 1.), ("production_mapping_cache_unchanged", 1.),
                             ("trial_preparations", 0.), ("volume_direction_actions", 1.)):
        compare(0, "invariant", metric, expected)
    for metric in ("maximum_scatter_roundtrip_error", "maximum_surface_trace_error"):
        key = (0, "reconstruction", metric)
        expected_summary.add(key)
        audit.check("physical:"+metric, 0 <= observed[key] <= 512*EPS*max(1., float(np.max(np.abs(physical))), float(np.max(np.abs(physical_direction)))))
    audit.check("physical:summary_coverage", set(observed) == expected_summary)
    sensitivity = {}
    for field in FIELDS:
        sensitivity[field] = {}
        for key in (*SQUARES, "pullback_transport_cross_integral"):
            a, b = (out["fields"][field]["integrals"][key] for out in outputs)
            sensitivity[field][key] = {"q10": a, "q14": b, "difference": b-a,
                                       "relative_to_q14": ratio(abs(b-a), abs(b))}
    return {"metadata": meta, "quadratures": outputs, "quadrature_sensitivity": sensitivity,
            "surface": surface_fit(surface, physical, physical_direction, blocks, audit),
            "scope": "recorded stellar broken-scalar integrals; no interface jumps, RT transport, exact shell integrals or continuum gauge test"}


def producer_binding(direction_path, folder, audit):
    producer = direction_path.parent
    mapping = {"nonlinear_step_physical_base.txt":"physical_base", "nonlinear_step_normalized_base.txt":"base",
               "nonlinear_step_state_factors.txt":"state_factors", "nonlinear_step_residual_factors.txt":"residual_factors",
               "nonlinear_step_residual.txt":"residual", "nonlinear_step_action_full.txt":"action"}
    present = [name for name in mapping if (producer/name).exists()]
    if not present:
        return {"status":"producer_neighbors_unavailable", "bound":False}
    if len(present) != len(mapping) or not (producer/"nonlinear_step_layout.csv").exists():
        raise ValueError("Incomplete nonlinear-step direction producer artifacts")
    saved_hashes, equal = {}, {}
    for name, current in mapping.items():
        path=producer/name
        saved_hashes[name]=digest(path)
        a,b=vector(path),vector(folder/(P+current+".txt"))
        exact=np.array_equal(a,b)
        equal[name]=exact
        audit.check("producer:vector:"+name, exact if "factors" in name else vector_near(a,b))
    path=producer/"nonlinear_step_layout.csv"
    saved_hashes[path.name]=digest(path)
    audit.check("producer:layout", rows(path)==rows(folder/(P+"layout.csv")))
    run_path=producer.parent/(producer.name+"_run.json")
    run_bound=False
    if run_path.exists():
        run=json.loads(run_path.read_text())
        for name, actual in saved_hashes.items():
            audit.check("producer:recorded_output_hash:"+name, run["output_hashes"].get(name)==actual)
        audit.check("producer:direction_output_hash", run["output_hashes"].get(direction_path.name)==digest(direction_path))
        saved_hashes[run_path.name]=digest(run_path)
        run_bound=True
    return {"status":"producer_neighbors_compared", "bound":run_bound, "producer":str(producer),
            "hashes":saved_hashes,"vectors_bitwise_equal":equal,
            "qualification":"neighbor hashes checked against producer run when present; no FEM rerun"}


def provenance(folder, required, audit):
    root = folder.parent
    source_file, run_file = root / "provenance.json", root / (folder.name + "_run.json")
    if not source_file.exists() and not run_file.exists():
        return {"status": "unavailable", "bound": False}
    if not source_file.exists() or not run_file.exists():
        raise ValueError("Partially present campaign provenance")
    source, run = (json.loads(p.read_text()) for p in (source_file, run_file))
    audit.check("provenance:case_id", run["id"] == folder.name)
    command = run["command"]
    final_mode = "--weak-shape-final" in command
    audit.check("provenance:completed_probe_exit", run["exit_code"] == 0 or (final_mode and run["exit_code"] == 3))
    audit.check("provenance:required_output_hashes", required <= set(run["output_hashes"]))
    for name, expected in run["output_hashes"].items():
        path = folder / name
        audit.check("provenance:contained_output:"+name, path.resolve().is_relative_to(folder.resolve()))
        audit.check("provenance:output_hash:"+name, digest(path) == expected)
    for name, expected in run.get("input_state_hashes", {}).items():
        audit.check("provenance:input_hash:"+name, digest(Path(name)) == expected)
    case = next(c for c in source["cases"] if c["id"] == folder.name)
    expected = [str(root / "plateau_diagnostic"), *case["arguments"], "--mesh", source["mesh"], "--output", str(folder)]
    audit.check("provenance:command", command == expected)
    audit.check("provenance:binary_hash", digest(Path(command[0])) == source["binary_sha256"])
    audit.check("provenance:mesh_hash", digest(Path(source["mesh"])) == source["mesh_sha256"])
    manifest = root / "cases.json"
    if manifest.exists():
        audit.check("provenance:case_manifest", json.loads(manifest.read_text()) == source["cases"])
    meta_path=folder/"metadata.txt"
    if meta_path.exists():
        meta=dict(line.split("=",1) for line in meta_path.read_text().splitlines() if "=" in line)
        audit.check("provenance:mode", ("--weak-shape-probe" in command) != final_mode)
        audit.check("provenance:metadata_mode", meta["weak_shape_probe"]==str(int(not final_mode)) and meta["weak_shape_final"]==str(int(final_mode)))
        audit.check("provenance:metadata_scope", meta["weak_shape_direction_space"]=="normalized_current_factors" and
                    meta["weak_shape_operator"]=="full_monolithic" and meta["weak_shape_trial_policy"]=="raw_prepare_only")
        for option,key in (("--linear-tol","linear_tol"),("--restart","restart"),("--max-linear","max_linear"),
                           ("--max-newton","max_newton"),("--index","n"),("--angular-momentum","J")):
            if option in command:
                audit.check("provenance:metadata_argument:"+option, finite(command[command.index(option)+1])==finite(meta[key]))
        if final_mode:
            audit.check("provenance:reported_convergence", integer(meta["converged"])==int(run["exit_code"]==0))
            audit.check("provenance:final_accepted_state", np.array_equal(vector(folder/"accepted_state.txt"),vector(folder/(P+"physical_base.txt"))))
    producer = None
    for option, saved in (("--probe-state", "physical_base"), ("--probe-direction", "direction")):
        if option in command:
            input_path = Path(command[command.index(option)+1])
            if not input_path.is_absolute():
                input_path = Path(source["source"]) / input_path
            input_path = input_path.resolve()
            audit.check("provenance:input_bound:"+option, str(input_path) in run.get("input_state_hashes", {}))
            audit.check("provenance:input_vector:"+option, np.array_equal(vector(input_path), vector(folder/(P+saved+".txt"))))
            if option == "--probe-direction":
                producer = producer_binding(input_path,folder,audit)
    return {"status": "runner_hashes_checked", "bound": True, "source_commit": source["commit"],
            "provenance_sha256": digest(source_file), "run_sha256": digest(run_file), "binary_sha256": source["binary_sha256"],
            "direction_producer":producer}


def inspect(folder):
    folder = Path(folder).resolve()
    audit, required, observations = Audit(), set(), {}
    try:
        def get(name):
            filename = P + name + ".txt"
            required.add(filename)
            return vector(folder / filename)
        required.update(P+n+".csv" for n in ("summary", "layout", "blocks", "trials"))
        summary = {k: finite(r["value"]) for k, r in table(folder/(P+"summary.csv"), lambda r:r["metric"]).items()}
        for key, value in {"schema_version": 1, "accepted_commits": 0, "probe_exception": 0,
                           "accepted_buffers_bitwise_identical": 1, "normalization_factors_bitwise_identical": 1,
                           "direction_buffers_bitwise_identical": 1, "maximum_restore_error": 0,
                           "restored_action_error": 0}.items():
            audit.check("summary:"+key, summary[key] == value)
        x, xp, f, p, pp, jp, defect, factors, rf = [get(n) for n in
            ("base", "physical_base", "residual", "direction", "physical_direction", "action", "defect", "state_factors", "residual_factors")]
        n = x.size
        audit.check("vector_sizes", all(v.size == n for v in (xp,f,p,pp,jp,defect,factors,rf)))
        audit.check("summary:state_size", summary["state_size"] == n)
        audit.check("positive_factors", np.all(factors > 0) and np.all(rf > 0))
        audit.check("denormalized_base", vector_near(xp, x/factors))
        audit.check("denormalized_direction", vector_near(pp, p/factors))
        audit.check("defect_identity", vector_near(defect, f+jp, norm(f)+norm(jp)))
        audit.check("restored_residual_exact", np.array_equal(f, get("restored_residual")))
        audit.check("restored_action_exact", np.array_equal(jp, get("restored_action")))
        fn, pn, jn, dn = map(norm, (f,p,jp,defect))
        if fn <= 0 or pn <= 0:
            raise ValueError("Probe requires nonzero base residual and direction")
        threshold = max(summary["outer_absolute_tolerance"], summary["outer_tolerance"]*fn)
        audit.check("valid_linear_controls", 0 < summary["outer_tolerance"] < 1 and summary["outer_absolute_tolerance"] >= 0 and
                    integer(summary["maximum_iterations"]) > 0 and integer(summary["restart"]) > 0)
        for key, value in {"base_norm":fn, "direction_norm":pn, "action_norm":jn, "verified_defect_norm":dn,
                           "verified_relative_defect":dn/fn, "verified_threshold":threshold}.items():
            audit.close("summary:"+key, summary[key], value)
        audit.check("verified_gate_truth", summary["verified_gate"] in (0,1) and bool(summary["verified_gate"]) == (dn <= threshold))
        audit.check("supplied_direction_flag", summary["supplied_direction"] in (0,1))
        if not summary["supplied_direction"]:
            audit.close("linear_reported_true_norm", summary["linear_reported_true_norm"], dn, fn+jn)
            audit.check("linear_iterations", 0 <= integer(summary["linear_iterations"]) <= integer(summary["maximum_iterations"]))
        blocks = layout(folder/(P+"layout.csv"), n)
        actions = [get("action_b"+str(i)) for i in range(len(blocks["value"]))]
        audit.check("split_column_count", summary["split_columns"] == len(actions))
        audit.check("direct_action_count", summary["direct_jacobian_applications"] == len(actions)+2)
        audit.check("preparation_count", summary["preparations"] == 6 + int(not summary["supplied_direction"]))
        audit.check("finite_nonnegative_probe_time", summary["total_probe_seconds"] >= 0)
        total = np.sum(actions, axis=0)
        action_scale = math.fsum(norm(a) for a in actions)
        audit.check("split_action_sum", vector_near(total, jp, action_scale))
        audit.close("reported_split_error", summary["split_action_sum_error"], norm(total-jp), action_scale)
        reported = table(folder/(P+"blocks.csv"), lambda r:(integer(r["column"]),r["residual_block"]))
        audit.check("block_table_coverage", set(reported) == {(i,r[0]) for i in range(len(actions)) for r in blocks["residual"]})
        for i,(name,start,count) in enumerate(blocks["value"]):
            for rname,offset,size in blocks["residual"]:
                row = reported[(i,rname)]
                audit.check(f"block_metadata:{i}:{rname}", row["value_block"]==name and integer(row["value_offset"])==start and
                            integer(row["value_size"])==count and integer(row["residual_offset"])==offset and integer(row["residual_size"])==size)
                audit.close(f"block_direction_norm:{i}:{rname}", finite(row["direction_norm"]), norm(p[start:start+count]))
                audit.close(f"block_action_norm:{i}:{rname}", finite(row["action_norm"]), norm(actions[i][offset:offset+size]))
        svd, smallest_action = restricted_svd(p, actions, blocks["value"], audit)
        svd["smallest_action_residual_blocks"] = {name:norm(smallest_action[o:o+s]) for name,o,s in blocks["residual"]}
        trial_rows = table(folder/(P+"trials.csv"), lambda r:integer(r["trial"]))
        audit.check("trial_coverage", set(trial_rows)=={0,1,2})
        trials = []
        for i,alpha in enumerate((1/128,1/512,1/2048)):
            row = trial_rows[i]
            audit.check(f"trial_alpha:{i}", finite(row["alpha"])==alpha)
            audit.check(f"trial_prepared_flag:{i}", row["prepared"] in ("0","1"))
            audit.check(f"trial_candidate:{i}", vector_near(get("candidate_a"+str(i)), x+alpha*p, norm(x)+alpha*pn))
            prediction = f+alpha*jp
            audit.close(f"trial_prediction_norm:{i}", finite(row["linear_prediction_norm"]), norm(prediction))
            prepared = row["prepared"]=="1"
            entry = {"trial":i, "alpha":alpha, "prepared":prepared, "reason":integer(row["reason"]), "stage":integer(row["stage"])}
            if prepared:
                audit.check(f"trial_prepared_reasons:{i}", entry["reason"] == entry["stage"] == -1)
                ft = get("residual_a"+str(i)); remainder = ft-prediction
                values = {"residual_norm":norm(ft), "remainder_norm":norm(remainder), "remainder_over_alpha":norm(remainder)/alpha}
                for key,value in values.items():
                    audit.close(f"trial:{i}:{key}", finite(row[key]), value, (fn+alpha*jn+norm(ft)) / (alpha if key=="remainder_over_alpha" else 1))
                entry.update(values)
                entry.update(residual_over_base=norm(ft)/fn, remainder_over_alpha_squared=norm(remainder)/alpha**2,
                             linear_dot_remainder=float(checked_einsum("i,i->",prediction,remainder)))
            else:
                audit.check(f"trial_failed_reasons:{i}", entry["reason"]>=0 and entry["stage"]>=0)
                audit.check(f"trial_failed_vector_absent:{i}", not (folder/(P+"residual_a"+str(i)+".txt")).exists())
                audit.check(f"trial_failed_metrics:{i}", all(math.isnan(float(row[k])) for k in ("residual_norm","remainder_norm","remainder_over_alpha")))
            trials.append(entry)
        physical_data = physical_analysis(folder,xp,pp,blocks["value"],audit,required)
        if (folder/"metadata.txt").exists():
            required.add("metadata.txt")
            meta=dict(line.split("=",1) for line in (folder/"metadata.txt").read_text().splitlines() if "=" in line)
            for key,source in (("outer_tolerance","linear_tol"),("restart","restart"),("maximum_iterations","max_linear")):
                audit.check("metadata:linear_control:"+key,summary[key]==finite(meta[source]))
            if meta.get("weak_shape_final")=="1":
                required.add("accepted_state.txt")
        binding = provenance(folder,required,audit)
        observations = {"summary":summary, "provenance":binding,
                        "linear_classification":"verified_newton_direction" if dn<=threshold else "finite_direction_failed_true_gate",
                        "verified_linear_gate":dn<=threshold, "restricted_block_span_svd":svd, "trials":trials,
                        "block_action_cancellation_factor":ratio(action_scale,jn), "physical":physical_data,
                        "input_vector_hashes":{name:digest(folder/(P+name+".txt")) for name in ("physical_base","base","direction","state_factors","residual_factors")},
                        "layout_sha256":digest(folder/(P+"layout.csv"))}
    except (OSError, ValueError, KeyError, TypeError, IndexError, StopIteration, OverflowError, ZeroDivisionError, np.linalg.LinAlgError) as error:
        audit.check("complete_consistent_artifacts: " + str(error),False)
    passed = not audit.failures
    return {"directory":str(folder), "artifact_pass":passed, "observations_valid":passed,
            "checks_count":audit.count, "failures":audit.failures,
            "required_artifacts":sorted(required), "observations":observations if passed else {}}


def comparisons(cases):
    output=[]
    for i,left in enumerate(cases):
        for right in cases[i+1:]:
            if not left["artifact_pass"] or not right["artifact_pass"]:
                continue
            a,b=left["observations"],right["observations"]
            same_layout=a["layout_sha256"]==b["layout_sha256"]
            same={k:a["input_vector_hashes"][k]==b["input_vector_hashes"][k] for k in a["input_vector_hashes"]}
            output.append({"left":left["directory"],"right":right["directory"],"same_layout":same_layout,"same_saved_vectors":same,
                           "same_state_metric_comparison":same_layout and all(same[k] for k in ("physical_base","base","state_factors","residual_factors")),
                           "qualification":"different discretizations or factors require physical comparisons; raw coefficient norms are not comparable"})
    return output


def tiny_fixture(folder, failed_gate=False):
    """Construct a complete eight-coordinate synthetic producer, with no FEM."""
    folder.mkdir(parents=True,exist_ok=True)
    def save(name,v):
        (folder/(P+name+".txt")).write_text(str(len(v))+"\n"+"\n".join(format(float(x),".17g") for x in v)+"\n")
    def write(name,records):
        with (folder/name).open("w",newline="") as stream:
            out=csv.DictWriter(stream,fieldnames=list(records[0]));out.writeheader();out.writerows(records)
    sf=np.arange(1.,9.);xp=np.array([2.,0.,0.,0.,0.,0.,0.,1.]);x=xp*sf
    p=np.arange(1.,9.)*1e-4;f=-p;jp=p*(.9 if failed_gate else 1.)
    for name,v in {"base":x,"physical_base":xp,"residual":f,"direction":p,"physical_direction":p/sf,
                   "action":jp,"defect":f+jp,"state_factors":sf,"residual_factors":np.ones(8),
                   "restored_residual":f,"restored_action":jp}.items():save(name,v)
    blocks=[("density",0,1),("surface_deformation",1,6),("constant",7,1)]
    write(P+"layout.csv",[{"kind":kind,"block":name,"offset":o,"size":n} for kind in ("value","residual") for name,o,n in blocks])
    br=[]
    for i,(name,o,n) in enumerate(blocks):
        a=np.zeros(8);a[o:o+n]=jp[o:o+n];save("action_b"+str(i),a)
        for rn,ro,rs in blocks:
            br.append(dict(column=i,value_block=name,value_offset=o,value_size=n,residual_block=rn,residual_offset=ro,
                           residual_size=rs,direction_norm=norm(p[o:o+n]),action_norm=norm(a[ro:ro+rs])))
    write(P+"blocks.csv",br)
    summary=dict(schema_version=1,state_size=8,supplied_direction=1,accepted_commits=0,outer_tolerance=.03,
                 outer_absolute_tolerance=0,maximum_iterations=200,restart=160,base_norm=norm(f),direction_norm=norm(p),
                 action_norm=norm(jp),verified_defect_norm=norm(f+jp),verified_relative_defect=norm(f+jp)/norm(f),
                 verified_threshold=.03*norm(f),verified_gate=int(not failed_gate),split_columns=3,split_action_sum_error=0,
                 restored_action_error=0,probe_exception=0,accepted_buffers_bitwise_identical=1,normalization_factors_bitwise_identical=1,
                 direction_buffers_bitwise_identical=1,maximum_restore_error=0,preparations=6,direct_jacobian_applications=5,total_probe_seconds=1.)
    write(P+"summary.csv",[dict(metric=k,value=v) for k,v in summary.items()])
    tr=[]
    for i,a in enumerate((1/128,1/512,1/2048)):
        save("candidate_a"+str(i),x+a*p);pred=f+a*jp;ft=pred+a*a*np.ones(8)*1e-5;save("residual_a"+str(i),ft)
        tr.append(dict(trial=i,alpha=a,prepared=1,reason=-1,stage=-1,residual_norm=norm(ft),remainder_norm=norm(ft-pred),
                       remainder_over_alpha=norm(ft-pred)/a,linear_prediction_norm=norm(pred)))
    write(P+"trials.csv",tr)
    meta=dict(schema="weak_shape_physical_v1",dg_derivatives="broken DG derivatives",gravity_vector="not evaluated",
              density_family_order="2",enthalpy_family_order="3",potential_family_order="2",displacement_family_order="3")
    write(PP+"metadata.csv",[dict(key=k,value=v) for k,v in meta.items()])
    elements=[];geometry=[];affine=[];ps=[]
    for q in (10,14):
        g=dict(quadrature=q,element=0,attribute=1,geometry=5,points=8,reference_mesh_geometry_order=4,density_element_order=2,
               enthalpy_element_order=3,potential_element_order=2,displacement_element_order=3,reference_radius_min=.1,
               reference_radius_max=.9,physical_radius_min=.1,physical_radius_max=.9,physical_volume=1.,
               base_displacement_squared_integral=0.,direction_displacement_squared_integral=.2)
        geometry.append(g)
        ar=dict(quadrature=q,element=0,points=8,reference_volume=1.,value_error_squared_integral=.01,
                reference_value_squared_integral=1.,gradient_error_squared_integral=.03,reference_gradient_squared_integral=3.)
        affine.append(ar)
        for field in FIELDS:
            vals=dict(zip(INTEGRALS,(2.,3.,2.,1.,4.,9.,4.,1.,6.)))
            elements.append(dict(quadrature=q,element=0,attribute=1,geometry=5,field=field,
                                 gradient_kind="continuous_H1" if field=="specific_enthalpy" else "broken_DG",points=8,physical_volume=1.,**vals))
            measured={k:vals[k] for k in (*SQUARES,"pullback_transport_cross_integral")}
            measured.update(pullback_l2_norm=3.,transport_l2_norm=2.,eulerian_l2_norm=1.,squared_norm_identity_error=0.)
            ps.extend(dict(quadrature=q,field=field,metric=k,value=v) for k,v in measured.items())
        ps.extend(dict(quadrature=q,field="geometry",metric=k,value=v) for k,v in
                  dict(stellar_elements=1,physical_volume=1.,base_displacement_squared_integral=0.,direction_displacement_squared_integral=.2).items())
        am={k:v for k,v in ar.items() if k not in ("quadrature","element","points")}
        am.update(relative_value_l2_error=.1,relative_gradient_l2_error=.1)
        ps.extend(dict(quadrature=q,field="affine",metric=k,value=v) for k,v in am.items())
    ps.extend(dict(quadrature=0,field="invariant",metric=k,value=v) for k,v in
              dict(accepted_trial_prepared_normalization_unchanged=1,production_mapping_cache_unchanged=1,trial_preparations=0,volume_direction_actions=1).items())
    ps.extend(dict(quadrature=0,field="reconstruction",metric=k,value=0.) for k in
              ("maximum_scatter_roundtrip_error","maximum_surface_trace_error"))
    for name,data in (("elements",elements),("geometry",geometry),("affine",affine),("summary",ps)):write(PP+name+".csv",data)
    xyz=[[1.,0.,0.],[-1.,0.,0.],[0.,1.,0.],[0.,-1.,0.],[0.,0.,1.],[0.,0.,-1.]]
    write(PP+"surface_nodes.csv",[dict(parameter=i,volume_scalar_true_dof=i,unit_radial_x=v[0],unit_radial_y=v[1],
          unit_radial_z=v[2],reference_radius=1.,physical_radius=1.,base_parameter=0.,physical_direction_parameter=p[i+1]/sf[i+1]) for i,v in enumerate(xyz)])


def self_test():
    """Tiny algebraic controls only; never opens campaign artifacts."""
    import tempfile
    count=0
    def require(condition):
        nonlocal count
        assert condition
        count+=1
    audit=Audit()
    p=np.array([2.,3.,0.])
    actions=[np.array([2.,0.,0.]),np.array([-3.,3e-8,0.]),np.zeros(3)]
    out,action=restricted_svd(p,actions,[("a",0,1),("b",1,1),("zero",2,1)],audit)
    require(not audit.failures)
    require(near(out["singular_values"][-1],1e-8/math.sqrt(2),1e-8))
    require(out["active_blocks"]==["a","b"])
    require(norm(action)<1e-7)
    bad=Audit();restricted_svd(p,actions[:2]+[np.ones(3)],[("a",0,1),("b",1,1),("zero",2,1)],bad)
    require(bool(bad.failures))
    require(not vector_near(np.array([1.,2.]),np.array([1.,2.001])))
    require(near(1e-20,0.,1.))
    # This size triggers spurious Accelerate matmul flags on the investigation
    # host. Use independent scalar sums, not another BLAS product, as the oracle.
    a=np.linspace(-1.,1.,257*6).reshape(257,6);v=np.arange(1.,7.)
    expected=np.array([math.fsum(float(x*y) for x,y in zip(row,v)) for row in a])
    with np.errstate(divide="raise",over="raise",invalid="raise"):
        require(vector_near(checked_einsum("ij,j->i",a,v),expected,norm(expected)))
        b=np.arange(18.).reshape(6,3)/18.
        expected_matrix=np.array([[math.fsum(float(a[i,k]*b[k,j]) for k in range(6)) for j in range(3)] for i in range(257)])
        require(vector_near(checked_einsum("ij,jk->ik",a,b),expected_matrix,norm(expected_matrix)))
    try:checked_einsum("i,i->",np.array([math.inf]),np.ones(1));raise AssertionError("nonfinite contraction accepted")
    except ValueError:require(True)
    with tempfile.TemporaryDirectory() as directory:
        d=Path(directory)
        (d/"v").write_text("3\n1\n2\n3\n")
        require(np.array_equal(vector(d/"v"),np.array([1.,2.,3.])))
        for text in ("2\n1\n", "2\n1\nnan\n", "2\n1\ninf\n"):
            (d/"v").write_text(text)
            try:vector(d/"v");raise AssertionError("invalid vector accepted")
            except ValueError:require(True)
        (d/"layout").write_text("kind,block,offset,size\nvalue,a,0,2\nresidual,b,0,2\n")
        require(layout(d/"layout",2)["value"]==[("a",0,2)])
        (d/"layout").write_text("kind,block,offset,size\nvalue,a,1,2\nresidual,b,0,2\n")
        try:layout(d/"layout",2);raise AssertionError("bad layout accepted")
        except ValueError:require(True)
    xyz=np.array([[1.,0.,0.],[-1.,0.,0.],[0.,1.,0.],[0.,-1.,0.],[0.,0.,1.],[0.,0.,-1.]])
    values=2.+3.*xyz[:,0]
    surf={i:{"volume_scalar_true_dof":str(i),"unit_radial_x":str(v[0]),"unit_radial_y":str(v[1]),"unit_radial_z":str(v[2]),
             "reference_radius":"1","physical_radius":"1","base_parameter":"0","physical_direction_parameter":str(values[i])} for i,v in enumerate(xyz)}
    audit=Audit();fits=surface_fit(surf,np.zeros(6),values,[("surface_deformation",0,6)],audit)
    require(not audit.failures)
    require(fits["fits"][1]["relative_nodal_l2_remainder"]<1e-14)
    require(fits["fits"][0]["relative_nodal_l2_remainder"]>.1)
    with tempfile.TemporaryDirectory() as directory:
        d=Path(directory).resolve()/"case"
        tiny_fixture(d)
        result=inspect(d)
        assert result["artifact_pass"],result["failures"]
        require(result["observations"]["verified_linear_gate"])
        require(result["observations"]["provenance"]["status"]=="unavailable")
        require(result["observations"]["physical"]["quadratures"][0]["fields"]["density"]["eulerian_over_pullback_plus_transport_l2"]==.2)
        tiny_fixture(d,failed_gate=True)
        result=inspect(d)
        require(result["artifact_pass"] and result["observations"]["linear_classification"]=="finite_direction_failed_true_gate")
        for name in ("defect.txt","restored_action.txt","candidate_a1.txt","physical_base.txt","action_b1.txt"):
            tiny_fixture(d)
            path=d/(P+name);lines=path.read_text().splitlines();lines[1]=str(float(lines[1])+.001);path.write_text("\n".join(lines)+"\n")
            require(not inspect(d)["artifact_pass"])
        for filename,key in (("elements","pullback_transport_cross_integral"),("summary","value"),("geometry","physical_volume")):
            tiny_fixture(d);path=d/(PP+filename+".csv");data=rows(path);data[0][key]=str(float(data[0][key])+.1)
            with path.open("w",newline="") as stream:
                writer=csv.DictWriter(stream,fieldnames=list(data[0]));writer.writeheader();writer.writerows(data)
            require(not inspect(d)["artifact_pass"])
        for missing in (False,True):
            tiny_fixture(d);path=d/(PP+"summary.csv");data=rows(path)
            if missing:data=[r for r in data if r["metric"]!="production_mapping_cache_unchanged"]
            else:
                next(r for r in data if r["metric"]=="production_mapping_cache_unchanged")["value"]="0"
            with path.open("w",newline="") as stream:
                writer=csv.DictWriter(stream,fieldnames=list(data[0]));writer.writeheader();writer.writerows(data)
            require(not inspect(d)["artifact_pass"])
        tiny_fixture(d)
        path=d/(P+"trials.csv");data=rows(path);data[0].update(prepared="0",reason="2",stage="3",
            residual_norm="nan",remainder_norm="nan",remainder_over_alpha="nan")
        with path.open("w",newline="") as stream:
            writer=csv.DictWriter(stream,fieldnames=list(data[0]));writer.writeheader();writer.writerows(data)
        (d/(P+"residual_a0.txt")).unlink()
        require(inspect(d)["artifact_pass"])
        tiny_fixture(d)
        root=d.parent
        executable=root/"plateau_diagnostic";executable.write_text("synthetic binary fixture")
        mesh=root/"mesh";mesh.write_text("synthetic mesh fixture")
        case={"id":"case","arguments":["--weak-shape-probe"]}
        prov={"commit":"0"*40,"source":str(root),"cases":[case],"mesh":str(mesh),"mesh_sha256":digest(mesh),"binary_sha256":digest(executable)}
        (root/"provenance.json").write_text(json.dumps(prov))
        run={"id":"case","exit_code":0,"input_state_hashes":{},"command":[str(executable),"--weak-shape-probe","--mesh",str(mesh),"--output",str(d)],
             "output_hashes":{p.name:digest(p) for p in d.iterdir()}}
        (root/"case_run.json").write_text(json.dumps(run))
        result=inspect(d)
        assert result["artifact_pass"],result["failures"]
        require(result["observations"]["provenance"]["bound"])
        run["output_hashes"][P+"action.txt"]="0"*64
        (root/"case_run.json").write_text(json.dumps(run))
        require(not inspect(d)["artifact_pass"])
        producer=root/"producer";producer.mkdir()
        for old,new in {"physical_base":"physical_base","normalized_base":"base","state_factors":"state_factors",
                        "residual_factors":"residual_factors","residual":"residual","action_full":"action","direction_full":"direction"}.items():
            (producer/("nonlinear_step_"+old+".txt")).write_bytes((d/(P+new+".txt")).read_bytes())
        (producer/"nonlinear_step_layout.csv").write_bytes((d/(P+"layout.csv")).read_bytes())
        (root/"producer_run.json").write_text(json.dumps({"output_hashes":{p.name:digest(p) for p in producer.iterdir()}}))
        audit=Audit();bound=producer_binding(producer/"nonlinear_step_direction_full.txt",d,audit)
        require(bound["bound"] and not audit.failures)
        (producer/"nonlinear_step_state_factors.txt").write_text("8\n"+"2\n"*8)
        audit=Audit();producer_binding(producer/"nonlinear_step_direction_full.txt",d,audit)
        require(bool(audit.failures))
    print(f"{count} tiny synthetic checks passed")


def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument("directories",nargs="*",type=Path)
    parser.add_argument("--output",type=Path)
    parser.add_argument("--self-test",action="store_true")
    args=parser.parse_args()
    if args.self_test:
        self_test();return 0
    if not args.directories or args.output is None:
        parser.error("Provide one or more probe directories and --output")
    destination=args.output.resolve()
    if destination.exists() or any(destination.is_relative_to(d.resolve()) for d in args.directories):
        parser.error("Use a fresh output JSON outside the immutable probe directories")
    cases=[inspect(d) for d in args.directories]
    result={"schema_version":1,"passed":all(c["artifact_pass"] for c in cases),"checker_sha256":digest(Path(__file__)),
            "numpy_version":np.__version__,"tolerances":{"relative":RTOL,"cancellation_roundoff_multiplier":2048,"absolute_floor":1e-28},
            "scope":"Independent recorded-vector algebra/hash and physical-integral identity audit; no FEM rerun, quadrature rerun, full-rank, gauge, convergence or physical-model certificate",
            "cases":cases,"comparisons":comparisons(cases)}
    destination.parent.mkdir(parents=True,exist_ok=True)
    destination.write_text(json.dumps(result,indent=2,allow_nan=False)+"\n")
    print(json.dumps({"passed":result["passed"],"cases":len(cases),"output":str(destination)}))
    return 0 if result["passed"] else 1


if __name__=="__main__":
    sys.exit(main())