"""
qaoa_adapter.py — Bidirectional QAOA Conversion Adapter Set

Borrows existing Lean-native problem representations and makes them
interactable via QAOA (Quantum Approximate Optimization Algorithm).

Forward:   Problem → QUBO → Ising → Pauli strings → Cirq/Qiskit circuit
Backward:  Measurements → bitstring → solution → problem update

Hooks into existing pipeline:
  - EntropyMeasures.QUBOFormulation (Array Array Q16_16)  → QUBO
  - RotationQUBO.QUBOField (frustration, energyScale)      → QUBO
  - braid_search.build_qubo_matrix (bracket QUBO)          → QUBO
  - qubo_highs.solve_qubo_highs (SA / HiGHS)               → stochastic solver
  - Lean #eval bridge (via eigensolid_lean_bridge)         → live Lean eval

Stochastic abuse: runs QAOA and classical stochastic solvers side-by-side
on the same QUBO for comparison.

All Q16_16 boundary conversion is explicit. No decision or gating logic.
"""

from __future__ import annotations

import json
import math
import time
from dataclasses import dataclass, field, asdict
from pathlib import Path
from typing import Any, Optional

import sys as _sys
_sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "lib"))
_sys.path.insert(0, str(Path(__file__).resolve().parent))
from q16 import Q16_SCALE, from_q16, to_q16
from braid_search import _q16_signed

# ── optional backend imports ──────────────────────────────────────────

try:
    import cirq
    _HAS_CIRQ = True
except ImportError:
    _HAS_CIRQ = False

try:
    from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
    _HAS_QISKIT = True
except ImportError:
    _HAS_QISKIT = False

try:
    import numpy as np
    _HAS_NUMPY = True
except ImportError:
    _HAS_NUMPY = False


# =========================================================================
# I. Data Models
# =========================================================================

@dataclass
class QUBO:
    """Quadratic Unconstrained Binary Optimization problem.

    Minimize  x^T Q x  where x_i ∈ {0, 1}.
    matrix[(i,j)] = coefficient for x_i * x_j (i <= j, upper triangular).
    """
    n: int
    matrix: dict[tuple[int, int], float] = field(default_factory=dict)
    offset: float = 0.0

    def __post_init__(self) -> None:
        keys = list(self.matrix.keys())
        for i, j in keys:
            if i > j:
                self.matrix[(j, i)] = self.matrix.pop((i, j))

    def energy(self, x: list[int]) -> float:
        e = self.offset
        for (i, j), qij in self.matrix.items():
            e += qij * x[i] * x[j]
        return e


@dataclass
class Ising:
    """Ising Hamiltonian: H = Σ h_i s_i + Σ J_ij s_i s_j + offset.

    s_i ∈ {+1, -1}.  J uses upper triangular (i < j).
    """
    n: int
    h: list[float] = field(default_factory=list)
    J: dict[tuple[int, int], float] = field(default_factory=dict)
    offset: float = 0.0

    def energy(self, s: list[int]) -> float:
        if len(s) != self.n:
            raise ValueError(f"expected {self.n} spins, got {len(s)}")
        e = self.offset
        for i in range(self.n):
            e += self.h[i] * s[i]
        for (i, j), Jij in self.J.items():
            e += Jij * s[i] * s[j]
        return e


@dataclass
class PauliSum:
    """Pauli string representation of an Ising Hamiltonian.

    terms: list of (pauli_string, coefficient)
           e.g. ("ZZII", 0.5), ("ZIII", -1.0)
    offset: constant (identity) term
    """
    n: int
    terms: list[tuple[str, float]] = field(default_factory=list)
    offset: float = 0.0


# =========================================================================
# II. Q16_16 Boundary Utilities
# =========================================================================

def q16_to_float(raw: int) -> float:
    """Convert Q16_16 raw integer to float (boundary only)."""
    return from_q16(raw)


def float_to_q16(value: float) -> int:
    """Convert float to Q16_16 raw integer (boundary only)."""
    return to_q16(value)


# =========================================================================
# III-A. Lean QUBOFormulation → QUBO
#   Source: EntropyMeasures.lean — QUBOFormulation { matrix, numVariables }
#   matrix is Array (Array Q16_16): outer = rows, inner = cols
#   objective: E = Σ_i Σ_j Q_ij * x_i * x_j  (both true → add)
# =========================================================================

def lean_qubo_formulation_to_qubo(
    matrix: list[list[int]],
    num_variables: Optional[int] = None,
    offset_q16: int = 0,
) -> QUBO:
    """Convert a Lean QUBOFormulation matrix to a QUBO.

    Args:
        matrix: N×N array of Q16_16 raw integers (from Lean's
                QUBOFormulation.matrix : Array (Array Q16_16))
        num_variables: Number of binary variables (defaults to len(matrix))
        offset_q16: Q16_16 constant term

    Returns:
        QUBO with float coefficients converted at the boundary.
    """
    n = num_variables or len(matrix)
    Q: dict[tuple[int, int], float] = {}
    for i in range(min(n, len(matrix))):
        row = matrix[i] if i < len(matrix) else []
        for j in range(min(n, len(row))):
            qij = row[j]
            if qij != 0:
                key = (min(i, j), max(i, j))
                Q[key] = Q.get(key, 0.0) + q16_to_float(qij)
    return QUBO(n=n, matrix=Q, offset=q16_to_float(offset_q16))


def qubo_to_lean_qubo_formulation(qubo: QUBO) -> list[list[int]]:
    """Convert a QUBO back to Lean QUBOFormulation matrix (Q16_16 ints).

    Returns N×N list of Q16_16 raw integers, suitable for formatting
    into Lean syntax::

        QUBOFormulation.mk #[#[q00, q01, ...], #[q10, q11, ...], ...]
    """
    n = qubo.n
    matrix: list[list[int]] = [[0] * n for _ in range(n)]
    for (i, j), qij in qubo.matrix.items():
        raw = float_to_q16(qij)
        matrix[i][j] = raw
        if i != j:
            matrix[j][i] = raw
    return matrix


# =========================================================================
# III-B. Lean RotationQUBO.QUBOField → QUBO
#   Source: RotationQUBO.lean — QUBOField { frustration, energyScale }
#   fieldEnergy(x) = x² / (1 + δ²) - energyScale
#   isFrustrated(x) = fieldEnergy > 0
#   We discretize x across [-range, +range] into binary variables.
# =========================================================================

def lean_qubo_field_to_qubo(
    frustration_raw: int,
    energy_scale_raw: int,
    n_vars: int = 8,
    field_range: float = 2.0,
) -> QUBO:
    """Discretize a Lean RotationQUBO.QUBOField into a binary QUBO.

    Maps the continuous field energy E(x) = x²/(1+δ²) - energyScale
    onto n_vars binary variables by discretizing x ∈ [-range, +range].

    The QUBO encodes: x = Σ_k s_k * Δ where s_k ∈ {0,1} and
    Δ = 2*range / n_vars.  The energy becomes a quadratic in the s_k.

    Variables: x_i = 1 if discretization point i is selected (one-hot).
    """
    δ = q16_to_float(frustration_raw)
    ε = q16_to_float(energy_scale_raw)
    denom = 1.0 + δ * δ

    Q: dict[tuple[int, int], float] = {}
    for i in range(n_vars):
        xi = -field_range + (2.0 * field_range * i) / max(n_vars - 1, 1)
        Ei = (xi * xi) / denom - ε

        # One-hot diagonal: energy at this point
        Q[(i, i)] = Ei

    return QUBO(n=n_vars, matrix=Q)


# =========================================================================
# III-C. braid_search.build_qubo_matrix → QUBO
#   Source: braid_search.py — build_qubo_matrix(brackets) → dict[(i,j), int]
#   Each bracket has: lower, upper, gap, kappa, phi, admissible
#   bracket_cost = -base + |gap|/10  (Q16_16)
#   crossing_penalty = overlap * 2*Q16_ONE + diversity_reward
# =========================================================================

def braid_search_brackets_to_qubo(
    brackets: list[dict],
    as_q16: bool = False,
) -> QUBO:
    """Convert braid brackets to QUBO using the existing braid_search pipeline.

    Args:
        brackets: List of bracket dicts with lower/upper/gap/kappa/phi/admissible.
        as_q16: If True, keep raw Q16_16 ints; if False, convert to float.

    Returns:
        QUBO problem (8-strand braid crossing selection).
    """
    try:
        _sys.path.insert(0, str(Path(__file__).resolve().parent))
        from braid_search import build_qubo_matrix
    except ImportError:
        return _fallback_bracket_qubo(brackets)

    q16_matrix = build_qubo_matrix(brackets)
    n = max(max(i, j) for (i, j) in q16_matrix) + 1 if q16_matrix else 0

    Q: dict[tuple[int, int], float] = {}
    for (i, j), raw in q16_matrix.items():
        if as_q16:
            Q[(i, j)] = float(raw)
        else:
            Q[(i, j)] = q16_to_float(raw)

    return QUBO(n=n, matrix=Q)


def _fallback_bracket_qubo(brackets: list[dict]) -> QUBO:
    """Fallback if braid_search is not importable."""
    n = len(brackets)
    Q: dict[tuple[int, int], float] = {}
    for i, b in enumerate(brackets):
        gap = q16_to_float(b.get("gap", 32768))
        admissible = b.get("admissible", True)
        base = 1.0 if admissible else 2.0
        Q[(i, i)] = -base + abs(gap) * 0.1
        for j in range(i + 1, n):
            bj = brackets[j]
            overlap = min(b.get("upper", 0), bj.get("upper", 0)) - \
                      max(b.get("lower", 0), bj.get("lower", 0))
            if overlap > 0:
                Q[(i, j)] = 2.0
    return QUBO(n=n, matrix=Q)


# =========================================================================
# III-D. BraidReceipt → QUBO  (8-strand braid)
# =========================================================================

def braid_receipt_to_qubo(receipt: dict) -> QUBO:
    """Convert a BraidReceipt JSON dict to a QUBO over 8 strand variables.

    Variable mapping:
      x_i = 1 if strand i is "active" (selected in crossing), 0 otherwise.

    Cost terms:
      - Diagonal: kappa_i (phase accumulation — lower is better)
      - Off-diagonal: slot collision penalty if two strands share a slot
      - Scar penalty: if scar_absent is false, each scarred strand adds cost
      - Sidon slack bonus: want high slack = low max label used
      - Residual penalty: non-converged strands increase cost
    """
    braid = receipt.get("braid", receipt)
    strands = braid.get("strands", [])
    if not strands:
        strands = _receipt_strands_from_bracket(braid)

    n = len(strands)
    Q: dict[tuple[int, int], float] = {}
    offset = 0.0

    slot_map: dict[int, list[int]] = {}

    for i, s in enumerate(strands):
        kappa = q16_to_float(s.get("kappa", 0))
        phi = q16_to_float(s.get("phi", 0))
        slot = s.get("slot", 0)
        converged = s.get("converged", True)
        residual = q16_to_float(s.get("residual", 0))

        Q[(i, i)] = kappa + 0.1 * abs(phi)

        if not converged and residual > 0:
            Q[(i, i)] = Q[(i, i)] + 10.0 * residual

        slot_map.setdefault(slot, []).append(i)

    for slot, indices in slot_map.items():
        if len(indices) > 1:
            penalty = 50.0
            for a in range(len(indices)):
                for b in range(a + 1, len(indices)):
                    i, j = indices[a], indices[b]
                    Q[(i, j)] = Q.get((i, j), 0.0) + penalty

    scar_absent = receipt.get("scar_absent", braid.get("scar_absent", True))
    if not scar_absent:
        for i in range(n):
            Q[(i, i)] = Q.get((i, i), 0.0) + 100.0

    sidon_slack = receipt.get("sidon_slack", 0)
    max_label = 128 - sidon_slack
    offset += max_label / 128.0

    return QUBO(n=n, matrix=Q, offset=offset)


