#!/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,
reference-structure displacement, and (optionally) agreement with a declared
prototype. It is deliberately not a stability, synthesis, or property
predictor.

v4.4 (calibrated against 88 wild COD CIFs, post 019fdf2b): the min-pair gate is
occupancy-aware (split/partial sites are alternatives in a probability
distribution, not simultaneous residents), the reference matcher standardizes
the input setting before comparing and bows out on triclinic cells, the
formula check compares reduced compositions against the raw parse (allowing
declared-but-unlocated H), and space-group labels are compared by IT number,
not by string (Pbnm == Pnma == #62).

The card catches what symmetry analysis alone cannot: coordinate permutation,
overlapping atoms, systematic geometric corruption that hides behind a
valid-looking space group label, and species label swaps that preserve every
distance. v4.3 also reads the file's claims about itself (declared space
group, symmetry operations, formula) and checks composition charge balance —
a missing anion is invisible to geometry gates but not to oxidation states. The reference-structure matcher compares each atom to its ideal
position in the symmetry-refined cell. The prototype gate compares the
candidate against a species-carrying template built from an independent
Wyckoff implementation (ASE spacegroup tables) — the one check that catches
label swaps, and the one that must never share the generator's bugs.

Usage:
    python structure_sanity_card.py <cif> --output card.md [--json report.json]
                                        [--prototype rocksalt|cesium_chloride|
                                          zincblende|fluorite|spinel]

Dependencies: pymatgen, numpy, spglib. The prototype gate additionally needs
scipy and ase (imported only when --prototype is used).

