"""图论/动态规划工具箱:开放路径排序、下界与路由度量(S1/S5)。 设计原则(docs/第三四问图论与动态规划优化方案.md §1、§3.2、§3.4): - **开放路径**(起点固定为原点、终点自由)是本题的真实图问题:一次作业回合 = 依次停靠若干作业点,代价 = Σ 移动时间 + Σ 动作耗时; - 小规模(n ≤ ``EXACT_LIMIT``)直接用 Held-Karp 精确 DP,规模更大时用 最近邻打底 + 2-opt / Or-opt 改进,并始终给出 **MST 下界**用于在线后悔值审计; - 所有函数只做"排序与度量",不承担任何安全责任,也不改变覆盖/清除语义。 约定:边权为**欧氏距离**(米);时间折算由调用方按 ``speed_mps`` 处理。 2-opt / Or-opt 用 O(1) 增量代价评估,使 n ≈ 10² 量级仍可在毫秒内完成。 """ from __future__ import annotations import math from .config import Vec2 EXACT_LIMIT = 13 # Held-Karp 状态数 2^n·n:n=13 → 106 496,纯 Python 可接受 OR_OPT_LIMIT = 60 # 作业点数超过此值时跳过高阶 Or-opt(见 improve_open_path) def path_length(order: list[Vec2], start: Vec2 = (0.0, 0.0)) -> float: """开放路径长度(起点 start,终点自由)。""" total = 0.0 position = start for point in order: total += math.dist(position, point) position = point return total def nearest_neighbor_order(points: list[Vec2], start: Vec2 = (0.0, 0.0)) -> list[Vec2]: """最近邻开放路径(O(n²),启发式上界)。""" remaining = list(points) order: list[Vec2] = [] position = start while remaining: nxt = min(remaining, key=lambda q: math.dist(position, q)) remaining.remove(nxt) order.append(nxt) position = nxt return order def two_opt_open(order: list[Vec2], start: Vec2 = (0.0, 0.0), max_passes: int = 8) -> list[Vec2]: """开放路径 2-opt(终点自由:只反转内部片段,不做首尾相接)。 增量代价 O(1):反转 [i..j] 只影响 (prev, i) 与 (j, next) 两条边。 """ best = list(order) n = len(best) if n < 3: return best for _ in range(max_passes): improved = False for i in range(n - 1): prev = best[i - 1] if i > 0 else start pi = best[i] d_prev_i = math.dist(prev, pi) for j in range(i + 1, n): pj = best[j] nxt = best[j + 1] if j + 1 < n else None d_prev_j = math.dist(prev, pj) removed = d_prev_i added = d_prev_j if nxt is not None: removed += math.dist(pj, nxt) added += math.dist(pi, nxt) if added < removed - 1e-9: best[i:j + 1] = best[i:j + 1][::-1] improved = True pi = best[i] d_prev_i = math.dist(prev, pi) if not improved: break return best def or_opt_open(order: list[Vec2], start: Vec2 = (0.0, 0.0), max_passes: int = 4, max_segment: int = 3) -> list[Vec2]: """Or-opt:把长度 ≤ max_segment 的片段搬家(开放路径,增量代价 O(1))。""" best = list(order) n = len(best) if n < 4: return best for _ in range(max_passes): improved = False for length in range(1, min(max_segment, n - 1) + 1): for i in range(0, n - length + 1): segment = best[i:i + length] rest = best[:i] + best[i + length:] for j in range(len(rest) + 1): if j == i: continue candidate = rest[:j] + segment + rest[j:] if _segment_move_gain(best, i, length, j, start) > 1e-9: best = candidate improved = True break if improved: break if improved: break if not improved: break return best def _segment_move_gain(order: list[Vec2], i: int, length: int, j: int, start: Vec2) -> float: """把 order[i:i+length] 插到去掉该段后的第 j 位所省下的长度(正=更省)。""" rest = order[:i] + order[i + length:] segment = order[i:i + length] def length_of(seq: list[Vec2]) -> float: return path_length(seq, start) return length_of(order) - length_of(rest[:j] + segment + rest[j:]) def improve_open_path(points: list[Vec2], start: Vec2 = (0.0, 0.0), max_passes: int = 8) -> list[Vec2]: """最近邻 + 2-opt + Or-opt 的开放路径改进(大规模兜底)。 Or-opt 的搬移评估是 O(n³),作业点数超过 :data:`OR_OPT_LIMIT` 时只做 O(n²) 的 2-opt(度量口径下 Or-opt 的相对收益远小于其墙钟代价)。 """ order = nearest_neighbor_order(points, start) order = two_opt_open(order, start, max_passes) if len(order) > OR_OPT_LIMIT: return order return or_opt_open(order, start, max_passes) def held_karp_open(points: list[Vec2], start: Vec2 = (0.0, 0.0)) -> list[Vec2]: """Held-Karp 精确开放路径(起点固定、终点自由)。 状态 ``dp[mask][j]`` = 访问集合 mask、当前停在 j 的最小长度; 规模 2^n·n,仅用于 n ≤ :data:`EXACT_LIMIT`。返回最优访问顺序。 """ n = len(points) if n == 0: return [] if n == 1: return list(points) dist = [[math.dist(points[i], points[j]) for j in range(n)] for i in range(n)] from_start = [math.dist(start, points[j]) for j in range(n)] full = (1 << n) - 1 inf = float("inf") dp: list[list[float]] = [[inf] * n for _ in range(1 << n)] parent: list[list[int]] = [[-1] * n for _ in range(1 << n)] for j in range(n): dp[1 << j][j] = from_start[j] for mask in range(1 << n): row = dp[mask] for j in range(n): current = row[j] if current == inf or not (mask >> j) & 1: continue rest = full ^ mask while rest: bit = rest & -rest k = bit.bit_length() - 1 rest ^= bit candidate = current + dist[j][k] nxt = mask | bit if candidate < dp[nxt][k]: dp[nxt][k] = candidate parent[nxt][k] = j best_end = min(range(n), key=lambda j: dp[full][j]) order: list[Vec2] = [] mask, j = full, best_end while j != -1: order.append(points[j]) previous = parent[mask][j] mask ^= 1 << j j = previous order.reverse() return order def optimal_open_path(points: list[Vec2], start: Vec2 = (0.0, 0.0), exact_limit: int = EXACT_LIMIT) -> list[Vec2]: """离线最优(或近优)开放路径:小规模精确,大规模 2-opt + Or-opt。""" unique = _dedupe(points) if len(unique) <= exact_limit: return held_karp_open(unique, start) return improve_open_path(unique, start) def mst_weight(points: list[Vec2], start: Vec2 | None = (0.0, 0.0)) -> float: """最小生成树权重(Prim,O(n²));含 start 时给出开放路径最优值的下界。 开放路径删除其最后一条边即得一条生成树,故 MST 权重 ≤ 最优路径长度。 """ nodes = ([start] if start is not None else []) + list(points) n = len(nodes) if n <= 1: return 0.0 in_tree = [False] * n best = [float("inf")] * n best[0] = 0.0 total = 0.0 for _ in range(n): u, u_cost = -1, float("inf") for i in range(n): if not in_tree[i] and best[i] < u_cost: u, u_cost = i, best[i] if u < 0: break in_tree[u] = True total += u_cost ux, uy = nodes[u] for v in range(n): if in_tree[v]: continue distance = math.hypot(nodes[v][0] - ux, nodes[v][1] - uy) if distance < best[v]: best[v] = distance return total def _dedupe(points, cell_m: float = 0.5) -> list[Vec2]: """按 ``cell_m`` 网格去重(同一物理点只算一次作业)。""" seen: set[tuple[int, int]] = set() unique: list[Vec2] = [] for point in points: key = (int(round(point[0] / cell_m)), int(round(point[1] / cell_m))) if key in seen: continue seen.add(key) unique.append((float(point[0]), float(point[1]))) return unique def route_metrics(service_points, actual_move_m: float, speed_mps: float = 5.0) -> dict: """路由度量:离线最优开放回路、MST 下界与相对 gap(S1 必报列)。""" unique = _dedupe(service_points) offline = path_length(optimal_open_path(unique)) if unique else 0.0 lower = mst_weight(unique) gap = ((actual_move_m - offline) / actual_move_m) if actual_move_m > 0 else 0.0 return { "distinct_service_points": len(unique), "offline_opt_open_path_m": round(offline, 1), "offline_opt_open_path_s": round(offline / speed_mps, 1), "route_lower_bound_s": round(lower / speed_mps, 1), "route_gap_vs_offline": round(gap, 4), }