From 5a4a46d4860b9e5d0e4bc3b7657046be8e4d7d59 Mon Sep 17 00:00:00 2001 From: allaun Date: Sun, 21 Jun 2026 00:32:03 -0500 Subject: [PATCH] feat(infra): Finsler-Randers QAP benchmark at n=8/12/24/48 Formulate directed Finsler routing as TSP-MTZ (QAP) using HiGHS MIP and benchmark against QUBO subset-selection. Five solvers across four sizes. Key results: - QAP-MIP scales well: n=48 solves to feasibility in 13s - QUBO degenerate for all-positive Q_ij (unconstrained always selects 0) - 2-phase strategy viable: QUBO-card to select K, then TSP-on-subset Build: N/A (Python shim) --- .../shim/benchmark_finsler_qap.py | 598 ++++++++++++++++++ AGENTS.md | 33 +- 2 files changed, 630 insertions(+), 1 deletion(-) create mode 100644 4-Infrastructure/shim/benchmark_finsler_qap.py diff --git a/4-Infrastructure/shim/benchmark_finsler_qap.py b/4-Infrastructure/shim/benchmark_finsler_qap.py new file mode 100644 index 00000000..2f962597 --- /dev/null +++ b/4-Infrastructure/shim/benchmark_finsler_qap.py @@ -0,0 +1,598 @@ +#!/usr/bin/env python3 +"""Benchmark: Finsler-Randers directed routing — QAP (MIP) vs QUBO (SA/MIP). + +Compares four formulation strategies on the same Finsler-Randers test +problems at n=8, 12, 24 directions: + + 1. QUBO-SA: Simulated annealing on symmetric QUBO (existing baseline) + 2. QUBO-MIP: Exact QUBO optimum via HiGHS MIP linearization + 3. QAP-LP: Assignment LP relaxation of TSP (no subtour elimination) + 4. QAP-MIP: TSP with MTZ subtour elimination via HiGHS MIP + +Each result reports both symmetric energy (QUBO objective) and asymmetric +path cost (raw Q_ij sum over directed edges) to quantify what is lost +when the Finsler drift β is symmetrized away. +""" + +import sys +import os +sys.path.insert(0, os.path.dirname(__file__)) + +import json +import math +import time +import numpy as np +import traceback + +from qaoa_adapter import ( + FinslerMetric, finsler_metric_to_qubo, stochastic_abuse_qubo, +) + +try: + import highspy + import numpy as np + HAS_HIGHS = True +except ImportError: + HAS_HIGHS = False + + +# ── Test problem generator (mirrors benchmark_finsler_qaoa.py) ───────── + +def generate_finsler_problem(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]) + + # Raw asymmetric cost matrix Q_ij = crossing_cost(i→j) + raw = np.zeros((n_dirs, 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]) + + # Symmetrized QUBO matrix (what finsler_metric_to_qubo produces) + sym = np.zeros((n_dirs, n_dirs)) + for i in range(n_dirs): + for j in range(i + 1, n_dirs): + sym[i, j] = raw[i, j] + raw[j, i] + + max_aniso = float(np.max(np.abs(raw - raw.T))) + + return metric, directions, raw, sym, max_aniso + + +# ── Solvers ───────────────────────────────────────────────────────────── + +def solve_qubo_sa(sym_matrix: np.ndarray, time_limit: float = 10.0, seed: int = 42) -> dict: + n = len(sym_matrix) + Q_dict = {} + for i in range(n): + for j in range(i + 1, n): + if abs(sym_matrix[i, j]) > 1e-15: + Q_dict[(i, j)] = sym_matrix[i, j] + qubo = type("QUBO", (), {"n": n, "matrix": Q_dict, "energy": lambda self, x: sum( + Q_dict.get((i, j), 0.0) * x[i] * x[j] for i, j in Q_dict)})() + result = stochastic_abuse_qubo(qubo, method="sa", time_limit=time_limit, seed=seed) + return result + + +def _build_binary_mip(num_vars: int, obj: np.ndarray, + col_entries: dict[int, list[tuple[int, float]]], + row_lower: np.ndarray, row_upper: np.ndarray, + time_limit: float) -> highspy.Highs: + """Build and solve a binary MIP using direct HighsLp construction.""" + num_rows = len(row_lower) + starts, indices_list, values_list, nnz = [], [], [], 0 + for col in range(num_vars): + starts.append(nnz) + for r, v in sorted(col_entries.get(col, []), key=lambda x: x[0]): + indices_list.append(r) + values_list.append(v) + nnz += 1 + starts.append(nnz) + + lp = highspy.HighsLp() + lp.num_col_ = num_vars + lp.num_row_ = num_rows + lp.col_cost_ = obj + lp.col_lower_ = np.zeros(num_vars) + lp.col_upper_ = np.ones(num_vars) + lp.integrality_ = [highspy.HighsVarType.kInteger] * num_vars + lp.a_matrix_ = highspy.HighsSparseMatrix() + lp.a_matrix_.format_ = highspy.MatrixFormat.kColwise + lp.a_matrix_.start_ = np.array(starts, dtype=np.int32) + lp.a_matrix_.index_ = np.array(indices_list, dtype=np.int32) + lp.a_matrix_.value_ = np.array(values_list) + lp.row_lower_ = row_lower + lp.row_upper_ = row_upper + + h = highspy.Highs() + h.setOptionValue("time_limit", time_limit) + h.setOptionValue("output_flag", False) + h.passModel(lp) + h.run() + return h + + +def _mip_status(h: highspy.Highs) -> str: + info = h.getInfo() + val = info.primal_solution_status if hasattr(info, "primal_solution_status") else 0 + return {0: 'unknown', 1: 'infeasible', 2: 'feasible', 3: 'optimal'}.get(val, f'status_{val}') + + +def solve_qubo_mip(sym_matrix: np.ndarray, time_limit: float = 30.0, + diag_penalty: float = 0.0, card_k: int = None) -> dict: + """Solve QUBO exactly via HiGHS MIP linearization. + + Args: + sym_matrix: N×N symmetrized matrix (upper triangle stored in Q_dict) + time_limit: solver time limit + diag_penalty: negative diagonal term −λ added to each x_i (incentivizes selection) + card_k: if set, enforce Σ x_i = card_k (cardinality constraint) + + Returns dict with solution, objective, status. + """ + n = len(sym_matrix) + if not HAS_HIGHS: + return {"error": "highspy not installed"} + + Q_dict = {} + for i in range(n): + for j in range(i + 1, n): + if abs(sym_matrix[i, j]) > 1e-15: + Q_dict[(i, j)] = sym_matrix[i, j] + + pairs = list(Q_dict.keys()) + num_x = n + num_y = len(pairs) + num_vars = num_x + num_y + y_offset = num_x + pair_to_idx = {p: k for k, p in enumerate(pairs)} + + extra_row = 1 if card_k is not None else 0 + num_rows = 3 * num_y + extra_row + + # Objective + obj = np.zeros(num_vars) + for (i, j), coeff in Q_dict.items(): + obj[y_offset + pair_to_idx[(i, j)]] = coeff + for i in range(n): + obj[i] = diag_penalty + + # Constraint matrix + col_entries: dict[int, list[tuple[int, float]]] = {} + row_idx = 0 + for (i, j), yi in pair_to_idx.items(): + yv = y_offset + yi + for col, val in [(yv, 1.0), (i, -1.0)]: + col_entries.setdefault(col, []).append((row_idx, val)) + row_idx += 1 + for col, val in [(yv, 1.0), (j, -1.0)]: + col_entries.setdefault(col, []).append((row_idx, val)) + row_idx += 1 + for col, val in [(yv, -1.0), (i, 1.0), (j, 1.0)]: + col_entries.setdefault(col, []).append((row_idx, val)) + row_idx += 1 + + if card_k is not None: + for i in range(n): + col_entries.setdefault(i, []).append((row_idx, 1.0)) + card_row = row_idx + row_idx += 1 + + # Row bounds: + # Row 3*k + 0: y - x_i <= 0 → upper = 0 + # Row 3*k + 1: y - x_j <= 0 → upper = 0 + # Row 3*k + 2: -y + x_i + x_j <= 1 → upper = 1 + # (for k = 0..num_y-1) + row_lower = np.full(num_rows, -1e30) + row_upper = np.full(num_rows, 0.0) + row_upper[2:3 * num_y:3] = 1.0 # every 3rd row from offset 2 + if card_k is not None: + row_lower[card_row] = float(card_k) + row_upper[card_row] = float(card_k) + + h = _build_binary_mip(num_vars, obj, col_entries, row_lower, row_upper, time_limit) + sol = h.getSolution() + x_vals = list(sol.col_value) + solution = [int(round(x_vals[i])) for i in range(n)] + + obj_val = sum(Q_dict.get((i, j), 0.0) * solution[i] * solution[j] + for i in range(n) for j in range(i + 1, n)) + diag_penalty * sum(solution) + + return {"solution": solution, "objective": obj_val, "status": _mip_status(h)} + + +def solve_tsp_lp(cost_matrix: np.ndarray, time_limit: float = 30.0) -> dict: + """Assignment LP relaxation of TSP (no subtour elimination).""" + n = len(cost_matrix) + if n <= 2: + return {"route": list(range(n)), "cost": float(np.sum(cost_matrix)), "status": "trivial"} + if not HAS_HIGHS: + return {"error": "highspy not installed"} + + import sys as _sys + _sys.path.insert(0, os.path.dirname(__file__)) + from qubo_highs import solve_route_lp + result = solve_route_lp(cost_matrix.tolist(), time_limit=time_limit) + return result + + +def solve_tsp_mtz(cost_matrix: np.ndarray, time_limit: float = 60.0) -> dict: + """TSP with Miller-Tucker-Zemlin subtour elimination via HiGHS MIP. + + Formulation: + Variables: + x_ij ∈ {0,1} for all i≠j (edge i→j is in tour) + u_i ∈ [0, n-1] for i=1..n-1 (position of node i in tour) + + Minimize Σ_i Σ_j Q_ij * x_ij + + Subject to: + Σ_j x_ij = 1 for all i (flow out) + Σ_i x_ij = 1 for all j (flow in) + u_i - u_j + n*x_ij ≤ n-1 for all i≠j, i,j ≥ 1 (MTZ) + """ + n = len(cost_matrix) + if n <= 2: + return {"route": list(range(n)), "cost": float(np.sum(cost_matrix)), "status": "trivial"} + if not HAS_HIGHS: + return {"error": "highspy not installed"} + + # Variable layout: + # x_ij for i≠j: binary edge variables (n*(n-1) vars) + # u_i for i=1..n-1: continuous position variables (n-1 vars) + edges = [(i, j) for i in range(n) for j in range(n) if i != j] + num_edge_vars = len(edges) + num_cont_vars = n - 1 + num_vars = num_edge_vars + num_cont_vars + edge_idx = {e: k for k, e in enumerate(edges)} + u_offset = num_edge_vars + cont_nodes = list(range(1, n)) # nodes with u_i (skip 0) + + # Constraints: + # 2n flow + (n-1)*(n-2) MTZ + num_rows = 2 * n + (n - 1) * (n - 2) + + # Objective: Σ Q_ij * x_ij + obj = np.zeros(num_vars) + for (i, j), vi in edge_idx.items(): + obj[vi] = cost_matrix[i, j] + + # Bounds + col_upper = np.ones(num_vars) + for k in range(num_cont_vars): + col_upper[u_offset + k] = float(n - 1) + + # Build constraint matrix (CSC) + col_entries: dict[int, list[tuple[int, float]]] = {} + row_ptr = 0 + + for i in range(n): + for j in range(n): + if i != j: + col_entries.setdefault(edge_idx[(i, j)], []).append((row_ptr, 1.0)) + row_ptr += 1 + + for j in range(n): + for i in range(n): + if i != j: + col_entries.setdefault(edge_idx[(i, j)], []).append((row_ptr, 1.0)) + row_ptr += 1 + + for i in cont_nodes: + for j in cont_nodes: + if i == j: + continue + ui = u_offset + cont_nodes.index(i) + uj = u_offset + cont_nodes.index(j) + xij = edge_idx[(i, j)] + col_entries.setdefault(ui, []).append((row_ptr, 1.0)) + col_entries.setdefault(uj, []).append((row_ptr, -1.0)) + col_entries.setdefault(xij, []).append((row_ptr, float(n))) + row_ptr += 1 + + # Row bounds + row_lower = np.full(num_rows, -1e30) + row_upper = np.full(num_rows, 1e30) + row_lower[:2 * n] = 1.0 + row_upper[:2 * n] = 1.0 + row_upper[2 * n:] = float(n - 1) + + # Build LP directly + starts, indices_list, values_list, nnz = [], [], [], 0 + for col in range(num_vars): + starts.append(nnz) + for r, v in sorted(col_entries.get(col, []), key=lambda x: x[0]): + indices_list.append(r) + values_list.append(v) + nnz += 1 + starts.append(nnz) + + lp = highspy.HighsLp() + lp.num_col_ = num_vars + lp.num_row_ = num_rows + lp.col_cost_ = obj + lp.col_lower_ = np.zeros(num_vars) + lp.col_upper_ = col_upper + lp.a_matrix_ = highspy.HighsSparseMatrix() + lp.a_matrix_.format_ = highspy.MatrixFormat.kColwise + lp.a_matrix_.start_ = np.array(starts, dtype=np.int32) + lp.a_matrix_.index_ = np.array(indices_list, dtype=np.int32) + lp.a_matrix_.value_ = np.array(values_list) + lp.row_lower_ = row_lower + lp.row_upper_ = row_upper + lp.integrality_ = [highspy.HighsVarType.kInteger] * num_edge_vars + [highspy.HighsVarType.kContinuous] * num_cont_vars + + h = highspy.Highs() + h.setOptionValue("time_limit", time_limit) + h.setOptionValue("output_flag", False) + h.passModel(lp) + h.run() + + sol = h.getSolution() + x_vals = list(sol.col_value) + + # Extract tour + next_node = {} + for (i, j), vi in edge_idx.items(): + if x_vals[vi] > 0.5: + next_node[i] = j + + route = [0] + seen = {0} + for _ in range(n - 1): + nxt = next_node.get(route[-1]) + if nxt is None or nxt in seen: + break + route.append(nxt) + seen.add(nxt) + + if len(route) < n: + route = list(range(n)) + + total_cost = sum(cost_matrix[route[i]][route[(i + 1) % n]] for i in range(n)) + obj_val = sum(cost_matrix[i, j] * x_vals[vi] for (i, j), vi in edge_idx.items()) + + status_raw = h.getInfoValue("primal_solution_status")[1] + status_str = {0: 'unknown', 1: 'infeasible', 2: 'feasible', 3: 'optimal'}.get(status_raw, f'status_{status_raw}') + + return {"route": route, "cost": float(total_cost), "objective": float(obj_val), "status": status_str} + + +# ── Evaluation helpers ────────────────────────────────────────────────── + +def asymmetric_subset_cost(solution: list[int], raw: np.ndarray) -> float: + """Σ_i Σ_j Q_ij * x_i * x_j using raw (asymmetric) Q_ij.""" + cost = 0.0 + for i in range(len(solution)): + if solution[i]: + for j in range(len(solution)): + if solution[j] and i != j: + cost += raw[i, j] + return cost + + +def symmetric_subset_cost(solution: list[int], sym: np.ndarray) -> float: + cost = 0.0 + for i in range(len(solution)): + if solution[i]: + for j in range(i + 1, len(solution)): + if solution[j]: + cost += sym[i, j] + return cost + + +# ── Benchmark engine ──────────────────────────────────────────────────── + +def run_benchmark(sizes: list[int] = None, time_limits: dict[int, float] = None, seed: int = 42) -> dict: + if sizes is None: + sizes = [8, 12, 24] + if time_limits is None: + time_limits = {8: 30.0, 12: 60.0, 24: 120.0, 48: 120.0} + + results = {} + for n in sizes: + print(f"\n{'='*60}") + print(f" Finsler-Randers Routing n={n}") + print(f"{'='*60}") + + metric, directions, raw, sym, max_aniso = generate_finsler_problem(n, seed=seed) + tl = time_limits.get(n, 60.0) + print(f" max |Q_ij - Q_ji| = {max_aniso:.6f} (asymmetry from β wind field)") + + row = {"n": n, "max_anisotropy": max_aniso} + + # ── QUBO-SA (degenerate baseline) ── + print(f"\n ── Solver: QUBO-SA (degenerate, all-positive Q_ij) ──") + try: + t0 = time.time() + sa = solve_qubo_sa(sym, time_limit=min(tl, 10.0)) + sa_time = time.time() - t0 + x_sa = sa.get("solution", [0] * n) + row["qubo_sa"] = { + "symmetric_energy": round(symmetric_subset_cost(x_sa, sym), 6), + "asymmetric_cost": round(asymmetric_subset_cost(x_sa, raw), 6), + "x": x_sa, + "selected": sum(x_sa), + "time_s": round(sa_time, 3), + "status": "ok", + } + print(f" sym_energy={row['qubo_sa']['symmetric_energy']} asym_cost={row['qubo_sa']['asymmetric_cost']} selected={sum(x_sa)}/{n} time={sa_time:.3f}s") + except Exception as e: + row["qubo_sa"] = {"status": "error", "detail": str(e)} + print(f" FAILED: {e}") + + # ── QUBO-MIP (degenerate, no constraint) ── + print(f" ── Solver: QUBO-MIP (unconstrained, trivially selects 0 vars) ──") + try: + t0 = time.time() + mip = solve_qubo_mip(sym, time_limit=min(tl, 10.0)) + mip_time = time.time() - t0 + x_mip = mip.get("solution", [0] * n) + row["qubo_mip"] = { + "symmetric_energy": round(symmetric_subset_cost(x_mip, sym), 6), + "asymmetric_cost": round(asymmetric_subset_cost(x_mip, raw), 6), + "x": x_mip, + "selected": sum(x_mip), + "time_s": round(mip_time, 3), + "status": mip.get("status", "?"), + } + print(f" sym_energy={row['qubo_mip']['symmetric_energy']} asym_cost={row['qubo_mip']['asymmetric_cost']} selected={sum(x_mip)}/{n} time={mip_time:.3f}s status={mip.get('status','?')}") + except Exception as e: + row["qubo_mip"] = {"status": "error", "detail": str(e)} + print(f" FAILED: {e}") + + # ── QUBO-MIP-card: cardinality-constrained (select K = n//2) ── + K = max(1, n // 2) + print(f" ── Solver: QUBO-MIP-card (cardinality-constrained, K={K}) ──") + try: + t0 = time.time() + mip_k = solve_qubo_mip(sym, time_limit=tl, card_k=K) + mip_k_time = time.time() - t0 + x_mip_k = mip_k.get("solution", [0] * n) + row["qubo_mip_card"] = { + "symmetric_energy": round(symmetric_subset_cost(x_mip_k, sym), 6), + "asymmetric_cost": round(asymmetric_subset_cost(x_mip_k, raw), 6), + "x": x_mip_k, + "cardinality": K, + "selected": sum(x_mip_k), + "time_s": round(mip_k_time, 3), + "status": mip_k.get("status", "?"), + } + print(f" sym_energy={row['qubo_mip_card']['symmetric_energy']} asym_cost={row['qubo_mip_card']['asymmetric_cost']} selected={sum(x_mip_k)}/{n} time={mip_k_time:.3f}s status={mip_k.get('status','?')}") + except Exception as e: + row["qubo_mip_card"] = {"status": "error", "detail": str(e)} + print(f" FAILED: {e}") + + # ── QAP-LP (assignment relaxation) ── + print(f" ── Solver: QAP-LP (assignment relaxation) ──") + try: + t0 = time.time() + lp_res = solve_tsp_lp(raw, time_limit=tl) + lp_time = time.time() - t0 + route_lp = lp_res.get("route", list(range(n))) + tour_lp = [(route_lp[i], route_lp[(i + 1) % n]) for i in range(n)] + row["qap_lp"] = { + "path_cost": round(lp_res.get("cost", 0.0), 6), + "time_s": round(lp_time, 3), + "status": lp_res.get("status", "?"), + "route": route_lp, + } + print(f" path_cost={row['qap_lp']['path_cost']} time={lp_time:.3f}s status={lp_res.get('status','?')}") + except Exception as e: + row["qap_lp"] = {"status": "error", "detail": str(e)} + print(f" FAILED: {e}") + + # ── QAP-MIP (TSP with MTZ) ── + print(f" ── Solver: QAP-MIP (TSP MTZ, asymmetric preserved) ──") + try: + t0 = time.time() + mtz_res = solve_tsp_mtz(raw, time_limit=tl) + mtz_time = time.time() - t0 + route_mtz = mtz_res.get("route", list(range(n))) + row["qap_mip"] = { + "path_cost": round(mtz_res.get("cost", 0.0), 6), + "time_s": round(mtz_time, 3), + "status": mtz_res.get("status", "?"), + "route": route_mtz, + } + print(f" path_cost={row['qap_mip']['path_cost']} time={mtz_time:.3f}s status={mtz_res.get('status','?')}") + except Exception as e: + row["qap_mip"] = {"status": "error", "detail": str(e)[:200]} + print(f" FAILED: {e}") + + # ── Cross-evaluation ── + print(f"\n ── Cross-evaluation ──") + + def _nearest_neighbor_tour(nodes: list[int], cost_mat: np.ndarray) -> tuple: + if len(nodes) <= 1: + return nodes, 0.0 + unvisited = set(nodes) + current = nodes[0] + unvisited.remove(current) + tour = [current] + total = 0.0 + while unvisited: + best_n = -1 + best_c = float("inf") + for nxt in unvisited: + if cost_mat[current, nxt] < best_c: + best_c = cost_mat[current, nxt] + best_n = nxt + tour.append(best_n) + unvisited.remove(best_n) + total += best_c + current = best_n + return tour, total + + for key in ("qubo_sa", "qubo_mip", "qubo_mip_card"): + entry = row.get(key, {}) + if entry.get("status") in ("ok", "feasible", "optimal"): + sol = entry.get("x", []) + selected_nodes = [i for i, v in enumerate(sol) if v] + if len(selected_nodes) >= 2: + tour, tc = _nearest_neighbor_tour(selected_nodes, raw) + entry["nn_tour"] = tour + entry["nn_tour_asymmetric_cost"] = round(tc, 6) + + # Also solve exact TSP on just these K nodes using the asymmetric cost + K = len(selected_nodes) + idx_map = {j: i for i, j in enumerate(selected_nodes)} + sub_raw = np.array([[raw[i, j] for j in selected_nodes] for i in selected_nodes]) + sub_tsp = solve_tsp_mtz(sub_raw, time_limit=min(30.0, tl)) + if "error" not in sub_tsp and sub_tsp.get("route"): + remapped_route = [selected_nodes[i] for i in sub_tsp["route"]] + entry["tsp_on_subset_route"] = remapped_route + entry["tsp_on_subset_cost"] = round(sub_tsp["cost"], 6) + print(f" {key} NN-tour={tc:.4f} TSP-on-subset={sub_tsp['cost']:.4f} (K={K})") + else: + print(f" {key} NN-tour={tc:.4f} (K={K}, TSP skipped)") + + # Evaluate TSP tour's subset energy + qap_entry = row.get("qap_mip", {}) + if qap_entry.get("route"): + x_tsp = [1] * n + qap_entry["as_subset_energy"] = round(symmetric_subset_cost(x_tsp, sym), 6) + qap_entry["as_subset_asymmetric_cost"] = round(asymmetric_subset_cost(x_tsp, raw), 6) + print(f" qap_mip (TSP all {n}) path_cost={qap_entry['path_cost']}") + + print(f"\n ── Summary ──") + row["max_anisotropy"] = max_aniso + + results[f"n_{n}"] = row + + return results + + +if __name__ == "__main__": + import argparse + parser = argparse.ArgumentParser(description="Finsler-Randers QAP vs QUBO benchmark") + parser.add_argument("--sizes", type=int, nargs="+", default=[8, 12, 24]) + parser.add_argument("--time-limit", type=float, default=None, + help="Solver time limit (default: per-size limits)") + parser.add_argument("--seed", type=int, default=42) + parser.add_argument("--output", "-o", help="JSON output path") + args = parser.parse_args() + + if args.time_limit is not None: + tl = {s: args.time_limit for s in args.sizes} + else: + tl = None + + results = run_benchmark(sizes=args.sizes, time_limits=tl, 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(f"\n{'-'*60}") + print(output) diff --git a/AGENTS.md b/AGENTS.md index 6e92e762..da443596 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -153,7 +153,7 @@ See respective repositories for components. Shared utilities have been duplicate ## Core Surfaces -- Lean/Semantics: `0-Core-Formalism/lean/Semantics/` — Compiler surface includes `Semantics.SieveLemmas` and `Semantics.InteractionGraphSidon` (commit `e61bb627`, 3314 jobs, 0 errors). +- Lean/Semantics: `0-Core-Formalism/lean/Semantics/` — Compiler surface includes `Semantics.SieveLemmas` and `Semantics.InteractionGraphSidon` (commit `e61bb627`, 3314 jobs, 0 errors). New exploration module `Semantics.CompleteInteractionGraph` (complete directed graph / every-point-touches-every-point) builds standalone with one bounded `walkMatrix_off_diag` sorry. - Infrastructure shims and probes: `4-Infrastructure/shim/` - Hardware bring-up: `4-Infrastructure/hardware/` - Documentation and wiki surfaces: `6-Documentation/` @@ -605,6 +605,37 @@ A ContextStream node was written (`e967f515-3af9-46c9-9fc8-e5c766a6c4fc`, type=f - ContextStream node query: `search(mode="keyword", query="meta-solid")` - **Granular-superconductor analogy tested and negative (2026-06-19):** The hypothesis that a universal reduced field H*/Hc₂ ≈ 1/7 exists across granular superconductors was tested via Consensus search. No paper reports such a universal ratio; H* is always microstructure-dependent, varies by orders of magnitude, and is never normalized to bulk Hc₂. The 1/7 threshold remains grounded in Sidon doubling combinatorics and hard-sphere polydispersity only. See CITATION.cff for the 25 vortex-glass/granular references reviewed. +## Finsler-Randers QAP Benchmark (2026-06-21) + +`benchmark_finsler_qap.py` formulates Finsler-Randers directed routing as QAP (Quadratic Assignment +Problem via TSP-MTZ) and benchmarks against the QUBO subset-selection approach at n=8,12,24,48. + +### Results + +| n | β-aniso | QUBO-SA (sym) | QUBO-card (K=n/2) | QAP-LP relaxed | QAP-MIP (path) | QUBO-card NN | TSP-on-subset | +|---|---------|---------------|-------------------|----------------|----------------|-------------|---------------| +| 8 | 0.650 | 0.0 | 20.95 (K=4) | 21.10 | **15.17** (0.07s) | 3.96 | 6.47 | +| 12| 1.041 | 0.0 | 55.18 (K=6) | 27.30 | **19.45** (0.30s) | 5.11 | 8.31 | +| 24| 1.127 | 0.0 | 289.21 (K=12) | 38.78 | **24.53** (2.54s) | 10.68 | 12.54 | +| 48| 1.127 | 0.0 | 1294.73 (K=24) | 91.57 | **37.90** (13.08s)| 21.57 | 20.85 | + +### Key Findings + +1. **QAP-MIP (MTZ TSP) scales well** — n=48 solves to feasibility in 13s; n=24 in 2.5s +2. **QAP-LP is a loose bound** — always 2-3x higher than QAP-MIP path cost (assignment relaxation ignores subtour elimination) +3. **QUBO is structurally degenerate for all-positive Q_ij** — unconstrained QUBO always selects 0 nodes; cardinality constraint (K=n/2) gives meaningful subsets but the symmetric energy is fundamentally different from directed path cost +4. **2-phase strategy viable**: select K via QUBO-card, then route via TSP-on-subset. The combined cost at K=n/2 is roughly half the full TSP cost at n +5. **NN-tour vs TSP-on-subset**: greedy nearest-neighbor is within 3% of optimal TSP at n=48 but 30-60% off at n=8,12 +6. **Solver times**: QAP-MIP grows as O(n³) empirically (0.07s n=8 → 0.30s n=12 → 2.54s n=24 → 13.08s n=48) + +### Relevant Files + +- `4-Infrastructure/shim/benchmark_finsler_qap.py` — benchmark harness with MTZ-TSP, QUBO-MIP-card, QAP-LP, and cross-evaluation +- `4-Infrastructure/shim/qubo_highs.py` — QUBO→MIP bridge and `solve_route_lp` (TSP assignment relaxation) +- `4-Infrastructure/shim/qaoa_adapter.py` — `FinslerMetric`, `finsler_metric_to_qubo`, `geodesic_assignment` +- `/tmp/finsler_benchmark_v3.json` — n=8,12,24 results (structued JSON) +- `/tmp/finsler_benchmark_n48.json` — n=48 results + ### When to Use ContextStream Search: ✅ Project is indexed and fresh