Skip to content

Instantly share code, notes, and snippets.

@gukoff

gukoff/exchange_solver.py

Created August 2, 2026 17:24
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save gukoff/e62a157d51c86f7f2dab4dfa04aa0f9e to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save gukoff/e62a157d51c86f7f2dab4dfa04aa0f9e to your computer and use it in GitHub Desktop.
Download ZIP
Magic Hexagon search algorithm using simulated annealing. Optimized for performance
Raw
exchange_solver.py
#!/usr/bin/env python3
"""Permutation-preserving solver for balanced abnormal magic hexagons.
This solver works in the half-turn antisymmetric branch H(-c)=-H(c), H(0)=0.
It never leaves the consecutive shell: the state is always a signed permutation
of 1..m on antipodal cell pairs. The only objective is the 3(n-1) positive-line
residual vector.
The model is equivalent to three signed triangular arrays A,B,C indexed by
0 <= a < b <= R=n-1. See the accompanying research note for the formulas.
Typical use inside the project repository:
python exchange_solver.py \
--order 13 \
--hint known_solutions/MagicHexagon-Order12-sum_zero.mhx \
--allow-cross-order \
--workers 32 --until-found \
--steps 3000000 --finisher-energy 100 \
--finisher-size 72 --finisher-attempts 12 \
--out artifacts/order13_exchange.mhx \
--best-out artifacts/order13_exchange_best.mhx
"""
from __future__ import annotations
# R3_NUMBA_PERF_BOOTSTRAP_BEGIN
# Numba reads its profiling/debug environment during import. Re-exec before
# importing numpy/numba so --numba-perf works as a normal command-line switch.
import os
import sys
if (
"--numba-perf" in sys.argv
and os.environ.get("R3_NUMBA_PERF_ACTIVE") != "1"
):
_perf_env = os.environ.copy()
_perf_env["R3_NUMBA_PERF_ACTIVE"] = "1"
_perf_env["NUMBA_ENABLE_PROFILING"] = "1"
_perf_env["NUMBA_DEBUGINFO"] = "1"
# A unique cache directory forces a fresh debug/profiling compilation.
# Numba cannot expose object bytes for some cache-loaded specializations,
# which would prevent construction of the live perf map.
_perf_env.setdefault(
"NUMBA_CACHE_DIR",
f"/tmp/r3-numba-perf-cache-{os.getpid()}",
)
os.execvpe(sys.executable, [sys.executable, *sys.argv], _perf_env)
# R3_NUMBA_PERF_BOOTSTRAP_END
import argparse
import concurrent.futures
import math
import multiprocessing as mp
import os
import re
import time
from pathlib import Path
import numpy as np
from numba import njit
def cells(order: int) -> list[tuple[int, int, int]]:
R = order - 1
return [
(q, r, -q-r)
for r in range(-R, R+1)
for q in range(max(-R, -r-R), min(R, -r+R)+1)
]
def parse_mhx(path: Path) -> tuple[int, dict[tuple[int, int, int], int]]:
text = path.read_text(encoding="utf-8")
mo = re.search(r"^order:\s*(\d+)\s*$", text, re.MULTILINE)
if mo is None:
raise ValueError(f"{path}: missing order header")
n = int(mo.group(1))
if "rows:" not in text:
raise ValueError(f"{path}: missing rows section")
row_lines = [ln.strip() for ln in text.split("rows:", 1)[1].splitlines() if ln.strip()]
rows = [[int(x) for x in ln.split()] for ln in row_lines]
expected = [n+i for i in range(n)] + [n+i for i in range(n-2, -1, -1)]
if [len(row) for row in rows] != expected:
raise ValueError(f"{path}: wrong row shape")
R = n - 1
out: dict[tuple[int, int, int], int] = {}
for row_i, row in enumerate(rows):
r = row_i - R
q0 = max(-R, -r-R)
for j, value in enumerate(row):
q = q0 + j
out[(q, r, -q-r)] = value
return n, out
def representative_data(order: int):
cs = cells(order)
reps: list[tuple[int, int, int]] = []
for c in cs:
if c == (0, 0, 0):
continue
nc = (-c[0], -c[1], -c[2])
if c > nc:
reps.append(c)
rid = {c: i for i, c in enumerate(reps)}
R = order - 1
B = np.zeros((3*R, len(reps)), dtype=np.int8)
for i, c in enumerate(reps):
for axis in range(3):
z = c[axis]
if z > 0:
B[axis*R + z-1, i] = 1
elif z < 0:
B[axis*R + (-z)-1, i] = -1
return cs, reps, rid, B
def psi0(c: tuple[int, int, int]) -> int:
q, r, s = c
return (q-r)*(r-s)*(s-q)
def exact_rank_seed(order: int, reps: list[tuple[int, int, int]]) -> np.ndarray:
"""Consecutive seed from |psi0| ranks; ties are broken deterministically."""
m = len(reps)
raw = np.array([psi0(c) for c in reps], dtype=np.int64)
signs = np.where(raw >= 0, 1, -1)
# Coordinates break the large Weyl-orbit ties which make f(psi0) non-injective.
order_idx = sorted(
range(m),
key=lambda i: (abs(int(raw[i])), reps[i][0], reps[i][1], reps[i][2]),
)
mags = np.empty(m, dtype=np.int64)
mags[np.array(order_idx, dtype=int)] = np.arange(1, m+1, dtype=np.int64)
return signs*mags
def cross_order_seed(
order: int,
reps: list[tuple[int, int, int]],
hint_order: int,
hint_map: dict[tuple[int, int, int], int],
) -> np.ndarray:
R = order - 1
Rh = hint_order - 1
hint_cells = np.array(list(hint_map.keys()), dtype=float)
hint_vals = np.array([hint_map[tuple(map(int, c))] for c in hint_cells], dtype=np.int64)
hint_norm = hint_cells / max(1, Rh)
vals = np.empty(len(reps), dtype=np.int64)
for i, c in enumerate(reps):
p = np.array(c, dtype=float) / max(1, R)
j = int(np.argmin(np.sum((hint_norm-p)**2, axis=1)))
vals[i] = hint_vals[j]
signs = np.sign(vals)
for i in range(len(signs)):
if signs[i] == 0:
signs[i] = 1 if psi0(reps[i]) >= 0 else -1
order_idx = sorted(
range(len(reps)),
key=lambda i: (abs(int(vals[i])), reps[i][0], reps[i][1], reps[i][2]),
)
mags = np.empty(len(reps), dtype=np.int64)
mags[np.array(order_idx, dtype=int)] = np.arange(1, len(reps)+1, dtype=np.int64)
return signs.astype(np.int64)*mags
# R3_POTENTIAL_GUIDANCE_V1_BEGIN
# R3_POTENTIAL_GUIDANCE_PATCH_APPLIED
# Local-cycle ordering matches tools/hypothesis_probes.py and
# tools/visualize_mhx.py. The alternating stencil has Fourier symbol
# proportional to (e^{iq}-1)(e^{ir}-1)(e^{iq}-e^{ir}), so its inverse strongly
# amplifies broad terrain modes. We exploit that structure only at startup.
_POTENTIAL_CYCLE_OFFSETS = (
(0, -1, 1),
(1, -1, 0),
(-1, 0, 1),
(1, 0, -1),
(-1, 1, 0),
(0, 1, -1),
)
# D6 on cube coordinates is every coordinate permutation, optionally followed
# by global negation. Reflection can reverse the oriented-cycle coefficient;
# pairwise correlation alignment below absorbs that harmless global sign.
_POTENTIAL_D6_TRANSFORMS = (
((0, 1, 2), 1), ((0, 2, 1), 1),
((1, 0, 2), 1), ((1, 2, 0), 1),
((2, 0, 1), 1), ((2, 1, 0), 1),
((0, 1, 2), -1), ((0, 2, 1), -1),
((1, 0, 2), -1), ((1, 2, 0), -1),
((2, 0, 1), -1), ((2, 1, 0), -1),
)
def _potential_cycle_data(order: int):
"""Return full cells, valid cycle centers, and six cell indices per cycle."""
full_cells = cells(order)
index = {c: i for i, c in enumerate(full_cells)}
centers: list[tuple[int, int, int]] = []
cycle_rows: list[list[int]] = []
for c in full_cells:
neighbors = [
(c[0] + dq, c[1] + dr, c[2] + ds)
for dq, dr, ds in _POTENTIAL_CYCLE_OFFSETS
]
if all(nn in index for nn in neighbors):
centers.append(c)
cycle_rows.append([index[nn] for nn in neighbors])
return (
full_cells,
np.asarray(centers, dtype=np.int16),
np.asarray(cycle_rows, dtype=np.int32),
)
@njit(cache=True)
def _potential_recover_cg(
cycle_cells: np.ndarray,
target: np.ndarray,
tolerance: float,
max_iterations: int,
):
"""Solve C.T w = target through CG on C C.T, starting at zero.
C has two gauge dependencies. Starting in the range of C and applying CG
to the consistent positive-semidefinite normal system selects the
minimum-norm solution without forming a dense pseudoinverse. Every matrix
application touches exactly twelve scalar stencil entries per cycle.
"""
cycle_count = cycle_cells.shape[0]
cell_count = target.shape[0]
rhs = np.empty(cycle_count, dtype=np.float64)
for i in range(cycle_count):
c = cycle_cells[i]
rhs[i] = (
-target[c[0]] + target[c[1]] + target[c[2]]
- target[c[3]] - target[c[4]] + target[c[5]]
)
potential = np.zeros(cycle_count, dtype=np.float64)
residual = rhs.copy()
direction = residual.copy()
applied = np.empty(cycle_count, dtype=np.float64)
cell_scratch = np.empty(cell_count, dtype=np.float64)
rr = 0.0
for i in range(cycle_count):
rr += residual[i] * residual[i]
rr_initial = rr
if rr_initial == 0.0:
return potential, 0, 0.0
threshold = tolerance * tolerance * rr_initial
iteration = 0
while iteration < max_iterations and rr > threshold:
for j in range(cell_count):
cell_scratch[j] = 0.0
# cell_scratch = C.T @ direction
for i in range(cycle_count):
value = direction[i]
c = cycle_cells[i]
cell_scratch[c[0]] -= value
cell_scratch[c[1]] += value
cell_scratch[c[2]] += value
cell_scratch[c[3]] -= value
cell_scratch[c[4]] -= value
cell_scratch[c[5]] += value
# applied = C @ cell_scratch
denominator = 0.0
for i in range(cycle_count):
c = cycle_cells[i]
value = (
-cell_scratch[c[0]] + cell_scratch[c[1]] + cell_scratch[c[2]]
- cell_scratch[c[3]] - cell_scratch[c[4]] + cell_scratch[c[5]]
)
applied[i] = value
denominator += direction[i] * value
if denominator <= 1.0e-30:
break
alpha = rr / denominator
rr_new = 0.0
for i in range(cycle_count):
potential[i] += alpha * direction[i]
residual[i] -= alpha * applied[i]
rr_new += residual[i] * residual[i]
beta = rr_new / rr
for i in range(cycle_count):
direction[i] = residual[i] + beta * direction[i]
rr = rr_new
iteration += 1
return potential, iteration, math.sqrt(rr / rr_initial)
@njit(cache=True)
def _potential_interpolate_idw(
source_coordinates: np.ndarray,
source_values: np.ndarray,
destination_coordinates: np.ndarray,
neighbors: int,
power: float,
) -> np.ndarray:
"""K-nearest IDW on normalized cube coordinates, allocation once per call."""
source_count = source_coordinates.shape[0]
destination_count = destination_coordinates.shape[0]
neighbor_count = min(neighbors, source_count)
output = np.empty(destination_count, dtype=np.float64)
best_distance = np.empty(neighbor_count, dtype=np.float64)
best_index = np.empty(neighbor_count, dtype=np.int64)
exponent = 0.5 * power
for i in range(destination_count):
for k in range(neighbor_count):
best_distance[k] = 1.0e300
best_index[k] = -1
for j in range(source_count):
dq = destination_coordinates[i, 0] - source_coordinates[j, 0]
dr = destination_coordinates[i, 1] - source_coordinates[j, 1]
ds = destination_coordinates[i, 2] - source_coordinates[j, 2]
distance = dq * dq + dr * dr + ds * ds
if distance < best_distance[neighbor_count - 1]:
position = neighbor_count - 1
while position > 0 and distance < best_distance[position - 1]:
best_distance[position] = best_distance[position - 1]
best_index[position] = best_index[position - 1]
position -= 1
best_distance[position] = distance
best_index[position] = j
if best_distance[0] < 1.0e-28:
output[i] = source_values[best_index[0]]
else:
total_weight = 0.0
total_value = 0.0
for k in range(neighbor_count):
weight = 1.0 / (
best_distance[k] ** exponent + 1.0e-30
)
total_weight += weight
total_value += weight * source_values[best_index[k]]
output[i] = total_value / total_weight
return output
@njit(cache=True)
def _potential_cycle_transpose(
cycle_cells: np.ndarray,
potential: np.ndarray,
cell_count: int,
) -> np.ndarray:
"""Compute the value-field score C.T @ potential with the six-point stencil."""
output = np.zeros(cell_count, dtype=np.float64)
for i in range(cycle_cells.shape[0]):
value = potential[i]
c = cycle_cells[i]
output[c[0]] -= value
output[c[1]] += value
output[c[2]] += value
output[c[3]] -= value
output[c[4]] -= value
output[c[5]] += value
return output
def _potential_signed_rank(
reps: list[tuple[int, int, int]],
scores: np.ndarray,
) -> np.ndarray:
"""Turn arbitrary real scores into a signed permutation of 1..m."""
score = np.asarray(scores, dtype=np.float64)
signs = np.sign(score).astype(np.int64)
for raw_index in np.flatnonzero(signs == 0):
i = int(raw_index)
signs[i] = 1 if psi0(reps[i]) >= 0 else -1
coordinates = np.asarray(reps, dtype=np.int64)
# np.lexsort uses its final key as primary: |score|, then q, r, s.
rank_order = np.lexsort(
(
coordinates[:, 2],
coordinates[:, 1],
coordinates[:, 0],
np.abs(score),
)
)
magnitudes = np.empty(len(reps), dtype=np.int64)
magnitudes[rank_order] = np.arange(1, len(reps) + 1, dtype=np.int64)
return np.ascontiguousarray(signs * magnitudes)
def _potential_correlation(a: np.ndarray, b: np.ndarray) -> float:
aa = a - float(np.mean(a))
bb = b - float(np.mean(b))
denominator = float(np.linalg.norm(aa) * np.linalg.norm(bb))
if denominator == 0.0:
return 0.0
return float(np.dot(aa, bb) / denominator)
def _potential_parse_weights(raw: str) -> tuple[float, ...]:
weights: list[float] = []
for token in raw.split(","):
token = token.strip()
if not token:
continue
value = float(token)
if not 0.0 <= value <= 1.0:
raise ValueError("potential blend weights must lie in [0,1]")
if 0.0 < value < 1.0 and value not in weights:
weights.append(value)
if not weights:
raise ValueError("--potential-blend-weights contains no interior weights")
return tuple(weights)
def _potential_load_or_recover(
path: Path,
cache_dir: Path | None,
use_cache: bool,
tolerance: float,
max_iterations: int,
):
"""Load a cached potential or recover it from an exact MHX witness."""
import hashlib
hint_order, hint_map = parse_mhx(path)
source_bytes = path.read_bytes()
digest = hashlib.sha256(source_bytes).hexdigest()[:20]
cache_path = None
if use_cache and cache_dir is not None:
cache_dir.mkdir(parents=True, exist_ok=True)
cache_path = cache_dir / f"{path.stem}-n{hint_order}-{digest}.npz"
if cache_path.exists():
with np.load(cache_path, allow_pickle=False) as data:
centers = np.ascontiguousarray(data["centers"], dtype=np.int16)
potential = np.ascontiguousarray(data["potential"], dtype=np.float64)
iterations = int(data["iterations"])
relative_residual = float(data["relative_residual"])
return (
hint_order,
hint_map,
centers,
potential,
iterations,
relative_residual,
True,
)
hint_cells, centers, cycle_cells = _potential_cycle_data(hint_order)
scale = float(hint_order * (hint_order - 1))
target = np.asarray(
[scale * float(hint_map[c]) for c in hint_cells],
dtype=np.float64,
)
potential, iterations, relative_residual = _potential_recover_cg(
cycle_cells,
target,
tolerance,
max_iterations,
)
potential = np.ascontiguousarray(potential)
if relative_residual > max(1.0e-7, 100.0 * tolerance):
raise RuntimeError(
f"potential recovery did not converge for {path}: "
f"relative normal residual={relative_residual:.3e} "
f"after {iterations} iterations"
)
if cache_path is not None:
np.savez(
cache_path,
centers=centers,
potential=potential,
iterations=np.asarray(iterations, dtype=np.int64),
relative_residual=np.asarray(relative_residual, dtype=np.float64),
)
return (
hint_order,
hint_map,
centers,
potential,
iterations,
relative_residual,
False,
)
def build_potential_seed_bank(
*,
order: int,
reps: list[tuple[int, int, int]],
B: np.ndarray,
source_paths: list[Path],
bank_size: int,
idw_neighbors: int,
idw_power: float,
blend_weights: tuple[float, ...],
hybrid_weights: tuple[float, ...],
pair_source_limit: int,
alignment_limit: int,
use_pair_blends: bool,
minimum_distance: float,
cache_dir: Path | None,
use_cache: bool,
cg_tolerance: float,
cg_iterations: int,
):
"""Build ordinary signed-permutation seeds from smooth cross-order terrains.
Candidate generation is intentionally generous because it occurs once. The
target magic residual itself selects the best candidates; no heuristic
terrain score substitutes for the actual objective.
"""
if bank_size < 1:
raise ValueError("--potential-bank-size must be >= 1")
if idw_neighbors < 1:
raise ValueError("--potential-idw-neighbors must be >= 1")
if idw_power <= 0.0:
raise ValueError("--potential-idw-power must be > 0")
if pair_source_limit < 1:
raise ValueError("--potential-pair-sources must be >= 1")
if alignment_limit < 1:
raise ValueError("--potential-alignments must be >= 1")
if not 0.0 <= minimum_distance <= 1.0:
raise ValueError("--potential-min-seed-distance must lie in [0,1]")
# De-duplicate paths while preserving user order.
unique_paths: list[Path] = []
seen_paths: set[Path] = set()
for path in source_paths:
resolved = path.resolve()
if resolved not in seen_paths:
seen_paths.add(resolved)
unique_paths.append(path)
if not unique_paths:
raise ValueError("potential guidance requires at least one source")
target_cells, target_centers, target_cycles = _potential_cycle_data(order)
target_center_coordinates = (
target_centers.astype(np.float64) / float(max(1, order - 1))
)
target_cell_index = {c: i for i, c in enumerate(target_cells)}
representative_indices = np.asarray(
[target_cell_index[c] for c in reps],
dtype=np.int64,
)
B64 = B.astype(np.int64)
# Each source stores twelve target-grid terrains, one per D6 orientation.
sources: list[dict] = []
source_maps: list[tuple[int, dict[tuple[int, int, int], int]]] = []
for path in unique_paths:
(
hint_order,
hint_map,
source_centers,
potential,
iterations,
relative_residual,
cache_hit,
) = _potential_load_or_recover(
path,
cache_dir,
use_cache,
cg_tolerance,
cg_iterations,
)
source_maps.append((hint_order, hint_map))
rms = math.sqrt(float(np.mean(potential * potential)))
if rms == 0.0:
raise RuntimeError(f"zero recovered potential for {path}")
normalized_potential = potential / rms
transformed_terrains: list[np.ndarray] = []
for permutation, coordinate_sign in _POTENTIAL_D6_TRANSFORMS:
transformed_coordinates = np.ascontiguousarray(
coordinate_sign
* source_centers[:, permutation].astype(np.float64)
/ float(max(1, hint_order - 1))
)
terrain = _potential_interpolate_idw(
transformed_coordinates,
normalized_potential,
target_center_coordinates,
idw_neighbors,
idw_power,
)
target_rms = math.sqrt(float(np.mean(terrain * terrain)))
if target_rms > 0.0:
terrain /= target_rms
transformed_terrains.append(np.ascontiguousarray(terrain))
print(
f"potential-source path={path} order={hint_order} "
f"cycles={len(source_centers)} cg_iterations={iterations} "
f"relative_residual={relative_residual:.3e} "
f"cache={'hit' if cache_hit else 'miss'}",
flush=True,
)
sources.append(
{
"path": path,
"order": hint_order,
"terrains": transformed_terrains,
}
)
candidates: list[tuple[int, int, str, np.ndarray]] = []
seen_assignments: set[bytes] = set()
psi_seed = np.ascontiguousarray(exact_rank_seed(order, reps), dtype=np.int64)
psi_guide_score = psi_seed.astype(np.float64)
psi_guide_rms = math.sqrt(float(np.mean(psi_guide_score * psi_guide_score)))
psi_guide_score /= psi_guide_rms
def add_seed(seed: np.ndarray, label: str) -> None:
x = np.ascontiguousarray(seed, dtype=np.int64)
key = x.tobytes()
if key in seen_assignments:
return
seen_assignments.add(key)
residual = B64 @ x
energy = int(residual @ residual)
maximum_residual = int(np.max(np.abs(residual)))
candidates.append((energy, maximum_residual, label, x))
def add_terrain(terrain: np.ndarray, label: str) -> None:
full_scores = _potential_cycle_transpose(
target_cycles,
np.ascontiguousarray(terrain, dtype=np.float64),
len(target_cells),
)
representative_scores = np.ascontiguousarray(
full_scores[representative_indices],
dtype=np.float64,
)
add_seed(
_potential_signed_rank(reps, representative_scores),
label,
)
# Homotopy seeds inject a controlled amount of the learned terrain into
# the very low-residual psi0-rank score. This yields many structurally
# distinct basins without surrendering the excellent initial balance of
# psi0. It remains startup-only; the annealer sees an ordinary seed.
terrain_rms = math.sqrt(
float(np.mean(representative_scores * representative_scores))
)
if terrain_rms > 0.0:
normalized_terrain_score = representative_scores / terrain_rms
for terrain_weight in hybrid_weights:
hybrid_score = (
(1.0 - terrain_weight) * psi_guide_score
+ terrain_weight * normalized_terrain_score
)
add_seed(
_potential_signed_rank(reps, hybrid_score),
f"hybrid:psi0@{1.0-terrain_weight:.6g}+"
f"{label}@{terrain_weight:.6g}",
)
# Always retain legacy baselines in the competition. Potential guidance
# cannot make startup worse merely because a source landscape is atypical.
add_seed(psi_seed, "legacy:psi0-rank")
for source, (hint_order, hint_map) in zip(sources, source_maps):
add_seed(
cross_order_seed(order, reps, hint_order, hint_map),
f"legacy:nearest-cell:{source['path'].name}",
)
# Individual potential extrapolations provide robust fallbacks and a full
# D6 trajectory portfolio even when only one neighboring order is known.
for source in sources:
for transform_index, terrain in enumerate(source["terrains"]):
add_terrain(
terrain,
f"potential:{source['path'].name}:d6={transform_index}",
)
# Use the closest few source orders for pairwise interpolation. For every
# pair, retain the strongest D6 alignments and sweep asymmetric weights;
# this is especially useful when the higher-order terrain predicts the
# missing order better than the lower-order terrain.
closest_indices = sorted(
range(len(sources)),
key=lambda i: (abs(int(sources[i]["order"]) - order), i),
)[:pair_source_limit]
if use_pair_blends:
for left_position in range(len(closest_indices)):
for right_position in range(left_position + 1, len(closest_indices)):
left = sources[closest_indices[left_position]]
right = sources[closest_indices[right_position]]
alignments: list[tuple[float, float, int, int]] = []
for left_transform, left_terrain in enumerate(left["terrains"]):
for right_transform, right_terrain in enumerate(right["terrains"]):
correlation = _potential_correlation(
left_terrain,
right_terrain,
)
alignments.append(
(
abs(correlation),
correlation,
left_transform,
right_transform,
)
)
alignments.sort(reverse=True)
for (
_absolute_correlation,
correlation,
left_transform,
right_transform,
) in alignments[:alignment_limit]:
left_terrain = left["terrains"][left_transform]
right_terrain = right["terrains"][right_transform]
if correlation < 0.0:
right_terrain = -right_terrain
for weight in blend_weights:
terrain = (
weight * left_terrain
+ (1.0 - weight) * right_terrain
)
add_terrain(
terrain,
"blend:"
f"{left['path'].name}@{weight:.6g}:d6={left_transform}+"
f"{right['path'].name}@{1.0-weight:.6g}:d6={right_transform}",
)
# Multi-source consensus: use each orientation of the closest source as an
# anchor, align every other source to it, then average with inverse order
# distance. This costs only twelve additional candidates.
if len(sources) >= 2:
anchor_index = closest_indices[0]
anchor_source = sources[anchor_index]
for anchor_transform, anchor_terrain in enumerate(anchor_source["terrains"]):
weighted_sum = np.zeros_like(anchor_terrain)
total_weight = 0.0
anchor_weight = 1.0 / (
1.0 + abs(int(anchor_source["order"]) - order)
)
weighted_sum += anchor_weight * anchor_terrain
total_weight += anchor_weight
labels = [
f"{anchor_source['path'].name}:d6={anchor_transform}"
]
for source_index, source in enumerate(sources):
if source_index == anchor_index:
continue
best_correlation = -1.0
best_signed_terrain = source["terrains"][0]
best_transform = 0
for transform_index, terrain in enumerate(source["terrains"]):
correlation = _potential_correlation(anchor_terrain, terrain)
absolute_correlation = abs(correlation)
if absolute_correlation > best_correlation:
best_correlation = absolute_correlation
best_signed_terrain = terrain if correlation >= 0.0 else -terrain
best_transform = transform_index
source_weight = 1.0 / (
1.0 + abs(int(source["order"]) - order)
)
weighted_sum += source_weight * best_signed_terrain
total_weight += source_weight
labels.append(f"{source['path'].name}:d6={best_transform}")
add_terrain(
weighted_sum / total_weight,
"consensus:" + "+".join(labels),
)
candidates.sort(key=lambda item: (item[0], item[1], item[2]))
if not candidates:
raise RuntimeError("potential guidance generated no seed candidates")
# Greedy diversity filtering prevents a seed bank consisting solely of tiny
# rank perturbations of one blend. Equality compares complete signed
# placements; the bank is small, so this startup-only O(bank^2*m) pass is
# negligible.
selected: list[tuple[int, int, str, np.ndarray]] = []
for candidate in candidates:
x = candidate[3]
sufficiently_different = True
for prior in selected:
distance = float(np.mean(x != prior[3]))
if distance < minimum_distance:
sufficiently_different = False
break
if sufficiently_different:
selected.append(candidate)
if len(selected) >= bank_size:
break
if len(selected) < bank_size:
selected_keys = {item[3].tobytes() for item in selected}
for candidate in candidates:
key = candidate[3].tobytes()
if key not in selected_keys:
selected.append(candidate)
selected_keys.add(key)
if len(selected) >= bank_size:
break
print(
f"potential-seed-bank generated={len(candidates)} "
f"selected={len(selected)} target_order={order}",
flush=True,
)
for rank, (energy, maximum_residual, label, _x) in enumerate(selected):
print(
f"potential-seed rank={rank} energy={energy} "
f"max_residual={maximum_residual} source={label}",
flush=True,
)
return (
[np.ascontiguousarray(item[3]) for item in selected],
selected,
)
# R3_POTENTIAL_GUIDANCE_V1_END
def sparse_columns(B: np.ndarray):
m = B.shape[1]
rows = np.full((m, 3), -1, dtype=np.int16)
coeff = np.zeros((m, 3), dtype=np.int8)
degree = np.zeros(m, dtype=np.int8)
for i in range(m):
nz = np.flatnonzero(B[:, i])
degree[i] = len(nz)
rows[i, :len(nz)] = nz
coeff[i, :len(nz)] = B[nz, i]
return rows, coeff, degree
# HOTLOOP_OPTIMIZED_V2_BEGIN
# R3_OBSERVABILITY_V3_APPLIED
@njit(cache=True, inline="always")
def _rotl64(x: np.uint64, k: int) -> np.uint64:
return (x << k) | (x >> (64 - k))
@njit(cache=True, inline="always")
def _splitmix64_next(state: np.ndarray) -> np.uint64:
"""Advance one-word SplitMix64 state; used only to seed xoroshiro."""
state[0] += np.uint64(0x9E3779B97F4A7C15)
z = state[0]
z = (z ^ (z >> 30)) * np.uint64(0xBF58476D1CE4E5B9)
z = (z ^ (z >> 27)) * np.uint64(0x94D049BB133111EB)
return z ^ (z >> 31)
@njit(cache=True, inline="always")
def _rng_next_u64(state: np.ndarray) -> np.uint64:
"""xoroshiro128+; state is private to one annealing worker."""
s0 = state[0]
s1 = state[1]
result = s0 + s1
s1 ^= s0
state[0] = _rotl64(s0, 55) ^ s1 ^ (s1 << 14)
state[1] = _rotl64(s1, 36)
return result
@njit(cache=True, inline="always")
def _rng_bounded(state: np.ndarray, bound: int) -> int:
# Multiply-high reduction avoids the hardware integer division generated by
# ``random % bound``. For bound << 2^32 the bias is tiny and harmless for
# this stochastic search, just as with the previous modulo reduction.
random32 = _rng_next_u64(state) >> np.uint64(32)
return np.int64((random32 * np.uint64(bound)) >> np.uint64(32))
@njit(cache=True, inline="always")
def _rng_float(state: np.ndarray) -> float:
# Uniform double in [0,1), using the high 53 random bits.
return float(_rng_next_u64(state) >> 11) * 1.1102230246251565e-16
@njit(cache=True)
def _anneal(
x0: np.ndarray,
rows: np.ndarray,
coeff: np.ndarray,
degree: np.ndarray,
steps: int,
seed: int,
t0: float,
t1: float,
epoch: int,
triple_probability: float,
quad_probability: float,
progress_steps: int = 0,
worker_id: int = -1,
acceptance_cache_size: int = 0,
):
"""Allocation-free proposal loop with optional low-rate progress."""
x = x0.copy()
m = len(x)
line_count = 0
for i in range(m):
for k in range(degree[i]):
candidate = int(rows[i, k]) + 1
if candidate > line_count:
line_count = candidate
residual = np.zeros(line_count, dtype=np.int64)
for i in range(m):
value = x[i]
for k in range(degree[i]):
residual[rows[i, k]] += coeff[i, k] * value
energy = np.int64(0)
for rr in range(line_count):
energy += residual[rr] * residual[rr]
best_energy = energy
best_x = np.empty_like(x)
best_residual = np.empty_like(residual)
for i in range(m):
best_x[i] = x[i]
for rr in range(line_count):
best_residual[rr] = residual[rr]
# These buffers are allocated once per worker invocation and reused for
# every proposal. The old loop allocated 5-8 NumPy arrays per proposal.
move_idx = np.empty(4, dtype=np.int32)
magnitudes = np.empty(4, dtype=np.int64)
new_values = np.empty(4, dtype=np.int64)
permutation = np.empty(4, dtype=np.int8)
# A direct row-indexed accumulator replaces the old linear search through
# an allocated list of affected rows. At most 12 rows are touched.
row_delta = np.zeros(line_count, dtype=np.int64)
touched_rows = np.empty(12, dtype=np.int16)
row_touched = np.zeros(line_count, dtype=np.uint8)
# Seed a private xoroshiro128+ stream. No TLS lookup and no
# numba_get_np_random_state()/numba_rnd_shuffle() calls remain in the loop.
splitmix_state = np.empty(1, dtype=np.uint64)
splitmix_state[0] = np.uint64(seed) ^ np.uint64(0xD1B54A32D192ED03)
rng_state = np.empty(2, dtype=np.uint64)
rng_state[0] = _splitmix64_next(splitmix_state)
rng_state[1] = _splitmix64_next(splitmix_state)
if rng_state[0] == 0 and rng_state[1] == 0:
rng_state[1] = np.uint64(1)
# Optional RAM-for-speed trade. If enabled, precompute Exp(1) variates.
# U < exp(-dE/T) is exactly equivalent to dE < T*(-log(U)). Reusing a
# large, worker-private table removes both libm exp() and one RNG draw from
# every uphill Metropolis test. A power-of-two size permits mask indexing.
cache_slots = acceptance_cache_size if acceptance_cache_size > 0 else 1
acceptance_thresholds = np.empty(cache_slots, dtype=np.float64)
acceptance_mask = cache_slots - 1
acceptance_position = 0
acceptance_stride = 1
if acceptance_cache_size > 0:
for cache_i in range(acceptance_cache_size):
cache_u = _rng_float(rng_state)
if cache_u <= 0.0:
cache_u = 1.1102230246251565e-16
acceptance_thresholds[cache_i] = -math.log(cache_u)
acceptance_position = _rng_bounded(rng_state, acceptance_cache_size)
acceptance_stride = _rng_bounded(rng_state, acceptance_cache_size) | 1
accepted_moves = np.int64(0)
next_progress = progress_steps
if best_energy == 0:
return best_energy, best_x, best_residual, 0
if epoch < 1:
epoch = 1
if epoch <= 1 or t0 <= 0.0 or t1 <= 0.0:
cooling = 1.0
else:
cooling = math.exp(math.log(t1 / t0) / float(epoch - 1))
temperature = t0
epoch_position = 0
for step in range(steps):
move_draw = _rng_float(rng_state)
if m >= 4 and move_draw < quad_probability:
move_size = 4
elif m >= 3 and move_draw < quad_probability + triple_probability:
move_size = 3
else:
move_size = 2
# Draw distinct target positions. Rejection is very cheap because
# move_size <= 4 while m=234 already at order 13.
for u in range(move_size):
while True:
candidate = _rng_bounded(rng_state, m)
clash = False
for w in range(u):
if move_idx[w] == candidate:
clash = True
break
if not clash:
move_idx[u] = candidate
break
old_value = x[move_idx[u]]
magnitudes[u] = old_value if old_value >= 0 else -old_value
permutation[u] = u
# Manual Fisher-Yates using the local generator. No temporary array
# and no Numba np.random.shuffle helper.
for u in range(move_size - 1, 0, -1):
v = _rng_bounded(rng_state, u + 1)
tmp = permutation[u]
permutation[u] = permutation[v]
permutation[v] = tmp
# One 64-bit draw supplies all independent sign bits for the move.
sign_bits = _rng_next_u64(rng_state)
changed = False
for u in range(move_size):
value = magnitudes[permutation[u]]
if ((sign_bits >> u) & np.uint64(1)) == 0:
value = -value
new_values[u] = value
if value != x[move_idx[u]]:
changed = True
if changed:
touched_count = 0
for u in range(move_size):
col = move_idx[u]
value_delta = new_values[u] - x[col]
if value_delta == 0:
continue
for z in range(degree[col]):
rr = rows[col, z]
if row_touched[rr] == 0:
row_touched[rr] = 1
touched_rows[touched_count] = rr
touched_count += 1
row_delta[rr] += coeff[col, z] * value_delta
delta_energy = np.int64(0)
for w in range(touched_count):
rr = touched_rows[w]
delta = row_delta[rr]
delta_energy += 2 * residual[rr] * delta + delta * delta
accept = delta_energy <= 0
if not accept and temperature > 0.0:
if acceptance_cache_size > 0:
threshold = acceptance_thresholds[acceptance_position]
acceptance_position = (
acceptance_position + acceptance_stride
) & acceptance_mask
accept = float(delta_energy) < temperature * threshold
else:
accept = _rng_float(rng_state) < math.exp(
-float(delta_energy) / temperature
)
if accept:
accepted_moves += 1
for u in range(move_size):
x[move_idx[u]] = new_values[u]
for w in range(touched_count):
rr = touched_rows[w]
residual[rr] += row_delta[rr]
energy += delta_energy
if energy < best_energy:
best_energy = energy
# Copy into already-owned output buffers. The old version
# allocated two new arrays on every incumbent improvement.
for i in range(m):
best_x[i] = x[i]
for rr in range(line_count):
best_residual[rr] = residual[rr]
# Clear only the rows touched by this proposal.
for w in range(touched_count):
rr = touched_rows[w]
row_delta[rr] = 0
row_touched[rr] = 0
if best_energy == 0:
return best_energy, best_x, best_residual, step + 1
# Multiplicative schedule: one exp/log pair per anneal call instead of
# a power operation on every proposal.
epoch_position += 1
if epoch_position >= epoch:
epoch_position = 0
temperature = t0
else:
temperature *= cooling
# One predictable comparison per proposal when progress is enabled.
# Printing occurs only at the requested large interval.
if progress_steps > 0 and step + 1 >= next_progress:
print(
"r3-progress worker", worker_id,
"step", step + 1,
"of", steps,
"best", best_energy,
"current", energy,
"accepted", accepted_moves,
)
next_progress += progress_steps
return best_energy, best_x, best_residual, steps
# HOTLOOP_OPTIMIZED_V2_END
# R3_NUMBA_PERF_MAP_BEGIN
_NUMBA_PERF_MAP_PIDS: set[int] = set()
def _write_numba_perf_map(dispatchers) -> None:
"""Write a Linux perf map for the current process's Numba code.
Numba/LLVM profiling events are useful for ``perf record`` +
``perf inject --jit``. ``perf top`` also understands the conventional
/tmp/perf-<pid>.map format, so this helper publishes live symbol ranges.
This is debug-only code. It deliberately uses Numba's private compiled
object accessor and the system ``nm`` command, and therefore fails softly.
"""
import subprocess
import tempfile
pid = os.getpid()
if pid in _NUMBA_PERF_MAP_PIDS:
return
output_lines: list[str] = []
try:
for label, dispatcher in dispatchers:
for signature, cres in dispatcher.overloads.items():
obj = cres.library._get_compiled_object()
if not obj:
continue
tmp_name = None
try:
with tempfile.NamedTemporaryFile(prefix="r3-numba-", suffix=".o", delete=False) as tmp:
tmp.write(obj)
tmp_name = tmp.name
completed = subprocess.run(
["nm", "-S", "--defined-only", tmp_name],
check=True,
text=True,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
)
finally:
if tmp_name is not None:
try:
os.unlink(tmp_name)
except OSError:
pass
native_name = cres.fndesc.llvm_func_name
cpython_name = cres.fndesc.llvm_cpython_wrapper_name
cfunc_name = cres.fndesc.llvm_cfunc_wrapper_name
signature_text = str(signature).replace(" ", "")
for line in completed.stdout.splitlines():
fields = line.split(maxsplit=3)
if len(fields) != 4:
continue
_offset_hex, size_hex, symbol_type, symbol_name = fields
if symbol_type not in ("t", "T"):
continue
try:
size = int(size_hex, 16)
address = int(cres.library.get_pointer_to_function(symbol_name))
except Exception:
continue
if address == 0 or size <= 0:
continue
if symbol_name == native_name:
friendly = f"numba::{label}[{signature_text}]"
elif symbol_name == cpython_name:
friendly = f"numba::cpython::{label}[{signature_text}]"
elif symbol_name == cfunc_name:
friendly = f"numba::cfunc::{label}[{signature_text}]"
else:
friendly = f"numba::{label}::{symbol_name}"
output_lines.append(f"{address:x} {size:x} {friendly}")
if output_lines:
map_path = Path(f"/tmp/perf-{pid}.map")
map_path.write_text("\n".join(output_lines) + "\n", encoding="utf-8")
_NUMBA_PERF_MAP_PIDS.add(pid)
except Exception as exc:
print(f"warning: could not create Numba perf map for pid={pid}: {exc}", flush=True)
# R3_NUMBA_PERF_MAP_END
def _anneal_job(payload):
"""Pickle-friendly wrapper used by ProcessPoolExecutor."""
(
x0,
rows,
coeff,
degree,
steps,
seed,
t0,
t1,
epoch,
triple_probability,
quad_probability,
progress_steps,
worker_id,
acceptance_cache_size,
numba_perf,
) = payload
if progress_steps > 0:
try:
sys.stdout.reconfigure(line_buffering=True, write_through=True)
except Exception:
pass
# Under Linux/fork the dispatcher is already warm. Under spawn, compile a
# one-step specialization before publishing the perf map.
if numba_perf and not _anneal.signatures:
_anneal(
x0,
rows,
coeff,
degree,
1,
seed,
1.0,
1.0,
1,
triple_probability,
quad_probability,
0,
worker_id,
0,
)
if numba_perf:
_write_numba_perf_map((("r3_anneal", _anneal),))
return _anneal(
x0,
rows,
coeff,
degree,
steps,
seed,
t0,
t1,
epoch,
triple_probability,
quad_probability,
progress_steps,
worker_id,
acceptance_cache_size,
)
def _stop_executor_now(
pool: concurrent.futures.ProcessPoolExecutor,
futures,
) -> None:
"""Cancel queued jobs and terminate running worker processes immediately.
``Future.cancel()`` only affects jobs that have not started. Python 3.13's
ProcessPoolExecutor also has no public API for killing running workers, so
on that version we use its process table as a compatibility fallback.
Python 3.14+ is handled through ``terminate_workers()`` when available.
"""
for future in futures:
future.cancel()
terminate_workers = getattr(pool, "terminate_workers", None)
if terminate_workers is not None:
terminate_workers()
return
# Python 3.13 compatibility. Capture the processes before shutdown clears
# the executor's private process table.
processes = list(getattr(pool, "_processes", {}).values())
for process in processes:
if process.is_alive():
process.terminate()
# Give SIGTERM/TerminateProcess a brief chance, then escalate. The short
# joins reap the children without waiting for their configured anneal steps.
deadline = time.monotonic() + 1.0
for process in processes:
process.join(timeout=max(0.0, deadline - time.monotonic()))
for process in processes:
if process.is_alive():
kill = getattr(process, "kill", None)
if kill is not None:
kill()
else:
process.terminate()
for process in processes:
process.join(timeout=0.2)
pool.shutdown(wait=False, cancel_futures=True)
def perturb_state(x: np.ndarray, seed: int, moves: int) -> np.ndarray:
"""Diversify an incumbent while preserving the signed absolute permutation."""
if moves <= 0:
return x.copy()
rng = np.random.default_rng(seed)
y = x.copy()
m = len(y)
for _ in range(moves):
k = 3 if rng.random() < 0.7 else 4
idx = rng.choice(m, size=k, replace=False)
mags = np.abs(y[idx]).copy()
rng.shuffle(mags)
signs = rng.choice(np.array([-1, 1], dtype=np.int64), size=k)
y[idx] = signs * mags
return y
def _select_active_columns(
B: np.ndarray,
x: np.ndarray,
residual: np.ndarray,
size: int,
closure_layers: int,
seed: int,
) -> np.ndarray:
"""Choose a sparse exact-LNS neighborhood around the defective rows.
The selection deliberately mixes deterministic min-conflict columns with
randomized closure columns. Repeated finisher attempts therefore explore
genuinely different exact neighborhoods.
"""
rng = np.random.default_rng(seed)
m = B.shape[1]
size = max(4, min(int(size), m))
defect_rows = np.flatnonzero(residual)
if len(defect_rows) == 0:
return np.arange(min(size, m), dtype=np.int64)
# Correlation with the current residual is the first-order energy gradient.
gradient = np.abs(B.astype(np.int64).T @ residual.astype(np.int64)).astype(float)
incidence = np.count_nonzero(B[defect_rows, :], axis=0).astype(float)
jitter = rng.random(m) * 1e-3
score = gradient + 4.0 * incidence + jitter
active: set[int] = set()
# Ensure every defective row has several adjustable columns.
quota = max(4, min(12, size // max(1, len(defect_rows))))
for rr in defect_rows[np.argsort(-np.abs(residual[defect_rows]))]:
cols = np.flatnonzero(B[rr, :])
if len(cols) == 0:
continue
ranked = cols[np.argsort(-score[cols])]
for col in ranked[:quota]:
active.add(int(col))
if len(active) >= size:
break
if len(active) >= size:
break
# Add the globally strongest columns incident to a defect.
defect_cols = np.flatnonzero(incidence > 0)
for col in defect_cols[np.argsort(-score[defect_cols])]:
active.add(int(col))
if len(active) >= size:
break
# One or more incidence-closure layers let the exact solver compensate on
# rows that would otherwise become newly defective.
for _ in range(max(0, int(closure_layers))):
if len(active) >= size or not active:
break
active_arr = np.fromiter(active, dtype=np.int64)
touched_rows = np.flatnonzero(np.any(B[:, active_arr] != 0, axis=1))
closure_cols = np.flatnonzero(np.any(B[touched_rows, :] != 0, axis=0))
# Randomized tie breaking is important when the incumbent residual is tiny.
closure_score = score[closure_cols] + rng.random(len(closure_cols))
for col in closure_cols[np.argsort(-closure_score)]:
active.add(int(col))
if len(active) >= size:
break
# Fill any remaining slots with a mixture of strong and random columns.
if len(active) < size:
remaining = np.array([i for i in range(m) if i not in active], dtype=np.int64)
if len(remaining):
mixed = score[remaining] + 0.5 * rng.random(len(remaining))
for col in remaining[np.argsort(-mixed)]:
active.add(int(col))
if len(active) >= size:
break
return np.array(sorted(active), dtype=np.int64)
def exact_lns_finisher(
B: np.ndarray,
x: np.ndarray,
residual: np.ndarray,
*,
attempts: int,
active_size: int,
closure_layers: int,
time_limit: float,
workers: int,
seed: int,
) -> tuple[bool, np.ndarray, np.ndarray, dict]:
"""Try exact signed-permutation repair on small CP-SAT neighborhoods.
Outside the active set all labels are fixed. Inside it, CP-SAT may apply
any signed permutation of the magnitudes currently present there. Thus the
global consecutive shell is preserved exactly, while every line equation is
imposed as an equality.
"""
from ortools.sat.python import cp_model
m = len(x)
best_x = x.copy()
best_residual = residual.copy()
best_energy = int(residual @ residual)
stats = {"attempts": 0, "best_energy": best_energy, "active_sizes": []}
for attempt in range(max(0, int(attempts))):
active = _select_active_columns(
B,
best_x,
best_residual,
size=active_size,
closure_layers=closure_layers,
seed=seed + 7919 * attempt,
)
stats["attempts"] += 1
stats["active_sizes"].append(int(len(active)))
active_set = set(map(int, active.tolist()))
fixed = np.array([i for i in range(m) if i not in active_set], dtype=np.int64)
fixed_sum = (
B[:, fixed].astype(np.int64) @ best_x[fixed]
if len(fixed)
else np.zeros(B.shape[0], dtype=np.int64)
)
target = -fixed_sum
mags = np.abs(best_x[active]).astype(np.int64)
mag_list = [int(v) for v in mags.tolist()]
mag_to_slot = {v: i for i, v in enumerate(mag_list)}
k = len(active)
model = cp_model.CpModel()
slot_vars = [model.NewIntVar(0, k - 1, f"slot_{j}") for j in range(k)]
model.AddAllDifferent(slot_vars)
mag_vars = [model.NewIntVar(1, m, f"mag_{j}") for j in range(k)]
sign_vars = [model.NewBoolVar(f"pos_{j}") for j in range(k)]
val_vars = [model.NewIntVar(-m, m, f"val_{j}") for j in range(k)]
for j in range(k):
model.AddElement(slot_vars[j], mag_list, mag_vars[j])
model.Add(val_vars[j] == mag_vars[j]).OnlyEnforceIf(sign_vars[j])
model.Add(val_vars[j] == -mag_vars[j]).OnlyEnforceIf(sign_vars[j].Not())
old = int(best_x[active[j]])
model.AddHint(slot_vars[j], mag_to_slot[abs(old)])
model.AddHint(mag_vars[j], abs(old))
model.AddHint(sign_vars[j], 1 if old > 0 else 0)
model.AddHint(val_vars[j], old)
# Every row is exact. Rows not touched by the active set reduce to a
# constant equality and cheaply reject a bad neighborhood.
for rr in range(B.shape[0]):
terms = [
int(B[rr, active[j]]) * val_vars[j]
for j in range(k)
if B[rr, active[j]] != 0
]
if terms:
model.Add(sum(terms) == int(target[rr]))
elif int(target[rr]) != 0:
# The neighborhood cannot repair this row.
model.Add(0 == 1)
model.AddDecisionStrategy(
slot_vars,
cp_model.CHOOSE_MIN_DOMAIN_SIZE,
cp_model.SELECT_MIN_VALUE,
)
solver = cp_model.CpSolver()
solver.parameters.max_time_in_seconds = float(time_limit)
solver.parameters.num_search_workers = max(1, int(workers))
solver.parameters.random_seed = int(seed + 104729 * attempt)
solver.parameters.log_search_progress = False
solver.parameters.cp_model_presolve = True
status = solver.Solve(model)
if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
continue
candidate = best_x.copy()
for j, col in enumerate(active):
candidate[col] = solver.Value(val_vars[j])
candidate_residual = B.astype(np.int64) @ candidate
candidate_energy = int(candidate_residual @ candidate_residual)
if candidate_energy < best_energy:
best_energy = candidate_energy
best_x = candidate
best_residual = candidate_residual
stats["best_energy"] = best_energy
if candidate_energy == 0:
return True, candidate, candidate_residual, stats
return False, best_x, best_residual, stats
def write_mhx(path: Path, order: int, reps, x: np.ndarray, exact: bool, residual: np.ndarray):
cs = cells(order)
cmap: dict[tuple[int, int, int], int] = {(0, 0, 0): 0}
x_by_rep = {c: int(x[i]) for i, c in enumerate(reps)}
for c in cs:
if c == (0, 0, 0):
continue
nc = (-c[0], -c[1], -c[2])
rep = c if c > nc else nc
cmap[c] = x_by_rep[rep] if c == rep else -x_by_rep[rep]
R = order - 1
rows_out = []
for r in range(-R, R+1):
q0 = max(-R, -r-R)
q1 = min(R, -r+R)
rows_out.append([cmap[(q, r, -q-r)] for q in range(q0, q1+1)])
m = len(reps)
source = "r3 signed-permutation exchange solver"
if not exact:
source += f"; NONMAGIC BEST CANDIDATE energy={int(residual@residual)} max_residual={int(np.max(np.abs(residual)))}"
lines = [
"MHX 1", "kind: filled", f"order: {order}", "magic: 0",
f"domain: {-m}..{m}", f"source: {source}",
"layout: rows are top-to-bottom; values in each row are left-to-right", "", "rows:",
]
lines.extend(" ".join(map(str, row)) for row in rows_out)
lines.append("")
path.parent.mkdir(parents=True, exist_ok=True)
path.write_text("\n".join(lines), encoding="utf-8")
def main() -> int:
ap = argparse.ArgumentParser()
ap.add_argument("--order", type=int, required=True)
ap.add_argument("--hint", type=Path)
ap.add_argument("--allow-cross-order", action="store_true")
ap.add_argument(
"--potential-guide",
action="store_true",
help=(
"Recover the local-cycle potential of --hint and seed from its "
"interpolated terrain instead of nearest cell values."
),
)
ap.add_argument(
"--potential-hint",
action="append",
type=Path,
default=[],
help=(
"Additional exact MHX source for potential-field guidance. "
"Repeat for multiple orders; pairwise D6-aligned blends are tested."
),
)
ap.add_argument(
"--no-potential-include-hint",
action="store_true",
help="Do not add --hint to a nonempty --potential-hint source list.",
)
ap.add_argument(
"--potential-bank-size",
type=int,
default=32,
help="Number of diverse terrain-derived signed-permutation seeds to retain.",
)
ap.add_argument(
"--potential-reset-every",
type=int,
default=4,
help=(
"Start every Kth restart from the next potential seed; 0 disables "
"periodic resets after the initial seed-bank launch."
),
)
ap.add_argument(
"--potential-idw-neighbors",
type=int,
default=12,
help="Nearest source cycle centers used for normalized-coordinate IDW.",
)
ap.add_argument(
"--potential-idw-power",
type=float,
default=2.0,
help="Inverse-distance interpolation power for potential terrains.",
)
ap.add_argument(
"--potential-blend-weights",
default=(
"0.025,0.05,0.075,0.1,0.125,0.15,0.2,0.3,0.5,"
"0.7,0.8,0.85,0.9,0.925,0.95,0.975"
),
help="Comma-separated first-source weights used in pairwise blends.",
)
ap.add_argument(
"--potential-hybrid-weights",
default="0.005,0.01,0.02,0.03,0.05,0.075,0.1,0.15,0.2",
help=(
"Comma-separated terrain fractions mixed into the psi0-rank score "
"before reranking; creates low-energy terrain-aware homotopy seeds."
),
)
ap.add_argument(
"--no-potential-hybrids",
action="store_true",
help="Skip psi0/terrain hybrid seeds and retain only pure terrains.",
)
ap.add_argument(
"--potential-pair-sources",
type=int,
default=4,
help="Use at most this many closest-order sources for pairwise blends.",
)
ap.add_argument(
"--potential-alignments",
type=int,
default=48,
help="Best D6 alignment pairs retained per source pair.",
)
ap.add_argument(
"--no-potential-pair-blends",
action="store_true",
help="Generate individual/consensus potential seeds but skip pair blends.",
)
ap.add_argument(
"--potential-min-seed-distance",
type=float,
default=0.05,
help="Minimum fraction of differing signed placements in the selected bank.",
)
ap.add_argument(
"--potential-cache-dir",
type=Path,
default=Path("artifacts/r3_potential_cache"),
help="Directory for recovered-potential .npz caches.",
)
ap.add_argument(
"--no-potential-cache",
action="store_true",
help="Recompute source potentials instead of reading/writing cache files.",
)
ap.add_argument(
"--potential-cg-tolerance",
type=float,
default=1.0e-10,
help="Relative normal-residual tolerance for sparse potential recovery.",
)
ap.add_argument(
"--potential-cg-iterations",
type=int,
default=20000,
help="Maximum sparse-CG iterations per uncached source potential.",
)
ap.add_argument(
"--time-limit",
type=float,
default=300.0,
help="Finite wall-clock scheduling budget in seconds. Ignored with --until-found.",
)
ap.add_argument(
"--until-found",
action="store_true",
help="Run unbounded batches of parallel restarts until an exact solution is found.",
)
ap.add_argument(
"--workers",
type=int,
default=max(1, min(16, os.cpu_count() or 1)),
help="Parallel annealing processes (default: min(16, CPU count)).",
)
ap.add_argument(
"--restarts",
type=int,
default=32,
help="Total restart count in finite mode; ignored as a cap with --until-found.",
)
ap.add_argument("--steps", type=int, default=5_000_000)
ap.add_argument("--epoch", type=int, default=500_000)
ap.add_argument("--temperature-start", type=float, default=100.0)
ap.add_argument("--temperature-end", type=float, default=0.005)
ap.add_argument("--triple-probability", type=float, default=0.7)
ap.add_argument("--quad-probability", type=float, default=0.05)
ap.add_argument(
"--perturbations",
type=int,
default=24,
help="Signed 3/4-position perturbations applied before most restarts.",
)
ap.add_argument("--seed", type=int, default=1)
progress_group = ap.add_mutually_exclusive_group()
progress_group.add_argument(
"--progress-steps",
type=int,
default=0,
metavar="K",
help="Emit worker progress every K annealing proposals (0 disables).",
)
progress_group.add_argument(
"--progress-epoch",
action="store_true",
help="Emit worker progress at the end of every annealing epoch.",
)
ap.add_argument(
"--acceptance-cache-size",
type=int,
default=0,
metavar="N",
help=(
"Power-of-two count of cached Exp(1) thresholds per worker; "
"removes exp() from uphill acceptance tests. 0 disables."
),
)
ap.add_argument(
"--numba-perf",
action="store_true",
help=(
"Enable Numba profiling/debug metadata and write live "
"/tmp/perf-PID.map symbols for Linux perf."
),
)
ap.add_argument("--out", type=Path)
ap.add_argument("--best-out", type=Path)
ap.add_argument(
"--finisher",
action=argparse.BooleanOptionalAction,
default=True,
help="Enable the exact CP-SAT terminal LNS finisher.",
)
ap.add_argument(
"--finisher-energy",
type=int,
default=100,
help="Invoke the exact finisher when the incumbent energy is at most this value.",
)
ap.add_argument("--finisher-size", type=int, default=72)
ap.add_argument("--finisher-attempts", type=int, default=8)
ap.add_argument("--finisher-closure", type=int, default=1)
ap.add_argument(
"--finisher-time-limit",
type=float,
default=20.0,
help="CP-SAT time limit per exact-neighborhood attempt.",
)
ap.add_argument(
"--finisher-workers",
type=int,
default=0,
help="CP-SAT workers for the finisher; 0 means use --workers.",
)
ap.add_argument(
"--finisher-every",
type=int,
default=4,
help="Retry the finisher every N parallel batches even without an energy improvement.",
)
args = ap.parse_args()
if args.acceptance_cache_size < 0:
raise SystemExit("--acceptance-cache-size must be >= 0")
if (
args.acceptance_cache_size > 0
and args.acceptance_cache_size & (args.acceptance_cache_size - 1)
):
raise SystemExit("--acceptance-cache-size must be a power of two")
if args.progress_steps < 0:
raise SystemExit("--progress-steps must be >= 0")
if args.progress_steps > 0 or args.progress_epoch:
try:
sys.stdout.reconfigure(line_buffering=True, write_through=True)
except Exception:
pass
if args.order < 2:
raise SystemExit("--order must be >= 2")
if args.workers < 1:
raise SystemExit("--workers must be >= 1")
if not args.until_found and args.restarts < 1:
raise SystemExit("--restarts must be >= 1 in finite mode")
if not args.until_found and args.time_limit <= 0:
raise SystemExit("--time-limit must be positive unless --until-found is used")
_cs, reps, _rid, B = representative_data(args.order)
B64 = B.astype(np.int64)
rows, coeff, degree = sparse_columns(B)
potential_paths = list(args.potential_hint or [])
if args.hint is not None and (
args.potential_guide
or (potential_paths and not args.no_potential_include_hint)
):
hint_resolved = args.hint.resolve()
if all(path.resolve() != hint_resolved for path in potential_paths):
potential_paths.insert(0, args.hint)
if potential_paths:
try:
potential_weights = _potential_parse_weights(
args.potential_blend_weights
)
potential_hybrid_weights = (
()
if args.no_potential_hybrids
else _potential_parse_weights(args.potential_hybrid_weights)
)
seed_bank, potential_seed_records = build_potential_seed_bank(
order=args.order,
reps=reps,
B=B,
source_paths=potential_paths,
bank_size=args.potential_bank_size,
idw_neighbors=args.potential_idw_neighbors,
idw_power=args.potential_idw_power,
blend_weights=potential_weights,
hybrid_weights=potential_hybrid_weights,
pair_source_limit=args.potential_pair_sources,
alignment_limit=args.potential_alignments,
use_pair_blends=not args.no_potential_pair_blends,
minimum_distance=args.potential_min_seed_distance,
cache_dir=args.potential_cache_dir,
use_cache=not args.no_potential_cache,
cg_tolerance=args.potential_cg_tolerance,
cg_iterations=args.potential_cg_iterations,
)
except (ValueError, RuntimeError, KeyError, OSError) as exc:
raise SystemExit(f"potential guidance failed: {exc}") from exc
seed_x = seed_bank[0].copy()
seed_desc = (
f"potential bank ({len(seed_bank)} seeds from "
f"{len(potential_paths)} source(s))"
)
elif args.hint is None:
seed_x = exact_rank_seed(args.order, reps)
seed_desc = "psi0 rank"
seed_bank = [np.ascontiguousarray(seed_x)]
else:
hn, hmap = parse_mhx(args.hint)
if hn != args.order and not args.allow_cross_order:
raise SystemExit("hint order differs; add --allow-cross-order")
seed_x = cross_order_seed(args.order, reps, hn, hmap)
seed_desc = f"hint {args.hint} (order {hn})"
seed_bank = [np.ascontiguousarray(seed_x)]
seed_residual = B64 @ seed_x
print(
f"order={args.order} pairs={len(reps)} independent_lines={np.linalg.matrix_rank(B.astype(float))} "
f"seed={seed_desc} energy={int(seed_residual@seed_residual)} "
f"max_residual={int(np.max(np.abs(seed_residual)))} workers={args.workers} "
f"mode={'until-found' if args.until_found else 'finite'}"
)
if not np.array_equal(np.sort(np.abs(seed_x)), np.arange(1, len(reps)+1)):
raise RuntimeError("seed is not a signed permutation")
best_energy = int(seed_residual @ seed_residual)
best_x = seed_x.copy()
best_residual = seed_residual.copy()
deadline = math.inf if args.until_found else time.monotonic() + args.time_limit
# Compile the Numba kernel before forking workers.
_anneal(
seed_x,
rows,
coeff,
degree,
1,
args.seed,
1.0,
1.0,
1,
args.triple_probability,
args.quad_probability,
0,
-1,
0,
)
if not args.until_found:
deadline = time.monotonic() + args.time_limit
progress_interval = args.epoch if args.progress_epoch else args.progress_steps
if progress_interval < 0:
raise SystemExit("--progress-steps must be >= 0")
restart_index = 0
batch_index = 0
last_finisher_energy: int | None = None
finite_remaining = args.restarts
# Fork preserves the warmed Numba dispatcher on Linux. Fall back to the
# platform default context where fork is unavailable.
try:
context = mp.get_context("fork")
except ValueError:
context = mp.get_context()
pool = concurrent.futures.ProcessPoolExecutor(
max_workers=args.workers,
mp_context=context,
)
pool_stopped_hard = False
current_futures = {}
try:
while best_energy != 0:
if time.monotonic() >= deadline:
break
if not args.until_found and finite_remaining <= 0:
break
jobs = args.workers if args.until_found else min(args.workers, finite_remaining)
payloads = []
job_ids = []
incumbent_snapshot = best_x.copy()
for _ in range(jobs):
ridx = restart_index
restart_index += 1
if len(seed_bank) > 1 and (
ridx < len(seed_bank)
or (
args.potential_reset_every > 0
and ridx % args.potential_reset_every == 0
)
):
if ridx < len(seed_bank):
bank_index = ridx
else:
bank_index = (
ridx // args.potential_reset_every
) % len(seed_bank)
start_x = seed_bank[bank_index].copy()
elif ridx % 7 == 0:
start_x = seed_x.copy()
else:
start_x = perturb_state(
incumbent_snapshot,
args.seed + 65537 * ridx,
args.perturbations,
)
payloads.append(
(
start_x.astype(np.int64),
rows,
coeff,
degree,
args.steps,
args.seed + 1009 * ridx,
args.temperature_start,
args.temperature_end,
args.epoch,
args.triple_probability,
args.quad_probability,
progress_interval,
ridx,
args.acceptance_cache_size,
args.numba_perf,
)
)
job_ids.append(ridx)
# as_completed() is essential: submission-order future.result()
# can hide a solution already returned by a later worker.
future_to_ridx = {
pool.submit(_anneal_job, payload): ridx
for ridx, payload in zip(job_ids, payloads)
}
current_futures = future_to_ridx
for future in concurrent.futures.as_completed(future_to_ridx):
ridx = future_to_ridx[future]
energy, candidate_x, candidate_residual, used = future.result()
energy = int(energy)
if energy < best_energy:
best_energy = energy
best_x = candidate_x.copy()
best_residual = candidate_residual.copy()
if args.best_out is not None and best_energy != 0:
write_mhx(
args.best_out,
args.order,
reps,
best_x,
False,
best_residual,
)
print(
f"restart={ridx:6d} run_best={energy:8d} global_best={best_energy:8d} "
f"max_residual={int(np.max(np.abs(best_residual))):4d} steps={used}",
flush=True,
)
if energy == 0:
# Do not leave the executor through its context manager:
# that would perform shutdown(wait=True) and wait for every
# other worker to consume all of --steps.
print(
f"worker restart={ridx} found an exact solution after {used} steps; "
"terminating the remaining workers",
flush=True,
)
_stop_executor_now(pool, future_to_ridx)
pool_stopped_hard = True
break
current_futures = {}
if pool_stopped_hard:
break
if not args.until_found:
finite_remaining -= jobs
batch_index += 1
if best_energy == 0:
break
should_finish = (
args.finisher
and best_energy <= args.finisher_energy
and (
last_finisher_energy is None
or best_energy < last_finisher_energy
or batch_index % max(1, args.finisher_every) == 0
)
)
if should_finish:
print(
f"finisher: energy={best_energy} attempts={args.finisher_attempts} "
f"active_size={args.finisher_size} closure={args.finisher_closure}"
)
found, fx, fr, fstats = exact_lns_finisher(
B,
best_x,
best_residual,
attempts=args.finisher_attempts,
active_size=args.finisher_size,
closure_layers=args.finisher_closure,
time_limit=args.finisher_time_limit,
workers=args.finisher_workers or args.workers,
seed=args.seed + 1_000_003 * batch_index,
)
last_finisher_energy = best_energy
fenergy = int(fr @ fr)
if fenergy < best_energy:
best_energy = fenergy
best_x = fx.copy()
best_residual = fr.copy()
if args.best_out is not None and best_energy != 0:
write_mhx(
args.best_out,
args.order,
reps,
best_x,
False,
best_residual,
)
print(
f"finisher_result found={found} best_energy={best_energy} "
f"attempts={fstats['attempts']} active_sizes={fstats['active_sizes']}"
)
if found:
break
if args.until_found and batch_index % 10 == 0:
print(
f"campaign_checkpoint batches={batch_index} restarts={restart_index} "
f"best_energy={best_energy} residual={best_residual.tolist()}"
)
except BaseException:
# KeyboardInterrupt and unexpected worker failures must not strand a
# pool of long-running Numba jobs either.
if not pool_stopped_hard:
_stop_executor_now(pool, current_futures)
pool_stopped_hard = True
raise
finally:
if not pool_stopped_hard:
pool.shutdown(wait=True, cancel_futures=True)
exact = best_energy == 0
if exact:
out = args.out or Path(f"order{args.order}_r3_exchange.mhx")
write_mhx(out, args.order, reps, best_x, True, best_residual)
print(f"FOUND exact consecutive magic hexagon: {out}")
return 0
if args.best_out is not None:
write_mhx(args.best_out, args.order, reps, best_x, False, best_residual)
print(f"wrote best nonmagic candidate: {args.best_out}")
print(
f"no exact solution in budget; best energy={best_energy}, "
f"max_residual={int(np.max(np.abs(best_residual)))}, residual={best_residual.tolist()}"
)
return 1
if __name__ == "__main__":
raise SystemExit(main())
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment