mirror of
https://github.com/allaunthefox/Research-Stack.git
synced 2026-07-31 03:05:21 +00:00
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)
This commit is contained in:
parent
7b498b95e4
commit
5a4a46d486
2 changed files with 630 additions and 1 deletions
598
4-Infrastructure/shim/benchmark_finsler_qap.py
Normal file
598
4-Infrastructure/shim/benchmark_finsler_qap.py
Normal file
|
|
@ -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)
|
||||
33
AGENTS.md
33
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
|
||||
|
||||
<!-- BEGIN ContextStream -->
|
||||
### When to Use ContextStream Search:
|
||||
✅ Project is indexed and fresh
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue