#!/usr/bin/env python3
"""Generate a structure sanity card — a pre-interpretation CIF sanity check.

A space group label is a question, not a verdict. This card checks what is
encoded in a CIF: geometry, symmetry sensitivity, coordination, bond lengths,
and reference-structure displacement. It is deliberately not a stability,
synthesis, or property predictor.

The card catches what symmetry analysis alone cannot: coordinate permutation,
overlapping atoms, and systematic geometric corruption that hides behind a
valid-looking space group label. The reference-structure matcher compares each
atom to its ideal position in the symmetry-refined cell — atoms at wrong
symmetry-equivalent positions show up as large displacements even when the
space group is correct.

Usage:
    python structure_sanity_card.py <cif> [--output card.md] [--json report.json]
"""
from __future__ import annotations

import argparse
import json
from pathlib import Path
from itertools import combinations

import numpy as np
from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer

DEFAULT_SYMPRECS = [0.01, 0.02, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5, 0.75, 1.0]


def fine_symmetry_sweep(structure: Structure, symprecs: list[float]) -> list[dict]:
    """Sweep symprec from tight to loose, recording space group at each point."""
    results = []
    for sp in symprecs:
        try:
            sga = SpacegroupAnalyzer(structure, symprec=sp)
            results.append({
                "symprec": sp,
                "spacegroup": sga.get_space_group_symbol(),
                "number": sga.get_space_group_number(),
                "crystal_system": sga.get_crystal_system(),
            })
        except Exception:
            results.append({"symprec": sp, "spacegroup": "ERROR", "number": -1, "crystal_system": "unknown"})
    return results


def min_pair_distance(structure: Structure) -> dict:
    """Find the shortest interatomic distance and identify the pair."""
    if len(structure) < 2:
        return {"distance": None, "sites": None}
    d_min, pair = float("inf"), None
    for i, j in combinations(range(len(structure)), 2):
        d = structure.get_distance(i, j)
        if d < d_min:
            d_min = d
            pair = (i, j, structure[i].specie.symbol, structure[j].specie.symbol)
    return {"distance": round(float(d_min), 4), "sites": pair}


def coordination_fingerprint(structure: Structure, cutoff: float = 3.0) -> list[dict]:
    """Per-site coordination with nearest-neighbor shell detection."""
    fingerprints = []
    for i, site in enumerate(structure):
        neighbors = structure.get_neighbors(site, cutoff)
        neighbor_list = sorted([(n.specie.symbol, n.nn_distance) for n in neighbors], key=lambda x: x[1])
        if neighbor_list:
            d_min = neighbor_list[0][1]
            coord_shell = [(e, d) for e, d in neighbor_list if d <= max(d_min * 1.3, d_min + 0.8)]
        else:
            coord_shell = []
        elem_counts = {}
        for elem, dist in coord_shell:
            elem_counts.setdefault(elem, []).append(round(float(dist), 4))
        fingerprints.append({"site": i, "element": site.specie.symbol, "cn": len(coord_shell), "neighbors": elem_counts})
    return fingerprints


def bond_stats(structure: Structure, cutoff: float = 3.0) -> dict:
    """Collect bond length statistics per element pair."""
    pairs = {}
    for site_i in structure:
        for n in structure.get_neighbors(site_i, cutoff):
            pair = tuple(sorted([site_i.specie.symbol, n.specie.symbol]))
            pairs.setdefault(pair, []).append(round(float(n.nn_distance), 4))
    stats = {}
    for pair, dists in pairs.items():
        arr = np.array(dists)
        label = f"{pair[0]}-{pair[1]}"
        stats[label] = {"count": len(dists), "min": float(arr.min()), "max": float(arr.max()),
                        "mean": float(arr.mean()), "std": float(arr.std())}
    return stats