def _receipt_strands_from_bracket(bracket: dict) -> list[dict]:
    """Synthesize strand array from a minimal bracket-only receipt."""
    default = bracket.get("bracket", bracket)
    gaps = [default.get("gap", 32768)] * 8
    kappas = [default.get("kappa", 0)] * 8
    return [
        {"id": chr(ord("a") + i), "slot": 1 << i, "kappa": kappas[i],
         "phi": 0, "residual": 0, "converged": True}
        for i in range(8)
    ]


# =========================================================================
# III-F. Bidirectional QAOA ↔ Greek HachimojiSubstitution Decoder
#
# Each of the 8 QUBO variables maps to a Greek Hachimoji state with a
# defined phase (45° steps), chirality, and flow direction per the
# Omindirection Logogram contract (OMINDIRECTION_LOGOGRAM_DESIGN_AND_COMPILER.md).
#
# Forward  (LTR, phases 0-135°): Φ Λ Ρ Κ — normal Baker regime
# Backward (RTL, phases 180-315°): Ω Σ Π Ζ — quarantine/tearing regime
#
# Lean source: Semantics/HachimojiSubstitution.lean §5-6
# =========================================================================

# Ordered by bit index 0-7 (matches braid_receipt_to_qubo variable order)
GREEK_STATES: list[str] = ["Φ", "Λ", "Ρ", "Κ", "Ω", "Σ", "Π", "Ζ"]

GREEK_PHASE: dict[str, int] = {
    "Φ": 0, "Λ": 45, "Ρ": 90, "Κ": 135,
    "Ω": 180, "Σ": 225, "Π": 270, "Ζ": 315,
}

def _phase_to_chirality(phase: int) -> str:
    """Omindirection Principle 3: chirality is a projection of phase."""
    if phase in (0, 180):
        return "ambidextrous"
    elif phase < 180:
        return "left"
    else:
        return "right"

def _phase_to_direction(phase: int) -> str:
    """Phases 0-135 = forward (LTR); 180-315 = reverse (RTL)."""
    return "forward" if phase < 180 else "reverse"

# Receipt Bool fields controlled by each Greek state
# (payloadBound, contradictionWitness, tearBoundary, detachedMass, residualLane)
_GREEK_RECEIPT_BITS: dict[str, tuple[bool, bool, bool, bool, bool]] = {
    "Φ": (True,  False, False, False, False),
    "Λ": (True,  False, False, False, False),
    "Ρ": (False, False, False, False, True),
    "Κ": (False, False, False, True,  False),
    "Ω": (True,  True,  True,  True,  True),
    "Σ": (False, False, True,  False, False),
    "Π": (False, False, False, False, False),
    "Ζ": (False, False, False, True,  False),
}

def _greek_to_regime(state: str) -> str:
    """SemanticRegime for a Greek state."""
    if state in ("Φ", "Λ"):
        return "beautifulTopologicalFolding"
    elif state in ("Ρ", "Κ"):
        return "uglyAsymmetricPruning"
    else:
        return "horribleManifoldTearing"


def decode_bitstring_to_greek_states(bits: list[bool]) -> list[dict]:
    """Backward QAOA decoder: measurement bitstring → Greek state atoms.

    Each active bit (True) becomes an omindirectional atom with:
      symbol, phase, chirality, direction, regime.

    Matches Lean: HachimojiSubstitution.bitToGreek / Greek.HachimojiBase.phase
    """
    atoms = []
    for i, active in enumerate(bits[:8]):
        if active:
            sym = GREEK_STATES[i]
            phase = GREEK_PHASE[sym]
            atoms.append({
                "bit": i,
                "symbol": sym,
                "phase": phase,
                "chirality": _phase_to_chirality(phase),
                "direction": _phase_to_direction(phase),
                "regime": _greek_to_regime(sym),
            })
    return atoms


def greek_states_to_logogram_receipt(active_states: list[dict]) -> dict:
    """Backward QAOA decoder: Greek atom list → LogogramReceipt dict.

    Dominant state = lowest phase (most stable).
    Bool fields accumulated from ALL active states (OR-fold).
    Matches Lean: HachimojiSubstitution.fromQAOABitstring

    The receipt dict is compatible with RRCLogogramProjection.LogogramReceipt.
    """
    if not active_states:
        return {
            "shape": "logogramProjection",
            "status": "hold",
            "regime": "uglyAsymmetricPruning",
            "payloadBound": False,
            "contradictionWitness": False,
            "tearBoundary": False,
            "detachedMass": False,
            "residualLane": False,
        }

    # Dominant = lowest phase
    dominant = min(active_states, key=lambda a: a["phase"])

    # Accumulate Bool fields from all active states
    pb  = any(_GREEK_RECEIPT_BITS[a["symbol"]][0] for a in active_states)
    cw  = any(_GREEK_RECEIPT_BITS[a["symbol"]][1] for a in active_states)
    tb  = any(_GREEK_RECEIPT_BITS[a["symbol"]][2] for a in active_states)
    dm  = any(_GREEK_RECEIPT_BITS[a["symbol"]][3] for a in active_states)
    rl  = any(_GREEK_RECEIPT_BITS[a["symbol"]][4] for a in active_states)

    # status = hold if Ζ (bit 7) is active
    zeta_active = any(a["symbol"] == "Ζ" for a in active_states)

    return {
        "shape": "logogramProjection",
        "status": "hold" if zeta_active else "candidate",
        "regime": dominant["regime"],
        "payloadBound": pb,
        "contradictionWitness": cw,
        "tearBoundary": tb,
        "detachedMass": dm,
        "residualLane": rl,
        # Omindirection metadata
        "chirality_witness": {
            "direction": dominant["direction"],
            "handedness": dominant["chirality"],
            "phase": dominant["phase"],
            "placement": "quarantine" if dominant["direction"] == "reverse" else "row",
            "mode": "greek_hachimoji_v1",
        },
    }


def logogram_receipt_to_bitstring(receipt: dict) -> list[bool]:
    """RTL direction: LogogramReceipt → 8-bit Greek signature.

    Matches Lean: HachimojiSubstitution.toQAOABitstring
    This is the backward path for re-encoding into a QUBO update.
    """
    regime = receipt.get("regime", "")
    return [
        bool(receipt.get("payloadBound", False)),                        # Φ
        regime == "beautifulTopologicalFolding",                         # Λ
        bool(receipt.get("residualLane", False)),                        # Ρ
        bool(receipt.get("detachedMass", False)),                        # Κ
        bool(receipt.get("contradictionWitness", False)),                # Ω
        bool(receipt.get("tearBoundary", False)),                        # Σ
        regime == "horribleManifoldTearing",                             # Π
        receipt.get("status", "") == "hold",                             # Ζ
    ]


def qaoa_measurement_to_receipt(bits: list[bool]) -> dict:
    """Full bidi decoder: QAOA bitstring → LogogramReceipt.

    Combines decode_bitstring_to_greek_states + greek_states_to_logogram_receipt.
    This is the primary backward path entry point.
    """
    atoms = decode_bitstring_to_greek_states(bits)
    return greek_states_to_logogram_receipt(atoms)


# =========================================================================
# IV. Forward: LonelyRunner β₀/Scar → QUBO
# =========================================================================

def lonely_runner_scar_to_qubo(
    scar_field: Any,
    *,
    beta0_weight: float = 10.0,
    coverage_weight: float = 1.0,
) -> QUBO:
    """Convert a LonelyRunner scar field to a QUBO over N circle points.

    The scar field μ(t,θ) ∈ {0, 1} indicates uncovered (lonely) regions.

    Variables: x_i = scarred[i] (1 = uncovered, 0 = covered), N points.

    Cost: maximize β₀ (connected components of uncovered region)
          = minimize -β₀ + coverage_penalty

    β₀ = Σ_i x_i * (1 - x_{i-1})   (rising edges, cyclically)

    For an end-to-end QAOA that finds a lonely time:
        H_cost = -β₀ + λ * coverage
        The ground state = max β₀ at minimal coverage = a lonely time.
    """
    if isinstance(scar_field, list) and _HAS_NUMPY:
        import numpy as _np
        scar_field = _np.array(scar_field)

    if _HAS_NUMPY:
        import numpy as _np
        T, N = scar_field.shape
        avg_scar = _np.mean(scar_field, axis=0) > 0.5
    else:
        T = len(scar_field)
        N = len(scar_field[0]) if T > 0 else 0
        avg_scar = [sum(scar_field[t][i] for t in range(T)) / T > 0.5
                    for i in range(N)]

    Q: dict[tuple[int, int], float] = {}

    for i in range(N):
        prev = (i - 1) % N
        Q[(i, i)] = Q.get((i, i), 0.0) - beta0_weight
        Q[(prev, i)] = Q.get((prev, i), 0.0) + beta0_weight

    for i in range(N):
        bias = -coverage_weight if avg_scar[i] > 0.5 else coverage_weight
        Q[(i, i)] = Q.get((i, i), 0.0) + bias

    return QUBO(n=N, matrix=Q)


# =========================================================================
# V. Forward: Goormaghtigh Cost → QUBO
# =========================================================================

def goormaghtigh_cost_to_qubo(
    x_range: tuple[int, int] = (2, 90),
    m_range: tuple[int, int] = (3, 13),
    penalty_weight: float = 1000.0,
) -> QUBO:
    """Convert the Goormaghtigh repunit collision problem to a QUBO.

    Repunit: (x^m - 1)/(x - 1) = 1 + x + ... + x^{m-1}

    Goal: find (x₁,m₁) ≠ (x₂,m₂) s.t. repunit(x₁,m₁) = repunit(x₂,m₂).

    QUBO encoding:
      - One-hot encoding for x and m in each repunit
      - Collision penalty proportional to repunit difference
    """
    x_vals = list(range(x_range[0], x_range[1] + 1))
    m_vals = list(range(m_range[0], m_range[1] + 1))
    nx, nm = len(x_vals), len(m_vals)

    off_x1 = 0
    off_m1 = nx
    off_x2 = nx + nm
    off_m2 = nx + nm + nx

    n = nx + nm + nx + nm
    Q: dict[tuple[int, int], float] = {}

    def _idx(offset: int, k: int) -> int:
        return offset + k

    for offset in (off_x1, off_x2):
        for i in range(nx):
            Q[(_idx(offset, i), _idx(offset, i))] = \
                Q.get((_idx(offset, i), _idx(offset, i)), 0.0) - penalty_weight
        for i in range(nx):
            for j in range(i + 1, nx):
                Q[(_idx(offset, i), _idx(offset, j))] = \
                    Q.get((_idx(offset, i), _idx(offset, j)), 0.0) + 2.0 * penalty_weight

    for offset in (off_m1, off_m2):
        for i in range(nm):
            Q[(_idx(offset, i), _idx(offset, i))] = \
                Q.get((_idx(offset, i), _idx(offset, i)), 0.0) - penalty_weight
        for i in range(nm):
            for j in range(i + 1, nm):
                Q[(_idx(offset, i), _idx(offset, j))] = \
                    Q.get((_idx(offset, i), _idx(offset, j)), 0.0) + 2.0 * penalty_weight

    for xi in range(nx):
        for mi in range(nm):
            r1 = _repunit(x_vals[xi], m_vals[mi])
            for xj in range(nx):
                for mj in range(nm):
                    if xi == xj and mi == mj:
                        continue
                    r2 = _repunit(x_vals[xj], m_vals[mj])
                    diff = abs(r1 - r2)
                    i1 = _idx(off_x1, xi)
                    i2 = _idx(off_m1, mi)
                    i3 = _idx(off_x2, xj)
                    i4 = _idx(off_m2, mj)
                    w = penalty_weight / max(diff, 1.0)
                    Q[(i1, i3)] = Q.get((i1, i3), 0.0) + w * 0.25
                    Q[(i2, i4)] = Q.get((i2, i4), 0.0) + w * 0.25

    return QUBO(n=n, matrix=Q)


