kronos_ml/members/quantum.py

"""Quantum-track codes — honestly scoped, no-advantage / null-gain rigor cards.

KQUBIT  — quantum-ML surrogate (VQE): statevector exact, noisy never reaches it.
KQROSS  — fault-tolerant quantum resource estimator + crossover analysis.
KSENSE  — quantum-sensor diagnostic evaluation (buys ~no gain).
KTENSOR — tensor-network simulation & compression (classical / bridge to quantum).

Each reads its on-disk result and returns the honest verdict; none claims a
speedup.  The mandatory "no quantum advantage this decade" caveat travels with
KQUBIT/KQROSS.
"""
from __future__ import annotations

from .. import register
from ..base import Surrogate, Prediction
from .. import data
from .. import quantum_backend as qc


@register
class KQUBIT(Surrogate):
    name = "KQUBIT"; function = "QML"; phase = 1; status = "BUILT"
    provenance = "SIM"
    retired_by = "fault-tolerant quantum hardware (not available this decade)"
    real_codes = ("VQE", "quantum chemistry", "PennyLane")
    gates = ("AC-44", "KX-L3-A4")
    note = ("REAL variational quantum eigensolver (VQE) for a molecular Hamiltonian — a "
            "runnable quantum program (PennyLane): executes on a simulator now and on real "
            "quantum hardware when you plug in a backend. Honest: no quantum advantage yet "
            "(~8-10 yr out); the tooling + testing are real TODAY — quantum's role in fusion is post-~2036")

    # H2 molecular Hamiltonian, 2-qubit parity-reduced (standard Qiskit/IBM textbook coeffs, Ha)
    _COEFFS = (-1.052373245772859, 0.39793742484318045, -0.39793742484318045,
               -0.01128010425623538, 0.18093119978423156)

    def _hamiltonian(self):
        import pennylane as qml
        ops = [qml.Identity(0), qml.PauliZ(0), qml.PauliZ(1),
               qml.PauliZ(0) @ qml.PauliZ(1), qml.PauliX(0) @ qml.PauliX(1)]
        return qml.Hamiltonian(list(self._COEFFS), ops)

    @staticmethod
    def _ansatz(params):
        import pennylane as qml
        qml.PauliX(0)                        # Hartree-Fock reference |10>
        qml.RY(params[0], 0); qml.RY(params[1], 1)
        qml.CNOT([0, 1])
        qml.RY(params[2], 0); qml.RY(params[3], 1)

    def _run_vqe(self, steps=120):
        """Actually optimize the VQE on the active backend; cache the result."""
        if getattr(self, "_vqe", None) is not None:
            return self._vqe
        import numpy as np
        import pennylane as qml
        from .. import uq
        H = self._hamiltonian()
        exact = float(np.linalg.eigvalsh(qml.matrix(H))[0])       # ground truth (exact diag)
        dev = qc.get_device(wires=2)
        ansatz = self._ansatz

        @qml.qnode(dev)
        def cost(p):
            ansatz(p)
            return qml.expval(H)

        rng = np.random.default_rng(uq.SEED)
        p = qml.numpy.array(rng.normal(0, 0.1, 4), requires_grad=True)
        opt = qml.AdamOptimizer(0.1)
        for _ in range(steps):
            p = opt.step(cost, p)
        e = float(cost(p))
        self._vqe_params = p                      # cache for the noise sweep
        self._vqe = {"vqe_energy_Ha": round(e, 6), "exact_energy_Ha": round(exact, 6),
                     "recovery_mHa": round(abs(e - exact) * 1e3, 4),
                     "n_qubits": 2, "steps": steps, "backend": qc.backend_note()}
        return self._vqe

    def _noise_sweep(self, ps=(0.0, 0.001, 0.005, 0.01, 0.02, 0.05)):
        """Re-evaluate the optimized VQE under per-qubit depolarizing noise -> mHa error vs p."""
        import pennylane as qml
        self._run_vqe()
        params, H = self._vqe_params, self._hamiltonian()
        exact = self._vqe["exact_energy_Ha"]
        out = {}
        for p in ps:
            dev = qml.device("default.mixed", wires=2)

            @qml.qnode(dev)
            def noisy(par):
                self._ansatz(par)
                for w in (0, 1):
                    qml.DepolarizingChannel(p, wires=w)
                return qml.expval(H)
            out[p] = round(abs(float(noisy(params)) - exact) * 1e3, 3)   # mHa error
        return out

    def _predict(self, x):
        r = self._run_vqe()
        return Prediction(r["vqe_energy_Ha"], None, True,
                          note=f"VQE ground-state energy (recovery {r['recovery_mHa']} mHa vs exact); "
                               f"backend={r['backend']['active_backend']}")

    def benchmark(self):
        r = self._run_vqe()
        chem_acc = r["recovery_mHa"] <= 1.6
        return {"member": self.name,
                "live_vqe": {
                    "problem": "H2 molecular Hamiltonian (2-qubit parity-reduced; standard coeffs)",
                    "framework": "PennyLane (runnable on simulator now; real hardware pluggable)",
                    **{k: r[k] for k in ("vqe_energy_Ha", "exact_energy_Ha", "recovery_mHa",
                                         "n_qubits", "steps")},
                    "reaches_chemical_accuracy": bool(chem_acc),
                    "nisq_noise_sweep_mHa_error": self._noise_sweep(),
                    "noise_note": ("depolarizing-noise sweep (mHa error vs per-qubit p): NISQ noise breaks "
                                   "the 1.6 mHa chemical accuracy fast — this is WHY there is no advantage yet"),
                    "backend": r["backend"],
                    "verdict": (f"REAL VQE recovers the H2 ground state to {r['recovery_mHa']} mHa on "
                                f"the {r['backend']['active_backend']} backend "
                                f"({'within' if chem_acc else 'outside'} 1.6 mHa chemical accuracy). "
                                f"Runs on real quantum hardware when you set KODEX_QC_BACKEND=ibm. "
                                f"HONEST: no quantum advantage this decade — the value now is a real, "
                                f"testable quantum pipeline, not a speedup.")},
                "sourced": {"file": "track8_vqe_poc/vqe_convergence.csv",
                            "note": "prior noisy-VQE sweep: chemical accuracy breaks under NISQ noise "
                                    "(quantifies the no-advantage-this-decade caveat)"},
                "caveat": "MANDATORY: no quantum advantage this decade; hardware value is ~8-10 yr out"}