def reference_match(structure: Structure, symprec: float = 0.1) -> dict:
    """Compare each atom to its ideal position in the symmetry-refined cell.

    Uses spglib.refine_cell to build the ideal structure from the detected
    space group, then measures per-atom displacement (minimum-image, PBC).
    Large displacements mean atoms sit at wrong symmetry-equivalent positions
    even though the space group label looks correct — the signature of
    coordinate permutation or systematic corruption that symmetry alone
    cannot catch.

    Returns trivial zeros for P1 (no symmetry to refine against).
    """
    try:
        import spglib
    except ImportError:
        return {"available": False}

    lattice = np.array(structure.lattice.matrix)
    positions = np.array(structure.frac_coords)
    numbers = np.array([site.specie.Z for site in structure])
    cell = (lattice, positions, numbers)

    ds = spglib.get_symmetry_dataset(cell, symprec=symprec)
    sg_symbol = ds.international if ds else "unknown"
    sg_number = ds.number if ds else -1

    refined = spglib.refine_cell(cell, symprec=symprec)
    if refined is None:
        return {"available": True, "spacegroup": sg_symbol, "number": sg_number, "error": "refine_cell failed"}

    lat_r, pos_r, num_r = refined
    if len(num_r) != len(structure):
        return {"available": True, "spacegroup": sg_symbol, "number": sg_number,
                "error": f"refined cell has {len(num_r)} atoms, expected {len(structure)}"}

    displacements = []
    for i in range(len(structure)):
        zi = structure[i].specie.Z
        best_d = float("inf")
        for j in range(len(num_r)):
            if num_r[j] == zi:
                d_frac = structure[i].frac_coords - pos_r[j]
                d_frac = d_frac - np.round(d_frac)
                d_cart = d_frac @ lattice
                d = float(np.linalg.norm(d_cart))
                if d < best_d:
                    best_d = d
        displacements.append(best_d)

    disps = np.array(displacements)
    per_element = {}
    for i, site in enumerate(structure):
        per_element.setdefault(site.specie.symbol, []).append(disps[i])

    elem_stats = {}
    for elem, vals in sorted(per_element.items()):
        arr = np.array(vals)
        elem_stats[elem] = {
            "count": len(arr),
            "max_disp": round(float(arr.max()), 4),
            "mean_disp": round(float(arr.mean()), 4),
            "n_above_0p5": int(np.sum(arr > 0.5)),
        }

    return {
        "available": True,
        "spacegroup": sg_symbol,
        "number": sg_number,
        "max_displacement": round(float(disps.max()), 4),
        "mean_displacement": round(float(disps.mean()), 4),
        "median_displacement": round(float(np.median(disps)), 4),
        "n_above_0p5": int(np.sum(disps > 0.5)),
        "n_above_1p0": int(np.sum(disps > 1.0)),
        "n_above_2p0": int(np.sum(disps > 2.0)),
        "n_atoms": len(disps),
        "per_element": elem_stats,
        "is_p1": sg_number == 1,
    }