Provenance: built during the August 2026 curiosity-window series on structure
validation (Ouro, #materials-science). The scars that shaped it — including
the pymatgen origin-choice spinel trap — are documented in the accompanying
posts and in lessons-structure-validation.md.
"""
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, Element
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.

    v4.4: occupancy-aware. Partially occupied sites are alternatives in a
    probability distribution, not simultaneous residents. A short pair only
    counts as a physical overlap when the two sites can coexist in the same
    cell, i.e. occ_i + occ_j > 1. The wild-sample calibration (post 019fdf2b)
    showed 9 of 10 min-pair FAILs were disorder representations: split oxygen
    sites at 0.15 A with 50% occupancy each. The one genuine survivor (a
    partial W 0.33 A from a fully occupied O) still trips the gate because
    the occupancies sum above 1.

    Returns the shortest *simultaneously occupied* pair distance as
    "distance" (the gate verdict uses this), plus "distance_any" — the
    shortest pair ignoring occupancy — so the card can note excluded
    split-site alternatives instead of screaming about them.
    """
    if len(structure) < 2:
        return {"distance": None, "sites": None}
    occs = [float(sum(site.species.values())) for site in structure]
    labels = [site.species_string for site in structure]
    d_min, pair = float("inf"), None          # shortest coexisting pair
    d_any, pair_any = float("inf"), None      # shortest pair, occupancy-blind
    for i, j in combinations(range(len(structure)), 2):
        d = structure.get_distance(i, j)
        if d < d_any:
            d_any, pair_any = d, (i, j, labels[i], labels[j], round(occs[i], 3), round(occs[j], 3))
        if occs[i] + occs[j] > 1.0 + 1e-6 and d < d_min:
            d_min, pair = d, (i, j, labels[i], labels[j], round(occs[i], 3), round(occs[j], 3))
    out = {
        "distance": round(float(d_min), 4) if pair else None,
        "sites": pair,
        "distance_any": round(float(d_any), 4) if pair_any else None,
        "sites_any": pair_any,
    }
    if pair is None and pair_any is not None:
        out["all_mutually_exclusive"] = True
    elif pair_any and pair and d_any < d_min - 1e-9:
        out["excluded_alternatives"] = True
    return out


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.

    v4.4: the input is first mapped into the standard setting with
    spglib.standardize_cell (no idealization), so non-standard space group
    settings (Pcab vs Pbca, P 21/n vs P 21/c) compare in the same frame as
    the refined ideal. v4.2's anchor origin-shift search then absorbs any
    remaining origin-choice difference. Triclinic cells (IT #1-2) bow out:
    with no meaningful unique setting to refine against, the comparison only
    manufactures displacement storms (the Fe16Sb P-1 false positive).

    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.
    """
    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

    base = {"available": True, "spacegroup": sg_symbol, "number": sg_number}

    # v4.4: triclinic bow-out. P1 was already trivial (every atom is its own
    # Wyckoff position); P-1 has centrosymmetry but no unique axis setting,
    # so refine_cell's frame choice manufactures false displacement storms.
    if sg_number <= 2:
        return {**base, "is_p1": sg_number == 1,
                "uninformative": ("triclinic" if sg_number == 2 else "p1"),
                "max_displacement": 0.0, "mean_displacement": 0.0,
                "median_displacement": 0.0, "n_above_0p5": 0, "n_above_1p0": 0,
                "n_above_2p0": 0, "n_atoms": len(structure), "per_element": {}}

    # v4.4: normalize the input into the standard setting BEFORE comparing.
    # standardize_cell(no_idealize=True) keeps measured positions but rotates
    # non-standard settings (Pcab, P 21/n) into the conventional frame that
    # refine_cell also reports, so only an origin choice can still differ.
    std = spglib.standardize_cell(cell, to_primitive=False, no_idealize=True, symprec=symprec)
    if std is None:
        return {**base, "error": "standardize_cell failed"}
    lat_s, pos_s, num_s = std

    refined = spglib.refine_cell((lat_s, pos_s, num_s), symprec=symprec)
    if refined is None:
        return {**base, "error": "refine_cell failed"}

    lat_r, pos_r, num_r = refined
    if len(num_r) != len(num_s):
        return {**base, "error": f"refined cell has {len(num_r)} atoms, expected {len(num_s)}"}

    # v4.2 anchor shift search (origin choice), now in the shared standard
    # frame: find the rigid shift that best aligns the standardized input to
    # the refined positions (anchor: each same-species candidate for atom 0).
    pos_in = np.array(pos_s)
    best_shift = np.zeros(3)
    best_max = float("inf")
    anchor_z = int(num_s[0])
    for j in range(len(num_r)):
        if int(num_r[j]) != anchor_z:
            continue
        shift = (pos_r[j] - pos_in[0]) % 1.0
        pos_sh = (pos_in + shift) % 1.0
        mx = 0.0
        for i in range(len(pos_sh)):
            zi = int(num_s[i])
            dmin = float("inf")
            for k in range(len(num_r)):
                if int(num_r[k]) == zi:
                    d_frac = pos_sh[i] - pos_r[k]
                    d_frac = d_frac - np.round(d_frac)
                    d = float(np.linalg.norm(d_frac @ lat_s))
                    if d < dmin:
                        dmin = d
            if dmin > mx:
                mx = dmin
            if mx >= best_max:
                break
        if mx < best_max:
            best_max = mx
            best_shift = shift
    pos_in = (pos_in + best_shift) % 1.0

    displacements = []
    for i in range(len(num_s)):
        zi = int(num_s[i])
        best_d = float("inf")
        for j in range(len(num_r)):
            if int(num_r[j]) == zi:
                d_frac = pos_in[i] - pos_r[j]
                d_frac = d_frac - np.round(d_frac)
                d_cart = d_frac @ lat_s
                d = float(np.linalg.norm(d_cart))
                if d < best_d:
                    best_d = d
        displacements.append(best_d)

    disps = np.array(displacements)
    per_element = {}
    from pymatgen.core import Element as _El
    for i in range(len(num_s)):
        per_element.setdefault(_El.from_Z(int(num_s[i])).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 {
        **base,
        "is_p1": False,
        "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,
        "setting_shift": [round(float(x), 4) for x in best_shift],
    }


# ---------------------------------------------------------------------------
# v4.3 gates: the file's claims about itself, and composition plausibility
#
# The battery (post 019fdebb) showed three blind spots: a missing atom passes
# every geometry gate, the declared SG label is consumed by the parser but
# never checked, and species swaps need a declared prototype. The first two
# close here; the third is narrowed by matching species-strict FIRST in the
# prototype gate (anonymized point-set matching on dense cells absorbs real
# positional corruption — see LESSONS.md).
# ---------------------------------------------------------------------------

ANIONS = {"H", "B", "C", "N", "O", "F", "Si", "P", "S", "Cl", "As", "Se", "Te", "Br", "I"}

# Common oxidation states, SMACT-style. Plausibility only — not a chemistry claim.
COMMON_OXI = {
    "H": [1, -1], "Li": [1], "Na": [1], "K": [1], "Rb": [1], "Cs": [1],
    "Be": [2], "Mg": [2], "Ca": [2], "Sr": [2], "Ba": [2],
    "B": [3], "Al": [3], "Ga": [3], "In": [3], "Tl": [1, 3],
    "C": [4, -4], "Si": [4, -4], "Ge": [2, 4], "Sn": [2, 4], "Pb": [2, 4],
    "N": [-3, 3, 5], "P": [-3, 3, 5], "As": [-3, 3, 5], "Sb": [-3, 3, 5], "Bi": [3, 5],
    "O": [-2, -1], "S": [-2, 4, 6], "Se": [-2, 4, 6], "Te": [-2, 4, 6],
    "F": [-1], "Cl": [-1, 1, 3, 5, 7], "Br": [-1, 1, 5], "I": [-1, 1, 5, 7],
    "Sc": [3], "Ti": [2, 3, 4], "V": [2, 3, 4, 5], "Cr": [2, 3, 6],
    "Mn": [2, 3, 4, 7], "Fe": [2, 3], "Co": [2, 3], "Ni": [2, 3],
    "Cu": [1, 2], "Zn": [2], "Y": [3], "Zr": [4], "Nb": [3, 5], "Mo": [4, 6],
    "Ru": [3, 4], "Rh": [3], "Pd": [2, 4], "Ag": [1], "Cd": [2],
    "Hf": [4], "Ta": [5], "W": [4, 6], "Re": [4, 7], "Os": [4], "Ir": [3, 4],
    "Pt": [2, 4], "Au": [1, 3], "Hg": [1, 2],
    "La": [3], "Ce": [3, 4], "Pr": [3, 4], "Nd": [3], "Sm": [2, 3], "Eu": [2, 3],
    "Gd": [3], "Tb": [3, 4], "Dy": [3], "Ho": [3], "Er": [3], "Tm": [3],
    "Yb": [2, 3], "Lu": [3], "Th": [4], "U": [4, 6],
}


def _norm_hm(raw: str | None) -> str | None:
    """'F m -3 m' / "'Fm-3m'" -> 'Fm-3m'."""
    if raw is None:
        return None
    s = raw.strip().strip("'").strip('"')
    return "".join(s.split())


# Non-standard / slashed Hermann-Mauguin settings pymatgen's SpaceGroup does
# not parse. Keys are space-stripped, uppercased. Values are IT numbers — the
# whole point of v4.4's label check is that Pbnm and Pnma are both #62.
_HM_IT_ALIASES = {
    "P21/N": 14, "P21/C": 14, "P21/A": 14, "P21/B": 14, "P121/N1": 14,
    "P121/C1": 14, "P121/A1": 14, "P121/B1": 14, "P21/11": 14,
    "C2/N": 15, "I2/A": 15, "I2/C": 15, "C12/C1": 15, "C12/N1": 15,
    "C2/M": 12, "I2/M": 12, "P2/M": 10, "P2/C": 13, "P2/N": 13, "P2/A": 13,
    "PNNA": 52, "PBNN": 52, "PCCN": 49, "PBAA": 54, "PMAA": 49,
    "AMAM": 63, "BMBM": 63, "CM CM": 63, "CMCM": 63, "AMMA": 67,
    "FDD2": 43, "I41/AMD": 141, "I41/ACD": 142, "R3C": 161, "R3M": 160,
}


def _hm_to_it(hm: str | None) -> int | None:
    """Map a Hermann-Mauguin label to its IT number, surviving non-standard
    axis settings. Returns None when the label cannot be resolved."""
    if not hm:
        return None
    key = "".join(str(hm).split()).upper()
    try:
        from pymatgen.symmetry.groups import SpaceGroup
        return int(SpaceGroup(key).int_number)
    except Exception:
        pass
    return _HM_IT_ALIASES.get(key)


def declared_metadata(cif_path: Path) -> dict:
    """Read what the CIF declares about itself, raw from the text.

    pymatgen's parser consumes these tags but never checks them against what
    it built. Returns declared H-M label, IT number, ops-loop length,
    formula_sum, and the number of listed atom sites.
    """
    import re
    text = Path(cif_path).read_text()
    lines = text.splitlines()

    def tag(name):
        for ln in lines:
            if ln.strip().startswith(name):
                rest = ln.strip()[len(name):].strip()
                return rest if rest else None
        return None

    hm = tag("_symmetry_space_group_name_H-M") or tag("_space_group_name_H-M_alt")
    it_raw = tag("_symmetry_Int_Tables_number") or tag("_space_group_IT_number")
    try:
        it_num = int(str(it_raw).strip("'").strip('"')) if it_raw else None
    except ValueError:
        it_num = None
    formula_sum = tag("_chemical_formula_sum")

    # count ops in the symmetry-equiv loop and sites in the atom_site loop
    ops_n, sites_n = None, None
    i = 0
    while i < len(lines):
        if lines[i].strip() == "loop_":
            j = i + 1
            headers = []
            while j < len(lines) and lines[j].strip().startswith("_"):
                headers.append(lines[j].strip().split()[0])
                j += 1
            body = []
            while j < len(lines) and lines[j].strip() and not lines[j].strip().startswith(("loop_", "_", "data_")):
                body.append(lines[j])
                j += 1
            if any(h in ("_symmetry_equiv_pos_as_xyz", "_space_group_symop_operation_xyz") for h in headers):
                ops_n = len(body)
            if any(h.startswith("_atom_site_fract") for h in headers):
                sites_n = len(body)
            i = j
        else:
            i += 1
    return {"declared_hm": _norm_hm(hm), "declared_it": it_num,
            "ops_listed": ops_n, "formula_sum": formula_sum,
            "sites_listed": sites_n}


def metadata_check(structure: Structure, cif_path: Path, sweep: list[dict]) -> dict:
    """Consistency of the file's declared metadata with what was parsed."""
    md = declared_metadata(cif_path)
    out = dict(md)
    out["label_ops_consistent"] = None
    out["formula_sum_matches"] = None
    out["declared_vs_detected"] = None
    out["redundant_listing"] = False

    # 1. declared label vs ops-loop length (centering conventions accepted)
    if md["declared_hm"] and md["ops_listed"] is not None:
        try:
            from pymatgen.symmetry.groups import SpaceGroup
            n_full = len(SpaceGroup(md["declared_hm"]).symmetry_ops)
            out["ops_expected"] = n_full
            out["label_ops_consistent"] = md["ops_listed"] in {n_full, max(n_full // 2, 1), max(n_full // 4, 1)}
        except Exception:
            out["ops_expected"] = None

    # 2. declared formula vs parsed composition (v4.4)
    # Compare REDUCED compositions, so formula-unit vs full-cell-contents
    # declarations both match. Compare against the raw parse (occupancies
    # included), not the ordered surrogate — the surrogate deletes minority
    # species from mixed sites and manufactures mismatches. Hydrogen declared
    # in the formula but never given coordinates is honest standard practice
    # (H not located in the refinement), not a stale header.
    if md["formula_sum"]:
        try:
            from pymatgen.core import Composition
            decl = Composition(md["formula_sum"].strip("'").strip('"')).get_el_amt_dict()
            parsed = structure.composition.get_el_amt_dict()
            h_unlocated = False
            if "H" in decl and "H" not in parsed:
                decl = {k: v for k, v in decl.items() if k != "H"}
                h_unlocated = True
            # Scale-invariant comparison: mole fractions, not absolute
            # amounts. (reduced_composition cannot reduce float compositions
            # with rounding noise — Ba11.98Pr4.02Co8O30 stays "reduced" and
            # false-mismatches Ba6Pr2Co4O15. Scar from the v4.4 re-run.)
            def _fracs(d):
                tot = sum(d.values())
                return {el: amt / tot for el, amt in d.items()} if tot > 0 else {}
            rd, rp = _fracs(decl), _fracs(parsed)
            ok = bool(rd) and all(
                abs(rd.get(el, 0.0) - rp.get(el, 0.0)) <= max(0.03 * max(rd.get(el, 0.0), rp.get(el, 0.0)), 0.005)
                for el in set(rd) | set(rp))
            out["formula_sum_matches"] = ok
            if h_unlocated:
                out["formula_h_unlocated"] = True
        except Exception:
            pass

    # 3. declared SG vs detected, by IT number (v4.4). Pbnm and Pnma are the
    # same group (#62) in different axis settings; a string comparison flags
    # honest setting choices as corruption. Compare the declared IT number
    # (from the header tag, else resolved from the label) against every number
    # the symmetry sweep detects.
    valid = [e for e in sweep if e["spacegroup"] != "ERROR"]
    if md["declared_hm"] and valid:
        tight, loose = valid[0]["spacegroup"], valid[-1]["spacegroup"]
        decl_it = md["declared_it"] or _hm_to_it(md["declared_hm"])
        out["declared_it_resolved"] = decl_it
        detected_numbers = {e["number"] for e in valid if e.get("number", -1) > 0}
        if md["declared_hm"] in ("P1", "P 1") or decl_it == 1:
            out["declared_vs_detected"] = "p1_declared"
        elif decl_it is not None:
            out["declared_vs_detected"] = "match" if decl_it in detected_numbers else "mismatch"
        elif md["declared_hm"] in (tight, loose):
            out["declared_vs_detected"] = "match"   # unresolvable label: string fallback
        else:
            out["declared_vs_detected"] = "mismatch"
        out["detected_tight"], out["detected_loose"] = tight, loose

    # 4. redundant full-cell listing under non-P1 ops (parser-fragile style)
    if md["sites_listed"] is not None and md["ops_listed"] and md["ops_listed"] > 1:
        out["redundant_listing"] = md["sites_listed"] >= len(structure)
    return out


def charge_balance(el_amt: dict) -> tuple[bool | None, dict | None]:
    """SMACT-style plausibility: does any common oxidation-state assignment
    balance charge? Returns (valid, example_assignment). (None, None) when an
    element is outside the table."""
    from fractions import Fraction
    counts = {}
    for el, amt in el_amt.items():
        if el not in COMMON_OXI:
            return None, None
        counts[el] = Fraction(amt).limit_denominator(24)
    scale = 1
    for f in counts.values():
        scale = scale * f.denominator // __import__("math").gcd(scale, f.denominator)
    int_counts = {el: int(f * scale) for el, f in counts.items()}
    # DP over achievable total charges
    reach = {0: {}}
    for el, n in int_counts.items():
        new = {}
        for total, assign in reach.items():
            for ox in COMMON_OXI[el]:
                t = total + n * ox
                if t not in new:
                    new[t] = {**assign, el: ox}
        reach = new
        if len(reach) > 20000:  # safety valve
            return None, None
    if 0 in reach:
        return True, reach[0]
    return False, None


def stoichiometry_check(structure: Structure) -> dict:
    """Composition plausibility via oxidation-state charge balance.

    v4.4: skipped on disordered cells. Per-cell charge balance assumes one
    representative configuration; fractional occupancies violate that (and
    disorder often exists precisely to balance charge non-locally). The
    occupancy gate already flags disorder for inspection."""
    el_amt = {str(el): amt for el, amt in structure.composition.get_el_amt_dict().items()}
    out = {"composition": structure.composition.formula.replace(" ", "")}
    if not structure.is_ordered:
        out.update({"skipped": True, "charge_balance": None,
                    "reason": "disordered composition; per-cell charge balance is ill-posed "
                              "for fractional occupancies"})
        return out
    if not any(el in ANIONS for el in el_amt):
        out.update({"skipped": True, "charge_balance": None,
                    "reason": "no anion species (intermetallic/elemental); charge balance not applicable"})
        return out
    valid, assignment = charge_balance(el_amt)
    out.update({"skipped": False, "charge_balance": valid, "assignment": assignment})
    if valid is None:
        out["reason"] = "element outside the common-oxidation-state table"
    return out


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

    # 1. Minimum pair distance — occupancy-aware (v4.4). Split/partial sites
    # are alternatives in a disorder model, not simultaneous residents; only
    # pairs whose occupancies sum above 1 can physically overlap.
    mp = report["min_pair"]
    mpd, sites = mp["distance"], mp["sites"]
    d_any, sites_any = mp.get("distance_any"), mp.get("sites_any")
    if mpd is not None and mpd < 0.5:
        notes.append(("FAIL", f"Minimum pair distance {mpd:.3f} Å between simultaneously "
                      f"occupied sites {sites} is physically impossible. Atoms are "
                      "overlapping. Structure is corrupted."))
    elif mpd is not None and mpd < 1.0:
        notes.append(("CHECK", f"Minimum pair distance {mpd:.3f} Å between simultaneously "
                      f"occupied sites {sites} is suspiciously short. Atoms 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."))
    if d_any is not None and d_any < 1.5 and (mpd is None or d_any < mpd - 1e-9):
        notes.append(("NOTE", f"Minimum pair distance {d_any:.3f} Å between sites {sites_any} "
                      "is excluded from the verdict: their occupancies sum to 1.0 or less, "
                      "so they are alternatives in a disorder model, not simultaneous "
                      "residents. Confirm the disorder representation is intended."))
    if mpd is None and d_any is not None:
        notes.append(("NOTE", "Minimum pair distance: every listed pair is mutually exclusive "
                      "(occupancies sum to 1.0 or less), so no simultaneously-occupied pair "
                      "exists to judge. This is a fully disorder-modelled cell."))

    # 1b. Declared metadata — what the file claims about itself (v4.3)
    md = report.get("metadata", {})
    if md:
        decl = md.get("declared_hm")
        if md.get("label_ops_consistent") is False:
            notes.append(("CHECK", f"Declared {decl} (#{md.get('declared_it')}) but the symmetry "
                          f"operations loop lists {md.get('ops_listed')} operations "
                          f"(expected ~{md.get('ops_expected')}). The label and the symmetry "
                          "data disagree — a stale or mis-stamped header."))
        if md.get("formula_sum_matches") is False:
            notes.append(("CHECK", f"Declared _chemical_formula_sum {md.get('formula_sum')} does "
                          f"not match the listed sites ({report['formula']}). "
                          "Stale header, or a site was added/dropped without updating the metadata."))
        elif md.get("formula_h_unlocated") and md.get("formula_sum_matches"):
            notes.append(("NOTE", f"Declared formula {md.get('formula_sum')} includes H atoms "
                          "with no coordinates — H not located in the refinement, which is "
                          "standard practice. Counted as a match."))
        dvt = md.get("declared_vs_detected")
        if dvt == "mismatch":
            notes.append(("CHECK", f"Declared space group {decl} is not recovered from the geometry "
                          f"(sweep detects {md.get('detected_tight')} to {md.get('detected_loose')}). "
                          "The label describes a different structure than the coordinates do."))
        elif dvt == "match" and md.get("label_ops_consistent") is not False:
            notes.append(("PASS", f"Declared space group {decl} is consistent with the "
                          "detected symmetry and the symmetry data."))
        elif dvt == "p1_declared" and md.get("detected_loose") not in (None, "P1"):
            notes.append(("NOTE", f"No symmetry declared (P1-style full-cell listing); the geometry "
                          f"itself has {md.get('detected_loose')}. Fine for pipelines, but the "
                          "file carries less information than it could."))
        if md.get("redundant_listing"):
            notes.append(("NOTE", f"All {md.get('sites_listed')} cell sites are listed explicitly under "
                          f"{md.get('ops_listed')} symmetry operations. Parsers that do not merge "
                          "equivalent sites will build a pile of overlapping atoms from this file."))

    # 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."))

    # 3b. Stoichiometry plausibility (v4.3)
    st = report.get("stoichiometry", {})
    if st:
        if st.get("skipped"):
            notes.append(("NOTE", "Oxidation-state plausibility skipped: "
                          + st.get("reason", "not applicable") + "."))
        elif st.get("charge_balance") is True:
            assign = ", ".join(f"{el}{ox:+d}" for el, ox in sorted((st.get("assignment") or {}).items()))
            notes.append(("PASS", f"Charge balances with common oxidation states ({assign})."))
        elif st.get("charge_balance") is False:
            notes.append(("CHECK", f"No common oxidation-state assignment balances charge for "
                          f"{st.get('composition')}. A missing or extra site reads exactly like "
                          "this — or the chemistry is genuinely unusual (mixed valence, "
                          "intercalation). Inspect before trusting."))
        else:
            notes.append(("NOTE", "Oxidation-state plausibility inconclusive: "
                          + st.get("reason", "unknown") + "."))

    # 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("uninformative") == "triclinic":
            notes.append(("NOTE", "Reference-structure match bows out: the cell is triclinic "
                          f"({rm['spacegroup']}), where no unique standard setting exists to "
                          "refine against. Displacement numbers here would say more about the "
                          "matcher than about the structure — rely on the geometry gates."))
        elif rm.get("uninformative") == "p1" or rm.get("is_p1"):
            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."))


    # 7. Prototype gate (only when a prototype was declared)
    pg = report.get("prototype_gate")
    if pg:
        verdict = pg.get("verdict", "ERROR")
        proto = pg.get("prototype", "?")
        if verdict == "PASS":
            notes.append(("PASS", f"Prototype gate: geometry and species roles both match "
                          f"declared {proto}."))
        elif verdict == "PASS(inverted roles)":
            notes.append(("NOTE", f"Prototype gate ({proto}): {pg.get('note', '')}"))
        elif verdict == "INSPECT(species)":
            n = pg.get("mismatched_atoms")
            tot = pg.get("total_atoms")
            frac = f" ({n}/{tot} sites)" if n is not None and tot else ""
            notes.append(("CHECK", f"Prototype gate ({proto}): geometry matches but species "
                          f"occupy each other's roles{frac}. Label swap (corruption) or genuine "
                          f"antisite disorder (physics) — inspect, do not auto-reject."))
        elif verdict == "INSPECT(geometry)":
            extra = (f" Mean candidate-template distance {pg.get('mean_assigned_distance_A')} Å."
                     if pg.get("mean_assigned_distance_A") is not None else "")
            notes.append(("CHECK", f"Prototype gate: point set does not match declared {proto}.{extra} "
                          "Wrong prototype, wrong setting, free-parameter corruption, or severe distortion."))
        elif verdict == "INSPECT(chemistry)":
            notes.append(("CHECK", f"Prototype gate ({proto}): {pg.get('detail', '')}"))
        else:
            notes.append(("NOTE", f"Prototype gate ({proto}) unavailable: "
                          f"{pg.get('detail', 'unknown error')}"))

    return notes


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 (no coexisting pair)"
    if mpd.get("distance_any") is not None and (
            mpd.get("excluded_alternatives") or mpd.get("all_mutually_exclusive")):
        mpd_str += f"; shortest raw pair {mpd['distance_any']:.4f} Å excluded (mutually exclusive occupancies)"

    rm0 = report.get("reference_match", {})
    shift = rm0.get("setting_shift") if rm0.get("available") else None
    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),
    ]

    md = report.get("metadata", {})
    if md.get("declared_hm"):
        rows.append(("Declared SG (file)", f"{md['declared_hm']} (#{md.get('declared_it')}), "
                     f"{md.get('ops_listed')} ops listed"))
    if md.get("formula_sum"):
        rows.append(("Declared formula (file)", str(md["formula_sum"])))
    if shift and max(abs(x) for x in shift) > 0.01:
        rows.append(("Origin choice", "non-standard vs spglib; matched after shift "
                     + f"({shift[0]:.3f}, {shift[1]:.3f}, {shift[2]:.3f})"))
    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") and rm.get("uninformative"):
        rm_section = (f"\n## Reference-structure match\n\nBow-out: the cell is "
                      f"{'triclinic (' + rm['spacegroup'] + ')' if rm['uninformative'] == 'triclinic' else 'P1'} — "
                      "no unique standard setting to refine against, so the displacement "
                      "comparison is uninformative here.\n")
    elif 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"

    # Prototype gate section
    pg = report.get("prototype_gate")
    pg_section = ""
    if pg:
        pg_rows = [
            ("Declared prototype", str(pg.get("prototype"))),
            ("Verdict", str(pg.get("verdict"))),
            ("Best role assignment", str(pg.get("best_assignment"))),
            ("Mismatched atoms", f"{pg.get('mismatched_atoms')} / {pg.get('total_atoms')}"),
        ]
        pg_table = "\n".join(f"| {k} | {v} |" for k, v in pg_rows)
        mism = pg.get("mismatches") or []
        mism_lines = "\n".join(
            f"| {m['candidate_site']} | {m['candidate_species']} | {m['template_site']} | {m['template_species']} |"
            for m in mism
        )
        mism_table = (f"\n| Candidate site | Candidate species | Template site | Template species |\n"
                      f"| --- | --- | --- | --- |\n{mism_lines}\n") if mism else ""
        note_line = f"\n{pg['note']}\n" if pg.get("note") else ""
        pg_section = (
            "\n## Prototype gate\n\n"
            "Compares the candidate against a species-carrying template built from ASE's\n"
            "independent Wyckoff tables, scaled to the candidate's volume/atom. Species-strict\n"
            "matching runs first (v4.3); anonymized geometry only discriminates label swaps\n"
            "from parameter corruption. A species mismatch is always INSPECT, never BROKEN —\n"
            "antisite disorder is real physics.\n\n"
            "| Field | Value |\n| --- | --- |\n"
            f"{pg_table}\n{note_line}{mism_table}"
        )

    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}{pg_section}
## Bond length statistics

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

def generate_report(cif_path: Path, symprecs: list[float] | None = None,
                    prototype: str | None = None) -> dict:
    if symprecs is None:
        symprecs = DEFAULT_SYMPRECS
    structure_raw = Structure.from_file(str(cif_path))
    # v4.2: disordered sites crash .specie-based gates. Geometry gates card the
    # dominant-species surrogate; the occupancy gate reports disorder
    # separately (flag captured BEFORE the surrogate swap).
    # v4.4: gates that compare the file against itself — min-pair distance
    # (now occupancy-aware), declared-formula, and charge balance — run on the
    # RAW parse. The surrogate deletes minority species from mixed sites and
    # manufactured false formula/pair alarms in the wild calibration.
    _ordered_input = all(st.is_ordered for st in structure_raw)
    structure = structure_raw
    if not _ordered_input:
        geom = structure_raw.copy()
        for i, st in enumerate(geom):
            if not st.is_ordered:
                dominant = max(st.species.items(), key=lambda kv: kv[1])[0]
                geom[i] = dominant, st.frac_coords
        structure = geom
    report = {
        "source": str(cif_path),
        "formula": structure_raw.composition.reduced_formula,
        "sites": len(structure_raw),
        "volume_per_atom": round(float(structure_raw.volume / len(structure_raw)), 4),
        "density": round(float(structure_raw.density), 4),
        "ordered": _ordered_input,
        "lattice": {
            "a": round(float(structure_raw.lattice.a), 4),
            "b": round(float(structure_raw.lattice.b), 4),
            "c": round(float(structure_raw.lattice.c), 4),
            "alpha": round(float(structure_raw.lattice.alpha), 2),
            "beta": round(float(structure_raw.lattice.beta), 2),
            "gamma": round(float(structure_raw.lattice.gamma), 2),
        },
        "min_pair": min_pair_distance(structure_raw),
        "symmetry_sweep": fine_symmetry_sweep(structure, symprecs),
        "coordination": coordination_fingerprint(structure),
        "bond_stats": bond_stats(structure),
        "reference_match": reference_match(structure),
    }
    report["metadata"] = metadata_check(structure_raw, cif_path, report["symmetry_sweep"])
    report["stoichiometry"] = stoichiometry_check(structure_raw)
    if prototype is not None:
        report["prototype_gate"] = prototype_gate(structure, prototype)
    return report

# Prototype-aware gate (added v4, 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.
# ---------------------------------------------------------------------------


# 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."""
    from ase.spacegroup import crystal as _ase_crystal
    from pymatgen.io.ase import AseAtomsAdaptor
    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).
    """
    from scipy.optimize import linear_sum_assignment
    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(chemistry)",
                "detail": f"candidate species {elements} cannot fill the roles of {proto_name}. "
                          "The declaration may be wrong — or a species label may be corrupted."}

    from pymatgen.analysis.structure_matcher import StructureMatcher
    vol_per_atom = structure.volume / len(structure)
    matcher = StructureMatcher(primitive_cell=False, attempt_supercell=True,
                               ltol=0.1, stol=stol, angle_tol=5)
    # v4.3: species-strict FIRST. Anonymized union point-set matching on dense
    # cells absorbs 0.2-0.5 Å of positional corruption (LESSONS.md), so a
    # geometry-first order can certify a parameter-corrupted file. When
    # species-strict fails but anonymized geometry passes, the Hungarian mean
    # assigned distance decides what we are looking at: ~0 Å means a clean
    # label swap (INSPECT(species)); large means free-parameter corruption or
    # distortion that the anonymized match hid (INSPECT(geometry)).
    results = []
    for assign in assignments:
        tmpl = _build_template(proto_name, assign, vol_per_atom)
        species_ok = matcher.fit(structure, tmpl)
        if species_ok:
            results.append({"assignment": assign, "verdict": "PASS",
                            "mismatched_atoms": 0, "total_atoms": len(structure)})
            continue
        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
        if len(tmpl) == len(structure):
            n_mism, mean_d, details = _mismatch_count(structure, tmpl)
        else:
            n_mism, mean_d, details = None, None, []
        if mean_d is not None and mean_d > 0.1:
            results.append({"assignment": assign, "verdict": "INSPECT(geometry)",
                            "mean_assigned_distance_A": mean_d,
                            "detail": f"species-strict match fails and atoms sit on average "
                                      f"{mean_d:.2f} Å from template positions. Anonymized "
                                      "geometry matching absorbed this — the pattern is "
                                      "free-parameter corruption or real distortion, not a "
                                      "clean label swap."})
            continue
        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

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)
    parser.add_argument("--prototype", choices=sorted(PROTOTYPES), default=None,
                        help="Declare the expected prototype to enable the species-aware "
                             "template gate (catches label swaps geometry cannot see).")
    args = parser.parse_args()

    report = generate_report(args.cif, args.symprec, prototype=args.prototype)
    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()