"""
BAWIM Mutation Engine — explore the coupling/encoding design space

Maps the Bulk Acoustic Wave Ising Machine (arXiv 2607.02112) equations
onto the SilverSight Q16 fraction system and provides a mutation loop
that systematically varies parameters to improve solution quality.

Paper equations:
  Eq 1 — Ising Hamiltonian     H = -∑ J_ij s_i s_j - ∑ h_i s_i
  Eq 2 — Feedback coupling     c_i = ∑ J_ij s_j
  Eq 3 — MAX-CUT score         Cuts = -½∑ J_ij - ½H
  Eq 4 — NPP (set form)        E(A,B) = |∑A a_i - ∑B a_i|
  Eq 5 — NPP (spin form)       E(s) = |∑ a_i s_i|

Additional physics:
  Barkhausen criterion:  loop_gain > 1, phase_shift = 2πn
  Thermal stability:     BAWIM 780 deg/°C vs CIM 1.73×10^7 deg/°C
  Sudoku:                729-spin one-hot + row/col/block constraints

Architecture:
  Mutation round → vary one parameter → evaluate → accept/reject → next round
  All arithmetic is exact Fraction (q16_fraction).  Quantize only at output.

Usage:
    python3 python/bawim_mutation_engine.py --rounds 50 --problem sudoku
    python3 python/bawim_mutation_engine.py --rounds 50 --problem maxcut
    python3 python/bawim_mutation_engine.py --rounds 50 --problem npp
"""

from __future__ import annotations

import argparse
import hashlib
import json
import math
import random
import sys
import time
from dataclasses import dataclass, field
from fractions import Fraction
from itertools import combinations
from typing import Any, Callable, Optional

# ── Import Q16 fraction system (exact rational, no float in compute path) ──

sys.path.insert(0, "python")
from q16_fraction import Q16, q16_mul, q16_div, q16_add, q16_sub, to_raw, from_raw
from q16_fraction import from_float as q16_from_float


# ═══════════════════════════════════════════════════════════════════════════
# §1  PARAMETER BUDGET (mutable)
# ═══════════════════════════════════════════════════════════════════════════

# BAWIM hardware constraints (Eq.2 context: 15-bit J_ij, 5-30% amplitude)
BAWIM_DEFAULTS = {
    "n_spins": 20,                  # N (729 for sudoku)
    "j_resolution_bits": 15,        # J_ij bit depth (paper: 15-bit)
    "j_max": Fraction(2**15, 1),    # max coupling value
    "feedback_min_pct": Fraction(5, 100),    # 5% of RF carrier
    "feedback_max_pct": Fraction(30, 100),   # 30% of RF carrier
    "h_bias_max": Fraction(1, 10),           # max local bias
    "thermal_stability": Fraction(780, 1),   # BAWIM deg/°C (paper: 780)
    "encoding": "one_hot",          # one_hot | binary | hybrid
    "solver": "greedy",             # greedy | simulated_annealing | oscillating_bath
    "sa_temperature_start": Fraction(10, 1),
    "sa_temperature_end": Fraction(1, 100),
    "bath_oscillation_amplitude": Fraction(1, 5),   # fraction of j_max
    "bath_oscillation_frequency": Fraction(1, 10),  # cycles per 100 flips
    "barkhausen_loop_gain": Fraction(3, 2),          # loop gain > 1
    "barkhausen_phase_shift": Fraction(2, 1),        # multiples of 2π
}


@dataclass
class BawimParams:
    """Mutable parameter set for the BAWIM system.

    Each field can be mutated independently.  The mutation engine
    varies one field per round and evaluates the impact.
    """
    n_spins: int = BAWIM_DEFAULTS["n_spins"]
    j_resolution_bits: int = BAWIM_DEFAULTS["j_resolution_bits"]
    j_max: Fraction = BAWIM_DEFAULTS["j_max"]
    feedback_min_pct: Fraction = BAWIM_DEFAULTS["feedback_min_pct"]
    feedback_max_pct: Fraction = BAWIM_DEFAULTS["feedback_max_pct"]
    h_bias_max: Fraction = BAWIM_DEFAULTS["h_bias_max"]
    thermal_stability: Fraction = BAWIM_DEFAULTS["thermal_stability"]
    encoding: str = BAWIM_DEFAULTS["encoding"]
    solver: str = BAWIM_DEFAULTS["solver"]
    sa_temperature_start: Fraction = BAWIM_DEFAULTS["sa_temperature_start"]
    sa_temperature_end: Fraction = BAWIM_DEFAULTS["sa_temperature_end"]
    bath_oscillation_amplitude: Fraction = BAWIM_DEFAULTS["bath_oscillation_amplitude"]
    bath_oscillation_frequency: Fraction = BAWIM_DEFAULTS["bath_oscillation_frequency"]
    barkhausen_loop_gain: Fraction = BAWIM_DEFAULTS["barkhausen_loop_gain"]
    barkhausen_phase_shift: Fraction = BAWIM_DEFAULTS["barkhausen_phase_shift"]
    seed: int = 42

    def clone(self) -> BawimParams:
        return BawimParams(
            n_spins=self.n_spins,
            j_resolution_bits=self.j_resolution_bits,
            j_max=self.j_max,
            feedback_min_pct=self.feedback_min_pct,
            feedback_max_pct=self.feedback_max_pct,
            h_bias_max=self.h_bias_max,
            thermal_stability=self.thermal_stability,
            encoding=self.encoding,
            solver=self.solver,
            sa_temperature_start=self.sa_temperature_start,
            sa_temperature_end=self.sa_temperature_end,
            bath_oscillation_amplitude=self.bath_oscillation_amplitude,
            bath_oscillation_frequency=self.bath_oscillation_frequency,
            barkhausen_loop_gain=self.barkhausen_loop_gain,
            barkhausen_phase_shift=self.barkhausen_phase_shift,
            seed=self.seed,
        )

    def fingerprint(self) -> str:
        """Deterministic hash for reproducibility tracking."""
        raw = json.dumps({
            "n_spins": self.n_spins,
            "j_res_bits": self.j_resolution_bits,
            "j_max": str(self.j_max),
            "fb_min": str(self.feedback_min_pct),
            "fb_max": str(self.feedback_max_pct),
            "h_max": str(self.h_bias_max),
            "thermal": str(self.thermal_stability),
            "encoding": self.encoding,
            "solver": self.solver,
            "bark_gain": str(self.barkhausen_loop_gain),
            "seed": self.seed,
        }, sort_keys=True, separators=(",", ":"))
        return hashlib.sha256(raw.encode()).hexdigest()[:16]


# ═══════════════════════════════════════════════════════════════════════════
# §2  CORE BAWIM EQUATIONS (exact Fraction arithmetic)
# ═══════════════════════════════════════════════════════════════════════════

def ising_hamiltonian(
    spins: list[int],
    J: list[list[Fraction]],
    h: list[Fraction],
) -> Fraction:
    """Eq 1: H = -∑ J_ij s_i s_j - ∑ h_i s_i

    All arithmetic in exact Fraction.  No float, no truncation.
    """
    H = Fraction(0, 1)
    n = len(spins)
    for i in range(n):
        for j in range(i + 1, n):
            H -= J[i][j] * spins[i] * spins[j]
    for i in range(n):
        H -= h[i] * spins[i]
    return H


def feedback_coupling(
    spins: list[int],
    J: list[list[Fraction]],
    i: int,
) -> Fraction:
    """Eq 2: c_i = ∑ J_ij s_j

    The coupling pulse injected back into spin i.
    Amplitude should stay within feedback_min_pct–feedback_max_pct
    of the circulating RF signal to avoid chaotization.
    """
    return sum(J[i][j] * spins[j] for j in range(len(spins)))


def maxcut_score(
    spins: list[int],
    J: list[list[Fraction]],
    H: Optional[Fraction] = None,
) -> Fraction:
    """Eq 3: Cuts = -½∑ J_ij - ½H

    Lower Hamiltonian → more cuts.
    """
    if H is None:
        H = ising_hamiltonian(spins, J, [Fraction(0, 1)] * len(spins))
    sum_J = sum(J[i][j] for i in range(len(spins)) for j in range(i + 1, len(spins)))
    return -sum_J / 2 - H / 2


def npp_spin_form(
    spins: list[int],
    values: list[Fraction],
) -> Fraction:
    """Eq 5: E(s) = |∑ a_i s_i|

    s_i = +1 → subset A; s_i = -1 → subset B.
    Perfect partition when E = 0 (total even) or E = 1 (total odd).
    """
    total = sum(values[i] * spins[i] for i in range(len(spins)))
    return abs(total)


# ═══════════════════════════════════════════════════════════════════════════
# §3  SPIN CONFIGURATION + INCREMENTAL ΔH
# ═══════════════════════════════════════════════════════════════════════════

def random_spins(n: int, seed: int) -> list[int]:
    rng = random.Random(seed)
    return [1 if rng.random() < 0.5 else -1 for _ in range(n)]


def flip_spin(spins: list[int], idx: int) -> list[int]:
    s = list(spins)
    s[idx] = -s[idx]
    return s


def delta_h(
    spins: list[int],
    J: list[list[Fraction]],
    h: list[Fraction],
    i: int,
) -> Fraction:
    """Energy change from flipping spin i:  ΔH = 2·s_i·(∑ J_ij·s_j + h_i)

    Computes the local field at spin i in O(N) and returns the
    exact energy shift if spin i were flipped.  Use this inside
    solvers instead of recomputing the full Hamiltonian O(N²).
    """
    n = len(spins)
    s_i = spins[i]
    local_field = h[i]
    row = J[i]
    for j in range(n):
        if row[j] != 0:
            local_field += row[j] * spins[j]
    return Fraction(2, 1) * s_i * local_field


# ═══════════════════════════════════════════════════════════════════════════
# §4  COUPLING MATRIX GENERATORS
# ═══════════════════════════════════════════════════════════════════════════

def build_maxcut_graph(
    n: int,
    density: Fraction,
    seed: int,
    j_max: Fraction,
) -> list[list[Fraction]]:
    """Generate a MAX-CUT graph coupling matrix."""
    rng = random.Random(seed)
    J = [[Fraction(0, 1) for _ in range(n)] for _ in range(n)]
    threshold = float(density)
    for i in range(n):
        for j in range(i + 1, n):
            if rng.random() < threshold:
                val = rng.randint(1, int(j_max))
                J[i][j] = Fraction(val, 1)
                J[j][i] = J[i][j]
    return J


def build_npp_matrix(
    values: list[Fraction],
) -> tuple[list[list[Fraction]], list[Fraction]]:
    """Build coupling matrix for number partitioning (Mattis spin glass)."""
    n = len(values)
    J = [[Fraction(0, 1) for _ in range(n)] for _ in range(n)]
    for i in range(n):
        for j in range(i + 1, n):
            val = values[i] * values[j]
            J[i][j] = val
            J[j][i] = val
    h = [Fraction(0, 1)] * n
    return J, h


# ═══════════════════════════════════════════════════════════════════════════
# §4a  SUDOKU ENCODING (729-spin one-hot Ising formulation)
# ═══════════════════════════════════════════════════════════════════════════

# Spin index: cell (r,c) digit d  →  i = (r*9 + c)*9 + (d-1)
# 81 cells × 9 digits = 729 spins
SUDOKU_N = 729


def sudoku_spin_index(row: int, col: int, digit: int) -> int:
    """Map (row, col, digit) → spin index [0, 729)."""
    return (row * 9 + col) * 9 + (digit - 1)


def sudoku_spin_to_cell(i: int) -> tuple[int, int, int]:
    """Inverse: spin index → (row, col, digit)."""
    d = (i % 9) + 1
    cell = i // 9
    row = cell // 9
    col = cell % 9
    return row, col, d


