diff --git a/python/bawim_mutation_engine.py b/python/bawim_mutation_engine.py index fe4427ff..9644ba29 100644 --- a/python/bawim_mutation_engine.py +++ b/python/bawim_mutation_engine.py @@ -5,20 +5,26 @@ 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. -Papers equations: +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 100 --problem maxcut - python3 python/bawim_mutation_engine.py --rounds 100 --problem npp + 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 @@ -48,7 +54,7 @@ from q16_fraction import from_float as q16_from_float # BAWIM hardware constraints (Eq.2 context: 15-bit J_ij, 5-30% amplitude) BAWIM_DEFAULTS = { - "n_spins": 20, # N + "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 @@ -56,6 +62,13 @@ BAWIM_DEFAULTS = { "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π } @@ -74,6 +87,13 @@ class BawimParams: 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: @@ -86,6 +106,13 @@ class BawimParams: 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, ) @@ -100,6 +127,8 @@ class BawimParams: "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] @@ -171,7 +200,7 @@ def npp_spin_form( # ═══════════════════════════════════════════════════════════════════════════ -# §3 SPIN CONFIGURATION +# §3 SPIN CONFIGURATION + INCREMENTAL ΔH # ═══════════════════════════════════════════════════════════════════════════ def random_spins(n: int, seed: int) -> list[int]: @@ -185,8 +214,30 @@ def flip_spin(spins: list[int], idx: int) -> list[int]: 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 (mutation targets) +# §4 COUPLING MATRIX GENERATORS # ═══════════════════════════════════════════════════════════════════════════ def build_maxcut_graph( @@ -195,11 +246,7 @@ def build_maxcut_graph( seed: int, j_max: Fraction, ) -> list[list[Fraction]]: - """Generate a MAX-CUT graph coupling matrix. - - density is the fraction of edges present. - Each edge gets a random coupling in [0, j_max]. - """ + """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) @@ -215,11 +262,7 @@ def build_maxcut_graph( def build_npp_matrix( values: list[Fraction], ) -> tuple[list[list[Fraction]], list[Fraction]]: - """Build coupling matrix for number partitioning. - - J_ij = a_i * a_j (Mattis spin glass form) - h = 0 - """ + """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): @@ -232,38 +275,370 @@ def build_npp_matrix( # ═══════════════════════════════════════════════════════════════════════════ -# §5 SOLVER (greedy descent — mutation target) +# §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 = 100, + max_flips: int = 500, ) -> tuple[list[int], Fraction]: - """Greedy energy minimization by single-spin flips. - - Each flip that lowers H is accepted. Stops when no flip improves. - — MUTATION TARGET: replace with simulated annealing, QAOA, etc. - """ - current = list(spins) - H = ising_hamiltonian(current, J, h) - improved = True - n = len(current) - while improved and max_flips > 0: - improved = False + """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): - candidate = flip_spin(current, i) - H_candidate = ising_hamiltonian(candidate, J, h) - if H_candidate < H: - current = candidate - H = H_candidate - improved = True - max_flips -= 1 + 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 = 10, + 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 # ═══════════════════════════════════════════════════════════════════════════ @@ -271,44 +646,44 @@ def greedy_descent( @dataclass class Mutation: """A single mutation: which parameter changed, from what to what.""" - target: str # parameter name + target: str old_value: Any new_value: Any - delta_energy: Optional[Fraction] = None # improvement (negative = better) + delta_energy: Optional[Fraction] = None MUTATION_TARGETS = [ - # n_spins is excluded: changing N changes the graph, making - # cross-round energy comparisons meaningless. "j_resolution_bits", "j_max", "feedback_min_pct", "feedback_max_pct", "h_bias_max", "encoding", - "solver", # swap solver strategy + "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. - - Returns (new_params, description_of_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]) - new = max(4, min(200, p.n_spins + delta)) - p.n_spins = new + p.n_spins = max(4, min(200, p.n_spins + delta)) elif target == "j_resolution_bits": delta = rng.choice([-2, -1, 1, 2]) - new = max(4, min(20, p.j_resolution_bits + delta)) - p.j_resolution_bits = new - p.j_max = Fraction(2**new, 1) + 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), @@ -331,7 +706,7 @@ def mutate_param(params: BawimParams, rng: random.Random) -> tuple[BawimParams, elif target == "h_bias_max": scale = rng.choice([Fraction(1, 2), Fraction(3, 4), Fraction(1, 1), - Fraction(3, 2), Fraction(2, 1)]) + Fraction(5, 4), Fraction(3, 2), Fraction(2, 1)]) p.h_bias_max = max(Fraction(1, 100), min(Fraction(10, 1), (p.h_bias_max * scale).limit_denominator(100))) @@ -339,8 +714,41 @@ def mutate_param(params: BawimParams, rng: random.Random) -> tuple[BawimParams, p.encoding = rng.choice(["one_hot", "binary", "hybrid"]) elif target == "solver": - # Future: swap between greedy_descent, simulated_annealing, QAOA - pass + 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(1, 2), + p.bath_oscillation_amplitude * rng.choice([Fraction(1, 2), Fraction(3, 4), + Fraction(1, 1), Fraction(3, 2)]))) + + 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) @@ -350,6 +758,90 @@ def mutate_param(params: BawimParams, rng: random.Random) -> tuple[BawimParams, # §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=2000, + 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=10, + amplitude=params.bath_oscillation_amplitude, + frequency=params.bath_oscillation_frequency, + seed=params.seed, + ) + else: + return greedy_descent(spins, J, h, max_flips=500) + + +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. + """ + # Standard hard Sudoku puzzle (Wikipedia) + puzzle_str = ( + "530070000" + "600195000" + "098000060" + "800060003" + "400803001" + "700020006" + "060000280" + "000419005" + "000080079" + ) + clues = parse_sudoku_puzzle(puzzle_str) + + J = build_sudoku_couplings(penalty=Fraction(10, 1)) + h = build_sudoku_biases(clues, h_clue=Fraction(50, 1), h_anti=Fraction(10, 1)) + + 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.""" rng = random.Random(params.seed) @@ -361,14 +853,21 @@ def evaluate_maxcut(params: BawimParams) -> dict[str, Any]: ) h = [Fraction(0, 1)] * params.n_spins spins = random_spins(params.n_spins, params.seed) - final_spins, H = greedy_descent(spins, J, h) + 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, } @@ -379,14 +878,20 @@ def evaluate_npp(params: BawimParams) -> dict[str, Any]: 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 = greedy_descent(spins, J, h) + 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, } @@ -406,19 +911,11 @@ class MutationRound: class BawimMutationEngine: - """Run mutation rounds to improve BAWIM solution quality. - - Each round: - 1. Clone current params - 2. Mutate one parameter - 3. Evaluate on the problem - 4. Accept if energy improved (or with probability for exploration) - 5. Record the round - """ + """Run mutation rounds to improve BAWIM solution quality.""" def __init__( self, - problem: str = "maxcut", + problem: str = "sudoku", initial_params: Optional[BawimParams] = None, temperature: Fraction = Fraction(1, 10), ): @@ -431,7 +928,9 @@ class BawimMutationEngine: self.rng = random.Random(self.params.seed) def evaluate(self, params: BawimParams) -> dict[str, Any]: - if self.problem == "maxcut": + if self.problem == "sudoku": + return evaluate_sudoku(params) + elif self.problem == "maxcut": return evaluate_maxcut(params) elif self.problem == "npp": return evaluate_npp(params) @@ -439,10 +938,15 @@ class BawimMutationEngine: raise ValueError(f"Unknown problem: {self.problem}") def energy_key(self, result: dict[str, Any]) -> Fraction: - if self.problem == "maxcut": - return -result["cuts"] # more cuts = better (lower energy) + if self.problem == "sudoku": + # Lower violations + lower H = better. + # Strongly weight violations (each violation adds 100 to energy) + v = Fraction(result["total_violations"], 1) + return v * 100 + result["H"] + elif self.problem == "maxcut": + return -result["cuts"] else: - return result["E"] # lower E = better + return result["E"] def step(self) -> MutationRound: new_params, mutation = mutate_param(self.params, self.rng) @@ -453,12 +957,10 @@ class BawimMutationEngine: if self.best_energy is None: accepted = True elif new_energy < self.best_energy: - accepted = True # strict improvement + accepted = True else: - # Probabilistic acceptance (simulated annealing style) delta = new_energy - self.best_energy if self.temperature > 0: - # Convert Fraction to float for exp (I/O boundary — display only) p = math.exp(-float(delta) / float(self.temperature)) accepted = self.rng.random() < p @@ -490,7 +992,7 @@ class BawimMutationEngine: energy = float(self.energy_key(record.result)) print( f" [{r:4d}] {sign} " - f"{record.mutation.target:20s} " + f"{record.mutation.target:25s} " f"Δ={float(delta):+.6f} " f"E={energy:.6f}" ) @@ -499,7 +1001,7 @@ class BawimMutationEngine: def summary(self) -> dict[str, Any]: best = self.best_result or {} accepted = sum(1 for h in self.history if h.accepted) - return { + s = { "problem": self.problem, "rounds": len(self.history), "accepted": accepted, @@ -511,12 +1013,33 @@ class BawimMutationEngine: "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" +) # ═══════════════════════════════════════════════════════════════════════════ @@ -528,15 +1051,35 @@ def main() -> int: 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=["maxcut", "npp"], default="maxcut") + 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("--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, @@ -545,7 +1088,7 @@ def main() -> int: print(f"BAWIM Mutation Engine — {args.problem.upper()}") print(f" Initial params: n={params.n_spins}, J_bits={params.j_resolution_bits}") - print(f" Temperature: {args.temperature}") + print(f" Solver: {params.solver}, Temperature: {args.temperature}") print(f" Running {args.rounds} rounds...\n") start = time.time() @@ -559,6 +1102,9 @@ def main() -> int: 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']}")