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.

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)

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):
+        for j in range(i + 1, n):
-            if i == j:
+            val = (raw[i][j] + raw[j][i]) / scale
-                continue
-            val = raw[i][j] / 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
+    # Check anisotropy via raw matrix (Q_ij vs Q_ji before summation)
-    anisotropic = is_anisotropic(qubo)
+    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