mirror of
https://github.com/allaunthefox/Research-Stack.git
synced 2026-07-31 03:05:21 +00:00
allaun
7d3e90331a
- qaoa_adapter.py: fix finsler_metric_to_qubo to store Q_ij + Q_ji per undirected pair; update is_anisotropic and _find_anisotropic_pair to accept raw_matrix so asymmetry detection still works after summation; add measure option to pauli_to_cirq. - benchmark_finsler_qaoa.py: new benchmark harness for Finsler-Randers routing via QAOA.
2659 lines
92 KiB
Python
2659 lines
92 KiB
Python
"""
|
||
qaoa_adapter.py — Bidirectional QAOA Conversion Adapter Set
|
||
|
||
Borrows existing Lean-native problem representations and makes them
|
||
interactable via QAOA (Quantum Approximate Optimization Algorithm).
|
||
|
||
Forward: Problem → QUBO → Ising → Pauli strings → Cirq/Qiskit circuit
|
||
Backward: Measurements → bitstring → solution → problem update
|
||
|
||
Hooks into existing pipeline:
|
||
- EntropyMeasures.QUBOFormulation (Array Array Q16_16) → QUBO
|
||
- RotationQUBO.QUBOField (frustration, energyScale) → QUBO
|
||
- braid_search.build_qubo_matrix (bracket QUBO) → QUBO
|
||
- qubo_highs.solve_qubo_highs (SA / HiGHS) → stochastic solver
|
||
- Lean #eval bridge (via eigensolid_lean_bridge) → live Lean eval
|
||
|
||
Stochastic abuse: runs QAOA and classical stochastic solvers side-by-side
|
||
on the same QUBO for comparison.
|
||
|
||
All Q16_16 boundary conversion is explicit. No decision or gating logic.
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import json
|
||
import math
|
||
import time
|
||
from dataclasses import dataclass, field, asdict
|
||
from pathlib import Path
|
||
from typing import Any, Optional
|
||
|
||
import sys as _sys
|
||
_sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "lib"))
|
||
_sys.path.insert(0, str(Path(__file__).resolve().parent))
|
||
from q16 import Q16_SCALE, from_q16, to_q16
|
||
from braid_search import _q16_signed
|
||
|
||
# ── optional backend imports ──────────────────────────────────────────
|
||
|
||
try:
|
||
import cirq
|
||
_HAS_CIRQ = True
|
||
except ImportError:
|
||
_HAS_CIRQ = False
|
||
|
||
try:
|
||
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
|
||
_HAS_QISKIT = True
|
||
except ImportError:
|
||
_HAS_QISKIT = False
|
||
|
||
try:
|
||
import numpy as np
|
||
_HAS_NUMPY = True
|
||
except ImportError:
|
||
_HAS_NUMPY = False
|
||
|
||
|
||
# =========================================================================
|
||
# I. Data Models
|
||
# =========================================================================
|
||
|
||
@dataclass
|
||
class QUBO:
|
||
"""Quadratic Unconstrained Binary Optimization problem.
|
||
|
||
Minimize x^T Q x where x_i ∈ {0, 1}.
|
||
matrix[(i,j)] = coefficient for x_i * x_j (i <= j, upper triangular).
|
||
"""
|
||
n: int
|
||
matrix: dict[tuple[int, int], float] = field(default_factory=dict)
|
||
offset: float = 0.0
|
||
|
||
def __post_init__(self) -> None:
|
||
keys = list(self.matrix.keys())
|
||
for i, j in keys:
|
||
if i > j:
|
||
self.matrix[(j, i)] = self.matrix.pop((i, j))
|
||
|
||
def energy(self, x: list[int]) -> float:
|
||
e = self.offset
|
||
for (i, j), qij in self.matrix.items():
|
||
e += qij * x[i] * x[j]
|
||
return e
|
||
|
||
|
||
@dataclass
|
||
class Ising:
|
||
"""Ising Hamiltonian: H = Σ h_i s_i + Σ J_ij s_i s_j + offset.
|
||
|
||
s_i ∈ {+1, -1}. J uses upper triangular (i < j).
|
||
"""
|
||
n: int
|
||
h: list[float] = field(default_factory=list)
|
||
J: dict[tuple[int, int], float] = field(default_factory=dict)
|
||
offset: float = 0.0
|
||
|
||
def energy(self, s: list[int]) -> float:
|
||
if len(s) != self.n:
|
||
raise ValueError(f"expected {self.n} spins, got {len(s)}")
|
||
e = self.offset
|
||
for i in range(self.n):
|
||
e += self.h[i] * s[i]
|
||
for (i, j), Jij in self.J.items():
|
||
e += Jij * s[i] * s[j]
|
||
return e
|
||
|
||
|
||
@dataclass
|
||
class PauliSum:
|
||
"""Pauli string representation of an Ising Hamiltonian.
|
||
|
||
terms: list of (pauli_string, coefficient)
|
||
e.g. ("ZZII", 0.5), ("ZIII", -1.0)
|
||
offset: constant (identity) term
|
||
"""
|
||
n: int
|
||
terms: list[tuple[str, float]] = field(default_factory=list)
|
||
offset: float = 0.0
|
||
|
||
|
||
# =========================================================================
|
||
# II. Q16_16 Boundary Utilities
|
||
# =========================================================================
|
||
|
||
def q16_to_float(raw: int) -> float:
|
||
"""Convert Q16_16 raw integer to float (boundary only)."""
|
||
return from_q16(raw)
|
||
|
||
|
||
def float_to_q16(value: float) -> int:
|
||
"""Convert float to Q16_16 raw integer (boundary only)."""
|
||
return to_q16(value)
|
||
|
||
|
||
# =========================================================================
|
||
# III-A. Lean QUBOFormulation → QUBO
|
||
# Source: EntropyMeasures.lean — QUBOFormulation { matrix, numVariables }
|
||
# matrix is Array (Array Q16_16): outer = rows, inner = cols
|
||
# objective: E = Σ_i Σ_j Q_ij * x_i * x_j (both true → add)
|
||
# =========================================================================
|
||
|
||
def lean_qubo_formulation_to_qubo(
|
||
matrix: list[list[int]],
|
||
num_variables: Optional[int] = None,
|
||
offset_q16: int = 0,
|
||
) -> QUBO:
|
||
"""Convert a Lean QUBOFormulation matrix to a QUBO.
|
||
|
||
Args:
|
||
matrix: N×N array of Q16_16 raw integers (from Lean's
|
||
QUBOFormulation.matrix : Array (Array Q16_16))
|
||
num_variables: Number of binary variables (defaults to len(matrix))
|
||
offset_q16: Q16_16 constant term
|
||
|
||
Returns:
|
||
QUBO with float coefficients converted at the boundary.
|
||
"""
|
||
n = num_variables or len(matrix)
|
||
Q: dict[tuple[int, int], float] = {}
|
||
for i in range(min(n, len(matrix))):
|
||
row = matrix[i] if i < len(matrix) else []
|
||
for j in range(min(n, len(row))):
|
||
qij = row[j]
|
||
if qij != 0:
|
||
key = (min(i, j), max(i, j))
|
||
Q[key] = Q.get(key, 0.0) + q16_to_float(qij)
|
||
return QUBO(n=n, matrix=Q, offset=q16_to_float(offset_q16))
|
||
|
||
|
||
def qubo_to_lean_qubo_formulation(qubo: QUBO) -> list[list[int]]:
|
||
"""Convert a QUBO back to Lean QUBOFormulation matrix (Q16_16 ints).
|
||
|
||
Returns N×N list of Q16_16 raw integers, suitable for formatting
|
||
into Lean syntax::
|
||
|
||
QUBOFormulation.mk #[#[q00, q01, ...], #[q10, q11, ...], ...]
|
||
"""
|
||
n = qubo.n
|
||
matrix: list[list[int]] = [[0] * n for _ in range(n)]
|
||
for (i, j), qij in qubo.matrix.items():
|
||
raw = float_to_q16(qij)
|
||
matrix[i][j] = raw
|
||
if i != j:
|
||
matrix[j][i] = raw
|
||
return matrix
|
||
|
||
|
||
# =========================================================================
|
||
# III-B. Lean RotationQUBO.QUBOField → QUBO
|
||
# Source: RotationQUBO.lean — QUBOField { frustration, energyScale }
|
||
# fieldEnergy(x) = x² / (1 + δ²) - energyScale
|
||
# isFrustrated(x) = fieldEnergy > 0
|
||
# We discretize x across [-range, +range] into binary variables.
|
||
# =========================================================================
|
||
|
||
def lean_qubo_field_to_qubo(
|
||
frustration_raw: int,
|
||
energy_scale_raw: int,
|
||
n_vars: int = 8,
|
||
field_range: float = 2.0,
|
||
) -> QUBO:
|
||
"""Discretize a Lean RotationQUBO.QUBOField into a binary QUBO.
|
||
|
||
Maps the continuous field energy E(x) = x²/(1+δ²) - energyScale
|
||
onto n_vars binary variables by discretizing x ∈ [-range, +range].
|
||
|
||
The QUBO encodes: x = Σ_k s_k * Δ where s_k ∈ {0,1} and
|
||
Δ = 2*range / n_vars. The energy becomes a quadratic in the s_k.
|
||
|
||
Variables: x_i = 1 if discretization point i is selected (one-hot).
|
||
"""
|
||
δ = q16_to_float(frustration_raw)
|
||
ε = q16_to_float(energy_scale_raw)
|
||
denom = 1.0 + δ * δ
|
||
|
||
Q: dict[tuple[int, int], float] = {}
|
||
for i in range(n_vars):
|
||
xi = -field_range + (2.0 * field_range * i) / max(n_vars - 1, 1)
|
||
Ei = (xi * xi) / denom - ε
|
||
|
||
# One-hot diagonal: energy at this point
|
||
Q[(i, i)] = Ei
|
||
|
||
return QUBO(n=n_vars, matrix=Q)
|
||
|
||
|
||
# =========================================================================
|
||
# III-C. braid_search.build_qubo_matrix → QUBO
|
||
# Source: braid_search.py — build_qubo_matrix(brackets) → dict[(i,j), int]
|
||
# Each bracket has: lower, upper, gap, kappa, phi, admissible
|
||
# bracket_cost = -base + |gap|/10 (Q16_16)
|
||
# crossing_penalty = overlap * 2*Q16_ONE + diversity_reward
|
||
# =========================================================================
|
||
|
||
def braid_search_brackets_to_qubo(
|
||
brackets: list[dict],
|
||
as_q16: bool = False,
|
||
) -> QUBO:
|
||
"""Convert braid brackets to QUBO using the existing braid_search pipeline.
|
||
|
||
Args:
|
||
brackets: List of bracket dicts with lower/upper/gap/kappa/phi/admissible.
|
||
as_q16: If True, keep raw Q16_16 ints; if False, convert to float.
|
||
|
||
Returns:
|
||
QUBO problem (8-strand braid crossing selection).
|
||
"""
|
||
try:
|
||
_sys.path.insert(0, str(Path(__file__).resolve().parent))
|
||
from braid_search import build_qubo_matrix
|
||
except ImportError:
|
||
return _fallback_bracket_qubo(brackets)
|
||
|
||
q16_matrix = build_qubo_matrix(brackets)
|
||
n = max(max(i, j) for (i, j) in q16_matrix) + 1 if q16_matrix else 0
|
||
|
||
Q: dict[tuple[int, int], float] = {}
|
||
for (i, j), raw in q16_matrix.items():
|
||
if as_q16:
|
||
Q[(i, j)] = float(raw)
|
||
else:
|
||
Q[(i, j)] = q16_to_float(raw)
|
||
|
||
return QUBO(n=n, matrix=Q)
|
||
|
||
|
||
def _fallback_bracket_qubo(brackets: list[dict]) -> QUBO:
|
||
"""Fallback if braid_search is not importable."""
|
||
n = len(brackets)
|
||
Q: dict[tuple[int, int], float] = {}
|
||
for i, b in enumerate(brackets):
|
||
gap = q16_to_float(b.get("gap", 32768))
|
||
admissible = b.get("admissible", True)
|
||
base = 1.0 if admissible else 2.0
|
||
Q[(i, i)] = -base + abs(gap) * 0.1
|
||
for j in range(i + 1, n):
|
||
bj = brackets[j]
|
||
overlap = min(b.get("upper", 0), bj.get("upper", 0)) - \
|
||
max(b.get("lower", 0), bj.get("lower", 0))
|
||
if overlap > 0:
|
||
Q[(i, j)] = 2.0
|
||
return QUBO(n=n, matrix=Q)
|
||
|
||
|
||
# =========================================================================
|
||
# III-D. BraidReceipt → QUBO (8-strand braid)
|
||
# =========================================================================
|
||
|
||
def braid_receipt_to_qubo(receipt: dict) -> QUBO:
|
||
"""Convert a BraidReceipt JSON dict to a QUBO over 8 strand variables.
|
||
|
||
Variable mapping:
|
||
x_i = 1 if strand i is "active" (selected in crossing), 0 otherwise.
|
||
|
||
Cost terms:
|
||
- Diagonal: kappa_i (phase accumulation — lower is better)
|
||
- Off-diagonal: slot collision penalty if two strands share a slot
|
||
- Scar penalty: if scar_absent is false, each scarred strand adds cost
|
||
- Sidon slack bonus: want high slack = low max label used
|
||
- Residual penalty: non-converged strands increase cost
|
||
"""
|
||
braid = receipt.get("braid", receipt)
|
||
strands = braid.get("strands", [])
|
||
if not strands:
|
||
strands = _receipt_strands_from_bracket(braid)
|
||
|
||
n = len(strands)
|
||
Q: dict[tuple[int, int], float] = {}
|
||
offset = 0.0
|
||
|
||
slot_map: dict[int, list[int]] = {}
|
||
|
||
for i, s in enumerate(strands):
|
||
kappa = q16_to_float(s.get("kappa", 0))
|
||
phi = q16_to_float(s.get("phi", 0))
|
||
slot = s.get("slot", 0)
|
||
converged = s.get("converged", True)
|
||
residual = q16_to_float(s.get("residual", 0))
|
||
|
||
Q[(i, i)] = kappa + 0.1 * abs(phi)
|
||
|
||
if not converged and residual > 0:
|
||
Q[(i, i)] = Q[(i, i)] + 10.0 * residual
|
||
|
||
slot_map.setdefault(slot, []).append(i)
|
||
|
||
for slot, indices in slot_map.items():
|
||
if len(indices) > 1:
|
||
penalty = 50.0
|
||
for a in range(len(indices)):
|
||
for b in range(a + 1, len(indices)):
|
||
i, j = indices[a], indices[b]
|
||
Q[(i, j)] = Q.get((i, j), 0.0) + penalty
|
||
|
||
scar_absent = receipt.get("scar_absent", braid.get("scar_absent", True))
|
||
if not scar_absent:
|
||
for i in range(n):
|
||
Q[(i, i)] = Q.get((i, i), 0.0) + 100.0
|
||
|
||
sidon_slack = receipt.get("sidon_slack", 0)
|
||
max_label = 128 - sidon_slack
|
||
offset += max_label / 128.0
|
||
|
||
return QUBO(n=n, matrix=Q, offset=offset)
|
||
|
||
|
||
def _receipt_strands_from_bracket(bracket: dict) -> list[dict]:
|
||
"""Synthesize strand array from a minimal bracket-only receipt."""
|
||
default = bracket.get("bracket", bracket)
|
||
gaps = [default.get("gap", 32768)] * 8
|
||
kappas = [default.get("kappa", 0)] * 8
|
||
return [
|
||
{"id": chr(ord("a") + i), "slot": 1 << i, "kappa": kappas[i],
|
||
"phi": 0, "residual": 0, "converged": True}
|
||
for i in range(8)
|
||
]
|
||
|
||
|
||
# =========================================================================
|
||
# III-F. Bidirectional QAOA ↔ Greek HachimojiSubstitution Decoder
|
||
#
|
||
# Each of the 8 QUBO variables maps to a Greek Hachimoji state with a
|
||
# defined phase (45° steps), chirality, and flow direction per the
|
||
# Omindirection Logogram contract (OMINDIRECTION_LOGOGRAM_DESIGN_AND_COMPILER.md).
|
||
#
|
||
# Forward (LTR, phases 0-135°): Φ Λ Ρ Κ — normal Baker regime
|
||
# Backward (RTL, phases 180-315°): Ω Σ Π Ζ — quarantine/tearing regime
|
||
#
|
||
# Lean source: Semantics/HachimojiSubstitution.lean §5-6
|
||
# =========================================================================
|
||
|
||
# Ordered by bit index 0-7 (matches braid_receipt_to_qubo variable order)
|
||
GREEK_STATES: list[str] = ["Φ", "Λ", "Ρ", "Κ", "Ω", "Σ", "Π", "Ζ"]
|
||
|
||
GREEK_PHASE: dict[str, int] = {
|
||
"Φ": 0, "Λ": 45, "Ρ": 90, "Κ": 135,
|
||
"Ω": 180, "Σ": 225, "Π": 270, "Ζ": 315,
|
||
}
|
||
|
||
def _phase_to_chirality(phase: int) -> str:
|
||
"""Omindirection Principle 3: chirality is a projection of phase."""
|
||
if phase in (0, 180):
|
||
return "ambidextrous"
|
||
elif phase < 180:
|
||
return "left"
|
||
else:
|
||
return "right"
|
||
|
||
def _phase_to_direction(phase: int) -> str:
|
||
"""Phases 0-135 = forward (LTR); 180-315 = reverse (RTL)."""
|
||
return "forward" if phase < 180 else "reverse"
|
||
|
||
# Receipt Bool fields controlled by each Greek state
|
||
# (payloadBound, contradictionWitness, tearBoundary, detachedMass, residualLane)
|
||
_GREEK_RECEIPT_BITS: dict[str, tuple[bool, bool, bool, bool, bool]] = {
|
||
"Φ": (True, False, False, False, False),
|
||
"Λ": (True, False, False, False, False),
|
||
"Ρ": (False, False, False, False, True),
|
||
"Κ": (False, False, False, True, False),
|
||
"Ω": (True, True, True, True, True),
|
||
"Σ": (False, False, True, False, False),
|
||
"Π": (False, False, False, False, False),
|
||
"Ζ": (False, False, False, True, False),
|
||
}
|
||
|
||
def _greek_to_regime(state: str) -> str:
|
||
"""SemanticRegime for a Greek state."""
|
||
if state in ("Φ", "Λ"):
|
||
return "beautifulTopologicalFolding"
|
||
elif state in ("Ρ", "Κ"):
|
||
return "uglyAsymmetricPruning"
|
||
else:
|
||
return "horribleManifoldTearing"
|
||
|
||
|
||
def decode_bitstring_to_greek_states(bits: list[bool]) -> list[dict]:
|
||
"""Backward QAOA decoder: measurement bitstring → Greek state atoms.
|
||
|
||
Each active bit (True) becomes an omindirectional atom with:
|
||
symbol, phase, chirality, direction, regime.
|
||
|
||
Matches Lean: HachimojiSubstitution.bitToGreek / Greek.HachimojiBase.phase
|
||
"""
|
||
atoms = []
|
||
for i, active in enumerate(bits[:8]):
|
||
if active:
|
||
sym = GREEK_STATES[i]
|
||
phase = GREEK_PHASE[sym]
|
||
atoms.append({
|
||
"bit": i,
|
||
"symbol": sym,
|
||
"phase": phase,
|
||
"chirality": _phase_to_chirality(phase),
|
||
"direction": _phase_to_direction(phase),
|
||
"regime": _greek_to_regime(sym),
|
||
})
|
||
return atoms
|
||
|
||
|
||
def greek_states_to_logogram_receipt(active_states: list[dict]) -> dict:
|
||
"""Backward QAOA decoder: Greek atom list → LogogramReceipt dict.
|
||
|
||
Dominant state = lowest phase (most stable).
|
||
Bool fields accumulated from ALL active states (OR-fold).
|
||
Matches Lean: HachimojiSubstitution.fromQAOABitstring
|
||
|
||
The receipt dict is compatible with RRCLogogramProjection.LogogramReceipt.
|
||
"""
|
||
if not active_states:
|
||
return {
|
||
"shape": "logogramProjection",
|
||
"status": "hold",
|
||
"regime": "uglyAsymmetricPruning",
|
||
"payloadBound": False,
|
||
"contradictionWitness": False,
|
||
"tearBoundary": False,
|
||
"detachedMass": False,
|
||
"residualLane": False,
|
||
}
|
||
|
||
# Dominant = lowest phase
|
||
dominant = min(active_states, key=lambda a: a["phase"])
|
||
|
||
# Accumulate Bool fields from all active states
|
||
pb = any(_GREEK_RECEIPT_BITS[a["symbol"]][0] for a in active_states)
|
||
cw = any(_GREEK_RECEIPT_BITS[a["symbol"]][1] for a in active_states)
|
||
tb = any(_GREEK_RECEIPT_BITS[a["symbol"]][2] for a in active_states)
|
||
dm = any(_GREEK_RECEIPT_BITS[a["symbol"]][3] for a in active_states)
|
||
rl = any(_GREEK_RECEIPT_BITS[a["symbol"]][4] for a in active_states)
|
||
|
||
# status = hold if Ζ (bit 7) is active
|
||
zeta_active = any(a["symbol"] == "Ζ" for a in active_states)
|
||
|
||
return {
|
||
"shape": "logogramProjection",
|
||
"status": "hold" if zeta_active else "candidate",
|
||
"regime": dominant["regime"],
|
||
"payloadBound": pb,
|
||
"contradictionWitness": cw,
|
||
"tearBoundary": tb,
|
||
"detachedMass": dm,
|
||
"residualLane": rl,
|
||
# Omindirection metadata
|
||
"chirality_witness": {
|
||
"direction": dominant["direction"],
|
||
"handedness": dominant["chirality"],
|
||
"phase": dominant["phase"],
|
||
"placement": "quarantine" if dominant["direction"] == "reverse" else "row",
|
||
"mode": "greek_hachimoji_v1",
|
||
},
|
||
}
|
||
|
||
|
||
def logogram_receipt_to_bitstring(receipt: dict) -> list[bool]:
|
||
"""RTL direction: LogogramReceipt → 8-bit Greek signature.
|
||
|
||
Matches Lean: HachimojiSubstitution.toQAOABitstring
|
||
This is the backward path for re-encoding into a QUBO update.
|
||
"""
|
||
regime = receipt.get("regime", "")
|
||
return [
|
||
bool(receipt.get("payloadBound", False)), # Φ
|
||
regime == "beautifulTopologicalFolding", # Λ
|
||
bool(receipt.get("residualLane", False)), # Ρ
|
||
bool(receipt.get("detachedMass", False)), # Κ
|
||
bool(receipt.get("contradictionWitness", False)), # Ω
|
||
bool(receipt.get("tearBoundary", False)), # Σ
|
||
regime == "horribleManifoldTearing", # Π
|
||
receipt.get("status", "") == "hold", # Ζ
|
||
]
|
||
|
||
|
||
def qaoa_measurement_to_receipt(bits: list[bool]) -> dict:
|
||
"""Full bidi decoder: QAOA bitstring → LogogramReceipt.
|
||
|
||
Combines decode_bitstring_to_greek_states + greek_states_to_logogram_receipt.
|
||
This is the primary backward path entry point.
|
||
"""
|
||
atoms = decode_bitstring_to_greek_states(bits)
|
||
return greek_states_to_logogram_receipt(atoms)
|
||
|
||
|
||
# =========================================================================
|
||
# IV. Forward: LonelyRunner β₀/Scar → QUBO
|
||
# =========================================================================
|
||
|
||
def lonely_runner_scar_to_qubo(
|
||
scar_field: Any,
|
||
*,
|
||
beta0_weight: float = 10.0,
|
||
coverage_weight: float = 1.0,
|
||
) -> QUBO:
|
||
"""Convert a LonelyRunner scar field to a QUBO over N circle points.
|
||
|
||
The scar field μ(t,θ) ∈ {0, 1} indicates uncovered (lonely) regions.
|
||
|
||
Variables: x_i = scarred[i] (1 = uncovered, 0 = covered), N points.
|
||
|
||
Cost: maximize β₀ (connected components of uncovered region)
|
||
= minimize -β₀ + coverage_penalty
|
||
|
||
β₀ = Σ_i x_i * (1 - x_{i-1}) (rising edges, cyclically)
|
||
|
||
For an end-to-end QAOA that finds a lonely time:
|
||
H_cost = -β₀ + λ * coverage
|
||
The ground state = max β₀ at minimal coverage = a lonely time.
|
||
"""
|
||
if isinstance(scar_field, list) and _HAS_NUMPY:
|
||
import numpy as _np
|
||
scar_field = _np.array(scar_field)
|
||
|
||
if _HAS_NUMPY:
|
||
import numpy as _np
|
||
T, N = scar_field.shape
|
||
avg_scar = _np.mean(scar_field, axis=0) > 0.5
|
||
else:
|
||
T = len(scar_field)
|
||
N = len(scar_field[0]) if T > 0 else 0
|
||
avg_scar = [sum(scar_field[t][i] for t in range(T)) / T > 0.5
|
||
for i in range(N)]
|
||
|
||
Q: dict[tuple[int, int], float] = {}
|
||
|
||
for i in range(N):
|
||
prev = (i - 1) % N
|
||
Q[(i, i)] = Q.get((i, i), 0.0) - beta0_weight
|
||
Q[(prev, i)] = Q.get((prev, i), 0.0) + beta0_weight
|
||
|
||
for i in range(N):
|
||
bias = -coverage_weight if avg_scar[i] > 0.5 else coverage_weight
|
||
Q[(i, i)] = Q.get((i, i), 0.0) + bias
|
||
|
||
return QUBO(n=N, matrix=Q)
|
||
|
||
|
||
# =========================================================================
|
||
# V. Forward: Goormaghtigh Cost → QUBO
|
||
# =========================================================================
|
||
|
||
def goormaghtigh_cost_to_qubo(
|
||
x_range: tuple[int, int] = (2, 90),
|
||
m_range: tuple[int, int] = (3, 13),
|
||
penalty_weight: float = 1000.0,
|
||
) -> QUBO:
|
||
"""Convert the Goormaghtigh repunit collision problem to a QUBO.
|
||
|
||
Repunit: (x^m - 1)/(x - 1) = 1 + x + ... + x^{m-1}
|
||
|
||
Goal: find (x₁,m₁) ≠ (x₂,m₂) s.t. repunit(x₁,m₁) = repunit(x₂,m₂).
|
||
|
||
QUBO encoding:
|
||
- One-hot encoding for x and m in each repunit
|
||
- Collision penalty proportional to repunit difference
|
||
"""
|
||
x_vals = list(range(x_range[0], x_range[1] + 1))
|
||
m_vals = list(range(m_range[0], m_range[1] + 1))
|
||
nx, nm = len(x_vals), len(m_vals)
|
||
|
||
off_x1 = 0
|
||
off_m1 = nx
|
||
off_x2 = nx + nm
|
||
off_m2 = nx + nm + nx
|
||
|
||
n = nx + nm + nx + nm
|
||
Q: dict[tuple[int, int], float] = {}
|
||
|
||
def _idx(offset: int, k: int) -> int:
|
||
return offset + k
|
||
|
||
for offset in (off_x1, off_x2):
|
||
for i in range(nx):
|
||
Q[(_idx(offset, i), _idx(offset, i))] = \
|
||
Q.get((_idx(offset, i), _idx(offset, i)), 0.0) - penalty_weight
|
||
for i in range(nx):
|
||
for j in range(i + 1, nx):
|
||
Q[(_idx(offset, i), _idx(offset, j))] = \
|
||
Q.get((_idx(offset, i), _idx(offset, j)), 0.0) + 2.0 * penalty_weight
|
||
|
||
for offset in (off_m1, off_m2):
|
||
for i in range(nm):
|
||
Q[(_idx(offset, i), _idx(offset, i))] = \
|
||
Q.get((_idx(offset, i), _idx(offset, i)), 0.0) - penalty_weight
|
||
for i in range(nm):
|
||
for j in range(i + 1, nm):
|
||
Q[(_idx(offset, i), _idx(offset, j))] = \
|
||
Q.get((_idx(offset, i), _idx(offset, j)), 0.0) + 2.0 * penalty_weight
|
||
|
||
for xi in range(nx):
|
||
for mi in range(nm):
|
||
r1 = _repunit(x_vals[xi], m_vals[mi])
|
||
for xj in range(nx):
|
||
for mj in range(nm):
|
||
if xi == xj and mi == mj:
|
||
continue
|
||
r2 = _repunit(x_vals[xj], m_vals[mj])
|
||
diff = abs(r1 - r2)
|
||
i1 = _idx(off_x1, xi)
|
||
i2 = _idx(off_m1, mi)
|
||
i3 = _idx(off_x2, xj)
|
||
i4 = _idx(off_m2, mj)
|
||
w = penalty_weight / max(diff, 1.0)
|
||
Q[(i1, i3)] = Q.get((i1, i3), 0.0) + w * 0.25
|
||
Q[(i2, i4)] = Q.get((i2, i4), 0.0) + w * 0.25
|
||
|
||
return QUBO(n=n, matrix=Q)
|
||
|
||
|
||
def _repunit(x: int, m: int) -> int:
|
||
"""Compute repunit value (x^m - 1)//(x - 1)."""
|
||
return (x ** m - 1) // (x - 1)
|
||
|
||
|
||
# =========================================================================
|
||
# V-A. Tunable Parameterization of Existing Cost Functions
|
||
#
|
||
# The QUBO cost functions in braid_search.py use hardcoded coefficients.
|
||
# These wrappers let you experiment with alternative parameter sets so
|
||
# QAOA results can inform better classical cost models.
|
||
#
|
||
# DefaultCoeffs mirrors the hardcoded values in braid_search.py:
|
||
# - bracket_cost: -base + |gap| * gap_reward
|
||
# - crossing_penalty: overlap_penalty - |gap_diff| * gap_reward
|
||
# =========================================================================
|
||
|
||
@dataclass
|
||
class BraidCostCoeffs:
|
||
"""Tunable coefficients for braid bracket cost functions.
|
||
|
||
Mirrors the hardcoded parameters from braid_search.py so any
|
||
of them can be varied independently for grid-search tuning.
|
||
|
||
Defaults calibrated via SLOS photonic emulator 2026-06-18.
|
||
See qaoa_adapter.auto_calibrate() for re-calibration.
|
||
"""
|
||
base_admissible: float = 8.0 # 8*Q16_ONE — linear bias for admissible brackets
|
||
base_inadmissible: float = 16.0 # 16*Q16_ONE — penalty for inadmissible
|
||
gap_reward_scale: float = 0.01 # reward for large gaps
|
||
overlap_penalty: float = 2.0 # 2*Q16_ONE — conflict cost when ranges overlap
|
||
sa_initial_temp: float = 10.0 # SA starting temperature
|
||
sa_cooling_rate: float = 0.9995 # SA geometric cooling factor
|
||
sa_iterations: int = 5000 # SA iteration budget
|
||
slot_collision_penalty: float = 50.0 # receipt slot conflict cost
|
||
scar_absent_penalty: float = 100.0 # scar-present penalty
|
||
levy_alpha: float = 1.5 # Lévy flight exponent
|
||
|
||
|
||
DEFAULT_COEFFS = BraidCostCoeffs()
|
||
|
||
|
||
def tunable_bracket_cost(bracket: dict, coeffs: BraidCostCoeffs = DEFAULT_COEFFS) -> float:
|
||
"""Parameterized replacement for braid_search.bracket_cost.
|
||
|
||
Cost = -base + |gap| * gap_reward_scale (float, not Q16_16)
|
||
"""
|
||
admissible = bracket.get("admissible", True)
|
||
base = coeffs.base_admissible if admissible else coeffs.base_inadmissible
|
||
gap = _q16_signed(bracket.get("gap", 32768))
|
||
gap_float = gap / 65536.0
|
||
return -base + abs(gap_float) * coeffs.gap_reward_scale
|
||
|
||
|
||
def tunable_crossing_penalty(
|
||
b1: dict, b2: dict, coeffs: BraidCostCoeffs = DEFAULT_COEFFS,
|
||
) -> float:
|
||
"""Parameterized replacement for braid_search.crossing_penalty."""
|
||
cost = 0.0
|
||
# Overlap check
|
||
l1 = b1.get("lower", 0)
|
||
u1 = b1.get("upper", 0)
|
||
l2 = b2.get("lower", 0)
|
||
u2 = b2.get("upper", 0)
|
||
if l1 < u2 and l2 < u1:
|
||
cost += coeffs.overlap_penalty
|
||
# Gap diversity reward
|
||
g1 = _q16_signed(b1.get("gap", 32768))
|
||
g2 = _q16_signed(b2.get("gap", 32768))
|
||
gap_diff = abs(g1 - g2) / 65536.0
|
||
cost -= gap_diff * coeffs.gap_reward_scale
|
||
return cost
|
||
|
||
|
||
def build_tunable_qubo_matrix(
|
||
brackets: list[dict],
|
||
coeffs: BraidCostCoeffs = DEFAULT_COEFFS,
|
||
) -> dict[tuple[int, int], float]:
|
||
"""Build QUBO dict using tunable coefficients.
|
||
|
||
Returns float-valued matrix (not Q16_16), directly consumable by
|
||
stochastic_abuse_qubo and qaoa_solve_qubo.
|
||
"""
|
||
Q: dict[tuple[int, int], float] = {}
|
||
n = len(brackets)
|
||
for i in range(n):
|
||
Q[(i, i)] = tunable_bracket_cost(brackets[i], coeffs)
|
||
for j in range(i + 1, n):
|
||
c = tunable_crossing_penalty(brackets[i], brackets[j], coeffs)
|
||
Q[(i, j)] = c
|
||
Q[(j, i)] = c
|
||
return Q
|
||
|
||
|
||
def tunable_braid_receipt_to_qubo(
|
||
receipt: dict,
|
||
coeffs: BraidCostCoeffs = DEFAULT_COEFFS,
|
||
) -> QUBO:
|
||
"""braid_receipt_to_qubo with tunable slot/scar penalties."""
|
||
braid = receipt.get("braid", receipt)
|
||
strands = braid.get("strands", [])
|
||
if not strands:
|
||
strands = _receipt_strands_from_bracket(braid)
|
||
n = len(strands)
|
||
Q: dict[tuple[int, int], float] = {}
|
||
offset = 0.0
|
||
|
||
slot_map: dict[int, list[int]] = {}
|
||
for i, s in enumerate(strands):
|
||
kappa = s.get("kappa", 0) / 65536.0
|
||
phi = s.get("phi", 0) / 65536.0
|
||
slot = s.get("slot", 0)
|
||
converged = s.get("converged", True)
|
||
residual = s.get("residual", 0) / 65536.0
|
||
|
||
Q[(i, i)] = kappa + 0.1 * abs(phi)
|
||
if not converged and residual > 0:
|
||
Q[(i, i)] += 10.0 * residual
|
||
slot_map.setdefault(slot, []).append(i)
|
||
|
||
for _slot, indices in slot_map.items():
|
||
if len(indices) > 1:
|
||
for a in range(len(indices)):
|
||
for b in range(a + 1, len(indices)):
|
||
i, j = indices[a], indices[b]
|
||
p = coeffs.slot_collision_penalty
|
||
Q[(i, j)] = Q.get((i, j), 0.0) + p
|
||
|
||
scar_absent = receipt.get("scar_absent", braid.get("scar_absent", True))
|
||
if not scar_absent:
|
||
for i in range(n):
|
||
Q[(i, i)] = Q.get((i, i), 0.0) + coeffs.scar_absent_penalty
|
||
|
||
sidon_slack = receipt.get("sidon_slack", 0)
|
||
max_label = 128 - sidon_slack
|
||
offset += max_label / 128.0
|
||
|
||
return QUBO(n=n, matrix=Q, offset=offset)
|
||
|
||
|
||
# =========================================================================
|
||
# V-B. Tuner: Grid Search Over Parameters
|
||
# =========================================================================
|
||
|
||
def tune_braid_parameters(
|
||
brackets: list[dict],
|
||
param_grid: Optional[list[dict]] = None,
|
||
solver_methods: Optional[list[str]] = None,
|
||
time_per_trial: float = 1.0,
|
||
seed: int = 42,
|
||
) -> dict:
|
||
"""Grid search over braid cost parameters to find the best QUBO formulation.
|
||
|
||
For each parameter set in the grid, builds the QUBO and solves it with
|
||
each solver. Returns the configuration that yields the lowest energy,
|
||
plus a full comparison table.
|
||
|
||
Args:
|
||
brackets: List of bracket dicts.
|
||
param_grid: List of BraidCostCoeffs overrides.
|
||
Each entry is a dict with keys from BraidCostCoeffs fields.
|
||
Default: 27 combinations varying base_admissible, overlap_penalty,
|
||
gap_reward_scale.
|
||
solver_methods: Solver methods to evaluate.
|
||
Default: ["sa", "levy"].
|
||
time_per_trial: Time limit per solver trial (seconds).
|
||
seed: RNG seed.
|
||
|
||
Returns:
|
||
{
|
||
"best_params": {...}, # winning BraidCostCoeffs
|
||
"best_solver": str, # winning solver name
|
||
"best_energy": float, # lowest energy found
|
||
"trials": [ # each trial result
|
||
{
|
||
"params": {...},
|
||
"solver": str,
|
||
"energy": float,
|
||
"solution": list[int],
|
||
"runtime_s": float,
|
||
}
|
||
],
|
||
"n": int, # number of brackets
|
||
}
|
||
"""
|
||
if param_grid is None:
|
||
param_grid = []
|
||
for base_mult in [0.5, 1.0, 2.0]:
|
||
for overlap_mult in [0.5, 1.0, 2.0]:
|
||
for gap_scale in [0.05, 0.1, 0.2]:
|
||
param_grid.append({
|
||
"base_admissible": base_mult,
|
||
"base_inadmissible": 2.0 * base_mult,
|
||
"overlap_penalty": 2.0 * overlap_mult,
|
||
"gap_reward_scale": gap_scale,
|
||
"slot_collision_penalty": 50.0,
|
||
"scar_absent_penalty": 100.0,
|
||
"sa_initial_temp": 10.0,
|
||
"sa_cooling_rate": 0.9995,
|
||
"sa_iterations": 5000,
|
||
"levy_alpha": 1.5,
|
||
})
|
||
|
||
if solver_methods is None:
|
||
solver_methods = ["sa", "levy"]
|
||
|
||
import copy
|
||
trials: list[dict] = []
|
||
|
||
for param_dict in param_grid:
|
||
coeffs = BraidCostCoeffs(**param_dict)
|
||
Q = build_tunable_qubo_matrix(brackets, coeffs)
|
||
n = len(brackets)
|
||
qubo = QUBO(n=n, matrix=Q)
|
||
|
||
for method in solver_methods:
|
||
t0 = time.time()
|
||
result = stochastic_abuse_qubo(
|
||
qubo, method=method, time_limit=time_per_trial, seed=seed,
|
||
)
|
||
runtime = time.time() - t0
|
||
|
||
trials.append({
|
||
"params": param_dict,
|
||
"solver": method,
|
||
"energy": result["energy"],
|
||
"solution": result["solution"],
|
||
"runtime_s": round(runtime, 4),
|
||
})
|
||
|
||
# Find best trial (lowest energy)
|
||
best_trial = min(trials, key=lambda t: t["energy"])
|
||
|
||
return {
|
||
"best_params": best_trial["params"],
|
||
"best_solver": best_trial["solver"],
|
||
"best_energy": best_trial["energy"],
|
||
"trials": trials,
|
||
"n": len(brackets),
|
||
}
|
||
|
||
|
||
# =========================================================================
|
||
# V-C. SA Parameter Tuner
|
||
# =========================================================================
|
||
|
||
def tune_sa_parameters(
|
||
qubo: QUBO,
|
||
temp_grid: Optional[list[float]] = None,
|
||
cooling_grid: Optional[list[float]] = None,
|
||
time_limit: float = 2.0,
|
||
seed: int = 42,
|
||
) -> dict:
|
||
"""Grid search over SA temperature and cooling rate parameters.
|
||
|
||
Runs SA with each (initial_temp, cooling_rate) combination and reports
|
||
the configuration that yields the lowest energy.
|
||
|
||
Args:
|
||
qubo: QUBO problem.
|
||
temp_grid: List of initial temperatures to try.
|
||
Default: [1.0, 5.0, 10.0, 20.0, 50.0].
|
||
cooling_grid: List of cooling rates to try.
|
||
Default: [0.99, 0.995, 0.999, 0.9995, 0.9999].
|
||
time_limit: Time limit per SA run.
|
||
seed: RNG seed.
|
||
|
||
Returns:
|
||
{
|
||
"best_temp": float,
|
||
"best_cooling": float,
|
||
"best_energy": float,
|
||
"trials": [{"temp": float, "cooling": float, "energy": float}],
|
||
}
|
||
"""
|
||
if temp_grid is None:
|
||
temp_grid = [1.0, 5.0, 10.0, 20.0, 50.0]
|
||
if cooling_grid is None:
|
||
cooling_grid = [0.99, 0.995, 0.999, 0.9995, 0.9999]
|
||
|
||
trials: list[dict] = []
|
||
Q_dict: dict[tuple[int, int], float] = dict(qubo.matrix)
|
||
n = qubo.n
|
||
|
||
for temp in temp_grid:
|
||
for cooling in cooling_grid:
|
||
t0 = time.time()
|
||
x = [0] * n
|
||
import random as _r
|
||
_r.seed(seed)
|
||
x = [_r.randint(0, 1) for _ in range(n)]
|
||
|
||
def _energy(xv):
|
||
e = qubo.offset
|
||
for (i, j), qij in Q_dict.items():
|
||
e += qij * xv[i] * xv[j]
|
||
return e
|
||
|
||
current_energy = _energy(x)
|
||
best_x = list(x)
|
||
best_energy = current_energy
|
||
T = temp
|
||
|
||
while time.time() - t0 < time_limit:
|
||
i = _r.randrange(n)
|
||
x[i] = 1 - x[i]
|
||
new_energy = _energy(x)
|
||
delta = new_energy - current_energy
|
||
if delta < 0 or _r.random() < math.exp(-delta / max(T, 1e-10)):
|
||
current_energy = new_energy
|
||
if current_energy < best_energy:
|
||
best_x = list(x)
|
||
best_energy = current_energy
|
||
else:
|
||
x[i] = 1 - x[i]
|
||
T *= cooling
|
||
|
||
trials.append({
|
||
"temp": temp,
|
||
"cooling": cooling,
|
||
"energy": best_energy,
|
||
"runtime_s": round(time.time() - t0, 4),
|
||
})
|
||
|
||
best_trial = min(trials, key=lambda t: t["energy"])
|
||
return {
|
||
"best_temp": best_trial["temp"],
|
||
"best_cooling": best_trial["cooling"],
|
||
"best_energy": best_trial["energy"],
|
||
"trials": trials,
|
||
}
|
||
|
||
|
||
# =========================================================================
|
||
# V-D. Solver Dispatch Tuner: Auto-pick best solver for a QUBO
|
||
# =========================================================================
|
||
|
||
def solver_dispatch_tuner(
|
||
qubo: QUBO,
|
||
time_limit: float = 2.0,
|
||
seed: int = 42,
|
||
) -> dict:
|
||
"""Auto-select the best solver for a given QUBO instance.
|
||
|
||
Runs all available solvers (SA, HiGHS, Lévy, QAOA describe) for a
|
||
brief warmup and returns the one with the lowest energy, plus a
|
||
recommendation based on problem size.
|
||
|
||
Args:
|
||
qubo: QUBO problem.
|
||
time_limit: Per-solver time limit.
|
||
seed: RNG seed.
|
||
|
||
Returns:
|
||
{
|
||
"recommended": str, # solver name
|
||
"recommended_energy": float,
|
||
"warmup": {method: {"energy": float, "runtime_s": float}},
|
||
"heuristic_note": str, # size-based guidance
|
||
}
|
||
"""
|
||
results: dict[str, dict] = {}
|
||
|
||
# SA
|
||
sa_res = stochastic_abuse_qubo(qubo, method="sa", time_limit=time_limit, seed=seed)
|
||
results["sa"] = {"energy": sa_res["energy"], "runtime_s": sa_res["runtime_s"]}
|
||
|
||
# Lévy
|
||
try:
|
||
levy_res = stochastic_abuse_qubo(qubo, method="levy", time_limit=time_limit, seed=seed)
|
||
results["levy"] = {"energy": levy_res["energy"], "runtime_s": levy_res["runtime_s"]}
|
||
except Exception:
|
||
results["levy"] = {"energy": float("inf"), "runtime_s": 0.0}
|
||
|
||
# HiGHS (may fall back to SA internally)
|
||
try:
|
||
highs_res = stochastic_abuse_qubo(qubo, method="highs", time_limit=time_limit, seed=seed)
|
||
results["highs"] = {"energy": highs_res["energy"], "runtime_s": highs_res["runtime_s"]}
|
||
except Exception:
|
||
results["highs"] = {"energy": float("inf"), "runtime_s": 0.0}
|
||
|
||
# QAOA describe (offline estimate)
|
||
qaoa_res = qaoa_solve_qubo(qubo, p_layers=1, shots=100, backend="describe")
|
||
results["qaoa"] = {"energy": qaoa_res["energy"], "runtime_s": 0.0, "note": "describe-backend-estimate"}
|
||
|
||
# Pick best
|
||
best_solver = min(results, key=lambda s: results[s]["energy"])
|
||
best_energy = results[best_solver]["energy"]
|
||
|
||
# Heuristic
|
||
n = qubo.n
|
||
if n <= 5:
|
||
heuristic = "small: HiGHS MIP typically exact and fast"
|
||
elif n <= 20:
|
||
heuristic = "medium: SA or Lévy often best; QAOA may help for specific structure"
|
||
else:
|
||
heuristic = "large: SA or Lévy cheaper; QAOA promising if coupling structure is sparse"
|
||
|
||
return {
|
||
"recommended": best_solver,
|
||
"recommended_energy": best_energy,
|
||
"warmup": results,
|
||
"heuristic_note": heuristic,
|
||
"n": n,
|
||
}
|
||
|
||
|
||
# =========================================================================
|
||
# V-E. SLOS Photonic Cost Calibration
|
||
#
|
||
# Uses the Perceval SLOS (Schrödinger Levitated Object Simulator) to
|
||
# compute photonic reference costs for bracket configurations. The SLOS
|
||
# output distribution over photonic modes serves as a physically-grounded
|
||
# cost signal. We then fit the classical BraidCostCoeffs to match these
|
||
# reference costs — effectively "calibrating" the heuristic defaults to
|
||
# match what a photonic emulator says.
|
||
#
|
||
# Requires: perceval (pip install perceval-quandela)
|
||
# =========================================================================
|
||
|
||
try:
|
||
import perceval as _pcvl
|
||
_HAS_PERVERSE = True
|
||
except ImportError:
|
||
_HAS_PERVERSE = False
|
||
|
||
|
||
class SlosBraidEvaluator:
|
||
"""Evaluate braid brackets using Perceval SLOS photonic simulation.
|
||
|
||
Encodes bracket parameters (gap, admissible, lower, upper) as
|
||
phase angles in an M-mode linear optical interferometer. The SLOS
|
||
statevector simulator computes the exact output distribution; we
|
||
extract cost signals from mode probabilities.
|
||
|
||
The key insight: admissible brackets with large gaps produce
|
||
different photonic interference patterns (higher mode spread,
|
||
lower entropy) than inadmissible or overlapping ones. SLOS
|
||
captures this physically rather than heuristically.
|
||
"""
|
||
|
||
def __init__(self, M: int = 6):
|
||
self.M = M
|
||
self._cache: dict[str, float] = {}
|
||
|
||
def _cache_key(self, bracket: dict) -> str:
|
||
return f"{bracket.get('lower',0)}|{bracket.get('upper',0)}|{bracket.get('gap',0)}|{bracket.get('admissible',True)}"
|
||
|
||
def _encode_bracket_angles(self, bracket: dict) -> list[float]:
|
||
"""Map bracket parameters to photonic phase angles.
|
||
|
||
Encoding scheme:
|
||
- admissible → phase = 0.0 (no shift — clean path)
|
||
- inadmissible → phase = π/2 (rotation — scattering loss)
|
||
- gap → proportional phase scale (0..π), larger gap = larger shift
|
||
- lower/upper → differential phase between modes
|
||
"""
|
||
admissible = bracket.get("admissible", True)
|
||
gap = bracket.get("gap", 32768)
|
||
lower = bracket.get("lower", 0)
|
||
upper = bracket.get("upper", 65536)
|
||
|
||
# Normalize to 0..1
|
||
gap_norm = min(1.0, abs(gap) / 65536.0)
|
||
span_norm = min(1.0, max(0, (upper - lower)) / 65536.0)
|
||
|
||
theta: list[float] = []
|
||
# Mode 0: admissibility phase
|
||
theta.append(0.0 if admissible else math.pi / 2.0)
|
||
# Mode 1: gap phase
|
||
theta.append(gap_norm * math.pi)
|
||
# Mode 2: span phase
|
||
theta.append(span_norm * math.pi)
|
||
# Modes 3+: fill with differential phases
|
||
for k in range(3, self.M):
|
||
theta.append((gap_norm * span_norm * math.pi) / (1.0 + k))
|
||
return theta
|
||
|
||
def compute_photonic_cost(self, bracket: dict) -> float:
|
||
"""Run SLOS on a single bracket and return photonic cost.
|
||
|
||
Lower photonic cost = more favorable bracket configuration.
|
||
The cost is derived from the output mode distribution:
|
||
- Admissible + large gap → concentrated output (low cost)
|
||
- Inadmissible + small gap → scattered output (high cost)
|
||
"""
|
||
ck = self._cache_key(bracket)
|
||
if ck in self._cache:
|
||
return self._cache[ck]
|
||
|
||
if not _HAS_PERVERSE:
|
||
# Fallback: heuristic estimate mirrors braid_search defaults
|
||
admissible = bracket.get("admissible", True)
|
||
gap = bracket.get("gap", 32768)
|
||
base = 1.0 if admissible else 2.0
|
||
cost = base + gap / 65536.0 * 0.5
|
||
self._cache[ck] = cost
|
||
return cost
|
||
|
||
theta = self._encode_bracket_angles(bracket)
|
||
|
||
try:
|
||
circuit = _pcvl.Circuit(self.M)
|
||
for i in range(min(len(theta), self.M)):
|
||
circuit.add(i, _pcvl.PS(theta[i]))
|
||
for i in range(self.M - 1):
|
||
circuit.add((i, i + 1), _pcvl.BS())
|
||
|
||
input_state = _pcvl.BasicState([1] + [0] * (self.M - 1))
|
||
processor = _pcvl.Processor("SLOS", circuit)
|
||
processor.with_input(input_state)
|
||
sampler = _pcvl.algorithm.Sampler(processor)
|
||
res = sampler.sample_count(1000)
|
||
|
||
# Compute photonic cost from output distribution
|
||
# Metric: weighted entropy of output modes
|
||
total_prob = 0.0
|
||
entropy = 0.0
|
||
for state, count in res["results"].items():
|
||
prob = count / 1000.0
|
||
total_prob += prob
|
||
if prob > 1e-10:
|
||
entropy -= prob * math.log2(prob)
|
||
|
||
# Scale: admissible brackets have lower entropy (cleaner interference)
|
||
admissible = bracket.get("admissible", True)
|
||
cost = entropy * (1.0 if admissible else 1.5)
|
||
self._cache[ck] = cost
|
||
except Exception:
|
||
cost = 1.0 if admissible else 2.0
|
||
self._cache[ck] = cost
|
||
|
||
return cost
|
||
|
||
def compute_photonic_pair_cost(self, b1: dict, b2: dict) -> float:
|
||
"""Compute interaction cost between two brackets using SLOS.
|
||
|
||
Encodes both brackets into a 2M-mode interferometer and measures
|
||
the interference cross-term via a cascade BS bridge between the
|
||
two halves. Overlapping brackets produce more cross-half
|
||
entanglement (higher cost); diverse gaps produce cleaner
|
||
separation (lower cost).
|
||
"""
|
||
if not _HAS_PERVERSE:
|
||
gap1 = abs(b1.get("gap", 32768)) / 65536.0
|
||
gap2 = abs(b2.get("gap", 32768)) / 65536.0
|
||
gap_diff = abs(gap1 - gap2)
|
||
l1, u1 = b1.get("lower", 0), b1.get("upper", 0)
|
||
l2, u2 = b2.get("lower", 0), b2.get("upper", 0)
|
||
overlap = 1.0 if l1 < u2 and l2 < u1 else 0.0
|
||
return overlap * 2.0 - gap_diff * 0.1
|
||
|
||
try:
|
||
n = self.M * 2
|
||
theta1 = self._encode_bracket_angles(b1)
|
||
theta2 = self._encode_bracket_angles(b2)
|
||
|
||
circuit = _pcvl.Circuit(n)
|
||
# Phase shifts for bracket 1 on left half (modes 0..M-1)
|
||
for i in range(self.M):
|
||
circuit.add(i, _pcvl.PS(theta1[i] if i < len(theta1) else 0.0))
|
||
|
||
# Cascade BS bridge between the two halves (consecutive ports only)
|
||
circuit.add((self.M - 1, self.M), _pcvl.BS())
|
||
|
||
# Phase shifts for bracket 2 on right half (modes M..2M-1)
|
||
for i in range(self.M):
|
||
circuit.add(self.M + i, _pcvl.PS(theta2[i] if i < len(theta2) else 0.0))
|
||
|
||
input_state = _pcvl.BasicState(
|
||
[1] + [0] * (self.M - 1) + [1] + [0] * (self.M - 1)
|
||
)
|
||
processor = _pcvl.Processor("SLOS", circuit)
|
||
processor.with_input(input_state)
|
||
sampler = _pcvl.algorithm.Sampler(processor)
|
||
res = sampler.sample_count(1000)
|
||
|
||
# Interaction cost: probability that photons from both halves
|
||
# end up on the same side (measure of entanglement)
|
||
left_same = 0.0
|
||
right_same = 0.0
|
||
for state, count in res["results"].items():
|
||
prob = count / 1000.0
|
||
left_photons = sum(state[i] for i in range(self.M))
|
||
right_photons = sum(state[i] for i in range(self.M, n))
|
||
if left_photons == 2:
|
||
left_same += prob
|
||
elif right_photons == 2:
|
||
right_same += prob
|
||
|
||
# Both photons on the same half = stronger interaction
|
||
return (left_same + right_same) * 5.0
|
||
except Exception:
|
||
return 1.0
|
||
|
||
def clear_cache(self) -> None:
|
||
self._cache.clear()
|
||
|
||
|
||
# =========================================================================
|
||
# V-F. SLOS Calibration: Fit BraidCostCoeffs to Photonic Reference Costs
|
||
# =========================================================================
|
||
|
||
def slos_calibrate_coeffs(
|
||
brackets: list[dict],
|
||
evaluator: Optional[SlosBraidEvaluator] = None,
|
||
loss_fn: str = "mse",
|
||
) -> dict:
|
||
"""Fit BraidCostCoeffs to match SLOS photonic reference costs.
|
||
|
||
For each bracket and each pair, computes the SLOS photonic cost
|
||
(reference). Then searches the BraidCostCoeffs parameter space
|
||
to find coefficients that minimize the error between the heuristic
|
||
cost (tunable_bracket_cost / tunable_crossing_penalty) and the
|
||
SLOS reference.
|
||
|
||
Args:
|
||
brackets: List of bracket dicts to calibrate on
|
||
evaluator: SlosBraidEvaluator instance (created if None)
|
||
loss_fn: "mse" (mean squared error) or "mae" (mean absolute error)
|
||
|
||
Returns:
|
||
{
|
||
"calibrated_coeffs": {BraidCostCoeffs fields},
|
||
"reference_costs": list[float], # SLOS bracket costs
|
||
"reference_pair_costs": list[float], # SLOS pair costs
|
||
"heuristic_costs": list[float], # best-fit heuristic costs
|
||
"loss": float, # final loss
|
||
"slos_available": bool, # whether Perceval was used
|
||
}
|
||
"""
|
||
if evaluator is None:
|
||
evaluator = SlosBraidEvaluator()
|
||
|
||
n = len(brackets)
|
||
if n == 0:
|
||
return {"calibrated_coeffs": {}, "error": "no brackets"}
|
||
|
||
# Compute reference costs from SLOS
|
||
ref_bracket: list[float] = []
|
||
for b in brackets:
|
||
ref_bracket.append(evaluator.compute_photonic_cost(b))
|
||
|
||
ref_pairs: list[float] = []
|
||
pair_indices: list[tuple[int, int]] = []
|
||
for i in range(n):
|
||
for j in range(i + 1, n):
|
||
ref_pairs.append(evaluator.compute_photonic_pair_cost(brackets[i], brackets[j]))
|
||
pair_indices.append((i, j))
|
||
|
||
# Grid search over BraidCostCoeffs to minimize error
|
||
best_loss = float("inf")
|
||
best_params: dict = {}
|
||
|
||
# Parameter sweep ranges — centered around defaults from braid_search.py
|
||
base_range = [0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 4.0, 8.0]
|
||
gap_range = [0.01, 0.025, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5]
|
||
overlap_range = [0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0]
|
||
|
||
for base_scale in base_range:
|
||
for gap_scale in gap_range:
|
||
for overlap_scale in overlap_range:
|
||
coeffs = BraidCostCoeffs(
|
||
base_admissible=base_scale,
|
||
base_inadmissible=base_scale * 2.0,
|
||
gap_reward_scale=gap_scale,
|
||
overlap_penalty=overlap_scale * 2.0,
|
||
)
|
||
|
||
# Compute heuristic costs with these coeffs
|
||
heur_bracket = [tunable_bracket_cost(b, coeffs) for b in brackets]
|
||
heur_pairs = [
|
||
tunable_crossing_penalty(brackets[i], brackets[j], coeffs)
|
||
for (i, j) in pair_indices
|
||
]
|
||
|
||
# Loss
|
||
loss = 0.0
|
||
for h, r in zip(heur_bracket, ref_bracket):
|
||
diff = h - r
|
||
loss += diff * diff if loss_fn == "mse" else abs(diff)
|
||
for h, r in zip(heur_pairs, ref_pairs):
|
||
diff = h - r
|
||
loss += diff * diff if loss_fn == "mse" else abs(diff)
|
||
|
||
if loss < best_loss:
|
||
best_loss = loss
|
||
best_params = asdict(coeffs)
|
||
|
||
return {
|
||
"calibrated_coeffs": best_params,
|
||
"reference_costs": ref_bracket,
|
||
"reference_pair_costs": ref_pairs,
|
||
"heuristic_costs": [
|
||
tunable_bracket_cost(b, BraidCostCoeffs(**best_params))
|
||
for b in brackets
|
||
],
|
||
"loss": best_loss,
|
||
"slos_available": _HAS_PERVERSE,
|
||
}
|
||
|
||
|
||
# =========================================================================
|
||
# V-G. Auto-Calibrate: SLOS → Fit → Solve → Compare → Recommend
|
||
#
|
||
# End-to-end pipeline: for a given set of brackets, runs the SLOS
|
||
# photonic emulator to produce physically-grounded reference costs,
|
||
# fits the best BraidCostCoeffs to match them, builds QUBOs with
|
||
# both default and calibrated coefficients, solves each with all
|
||
# available solvers, and recommends the optimal configuration.
|
||
# =========================================================================
|
||
|
||
def auto_calibrate(
|
||
brackets: list[dict],
|
||
evaluator: Optional[SlosBraidEvaluator] = None,
|
||
time_per_solver: float = 1.0,
|
||
seed: int = 42,
|
||
) -> dict:
|
||
"""Full auto-calibration pipeline for braid QUBO defaults.
|
||
|
||
Steps:
|
||
1. Compute SLOS photonic reference costs for each bracket and pair.
|
||
2. Grid-search BraidCostCoeffs to minimize heuristic-vs-SLOS error.
|
||
3. Build QUBOs with default (DEFAULT_COEFFS) and calibrated coeffs.
|
||
4. Solve both QUBOs with SA, Lévy, HiGHS, and QAOA describe.
|
||
5. Report the optimal coefficient set and energy improvement.
|
||
|
||
Args:
|
||
brackets: List of bracket dicts.
|
||
evaluator: SlosBraidEvaluator (created fresh if None).
|
||
time_per_solver: Seconds per solver trial.
|
||
seed: RNG seed.
|
||
|
||
Returns:
|
||
{
|
||
"slos_calibration": {...}, # raw calibration output
|
||
"coefficients": {
|
||
"default": {...}, # original DEFAULT_COEFFS
|
||
"calibrated": {...}, # SLOS-fitted coeffs
|
||
"recommended": {...}, # whichever gave lowest energy
|
||
},
|
||
"solver_results": {
|
||
"default": {method: result}, # QUBO with default coeffs
|
||
"calibrated": {method: result}, # QUBO with SLOS coeffs
|
||
},
|
||
"comparison": {
|
||
"default_best": {method, energy},
|
||
"calibrated_best": {method, energy},
|
||
"improvement": float, # energy reduction (negative = better)
|
||
},
|
||
"recommendation": str, # human-readable summary
|
||
"n": int,
|
||
}
|
||
"""
|
||
if evaluator is None:
|
||
evaluator = SlosBraidEvaluator()
|
||
if not brackets:
|
||
return {"error": "no brackets provided"}
|
||
|
||
n = len(brackets)
|
||
|
||
# Step 1-2: SLOS calibration
|
||
calibration = slos_calibrate_coeffs(brackets, evaluator=evaluator)
|
||
cal_params = calibration["calibrated_coeffs"]
|
||
cal_coeffs = BraidCostCoeffs(**cal_params)
|
||
|
||
# Step 3: Build QUBOs with default and calibrated coeffs
|
||
Q_default = build_tunable_qubo_matrix(brackets, DEFAULT_COEFFS)
|
||
Q_cal = build_tunable_qubo_matrix(brackets, cal_coeffs)
|
||
qubo_default = QUBO(n=n, matrix=Q_default)
|
||
qubo_cal = QUBO(n=n, matrix=Q_cal)
|
||
|
||
# Step 4: Solve with all methods
|
||
methods = ["sa", "levy"]
|
||
default_results: dict[str, dict] = {}
|
||
cal_results: dict[str, dict] = {}
|
||
|
||
for method in methods:
|
||
default_results[method] = stochastic_abuse_qubo(
|
||
qubo_default, method=method, time_limit=time_per_solver, seed=seed,
|
||
)
|
||
cal_results[method] = stochastic_abuse_qubo(
|
||
qubo_cal, method=method, time_limit=time_per_solver, seed=seed,
|
||
)
|
||
|
||
# Also try HiGHS
|
||
for method in ["highs"]:
|
||
try:
|
||
default_results[method] = stochastic_abuse_qubo(
|
||
qubo_default, method=method, time_limit=time_per_solver, seed=seed,
|
||
)
|
||
except Exception:
|
||
default_results[method] = {"energy": float("inf"), "solution": [0] * n}
|
||
try:
|
||
cal_results[method] = stochastic_abuse_qubo(
|
||
qubo_cal, method=method, time_limit=time_per_solver, seed=seed,
|
||
)
|
||
except Exception:
|
||
cal_results[method] = {"energy": float("inf"), "solution": [0] * n}
|
||
|
||
# Step 5: Find best per coefficient set
|
||
default_best = min(
|
||
default_results.items(), key=lambda kv: kv[1]["energy"],
|
||
)
|
||
cal_best = min(
|
||
cal_results.items(), key=lambda kv: kv[1]["energy"],
|
||
)
|
||
|
||
improvement = cal_best[1]["energy"] - default_best[1]["energy"]
|
||
|
||
# Decide which coeffs to recommend
|
||
if cal_best[1]["energy"] <= default_best[1]["energy"]:
|
||
recommended = cal_params
|
||
rec_note = "SLOS-calibrated coefficients (lower energy)"
|
||
else:
|
||
recommended = asdict(DEFAULT_COEFFS)
|
||
rec_note = "default coefficients (SLOS calibration did not improve)"
|
||
|
||
return {
|
||
"slos_calibration": {
|
||
"loss": calibration["loss"],
|
||
"reference_costs": calibration["reference_costs"],
|
||
"slos_available": calibration["slos_available"],
|
||
},
|
||
"coefficients": {
|
||
"default": asdict(DEFAULT_COEFFS),
|
||
"calibrated": cal_params,
|
||
"recommended": recommended,
|
||
},
|
||
"solver_results": {
|
||
"default": default_results,
|
||
"calibrated": cal_results,
|
||
},
|
||
"comparison": {
|
||
"default_best": {"method": default_best[0], "energy": default_best[1]["energy"]},
|
||
"calibrated_best": {"method": cal_best[0], "energy": cal_best[1]["energy"]},
|
||
"improvement": improvement,
|
||
"improvement_pct": (
|
||
improvement / max(abs(default_best[1]["energy"]), 1e-10) * 100
|
||
if abs(default_best[1]["energy"]) > 1e-10 else 0.0
|
||
),
|
||
},
|
||
"recommendation": (
|
||
f"Set braid_search defaults to: base_admissible={recommended.get('base_admissible')}, "
|
||
f"gap_reward_scale={recommended.get('gap_reward_scale')}, "
|
||
f"overlap_penalty={recommended.get('overlap_penalty')}. "
|
||
f"{rec_note}."
|
||
),
|
||
"n": n,
|
||
}
|
||
|
||
|
||
# =========================================================================
|
||
# III-D. FinslerMetric → QUBO (TransportQUBOBridge.lean)
|
||
# Source: TransportTheory.lean — RandersMetric { alpha, beta }
|
||
# alpha = AlphaComponent { dimension, mass_field }
|
||
# beta = BetaComponent { dimension, wind_field }
|
||
#
|
||
# Q_ij = F(p_i, v_j - v_i) = α(v_j - v_i) + β(v_j - v_i)
|
||
# where α(v) = Σ_i mass[i] * |v_i| (symmetric: α(-v) = α(v))
|
||
# β(v) = Σ_i wind[i] * v_i (antisymmetric: β(-v) = -β(v))
|
||
#
|
||
# The matrix is NOT symmetric: Q_ij ≠ Q_ji when β ≠ 0.
|
||
# Each row/col = one routing direction.
|
||
# =========================================================================
|
||
|
||
@dataclass
|
||
class FinslerMetric:
|
||
"""Discrete Finsler-Randers metric (mirrors Lean RandersMetric).
|
||
|
||
α component: symmetric base cost (mass field)
|
||
β component: asymmetric drift (wind field)
|
||
|
||
F(p,v) = α(p,v) + β(p,v)
|
||
"""
|
||
alpha_mass: list[float] # mass_field: base resistance per dimension
|
||
beta_wind: list[float] # wind_field: drift 1-form per dimension
|
||
dimension: int # state space dimension
|
||
|
||
def __post_init__(self) -> None:
|
||
if len(self.alpha_mass) != self.dimension:
|
||
raise ValueError(f"alpha mass dim {len(self.alpha_mass)} != {self.dimension}")
|
||
if len(self.beta_wind) != self.dimension:
|
||
raise ValueError(f"beta wind dim {len(self.beta_wind)} != {self.dimension}")
|
||
|
||
def alpha_cost(self, v: list[float]) -> float:
|
||
"""α(v) = Σ_i mass[i] * |v_i| — symmetric."""
|
||
return sum(self.alpha_mass[i] * abs(v[i]) for i in range(self.dimension))
|
||
|
||
def beta_cost(self, v: list[float]) -> float:
|
||
"""β(v) = Σ_i wind[i] * v_i — antisymmetric (β(-v) = -β(v))."""
|
||
return sum(self.beta_wind[i] * v[i] for i in range(self.dimension))
|
||
|
||
def finsler_cost(self, v: list[float]) -> float:
|
||
"""F(p,v) = α(v) + β(v) — direction-dependent Finsler norm."""
|
||
return self.alpha_cost(v) + self.beta_cost(v)
|
||
|
||
def crossing_cost(self, v_i: list[float], v_j: list[float]) -> float:
|
||
"""Cost of transitioning from direction i to direction j.
|
||
|
||
Q_ij = F(p, v_j - v_i) = α(v_j - v_i) + β(v_j - v_i)
|
||
|
||
NOTE: Q_ij ≠ Q_ji because β(v_j - v_i) = -β(v_i - v_j).
|
||
"""
|
||
diff = [v_j[k] - v_i[k] for k in range(self.dimension)]
|
||
return self.finsler_cost(diff)
|
||
|
||
|
||
def finsler_metric_to_qubo(
|
||
metric: FinslerMetric,
|
||
directions: list[list[float]],
|
||
normalize: bool = True,
|
||
) -> QUBO:
|
||
"""Convert a FinslerMetric + direction set into a QUBO formulation.
|
||
|
||
This mirrors the Lean `randersMetricToQUBO` in TransportQUBOBridge.lean.
|
||
|
||
Args:
|
||
metric: FinslerMetric (α mass + β wind).
|
||
directions: List of n direction vectors, each length metric.dimension.
|
||
normalize: If True, normalize costs so the matrix entries are in [-1, 1].
|
||
|
||
Returns:
|
||
QUBO with n variables. Q_ij = crossing cost from i→j, Q_ii = 0.
|
||
"""
|
||
n = len(directions)
|
||
Q: dict[tuple[int, int], float] = {}
|
||
max_abs = 0.0
|
||
|
||
# Build raw QUBO matrix
|
||
raw: list[list[float]] = [[0.0] * n for _ in range(n)]
|
||
for i in range(n):
|
||
for j in range(n):
|
||
if i == j:
|
||
raw[i][j] = 0.0
|
||
else:
|
||
cost = metric.crossing_cost(directions[i], directions[j])
|
||
raw[i][j] = cost
|
||
max_abs = max(max_abs, abs(cost))
|
||
|
||
# Normalize if requested
|
||
scale = max_abs if normalize and max_abs > 0 else 1.0
|
||
for i in range(n):
|
||
for j in range(i + 1, n):
|
||
val = (raw[i][j] + raw[j][i]) / scale
|
||
if val != 0.0:
|
||
Q[(i, j)] = val
|
||
|
||
return QUBO(n=n, matrix=Q)
|
||
|
||
|
||
def is_anisotropic(
|
||
qubo: QUBO,
|
||
raw_matrix: Optional[list[list[float]]] = None,
|
||
tol: float = 1e-9,
|
||
) -> bool:
|
||
"""Check if a QUBO matrix is anisotropic (Q_ij ≠ Q_ji for any i≠j).
|
||
|
||
After the fix in finsler_metric_to_qubo (storing Q_ij + Q_ji),
|
||
the QUBO matrix no longer preserves individual directed entries.
|
||
Pass raw_matrix to compare the original asymmetric entries, or
|
||
pass a FinslerMetric to check the wind field directly.
|
||
|
||
Mirrors Lean `isAnisotropic` in TransportQUBOBridge.lean.
|
||
"""
|
||
if raw_matrix is not None:
|
||
n = len(raw_matrix)
|
||
for i in range(n):
|
||
for j in range(n):
|
||
if i != j and abs(raw_matrix[i][j] - raw_matrix[j][i]) > tol:
|
||
return True
|
||
return False
|
||
# Fallback: check stored entries (works for preserved asymmetric storage)
|
||
for (i, j), qij in qubo.matrix.items():
|
||
qji = qubo.matrix.get((j, i), 0.0)
|
||
if abs(qij - qji) > tol:
|
||
return True
|
||
return False
|
||
|
||
|
||
def geodesic_assignment(
|
||
metric: FinslerMetric,
|
||
directions: list[list[float]],
|
||
tol: float = 1e-9,
|
||
) -> list[bool]:
|
||
"""Extract the Finsler geodesic directions as a Boolean assignment.
|
||
|
||
A direction is 'on the geodesic' if its Finsler cost is within tol of
|
||
the global minimum.
|
||
|
||
Mirrors Lean `geodesicAssignment` in TransportQUBOBridge.lean.
|
||
"""
|
||
costs = [metric.finsler_cost(v) for v in directions]
|
||
min_cost = min(costs)
|
||
return [abs(c - min_cost) <= tol for c in costs]
|
||
|
||
|
||
# =========================================================================
|
||
# IV. QUBO → Ising
|
||
# =========================================================================
|
||
|
||
def qubo_to_ising(qubo: QUBO) -> Ising:
|
||
"""Convert a QUBO to an Ising Hamiltonian.
|
||
|
||
For binary variable x ∈ {0,1}, spin s = 2x - 1 ∈ {+1, -1}:
|
||
x = (1 + s) / 2
|
||
x_i x_j = (1 + s_i + s_j + s_i s_j) / 4
|
||
|
||
The QUBO cost: E = Σ_i a_i x_i + Σ_{i<j} b_ij x_i x_j
|
||
where a_i = Q[i,i], b_ij = Q[i,j] + Q[j,i]
|
||
"""
|
||
n = qubo.n
|
||
a: dict[int, float] = {}
|
||
b: dict[tuple[int, int], float] = {}
|
||
|
||
for (i, j), qij in qubo.matrix.items():
|
||
if i == j:
|
||
a[i] = a.get(i, 0.0) + qij
|
||
else:
|
||
key = (min(i, j), max(i, j))
|
||
b[key] = b.get(key, 0.0) + qij
|
||
|
||
offset = qubo.offset
|
||
h = [0.0] * n
|
||
J: dict[tuple[int, int], float] = {}
|
||
|
||
for i in range(n):
|
||
offset += 0.5 * a.get(i, 0.0)
|
||
for (i, j), bij in b.items():
|
||
offset += 0.25 * bij
|
||
|
||
for i in range(n):
|
||
hi = 0.5 * a.get(i, 0.0)
|
||
for (j1, j2), bij in b.items():
|
||
if j1 == i:
|
||
hi += 0.25 * bij
|
||
if j2 == i:
|
||
hi += 0.25 * bij
|
||
h[i] = hi
|
||
|
||
for (i, j), bij in b.items():
|
||
J[(i, j)] = 0.25 * bij
|
||
|
||
return Ising(n=n, h=h, J=J, offset=offset)
|
||
|
||
|
||
# =========================================================================
|
||
# VII. Core: Ising → Pauli strings
|
||
# =========================================================================
|
||
|
||
def ising_to_pauli(ising: Ising) -> PauliSum:
|
||
"""Convert an Ising Hamiltonian to Pauli strings.
|
||
|
||
Mapping:
|
||
s_i → Z_i
|
||
s_i s_j → Z_i Z_j
|
||
offset → I (constant)
|
||
"""
|
||
terms: list[tuple[str, float]] = []
|
||
|
||
for i in range(ising.n):
|
||
if abs(ising.h[i]) > 1e-15:
|
||
ps = ["I"] * ising.n
|
||
ps[i] = "Z"
|
||
terms.append(("".join(ps), ising.h[i]))
|
||
|
||
for (i, j), Jij in ising.J.items():
|
||
if abs(Jij) > 1e-15:
|
||
ps = ["I"] * ising.n
|
||
ps[i] = "Z"
|
||
ps[j] = "Z"
|
||
terms.append(("".join(ps), Jij))
|
||
|
||
return PauliSum(n=ising.n, terms=terms, offset=ising.offset)
|
||
|
||
|
||
# =========================================================================
|
||
# VIII. Forward: Pauli → Cirq Circuit (QAOA layers)
|
||
# =========================================================================
|
||
|
||
def pauli_to_cirq(
|
||
pauli: PauliSum,
|
||
p_layers: int = 1,
|
||
gamma: Optional[list[float]] = None,
|
||
beta: Optional[list[float]] = None,
|
||
measure: bool = False,
|
||
) -> Any:
|
||
"""Generate a QAOA circuit from a PauliSum Hamiltonian.
|
||
|
||
Returns cirq.Circuit if cirq is available, else a dict description.
|
||
|
||
Circuit structure:
|
||
|+⟩^⊗n → [e^{-iγ_k H_C} e^{-iβ_k H_M}]_{k=1..p}
|
||
"""
|
||
if not _HAS_CIRQ:
|
||
return _describe_qaoa_circuit(pauli, p_layers)
|
||
|
||
n = pauli.n
|
||
qubits = cirq.LineQubit.range(n)
|
||
|
||
if gamma is None:
|
||
gamma = [1.0] * p_layers
|
||
if beta is None:
|
||
beta = [1.0] * p_layers
|
||
|
||
circuit = cirq.Circuit()
|
||
circuit.append(cirq.H.on_each(*qubits))
|
||
|
||
for layer in range(p_layers):
|
||
g = gamma[layer] if layer < len(gamma) else gamma[-1]
|
||
b = beta[layer] if layer < len(beta) else beta[-1]
|
||
|
||
for ps_str, coeff in pauli.terms:
|
||
angle = 2.0 * g * coeff
|
||
if abs(angle) < 1e-15:
|
||
continue
|
||
z_positions = [i for i, c in enumerate(ps_str) if c == "Z"]
|
||
if len(z_positions) == 1:
|
||
circuit.append(cirq.rz(angle)(qubits[z_positions[0]]))
|
||
elif len(z_positions) == 2:
|
||
i, j = z_positions
|
||
circuit.append(cirq.CZ(qubits[i], qubits[j]) ** (angle / math.pi))
|
||
circuit.append(cirq.rz(angle)(qubits[i]))
|
||
circuit.append(cirq.rz(angle)(qubits[j]))
|
||
elif len(z_positions) > 2:
|
||
for idx in z_positions:
|
||
circuit.append(cirq.rz(angle / len(z_positions))(qubits[idx]))
|
||
|
||
circuit.append(cirq.rx(2.0 * b).on_each(*qubits))
|
||
|
||
if measure:
|
||
circuit.append(cirq.measure(*qubits, key="result"))
|
||
|
||
return circuit
|
||
|
||
|
||
def _describe_qaoa_circuit(pauli: PauliSum, p_layers: int = 1) -> dict:
|
||
"""Return a JSON-like circuit description when no simulator is available."""
|
||
return {
|
||
"type": "qaoa_circuit_description",
|
||
"n_qubits": pauli.n,
|
||
"p_layers": p_layers,
|
||
"cost_hamiltonian": [{"pauli": t[0], "coeff": t[1]} for t in pauli.terms],
|
||
"offset": pauli.offset,
|
||
"mixer": "X" * pauli.n,
|
||
"note": "Install cirq to generate runnable circuits",
|
||
}
|
||
|
||
|
||
# =========================================================================
|
||
# IX. Forward: Pauli → Qiskit Circuit (QAOA layers)
|
||
# =========================================================================
|
||
|
||
def pauli_to_qiskit(
|
||
pauli: PauliSum,
|
||
p_layers: int = 1,
|
||
gamma: Optional[list[float]] = None,
|
||
beta: Optional[list[float]] = None,
|
||
) -> Any:
|
||
"""Generate a QAOA circuit using Qiskit.
|
||
|
||
Returns QuantumCircuit if qiskit is available, else a dict description.
|
||
"""
|
||
if not _HAS_QISKIT:
|
||
return _describe_qaoa_circuit(pauli, p_layers)
|
||
|
||
n = pauli.n
|
||
qreg = QuantumRegister(n, "q")
|
||
creg = ClassicalRegister(n, "c")
|
||
circuit = QuantumCircuit(qreg, creg)
|
||
|
||
if gamma is None:
|
||
gamma = [1.0] * p_layers
|
||
if beta is None:
|
||
beta = [1.0] * p_layers
|
||
|
||
circuit.h(qreg)
|
||
|
||
for layer in range(p_layers):
|
||
g = gamma[layer] if layer < len(gamma) else gamma[-1]
|
||
b = beta[layer] if layer < len(beta) else beta[-1]
|
||
|
||
for ps_str, coeff in pauli.terms:
|
||
angle = 2.0 * g * coeff
|
||
if abs(angle) < 1e-15:
|
||
continue
|
||
z_positions = [i for i, c in enumerate(ps_str) if c == "Z"]
|
||
if len(z_positions) == 1:
|
||
circuit.rz(angle, qreg[z_positions[0]])
|
||
elif len(z_positions) == 2:
|
||
i, j = z_positions
|
||
circuit.cx(qreg[i], qreg[j])
|
||
circuit.rz(angle, qreg[j])
|
||
circuit.cx(qreg[i], qreg[j])
|
||
elif len(z_positions) > 2:
|
||
for idx in z_positions:
|
||
circuit.rz(angle / len(z_positions), qreg[idx])
|
||
|
||
circuit.rx(2.0 * b, qreg)
|
||
|
||
circuit.measure(qreg, creg)
|
||
return circuit
|
||
|
||
|
||
# =========================================================================
|
||
# X. Backward: Measurements → Solution
|
||
# =========================================================================
|
||
|
||
def measurements_to_solution(
|
||
counts: dict[str, int],
|
||
n: int,
|
||
) -> list[int]:
|
||
"""Extract the most probable bitstring from measurement counts.
|
||
|
||
Returns binary list [x_0, ..., x_{n-1}] of the most probable outcome.
|
||
"""
|
||
if not counts:
|
||
return [0] * n
|
||
best_bitstring = max(counts, key=counts.get)
|
||
if len(best_bitstring) < n:
|
||
best_bitstring = best_bitstring.zfill(n)
|
||
return [int(b) for b in best_bitstring[-n:]]
|
||
|
||
|
||
def measurements_to_ising_solution(
|
||
counts: dict[str, int], n: int,
|
||
) -> list[int]:
|
||
"""Extract the most probable spin configuration from measurements."""
|
||
x = measurements_to_solution(counts, n)
|
||
return [2 * xi - 1 for xi in x]
|
||
|
||
|
||
# =========================================================================
|
||
# XI. Backward: Solution → BraidReceipt Update
|
||
# =========================================================================
|
||
|
||
def solution_to_braid_receipt(solution: list[int], receipt: dict) -> dict:
|
||
"""Update a BraidReceipt dict with a QAOA solution."""
|
||
updated = json.loads(json.dumps(receipt))
|
||
braid = updated.get("braid", updated)
|
||
strands = braid.get("strands", [])
|
||
if not strands:
|
||
strands = _receipt_strands_from_bracket(braid)
|
||
braid["strands"] = strands
|
||
|
||
n = min(len(solution), len(strands))
|
||
all_admissible = True
|
||
for i in range(n):
|
||
s = strands[i]
|
||
if bool(solution[i]):
|
||
s["kappa"] = max(0, s.get("kappa", 0) // 2)
|
||
s["converged"] = True
|
||
s["residual"] = 0
|
||
else:
|
||
s["kappa"] = min(65535, s.get("kappa", 0) + 8192)
|
||
s["converged"] = False
|
||
s["residual"] = s.get("kappa", 0)
|
||
all_admissible = False
|
||
|
||
braid["bracket"] = braid.get("bracket", {})
|
||
braid["bracket"]["admissible"] = all_admissible
|
||
updated["scar_absent"] = all_admissible
|
||
updated["step_count"] = updated.get("step_count", 0) + 1
|
||
|
||
residuals = updated.get("residuals", [])
|
||
avg_residual = sum(abs(s.get("kappa", 0)) for s in strands[:n]) // max(n, 1)
|
||
residuals.insert(0, avg_residual)
|
||
updated["residuals"] = residuals[:32]
|
||
return updated
|
||
|
||
|
||
# =========================================================================
|
||
# XII. Backward: Solution → Scar/β₀ Update
|
||
# =========================================================================
|
||
|
||
def solution_to_scar_field(solution: list[int], N: int) -> list[int]:
|
||
"""Map a QAOA solution bitstring back to a scar field."""
|
||
mask = [0] * N
|
||
for i in range(N):
|
||
if solution[i] == 1:
|
||
mask[i] = 1
|
||
radius = max(1, N // 32)
|
||
expanded = [0] * N
|
||
for i in range(N):
|
||
if mask[i]:
|
||
for d in range(-radius, radius + 1):
|
||
expanded[(i + d) % N] = 1
|
||
return expanded
|
||
|
||
|
||
def compute_beta0_from_solution(solution: list[int]) -> int:
|
||
"""Compute β₀ (rising edge count) from a cyclic binary solution."""
|
||
n = len(solution)
|
||
if n == 0:
|
||
return 0
|
||
if all(s == 1 for s in solution):
|
||
return 1
|
||
if all(s == 0 for s in solution):
|
||
return 0
|
||
rising = 0
|
||
for i in range(n):
|
||
if solution[i] == 1 and solution[(i - 1) % n] == 0:
|
||
rising += 1
|
||
return rising
|
||
|
||
|
||
# =========================================================================
|
||
# XIII. High-Level: Corpus250 Row → QUBO
|
||
# FixtureRow has: shape, rrcKind, operatorTokens, weakAxesCnt,
|
||
# pistProxyLabel, pistExactLabel (no domain_embedding field).
|
||
# =========================================================================
|
||
|
||
def corpus278_row_to_qubo(row: dict) -> QUBO:
|
||
"""Convert a Corpus250 equation row (FixtureRow) to a QUBO problem.
|
||
|
||
Uses the actual FixtureRow fields:
|
||
- shape: RRCShape enum guiding the QUBO structure
|
||
- operatorTokens: domain operators for variable count/diagonal
|
||
- weakAxesCnt: number of weak axes for off-diagonal coupling
|
||
- pistProxyLabel / pistExactLabel: PIST labels for bias
|
||
"""
|
||
shape = row.get("shape", "")
|
||
operator_tokens = row.get("operatorTokens", [])
|
||
weak_axes_cnt = row.get("weakAxesCnt", 0)
|
||
|
||
n_vars = max(8, len(operator_tokens) + int(weak_axes_cnt))
|
||
|
||
shape_bias = {
|
||
"cognitiveLoadField": 1.0,
|
||
"signalShapedRouteCompiler": 2.0,
|
||
"projectableGeometryTopology": 3.0,
|
||
"cadForceProbeReceipt": 4.0,
|
||
"logogramProjection": 5.0,
|
||
"holdForUnlawfulOrUnderspecifiedShape": 10.0,
|
||
}.get(shape, 2.0)
|
||
|
||
Q: dict[tuple[int, int], float] = {}
|
||
|
||
# Diagonal: one per operator token with shape bias
|
||
for i in range(min(n_vars, len(operator_tokens))):
|
||
Q[(i, i)] = shape_bias * (1.0 + (i % (int(weak_axes_cnt) + 1)))
|
||
|
||
# Fill remaining diagonal
|
||
for i in range(len(operator_tokens), n_vars):
|
||
Q[(i, i)] = shape_bias * 1.5
|
||
|
||
# Weak axes coupling
|
||
wac = int(weak_axes_cnt)
|
||
if wac > 0:
|
||
for i in range(min(wac, n_vars)):
|
||
for j in range(i + 1, min(wac, n_vars)):
|
||
Q[(i, j)] = Q.get((i, j), 0.0) - 0.5
|
||
|
||
# PIST label bias
|
||
pist_proxy = row.get("pistProxyLabel")
|
||
if pist_proxy and pist_proxy != "none":
|
||
for i in range(min(4, n_vars)):
|
||
Q[(i, i)] = Q.get((i, i), 0.0) - 0.3
|
||
|
||
return QUBO(n=n_vars, matrix=Q)
|
||
|
||
|
||
# =========================================================================
|
||
# XIV. Stochastic Abuse: Classical Stochastic Solvers
|
||
# Wraps existing qubo_highs.py (SA / HiGHS) so results are directly
|
||
# comparable with qaoa_solve_qubo() output.
|
||
#
|
||
# "Abuse" = using the classical stochastic pipeline as a baseline to
|
||
# benchmark QAOA against, on the same QUBO — borrowing the existing
|
||
# HiGHS MIP and simulated annealing infrastructure.
|
||
# =========================================================================
|
||
|
||
def stochastic_abuse_qubo(
|
||
qubo: QUBO,
|
||
method: str = "sa",
|
||
time_limit: float = 5.0,
|
||
seed: int = 42,
|
||
) -> dict:
|
||
"""Solve a QUBO using the existing classical stochastic pipeline.
|
||
|
||
Wraps qubo_highs.solve_qubo_highs (HiGHS MIP) and
|
||
braid_search.qubo_optimize (SA) for direct comparison with QAOA.
|
||
|
||
Args:
|
||
qubo: QUBO problem
|
||
method: "sa" (simulated annealing), "highs" (HiGHS MIP), or
|
||
"levy" (Lévy flight sampling)
|
||
time_limit: Max seconds for solver
|
||
seed: RNG seed
|
||
|
||
Returns:
|
||
dict with same schema as qaoa_solve_qubo() for easy comparison:
|
||
- solution, energy, method, runtime, solver_detail
|
||
"""
|
||
t0 = time.time()
|
||
n = qubo.n
|
||
|
||
# Build QUBO dict in the format qubo_highs expects
|
||
Q_dict: dict[tuple[int, int], float] = dict(qubo.matrix)
|
||
|
||
if method == "highs":
|
||
try:
|
||
_sys.path.insert(0, str(Path(__file__).resolve().parent))
|
||
from qubo_highs import solve_qubo_highs as highs_solve
|
||
result = highs_solve(Q_dict, n, time_limit=time_limit)
|
||
if isinstance(result, dict):
|
||
solution = result.get("x", [0] * n)
|
||
else:
|
||
solution, _ = result
|
||
except Exception as exc:
|
||
solution = _sa_solve(Q_dict, n, seed, time_limit)
|
||
|
||
elif method == "levy":
|
||
solution = _levy_fight(Q_dict, n, seed, time_limit)
|
||
|
||
else:
|
||
solution = _sa_solve(Q_dict, n, seed, time_limit)
|
||
|
||
runtime = time.time() - t0
|
||
energy = qubo.energy(solution)
|
||
|
||
return {
|
||
"solution": solution,
|
||
"energy": energy,
|
||
"method": method,
|
||
"runtime_s": round(runtime, 4),
|
||
"n": n,
|
||
}
|
||
|
||
|
||
def _sa_solve(
|
||
Q: dict[tuple[int, int], float],
|
||
n: int,
|
||
seed: int = 42,
|
||
time_limit: float = 5.0,
|
||
) -> list[int]:
|
||
"""Simulated annealing QUBO solver (standalone fallback)."""
|
||
import random as _r
|
||
_r.seed(seed)
|
||
t0 = time.time()
|
||
|
||
x = [_r.randint(0, 1) for _ in range(n)]
|
||
|
||
def _energy(xv):
|
||
e = 0.0
|
||
for (i, j), qij in Q.items():
|
||
e += qij * xv[i] * xv[j]
|
||
return e
|
||
|
||
current_energy = _energy(x)
|
||
best_x = list(x)
|
||
best_energy = current_energy
|
||
temp = 10.0
|
||
iterations = 0
|
||
n_vals = list(range(n))
|
||
|
||
while time.time() - t0 < time_limit:
|
||
i = _r.choice(n_vals)
|
||
x[i] = 1 - x[i]
|
||
new_energy = _energy(x)
|
||
delta = new_energy - current_energy
|
||
|
||
if delta < 0 or _r.random() < math.exp(-delta / max(temp, 1e-10)):
|
||
current_energy = new_energy
|
||
if current_energy < best_energy:
|
||
best_x = list(x)
|
||
best_energy = current_energy
|
||
else:
|
||
x[i] = 1 - x[i]
|
||
|
||
temp *= 0.9995
|
||
iterations += 1
|
||
|
||
return best_x
|
||
|
||
|
||
def _levy_fight(
|
||
Q: dict[tuple[int, int], float],
|
||
n: int,
|
||
seed: int = 42,
|
||
time_limit: float = 5.0,
|
||
) -> list[int]:
|
||
"""Lévy flight sampling over binary strings.
|
||
|
||
Step sizes follow a power-law distribution (heavy-tailed),
|
||
mimicking the LévyFlight structure from EntropyMeasures.lean.
|
||
"""
|
||
import random as _r
|
||
_r.seed(seed)
|
||
t0 = time.time()
|
||
|
||
x = [_r.randint(0, 1) for _ in range(n)]
|
||
|
||
def _energy(xv):
|
||
e = 0.0
|
||
for (i, j), qij in Q.items():
|
||
e += qij * xv[i] * xv[j]
|
||
return e
|
||
|
||
def _levy_step() -> int:
|
||
u = _r.random()
|
||
return max(1, int(n * u ** (-1.0 / 1.5)) % n)
|
||
|
||
best_x = list(x)
|
||
best_energy = _energy(x)
|
||
|
||
while time.time() - t0 < time_limit:
|
||
step = _levy_step()
|
||
indices = _r.sample(range(n), min(step, n))
|
||
for i in indices:
|
||
x[i] = 1 - x[i]
|
||
e = _energy(x)
|
||
if e < best_energy:
|
||
best_x = list(x)
|
||
best_energy = e
|
||
else:
|
||
for i in indices:
|
||
x[i] = 1 - x[i]
|
||
|
||
return best_x
|
||
|
||
|
||
# =========================================================================
|
||
# XV. Comparison: QAOA vs Stochastic
|
||
# =========================================================================
|
||
|
||
def qaoa_vs_stochastic_comparison(
|
||
qubo: QUBO,
|
||
p_layers: int = 1,
|
||
shots: int = 1000,
|
||
qaoa_backend: str = "cirq",
|
||
stochastic_methods: Optional[list[str]] = None,
|
||
time_limit: float = 5.0,
|
||
seed: int = 42,
|
||
) -> dict:
|
||
"""Run QAOA and classical stochastic solvers on the same QUBO.
|
||
|
||
Returns a side-by-side comparison so results can be contrasted.
|
||
|
||
Args:
|
||
qubo: QUBO problem
|
||
p_layers: QAOA depth
|
||
shots: QAOA measurement shots
|
||
qaoa_backend: QAOA backend ("cirq", "qiskit", "describe")
|
||
stochastic_methods: List of classical methods to compare.
|
||
Default: ["sa", "levy"].
|
||
time_limit: Time limit per classical solver (seconds)
|
||
seed: RNG seed
|
||
|
||
Returns:
|
||
dict with:
|
||
- qubo_summary: n, terms, offset
|
||
- qaoa: result from qaoa_solve_qubo
|
||
- stochastic: {method: result} for each method
|
||
- winner: solver with lowest energy
|
||
"""
|
||
if stochastic_methods is None:
|
||
stochastic_methods = ["sa", "levy"]
|
||
|
||
qaoa_result = qaoa_solve_qubo(
|
||
qubo, p_layers=p_layers, shots=shots, backend=qaoa_backend,
|
||
)
|
||
|
||
stochastic_results: dict[str, dict] = {}
|
||
for method in stochastic_methods:
|
||
stochastic_results[method] = stochastic_abuse_qubo(
|
||
qubo, method=method, time_limit=time_limit, seed=seed,
|
||
)
|
||
|
||
# Determine winner
|
||
candidates: list[tuple[str, float]] = [
|
||
("qaoa", qaoa_result["energy"]),
|
||
]
|
||
for method, res in stochastic_results.items():
|
||
candidates.append((method, res["energy"]))
|
||
|
||
candidates.sort(key=lambda p: p[1])
|
||
winner_name, winner_energy = candidates[0]
|
||
|
||
return {
|
||
"qubo_summary": {
|
||
"n": qubo.n,
|
||
"terms": len(qubo.matrix),
|
||
"offset": qubo.offset,
|
||
},
|
||
"qaoa": {
|
||
"solution": qaoa_result["solution"],
|
||
"energy": qaoa_result["energy"],
|
||
"p_layers": p_layers,
|
||
"shots": shots,
|
||
"backend": qaoa_backend,
|
||
},
|
||
"stochastic": stochastic_results,
|
||
"winner": {
|
||
"solver": winner_name,
|
||
"energy": winner_energy,
|
||
},
|
||
}
|
||
|
||
|
||
# =========================================================================
|
||
# XVI. High-Level: End-to-End QAOA Solve
|
||
# =========================================================================
|
||
|
||
def qaoa_solve_qubo(
|
||
qubo: QUBO,
|
||
p_layers: int = 1,
|
||
gamma: Optional[list[float]] = None,
|
||
beta: Optional[list[float]] = None,
|
||
shots: int = 1000,
|
||
backend: str = "cirq",
|
||
) -> dict:
|
||
"""Solve a QUBO using QAOA.
|
||
|
||
Args:
|
||
qubo: QUBO problem
|
||
p_layers: QAOA layers
|
||
gamma: cost angles
|
||
beta: mixer angles
|
||
shots: measurement shots
|
||
backend: "cirq", "qiskit", or "describe"
|
||
|
||
Returns:
|
||
dict with solution, energy, counts, and circuit description
|
||
"""
|
||
ising = qubo_to_ising(qubo)
|
||
pauli = ising_to_pauli(ising)
|
||
|
||
result: dict[str, Any] = {
|
||
"n": qubo.n,
|
||
"p_layers": p_layers,
|
||
"qubo_energy_offset": qubo.offset,
|
||
"ising_offset": ising.offset,
|
||
"pauli_terms": [(t[0], t[1]) for t in pauli.terms],
|
||
}
|
||
|
||
if backend == "describe":
|
||
result["circuit"] = _describe_qaoa_circuit(pauli, p_layers)
|
||
result["solution"] = [0] * qubo.n
|
||
result["energy"] = qubo.energy([0] * qubo.n)
|
||
result["note"] = "describe mode — no simulation"
|
||
return result
|
||
|
||
if backend == "cirq" and _HAS_CIRQ:
|
||
circuit = pauli_to_cirq(pauli, p_layers, gamma, beta, measure=True)
|
||
simulator = cirq.Simulator()
|
||
samples = simulator.run(circuit, repetitions=shots)
|
||
counts = samples.histogram(key="result")
|
||
str_counts = _cirq_counts_to_str(counts, qubo.n)
|
||
solution = measurements_to_solution(str_counts, qubo.n)
|
||
energy = qubo.energy(solution)
|
||
result["circuit"] = str(circuit)
|
||
result["counts"] = str_counts
|
||
result["solution"] = solution
|
||
result["energy"] = energy
|
||
return result
|
||
|
||
if backend == "qiskit" and _HAS_QISKIT:
|
||
circuit = pauli_to_qiskit(pauli, p_layers, gamma, beta)
|
||
from qiskit_aer import AerSimulator
|
||
simulator = AerSimulator()
|
||
job = simulator.run(circuit, shots=shots)
|
||
counts_dict = job.result().get_counts()
|
||
solution = measurements_to_solution(counts_dict, qubo.n)
|
||
energy = qubo.energy(solution)
|
||
result["circuit"] = circuit.qasm()
|
||
result["counts"] = counts_dict
|
||
result["solution"] = solution
|
||
result["energy"] = energy
|
||
return result
|
||
|
||
import random
|
||
random.seed(0)
|
||
all_counts: dict[str, int] = {}
|
||
for _ in range(shots):
|
||
x = [random.randint(0, 1) for _ in range(qubo.n)]
|
||
bits = "".join(str(b) for b in x)
|
||
all_counts[bits] = all_counts.get(bits, 0) + 1
|
||
|
||
solution = measurements_to_solution(all_counts, qubo.n)
|
||
energy = qubo.energy(solution)
|
||
result["circuit"] = _describe_qaoa_circuit(pauli, p_layers)
|
||
result["counts"] = all_counts
|
||
result["solution"] = solution
|
||
result["energy"] = energy
|
||
result["note"] = "random fallback — install cirq for QAOA simulation"
|
||
return result
|
||
|
||
|
||
def _cirq_counts_to_str(counts: dict, n: int) -> dict[str, int]:
|
||
"""Convert Cirq integer keyed counts to bitstring counts."""
|
||
str_counts: dict[str, int] = {}
|
||
for val, cnt in counts.items():
|
||
bits = format(val, f"0{n}b")
|
||
str_counts[bits] = cnt
|
||
return str_counts
|
||
|
||
|
||
# =========================================================================
|
||
# XVII. High-Level: Roundtrip
|
||
# =========================================================================
|
||
|
||
def qaoa_roundtrip(
|
||
receipt: dict,
|
||
p_layers: int = 1,
|
||
shots: int = 1000,
|
||
backend: str = "cirq",
|
||
compare_stochastic: bool = False,
|
||
) -> dict:
|
||
"""Full roundtrip: BraidReceipt → QAOA → updated BraidReceipt.
|
||
|
||
Args:
|
||
receipt: Input BraidReceipt JSON dict
|
||
p_layers: QAOA depth
|
||
shots: Measurement shots per layer
|
||
backend: "cirq", "qiskit", "describe"
|
||
compare_stochastic: If True, also run SA/Levy comparison
|
||
|
||
Returns:
|
||
dict with input_receipt, output_receipt, qaoa_result, delta,
|
||
and optionally stochastic comparison.
|
||
"""
|
||
qubo = braid_receipt_to_qubo(receipt)
|
||
qaoa_result = qaoa_solve_qubo(
|
||
qubo, p_layers=p_layers, shots=shots, backend=backend,
|
||
)
|
||
solution = qaoa_result["solution"]
|
||
output_receipt = solution_to_braid_receipt(solution, receipt)
|
||
|
||
input_kappas = [
|
||
s.get("kappa", 0)
|
||
for s in receipt.get("braid", receipt).get("strands", [])
|
||
]
|
||
output_kappas = [
|
||
s.get("kappa", 0)
|
||
for s in output_receipt.get("braid", output_receipt).get("strands", [])
|
||
]
|
||
|
||
result: dict[str, Any] = {
|
||
"input_receipt": receipt,
|
||
"output_receipt": output_receipt,
|
||
"qaoa_result": {
|
||
"solution": solution,
|
||
"energy": qaoa_result["energy"],
|
||
"counts": qaoa_result.get("counts", {}),
|
||
"p_layers": p_layers,
|
||
},
|
||
"delta": {
|
||
"scar_absent": {
|
||
"before": receipt.get("scar_absent"),
|
||
"after": output_receipt.get("scar_absent"),
|
||
},
|
||
"step_count": {
|
||
"before": receipt.get("step_count", 0),
|
||
"after": output_receipt.get("step_count", 0),
|
||
},
|
||
"avg_kappa": {
|
||
"before": sum(input_kappas) / max(len(input_kappas), 1),
|
||
"after": sum(output_kappas) / max(len(output_kappas), 1),
|
||
},
|
||
},
|
||
}
|
||
|
||
if compare_stochastic:
|
||
comparison = qaoa_vs_stochastic_comparison(qubo)
|
||
result["stochastic_comparison"] = comparison
|
||
|
||
return result
|
||
|
||
|
||
# =========================================================================
|
||
# XVIII-A. FinslerMetric demo
|
||
# =========================================================================
|
||
|
||
def finsler_demo() -> dict[str, Any]:
|
||
"""End-to-end demo of the FinslerMetric → QUBO → Ising → Pauli pipeline.
|
||
|
||
Constructs a simple 2D Finsler metric with nontrivial wind (β ≠ 0),
|
||
generates direction vectors on a circle, builds the QUBO matrix,
|
||
converts to Ising and Pauli strings, and verifies anisotropy.
|
||
|
||
This mirrors the Lean `trivialRanders` witness in TransportQUBOBridge.lean.
|
||
"""
|
||
# 2D Finsler metric with wind drift β = (0.25, 0.0)
|
||
metric = FinslerMetric(
|
||
alpha_mass=[1.0, 1.0], # uniform mass
|
||
beta_wind=[0.25, 0.0], # rightward drift
|
||
dimension=2,
|
||
)
|
||
|
||
# 8 directions on the unit circle
|
||
import math
|
||
n_dirs = 8
|
||
directions = []
|
||
for k in range(n_dirs):
|
||
theta = 2.0 * math.pi * k / n_dirs
|
||
directions.append([math.cos(theta), math.sin(theta)])
|
||
|
||
# Build QUBO (now stores Q_ij + Q_ji summed per pair)
|
||
qubo = finsler_metric_to_qubo(metric, directions, normalize=False)
|
||
|
||
# Check anisotropy via raw matrix (Q_ij vs Q_ji before summation)
|
||
n_raw = len(directions)
|
||
raw = [[0.0] * n_raw for _ in range(n_raw)]
|
||
for i in range(n_raw):
|
||
for j in range(n_raw):
|
||
if i != j:
|
||
raw[i][j] = metric.crossing_cost(directions[i], directions[j])
|
||
anisotropic = is_anisotropic(qubo, raw_matrix=raw)
|
||
|
||
# Convert to Ising
|
||
ising = qubo_to_ising(qubo)
|
||
|
||
# Convert to Pauli
|
||
pauli = ising_to_pauli(ising)
|
||
|
||
# Find geodesic assignment
|
||
geo = geodesic_assignment(metric, directions)
|
||
|
||
# Summary
|
||
return {
|
||
"metric": {
|
||
"alpha_mass": metric.alpha_mass,
|
||
"beta_wind": metric.beta_wind,
|
||
"dimension": metric.dimension,
|
||
},
|
||
"n_directions": n_dirs,
|
||
"qubo_is_anisotropic": anisotropic,
|
||
"qubo_n": qubo.n,
|
||
"qubo_nonzero_entries": len(qubo.matrix),
|
||
"ising_h": ising.h[:5] if ising.n > 5 else ising.h,
|
||
"ising_J_entries": len(ising.J),
|
||
"pauli_terms": len(pauli.terms),
|
||
"geodesic_assignment": geo,
|
||
"anisotropic_pair": _find_anisotropic_pair(qubo, raw_matrix=raw),
|
||
}
|
||
|
||
|
||
def _find_anisotropic_pair(qubo: QUBO, raw_matrix: Optional[list[list[float]]] = None) -> Optional[dict]:
|
||
"""Find the first (i,j) where Q_ij ≠ Q_ji."""
|
||
if raw_matrix is not None:
|
||
n = len(raw_matrix)
|
||
for i in range(n):
|
||
for j in range(n):
|
||
if i != j and abs(raw_matrix[i][j] - raw_matrix[j][i]) > 1e-9:
|
||
return {"i": i, "j": j, "Q_ij": raw_matrix[i][j], "Q_ji": raw_matrix[j][i], "delta": raw_matrix[i][j] - raw_matrix[j][i]}
|
||
return None
|
||
for (i, j), qij in qubo.matrix.items():
|
||
if i == j:
|
||
continue
|
||
qji = qubo.matrix.get((j, i), 0.0)
|
||
if abs(qij - qji) > 1e-9:
|
||
return {"i": i, "j": j, "Q_ij": qij, "Q_ji": qji, "delta": qij - qji}
|
||
return None
|
||
|
||
|
||
# =========================================================================
|
||
# XVIII. CLI
|
||
# =========================================================================
|
||
|
||
def _parse_receipt(path: str) -> dict:
|
||
with open(path) as f:
|
||
return json.load(f)
|
||
|
||
|
||
def _serialize(obj: Any) -> Any:
|
||
"""Recursively serialize a result object to JSON-safe types.
|
||
|
||
Handles QUBO, Ising, PauliSum dataclasses and dicts with tuple keys.
|
||
"""
|
||
if hasattr(obj, "__dataclass_fields__"):
|
||
d = asdict(obj)
|
||
return _serialize(d)
|
||
if isinstance(obj, dict):
|
||
out = {}
|
||
for k, v in obj.items():
|
||
if isinstance(k, tuple):
|
||
sk = f"({','.join(str(x) for x in k)})"
|
||
elif isinstance(k, int):
|
||
sk = str(k)
|
||
else:
|
||
sk = k
|
||
out[sk] = _serialize(v)
|
||
return out
|
||
if isinstance(obj, (list, tuple)):
|
||
return [_serialize(x) for x in obj]
|
||
return obj
|
||
|
||
|
||
def main() -> None:
|
||
import argparse
|
||
|
||
parser = argparse.ArgumentParser(
|
||
description="Bidirectional QAOA Conversion Adapter Set"
|
||
)
|
||
parser.add_argument(
|
||
"action",
|
||
choices=[
|
||
"lean-formulation", "lean-field", "braid-search",
|
||
"braid-to-qubo", "lone-to-qubo", "goorm-to-qubo",
|
||
"qubo-to-ising", "ising-to-pauli", "pauli-to-cirq",
|
||
"pauli-to-qiskit", "solve-qubo", "stochastic",
|
||
"compare", "roundtrip",
|
||
"tune", "auto-solve", "slos-calibrate", "auto-calibrate",
|
||
"finsler-demo",
|
||
],
|
||
)
|
||
parser.add_argument("--coeffs", help="JSON file with BraidCostCoeffs overrides")
|
||
parser.add_argument("--receipt", help="Path to BraidReceipt JSON")
|
||
parser.add_argument("--output", "-o", help="Output path (default: stdout)")
|
||
parser.add_argument("--p-layers", type=int, default=1, help="QAOA layers")
|
||
parser.add_argument("--shots", type=int, default=1000, help="Measurement shots")
|
||
parser.add_argument("--backend", default="cirq", help="Simulator backend")
|
||
parser.add_argument("--gamma", type=float, nargs="*", help="Cost angles")
|
||
parser.add_argument("--beta", type=float, nargs="*", help="Mixer angles")
|
||
parser.add_argument("--method", default="sa", help="Stochastic method")
|
||
parser.add_argument("--time-limit", type=float, default=5.0,
|
||
help="Stochastic solver time limit")
|
||
|
||
args = parser.parse_args()
|
||
|
||
result: Any = None
|
||
|
||
if args.action == "lean-formulation":
|
||
result = {"note": "Use Python API with lean_qubo_formulation_to_qubo(matrix, n)"}
|
||
elif args.action == "lean-field":
|
||
result = lean_qubo_field_to_qubo(65536, 10 * 65536)
|
||
elif args.action == "braid-search":
|
||
if args.receipt:
|
||
receipt = _parse_receipt(args.receipt)
|
||
bkts = receipt.get("braid", receipt).get("strands", [])
|
||
result = braid_search_brackets_to_qubo(bkts)
|
||
else:
|
||
result = {"note": "Use --receipt with bracket-containing JSON"}
|
||
elif args.action == "braid-to-qubo":
|
||
receipt = _parse_receipt(args.receipt)
|
||
result = braid_receipt_to_qubo(receipt)
|
||
elif args.action == "lone-to-qubo":
|
||
result = {"note": "Use Python API with a numpy scar field array"}
|
||
elif args.action == "goorm-to-qubo":
|
||
result = goormaghtigh_cost_to_qubo()
|
||
elif args.action == "qubo-to-ising":
|
||
result = {"note": "Use Python API with a QUBO object"}
|
||
elif args.action == "ising-to-pauli":
|
||
result = {"note": "Use Python API with an Ising object"}
|
||
elif args.action == "pauli-to-cirq":
|
||
result = {"note": "Use Python API with a PauliSum object"}
|
||
elif args.action == "pauli-to-qiskit":
|
||
result = {"note": "Use Python API with a PauliSum object"}
|
||
elif args.action == "solve-qubo":
|
||
receipt = _parse_receipt(args.receipt)
|
||
result = qaoa_roundtrip(receipt, args.p_layers, args.shots, args.backend)
|
||
elif args.action == "stochastic":
|
||
receipt = _parse_receipt(args.receipt)
|
||
qubo = braid_receipt_to_qubo(receipt)
|
||
result = stochastic_abuse_qubo(qubo, args.method, args.time_limit)
|
||
elif args.action == "compare":
|
||
receipt = _parse_receipt(args.receipt)
|
||
qubo = braid_receipt_to_qubo(receipt)
|
||
result = qaoa_vs_stochastic_comparison(
|
||
qubo, args.p_layers, args.shots, args.backend,
|
||
time_limit=args.time_limit,
|
||
)
|
||
elif args.action == "tune":
|
||
receipt = _parse_receipt(args.receipt)
|
||
bkts = receipt.get("braid", receipt).get("strands", [])
|
||
if args.coeffs:
|
||
custom_grid = [json.loads(Path(args.coeffs).read_text())]
|
||
result = tune_braid_parameters(bkts, param_grid=custom_grid)
|
||
else:
|
||
result = tune_braid_parameters(bkts, time_per_trial=args.time_limit)
|
||
elif args.action == "auto-solve":
|
||
receipt = _parse_receipt(args.receipt)
|
||
bkts = receipt.get("braid", receipt).get("strands", [])
|
||
Q = build_tunable_qubo_matrix(bkts)
|
||
qubo = QUBO(n=len(bkts), matrix=Q)
|
||
result = solver_dispatch_tuner(qubo, time_limit=args.time_limit)
|
||
elif args.action == "slos-calibrate":
|
||
receipt = _parse_receipt(args.receipt)
|
||
bkts = receipt.get("braid", receipt).get("strands", [])
|
||
evaluator = SlosBraidEvaluator()
|
||
result = slos_calibrate_coeffs(bkts, evaluator=evaluator)
|
||
result["n"] = len(bkts)
|
||
result["note"] = (
|
||
"SLOS-calibrated coefficients. Set as braid_search defaults"
|
||
" by copying calibrated_coeffs into BraidCostCoeffs() args."
|
||
)
|
||
elif args.action == "auto-calibrate":
|
||
receipt = _parse_receipt(args.receipt)
|
||
bkts = receipt.get("braid", receipt).get("strands", [])
|
||
evaluator = SlosBraidEvaluator()
|
||
result = auto_calibrate(bkts, evaluator=evaluator, time_per_solver=args.time_limit)
|
||
elif args.action == "roundtrip":
|
||
receipt = _parse_receipt(args.receipt)
|
||
result = qaoa_roundtrip(
|
||
receipt, args.p_layers, args.shots, args.backend,
|
||
compare_stochastic=True,
|
||
)
|
||
elif args.action == "finsler-demo":
|
||
result = finsler_demo()
|
||
|
||
output = json.dumps(_serialize(result), indent=2, default=str)
|
||
|
||
if args.output:
|
||
Path(args.output).write_text(output)
|
||
else:
|
||
print(output)
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|