def _repunit(x: int, m: int) -> int:
    """Compute repunit value (x^m - 1)//(x - 1)."""
    return (x ** m - 1) // (x - 1)


# =========================================================================
# V-A. Tunable Parameterization of Existing Cost Functions
#
# The QUBO cost functions in braid_search.py use hardcoded coefficients.
# These wrappers let you experiment with alternative parameter sets so
# QAOA results can inform better classical cost models.
#
# DefaultCoeffs mirrors the hardcoded values in braid_search.py:
#   - bracket_cost:   -base + |gap| * gap_reward
#   - crossing_penalty: overlap_penalty - |gap_diff| * gap_reward
# =========================================================================

@dataclass
class BraidCostCoeffs:
    """Tunable coefficients for braid bracket cost functions.

    Mirrors the hardcoded parameters from braid_search.py so any
    of them can be varied independently for grid-search tuning.

    Defaults calibrated via SLOS photonic emulator 2026-06-18.
    See qaoa_adapter.auto_calibrate() for re-calibration.
    """
    base_admissible: float = 8.0      # 8*Q16_ONE — linear bias for admissible brackets
    base_inadmissible: float = 16.0   # 16*Q16_ONE — penalty for inadmissible
    gap_reward_scale: float = 0.01    # reward for large gaps
    overlap_penalty: float = 2.0      # 2*Q16_ONE — conflict cost when ranges overlap
    sa_initial_temp: float = 10.0     # SA starting temperature
    sa_cooling_rate: float = 0.9995   # SA geometric cooling factor
    sa_iterations: int = 5000         # SA iteration budget
    slot_collision_penalty: float = 50.0   # receipt slot conflict cost
    scar_absent_penalty: float = 100.0     # scar-present penalty
    levy_alpha: float = 1.5           # Lévy flight exponent


DEFAULT_COEFFS = BraidCostCoeffs()


def tunable_bracket_cost(bracket: dict, coeffs: BraidCostCoeffs = DEFAULT_COEFFS) -> float:
    """Parameterized replacement for braid_search.bracket_cost.

    Cost = -base + |gap| * gap_reward_scale   (float, not Q16_16)
    """
    admissible = bracket.get("admissible", True)
    base = coeffs.base_admissible if admissible else coeffs.base_inadmissible
    gap = _q16_signed(bracket.get("gap", 32768))
    gap_float = gap / 65536.0
    return -base + abs(gap_float) * coeffs.gap_reward_scale


def tunable_crossing_penalty(
    b1: dict, b2: dict, coeffs: BraidCostCoeffs = DEFAULT_COEFFS,
) -> float:
    """Parameterized replacement for braid_search.crossing_penalty."""
    cost = 0.0
    # Overlap check
    l1 = b1.get("lower", 0)
    u1 = b1.get("upper", 0)
    l2 = b2.get("lower", 0)
    u2 = b2.get("upper", 0)
    if l1 < u2 and l2 < u1:
        cost += coeffs.overlap_penalty
    # Gap diversity reward
    g1 = _q16_signed(b1.get("gap", 32768))
    g2 = _q16_signed(b2.get("gap", 32768))
    gap_diff = abs(g1 - g2) / 65536.0
    cost -= gap_diff * coeffs.gap_reward_scale
    return cost


def build_tunable_qubo_matrix(
    brackets: list[dict],
    coeffs: BraidCostCoeffs = DEFAULT_COEFFS,
) -> dict[tuple[int, int], float]:
    """Build QUBO dict using tunable coefficients.

    Returns float-valued matrix (not Q16_16), directly consumable by
    stochastic_abuse_qubo and qaoa_solve_qubo.
    """
    Q: dict[tuple[int, int], float] = {}
    n = len(brackets)
    for i in range(n):
        Q[(i, i)] = tunable_bracket_cost(brackets[i], coeffs)
        for j in range(i + 1, n):
            c = tunable_crossing_penalty(brackets[i], brackets[j], coeffs)
            Q[(i, j)] = c
            Q[(j, i)] = c
    return Q


def tunable_braid_receipt_to_qubo(
    receipt: dict,
    coeffs: BraidCostCoeffs = DEFAULT_COEFFS,
) -> QUBO:
    """braid_receipt_to_qubo with tunable slot/scar penalties."""
    braid = receipt.get("braid", receipt)
    strands = braid.get("strands", [])
    if not strands:
        strands = _receipt_strands_from_bracket(braid)
    n = len(strands)
    Q: dict[tuple[int, int], float] = {}
    offset = 0.0

    slot_map: dict[int, list[int]] = {}
    for i, s in enumerate(strands):
        kappa = s.get("kappa", 0) / 65536.0
        phi = s.get("phi", 0) / 65536.0
        slot = s.get("slot", 0)
        converged = s.get("converged", True)
        residual = s.get("residual", 0) / 65536.0

        Q[(i, i)] = kappa + 0.1 * abs(phi)
        if not converged and residual > 0:
            Q[(i, i)] += 10.0 * residual
        slot_map.setdefault(slot, []).append(i)

    for _slot, indices in slot_map.items():
        if len(indices) > 1:
            for a in range(len(indices)):
                for b in range(a + 1, len(indices)):
                    i, j = indices[a], indices[b]
                    p = coeffs.slot_collision_penalty
                    Q[(i, j)] = Q.get((i, j), 0.0) + p

    scar_absent = receipt.get("scar_absent", braid.get("scar_absent", True))
    if not scar_absent:
        for i in range(n):
            Q[(i, i)] = Q.get((i, i), 0.0) + coeffs.scar_absent_penalty

    sidon_slack = receipt.get("sidon_slack", 0)
    max_label = 128 - sidon_slack
    offset += max_label / 128.0

    return QUBO(n=n, matrix=Q, offset=offset)


# =========================================================================
# V-B. Tuner: Grid Search Over Parameters
# =========================================================================

def tune_braid_parameters(
    brackets: list[dict],
    param_grid: Optional[list[dict]] = None,
    solver_methods: Optional[list[str]] = None,
    time_per_trial: float = 1.0,
    seed: int = 42,
) -> dict:
    """Grid search over braid cost parameters to find the best QUBO formulation.

    For each parameter set in the grid, builds the QUBO and solves it with
    each solver. Returns the configuration that yields the lowest energy,
    plus a full comparison table.

    Args:
        brackets: List of bracket dicts.
        param_grid: List of BraidCostCoeffs overrides.
            Each entry is a dict with keys from BraidCostCoeffs fields.
            Default: 27 combinations varying base_admissible, overlap_penalty,
            gap_reward_scale.
        solver_methods: Solver methods to evaluate.
            Default: ["sa", "levy"].
        time_per_trial: Time limit per solver trial (seconds).
        seed: RNG seed.

    Returns:
        {
            "best_params": {...},         # winning BraidCostCoeffs
            "best_solver": str,            # winning solver name
            "best_energy": float,          # lowest energy found
            "trials": [                    # each trial result
                {
                    "params": {...},
                    "solver": str,
                    "energy": float,
                    "solution": list[int],
                    "runtime_s": float,
                }
            ],
            "n": int,                      # number of brackets
        }
    """
    if param_grid is None:
        param_grid = []
        for base_mult in [0.5, 1.0, 2.0]:
            for overlap_mult in [0.5, 1.0, 2.0]:
                for gap_scale in [0.05, 0.1, 0.2]:
                    param_grid.append({
                        "base_admissible": base_mult,
                        "base_inadmissible": 2.0 * base_mult,
                        "overlap_penalty": 2.0 * overlap_mult,
                        "gap_reward_scale": gap_scale,
                        "slot_collision_penalty": 50.0,
                        "scar_absent_penalty": 100.0,
                        "sa_initial_temp": 10.0,
                        "sa_cooling_rate": 0.9995,
                        "sa_iterations": 5000,
                        "levy_alpha": 1.5,
                    })

    if solver_methods is None:
        solver_methods = ["sa", "levy"]

    import copy
    trials: list[dict] = []

    for param_dict in param_grid:
        coeffs = BraidCostCoeffs(**param_dict)
        Q = build_tunable_qubo_matrix(brackets, coeffs)
        n = len(brackets)
        qubo = QUBO(n=n, matrix=Q)

        for method in solver_methods:
            t0 = time.time()
            result = stochastic_abuse_qubo(
                qubo, method=method, time_limit=time_per_trial, seed=seed,
            )
            runtime = time.time() - t0

            trials.append({
                "params": param_dict,
                "solver": method,
                "energy": result["energy"],
                "solution": result["solution"],
                "runtime_s": round(runtime, 4),
            })

    # Find best trial (lowest energy)
    best_trial = min(trials, key=lambda t: t["energy"])

    return {
        "best_params": best_trial["params"],
        "best_solver": best_trial["solver"],
        "best_energy": best_trial["energy"],
        "trials": trials,
        "n": len(brackets),
    }


# =========================================================================
# V-C. SA Parameter Tuner
# =========================================================================

def tune_sa_parameters(
    qubo: QUBO,
    temp_grid: Optional[list[float]] = None,
    cooling_grid: Optional[list[float]] = None,
    time_limit: float = 2.0,
    seed: int = 42,
) -> dict:
    """Grid search over SA temperature and cooling rate parameters.

    Runs SA with each (initial_temp, cooling_rate) combination and reports
    the configuration that yields the lowest energy.

    Args:
        qubo: QUBO problem.
        temp_grid: List of initial temperatures to try.
            Default: [1.0, 5.0, 10.0, 20.0, 50.0].
        cooling_grid: List of cooling rates to try.
            Default: [0.99, 0.995, 0.999, 0.9995, 0.9999].
        time_limit: Time limit per SA run.
        seed: RNG seed.

    Returns:
        {
            "best_temp": float,
            "best_cooling": float,
            "best_energy": float,
            "trials": [{"temp": float, "cooling": float, "energy": float}],
        }
    """
    if temp_grid is None:
        temp_grid = [1.0, 5.0, 10.0, 20.0, 50.0]
    if cooling_grid is None:
        cooling_grid = [0.99, 0.995, 0.999, 0.9995, 0.9999]

    trials: list[dict] = []
    Q_dict: dict[tuple[int, int], float] = dict(qubo.matrix)
    n = qubo.n

    for temp in temp_grid:
        for cooling in cooling_grid:
            t0 = time.time()
            x = [0] * n
            import random as _r
            _r.seed(seed)
            x = [_r.randint(0, 1) for _ in range(n)]

            def _energy(xv):
                e = qubo.offset
                for (i, j), qij in Q_dict.items():
                    e += qij * xv[i] * xv[j]
                return e

            current_energy = _energy(x)
            best_x = list(x)
            best_energy = current_energy
            T = temp

            while time.time() - t0 < time_limit:
                i = _r.randrange(n)
                x[i] = 1 - x[i]
                new_energy = _energy(x)
                delta = new_energy - current_energy
                if delta < 0 or _r.random() < math.exp(-delta / max(T, 1e-10)):
                    current_energy = new_energy
                    if current_energy < best_energy:
                        best_x = list(x)
                        best_energy = current_energy
                else:
                    x[i] = 1 - x[i]
                T *= cooling

            trials.append({
                "temp": temp,
                "cooling": cooling,
                "energy": best_energy,
                "runtime_s": round(time.time() - t0, 4),
            })

    best_trial = min(trials, key=lambda t: t["energy"])
    return {
        "best_temp": best_trial["temp"],
        "best_cooling": best_trial["cooling"],
        "best_energy": best_trial["energy"],
        "trials": trials,
    }