def assess(report: dict) -> list[tuple[str, str]]:
    """Generate gate notes — PASS / CHECK / FAIL / NOTE for each check."""
    notes = []

    # 1. Minimum pair distance — the most important geometry check
    mpd = report["min_pair"]["distance"]
    sites = report["min_pair"]["sites"]
    if mpd is not None and mpd < 0.5:
        notes.append(("FAIL", f"Minimum pair distance {mpd:.3f} Å is physically impossible. "
                      f"Atoms {sites} are overlapping. Structure is corrupted."))
    elif mpd is not None and mpd < 1.0:
        notes.append(("CHECK", f"Minimum pair distance {mpd:.3f} Å is suspiciously short. "
                      f"Atoms {sites} may be overlapping or misplaced."))
    elif mpd is not None and mpd < 1.5:
        notes.append(("CHECK", f"Minimum pair distance {mpd:.3f} Å is short. Inspect atoms {sites}."))
    elif mpd is not None:
        notes.append(("PASS", f"Minimum pair distance {mpd:.3f} Å is physically reasonable."))

    # 2. Symmetry stability
    sweep = report["symmetry_sweep"]
    valid = [e for e in sweep if e["spacegroup"] != "ERROR"]
    distinct_sgs = set(e["spacegroup"] for e in valid)
    tight_sg = valid[0]["spacegroup"] if valid else "unknown"
    loose_sg = valid[-1]["spacegroup"] if valid else "unknown"

    if len(distinct_sgs) == 1:
        notes.append(("PASS", f"Symmetry is ROBUST: {tight_sg} at all tolerances "
                      f"({sweep[0]['symprec']:.2f}-{valid[-1]['symprec']:.2f} Å). "
                      f"The space group is not a tolerance artifact."))
    else:
        transitions = []
        prev = None
        for e in valid:
            if e["spacegroup"] != prev:
                transitions.append(f"{e['symprec']:.2f}Å→{e['spacegroup']}")
                prev = e["spacegroup"]
        notes.append(("CHECK", f"Symmetry is FRAGILE: transitions at {', '.join(transitions)}. "
                      f"Tight: {tight_sg}, loose: {loose_sg}. "
                      f"The tight-tolerance result may reflect noise, not true symmetry."))

    # 3. Occupancy
    if report["ordered"]:
        notes.append(("PASS", "All sites are ordered."))
    else:
        notes.append(("CHECK", "Partial/disordered occupancies need an explicit intended-disorder record."))

    # 4. Bond length dispersion — flag high coefficient of variation
    for label, s in report["bond_stats"].items():
        if s["count"] > 2 and s["mean"] > 0:
            cv = s["std"] / s["mean"]
            if cv > 0.5:
                notes.append(("CHECK", f"{label} bond lengths are highly dispersed "
                              f"(CV={cv:.2f}, range {s['min']:.3f}-{s['max']:.3f} Å). "
                              f"May indicate mixed coordination or structural corruption."))

    # 5. Reference-structure match (displacement from symmetry-refined ideal)
    rm = report.get("reference_match", {})
    if rm.get("available") and not rm.get("error"):
        n_bad = rm["n_above_0p5"]
        max_d = rm["max_displacement"]
        if rm.get("is_p1"):
            if max_d < 0.01:
                notes.append(("NOTE", "Reference match is trivial: structure is P1, so every "
                              "atom is its own Wyckoff position. The displacement check is "
                              "uninformative for P1 — rely on geometry and bond stats instead."))
        elif n_bad > 0:
            frac = n_bad / rm["n_atoms"]
            worst_elem = max(rm["per_element"].items(), key=lambda x: x[1]["n_above_0p5"])
            if frac > 0.3 or max_d > 1.0:
                notes.append(("FAIL", f"Reference-structure match: {n_bad}/{rm['n_atoms']} atoms "
                              f"displaced >0.5 Å from symmetry-refined ideal (max {max_d:.3f} Å). "
                              f"Worst element: {worst_elem[0]} ({worst_elem[1]['n_above_0p5']}/{worst_elem[1]['count']} "
                              f"displaced). Atoms are at wrong symmetry-equivalent positions — "
                              f"the CIF is corrupted despite the {rm['spacegroup']} label."))
            elif n_bad > 0:
                notes.append(("CHECK", f"Reference-structure match: {n_bad}/{rm['n_atoms']} atoms "
                              f"displaced >0.5 Å from ideal (max {max_d:.3f} Å). "
                              f"May indicate a distorted variant of {rm['spacegroup']} or mild corruption."))
        else:
            notes.append(("PASS", f"Reference-structure match: all {rm['n_atoms']} atoms within "
                          f"0.5 Å of symmetry-refined ideal positions (max {max_d:.4f} Å). "
                          f"No coordinate corruption detected."))
    elif rm.get("error"):
        notes.append(("NOTE", f"Reference-structure match unavailable: {rm['error']}"))

    # 6. P1-specific guidance
    if tight_sg == "P1":
        if len(distinct_sgs) == 1:
            notes.append(("NOTE", "P1 is robust across all tolerances. This is a genuinely "
                          "low-symmetry structure, not a tolerance artifact. "
                          "P1 here is an encoded feature, not a diagnosis."))
        else:
            notes.append(("NOTE", "P1 at tight tolerance but higher symmetry emerges at "
                          f"looser tolerance. The structure may be a distorted version of {loose_sg}. "
                          "Investigate the intended prototype."))
    elif mpd is not None and mpd < 1.0:
        notes.append(("NOTE", f"Despite {tight_sg} symmetry, the geometry reveals corruption. "
                      "The space group label is necessary but not sufficient — "
                      "always check minimum pair distances."))

    return notes


def generate_report(cif_path: Path, symprecs: list[float] | None = None) -> dict:
    if symprecs is None:
        symprecs = DEFAULT_SYMPRECS
    structure = Structure.from_file(str(cif_path))
    return {
        "source": str(cif_path),
        "formula": structure.composition.reduced_formula,
        "sites": len(structure),
        "volume_per_atom": round(float(structure.volume / len(structure)), 4),
        "density": round(float(structure.density), 4),
        "ordered": all(s.is_ordered for s in structure),
        "lattice": {
            "a": round(float(structure.lattice.a), 4),
            "b": round(float(structure.lattice.b), 4),
            "c": round(float(structure.lattice.c), 4),
            "alpha": round(float(structure.lattice.alpha), 2),
            "beta": round(float(structure.lattice.beta), 2),
            "gamma": round(float(structure.lattice.gamma), 2),
        },
        "min_pair": min_pair_distance(structure),
        "symmetry_sweep": fine_symmetry_sweep(structure, symprecs),
        "coordination": coordination_fingerprint(structure),
        "bond_stats": bond_stats(structure),
        "reference_match": reference_match(structure),
    }


def markdown_card(report: dict, notes: list[tuple[str, str]]) -> str:
    sg_line = ", ".join(f"{e['symprec']:.2f}Å: {e['spacegroup']}" for e in report["symmetry_sweep"])
    mpd = report["min_pair"]
    mpd_str = f"{mpd['distance']:.4f} Å" if mpd["distance"] else "n/a"

    rows = [
        ("Source", report["source"]),
        ("Formula", report["formula"]),
        ("Sites", str(report["sites"])),
        ("Volume / atom", f"{report['volume_per_atom']:.4f} ų"),
        ("Density", f"{report['density']:.4f} g cm⁻³"),
        ("Min pair distance", mpd_str),
        ("Ordered", str(report["ordered"])),
        ("Symmetry sweep", sg_line),
    ]
    table = "\n".join(f"| {k} | {v} |" for k, v in rows)
    gate_lines = "\n".join(f"- **[{level}]** {text}" for level, text in notes)

    bond_lines = "\n".join(
        f"| {label} | {s['count']} | {s['min']:.3f} | {s['max']:.3f} | {s['mean']:.3f} | {s['std']:.3f} |"
        for label, s in sorted(report["bond_stats"].items())
    )

    # Reference-structure match section
    rm = report.get("reference_match", {})
    rm_section = ""
    if rm.get("available") and not rm.get("error"):
        rm_rows = [
            ("Detected space group", f"{rm['spacegroup']} (#{rm['number']})"),
            ("Max displacement", f"{rm['max_displacement']:.4f} Å"),
            ("Mean displacement", f"{rm['mean_displacement']:.4f} Å"),
            ("Median displacement", f"{rm['median_displacement']:.4f} Å"),
            ("Atoms >0.5 Å from ideal", f"{rm['n_above_0p5']} / {rm['n_atoms']}"),
            ("Atoms >1.0 Å from ideal", f"{rm['n_above_1p0']} / {rm['n_atoms']}"),
            ("Atoms >2.0 Å from ideal", f"{rm['n_above_2p0']} / {rm['n_atoms']}"),
        ]
        rm_table = "\n".join(f"| {k} | {v} |" for k, v in rm_rows)
        elem_lines = "\n".join(
            f"| {elem} | {s['count']} | {s['max_disp']:.4f} | {s['mean_disp']:.4f} | {s['n_above_0p5']} |"
            for elem, s in rm["per_element"].items()
        )
        rm_section = f"""
## Reference-structure match

Compares each atom to its ideal position in the symmetry-refined cell
(spglib `refine_cell` at symprec=0.1 Å). Large displacements mean atoms sit
at wrong symmetry-equivalent positions even though the space group label
looks correct. Trivially zero for P1 (no symmetry to refine against).

| Field | Value |
| --- | --- |
{rm_table}

| Element | Count | Max disp (Å) | Mean disp (Å) | >0.5 Å |
| --- | --- | --- | --- | --- |
{elem_lines}
"""
    elif rm.get("error"):
        rm_section = f"\n## Reference-structure match\n\nUnavailable: {rm['error']}\n"
    elif not rm.get("available"):
        rm_section = "\n## Reference-structure match\n\nspglib not available.\n"

    return f"""# Structure sanity card — {report['formula']}

A space group label is a question, not a verdict. This card checks what is
encoded in the CIF: geometry, symmetry sensitivity, bond lengths, and
reference-structure displacement. It is not a stability, synthesis, or
property claim.

| Field | Observation |
| --- | --- |
{table}

## Gate notes

{gate_lines}

{rm_section}
## Bond length statistics

| Pair | Count | Min (Å) | Max (Å) | Mean (Å) | Std (Å) |
| --- | --- | --- | --- | --- | --- |
{bond_lines}
"""


def main() -> None:
    parser = argparse.ArgumentParser(description="Generate a structure sanity card")
    parser.add_argument("cif", type=Path)
    parser.add_argument("--output", type=Path, required=True)
    parser.add_argument("--json", type=Path)
    parser.add_argument("--symprec", type=float, nargs="+", default=DEFAULT_SYMPRECS)
    args = parser.parse_args()

    report = generate_report(args.cif, args.symprec)
    notes = assess(report)
    args.output.write_text(markdown_card(report, notes))
    if args.json:
        args.json.write_text(json.dumps({"report": report, "gate_notes": notes}, indent=2))


