mirror of
https://github.com/allaunthefox/Research-Stack.git
synced 2026-07-31 03:05:21 +00:00
allaun
5a4a46d486
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)
598 lines
23 KiB
Python
598 lines
23 KiB
Python
#!/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)
|