# =========================================================================
# V-D. Solver Dispatch Tuner: Auto-pick best solver for a QUBO
# =========================================================================

def solver_dispatch_tuner(
    qubo: QUBO,
    time_limit: float = 2.0,
    seed: int = 42,
) -> dict:
    """Auto-select the best solver for a given QUBO instance.

    Runs all available solvers (SA, HiGHS, Lévy, QAOA describe) for a
    brief warmup and returns the one with the lowest energy, plus a
    recommendation based on problem size.

    Args:
        qubo: QUBO problem.
        time_limit: Per-solver time limit.
        seed: RNG seed.

    Returns:
        {
            "recommended": str,             # solver name
            "recommended_energy": float,
            "warmup": {method: {"energy": float, "runtime_s": float}},
            "heuristic_note": str,          # size-based guidance
        }
    """
    results: dict[str, dict] = {}

    # SA
    sa_res = stochastic_abuse_qubo(qubo, method="sa", time_limit=time_limit, seed=seed)
    results["sa"] = {"energy": sa_res["energy"], "runtime_s": sa_res["runtime_s"]}

    # Lévy
    try:
        levy_res = stochastic_abuse_qubo(qubo, method="levy", time_limit=time_limit, seed=seed)
        results["levy"] = {"energy": levy_res["energy"], "runtime_s": levy_res["runtime_s"]}
    except Exception:
        results["levy"] = {"energy": float("inf"), "runtime_s": 0.0}

    # HiGHS (may fall back to SA internally)
    try:
        highs_res = stochastic_abuse_qubo(qubo, method="highs", time_limit=time_limit, seed=seed)
        results["highs"] = {"energy": highs_res["energy"], "runtime_s": highs_res["runtime_s"]}
    except Exception:
        results["highs"] = {"energy": float("inf"), "runtime_s": 0.0}

    # QAOA describe (offline estimate)
    qaoa_res = qaoa_solve_qubo(qubo, p_layers=1, shots=100, backend="describe")
    results["qaoa"] = {"energy": qaoa_res["energy"], "runtime_s": 0.0, "note": "describe-backend-estimate"}

    # Pick best
    best_solver = min(results, key=lambda s: results[s]["energy"])
    best_energy = results[best_solver]["energy"]

    # Heuristic
    n = qubo.n
    if n <= 5:
        heuristic = "small: HiGHS MIP typically exact and fast"
    elif n <= 20:
        heuristic = "medium: SA or Lévy often best; QAOA may help for specific structure"
    else:
        heuristic = "large: SA or Lévy cheaper; QAOA promising if coupling structure is sparse"

    return {
        "recommended": best_solver,
        "recommended_energy": best_energy,
        "warmup": results,
        "heuristic_note": heuristic,
        "n": n,
    }


# =========================================================================
# V-E. SLOS Photonic Cost Calibration
#
# Uses the Perceval SLOS (Schrödinger Levitated Object Simulator) to
# compute photonic reference costs for bracket configurations. The SLOS
# output distribution over photonic modes serves as a physically-grounded
# cost signal.  We then fit the classical BraidCostCoeffs to match these
# reference costs — effectively "calibrating" the heuristic defaults to
# match what a photonic emulator says.
#
# Requires: perceval (pip install perceval-quandela)
# =========================================================================

try:
    import perceval as _pcvl
    _HAS_PERVERSE = True
except ImportError:
    _HAS_PERVERSE = False


class SlosBraidEvaluator:
    """Evaluate braid brackets using Perceval SLOS photonic simulation.

    Encodes bracket parameters (gap, admissible, lower, upper) as
    phase angles in an M-mode linear optical interferometer. The SLOS
    statevector simulator computes the exact output distribution; we
    extract cost signals from mode probabilities.

    The key insight: admissible brackets with large gaps produce
    different photonic interference patterns (higher mode spread,
    lower entropy) than inadmissible or overlapping ones. SLOS
    captures this physically rather than heuristically.
    """

    def __init__(self, M: int = 6):
        self.M = M
        self._cache: dict[str, float] = {}

    def _cache_key(self, bracket: dict) -> str:
        return f"{bracket.get('lower',0)}|{bracket.get('upper',0)}|{bracket.get('gap',0)}|{bracket.get('admissible',True)}"

    def _encode_bracket_angles(self, bracket: dict) -> list[float]:
        """Map bracket parameters to photonic phase angles.

        Encoding scheme:
          - admissible → phase = 0.0 (no shift — clean path)
          - inadmissible → phase = π/2 (rotation — scattering loss)
          - gap → proportional phase scale (0..π), larger gap = larger shift
          - lower/upper → differential phase between modes
        """
        admissible = bracket.get("admissible", True)
        gap = bracket.get("gap", 32768)
        lower = bracket.get("lower", 0)
        upper = bracket.get("upper", 65536)

        # Normalize to 0..1
        gap_norm = min(1.0, abs(gap) / 65536.0)
        span_norm = min(1.0, max(0, (upper - lower)) / 65536.0)

        theta: list[float] = []
        # Mode 0: admissibility phase
        theta.append(0.0 if admissible else math.pi / 2.0)
        # Mode 1: gap phase
        theta.append(gap_norm * math.pi)
        # Mode 2: span phase
        theta.append(span_norm * math.pi)
        # Modes 3+: fill with differential phases
        for k in range(3, self.M):
            theta.append((gap_norm * span_norm * math.pi) / (1.0 + k))
        return theta

    def compute_photonic_cost(self, bracket: dict) -> float:
        """Run SLOS on a single bracket and return photonic cost.

        Lower photonic cost = more favorable bracket configuration.
        The cost is derived from the output mode distribution:
        - Admissible + large gap → concentrated output (low cost)
        - Inadmissible + small gap → scattered output (high cost)
        """
        ck = self._cache_key(bracket)
        if ck in self._cache:
            return self._cache[ck]

        if not _HAS_PERVERSE:
            # Fallback: heuristic estimate mirrors braid_search defaults
            admissible = bracket.get("admissible", True)
            gap = bracket.get("gap", 32768)
            base = 1.0 if admissible else 2.0
            cost = base + gap / 65536.0 * 0.5
            self._cache[ck] = cost
            return cost

        theta = self._encode_bracket_angles(bracket)

        try:
            circuit = _pcvl.Circuit(self.M)
            for i in range(min(len(theta), self.M)):
                circuit.add(i, _pcvl.PS(theta[i]))
            for i in range(self.M - 1):
                circuit.add((i, i + 1), _pcvl.BS())

            input_state = _pcvl.BasicState([1] + [0] * (self.M - 1))
            processor = _pcvl.Processor("SLOS", circuit)
            processor.with_input(input_state)
            sampler = _pcvl.algorithm.Sampler(processor)
            res = sampler.sample_count(1000)

            # Compute photonic cost from output distribution
            # Metric: weighted entropy of output modes
            total_prob = 0.0
            entropy = 0.0
            for state, count in res["results"].items():
                prob = count / 1000.0
                total_prob += prob
                if prob > 1e-10:
                    entropy -= prob * math.log2(prob)

            # Scale: admissible brackets have lower entropy (cleaner interference)
            admissible = bracket.get("admissible", True)
            cost = entropy * (1.0 if admissible else 1.5)
            self._cache[ck] = cost
        except Exception:
            cost = 1.0 if admissible else 2.0
            self._cache[ck] = cost

        return cost

    def compute_photonic_pair_cost(self, b1: dict, b2: dict) -> float:
        """Compute interaction cost between two brackets using SLOS.

        Encodes both brackets into a 2M-mode interferometer and measures
        the interference cross-term via a cascade BS bridge between the
        two halves. Overlapping brackets produce more cross-half
        entanglement (higher cost); diverse gaps produce cleaner
        separation (lower cost).
        """
        if not _HAS_PERVERSE:
            gap1 = abs(b1.get("gap", 32768)) / 65536.0
            gap2 = abs(b2.get("gap", 32768)) / 65536.0
            gap_diff = abs(gap1 - gap2)
            l1, u1 = b1.get("lower", 0), b1.get("upper", 0)
            l2, u2 = b2.get("lower", 0), b2.get("upper", 0)
            overlap = 1.0 if l1 < u2 and l2 < u1 else 0.0
            return overlap * 2.0 - gap_diff * 0.1

        try:
            n = self.M * 2
            theta1 = self._encode_bracket_angles(b1)
            theta2 = self._encode_bracket_angles(b2)

            circuit = _pcvl.Circuit(n)
            # Phase shifts for bracket 1 on left half (modes 0..M-1)
            for i in range(self.M):
                circuit.add(i, _pcvl.PS(theta1[i] if i < len(theta1) else 0.0))

            # Cascade BS bridge between the two halves (consecutive ports only)
            circuit.add((self.M - 1, self.M), _pcvl.BS())

            # Phase shifts for bracket 2 on right half (modes M..2M-1)
            for i in range(self.M):
                circuit.add(self.M + i, _pcvl.PS(theta2[i] if i < len(theta2) else 0.0))

            input_state = _pcvl.BasicState(
                [1] + [0] * (self.M - 1) + [1] + [0] * (self.M - 1)
            )
            processor = _pcvl.Processor("SLOS", circuit)
            processor.with_input(input_state)
            sampler = _pcvl.algorithm.Sampler(processor)
            res = sampler.sample_count(1000)

            # Interaction cost: probability that photons from both halves
            # end up on the same side (measure of entanglement)
            left_same = 0.0
            right_same = 0.0
            for state, count in res["results"].items():
                prob = count / 1000.0
                left_photons = sum(state[i] for i in range(self.M))
                right_photons = sum(state[i] for i in range(self.M, n))
                if left_photons == 2:
                    left_same += prob
                elif right_photons == 2:
                    right_same += prob

            # Both photons on the same half = stronger interaction
            return (left_same + right_same) * 5.0
        except Exception:
            return 1.0

    def clear_cache(self) -> None:
        self._cache.clear()


# =========================================================================
# V-F. SLOS Calibration: Fit BraidCostCoeffs to Photonic Reference Costs
# =========================================================================

def slos_calibrate_coeffs(
    brackets: list[dict],
    evaluator: Optional[SlosBraidEvaluator] = None,
    loss_fn: str = "mse",
) -> dict:
    """Fit BraidCostCoeffs to match SLOS photonic reference costs.

    For each bracket and each pair, computes the SLOS photonic cost
    (reference).  Then searches the BraidCostCoeffs parameter space
    to find coefficients that minimize the error between the heuristic
    cost (tunable_bracket_cost / tunable_crossing_penalty) and the
    SLOS reference.

    Args:
        brackets: List of bracket dicts to calibrate on
        evaluator: SlosBraidEvaluator instance (created if None)
        loss_fn: "mse" (mean squared error) or "mae" (mean absolute error)

    Returns:
        {
            "calibrated_coeffs": {BraidCostCoeffs fields},
            "reference_costs": list[float],          # SLOS bracket costs
            "reference_pair_costs": list[float],      # SLOS pair costs
            "heuristic_costs": list[float],           # best-fit heuristic costs
            "loss": float,                            # final loss
            "slos_available": bool,                   # whether Perceval was used
        }
    """
    if evaluator is None:
        evaluator = SlosBraidEvaluator()

    n = len(brackets)
    if n == 0:
        return {"calibrated_coeffs": {}, "error": "no brackets"}

    # Compute reference costs from SLOS
    ref_bracket: list[float] = []
    for b in brackets:
        ref_bracket.append(evaluator.compute_photonic_cost(b))

    ref_pairs: list[float] = []
    pair_indices: list[tuple[int, int]] = []
    for i in range(n):
        for j in range(i + 1, n):
            ref_pairs.append(evaluator.compute_photonic_pair_cost(brackets[i], brackets[j]))
            pair_indices.append((i, j))

    # Grid search over BraidCostCoeffs to minimize error
    best_loss = float("inf")
    best_params: dict = {}

    # Parameter sweep ranges — centered around defaults from braid_search.py
    base_range = [0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 4.0, 8.0]
    gap_range = [0.01, 0.025, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5]
    overlap_range = [0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0]

    for base_scale in base_range:
        for gap_scale in gap_range:
            for overlap_scale in overlap_range:
                coeffs = BraidCostCoeffs(
                    base_admissible=base_scale,
                    base_inadmissible=base_scale * 2.0,
                    gap_reward_scale=gap_scale,
                    overlap_penalty=overlap_scale * 2.0,
                )

                # Compute heuristic costs with these coeffs
                heur_bracket = [tunable_bracket_cost(b, coeffs) for b in brackets]
                heur_pairs = [
                    tunable_crossing_penalty(brackets[i], brackets[j], coeffs)
                    for (i, j) in pair_indices
                ]

                # Loss
                loss = 0.0
                for h, r in zip(heur_bracket, ref_bracket):
                    diff = h - r
                    loss += diff * diff if loss_fn == "mse" else abs(diff)
                for h, r in zip(heur_pairs, ref_pairs):
                    diff = h - r
                    loss += diff * diff if loss_fn == "mse" else abs(diff)

                if loss < best_loss:
                    best_loss = loss
                    best_params = asdict(coeffs)

    return {
        "calibrated_coeffs": best_params,
        "reference_costs": ref_bracket,
        "reference_pair_costs": ref_pairs,
        "heuristic_costs": [
            tunable_bracket_cost(b, BraidCostCoeffs(**best_params))
            for b in brackets
        ],
        "loss": best_loss,
        "slos_available": _HAS_PERVERSE,
    }