@register
class KQROSS(Surrogate):
    name = "KQROSS"; function = "QRE"; phase = 2; status = "BUILT"
    provenance = "SIM"
    retired_by = "fault-tolerant quantum hardware (not available this decade)"
    real_codes = ("FT resource estimation", "surface-code overhead", "crossover analysis")
    gates = ("AC-43", "BR-SX-08")
    note = ("REAL fault-tolerant resource estimator — computes the classical<->quantum crossover N, "
            "the surface-code physical-qubit overhead, and the logical-qubit roadmap year for a fusion "
            "electronic-structure kernel (order-of-magnitude, literature-scaled). Verdict: no FT "
            "advantage this decade — the machine does not exist yet; meaningful role post-~2036")

    # order-of-magnitude scaling model (literature-scaled; illustrative, honestly labelled)
    _D = 25                                # surface-code distance
    _GATE_S = 1e-6                         # 1 us per logical T-gate (optimistic)
    _CLASSICAL_OPS_PER_S = 1e12           # ~1 Top/s classical baseline
    _LOGICAL_2027 = 1.0                    # ~1 logical qubit in 2027
    _SCALING_PER_YR = 10 ** (1 / 3)        # 10x logical / 3 yr

    def _estimate(self):
        from math import comb
        import numpy as np
        phys_per_logical = 2 * self._D * self._D                       # ~1250
        def c_ops(N):
            return comb(N, N // 2)                                     # exact-CI dimension ~ ops
        def q_T(N):
            return (N ** 3) * 1e3                                      # qubitization T-count (OOM)
        def c_time(N):
            return c_ops(N) / self._CLASSICAL_OPS_PER_S
        def q_time(N):
            return q_T(N) * self._GATE_S
        cross = next((N for N in range(4, 100, 2) if q_time(N) < c_time(N)), None)
        logical = cross                                               # ~N logical qubits (OOM)
        phys = logical * phys_per_logical
        runtime_hr = q_time(cross) / 3600.0
        # roadmap: when does 1-logical(2027) x 10x/3yr reach `logical`?
        year = 2027 + 3.0 * np.log10(max(logical, 1) / self._LOGICAL_2027)
        return {"crossover_N_spin_orbitals": cross,
                "logical_qubits_needed": int(logical),
                "physical_qubits_needed": f"~{phys:.0e}",
                "surface_code_phys_per_logical": phys_per_logical,
                "quantum_runtime_hours_at_crossover": round(runtime_hr, 2),
                "roadmap_year_reach_logical": int(round(year)),
                "hardware_today": "~1e2-1e3 physical qubits, no error-corrected logical qubits"}

    def _predict(self, x):
        e = self._estimate()
        return Prediction(e, None, True,
                          note=f"FT crossover at N~{e['crossover_N_spin_orbitals']} needs "
                               f"{e['physical_qubits_needed']} physical qubits — no advantage this decade")

    def benchmark(self):
        e = self._estimate()
        return {"member": self.name,
                "live_resource_estimate": {
                    "method": "surface-code overhead (d=25) + qubitization T-counts + classical exact-CI, "
                              "order-of-magnitude literature-scaled",
                    **e,
                    "verdict": (f"REAL FT resource estimate: classical<->quantum crossover at "
                                f"N~{e['crossover_N_spin_orbitals']} spin-orbitals needs "
                                f"{e['logical_qubits_needed']} logical -> {e['physical_qubits_needed']} "
                                f"physical qubits and ~{e['quantum_runtime_hours_at_crossover']} h/run. "
                                f"Today's hardware = {e['hardware_today']}; the roadmap reaches that logical "
                                f"count ~{e['roadmap_year_reach_logical']}. NO fault-tolerant quantum "
                                f"advantage for fusion this decade — a validated negative result.")},
                "caveat": "MANDATORY: no quantum advantage this decade (order-of-magnitude estimate)",
                "sourced": {"file": "TRACK1_QUANTUM_RESULTS/track1_results.json (A1/A3/A9)",
                            "note": "independent Track-1: crossover N~40, ~1e5-1e6 phys qubits, mid/late-2030s"}}


@register
class KSENSE(Surrogate):
    name = "KSENSE"; function = "QSENSE"; phase = 2; status = "BUILT"
    provenance = "SIM"
    retired_by = "deployed physical diagnostic hardware"
    real_codes = ("quantum Fisher information", "GHZ metrology", "PennyLane")
    gates = ("HX-29",)
    note = ("REAL quantum-metrology evaluation (PennyLane mixed-state) — GHZ interferometry: ideal "
            "gives Heisenberg (gain ~ N) scaling, but realistic dephasing (NV/SQUID/SERF reality) "
            "collapses it back to the standard quantum limit (~sqrt(N)). Honest null where fusion "
            "disruptions live; quantum's role in fusion is post-~2036")

    _P = 0.2                                  # per-qubit dephasing (realistic sensor decoherence)

    def _visibility(self, N, p):
        """GHZ interferometric visibility (X-parity contrast) under per-qubit dephasing."""
        import pennylane as qml
        dev = qml.device("default.mixed", wires=N)
        obs = qml.PauliX(0)
        for i in range(1, N):
            obs = obs @ qml.PauliX(i)

        @qml.qnode(dev)
        def sig():
            qml.Hadamard(0)
            for i in range(N - 1):
                qml.CNOT([0, i + 1])          # GHZ_N
            for i in range(N):
                qml.PhaseDamping(p, wires=i)  # dephasing channel
            return qml.expval(obs)
        return abs(float(sig()))

    def _scaling(self, Ns=(2, 3, 4, 5, 6, 7, 8)):
        import numpy as np
        rows, peak_N, peak_gain = [], None, -1.0
        for N in Ns:
            Vd = self._visibility(N, self._P)
            gain_ideal, gain_deph, gain_sql = float(N), N * Vd, float(np.sqrt(N))
            if gain_deph > peak_gain:
                peak_gain, peak_N = gain_deph, N
            rows.append({"N": N, "visibility_dephased": round(Vd, 3),
                         "GHZ_ideal_gain": round(gain_ideal, 2),
                         "GHZ_dephased_gain": round(gain_deph, 2),
                         "SQL_gain": round(gain_sql, 2)})
        return rows, peak_N

    def _predict(self, x):
        import numpy as np
        N = int(np.ravel(np.asarray(x, float))[0]) if x is not None else 4
        Vd = self._visibility(max(2, N), self._P)
        return Prediction(round(max(2, N) * Vd, 3), None, True,
                          note=f"GHZ dephased metrological gain at N={max(2,N)} (SQL={np.sqrt(N):.2f})")

    def benchmark(self):
        rows, peak_N = self._scaling()
        a, b = rows[0], rows[-1]
        g_ideal = round(b["GHZ_ideal_gain"] / a["GHZ_ideal_gain"], 2)     # ~ N growth (Heisenberg)
        g_deph = round(b["GHZ_dephased_gain"] / a["GHZ_dephased_gain"], 2)
        g_sql = round(b["SQL_gain"] / a["SQL_gain"], 2)                    # ~ sqrt(N)
        return {"member": self.name,
                "live_quantum_metrology": {
                    "method": "PennyLane mixed-state GHZ interferometry; Heisenberg vs SQL vs dephasing",
                    "dephasing_per_qubit": self._P,
                    "scaling": rows,
                    "gain_growth_N2_to_N8": {"GHZ_ideal": g_ideal, "GHZ_dephased": g_deph, "SQL": g_sql},
                    "analytic_turnover_N": peak_N,
                    "backend": qc.backend_note(),
                    "verdict": (f"REAL quantum-metrology computation: over N=2->8 the ideal GHZ gain grows "
                                f"{g_ideal}x (Heisenberg ~ N), but with per-qubit dephasing (p={self._P}) the "
                                f"dephased GHZ grows only {g_deph}x — essentially the SQL rate ({g_sql}x, "
                                f"~sqrt(N)). Dephasing ERASES the Heisenberg *scaling* to a bounded constant "
                                f"(gain turns over near N~{peak_N + 1}). HONEST NULL: no quantum-sensing "
                                f"advantage for fusion disruptions this decade — but a real, runnable "
                                f"metrology tool.")},
                "sourced": {"file": "TRACK1_QUANTUM_RESULTS/track1_results.json (A7_metrology)",
                            "note": "independent Track-1: correlated dephasing erases GHZ Heisenberg back to SQL"}}


@register
class KTENSOR(Surrogate):
    name = "KTENSOR"; function = "TN"; phase = 2; status = "BUILT"
    provenance = "SIM"
    retired_by = "exact many-body / quantum simulation"
    real_codes = ("tensor-network / low-rank SVD (MPS-style) ROM",)
    gates = ("KX-L3",)
    note = ("tensor-network / low-rank ROM compression of the real CGYRO A1e flux "
            "database (classical MPS-style SVD; a many-body method and bridge to "
            "quantum kernels — NO quantum-advantage claim; quantum's role in fusion is post-~2036)")

    _COLS = ("a_LT", "shear", "Q_i", "Q_e")   # operating grid + saturated fluxes

    def _load(self, turbulent_only=False):
        """Real CGYRO A1e flux database (the COMPLETE 16-pt map; NaN-free) as a
        (points x cols) matrix. The regen 9-pt DB is mostly NaN (a sparse linear scan),
        so we use the completed map, which is the honest, physically-complete dataset."""
        import numpy as np
        df = data.cgyro_flux_map_final()
        if turbulent_only:
            df = df[df["verdict"] == "turbulent"]
        A = df[list(self._COLS)].to_numpy(float)
        A = A[np.isfinite(A).all(1)]                       # drop any non-finite row (safety)
        src = "cgyro_flux_map_final" + (" (turbulent rows)" if turbulent_only else " (16 pts)")
        return A, list(self._COLS), src

    def _matrix(self, turbulent_only=False):
        """Per-column z-scored flux matrix so mixed-scale columns are comparable."""
        import numpy as np
        A, _, _ = self._load(turbulent_only)
        mu, sd = A.mean(0), A.std(0)
        sd = np.where(sd == 0, 1.0, sd)
        return (A - mu) / sd

    @staticmethod
    def _eff_rank(A):
        """Participation ratio of the singular-value spectrum (effective rank)."""
        import numpy as np
        s = np.linalg.svd(A, compute_uv=False)
        return float((s.sum() ** 2) / ((s ** 2).sum() + 1e-30))

    def _svd(self):
        import numpy as np
        A = self._matrix()
        U, s, Vt = np.linalg.svd(A, full_matrices=False)
        return A, U, s, Vt

    def _rel_error(self, rank):
        """rel-L2 (Frobenius) reconstruction error keeping `rank` singular values
        (= MPS bond dimension)."""
        import numpy as np
        A, U, s, Vt = self._svd()
        r = int(max(1, min(int(rank), len(s))))
        Ar = (U[:, :r] * s[:r]) @ Vt[:r]
        return float(np.linalg.norm(A - Ar) / (np.linalg.norm(A) + 1e-12))

    def _predict(self, x):
        """x = target bond dimension (retained rank) -> achievable reconstruction rel-error."""
        import numpy as np
        rank = int(np.ravel(np.asarray(x, dtype=float))[0])
        rmax = len(self._svd()[2])
        return Prediction(self._rel_error(rank), None, 1 <= rank <= rmax,
                          note=f"rel-L2 reconstruction of the CGYRO flux DB at bond-dim "
                               f"{rank} (of {rmax})")

    def benchmark(self):
        import numpy as np

        def curve_for(turb):
            A = self._matrix(turb)
            U, s, Vt = np.linalg.svd(A, full_matrices=False)

            def rel(r):
                r = int(max(1, min(r, len(s))))
                return float(np.linalg.norm(A - (U[:, :r] * s[:r]) @ Vt[:r]) / np.linalg.norm(A))
            var = (s ** 2) / (s ** 2).sum()
            return {"shape": list(A.shape),
                    "rel_error_vs_bond_dim": {r: round(rel(r), 4) for r in range(1, len(s) + 1)},
                    "cumulative_variance": {r: round(float(var[:r].sum()), 4)
                                            for r in range(1, len(s) + 1)},
                    "effective_rank": round(self._eff_rank(A), 2)}

        full = curve_for(False)
        turb = curve_for(True)
        v2 = full["cumulative_variance"][2]
        e2, e3 = full["rel_error_vs_bond_dim"][2], full["rel_error_vs_bond_dim"][3]
        return {"member": self.name,
                "live_low_rank_rom": {
                    "source": "real CGYRO A1e flux database (data.cgyro_flux_map_final)",
                    "columns": list(self._COLS),
                    "full_16pt_map": full, "turbulent_branch": turb,
                    "verdict": (f"MEASURED ROM characterization: the CGYRO flux database is only MODESTLY "
                                f"compressible — effective rank ~{full['effective_rank']} of 4; a bond-dim-2 "
                                f"MPS keeps {v2:.0%} of the variance but {e2:.0%} rel-L2 error (bond-dim 3 -> "
                                f"{e3:.0%}). The turbulent flux spans 3+ orders of magnitude with a sharp "
                                f"turbulent<->quiet transition, so it RESISTS dramatic low-rank compression. "
                                f"An honest counter to the naive 'flux is trivially low-rank' expectation. "
                                f"Classical SVD/MPS; NO quantum-advantage claim."),
                    "caveats": ["NOT a dramatic low-rank collapse: bond-dim-2 keeps ~91% variance but "
                                "~30% rel-L2 error",
                                "classical SVD/MPS ROM of real CGYRO data — a many-body / quantum-bridge "
                                "method, NOT a quantum-advantage claim",
                                "the 9-pt regen flux DB is mostly NaN (a sparse linear scan); the completed "
                                "16-pt map is the honest dataset used here"]},
                "method": "SVD low-rank / MPS bond-dimension compression (deterministic, CPU)",
                "sourced": {"file": "TRACK1_QUANTUM_RESULTS/track1_results.json (A5_MPS_flux)",
                            "note": "the Track-1 'rank-2 ~1%' figure was on the sparse 9x4 regen table; on "
                                    "the complete map the flux is not strongly low-rank (this card)"}}


@register
class KQERN(Surrogate):
    name = "KQERN"; function = "QKERNEL"; phase = 2; status = "BUILT"
    provenance = "SIM"
    retired_by = "fault-tolerant quantum hardware (not available this decade)"
    real_codes = ("quantum kernel", "PennyLane", "real MAST disruption features")
    gates = ("KX-L3",)
    note = ("REAL quantum-kernel classifier — embeds real MAST disruption features into a "
            "quantum feature map and classifies on the fidelity kernel (PennyLane); runs on a "
            "simulator now, real hardware pluggable. Honest: ties the classical kernel — a "
            "validated no-advantage result, but a real, runnable quantum-ML pipeline; role in fusion post-~2036")

    _N = 64
    _FEATS = (0, 6, 7, 8)          # ip_mean, beta_n, li, q95 (subset of KWARD's real features)

    def _data(self):
        import numpy as np
        from .. import get, uq
        X, y, _ = get("KWARD")._real_features()
        X = X[:, list(self._FEATS)]
        rng = np.random.default_rng(uq.SEED)
        pos = np.where(y == 1)[0]; neg = np.where(y == 0)[0]
        k = min(self._N // 2, len(pos), len(neg))
        idx = np.r_[rng.choice(pos, k, False), rng.choice(neg, k, False)]
        rng.shuffle(idx)
        Xs, ys = X[idx], y[idx]
        # robust scale to angles in [0, pi]
        lo, hi = np.nanpercentile(Xs, 5, 0), np.nanpercentile(Xs, 95, 0)
        Xs = np.clip((Xs - lo) / (hi - lo + 1e-9), 0, 1) * np.pi
        ntr = int(0.7 * len(ys))
        return Xs[:ntr], ys[:ntr], Xs[ntr:], ys[ntr:]

    def _kernel(self):
        import pennylane as qml
        nq = len(self._FEATS)
        dev = qc.get_device(wires=nq)

        @qml.qnode(dev)
        def overlap(a, b):
            qml.AngleEmbedding(a, wires=range(nq))
            qml.adjoint(qml.AngleEmbedding)(b, wires=range(nq))
            return qml.probs(wires=range(nq))
        return lambda a, b: float(overlap(a, b)[0])    # |<phi(a)|phi(b)>|^2

    def _gram(self, A, B, kfn):
        import numpy as np
        return np.array([[kfn(a, b) for b in B] for a in A])

    def _fit(self):
        if getattr(self, "_res", None) is not None:
            return self._res
        import numpy as np
        from sklearn.svm import SVC
        from sklearn.metrics import roc_auc_score
        Xtr, ytr, Xte, yte = self._data()
        kfn = self._kernel()
        Ktr = self._gram(Xtr, Xtr, kfn); Kte = self._gram(Xte, Xtr, kfn)
        qsvc = SVC(kernel="precomputed", probability=False).fit(Ktr, ytr)
        q_auc = float(roc_auc_score(yte, qsvc.decision_function(Kte)))
        c = SVC(kernel="rbf", probability=False).fit(Xtr, ytr)          # classical baseline
        c_auc = float(roc_auc_score(yte, c.decision_function(Xte)))
        self._res = {"quantum_kernel_auc": round(q_auc, 3), "classical_rbf_auc": round(c_auc, 3),
                     "n_train": len(ytr), "n_test": len(yte), "n_qubits": len(self._FEATS),
                     "features": ["ip_mean_MA", "beta_n", "li", "q95"], "backend": qc.backend_note()}
        return self._res

    def _predict(self, x):
        r = self._fit()
        return Prediction(r["quantum_kernel_auc"], None, True,
                          note=f"quantum-kernel disruption AUC (vs classical {r['classical_rbf_auc']}); "
                               f"backend={r['backend']['active_backend']}")

    def benchmark(self):
        r = self._fit()
        adv = r["quantum_kernel_auc"] - r["classical_rbf_auc"]
        return {"member": self.name,
                "live_quantum_kernel": {
                    "problem": "MAST disruption classification on a real quantum fidelity kernel",
                    "framework": "PennyLane AngleEmbedding kernel + precomputed-kernel SVM",
                    **{k: r[k] for k in ("quantum_kernel_auc", "classical_rbf_auc", "n_train",
                                         "n_test", "n_qubits", "features")},
                    "quantum_minus_classical_auc": round(adv, 3),
                    "backend": r["backend"],
                    "verdict": (f"REAL quantum kernel classifies MAST disruptions at AUC "
                                f"{r['quantum_kernel_auc']} vs classical RBF {r['classical_rbf_auc']} "
                                f"(Delta {adv:+.3f}) — {'no advantage' if adv <= 0.02 else 'marginal'}: a "
                                f"validated null, but a REAL runnable quantum-ML pipeline on real fusion "
                                f"data. Runs on hardware with KODEX_QC_BACKEND=ibm."),
                    "caveats": ["labels are the heuristic Ip-quench disruption labels (from KWARD)",
                                "quantum kernel ties classical here — no advantage; the value is a real, "
                                "hardware-ready quantum-ML tool, honest that advantage is ~8-10 yr out"]},
                "sourced": {"file": "TRACK1_QUANTUM_RESULTS/track1_results.json (A6_quantum_kernel_MAST)",
                            "note": "independent Track-1 run also found quantum-kernel AUC ~= classical (null)"}}


@register
class KQOPT(Surrogate):
    name = "KQOPT"; function = "QOPT"; phase = 3; status = "BUILT"
    provenance = "SIM"
    retired_by = "fault-tolerant quantum hardware (not available this decade)"
    real_codes = ("QAOA", "PennyLane")
    gates = ("KX-L3",)
    note = ("REAL QAOA quantum optimizer — solves a combinatorial design QUBO on an actual "
            "quantum circuit (PennyLane); runs on a simulator now, real hardware pluggable, and "
            "you can plug in your own QUBO. Honest: standard ~0.7-1.0 approx ratio, no speedup "
            "at this size (advantage post-~2036)")

    # a fixed 6-node design QUBO (illustrative MaxCut-style configuration problem)
    _EDGES = ((0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 0), (0, 2), (3, 5))
    _NQ = 6

    def _hamiltonians(self):
        import pennylane as qml
        cost_h = qml.Hamiltonian([0.5] * len(self._EDGES),
                                 [qml.PauliZ(i) @ qml.PauliZ(j) for (i, j) in self._EDGES])
        mixer_h = qml.Hamiltonian([1.0] * self._NQ, [qml.PauliX(i) for i in range(self._NQ)])
        return cost_h, mixer_h

    def _cut(self, bits):
        return sum(1 for (i, j) in self._EDGES if bits[i] != bits[j])

    def _optimum(self):
        import itertools
        return max(self._cut(b) for b in itertools.product((0, 1), repeat=self._NQ))

    def _run(self, p=2, steps=70):
        if getattr(self, "_res", None) is not None:
            return self._res
        import numpy as np
        import pennylane as qml
        from pennylane import qaoa
        from .. import uq
        cost_h, mixer_h = self._hamiltonians()
        dev = qc.get_device(wires=self._NQ)

        def qaoa_layer(gamma, beta):
            qaoa.cost_layer(gamma, cost_h)
            qaoa.mixer_layer(beta, mixer_h)

        def prep(params):
            for w in range(self._NQ):
                qml.Hadamard(w)
            qml.layer(qaoa_layer, p, params[0], params[1])

        @qml.qnode(dev)
        def energy(params):
            prep(params)
            return qml.expval(cost_h)

        @qml.qnode(dev)
        def probs(params):
            prep(params)
            return qml.probs(wires=range(self._NQ))

        rng = np.random.default_rng(uq.SEED)
        params = qml.numpy.array([rng.uniform(0, np.pi, p), rng.uniform(0, np.pi, p)],
                                 requires_grad=True)
        opt = qml.AdamOptimizer(0.1)
        for _ in range(steps):
            params = opt.step(energy, params)
        pr = np.array(probs(params))
        bits = [int(b) for b in format(int(pr.argmax()), f"0{self._NQ}b")]
        qaoa_cut, opt_cut = self._cut(bits), self._optimum()
        self._res = {"qaoa_cut": qaoa_cut, "optimal_cut": opt_cut,
                     "approx_ratio": round(qaoa_cut / opt_cut, 3), "p": p, "n_qubits": self._NQ,
                     "solution_bitstring": bits, "backend": qc.backend_note()}
        return self._res

    def _predict(self, x):
        r = self._run()
        return Prediction(r["approx_ratio"], None, True,
                          note=f"QAOA approx ratio {r['approx_ratio']} (cut {r['qaoa_cut']}/{r['optimal_cut']}); "
                               f"backend={r['backend']['active_backend']}")

    def benchmark(self):
        r = self._run()
        return {"member": self.name,
                "live_qaoa": {
                    "problem": "6-node design QUBO (MaxCut-style; plug in your own Q matrix)",
                    "framework": "PennyLane QAOA (p=2); runs on simulator now, real hardware pluggable",
                    **{k: r[k] for k in ("qaoa_cut", "optimal_cut", "approx_ratio", "p", "n_qubits",
                                         "solution_bitstring")},
                    "backend": r["backend"],
                    "verdict": (f"REAL QAOA reaches {r['approx_ratio']:.0%} of the brute-force optimum "
                                f"(cut {r['qaoa_cut']}/{r['optimal_cut']}, p={r['p']}, {r['n_qubits']} qubits) "
                                f"on the {r['backend']['active_backend']} backend — a working quantum "
                                f"optimizer, runnable on hardware with KODEX_QC_BACKEND=ibm. HONEST: no "
                                f"speedup vs classical at this size; advantage is ~8-10 yr out.")},
                "sourced": {"file": "TRACK1_QUANTUM_RESULTS/track1_results.json (A4_QAOA)",
                            "note": "independent Track-1 p=1 QAOA reached ~66% on an 8-node QUBO (consistent)"}}


@register
class KDYN(Surrogate):
    name = "KDYN"; function = "QDYN"; phase = 3; status = "BUILT"
    provenance = "SIM"
    retired_by = "fault-tolerant quantum hardware (not available this decade)"
    real_codes = ("Trotterized Hamiltonian simulation", "PennyLane")
    gates = ("KX-L3",)
    note = ("REAL quantum dynamics — Trotterized time-evolution of a transverse-field Ising "
            "'kinetic' Hamiltonian (PennyLane); the honest 1/n_steps cost curve for simulating "
            "plasma-like dynamics on a quantum computer. Runs on a simulator now, hardware pluggable; "
            "quantum's role in fusion is post-~2036")

    _NQ = 3; _J = 1.0; _HX = 0.8; _T = 1.0

    def _ham(self):
        import pennylane as qml
        coeffs = [-self._J] * (self._NQ - 1) + [-self._HX] * self._NQ
        ops = ([qml.PauliZ(i) @ qml.PauliZ(i + 1) for i in range(self._NQ - 1)]
               + [qml.PauliX(i) for i in range(self._NQ)])
        return qml.Hamiltonian(coeffs, ops)

    def _curve(self):
        import numpy as np
        import pennylane as qml
        from scipy.linalg import expm
        H = self._ham()
        coeffs, ops = H.terms()
        mats = [(float(c), qml.matrix(o, wire_order=range(self._NQ))) for c, o in zip(coeffs, ops)]
        Hm = sum(c * M for c, M in mats)
        Uex = expm(-1j * self._T * Hm)
        out = {}
        for n in (1, 2, 4, 8, 16, 32):
            dt = self._T / n
            Ustep = np.eye(2 ** self._NQ, dtype=complex)
            for c, M in mats:
                Ustep = expm(-1j * c * dt * M) @ Ustep
            Utr = np.linalg.matrix_power(Ustep, n)
            out[n] = round(float(np.linalg.norm(Utr - Uex, 2)), 5)
        return out

    def _backend_run(self):
        import pennylane as qml
        H = self._ham(); dev = qc.get_device(wires=self._NQ)

        @qml.qnode(dev)
        def circ(n):
            qml.ApproxTimeEvolution(H, self._T, n)
            return qml.expval(qml.PauliZ(0))
        return float(circ(16))

    def _predict(self, x):
        import numpy as np
        c = self._curve()
        n = int(np.ravel(np.asarray(x, float))[0]) if x is not None else 16
        key = min(c, key=lambda k: abs(k - n))
        return Prediction(c[key], None, True, note=f"Trotter spectral-norm error at {key} steps")

    def benchmark(self):
        c = self._curve()
        z0 = self._backend_run()
        st = sorted(c)
        return {"member": self.name,
                "live_trotter": {
                    "system": f"{self._NQ}-qubit transverse-field Ising (J={self._J}, h={self._HX}), t={self._T}",
                    "framework": "PennyLane Trotter (ApproxTimeEvolution); sim now, real hardware pluggable",
                    "spectral_norm_error_vs_steps": c,
                    "converges_as": "~1/n_steps (1st-order Trotter, as expected)",
                    "backend_expval_Z0_at_16_steps": round(z0, 4),
                    "backend": qc.backend_note(),
                    "verdict": (f"REAL Trotterized quantum dynamics: 1st-order error falls {c[st[0]]} -> "
                                f"{c[st[-1]]} from {st[0]} to {st[-1]} steps (~1/n) — the honest cost curve for "
                                f"simulating plasma-like Hamiltonian dynamics on a quantum computer. Runs on "
                                f"real hardware with KODEX_QC_BACKEND=ibm. No advantage at this size; a real, "
                                f"testable pipeline.")},
                "sourced": {"file": "TRACK1_QUANTUM_RESULTS/track1_results.json (A8_Trotter)",
                            "note": "independent Track-1: 1st-order Trotter error 1.257->0.032 with steps (consistent)"}}