if __name__ == "__main__":
    main()




# ---------------------------------------------------------------------------
# v4: prototype-aware gate (curiosity-window addition, 2026-08-06/07)
#
# The species-swap experiment showed geometry gates and BVS are blind to label
# exchanges. The only one-call catch is a declared prototype: match the
# candidate against a species-carrying template. Antisite disorder is real
# physics, so a species mismatch is always INSPECT, never BROKEN.
#
# SCAR: templates must come from an independent Wyckoff implementation.
# pymatgen Structure.from_spacegroup("Fd-3m", ...) silently produced an
# overlapping, broken 32e oxygen orbit here; ASE's spacegroup tables produced
# the textbook orbit. Every template below is built via ASE and verified
# (space group + no unphysical close pairs) before use.
# ---------------------------------------------------------------------------

import numpy as np
from scipy.optimize import linear_sum_assignment
from pymatgen.core import Element, Lattice
from pymatgen.analysis.structure_matcher import StructureMatcher
from ase.spacegroup import crystal as _ase_crystal
from pymatgen.io.ase import AseAtomsAdaptor

# Roles are listed cation-first; elements are assigned in increasing
# electronegativity (least electronegative -> first role). Each entry carries
# the ASE crystal() arguments for the conventional cell.
PROTOTYPES = {
    "rocksalt":  {"sg": "Fm-3m", "sg_num": 225, "roles": ["A", "B"],
                  "basis": [(0, 0, 0), (0.5, 0.5, 0.5)], "a": 5.6, "setting": 1},
    "cesium_chloride": {"sg": "Pm-3m", "sg_num": 221, "roles": ["A", "B"],
                        "basis": [(0, 0, 0), (0.5, 0.5, 0.5)], "a": 4.1, "setting": 1},
    "zincblende": {"sg": "F-43m", "sg_num": 216, "roles": ["A", "B"],
                   "basis": [(0, 0, 0), (0.25, 0.25, 0.25)], "a": 5.4, "setting": 1},
    "fluorite":  {"sg": "Fm-3m", "sg_num": 225, "roles": ["A", "B"],
                  "basis": [(0, 0, 0), (0.25, 0.25, 0.25)], "a": 5.5, "setting": 1},
    "spinel":    {"sg": "Fd-3m", "sg_num": 227, "roles": ["A", "B", "X"],
                  "basis": [(0.125, 0.125, 0.125), (0.5, 0.5, 0.5), (0.2625, 0.2625, 0.2625)],
                  "a": 8.4, "setting": 2},
}


def _build_template(proto_name, role_elements, vol_per_atom):
    """Build the prototype cell with real elements via ASE's Wyckoff tables,
    scaled to the candidate's volume/atom."""
    p = PROTOTYPES[proto_name]
    species = [role_elements[r] for r in p["roles"]]
    at = _ase_crystal(species, basis=p["basis"], spacegroup=p["sg_num"], setting=p["setting"],
                      cellpar=[p["a"], p["a"], p["a"], 90, 90, 90])
    tmpl = AseAtomsAdaptor.get_structure(at)
    scale = (vol_per_atom * len(tmpl) / tmpl.volume) ** (1 / 3)
    tmpl.scale_lattice(tmpl.volume * scale ** 3)
    return tmpl


def _role_assignments(elements, proto_name):
    """Map candidate elements to prototype roles by electronegativity.

    Spinel with two elements (Co3O4): one cation fills both cation roles, so
    normal-vs-inverse is species-invisible by construction. Spinel with three
    elements: try both cation permutations, because inversion is real physics.
    """
    roles = PROTOTYPES[proto_name]["roles"]
    ordered = sorted(elements, key=lambda e: Element(e).X)
    if proto_name == "spinel":
        if len(ordered) == 2:
            return [{"A": ordered[0], "B": ordered[0], "X": ordered[1]}]
        if len(ordered) == 3:
            c1, c2, anion = ordered
            return [{"A": c1, "B": c2, "X": anion}, {"A": c2, "B": c1, "X": anion}]
        return []
    if len(ordered) != len(roles):
        return []
    return [dict(zip(roles, ordered))]


def _anonymize(s):
    s2 = s.copy()
    for i in range(len(s2)):
        s2.replace(i, "Si")
    return s2


def _mismatch_count(cand, tmpl, shift_grid=(0.0, 0.25, 0.5, 0.75)):
    """Optimal atom-to-atom assignment (Hungarian) between candidate and
    template in their shared cell, over a small grid of origin shifts.

    Detail metric only: the boolean verdicts come from StructureMatcher.
    Assumes same cell setting (true for as-built templates vs conventional
    CIFs of these cubic prototypes).
    """
    a_c = cand.lattice.matrix
    cf, tf = cand.frac_coords, tmpl.frac_coords
    best = None
    for sx in shift_grid:
        for sy in shift_grid:
            for sz in shift_grid:
                shifted = (tf + np.array([sx, sy, sz])) % 1.0
                df = cf[:, None, :] - shifted[None, :, :]
                df -= np.round(df)
                D = np.linalg.norm(df @ a_c, axis=2)
                r, c = linear_sum_assignment(D)
                cost = D[r, c].mean()
                if best is None or cost < best[0]:
                    mism = sum(1 for i, j in zip(r, c)
                               if cand[i].specie.symbol != tmpl[j].specie.symbol)
                    best = (cost, mism, list(zip(r.tolist(), c.tolist())))
    details = [{"candidate_site": i, "candidate_species": cand[i].specie.symbol,
                "template_site": j, "template_species": tmpl[j].specie.symbol}
               for i, j in best[2] if cand[i].specie.symbol != tmpl[j].specie.symbol]
    return best[1], round(float(best[0]), 4), details[:10]


def prototype_gate(structure, proto_name, stol=0.35):
    """Compare a candidate against a declared prototype, species-aware.

    PASS                geometry and species roles both match
    PASS(inverted)      geometry matches; species match only with cation roles
                        inverted (inverse spinel - real physics)
    INSPECT(species)    geometry matches but atoms sit on sites the template
                        labels for another species: label swap (corruption)
                        or genuine antisite disorder (physics)
    INSPECT(geometry)   point set does not match the declared prototype
    """
    if proto_name not in PROTOTYPES:
        return {"prototype": proto_name, "verdict": "ERROR", "detail": "unknown prototype"}
    elements = sorted({sp.symbol for sp in structure.species})
    assignments = _role_assignments(elements, proto_name)
    if not assignments:
        return {"prototype": proto_name, "verdict": "INSPECT(geometry)",
                "detail": f"candidate chemistry {elements} does not fit the {proto_name} roles"}

    vol_per_atom = structure.volume / len(structure)
    matcher = StructureMatcher(primitive_cell=False, attempt_supercell=True,
                               ltol=0.1, stol=stol, angle_tol=5)
    results = []
    for assign in assignments:
        tmpl = _build_template(proto_name, assign, vol_per_atom)
        geo_ok = matcher.fit(_anonymize(structure), _anonymize(tmpl))
        if not geo_ok:
            results.append({"assignment": assign, "verdict": "INSPECT(geometry)",
                            "detail": "point set does not match declared prototype"})
            continue
        species_ok = matcher.fit(structure, tmpl)
        if species_ok:
            results.append({"assignment": assign, "verdict": "PASS",
                            "mismatched_atoms": 0, "total_atoms": len(structure)})
            continue
        if len(tmpl) == len(structure):
            n_mism, mean_d, details = _mismatch_count(structure, tmpl)
        else:
            n_mism, mean_d, details = None, None, []
        results.append({"assignment": assign, "verdict": "INSPECT(species)",
                        "mismatched_atoms": n_mism, "total_atoms": len(structure),
                        "mean_assigned_distance_A": mean_d, "mismatches": details})

    results.sort(key=lambda r: (0 if r["verdict"] == "PASS" else 1,
                                r.get("mismatched_atoms") if r.get("mismatched_atoms") is not None else 10**9))
    best = results[0]
    out = {"prototype": proto_name, "verdict": best["verdict"],
           "best_assignment": best.get("assignment"),
           "mismatched_atoms": best.get("mismatched_atoms"), "total_atoms": best.get("total_atoms"),
           "all_assignments": results}
    if best["verdict"] == "INSPECT(species)":
        out["note"] = ("Geometry matches the declared prototype but species occupy each other's roles. "
                       "Either a label swap (file corruption) or genuine antisite disorder (real physics, "
                       "e.g. spinel inversion). Inspect, do not auto-reject.")
        out["mismatches"] = best.get("mismatches")
        if len(results) > 1 and any(r["verdict"] == "PASS" for r in results[1:]):
            out["verdict"] = "PASS(inverted roles)"
            out["note"] = ("Matches only with cation roles inverted relative to the electronegativity "
                           "default -> consistent with an inverse spinel, which is real physics.")
    return out