# =========================================================================
# V-G. Auto-Calibrate: SLOS → Fit → Solve → Compare → Recommend
#
# End-to-end pipeline: for a given set of brackets, runs the SLOS
# photonic emulator to produce physically-grounded reference costs,
# fits the best BraidCostCoeffs to match them, builds QUBOs with
# both default and calibrated coefficients, solves each with all
# available solvers, and recommends the optimal configuration.
# =========================================================================

def auto_calibrate(
    brackets: list[dict],
    evaluator: Optional[SlosBraidEvaluator] = None,
    time_per_solver: float = 1.0,
    seed: int = 42,
) -> dict:
    """Full auto-calibration pipeline for braid QUBO defaults.

    Steps:
      1. Compute SLOS photonic reference costs for each bracket and pair.
      2. Grid-search BraidCostCoeffs to minimize heuristic-vs-SLOS error.
      3. Build QUBOs with default (DEFAULT_COEFFS) and calibrated coeffs.
      4. Solve both QUBOs with SA, Lévy, HiGHS, and QAOA describe.
      5. Report the optimal coefficient set and energy improvement.

    Args:
        brackets: List of bracket dicts.
        evaluator: SlosBraidEvaluator (created fresh if None).
        time_per_solver: Seconds per solver trial.
        seed: RNG seed.

    Returns:
        {
            "slos_calibration": {...},           # raw calibration output
            "coefficients": {
                "default": {...},                # original DEFAULT_COEFFS
                "calibrated": {...},             # SLOS-fitted coeffs
                "recommended": {...},            # whichever gave lowest energy
            },
            "solver_results": {
                "default": {method: result},     # QUBO with default coeffs
                "calibrated": {method: result},  # QUBO with SLOS coeffs
            },
            "comparison": {
                "default_best": {method, energy},
                "calibrated_best": {method, energy},
                "improvement": float,            # energy reduction (negative = better)
            },
            "recommendation": str,               # human-readable summary
            "n": int,
        }
    """
    if evaluator is None:
        evaluator = SlosBraidEvaluator()
    if not brackets:
        return {"error": "no brackets provided"}

    n = len(brackets)

    # Step 1-2: SLOS calibration
    calibration = slos_calibrate_coeffs(brackets, evaluator=evaluator)
    cal_params = calibration["calibrated_coeffs"]
    cal_coeffs = BraidCostCoeffs(**cal_params)

    # Step 3: Build QUBOs with default and calibrated coeffs
    Q_default = build_tunable_qubo_matrix(brackets, DEFAULT_COEFFS)
    Q_cal = build_tunable_qubo_matrix(brackets, cal_coeffs)
    qubo_default = QUBO(n=n, matrix=Q_default)
    qubo_cal = QUBO(n=n, matrix=Q_cal)

    # Step 4: Solve with all methods
    methods = ["sa", "levy"]
    default_results: dict[str, dict] = {}
    cal_results: dict[str, dict] = {}

    for method in methods:
        default_results[method] = stochastic_abuse_qubo(
            qubo_default, method=method, time_limit=time_per_solver, seed=seed,
        )
        cal_results[method] = stochastic_abuse_qubo(
            qubo_cal, method=method, time_limit=time_per_solver, seed=seed,
        )

    # Also try HiGHS
    for method in ["highs"]:
        try:
            default_results[method] = stochastic_abuse_qubo(
                qubo_default, method=method, time_limit=time_per_solver, seed=seed,
            )
        except Exception:
            default_results[method] = {"energy": float("inf"), "solution": [0] * n}
        try:
            cal_results[method] = stochastic_abuse_qubo(
                qubo_cal, method=method, time_limit=time_per_solver, seed=seed,
            )
        except Exception:
            cal_results[method] = {"energy": float("inf"), "solution": [0] * n}

    # Step 5: Find best per coefficient set
    default_best = min(
        default_results.items(), key=lambda kv: kv[1]["energy"],
    )
    cal_best = min(
        cal_results.items(), key=lambda kv: kv[1]["energy"],
    )

    improvement = cal_best[1]["energy"] - default_best[1]["energy"]

    # Decide which coeffs to recommend
    if cal_best[1]["energy"] <= default_best[1]["energy"]:
        recommended = cal_params
        rec_note = "SLOS-calibrated coefficients (lower energy)"
    else:
        recommended = asdict(DEFAULT_COEFFS)
        rec_note = "default coefficients (SLOS calibration did not improve)"

    return {
        "slos_calibration": {
            "loss": calibration["loss"],
            "reference_costs": calibration["reference_costs"],
            "slos_available": calibration["slos_available"],
        },
        "coefficients": {
            "default": asdict(DEFAULT_COEFFS),
            "calibrated": cal_params,
            "recommended": recommended,
        },
        "solver_results": {
            "default": default_results,
            "calibrated": cal_results,
        },
        "comparison": {
            "default_best": {"method": default_best[0], "energy": default_best[1]["energy"]},
            "calibrated_best": {"method": cal_best[0], "energy": cal_best[1]["energy"]},
            "improvement": improvement,
            "improvement_pct": (
                improvement / max(abs(default_best[1]["energy"]), 1e-10) * 100
                if abs(default_best[1]["energy"]) > 1e-10 else 0.0
            ),
        },
        "recommendation": (
            f"Set braid_search defaults to: base_admissible={recommended.get('base_admissible')}, "
            f"gap_reward_scale={recommended.get('gap_reward_scale')}, "
            f"overlap_penalty={recommended.get('overlap_penalty')}. "
            f"{rec_note}."
        ),
        "n": n,
    }


# =========================================================================
# III-D. FinslerMetric → QUBO  (TransportQUBOBridge.lean)
#   Source: TransportTheory.lean — RandersMetric { alpha, beta }
#   alpha = AlphaComponent { dimension, mass_field }
#   beta  = BetaComponent  { dimension, wind_field }
#
#   Q_ij = F(p_i, v_j - v_i) = α(v_j - v_i) + β(v_j - v_i)
#   where α(v) = Σ_i mass[i] * |v_i|  (symmetric: α(-v) = α(v))
#         β(v) = Σ_i wind[i] * v_i     (antisymmetric: β(-v) = -β(v))
#
#   The matrix is NOT symmetric: Q_ij ≠ Q_ji when β ≠ 0.
#   Each row/col = one routing direction.
# =========================================================================

@dataclass
class FinslerMetric:
    """Discrete Finsler-Randers metric (mirrors Lean RandersMetric).

    α component: symmetric base cost (mass field)
    β component: asymmetric drift (wind field)

    F(p,v) = α(p,v) + β(p,v)
    """
    alpha_mass: list[float]  # mass_field: base resistance per dimension
    beta_wind: list[float]   # wind_field: drift 1-form per dimension
    dimension: int           # state space dimension

    def __post_init__(self) -> None:
        if len(self.alpha_mass) != self.dimension:
            raise ValueError(f"alpha mass dim {len(self.alpha_mass)} != {self.dimension}")
        if len(self.beta_wind) != self.dimension:
            raise ValueError(f"beta wind dim {len(self.beta_wind)} != {self.dimension}")

    def alpha_cost(self, v: list[float]) -> float:
        """α(v) = Σ_i mass[i] * |v_i| — symmetric."""
        return sum(self.alpha_mass[i] * abs(v[i]) for i in range(self.dimension))

    def beta_cost(self, v: list[float]) -> float:
        """β(v) = Σ_i wind[i] * v_i — antisymmetric (β(-v) = -β(v))."""
        return sum(self.beta_wind[i] * v[i] for i in range(self.dimension))

    def finsler_cost(self, v: list[float]) -> float:
        """F(p,v) = α(v) + β(v) — direction-dependent Finsler norm."""
        return self.alpha_cost(v) + self.beta_cost(v)

    def crossing_cost(self, v_i: list[float], v_j: list[float]) -> float:
        """Cost of transitioning from direction i to direction j.

        Q_ij = F(p, v_j - v_i) = α(v_j - v_i) + β(v_j - v_i)

        NOTE: Q_ij ≠ Q_ji because β(v_j - v_i) = -β(v_i - v_j).
        """
        diff = [v_j[k] - v_i[k] for k in range(self.dimension)]
        return self.finsler_cost(diff)


def finsler_metric_to_qubo(
    metric: FinslerMetric,
    directions: list[list[float]],
    normalize: bool = True,
) -> QUBO:
    """Convert a FinslerMetric + direction set into a QUBO formulation.

    This mirrors the Lean `randersMetricToQUBO` in TransportQUBOBridge.lean.

    Args:
        metric: FinslerMetric (α mass + β wind).
        directions: List of n direction vectors, each length metric.dimension.
        normalize: If True, normalize costs so the matrix entries are in [-1, 1].

    Returns:
        QUBO with n variables.  Q_ij = crossing cost from i→j, Q_ii = 0.
    """
    n = len(directions)
    Q: dict[tuple[int, int], float] = {}
    max_abs = 0.0

    # Build raw QUBO matrix
    raw: list[list[float]] = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            if i == j:
                raw[i][j] = 0.0
            else:
                cost = metric.crossing_cost(directions[i], directions[j])
                raw[i][j] = cost
                max_abs = max(max_abs, abs(cost))

    # Normalize if requested
    scale = max_abs if normalize and max_abs > 0 else 1.0
    for i in range(n):
        for j in range(i + 1, n):
            val = (raw[i][j] + raw[j][i]) / scale
            if val != 0.0:
                Q[(i, j)] = val

    return QUBO(n=n, matrix=Q)


def is_anisotropic(
    qubo: QUBO,
    raw_matrix: Optional[list[list[float]]] = None,
    tol: float = 1e-9,
) -> bool:
    """Check if a QUBO matrix is anisotropic (Q_ij ≠ Q_ji for any i≠j).

    After the fix in finsler_metric_to_qubo (storing Q_ij + Q_ji),
    the QUBO matrix no longer preserves individual directed entries.
    Pass raw_matrix to compare the original asymmetric entries, or
    pass a FinslerMetric to check the wind field directly.

    Mirrors Lean `isAnisotropic` in TransportQUBOBridge.lean.
    """
    if raw_matrix is not None:
        n = len(raw_matrix)
        for i in range(n):
            for j in range(n):
                if i != j and abs(raw_matrix[i][j] - raw_matrix[j][i]) > tol:
                    return True
        return False
    # Fallback: check stored entries (works for preserved asymmetric storage)
    for (i, j), qij in qubo.matrix.items():
        qji = qubo.matrix.get((j, i), 0.0)
        if abs(qij - qji) > tol:
            return True
    return False


