mirror of
https://github.com/allaunthefox/Research-Stack.git
synced 2026-07-31 03:05:21 +00:00
feat(infra): Finsler QAOA adapter fixes and benchmark shim
- 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.
This commit is contained in:
parent
902191cfbb
commit
7d3e90331a
2 changed files with 197 additions and 11 deletions
154
4-Infrastructure/shim/benchmark_finsler_qaoa.py
Normal file
154
4-Infrastructure/shim/benchmark_finsler_qaoa.py
Normal file
|
|
@ -0,0 +1,154 @@
|
|||
#!/usr/bin/env python3
|
||||
"""Benchmark: Finsler QUBO — QAOA vs classical solvers across problem sizes.
|
||||
|
||||
Compares QAOA (Cirq), Simulated Annealing (SA), and Levy flights on
|
||||
TransportQUBOBridge Finsler-Randers QUBOs at n=8, 12, 24, 48 directions.
|
||||
"""
|
||||
|
||||
import sys
|
||||
import os
|
||||
sys.path.insert(0, os.path.dirname(__file__))
|
||||
|
||||
import json
|
||||
import math
|
||||
import time
|
||||
import numpy as np
|
||||
|
||||
from qaoa_adapter import (
|
||||
FinslerMetric, finsler_metric_to_qubo, qubo_to_ising,
|
||||
ising_to_pauli, pauli_to_cirq, qaoa_solve_qubo,
|
||||
stochastic_abuse_qubo, qaoa_vs_stochastic_comparison,
|
||||
)
|
||||
|
||||
def generate_finsler_qubo(n_dirs: int, dimension: int = 3, seed: int = 42) -> tuple:
|
||||
rng = np.random.RandomState(seed)
|
||||
alpha_mass = rng.uniform(0.5, 2.0, size=dimension).tolist()
|
||||
beta_wind = rng.uniform(-0.3, 0.3, size=dimension).tolist()
|
||||
metric = FinslerMetric(alpha_mass=alpha_mass, beta_wind=beta_wind, dimension=dimension)
|
||||
directions = []
|
||||
for k in range(n_dirs):
|
||||
vec = rng.randn(dimension)
|
||||
norm = math.sqrt(sum(v*v for v in vec))
|
||||
directions.append([v / norm for v in vec])
|
||||
qubo = finsler_metric_to_qubo(metric, directions, normalize=True)
|
||||
# Verify anisotropy
|
||||
raw = [[0.0]*n_dirs for _ in range(n_dirs)]
|
||||
for i in range(n_dirs):
|
||||
for j in range(n_dirs):
|
||||
if i != j:
|
||||
raw[i][j] = metric.crossing_cost(directions[i], directions[j])
|
||||
max_aniso = max(abs(raw[i][j] - raw[j][i]) for i in range(n_dirs) for j in range(n_dirs))
|
||||
return qubo, metric, directions, max_aniso
|
||||
|
||||
|
||||
def run_benchmark(sizes: list[int] = None, p_layers: int = 1, shots: int = 2000, time_limit: float = 5.0, seed: int = 42) -> dict:
|
||||
if sizes is None:
|
||||
sizes = [8, 12, 24, 48]
|
||||
results = {}
|
||||
for n in sizes:
|
||||
print(f"\n=== Finsler QUBO n={n} ===")
|
||||
qubo, metric, directions, max_aniso = generate_finsler_qubo(n, seed=seed)
|
||||
print(f" max |Q_ij - Q_ji| = {max_aniso:.6f} {'anisotropic' if max_aniso > 1e-9 else 'symmetric'}")
|
||||
row = {"n": n, "max_anisotropy": max_aniso}
|
||||
# QAOA
|
||||
try:
|
||||
t0 = time.time()
|
||||
qaoa_result = qaoa_solve_qubo(
|
||||
qubo, p_layers=p_layers, shots=shots, backend="cirq"
|
||||
)
|
||||
t_qaoa = time.time() - t0
|
||||
row["qaoa"] = {
|
||||
"energy": round(qaoa_result["energy"], 6),
|
||||
"time_s": round(t_qaoa, 3),
|
||||
"solution": qaoa_result["solution"],
|
||||
"p_layers": p_layers,
|
||||
}
|
||||
print(f" QAOA(p={p_layers}): energy={row['qaoa']['energy']}, time={t_qaoa:.3f}s")
|
||||
except Exception as e:
|
||||
row["qaoa"] = {"error": str(e)}
|
||||
print(f" QAOA: FAILED — {e}")
|
||||
# SA
|
||||
try:
|
||||
t0 = time.time()
|
||||
sa_result = stochastic_abuse_qubo(qubo, method="sa", time_limit=time_limit, seed=seed)
|
||||
t_sa = time.time() - t0
|
||||
row["sa"] = {
|
||||
"energy": round(sa_result["energy"], 6),
|
||||
"time_s": round(t_sa, 3),
|
||||
"solution": sa_result["solution"],
|
||||
}
|
||||
print(f" SA: energy={row['sa']['energy']}, time={t_sa:.3f}s")
|
||||
except Exception as e:
|
||||
row["sa"] = {"error": str(e)}
|
||||
print(f" SA: FAILED — {e}")
|
||||
# Levy flights
|
||||
try:
|
||||
t0 = time.time()
|
||||
levy_result = stochastic_abuse_qubo(qubo, method="levy", time_limit=time_limit, seed=seed)
|
||||
t_levy = time.time() - t0
|
||||
row["levy"] = {
|
||||
"energy": round(levy_result["energy"], 6),
|
||||
"time_s": round(t_levy, 3),
|
||||
"solution": levy_result["solution"],
|
||||
}
|
||||
print(f" Levy flights: energy={row['levy']['energy']}, time={t_levy:.3f}s")
|
||||
except Exception as e:
|
||||
row["levy"] = {"error": str(e)}
|
||||
print(f" Levy: FAILED — {e}")
|
||||
# Geodesic (ground truth)
|
||||
geodesic = geodesic_assignment(metric, directions)
|
||||
geo_energy = qubo.energy([1 if g else 0 for g in geodesic])
|
||||
all_false_energy = qubo.energy([0]*n)
|
||||
row["geodesic"] = {
|
||||
"energy": round(geo_energy, 6),
|
||||
"all_false_energy": round(all_false_energy, 6),
|
||||
"assignment": [int(g) for g in geodesic],
|
||||
}
|
||||
print(f" Geodesic truth: energy={geo_energy:.6f} (all-false: {all_false_energy:.6f})")
|
||||
# Determine winner
|
||||
energies = {}
|
||||
if "energy" in row.get("qaoa", {}):
|
||||
energies["qaoa"] = row["qaoa"]["energy"]
|
||||
if "energy" in row.get("sa", {}):
|
||||
energies["sa"] = row["sa"]["energy"]
|
||||
if "energy" in row.get("levy", {}):
|
||||
energies["levy"] = row["levy"]["energy"]
|
||||
if energies:
|
||||
best_name = min(energies, key=energies.get)
|
||||
row["winner"] = {"solver": best_name, "energy": energies[best_name]}
|
||||
row["gap_vs_geodesic"] = round(energies[best_name] - geo_energy, 6)
|
||||
print(f" Winner: {best_name} (gap vs geodesic: {row['gap_vs_geodesic']})")
|
||||
results[f"n_{n}"] = row
|
||||
return results
|
||||
|
||||
|
||||
def geodesic_assignment(metric, directions, tol=1e-9):
|
||||
costs = [metric.finsler_cost(v) for v in directions]
|
||||
min_cost = min(costs)
|
||||
return [abs(c - min_cost) <= tol for c in costs]
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
import argparse
|
||||
parser = argparse.ArgumentParser(description="Benchmark Finsler QUBO solvers")
|
||||
parser.add_argument("--sizes", type=int, nargs="+", default=[8, 12, 24, 48])
|
||||
parser.add_argument("--p-layers", type=int, default=1)
|
||||
parser.add_argument("--shots", type=int, default=2000)
|
||||
parser.add_argument("--time-limit", type=float, default=5.0)
|
||||
parser.add_argument("--seed", type=int, default=42)
|
||||
parser.add_argument("--output", "-o", help="JSON output path")
|
||||
args = parser.parse_args()
|
||||
results = run_benchmark(
|
||||
sizes=args.sizes,
|
||||
p_layers=args.p_layers,
|
||||
shots=args.shots,
|
||||
time_limit=args.time_limit,
|
||||
seed=args.seed,
|
||||
)
|
||||
output = json.dumps(results, indent=2)
|
||||
if args.output:
|
||||
with open(args.output, "w") as f:
|
||||
f.write(output)
|
||||
print(f"\nWrote results to {args.output}")
|
||||
else:
|
||||
print(output)
|
||||
|
|
@ -1589,21 +1589,36 @@ def finsler_metric_to_qubo(
|
|||
# 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, n):
|
||||
if i == j:
|
||||
continue
|
||||
val = raw[i][j] / scale
|
||||
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, tol: float = 1e-9) -> bool:
|
||||
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:
|
||||
|
|
@ -1716,6 +1731,7 @@ def pauli_to_cirq(
|
|||
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.
|
||||
|
||||
|
|
@ -1760,6 +1776,9 @@ def pauli_to_cirq(
|
|||
|
||||
circuit.append(cirq.rx(2.0 * b).on_each(*qubits))
|
||||
|
||||
if measure:
|
||||
circuit.append(cirq.measure(*qubits, key="result"))
|
||||
|
||||
return circuit
|
||||
|
||||
|
||||
|
|
@ -2271,7 +2290,7 @@ def qaoa_solve_qubo(
|
|||
return result
|
||||
|
||||
if backend == "cirq" and _HAS_CIRQ:
|
||||
circuit = pauli_to_cirq(pauli, p_layers, gamma, beta)
|
||||
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")
|
||||
|
|
@ -2425,11 +2444,17 @@ def finsler_demo() -> dict[str, Any]:
|
|||
theta = 2.0 * math.pi * k / n_dirs
|
||||
directions.append([math.cos(theta), math.sin(theta)])
|
||||
|
||||
# Build QUBO
|
||||
# Build QUBO (now stores Q_ij + Q_ji summed per pair)
|
||||
qubo = finsler_metric_to_qubo(metric, directions, normalize=False)
|
||||
|
||||
# Check anisotropy
|
||||
anisotropic = is_anisotropic(qubo)
|
||||
# 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)
|
||||
|
|
@ -2455,12 +2480,19 @@ def finsler_demo() -> dict[str, Any]:
|
|||
"ising_J_entries": len(ising.J),
|
||||
"pauli_terms": len(pauli.terms),
|
||||
"geodesic_assignment": geo,
|
||||
"anisotropic_pair": _find_anisotropic_pair(qubo),
|
||||
"anisotropic_pair": _find_anisotropic_pair(qubo, raw_matrix=raw),
|
||||
}
|
||||
|
||||
|
||||
def _find_anisotropic_pair(qubo: QUBO) -> Optional[dict]:
|
||||
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
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue