""" 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())