SilverSight/python/bawim_mutation_engine.py
allaun 03fe2b6c4e feat(bawim): parameter wiring + cached J matrix + grid scan results
Key changes:
  - Wired j_max, h_bias_max into Sudoku evaluation (params now affect result)
  - Cached _SUDOKU_J_BASE (729x729 rebuilt only once)
  - Widened mutation ranges for h_bias_max (0.01→100) and bath_amp (0.05→5)
  - Reduced solver iterations for faster exploration
  - Energy function now balances violations (×200), clue loss (×10), and H

Grid scan results across h_bias_max × bath_oscillation_amplitude:
  - 36 configs tested, 2 Pareto points found
  - h_bias=1: 81 violations, 0 clues (too weak)
  - h_bias=5: 285 violations, 1 clue (too strong — forces bad local minima)
  - h_bias≥10: 285 violations, 0 clues (bias dominates, solver can't explore)
  - bath_amp has negligible effect at this solver depth
2026-07-07 09:59:30 -05:00

1129 lines
42 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

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