def geodesic_assignment(
    metric: FinslerMetric,
    directions: list[list[float]],
    tol: float = 1e-9,
) -> list[bool]:
    """Extract the Finsler geodesic directions as a Boolean assignment.

    A direction is 'on the geodesic' if its Finsler cost is within tol of
    the global minimum.

    Mirrors Lean `geodesicAssignment` in TransportQUBOBridge.lean.
    """
    costs = [metric.finsler_cost(v) for v in directions]
    min_cost = min(costs)
    return [abs(c - min_cost) <= tol for c in costs]


# =========================================================================
# IV. QUBO → Ising
# =========================================================================

def qubo_to_ising(qubo: QUBO) -> Ising:
    """Convert a QUBO to an Ising Hamiltonian.

    For binary variable x ∈ {0,1}, spin s = 2x - 1 ∈ {+1, -1}:
        x = (1 + s) / 2
        x_i x_j = (1 + s_i + s_j + s_i s_j) / 4

    The QUBO cost:  E = Σ_i a_i x_i + Σ_{i<j} b_ij x_i x_j
    where a_i = Q[i,i], b_ij = Q[i,j] + Q[j,i]
    """
    n = qubo.n
    a: dict[int, float] = {}
    b: dict[tuple[int, int], float] = {}

    for (i, j), qij in qubo.matrix.items():
        if i == j:
            a[i] = a.get(i, 0.0) + qij
        else:
            key = (min(i, j), max(i, j))
            b[key] = b.get(key, 0.0) + qij

    offset = qubo.offset
    h = [0.0] * n
    J: dict[tuple[int, int], float] = {}

    for i in range(n):
        offset += 0.5 * a.get(i, 0.0)
    for (i, j), bij in b.items():
        offset += 0.25 * bij

    for i in range(n):
        hi = 0.5 * a.get(i, 0.0)
        for (j1, j2), bij in b.items():
            if j1 == i:
                hi += 0.25 * bij
            if j2 == i:
                hi += 0.25 * bij
        h[i] = hi

    for (i, j), bij in b.items():
        J[(i, j)] = 0.25 * bij

    return Ising(n=n, h=h, J=J, offset=offset)


# =========================================================================
# VII. Core: Ising → Pauli strings
# =========================================================================

def ising_to_pauli(ising: Ising) -> PauliSum:
    """Convert an Ising Hamiltonian to Pauli strings.

    Mapping:
        s_i    → Z_i
        s_i s_j → Z_i Z_j
        offset  → I (constant)
    """
    terms: list[tuple[str, float]] = []

    for i in range(ising.n):
        if abs(ising.h[i]) > 1e-15:
            ps = ["I"] * ising.n
            ps[i] = "Z"
            terms.append(("".join(ps), ising.h[i]))

    for (i, j), Jij in ising.J.items():
        if abs(Jij) > 1e-15:
            ps = ["I"] * ising.n
            ps[i] = "Z"
            ps[j] = "Z"
            terms.append(("".join(ps), Jij))

    return PauliSum(n=ising.n, terms=terms, offset=ising.offset)


# =========================================================================
# VIII. Forward: Pauli → Cirq Circuit  (QAOA layers)
# =========================================================================

def pauli_to_cirq(
    pauli: PauliSum,
    p_layers: int = 1,
    gamma: Optional[list[float]] = None,
    beta: Optional[list[float]] = None,
    measure: bool = False,
) -> Any:
    """Generate a QAOA circuit from a PauliSum Hamiltonian.

    Returns cirq.Circuit if cirq is available, else a dict description.

    Circuit structure:
        |+⟩^⊗n → [e^{-iγ_k H_C} e^{-iβ_k H_M}]_{k=1..p}
    """
    if not _HAS_CIRQ:
        return _describe_qaoa_circuit(pauli, p_layers)

    n = pauli.n
    qubits = cirq.LineQubit.range(n)

    if gamma is None:
        gamma = [1.0] * p_layers
    if beta is None:
        beta = [1.0] * p_layers

    circuit = cirq.Circuit()
    circuit.append(cirq.H.on_each(*qubits))

    for layer in range(p_layers):
        g = gamma[layer] if layer < len(gamma) else gamma[-1]
        b = beta[layer] if layer < len(beta) else beta[-1]

        for ps_str, coeff in pauli.terms:
            angle = 2.0 * g * coeff
            if abs(angle) < 1e-15:
                continue
            z_positions = [i for i, c in enumerate(ps_str) if c == "Z"]
            if len(z_positions) == 1:
                circuit.append(cirq.rz(angle)(qubits[z_positions[0]]))
            elif len(z_positions) == 2:
                i, j = z_positions
                circuit.append(cirq.CZ(qubits[i], qubits[j]) ** (angle / math.pi))
                circuit.append(cirq.rz(angle)(qubits[i]))
                circuit.append(cirq.rz(angle)(qubits[j]))
            elif len(z_positions) > 2:
                for idx in z_positions:
                    circuit.append(cirq.rz(angle / len(z_positions))(qubits[idx]))

        circuit.append(cirq.rx(2.0 * b).on_each(*qubits))

    if measure:
        circuit.append(cirq.measure(*qubits, key="result"))

    return circuit


def _describe_qaoa_circuit(pauli: PauliSum, p_layers: int = 1) -> dict:
    """Return a JSON-like circuit description when no simulator is available."""
    return {
        "type": "qaoa_circuit_description",
        "n_qubits": pauli.n,
        "p_layers": p_layers,
        "cost_hamiltonian": [{"pauli": t[0], "coeff": t[1]} for t in pauli.terms],
        "offset": pauli.offset,
        "mixer": "X" * pauli.n,
        "note": "Install cirq to generate runnable circuits",
    }


# =========================================================================
# IX. Forward: Pauli → Qiskit Circuit  (QAOA layers)
# =========================================================================

def pauli_to_qiskit(
    pauli: PauliSum,
    p_layers: int = 1,
    gamma: Optional[list[float]] = None,
    beta: Optional[list[float]] = None,
) -> Any:
    """Generate a QAOA circuit using Qiskit.

    Returns QuantumCircuit if qiskit is available, else a dict description.
    """
    if not _HAS_QISKIT:
        return _describe_qaoa_circuit(pauli, p_layers)

    n = pauli.n
    qreg = QuantumRegister(n, "q")
    creg = ClassicalRegister(n, "c")
    circuit = QuantumCircuit(qreg, creg)

    if gamma is None:
        gamma = [1.0] * p_layers
    if beta is None:
        beta = [1.0] * p_layers

    circuit.h(qreg)

    for layer in range(p_layers):
        g = gamma[layer] if layer < len(gamma) else gamma[-1]
        b = beta[layer] if layer < len(beta) else beta[-1]

        for ps_str, coeff in pauli.terms:
            angle = 2.0 * g * coeff
            if abs(angle) < 1e-15:
                continue
            z_positions = [i for i, c in enumerate(ps_str) if c == "Z"]
            if len(z_positions) == 1:
                circuit.rz(angle, qreg[z_positions[0]])
            elif len(z_positions) == 2:
                i, j = z_positions
                circuit.cx(qreg[i], qreg[j])
                circuit.rz(angle, qreg[j])
                circuit.cx(qreg[i], qreg[j])
            elif len(z_positions) > 2:
                for idx in z_positions:
                    circuit.rz(angle / len(z_positions), qreg[idx])

        circuit.rx(2.0 * b, qreg)

    circuit.measure(qreg, creg)
    return circuit


# =========================================================================
# X. Backward: Measurements → Solution
# =========================================================================

def measurements_to_solution(
    counts: dict[str, int],
    n: int,
) -> list[int]:
    """Extract the most probable bitstring from measurement counts.

    Returns binary list [x_0, ..., x_{n-1}] of the most probable outcome.
    """
    if not counts:
        return [0] * n
    best_bitstring = max(counts, key=counts.get)
    if len(best_bitstring) < n:
        best_bitstring = best_bitstring.zfill(n)
    return [int(b) for b in best_bitstring[-n:]]


def measurements_to_ising_solution(
    counts: dict[str, int], n: int,
) -> list[int]:
    """Extract the most probable spin configuration from measurements."""
    x = measurements_to_solution(counts, n)
    return [2 * xi - 1 for xi in x]


# =========================================================================
# XI. Backward: Solution → BraidReceipt Update
# =========================================================================

