From 7d3e90331a7f809a6194239b3a12ad848cf15ee0 Mon Sep 17 00:00:00 2001 From: allaun Date: Sun, 21 Jun 2026 01:02:01 -0500 Subject: [PATCH] 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. --- .../shim/benchmark_finsler_qaoa.py | 154 ++++++++++++++++++ 4-Infrastructure/shim/qaoa_adapter.py | 54 ++++-- 2 files changed, 197 insertions(+), 11 deletions(-) create mode 100644 4-Infrastructure/shim/benchmark_finsler_qaoa.py diff --git a/4-Infrastructure/shim/benchmark_finsler_qaoa.py b/4-Infrastructure/shim/benchmark_finsler_qaoa.py new file mode 100644 index 00000000..8eb2f1af --- /dev/null +++ b/4-Infrastructure/shim/benchmark_finsler_qaoa.py @@ -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) diff --git a/4-Infrastructure/shim/qaoa_adapter.py b/4-Infrastructure/shim/qaoa_adapter.py index 2ff8f597..c0e7a686 100644 --- a/4-Infrastructure/shim/qaoa_adapter.py +++ b/4-Infrastructure/shim/qaoa_adapter.py @@ -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