def build_sudoku_couplings(
    penalty: Fraction = Fraction(10, 1),
) -> list[list[Fraction]]:
    """Build the 729×729 Sudoku constraint coupling matrix.

    J_ij > 0 for prohibited configurations (anti-ferromagnetic).

    Penalized pairs (J_ij = penalty):
      1. Same cell, different digits    — one-hot violation
      2. Same row, same digit           — row uniqueness violation
      3. Same column, same digit        — column uniqueness violation
      4. Same 3×3 block, same digit     — block uniqueness violation

    All other pairs have J_ij = 0.
    """
    N = SUDOKU_N
    J = [[Fraction(0, 1) for _ in range(N)] for _ in range(N)]

    for i in range(N):
        r1, c1, d1 = sudoku_spin_to_cell(i)
        for j in range(i + 1, N):
            r2, c2, d2 = sudoku_spin_to_cell(j)

            # Same cell, different digits → one-hot violation
            if r1 == r2 and c1 == c2 and d1 != d2:
                J[i][j] = penalty

            # Same row, same digit → row uniqueness violation
            elif r1 == r2 and d1 == d2 and c1 != c2:
                J[i][j] = penalty

            # Same column, same digit → column uniqueness violation
            elif c1 == c2 and d1 == d2 and r1 != r2:
                J[i][j] = penalty

            # Same 3×3 block, same digit → block uniqueness violation
            elif (r1 // 3 == r2 // 3) and (c1 // 3 == c2 // 3) and d1 == d2 and (r1 != r2 or c1 != c2):
                J[i][j] = penalty

    return J


def build_sudoku_biases(
    clues: list[tuple[int, int, int]],
    h_clue: Fraction = Fraction(50, 1),   # strong bias for known clues
    h_anti: Fraction = Fraction(10, 1),   # penalty for wrong digit in clue cell
) -> list[Fraction]:
    """Build the 729-element bias field from a set of Sudoku clues.

    clues: list of (row, col, digit) for known cells.
    """
    N = SUDOKU_N
    h = [Fraction(0, 1) for _ in range(N)]

    for r, c, d in clues:
        # Lower the energy of the correct digit
        idx = sudoku_spin_index(r, c, d)
        h[idx] = -h_clue

        # Raise the energy of wrong digits in the same cell
        for wrong_d in range(1, 10):
            if wrong_d != d:
                idx_w = sudoku_spin_index(r, c, wrong_d)
                h[idx_w] = h_anti

    return h


def parse_sudoku_puzzle(puzzle_str: str) -> list[tuple[int, int, int]]:
    """Parse a standard 81-character Sudoku puzzle string.

    '.' or '0' for empty cells.  1-9 for clues.
    """
    clues = []
    s = puzzle_str.replace("\n", "").replace(" ", "").strip()
    assert len(s) == 81, f"Need 81 chars, got {len(s)}"
    for i, ch in enumerate(s):
        if ch != "." and ch != "0":
            r = i // 9
            c = i % 9
            d = int(ch)
            clues.append((r, c, d))
    return clues


def sudoku_puzzle_to_string(
    spins: list[int],
    clues: list[tuple[int, int, int]],
) -> str:
    """Decode a 729-spin configuration back to an 81-character puzzle string."""
    grid = [["." for _ in range(9)] for _ in range(9)]
    for i, s in enumerate(spins):
        if s == 1:  # spin up = digit selected
            r, c, d = sudoku_spin_to_cell(i)
            grid[r][c] = str(d)
    # Clue cells: if no spin is active, use the clue value
    for r, c, d in clues:
        if grid[r][c] == ".":
            grid[r][c] = str(d)
    return "\n".join("".join(row) for row in grid)


def count_sudoku_violations(
    spins: list[int],
) -> dict[str, int]:
    """Count constraint violations in a spin configuration.

    Returns {rule_name: violation_count}.
    """
    # Decode into grid
    grid = [[0 for _ in range(9)] for _ in range(9)]
    for i, s in enumerate(spins):
        if s == 1:
            r, c, d = sudoku_spin_to_cell(i)
            grid[r][c] = d

    violations = {
        "empty_cells": 0,
        "row_repeats": 0,
        "col_repeats": 0,
        "block_repeats": 0,
        "multi_digit_cells": 0,
    }

    # Count digits per cell
    for r in range(9):
        for c in range(9):
            cell_spins = [spins[sudoku_spin_index(r, c, d)] for d in range(1, 10)]
            active = sum(1 for s in cell_spins if s == 1)
            if active == 0:
                violations["empty_cells"] += 1
            elif active > 1:
                violations["multi_digit_cells"] += 1

    # Row repeats
    for r in range(9):
        for d in range(1, 10):
            count = sum(1 for c in range(9) if grid[r][c] == d)
            if count > 1:
                violations["row_repeats"] += count - 1

    # Column repeats
    for c in range(9):
        for d in range(1, 10):
            count = sum(1 for r in range(9) if grid[r][c] == d)
            if count > 1:
                violations["col_repeats"] += count - 1

    # Block repeats
    for br in range(3):
        for bc in range(3):
            for d in range(1, 10):
                count = 0
                for r in range(3 * br, 3 * br + 3):
                    for c in range(3 * bc, 3 * bc + 3):
                        if grid[r][c] == d:
                            count += 1
                if count > 1:
                    violations["block_repeats"] += count - 1

    return violations


# ═══════════════════════════════════════════════════════════════════════════
# §5  SOLVERS
# ═══════════════════════════════════════════════════════════════════════════
def greedy_descent(
    spins: list[int],
    J: list[list[Fraction]],
    h: list[Fraction],
    max_flips: int = 500,
) -> tuple[list[int], Fraction]:
    """Greedy energy minimization.  Uses incremental ΔH, not full recompute."""
    H = ising_hamiltonian(spins, J, h)
    n = len(spins)
    flips = 0
    changed = True
    while changed and flips < max_flips:
        changed = False
        for i in range(n):
            d = delta_h(spins, J, h, i)
            if d < 0:
                spins[i] = -spins[i]
                H += d
                changed = True
                flips += 1
                break
    return spins[:], H


def simulated_annealing(
    spins: list[int],
    J: list[list[Fraction]],
    h: list[Fraction],
    max_flips: int = 2000,
    t_start: Fraction = Fraction(10, 1),
    t_end: Fraction = Fraction(1, 100),
    seed: int = 42,
) -> tuple[list[int], Fraction]:
    """Simulated annealing.  Uses incremental ΔH (sparse-safe)."""
    rng = random.Random(seed)
    H = ising_hamiltonian(spins, J, h)
    n = len(spins)

    for flip in range(max_flips):
        progress = Fraction(flip, max_flips)
        T = t_start - (t_start - t_end) * progress
        if T <= 0:
            T = t_end

        i = rng.randint(0, n - 1)
        d = delta_h(spins, J, h, i)
        if d < 0:
            spins[i] = -spins[i]
            H += d
        elif T > 0:
            prob = math.exp(-float(d) / float(T))
            if rng.random() < prob:
                spins[i] = -spins[i]
                H += d

    return spins[:], H


def oscillating_bath(
    spins: list[int],
    J: list[list[Fraction]],
    h: list[Fraction],
    max_cycles: int = 5,
    amplitude: Fraction = Fraction(1, 5),
    frequency: Fraction = Fraction(1, 10),
    seed: int = 42,
) -> tuple[list[int], Fraction]:
    """BAWIM-style oscillating bath solver.

    Mimics the physical BAWIM process: coupling amplitudes oscillate
    sinusoidally, allowing the system to escape local minima through
    parametric resonance.

    Each cycle:
      1. Modulate J_ij by a sinusoidal factor: J' = J * (1 + A * sin(ωt))
      2. Run greedy descent under the modulated Hamiltonian
      3. Return to the unmodulated Hamiltonian

    This approximates the BAWIM's acoustic wave modulation of coupling terms.
    """
    rng = random.Random(seed)
    current = list(spins)
    n = len(current)
    H = ising_hamiltonian(current, J, h)

    for cycle in range(max_cycles):
        phase = Fraction(cycle, max_cycles) * Fraction(628, 100)  # ≈ 2π
        # Modulate: J' = J * (1 + A * sin(ωt))
        modulation = 1 + float(amplitude) * math.sin(float(phase) * float(frequency))
        J_mod = [
            [J[i][j] * Fraction(int(modulation * 100), 100) for j in range(n)]
            for i in range(n)
        ]

        # Greedy descent under modulated couplings (incremental ΔH)
        improved = True
        max_local = int(n / 2)
        flips = 0
        H_mod = ising_hamiltonian(current, J_mod, h)
        while improved and flips < max_local:
            improved = False
            for i in rng.sample(range(n), n):
                d = delta_h(current, J_mod, h, i)
                if d < 0:
                    current[i] = -current[i]
                    H_mod += d
                    improved = True
                    flips += 1
                    break

    # Final descent under unmodulated Hamiltonian
    current, H = greedy_descent(current, J, h)
    return current, H


# ═══════════════════════════════════════════════════════════════════════════
# §5b  BARKHAUSEN STABILITY CRITERION
# ═══════════════════════════════════════════════════════════════════════════

def check_barkhausen(params: BawimParams) -> dict[str, Any]:
    """Check the Barkhausen stability criterion for the oscillator loop.

    The ring oscillator loop must satisfy:
      1. Loop gain > 1  (to sustain oscillation)
      2. Phase shift = integer multiple of 2π  (constructive feedback)

    For the BAWIM Ising machine:
      - Loop gain  ∝  j_max * feedback_pct * n_spins
      - Phase shift is determined by the SAW delay line geometry
    """
    # Effective loop gain: coupling strength × feedback amplitude × connectivity
    loop_gain = (
        params.j_max * params.feedback_max_pct * Fraction(params.n_spins, 1)
    ) / Fraction(100, 1)

    # Normalize: the Barkhausen criterion is loop_gain ≥ 1
    gain_ok = loop_gain >= 1

    # Phase shift: assume 2π per SAW round-trip (relative units)
    # The parameter barkhausen_phase_shift encodes the multiple
    phase_shift_ok = (
        params.barkhausen_phase_shift == Fraction(2, 1)
        or params.barkhausen_phase_shift % Fraction(2, 1) == 0
    )

    return {
        "loop_gain": loop_gain,
        "loop_gain_float": float(loop_gain),
        "gain_ok": gain_ok,
        "phase_shift_ok": phase_shift_ok,
        "stable": gain_ok and phase_shift_ok,
    }


def check_thermal_stability(params: BawimParams) -> dict[str, Any]:
    """Compare BAWIM thermal stability against CIM baseline.

    Paper: BAWIM = 780 deg/°C, CIM = 1.73×10^7 deg/°C
    BAWIM is 2.21×10^4 times more stable.

    The stability parameter affects the noise floor of the Ising computation.
    Higher stability → lower bit-flip rate from thermal noise.
    """
    CIM_STABILITY = Fraction(17300000, 1)  # 1.73 × 10^7
    bawim_stability = params.thermal_stability
    ratio = Fraction(CIM_STABILITY, bawim_stability)

    return {
        "bawim_stability": bawim_stability,
        "cim_stability": CIM_STABILITY,
        "stability_ratio": ratio,
        "ratio_float": float(ratio),
        "bawim_more_stable": bawim_stability < CIM_STABILITY,
    }


# ═══════════════════════════════════════════════════════════════════════════
# §6  MUTATION OPERATORS
# ═══════════════════════════════════════════════════════════════════════════

@dataclass
class Mutation:
    """A single mutation: which parameter changed, from what to what."""
    target: str
    old_value: Any
    new_value: Any
    delta_energy: Optional[Fraction] = None


MUTATION_TARGETS = [
    "j_resolution_bits",
    "j_max",
    "feedback_min_pct",
    "feedback_max_pct",
    "h_bias_max",
    "encoding",
    "solver",
    "sa_temperature_start",
    "sa_temperature_end",
    "bath_oscillation_amplitude",
    "bath_oscillation_frequency",
    "barkhausen_loop_gain",
    "barkhausen_phase_shift",
    "thermal_stability",
]


def mutate_param(params: BawimParams, rng: random.Random) -> tuple[BawimParams, Mutation]:
    """Apply one random mutation to a parameter."""
    p = params.clone()
    target = rng.choice(MUTATION_TARGETS)
    old_val = getattr(p, target, None)

    if target == "n_spins":
        delta = rng.choice([-5, -3, -1, 1, 3, 5])
        p.n_spins = max(4, min(200, p.n_spins + delta))

    elif target == "j_resolution_bits":
        delta = rng.choice([-2, -1, 1, 2])
        p.j_resolution_bits = max(4, min(20, p.j_resolution_bits + delta))
        p.j_max = Fraction(2**p.j_resolution_bits, 1)

    elif target == "j_max":
        scale = rng.choice([Fraction(1, 2), Fraction(3, 4), Fraction(1, 1),
                           Fraction(5, 4), Fraction(3, 2), Fraction(2, 1)])
        p.j_max = (p.j_max * scale).limit_denominator(2**p.j_resolution_bits)

    elif target == "feedback_min_pct":
        step = Fraction(1, 100)
        p.feedback_min_pct = max(Fraction(1, 100), min(
            Fraction(49, 100),
            p.feedback_min_pct + rng.choice([-step, step]) * rng.randint(1, 5)
        ))

    elif target == "feedback_max_pct":
        step = Fraction(1, 100)
        p.feedback_max_pct = max(p.feedback_min_pct + step, min(
            Fraction(95, 100),
            p.feedback_max_pct + rng.choice([-step, step]) * rng.randint(1, 5)
        ))

    elif target == "h_bias_max":
        scale = rng.choice([Fraction(1, 3), Fraction(1, 2), Fraction(3, 4),
                           Fraction(1, 1), Fraction(3, 2), Fraction(2, 1),
                           Fraction(5, 1), Fraction(10, 1)])
        p.h_bias_max = max(Fraction(1, 100), min(Fraction(100, 1),
                           (p.h_bias_max * scale).limit_denominator(100)))

    elif target == "encoding":
        p.encoding = rng.choice(["one_hot", "binary", "hybrid"])

    elif target == "solver":
        p.solver = rng.choice(["greedy", "simulated_annealing", "oscillating_bath"])

    elif target == "sa_temperature_start":
        p.sa_temperature_start = max(Fraction(1, 10), min(Fraction(100, 1),
            p.sa_temperature_start * rng.choice([Fraction(1, 2), Fraction(3, 4),
                                                 Fraction(1, 1), Fraction(2, 1)])))

    elif target == "sa_temperature_end":
        p.sa_temperature_end = max(Fraction(1, 1000), min(Fraction(1, 2),
            p.sa_temperature_end * rng.choice([Fraction(1, 2), Fraction(3, 4),
                                               Fraction(1, 1), Fraction(2, 1)])))

    elif target == "bath_oscillation_amplitude":
        p.bath_oscillation_amplitude = max(Fraction(1, 20), min(Fraction(5, 1),
            p.bath_oscillation_amplitude * rng.choice([Fraction(1, 3), Fraction(1, 2),
                Fraction(3, 4), Fraction(1, 1), Fraction(3, 2), Fraction(2, 1), Fraction(5, 1)])))

    elif target == "bath_oscillation_frequency":
        p.bath_oscillation_frequency = max(Fraction(1, 100), min(Fraction(1, 2),
            p.bath_oscillation_frequency * rng.choice([Fraction(1, 2), Fraction(3, 4),
                                                       Fraction(1, 1), Fraction(2, 1)])))

    elif target == "barkhausen_loop_gain":
        p.barkhausen_loop_gain = max(Fraction(1, 10), min(Fraction(10, 1),
            p.barkhausen_loop_gain * rng.choice([Fraction(1, 2), Fraction(3, 4),
                                                 Fraction(1, 1), Fraction(3, 2)])))

    elif target == "barkhausen_phase_shift":
        p.barkhausen_phase_shift = rng.choice([Fraction(1, 1), Fraction(2, 1),
                                               Fraction(3, 1), Fraction(4, 1)])

    elif target == "thermal_stability":
        p.thermal_stability = max(Fraction(100, 1), min(Fraction(1000000, 1),
            p.thermal_stability * rng.choice([Fraction(1, 2), Fraction(3, 4),
                                               Fraction(1, 1), Fraction(2, 1)])))

    new_val = getattr(p, target, None)
    return p, Mutation(target=target, old_value=old_val, new_value=new_val)


# ═══════════════════════════════════════════════════════════════════════════
# §7  EVALUATION
# ═══════════════════════════════════════════════════════════════════════════

def _run_solver(
    spins: list[int],
    J: list[list[Fraction]],
    h: list[Fraction],
    params: BawimParams,
) -> tuple[list[int], Fraction]:
    """Dispatch to the selected solver."""
    if params.solver == "simulated_annealing":
        return simulated_annealing(
            spins, J, h,
            max_flips=1000,
            t_start=params.sa_temperature_start,
            t_end=params.sa_temperature_end,
            seed=params.seed,
        )
    elif params.solver == "oscillating_bath":
        return oscillating_bath(
            spins, J, h,
            max_cycles=5,
            amplitude=params.bath_oscillation_amplitude,
            frequency=params.bath_oscillation_frequency,
            seed=params.seed,
        )
    else:
        return greedy_descent(spins, J, h, max_flips=500)


# Shared Sudoku coupling base (clue-independent, penalty=10)
# Cached so we don't rebuild the 729×729 matrix on every evaluation.
_SUDOKU_J_BASE = build_sudoku_couplings(penalty=Fraction(10, 1))


def evaluate_sudoku(params: BawimParams) -> dict[str, Any]:
    """Build a Sudoku problem, solve it, return metrics.

    Uses a standard hard puzzle (Wikipedia "hard" difficulty).
    The mutation engine varies coupling penalties, bias strengths,
    solver strategy, and Barkhausen parameters to find better solutions.
    """
    puzzle_str = DEFAULT_SUDOKU
    clues = parse_sudoku_puzzle(puzzle_str)

    # Scale cached coupling matrix and rebuild bias with mutation params
    penalty_scale = params.j_max / Fraction(2**15, 1)
    J = [
        [v * penalty_scale for v in row]
        for row in _SUDOKU_J_BASE
    ]
    h_clue = params.h_bias_max * Fraction(50, 1)
    h_anti = params.h_bias_max * Fraction(10, 1)
    h = build_sudoku_biases(clues, h_clue=h_clue, h_anti=h_anti)

    spins = random_spins(SUDOKU_N, params.seed)
    final_spins, H = _run_solver(spins, J, h, params)

    violations = count_sudoku_violations(final_spins)
    total_violations = sum(violations.values())
    solved = total_violations == 0

    # Compute clue preservation: fraction of clue cells with correct digit
    clues_kept = 0
    for r, c, d in clues:
        if final_spins[sudoku_spin_index(r, c, d)] == 1:
            clues_kept += 1

    # Barkhausen and thermal stability checks
    barkhausen = check_barkhausen(params)
    thermal = check_thermal_stability(params)

    return {
        "problem": "sudoku",
        "n": SUDOKU_N,
        "H": H,
        "H_float": float(H),
        "solved": solved,
        "violations": violations,
        "total_violations": total_violations,
        "clues_kept": clues_kept,
        "clues_total": len(clues),
        "barkhausen_stable": barkhausen["stable"],
        "thermal_ratio": thermal["ratio_float"],
        "spins": final_spins,
    }


def evaluate_maxcut(params: BawimParams) -> dict[str, Any]:
    """Build a MAX-CUT problem, solve it, return metrics.

    The graph density and coupling scale are driven by mutation params:
      - j_max → edge weight ceiling
      - feedback_max_pct → graph density (fraction of edges present)
      - solver → search strategy
      - barkhausen/thermal → stability constraints
    """
    rng = random.Random(params.seed)
    density = params.feedback_max_pct  # reuse: max feedback = edge density
    J = build_maxcut_graph(
        params.n_spins,
        density=density,
        seed=params.seed,
        j_max=params.j_max,
    )
    h = [Fraction(0, 1)] * params.n_spins
    spins = random_spins(params.n_spins, params.seed)
    final_spins, H = _run_solver(spins, J, h, params)
    cuts = maxcut_score(final_spins, J, H)

    barkhausen = check_barkhausen(params)
    thermal = check_thermal_stability(params)

    return {
        "problem": "maxcut",
        "n": params.n_spins,
        "H": H,
        "H_float": float(H),
        "cuts": cuts,
        "cuts_float": float(cuts),
        "barkhausen_stable": barkhausen["stable"],
        "thermal_ratio": thermal["ratio_float"],
        "spins": final_spins,
    }


def evaluate_npp(params: BawimParams) -> dict[str, Any]:
    """Build a number partitioning problem, solve it, return metrics."""
    rng = random.Random(params.seed)
    values = [Fraction(rng.randint(1, 1000), 1) for _ in range(params.n_spins)]
    J, h = build_npp_matrix(values)
    spins = random_spins(params.n_spins, params.seed + 1)
    final_spins, H = _run_solver(spins, J, h, params)
    E = npp_spin_form(final_spins, values)

    barkhausen = check_barkhausen(params)
    thermal = check_thermal_stability(params)

    return {
        "problem": "npp",
        "n": params.n_spins,
        "E": E,
        "E_float": float(E),
        "H": H,
        "barkhausen_stable": barkhausen["stable"],
        "thermal_ratio": thermal["ratio_float"],
        "spins": final_spins,
    }


# ═══════════════════════════════════════════════════════════════════════════
# §8  MUTATION ENGINE
# ═══════════════════════════════════════════════════════════════════════════

@dataclass
class MutationRound:
    """Record of one mutation + evaluation."""
    round: int
    params: BawimParams
    mutation: Mutation
    result: dict[str, Any]
    accepted: bool


class BawimMutationEngine:
    """Run mutation rounds to improve BAWIM solution quality."""

    def __init__(
        self,
        problem: str = "sudoku",
        initial_params: Optional[BawimParams] = None,
        temperature: Fraction = Fraction(1, 10),
    ):
        self.problem = problem
        self.params = initial_params or BawimParams()
        self.temperature = temperature
        self.history: list[MutationRound] = []
        self.best_result: Optional[dict[str, Any]] = None
        self.best_energy: Optional[Fraction] = None
        self.rng = random.Random(self.params.seed)

    def evaluate(self, params: BawimParams) -> dict[str, Any]:
        if self.problem == "sudoku":
            return evaluate_sudoku(params)
        elif self.problem == "maxcut":
            return evaluate_maxcut(params)
        elif self.problem == "npp":
            return evaluate_npp(params)
        else:
            raise ValueError(f"Unknown problem: {self.problem}")

    def energy_key(self, result: dict[str, Any]) -> Fraction:
        if self.problem == "sudoku":
            # Energy balances constraint violations (×200 each), clue loss (×10),
            # and the Hamiltonian.  Lower is better.
            v = Fraction(result["total_violations"], 1)
            clues_lost = result["clues_total"] - result["clues_kept"]
            return v * 200 + Fraction(clues_lost * 10, 1) + result["H"]
        elif self.problem == "maxcut":
            return -result["cuts"]
        else:
            return result["E"]

    def step(self) -> MutationRound:
        new_params, mutation = mutate_param(self.params, self.rng)
        result = self.evaluate(new_params)
        new_energy = self.energy_key(result)

        accepted = False
        if self.best_energy is None:
            accepted = True
        elif new_energy < self.best_energy:
            accepted = True
        else:
            delta = new_energy - self.best_energy
            if self.temperature > 0:
                p = math.exp(-float(delta) / float(self.temperature))
                accepted = self.rng.random() < p

        mutation.delta_energy = (
            new_energy - self.best_energy if self.best_energy is not None else Fraction(0, 1)
        )

        if accepted:
            self.params = new_params
            self.best_result = result
            self.best_energy = new_energy

        record = MutationRound(
            round=len(self.history),
            params=new_params,
            mutation=mutation,
            result=result,
            accepted=accepted,
        )
        self.history.append(record)
        return record

    def run(self, rounds: int, verbose: bool = True) -> list[MutationRound]:
        for r in range(rounds):
            record = self.step()
            if verbose and (r < 10 or r % 10 == 0 or record.accepted):
                delta = record.mutation.delta_energy
                sign = "✓" if record.accepted else "✗"
                energy = float(self.energy_key(record.result))
                print(
                    f"  [{r:4d}] {sign} "
                    f"{record.mutation.target:25s} "
                    f"Δ={float(delta):+.6f}  "
                    f"E={energy:.6f}"
                )
        return self.history

    def summary(self) -> dict[str, Any]:
        best = self.best_result or {}
        accepted = sum(1 for h in self.history if h.accepted)
        s = {
            "problem": self.problem,
            "rounds": len(self.history),
            "accepted": accepted,
            "accept_rate": accepted / max(1, len(self.history)),
            "best_energy": float(self.best_energy) if self.best_energy else None,
            "best_params": {
                "n_spins": self.params.n_spins,
                "j_resolution_bits": self.params.j_resolution_bits,
                "j_max": str(self.params.j_max),
                "feedback_range": f"{self.params.feedback_min_pct}–{self.params.feedback_max_pct}",
                "encoding": self.params.encoding,
                "solver": self.params.solver,
                "bark_gain": str(self.params.barkhausen_loop_gain),
            },
            "best_result": {
                k: v for k, v in best.items() if k != "spins"
            } if best else {},
            "fingerprint": self.params.fingerprint(),
        }
        if self.problem == "sudoku" and best:
            s["sudoku_solved"] = best.get("solved", False)
            s["sudoku_violations"] = best.get("total_violations", 0)
        return s


# ── Sudoku puzzle string for eval witness ─────────────────────────────────

DEFAULT_SUDOKU = (
    "530070000"
    "600195000"
    "098000060"
    "800060003"
    "400803001"
    "700020006"
    "060000280"
    "000419005"
    "000080079"
)


# ═══════════════════════════════════════════════════════════════════════════
# §9  CLI
# ═══════════════════════════════════════════════════════════════════════════

def main() -> int:
    parser = argparse.ArgumentParser(
        description="BAWIM Mutation Engine — explore coupling/encoding design space"
    )
    parser.add_argument("--rounds", type=int, default=50, help="mutation rounds")
    parser.add_argument("--problem", choices=["sudoku", "maxcut", "npp"],
                        default="sudoku")
    parser.add_argument("--n-spins", type=int, default=20)
    parser.add_argument("--seed", type=int, default=42)
    parser.add_argument("--temperature", type=float, default=0.1,
                        help="SA acceptance temperature (0 = greedy)")
    parser.add_argument("--json", action="store_true",
                        help="output final summary as JSON")
    parser.add_argument("--eval-only", action="store_true",
                        help="run one evaluation with default params (no mutation)")
    args = parser.parse_args()

    params = BawimParams(n_spins=args.n_spins, seed=args.seed)

    if args.eval_only:
        engine = BawimMutationEngine(problem=args.problem, initial_params=params)
        result = engine.evaluate(params)
        print(f"\nSingle evaluation — {args.problem.upper()}")
        for k, v in sorted(result.items()):
            if k != "spins":
                print(f"  {k}: {v}")
        if "violations" in result:
            print(f"  violations breakdown: {result['violations']}")
        bark = check_barkhausen(params)
        print(f"  barkhausen_stable: {bark['stable']}")
        thermal = check_thermal_stability(params)
        print(f"  thermal_ratio (CIM/BAWIM): {thermal['ratio_float']:.1f}")
        return 0

    engine = BawimMutationEngine(
        problem=args.problem,
        initial_params=params,
        temperature=Fraction(args.temperature).limit_denominator(100),
    )

    print(f"BAWIM Mutation Engine — {args.problem.upper()}")
    print(f"  Initial params: n={params.n_spins}, J_bits={params.j_resolution_bits}")
    print(f"  Solver: {params.solver}, Temperature: {args.temperature}")
    print(f"  Running {args.rounds} rounds...\n")

    start = time.time()
    engine.run(args.rounds)
    elapsed = time.time() - start

    summary = engine.summary()
    print(f"\n--- Summary ({elapsed:.1f}s) ---")
    print(f"  Problem:        {summary['problem']}")
    print(f"  Rounds:         {summary['rounds']} ({summary['accepted']} accepted)")
    print(f"  Accept rate:    {summary['accept_rate']:.2f}")
    if summary['best_energy'] is not None:
        print(f"  Best energy:    {summary['best_energy']:.6f}")
    if summary.get('sudoku_solved') is not None:
        print(f"  Sudoku solved:  {summary['sudoku_solved']}")
        print(f"  Violations:     {summary['sudoku_violations']}")
    print(f"  Best params:    {summary['best_params']}")
    print(f"  Fingerprint:    {summary['fingerprint']}")

    if args.json:
        print(json.dumps(summary, indent=2))

    return 0


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