def solution_to_braid_receipt(solution: list[int], receipt: dict) -> dict:
    """Update a BraidReceipt dict with a QAOA solution."""
    updated = json.loads(json.dumps(receipt))
    braid = updated.get("braid", updated)
    strands = braid.get("strands", [])
    if not strands:
        strands = _receipt_strands_from_bracket(braid)
        braid["strands"] = strands

    n = min(len(solution), len(strands))
    all_admissible = True
    for i in range(n):
        s = strands[i]
        if bool(solution[i]):
            s["kappa"] = max(0, s.get("kappa", 0) // 2)
            s["converged"] = True
            s["residual"] = 0
        else:
            s["kappa"] = min(65535, s.get("kappa", 0) + 8192)
            s["converged"] = False
            s["residual"] = s.get("kappa", 0)
            all_admissible = False

    braid["bracket"] = braid.get("bracket", {})
    braid["bracket"]["admissible"] = all_admissible
    updated["scar_absent"] = all_admissible
    updated["step_count"] = updated.get("step_count", 0) + 1

    residuals = updated.get("residuals", [])
    avg_residual = sum(abs(s.get("kappa", 0)) for s in strands[:n]) // max(n, 1)
    residuals.insert(0, avg_residual)
    updated["residuals"] = residuals[:32]
    return updated


# =========================================================================
# XII. Backward: Solution → Scar/β₀ Update
# =========================================================================

def solution_to_scar_field(solution: list[int], N: int) -> list[int]:
    """Map a QAOA solution bitstring back to a scar field."""
    mask = [0] * N
    for i in range(N):
        if solution[i] == 1:
            mask[i] = 1
    radius = max(1, N // 32)
    expanded = [0] * N
    for i in range(N):
        if mask[i]:
            for d in range(-radius, radius + 1):
                expanded[(i + d) % N] = 1
    return expanded


def compute_beta0_from_solution(solution: list[int]) -> int:
    """Compute β₀ (rising edge count) from a cyclic binary solution."""
    n = len(solution)
    if n == 0:
        return 0
    if all(s == 1 for s in solution):
        return 1
    if all(s == 0 for s in solution):
        return 0
    rising = 0
    for i in range(n):
        if solution[i] == 1 and solution[(i - 1) % n] == 0:
            rising += 1
    return rising


# =========================================================================
# XIII. High-Level: Corpus250 Row → QUBO
#   FixtureRow has: shape, rrcKind, operatorTokens, weakAxesCnt,
#   pistProxyLabel, pistExactLabel (no domain_embedding field).
# =========================================================================

def corpus278_row_to_qubo(row: dict) -> QUBO:
    """Convert a Corpus250 equation row (FixtureRow) to a QUBO problem.

    Uses the actual FixtureRow fields:
      - shape: RRCShape enum guiding the QUBO structure
      - operatorTokens: domain operators for variable count/diagonal
      - weakAxesCnt: number of weak axes for off-diagonal coupling
      - pistProxyLabel / pistExactLabel: PIST labels for bias
    """
    shape = row.get("shape", "")
    operator_tokens = row.get("operatorTokens", [])
    weak_axes_cnt = row.get("weakAxesCnt", 0)

    n_vars = max(8, len(operator_tokens) + int(weak_axes_cnt))

    shape_bias = {
        "cognitiveLoadField": 1.0,
        "signalShapedRouteCompiler": 2.0,
        "projectableGeometryTopology": 3.0,
        "cadForceProbeReceipt": 4.0,
        "logogramProjection": 5.0,
        "holdForUnlawfulOrUnderspecifiedShape": 10.0,
    }.get(shape, 2.0)

    Q: dict[tuple[int, int], float] = {}

    # Diagonal: one per operator token with shape bias
    for i in range(min(n_vars, len(operator_tokens))):
        Q[(i, i)] = shape_bias * (1.0 + (i % (int(weak_axes_cnt) + 1)))

    # Fill remaining diagonal
    for i in range(len(operator_tokens), n_vars):
        Q[(i, i)] = shape_bias * 1.5

    # Weak axes coupling
    wac = int(weak_axes_cnt)
    if wac > 0:
        for i in range(min(wac, n_vars)):
            for j in range(i + 1, min(wac, n_vars)):
                Q[(i, j)] = Q.get((i, j), 0.0) - 0.5

    # PIST label bias
    pist_proxy = row.get("pistProxyLabel")
    if pist_proxy and pist_proxy != "none":
        for i in range(min(4, n_vars)):
            Q[(i, i)] = Q.get((i, i), 0.0) - 0.3

    return QUBO(n=n_vars, matrix=Q)


# =========================================================================
# XIV. Stochastic Abuse: Classical Stochastic Solvers
#   Wraps existing qubo_highs.py (SA / HiGHS) so results are directly
#   comparable with qaoa_solve_qubo() output.
#
#   "Abuse" = using the classical stochastic pipeline as a baseline to
#   benchmark QAOA against, on the same QUBO — borrowing the existing
#   HiGHS MIP and simulated annealing infrastructure.
# =========================================================================

def stochastic_abuse_qubo(
    qubo: QUBO,
    method: str = "sa",
    time_limit: float = 5.0,
    seed: int = 42,
) -> dict:
    """Solve a QUBO using the existing classical stochastic pipeline.

    Wraps qubo_highs.solve_qubo_highs (HiGHS MIP) and
    braid_search.qubo_optimize (SA) for direct comparison with QAOA.

    Args:
        qubo: QUBO problem
        method: "sa" (simulated annealing), "highs" (HiGHS MIP), or
                "levy" (Lévy flight sampling)
        time_limit: Max seconds for solver
        seed: RNG seed

    Returns:
        dict with same schema as qaoa_solve_qubo() for easy comparison:
          - solution, energy, method, runtime, solver_detail
    """
    t0 = time.time()
    n = qubo.n

    # Build QUBO dict in the format qubo_highs expects
    Q_dict: dict[tuple[int, int], float] = dict(qubo.matrix)

    if method == "highs":
        try:
            _sys.path.insert(0, str(Path(__file__).resolve().parent))
            from qubo_highs import solve_qubo_highs as highs_solve
            result = highs_solve(Q_dict, n, time_limit=time_limit)
            if isinstance(result, dict):
                solution = result.get("x", [0] * n)
            else:
                solution, _ = result
        except Exception as exc:
            solution = _sa_solve(Q_dict, n, seed, time_limit)

    elif method == "levy":
        solution = _levy_fight(Q_dict, n, seed, time_limit)

    else:
        solution = _sa_solve(Q_dict, n, seed, time_limit)

    runtime = time.time() - t0
    energy = qubo.energy(solution)

    return {
        "solution": solution,
        "energy": energy,
        "method": method,
        "runtime_s": round(runtime, 4),
        "n": n,
    }


def _sa_solve(
    Q: dict[tuple[int, int], float],
    n: int,
    seed: int = 42,
    time_limit: float = 5.0,
) -> list[int]:
    """Simulated annealing QUBO solver (standalone fallback)."""
    import random as _r
    _r.seed(seed)
    t0 = time.time()

    x = [_r.randint(0, 1) for _ in range(n)]

    def _energy(xv):
        e = 0.0
        for (i, j), qij in Q.items():
            e += qij * xv[i] * xv[j]
        return e

    current_energy = _energy(x)
    best_x = list(x)
    best_energy = current_energy
    temp = 10.0
    iterations = 0
    n_vals = list(range(n))

    while time.time() - t0 < time_limit:
        i = _r.choice(n_vals)
        x[i] = 1 - x[i]
        new_energy = _energy(x)
        delta = new_energy - current_energy

        if delta < 0 or _r.random() < math.exp(-delta / max(temp, 1e-10)):
            current_energy = new_energy
            if current_energy < best_energy:
                best_x = list(x)
                best_energy = current_energy
        else:
            x[i] = 1 - x[i]

        temp *= 0.9995
        iterations += 1

    return best_x


def _levy_fight(
    Q: dict[tuple[int, int], float],
    n: int,
    seed: int = 42,
    time_limit: float = 5.0,
) -> list[int]:
    """Lévy flight sampling over binary strings.

    Step sizes follow a power-law distribution (heavy-tailed),
    mimicking the LévyFlight structure from EntropyMeasures.lean.
    """
    import random as _r
    _r.seed(seed)
    t0 = time.time()

    x = [_r.randint(0, 1) for _ in range(n)]

    def _energy(xv):
        e = 0.0
        for (i, j), qij in Q.items():
            e += qij * xv[i] * xv[j]
        return e

    def _levy_step() -> int:
        u = _r.random()
        return max(1, int(n * u ** (-1.0 / 1.5)) % n)

    best_x = list(x)
    best_energy = _energy(x)

    while time.time() - t0 < time_limit:
        step = _levy_step()
        indices = _r.sample(range(n), min(step, n))
        for i in indices:
            x[i] = 1 - x[i]
        e = _energy(x)
        if e < best_energy:
            best_x = list(x)
            best_energy = e
        else:
            for i in indices:
                x[i] = 1 - x[i]

    return best_x


# =========================================================================
# XV. Comparison: QAOA vs Stochastic
# =========================================================================

def qaoa_vs_stochastic_comparison(
    qubo: QUBO,
    p_layers: int = 1,
    shots: int = 1000,
    qaoa_backend: str = "cirq",
    stochastic_methods: Optional[list[str]] = None,
    time_limit: float = 5.0,
    seed: int = 42,
) -> dict:
    """Run QAOA and classical stochastic solvers on the same QUBO.

    Returns a side-by-side comparison so results can be contrasted.

    Args:
        qubo: QUBO problem
        p_layers: QAOA depth
        shots: QAOA measurement shots
        qaoa_backend: QAOA backend ("cirq", "qiskit", "describe")
        stochastic_methods: List of classical methods to compare.
            Default: ["sa", "levy"].
        time_limit: Time limit per classical solver (seconds)
        seed: RNG seed

    Returns:
        dict with:
          - qubo_summary: n, terms, offset
          - qaoa: result from qaoa_solve_qubo
          - stochastic: {method: result} for each method
          - winner: solver with lowest energy
    """
    if stochastic_methods is None:
        stochastic_methods = ["sa", "levy"]

    qaoa_result = qaoa_solve_qubo(
        qubo, p_layers=p_layers, shots=shots, backend=qaoa_backend,
    )

    stochastic_results: dict[str, dict] = {}
    for method in stochastic_methods:
        stochastic_results[method] = stochastic_abuse_qubo(
            qubo, method=method, time_limit=time_limit, seed=seed,
        )

    # Determine winner
    candidates: list[tuple[str, float]] = [
        ("qaoa", qaoa_result["energy"]),
    ]
    for method, res in stochastic_results.items():
        candidates.append((method, res["energy"]))

    candidates.sort(key=lambda p: p[1])
    winner_name, winner_energy = candidates[0]

    return {
        "qubo_summary": {
            "n": qubo.n,
            "terms": len(qubo.matrix),
            "offset": qubo.offset,
        },
        "qaoa": {
            "solution": qaoa_result["solution"],
            "energy": qaoa_result["energy"],
            "p_layers": p_layers,
            "shots": shots,
            "backend": qaoa_backend,
        },
        "stochastic": stochastic_results,
        "winner": {
            "solver": winner_name,
            "energy": winner_energy,
        },
    }


# =========================================================================
# XVI. High-Level: End-to-End QAOA Solve
# =========================================================================

def qaoa_solve_qubo(
    qubo: QUBO,
    p_layers: int = 1,
    gamma: Optional[list[float]] = None,
    beta: Optional[list[float]] = None,
    shots: int = 1000,
    backend: str = "cirq",
) -> dict:
    """Solve a QUBO using QAOA.

    Args:
        qubo: QUBO problem
        p_layers: QAOA layers
        gamma: cost angles
        beta: mixer angles
        shots: measurement shots
        backend: "cirq", "qiskit", or "describe"

    Returns:
        dict with solution, energy, counts, and circuit description
    """
    ising = qubo_to_ising(qubo)
    pauli = ising_to_pauli(ising)

    result: dict[str, Any] = {
        "n": qubo.n,
        "p_layers": p_layers,
        "qubo_energy_offset": qubo.offset,
        "ising_offset": ising.offset,
        "pauli_terms": [(t[0], t[1]) for t in pauli.terms],
    }

    if backend == "describe":
        result["circuit"] = _describe_qaoa_circuit(pauli, p_layers)
        result["solution"] = [0] * qubo.n
        result["energy"] = qubo.energy([0] * qubo.n)
        result["note"] = "describe mode — no simulation"
        return result

    if backend == "cirq" and _HAS_CIRQ:
        circuit = pauli_to_cirq(pauli, p_layers, gamma, beta, measure=True)
        simulator = cirq.Simulator()
        samples = simulator.run(circuit, repetitions=shots)
        counts = samples.histogram(key="result")
        str_counts = _cirq_counts_to_str(counts, qubo.n)
        solution = measurements_to_solution(str_counts, qubo.n)
        energy = qubo.energy(solution)
        result["circuit"] = str(circuit)
        result["counts"] = str_counts
        result["solution"] = solution
        result["energy"] = energy
        return result

    if backend == "qiskit" and _HAS_QISKIT:
        circuit = pauli_to_qiskit(pauli, p_layers, gamma, beta)
        from qiskit_aer import AerSimulator
        simulator = AerSimulator()
        job = simulator.run(circuit, shots=shots)
        counts_dict = job.result().get_counts()
        solution = measurements_to_solution(counts_dict, qubo.n)
        energy = qubo.energy(solution)
        result["circuit"] = circuit.qasm()
        result["counts"] = counts_dict
        result["solution"] = solution
        result["energy"] = energy
        return result

    import random
    random.seed(0)
    all_counts: dict[str, int] = {}
    for _ in range(shots):
        x = [random.randint(0, 1) for _ in range(qubo.n)]
        bits = "".join(str(b) for b in x)
        all_counts[bits] = all_counts.get(bits, 0) + 1

    solution = measurements_to_solution(all_counts, qubo.n)
    energy = qubo.energy(solution)
    result["circuit"] = _describe_qaoa_circuit(pauli, p_layers)
    result["counts"] = all_counts
    result["solution"] = solution
    result["energy"] = energy
    result["note"] = "random fallback — install cirq for QAOA simulation"
    return result


def _cirq_counts_to_str(counts: dict, n: int) -> dict[str, int]:
    """Convert Cirq integer keyed counts to bitstring counts."""
    str_counts: dict[str, int] = {}
    for val, cnt in counts.items():
        bits = format(val, f"0{n}b")
        str_counts[bits] = cnt
    return str_counts


# =========================================================================
# XVII. High-Level: Roundtrip
# =========================================================================

def qaoa_roundtrip(
    receipt: dict,
    p_layers: int = 1,
    shots: int = 1000,
    backend: str = "cirq",
    compare_stochastic: bool = False,
) -> dict:
    """Full roundtrip: BraidReceipt → QAOA → updated BraidReceipt.

    Args:
        receipt: Input BraidReceipt JSON dict
        p_layers: QAOA depth
        shots: Measurement shots per layer
        backend: "cirq", "qiskit", "describe"
        compare_stochastic: If True, also run SA/Levy comparison

    Returns:
        dict with input_receipt, output_receipt, qaoa_result, delta,
        and optionally stochastic comparison.
    """
    qubo = braid_receipt_to_qubo(receipt)
    qaoa_result = qaoa_solve_qubo(
        qubo, p_layers=p_layers, shots=shots, backend=backend,
    )
    solution = qaoa_result["solution"]
    output_receipt = solution_to_braid_receipt(solution, receipt)

    input_kappas = [
        s.get("kappa", 0)
        for s in receipt.get("braid", receipt).get("strands", [])
    ]
    output_kappas = [
        s.get("kappa", 0)
        for s in output_receipt.get("braid", output_receipt).get("strands", [])
    ]

    result: dict[str, Any] = {
        "input_receipt": receipt,
        "output_receipt": output_receipt,
        "qaoa_result": {
            "solution": solution,
            "energy": qaoa_result["energy"],
            "counts": qaoa_result.get("counts", {}),
            "p_layers": p_layers,
        },
        "delta": {
            "scar_absent": {
                "before": receipt.get("scar_absent"),
                "after": output_receipt.get("scar_absent"),
            },
            "step_count": {
                "before": receipt.get("step_count", 0),
                "after": output_receipt.get("step_count", 0),
            },
            "avg_kappa": {
                "before": sum(input_kappas) / max(len(input_kappas), 1),
                "after": sum(output_kappas) / max(len(output_kappas), 1),
            },
        },
    }

    if compare_stochastic:
        comparison = qaoa_vs_stochastic_comparison(qubo)
        result["stochastic_comparison"] = comparison

    return result


# =========================================================================
# XVIII-A. FinslerMetric demo
# =========================================================================

def finsler_demo() -> dict[str, Any]:
    """End-to-end demo of the FinslerMetric → QUBO → Ising → Pauli pipeline.

    Constructs a simple 2D Finsler metric with nontrivial wind (β ≠ 0),
    generates direction vectors on a circle, builds the QUBO matrix,
    converts to Ising and Pauli strings, and verifies anisotropy.

    This mirrors the Lean `trivialRanders` witness in TransportQUBOBridge.lean.
    """
    # 2D Finsler metric with wind drift β = (0.25, 0.0)
    metric = FinslerMetric(
        alpha_mass=[1.0, 1.0],   # uniform mass
        beta_wind=[0.25, 0.0],    # rightward drift
        dimension=2,
    )

    # 8 directions on the unit circle
    import math
    n_dirs = 8
    directions = []
    for k in range(n_dirs):
        theta = 2.0 * math.pi * k / n_dirs
        directions.append([math.cos(theta), math.sin(theta)])

    # Build QUBO (now stores Q_ij + Q_ji summed per pair)
    qubo = finsler_metric_to_qubo(metric, directions, normalize=False)

    # Check anisotropy via raw matrix (Q_ij vs Q_ji before summation)
    n_raw = len(directions)
    raw = [[0.0] * n_raw for _ in range(n_raw)]
    for i in range(n_raw):
        for j in range(n_raw):
            if i != j:
                raw[i][j] = metric.crossing_cost(directions[i], directions[j])
    anisotropic = is_anisotropic(qubo, raw_matrix=raw)

    # Convert to Ising
    ising = qubo_to_ising(qubo)

    # Convert to Pauli
    pauli = ising_to_pauli(ising)

    # Find geodesic assignment
    geo = geodesic_assignment(metric, directions)

    # Summary
    return {
        "metric": {
            "alpha_mass": metric.alpha_mass,
            "beta_wind": metric.beta_wind,
            "dimension": metric.dimension,
        },
        "n_directions": n_dirs,
        "qubo_is_anisotropic": anisotropic,
        "qubo_n": qubo.n,
        "qubo_nonzero_entries": len(qubo.matrix),
        "ising_h": ising.h[:5] if ising.n > 5 else ising.h,
        "ising_J_entries": len(ising.J),
        "pauli_terms": len(pauli.terms),
        "geodesic_assignment": geo,
        "anisotropic_pair": _find_anisotropic_pair(qubo, raw_matrix=raw),
    }


def _find_anisotropic_pair(qubo: QUBO, raw_matrix: Optional[list[list[float]]] = None) -> Optional[dict]:
    """Find the first (i,j) where Q_ij ≠ Q_ji."""
    if raw_matrix is not None:
        n = len(raw_matrix)
        for i in range(n):
            for j in range(n):
                if i != j and abs(raw_matrix[i][j] - raw_matrix[j][i]) > 1e-9:
                    return {"i": i, "j": j, "Q_ij": raw_matrix[i][j], "Q_ji": raw_matrix[j][i], "delta": raw_matrix[i][j] - raw_matrix[j][i]}
        return None
    for (i, j), qij in qubo.matrix.items():
        if i == j:
            continue
        qji = qubo.matrix.get((j, i), 0.0)
        if abs(qij - qji) > 1e-9:
            return {"i": i, "j": j, "Q_ij": qij, "Q_ji": qji, "delta": qij - qji}
    return None


# =========================================================================
# XVIII. CLI
# =========================================================================

def _parse_receipt(path: str) -> dict:
    with open(path) as f:
        return json.load(f)


def _serialize(obj: Any) -> Any:
    """Recursively serialize a result object to JSON-safe types.

    Handles QUBO, Ising, PauliSum dataclasses and dicts with tuple keys.
    """
    if hasattr(obj, "__dataclass_fields__"):
        d = asdict(obj)
        return _serialize(d)
    if isinstance(obj, dict):
        out = {}
        for k, v in obj.items():
            if isinstance(k, tuple):
                sk = f"({','.join(str(x) for x in k)})"
            elif isinstance(k, int):
                sk = str(k)
            else:
                sk = k
            out[sk] = _serialize(v)
        return out
    if isinstance(obj, (list, tuple)):
        return [_serialize(x) for x in obj]
    return obj


def main() -> None:
    import argparse

    parser = argparse.ArgumentParser(
        description="Bidirectional QAOA Conversion Adapter Set"
    )
    parser.add_argument(
        "action",
        choices=[
            "lean-formulation", "lean-field", "braid-search",
            "braid-to-qubo", "lone-to-qubo", "goorm-to-qubo",
            "qubo-to-ising", "ising-to-pauli", "pauli-to-cirq",
            "pauli-to-qiskit", "solve-qubo", "stochastic",
            "compare", "roundtrip",
            "tune", "auto-solve", "slos-calibrate", "auto-calibrate",
            "finsler-demo",
        ],
    )
    parser.add_argument("--coeffs", help="JSON file with BraidCostCoeffs overrides")
    parser.add_argument("--receipt", help="Path to BraidReceipt JSON")
    parser.add_argument("--output", "-o", help="Output path (default: stdout)")
    parser.add_argument("--p-layers", type=int, default=1, help="QAOA layers")
    parser.add_argument("--shots", type=int, default=1000, help="Measurement shots")
    parser.add_argument("--backend", default="cirq", help="Simulator backend")
    parser.add_argument("--gamma", type=float, nargs="*", help="Cost angles")
    parser.add_argument("--beta", type=float, nargs="*", help="Mixer angles")
    parser.add_argument("--method", default="sa", help="Stochastic method")
    parser.add_argument("--time-limit", type=float, default=5.0,
                        help="Stochastic solver time limit")

    args = parser.parse_args()

    result: Any = None

    if args.action == "lean-formulation":
        result = {"note": "Use Python API with lean_qubo_formulation_to_qubo(matrix, n)"}
    elif args.action == "lean-field":
        result = lean_qubo_field_to_qubo(65536, 10 * 65536)
    elif args.action == "braid-search":
        if args.receipt:
            receipt = _parse_receipt(args.receipt)
            bkts = receipt.get("braid", receipt).get("strands", [])
            result = braid_search_brackets_to_qubo(bkts)
        else:
            result = {"note": "Use --receipt with bracket-containing JSON"}
    elif args.action == "braid-to-qubo":
        receipt = _parse_receipt(args.receipt)
        result = braid_receipt_to_qubo(receipt)
    elif args.action == "lone-to-qubo":
        result = {"note": "Use Python API with a numpy scar field array"}
    elif args.action == "goorm-to-qubo":
        result = goormaghtigh_cost_to_qubo()
    elif args.action == "qubo-to-ising":
        result = {"note": "Use Python API with a QUBO object"}
    elif args.action == "ising-to-pauli":
        result = {"note": "Use Python API with an Ising object"}
    elif args.action == "pauli-to-cirq":
        result = {"note": "Use Python API with a PauliSum object"}
    elif args.action == "pauli-to-qiskit":
        result = {"note": "Use Python API with a PauliSum object"}
    elif args.action == "solve-qubo":
        receipt = _parse_receipt(args.receipt)
        result = qaoa_roundtrip(receipt, args.p_layers, args.shots, args.backend)
    elif args.action == "stochastic":
        receipt = _parse_receipt(args.receipt)
        qubo = braid_receipt_to_qubo(receipt)
        result = stochastic_abuse_qubo(qubo, args.method, args.time_limit)
    elif args.action == "compare":
        receipt = _parse_receipt(args.receipt)
        qubo = braid_receipt_to_qubo(receipt)
        result = qaoa_vs_stochastic_comparison(
            qubo, args.p_layers, args.shots, args.backend,
            time_limit=args.time_limit,
        )
    elif args.action == "tune":
        receipt = _parse_receipt(args.receipt)
        bkts = receipt.get("braid", receipt).get("strands", [])
        if args.coeffs:
            custom_grid = [json.loads(Path(args.coeffs).read_text())]
            result = tune_braid_parameters(bkts, param_grid=custom_grid)
        else:
            result = tune_braid_parameters(bkts, time_per_trial=args.time_limit)
    elif args.action == "auto-solve":
        receipt = _parse_receipt(args.receipt)
        bkts = receipt.get("braid", receipt).get("strands", [])
        Q = build_tunable_qubo_matrix(bkts)
        qubo = QUBO(n=len(bkts), matrix=Q)
        result = solver_dispatch_tuner(qubo, time_limit=args.time_limit)
    elif args.action == "slos-calibrate":
        receipt = _parse_receipt(args.receipt)
        bkts = receipt.get("braid", receipt).get("strands", [])
        evaluator = SlosBraidEvaluator()
        result = slos_calibrate_coeffs(bkts, evaluator=evaluator)
        result["n"] = len(bkts)
        result["note"] = (
            "SLOS-calibrated coefficients. Set as braid_search defaults"
            " by copying calibrated_coeffs into BraidCostCoeffs() args."
        )
    elif args.action == "auto-calibrate":
        receipt = _parse_receipt(args.receipt)
        bkts = receipt.get("braid", receipt).get("strands", [])
        evaluator = SlosBraidEvaluator()
        result = auto_calibrate(bkts, evaluator=evaluator, time_per_solver=args.time_limit)
    elif args.action == "roundtrip":
        receipt = _parse_receipt(args.receipt)
        result = qaoa_roundtrip(
            receipt, args.p_layers, args.shots, args.backend,
            compare_stochastic=True,
        )
    elif args.action == "finsler-demo":
        result = finsler_demo()

    output = json.dumps(_serialize(result), indent=2, default=str)

    if args.output:
        Path(args.output).write_text(output)
    else:
        print(output)


if __name__ == "__main__":
    main()