Files
mcm-26-b/docs/paper/data/source-listings.tex
T
2026-09-13 11:45:01 +08:00

1463 lines
72 KiB
TeX
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
% Generated by scripts/build_code_appendix.py from mcm26b-核心代码.ipynb.
% Do not edit by hand: rerun the generator after the notebook changes.
\begingroup\linespread{1.0}\selectfont
\subsection*{动作计时与开放路径排序}
三段耗时按整数微秒累计;后备清除的走格顺序用最近邻打底并做 2-opt 改进。
\noindent\small \texttt{CostBreakdown}:一次合法动作的三段耗时:移动 + 换频道(仅检测)+ 动作本身。
\begin{lstlisting}[style=corecode]
class CostBreakdown:
"""一次合法动作的三段耗时 [题面]:移动 + 换频道(仅检测)+ 动作本身。"""
move_s: float
switch_s: float
action_s: float
@property
def total_s(self) -> float:
return self.total_us / 1_000_000
@property
def total_us(self) -> int:
"""逐段量化后按整数微秒累计,与官方虚拟时钟一致。"""
return sum(round(value * 1_000_000)
for value in (self.move_s, self.switch_s, self.action_s))
\end{lstlisting}
\noindent\small \texttt{measure\_cost}:合法 /measure 的耗时;移动按直线距离,换频道只看前后频道是否不同。
\begin{lstlisting}[style=corecode]
def measure_cost(prev_pos: Vec2, pos: Vec2, prev_channel: int, channel: int,
cost: CostConfig = CostConfig()) -> CostBreakdown:
"""合法 /measure 的耗时。移动按直线距离;换频道只看前后频道是否不同。"""
return CostBreakdown(
move_s=dist(prev_pos, pos) / cost.speed_mps,
switch_s=0.0 if channel == prev_channel else cost.switch_s,
action_s=cost.measure_s,
)
\end{lstlisting}
\noindent\small \texttt{clear\_cost}:合法 /clear 的耗时;清除不切换测向机频道,故不含换频项。
\begin{lstlisting}[style=corecode]
def clear_cost(prev_pos: Vec2, pos: Vec2, success: bool,
cost: CostConfig = CostConfig()) -> CostBreakdown:
"""合法 /clear 的耗时。"""
return CostBreakdown(
move_s=dist(prev_pos, pos) / cost.speed_mps,
switch_s=0.0,
action_s=cost.clear_success_s if success else cost.clear_fail_s,
)
\end{lstlisting}
\noindent\small \texttt{move\_time}:任意两点间的移动耗时,供规划使用。
\begin{lstlisting}[style=corecode]
def move_time(a: Vec2, b: Vec2, cost: CostConfig = CostConfig()) -> float:
"""任意两点间的移动耗时(规划用)。"""
return dist(a, b) / cost.speed_mps
\end{lstlisting}
\noindent\small \texttt{nearest\_neighbor\_order}:最近邻开放路径,作为路径排序的启发式上界。
\begin{lstlisting}[style=corecode]
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
\end{lstlisting}
\noindent\small \texttt{two\_opt\_open}:开放路径 2-opt 改进;终点自由,只反转内部片段。
\begin{lstlisting}[style=corecode]
def two_opt_open(order: list[Vec2], start: Vec2 = (0.0, 0.0),
max_passes: int = 8) -> list[Vec2]:
"""开放路径 2-opt(终点自由:只反转内部片段,不做首尾相接)。"""
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
\end{lstlisting}
\subsection*{问题一:角域交会、直径与最小包围圆}
半平面交求定位区域,顶点对枚举直径,Welzl 算法求最小包围圆,并输出直径圆覆盖判定与清除半径判据。
\noindent\small \texttt{Measurement}:一次示向度观测:测点坐标、示向度与误差半角(题设为 $1^\circ$)。
\begin{lstlisting}[style=corecode]
class Measurement:
x: float
y: float
bearing_deg: float
half_angle_deg: float = 1.0
\end{lstlisting}
\noindent\small \texttt{halfplanes}:一次观测的三个半平面约束:角域两侧边界与前向约束,法向量已单位化。
\begin{lstlisting}[style=corecode]
def halfplanes(m: Measurement) -> list[tuple[float, float, float]]:
"""Return a*x+b*y<=c with unit normals, including the forward…"""
theta = m.bearing_deg % 360
lower, upper = radians(theta - m.half_angle_deg), radians(theta + m.half_angle_deg)
middle = radians(theta)
normals = [(sin(lower), -cos(lower)), (-sin(upper), cos(upper)),
(-cos(middle), -sin(middle))]
return [(a, b, a * m.x + b * m.y) for a, b in normals]
\end{lstlisting}
\noindent\small \texttt{convex\_hull}:单调链凸包(逆时针、无重复末点),输入坐标先平移缩放到同一量级。
\begin{lstlisting}[style=corecode]
def convex_hull(points: list[Point], eps: float = EPS) -> list[Point]:
"""Monotone chain, CCW, no repeated final vertex; retain segm…"""
ordered = sorted(points)
unique = []
for p in ordered:
if not unique or distance(p, unique[-1]) > eps:
unique.append(p)
if len(unique) <= 2:
return unique
def chain(seq):
out = []
for p in seq:
while len(out) >= 2:
a = (out[-1][0] - out[-2][0], out[-1][1] - out[-2][1])
b = (p[0] - out[-1][0], p[1] - out[-1][1])
# Distance-to-line tolerance, not an unscaled area tolerance.
if cross(a, b) > eps * max(hypot(*a), hypot(*b), eps):
break
out.pop()
out.append(p)
return out
return chain(unique)[:-1] + chain(reversed(unique))[:-1]
\end{lstlisting}
\noindent\small \texttt{diameter}:顶点对枚举区域直径;有界凸集的最远点必为顶点。
\begin{lstlisting}[style=corecode]
def diameter(points: list[Point]):
if not points:
raise ValueError("空点集没有本程序定义下的直径")
a = b = points[0]
value = 0.0
for p, q in combinations(points, 2):
d = distance(p, q)
if d > value:
value, a, b = d, p, q
return value, (a, b)
\end{lstlisting}
\noindent\small \texttt{smallest\_enclosing\_circle}:Welzl 随机增量最小包围圆;固定随机种子保证结果可复现。
\begin{lstlisting}[style=corecode]
def smallest_enclosing_circle(points: list[Point], eps: float = EPS):
"""Welzl's randomized incremental algorithm, using an explici…"""
if not points:
raise ValueError("空集不计算最小包围圆")
if not isfinite(eps) or eps <= 0:
raise ValueError("eps必须为有限正数")
if any(len(p) != 2 or any(isinstance(x, bool) or not isinstance(x, (int, float))
or not isfinite(x) for x in p) for p in points):
raise ValueError("点集必须由有限二维坐标组成")
origin = points[0]
scale = max(distance(origin, p) for p in points)
if scale == 0:
return tuple(origin), 0.0
local = [((p[0]-origin[0])/scale, (p[1]-origin[1])/scale) for p in points]
tolerance = min(eps/scale, EPS)
Random(20260910).shuffle(local)
# Frame = (remaining prefix length, forced boundary, returning from child).
stack = [(len(local), (), False)]
circle = None
while stack:
n, boundary, returning = stack.pop()
if n == 0 or len(boundary) == 3:
# Forced support points must lie ON the boundary. For three points
# use their circumcircle, even when their unconstrained MEC is smaller.
circle = _circumcircle(*boundary) if len(boundary) == 3 else None
if circle is None and boundary:
circle = _enclosing_circle_reference(list(boundary), tolerance)
elif not returning:
stack.append((n, boundary, True))
stack.append((n-1, boundary, False))
elif not _circle_covers(circle, local[n-1], tolerance):
stack.append((n-1, boundary + (local[n-1],), False))
assert circle is not None
center = (origin[0]+circle[0][0]*scale, origin[1]+circle[0][1]*scale)
# Recompute after transforming back: returned circle contains all input points
# to distance() precision, including any roundoff from the large translation.
return center, max(distance(center, p) for p in points)
\end{lstlisting}
\noindent\small \texttt{solve}:问题一主入口:半平面交求定位区域,区分空集、无界与退化情形,再给出直径、直径圆覆盖判定与最小包围圆。
\begin{lstlisting}[style=corecode]
def solve(measurements: list[Measurement], eps: float = EPS) -> dict:
"""Solve Q1's pure angular intersection (no disk/range prior…"""
if not isfinite(eps) or eps <= 0:
raise ValueError("eps必须为有限正数")
if not measurements:
return {"status": "unbounded", "vertices": [], "diameter_m": None,
"diameter_infinite": True, "recession_direction": [1.0, 0.0],
"message": "没有观测,整个平面均可行"}
if len(measurements) > 100:
raise ValueError("参考实现最多接收100个检测点;更大规模应使用半平面交算法")
ox = sum(m.x for m in measurements) / len(measurements)
oy = sum(m.y for m in measurements) / len(measurements)
scale = max(1.0, max(hypot(m.x - ox, m.y - oy) for m in measurements))
local = [Measurement((m.x - ox) / scale, (m.y - oy) / scale,
m.bearing_deg, m.half_angle_deg) for m in measurements]
planes = [h for m in local for h in halfplanes(m)]
def inside(p):
tol = eps * max(1.0, abs(p[0]), abs(p[1]))
return all(a * p[0] + b * p[1] <= c + tol for a, b, c in planes)
feasible_vertices = []
for (a, b, c), (d, e, f) in combinations(planes, 2):
det = a * e - b * d
if det == 0.0:
continue # 精确平行:无交点(近平行对的交点是合法远端顶点,保留)
p = ((c * e - b * f) / det, (a * f - c * d) / det)
if isfinite(p[0]) and isfinite(p[1]) and inside(p):
feasible_vertices.append(p)
if not feasible_vertices:
return {"status": "empty", "vertices": [], "diameter_m": None,
"diameter_infinite": False, "tolerance_m": eps * scale,
"message": "测向角域没有公共点;检查频道、数据及误差界"}
# A nontrivial 2D recession cone has a ray tangent to an input halfplane.
# Use a much tighter angular threshold than the positional feasibility epsilon.
ray = None
for a, b, _ in planes:
for direction in ((b, -a), (-b, a)):
if all(d * direction[0] + e * direction[1] <= 2e-14
for d, e, _ in planes):
ray = direction
break
if ray is not None:
break
hull = convex_hull(feasible_vertices, eps)
def world(p):
return [ox + p[0] * scale, oy + p[1] * scale]
result = {"status": "unbounded" if ray else
("point" if len(hull) == 1 else "segment" if len(hull) == 2 else "polygon"),
"vertices": [world(p) for p in hull], "tolerance_m": eps * scale,
"diameter_infinite": ray is not None}
if ray is not None:
result.update(diameter_m=None, recession_direction=list(ray),
message="存在非零衰退方向,区域无界,直径为无穷大")
return result
d, pair = diameter(hull)
midpoint = ((pair[0][0] + pair[1][0]) / 2, (pair[0][1] + pair[1][1]) / 2)
distances = [distance(midpoint, p) for p in hull]
witness_index = max(range(len(hull)), key=distances.__getitem__)
max_distance = distances[witness_index]
center, radius = smallest_enclosing_circle(hull, eps)
area = abs(sum(cross(hull[i], hull[(i + 1) % len(hull)])
for i in range(len(hull)))) / 2
result.update(diameter_m=d * scale, diameter_pair=[world(p) for p in pair],
area_m2=area * scale ** 2,
diameter_circle={"center": world(midpoint), "radius_m": d * scale / 2,
"max_vertex_distance_m": max_distance * scale,
"covers": max_distance <= d / 2 + eps * 4,
"coverage_tolerance_m": 4 * eps * scale,
"max_excess_m": (max_distance-d/2) * scale,
"outside_vertex_indices": [i for i, value in enumerate(distances)
if value > d/2 + 4*eps],
"witness_vertex": world(hull[witness_index])
if max_distance > d/2 + 4*eps else None},
minimum_enclosing_circle={"center": world(center), "radius_m": radius * scale},
guaranteed_clearable_at_mec_center=radius * scale <= 20 - eps * scale * 4,
clearance_certificate={"model": "angular_intersection_only",
"clear_radius_m": 20.0,
"numerical_margin_m": 4*eps*scale,
"radius_with_margin_m": (radius+4*eps)*scale,
"remaining_margin_m": 20-(radius+4*eps)*scale,
"center": world(center)})
return result
\end{lstlisting}
\subsection*{问题二:候选集与保证接收的第二检测点}
候选集用圆盘外切多边形与角域半平面裁剪表示;安全盘半径取$1000-\rho$,盘内候选按最坏后验包围半径、同分按移动代价取优。
\noindent\small \texttt{CandidateRegion}:某频道候选集合的保守凸外包多边形(顶点按 CCW;允许退化为 1-2 点)。
\begin{lstlisting}[style=corecode]
class CandidateRegion:
"""某频道候选集合的保守凸外包多边形(顶点按 CCW;允许退化为 1-2 点)。"""
@classmethod
def disk(cls, center: Vec2, radius_m: float, sides: int,
slack_m: float | None = None) -> "CandidateRegion":
if slack_m is None:
slack_m = radius_m * (1 - math.cos(math.pi / sides))
return cls(circumscribed_polygon(center, radius_m, sides), slack_m)
@classmethod
def initial(cls, station: Vec2, config: ProblemConfig) -> "CandidateRegion":
"""首个 direction 后的基础集:Ω ∩ B(station, r_max) 的保守外包。"""
field_cfg, sides = config.field, config.search.polygon_sides
slack = field_cfg.omega_radius_m * (1 - math.cos(math.pi / sides))
region = cls.disk((0.0, 0.0), field_cfg.omega_radius_m, sides, slack)
region.clip_disk(station, field_cfg.r_eff_max_m, sides)
return region
def clip_disk(self, center: Vec2, radius_m: float, sides: int) -> None:
"""与圆盘 B(center, radius) 的外切多边形求交(近似误差入账)。"""
for a, b, c in polygon_halfplanes(circumscribed_polygon(center, radius_m, sides)):
self._clip(a, b, c)
self.slack_m += radius_m * (1 - math.cos(math.pi / sides))
def clip_wedge(self, station: Vec2, bearing_deg: float,
half_angle_deg: float) -> None:
"""与测向角域 W(station, bearing±half_angle) 求交(半角须 <90°)。
两个边界半平面与 localization.halfplanes 同口径,0°/360° 跨界由
sin/cos 直接处理,无需分支。
"""
if not 0 < half_angle_deg < 90:
raise ValueError("角域半角须在 (0, 90) 度")
lower = math.radians(bearing_deg - half_angle_deg)
upper = math.radians(bearing_deg + half_angle_deg)
sx, sy = float(station[0]), float(station[1])
self._clip(math.sin(lower), -math.cos(lower),
math.sin(lower) * sx - math.cos(lower) * sy)
self._clip(-math.sin(upper), math.cos(upper),
-math.sin(upper) * sx + math.cos(upper) * sy)
def mec(self) -> tuple[Vec2, float]:
"""最小包围圆(圆心, 半径)。半径是真实候选集 ρ 的上界。"""
if self._mec is None:
if not self.polygon:
raise InconsistentBelief("候选区域为空")
self._mec = smallest_enclosing_circle(self.polygon)
return self._mec
def witnesses(self, spacing_m: float) -> list[Vec2]:
"""顶点 + 边界按间距细分的采样点(最坏后验评估用)。"""
if len(self.polygon) < 3:
return list(self.polygon)
points = list(self.polygon)
count = len(self.polygon)
for index in range(count):
p, q = self.polygon[index], self.polygon[(index + 1) % count]
length = math.dist(p, q)
steps = int(length // spacing_m)
for k in range(1, steps + 1):
ratio = k / (steps + 1)
points.append((p[0] + ratio * (q[0] - p[0]),
p[1] + ratio * (q[1] - p[1])))
return _dedupe(points, tol=1e-9)
def contains(self, point: Vec2, tol_m: float = 0.05) -> bool:
"""点是否在候选集内(带容差;用于真源保留性检验,不做删点依据)。"""
if not self.polygon:
return False
if len(self.polygon) == 1:
return math.dist(point, self.polygon[0]) <= tol_m + self.slack_m
if len(self.polygon) == 2:
p, q = self.polygon
return _point_segment_distance(point, p, q) <= tol_m + self.slack_m
for index, p in enumerate(self.polygon):
q = self.polygon[(index + 1) % len(self.polygon)]
if (q[0] - p[0]) * (point[1] - p[1]) - (q[1] - p[1]) * (point[0] - p[0]) \
< -tol_m - self.slack_m:
return False
return True
\end{lstlisting}
\noindent\small \texttt{ClearCertificate}:清除充分条件:候选集外包半径与近似误差之和不超过清除半径。
\begin{lstlisting}[style=corecode]
class ClearCertificate:
"""清除证书:候选集外包半径(含近似误差)不超过 20m 阈值。"""
channel: int
center: Vec2
radius_m: float
limit_m: float
certified: bool
\end{lstlisting}
\noindent\small \texttt{clear\_certificate}:由候选集最小包围圆判定是否可以就地清除。
\begin{lstlisting}[style=corecode]
def clear_certificate(region: CandidateRegion, channel: int,
config: ProblemConfig = ProblemConfig()) -> ClearCertificate:
"""ρ_poly <= 20 - 外包近似误差 - 工程余量 才允许发起清除。"""
center, radius = region.mec()
limit = (config.field.clear_radius_m - region.slack_m
- config.search.clear_margin_m)
return ClearCertificate(channel, center, radius, limit, radius <= limit)
\end{lstlisting}
\noindent\small \texttt{candidate\_points}:保证接收盘内的候选点:盘心加各比例半径上的均匀方向点。
\begin{lstlisting}[style=corecode]
def candidate_points(center: Vec2, safe_radius_m: float, rings, angles: int) -> list[Vec2]:
"""保证接收盘内的候选点:盘心 + 各比例半径 × 均匀方向。"""
points = [tuple(center)]
for fraction in rings:
radius = fraction * safe_radius_m
if radius <= 1e-9:
continue
step = 2 * math.pi / angles
for k in range(angles):
angle = step * k + step / 2
points.append((center[0] + radius * math.cos(angle),
center[1] + radius * math.sin(angle)))
return points
\end{lstlisting}
\noindent\small \texttt{safe\_disk\_radius}:保证接收盘半径 $1000-\rho$;小于零表示没有保证接收的区域。
\begin{lstlisting}[style=corecode]
def safe_disk_radius(region: CandidateRegion,
config: ProblemConfig = ProblemConfig()) -> float:
"""保证接收盘 B(mec, safe) 的半径;<0 表示无保证(ρ 超过最小接收半径)。"""
return config.field.r_eff_min_m - region.mec()[1] - 1e-3
\end{lstlisting}
\noindent\small \texttt{select\_measurement\_point}:在保证接收盘内选下一检测点:主指标为最坏后验包围半径,同分取移动代价较小者。
\begin{lstlisting}[style=corecode]
def select_measurement_point(region: CandidateRegion, from_position: Vec2,
config: ProblemConfig = ProblemConfig()) -> SecondPointPlan:
"""在保证接收盘内选择下一检测点(P3 核心;对任意候选集通用)。"""
field, search = config.field, config.search
center, radius = region.mec()
half_angle = field.bearing_error_half_deg + search.bearing_rounding_margin_deg
def move_cost(point: Vec2) -> float:
return move_time(from_position, point, config.cost)
# 1mm 数值余量在 safe_disk_radius 内记账。
safe_radius = safe_disk_radius(region, config)
if safe_radius < 0:
return SecondPointPlan(tuple(center), radius, move_cost(center),
"无保证接收盘(ρ>1000):探索 MEC 圆心,可能 no_signal")
coarse = region.witnesses(search.witness_spacing_m)
samples = candidate_points(center, safe_radius,
search.candidate_rings, search.candidate_angles)
near_radius, sides = field.near_radius_m, search.polygon_sides
best_key, best_point, best_worst = None, None, None
for point in samples:
worst = _worst_posterior_radius(region, point, coarse, half_angle,
near_radius, sides)
key = (round(worst, 6), round(move_cost(point), 6))
if best_key is None or key < best_key:
best_key, best_point, best_worst = key, point, worst
dense = _worst_posterior_radius(region, best_point,
region.witnesses(search.witness_spacing_m / 4),
half_angle, near_radius, sides)
reason = (f"安全盘 B(centre,{safe_radius:.0f}m) 内 {len(samples)} 候选;"
f"粗采样最坏后验 {best_worst:.1f}m,1/4 间距加密复核 {dense:.1f}m")
return SecondPointPlan(tuple(best_point), dense, move_cost(best_point), reason)
\end{lstlisting}
\noindent\small \texttt{select\_second\_point}:问题二入口:由首次示向度构造候选集并给出第二检测点。
\begin{lstlisting}[style=corecode]
def select_second_point(client: SimulatorClient, first_position: Vec2,
first_bearing_deg: float,
config: ProblemConfig = ProblemConfig()) -> SecondPointPlan:
"""P3 入口:由首次 direction 构造 K1,返回保证接收的第二检测点方案。"""
region = CandidateRegion.initial(first_position, config)
half_angle = config.field.bearing_error_half_deg \
+ config.search.bearing_rounding_margin_deg
region.clip_wedge(first_position, first_bearing_deg, half_angle)
return select_measurement_point(region, client.state.position, config)
\end{lstlisting}
\subsection*{问题三:覆盖计划、检测记录与清除证书}
覆盖计划与运行期抽样校验、逐频道检测记录,以及七点解析覆盖与全向/定向两种覆盖语义的证书判据。
\noindent\small \texttt{CoverageCertificate}:覆盖证书:采样分辨率、位置与朝向目标数、缺口数与最坏命中距离。
\begin{lstlisting}[style=corecode]
class CoverageCertificate:
"""覆盖证书:采样语义 + 缺口 + 最坏命中距离。"""
kind: str # "full_disk" | "forward_half_disk"
target_step_m: float
headings: int
r_eff_min_m: float
targets: int
uncovered: int
worst_hit_distance_m: float
nn_tour_m: float
generator_sha256: str = ""
verify_heading_margin_deg: float = 0.0
@property
def valid(self) -> bool:
return self.uncovered == 0 and self.targets > 0
def as_dict(self) -> dict:
return {"kind": self.kind, "target_step_m": self.target_step_m,
"headings": self.headings, "r_eff_min_m": self.r_eff_min_m,
"targets": self.targets, "uncovered": self.uncovered,
"worst_hit_distance_m": round(self.worst_hit_distance_m, 3),
"nn_tour_m": round(self.nn_tour_m, 1),
"nn_tour_virtual_s": round(self.nn_tour_m / SPEED_MPS, 1),
"generator_sha256": self.generator_sha256,
"verify_heading_margin_deg": self.verify_heading_margin_deg}
\end{lstlisting}
\noindent\small \texttt{CoveragePlan}:覆盖计划:点集、覆盖语义与证书。
\begin{lstlisting}[style=corecode]
class CoveragePlan:
"""一份覆盖计划:点集 + 语义 + 证书(可为解析证明)。"""
name: str
problem: int
points: tuple[Vec2, ...]
kind: str # "analytic" | "certified"
coverage: str # "full_disk" | "forward_half_disk"
note: str = ""
certificate: CoverageCertificate | None = None
generator: dict = field(default_factory=dict)
def as_dict(self) -> dict:
return {"name": self.name, "problem": self.problem, "kind": self.kind,
"coverage": self.coverage, "note": self.note,
"generator": self.generator,
"certificate": (self.certificate.as_dict()
if self.certificate is not None else None),
"points": [[p[0], p[1]] for p in self.points]}
\end{lstlisting}
\noindent\small \texttt{coverage\_points}:49 点覆盖网格,单元对角 $700\sqrt2\approx989.95\le1000$ m。
\begin{lstlisting}[style=corecode]
def coverage_points(config: ProblemConfig = ProblemConfig()) -> list[Vec2]:
"""49 点覆盖网格生成(纯几何,可独立测试)。"""
step = config.search.grid_step_m
half = config.search.grid_half_cells
return [(step * i, step * j)
for j in range(half, -half - 1, -1)
for i in range(-half, half + 1)]
\end{lstlisting}
\noindent\small \texttt{CoverageLedger}:逐频道检测记录,按(频道, 网格点)二维记账。
\begin{lstlisting}[style=corecode]
class CoverageLedger:
"""逐频道覆盖账本:(频道, 网格点) -> 是否已测。"""
config: ProblemConfig = field(default_factory=ProblemConfig)
points: list[Vec2] = field(default_factory=coverage_points)
measured: set = field(default_factory=set) # {(channel, point_index)}
skipped: dict = field(default_factory=dict) # channel -> 原因
def mark_measured(self, channel: int, point_index: int) -> None:
self.measured.add((channel, point_index))
def measured_at(self, channel: int, point_index: int) -> bool:
return (channel, point_index) in self.measured
def skip_channel(self, channel: int, reason: str) -> None:
"""记录整频道跳过原因(如 already_cleared / located_pending_clear)。"""
self.skipped[channel] = reason
def pending_points(self, channel: int) -> list[int]:
"""该频道尚未测量的网格点下标。扫描调度策略未实现(P4)。"""
return [index for index in range(len(self.points))
if not self.measured_at(channel, index)]
def coverage_complete(self, channel: int) -> bool:
"""该频道覆盖是否完成(含显式跳过)。"""
return channel in self.skipped or not self.pending_points(channel)
\end{lstlisting}
\noindent\small \texttt{seven\_points}:原点与环上六个等分点的七站覆盖集。
\begin{lstlisting}[style=corecode]
def seven_points(ring_radius_m: float = RING_RADIUS_M) -> list[tuple[float, float]]:
"""原点 + 环上 6 等分点的七站覆盖集。"""
return [(0.0, 0.0)] + [(ring_radius_m * math.cos(k * math.pi / 3),
ring_radius_m * math.sin(k * math.pi / 3)) for k in range(6)]
\end{lstlisting}
\noindent\small \texttt{seven\_plan}:第三问解析覆盖计划;环半径 $1159.7$ m 时最坏覆盖距离不超过 $984.539$ m。
\begin{lstlisting}[style=corecode]
def seven_plan(ring_radius_m: float) -> CoveragePlan:
"""第三问解析计划(覆盖证明由 validate_seven_coverage 给出)。"""
return CoveragePlan(
name=f"seven-r{ring_radius_m:.1f}", problem=3, kind="analytic",
coverage="full_disk",
note=("解析证明 max(r/√3, sqrt(R²+r²−√3Rr)) < 1000;"
"r∈[1122.9, 1732] 时成立,巡检半径越小巡回越短"),
points=tuple(seven_points(ring_radius_m)))
\end{lstlisting}
\noindent\small \texttt{verify\_full\_disk}:第三问覆盖语义:圆域内每点都落在某个被访点的最小接收圆盘内。
\begin{lstlisting}[style=corecode]
def verify_full_disk(points: list[Vec2], step_m: float = VERIFY_STEP_M,
r_eff: float = R_EFF_MIN_M) -> CoverageCertificate:
"""第三问语义:Ω 内每点都在某个被访点的 r_eff 圆盘内。"""
points = list(points)
positions = omega_positions(step_m)
table = _bucket(points, r_eff)
missing, worst = 0, 0.0
for (px, py) in positions:
ci, cj = int(px // r_eff), int(py // r_eff)
best = None
for di in (-1, 0, 1):
for dj in (-1, 0, 1):
for c_index in table.get((ci + di, cj + dj), ()):
cx, cy = points[c_index]
distance = math.hypot(cx - px, cy - py)
if distance <= r_eff and (best is None or distance < best):
best = distance
if best is None:
missing += 1
else:
worst = max(worst, best)
return CoverageCertificate("full_disk", step_m, 1, r_eff, len(positions),
missing, worst, nn_tour_length(points))
\end{lstlisting}
\noindent\small \texttt{verify\_forward\_half\_disk}:第四问覆盖语义:任意位置与朝向都存在前向半盘内的被访点。
\begin{lstlisting}[style=corecode]
def verify_forward_half_disk(points: list[Vec2], step_m: float = VERIFY_STEP_M,
headings: int = VERIFY_HEADINGS,
r_eff: float = R_EFF_MIN_M,
heading_margin_deg: float = 0.0) -> CoverageCertificate:
"""第四问语义:∀(p,h) 存在前向半盘内的被访点。"""
points = list(points)
targets = forward_targets(step_m, headings)
table = _bucket(points, r_eff)
missing, worst = 0, 0.0
sin_margin = math.sin(math.radians(max(heading_margin_deg, 0.0)))
for (px, py, heading) in targets:
ux, uy = math.cos(heading), math.sin(heading)
ci, cj = int(px // r_eff), int(py // r_eff)
best = None
for di in (-1, 0, 1):
for dj in (-1, 0, 1):
for c_index in table.get((ci + di, cj + dj), ()):
cx, cy = points[c_index]
dx, dy = cx - px, cy - py
distance_sq = dx * dx + dy * dy
if distance_sq > r_eff * r_eff:
continue
if dx * ux + dy * uy < sin_margin * math.sqrt(distance_sq) - 1e-9:
continue
distance = math.sqrt(distance_sq)
if best is None or distance < best:
best = distance
if best is None:
missing += 1
else:
worst = max(worst, best)
return CoverageCertificate("forward_half_disk", step_m, headings, r_eff,
len(targets), missing, worst,
nn_tour_length(points))
\end{lstlisting}
\noindent\small \texttt{plan\_instance}:按名称取覆盖计划;证书计划在运行期按同一分辨率复核。
\begin{lstlisting}[style=corecode]
def plan_instance(name: str, verify: bool = True) -> CoveragePlan:
"""统一入口:按名称取覆盖计划(内置解析计划 / 内嵌证书计划)。"""
if name == "square49":
return square_49_plan()
if name == "seven":
return seven_plan(1200.0)
if name not in EMBEDDED_PLANS:
raise FileNotFoundError(f"未找到覆盖计划 {name!r};现有:{plan_names()}")
plan = embedded_plan(name)
if verify and plan.kind == "certified":
certificate = verify_plan(plan)
if not certificate.valid:
raise ValueError(
f"覆盖计划 {name!r} 运行期校验失败:"
f"{certificate.uncovered}/{certificate.targets} 个 (p,h) 未被覆盖")
plan = CoveragePlan(
name=plan.name, problem=plan.problem, kind=plan.kind,
coverage=plan.coverage, note=plan.note, generator=plan.generator,
certificate=CoverageCertificate(
kind=certificate.kind, target_step_m=certificate.target_step_m,
headings=certificate.headings, r_eff_min_m=certificate.r_eff_min_m,
targets=certificate.targets, uncovered=certificate.uncovered,
worst_hit_distance_m=certificate.worst_hit_distance_m,
nn_tour_m=certificate.nn_tour_m,
generator_sha256=(plan.certificate.generator_sha256
if plan.certificate else "")),
points=plan.points)
return plan
\end{lstlisting}
\subsection*{问题三:网格候选集}
逐频道的网格幸存集,按格心裕量 $\ell/\sqrt2$ 内缩后再删格,据此支持负观测更新、跳测规则与缺席证明。
\noindent\small \texttt{ChannelParticles}:单个频道的网格候选集:格心裕量 $\ell/\sqrt2$,所有删除判定内缩该裕量,故真源所在格永不被删。
\begin{lstlisting}[style=corecode]
class ChannelParticles:
"""单个频道的网格幸存集。"""
def cell_center(self, index: int) -> Vec2:
i, j = divmod(index, self.width)
return (self.x0 + (i + 0.5) * self.spacing_m,
self.x0 + (j + 0.5) * self.spacing_m)
def discard(self, indices: list[int]) -> None:
"""公开删除入口:联合信念把"朝向不可能"的格一并删除时使用。"""
self._discard(indices)
def _discard(self, indices: list[int]) -> None:
if not indices:
return
coarse_w = self._coarse_width()
for index in indices:
if index not in self.alive:
continue
self.alive.discard(index)
key = (index // self.width // _COARSE_FACTOR) * coarse_w \
+ (index % self.width) // _COARSE_FACTOR
count = self._coarse.get(key, 0) - 1
if count <= 0:
self._coarse.pop(key, None)
else:
self._coarse[key] = count
self.version += 1
self._mec = None
def _kill_disk(self, center: Vec2, radius_m: float) -> None:
"""删除格心距 center <= radius_m 的全部幸存格(整格判据见模块注释)。"""
if not self.alive:
return
reach = radius_m + _DIST_EPS_M
ilo = max(0, int(math.floor((center[0] - reach - self.x0) / self.spacing_m)))
ihi = min(self.width - 1,
int(math.ceil((center[0] + reach - self.x0) / self.spacing_m)))
jlo = max(0, int(math.floor((center[1] - reach - self.x0) / self.spacing_m)))
jhi = min(self.width - 1,
int(math.ceil((center[1] + reach - self.x0) / self.spacing_m)))
limit_sq = reach * reach
doomed = []
for i in range(ilo, ihi + 1):
dx = self.x0 + (i + 0.5) * self.spacing_m - center[0]
base = i * self.width
for j in range(jlo, jhi + 1):
index = base + j
if index not in self.alive:
continue
dy = self.x0 + (j + 0.5) * self.spacing_m - center[1]
if dx * dx + dy * dy <= limit_sq:
doomed.append(index)
self._discard(doomed)
def _kill_outside_disk(self, center: Vec2, radius_m: float) -> None:
"""删除格心距 center > radius_m 的全部幸存格(保留盘内)。"""
if not self.alive:
return
reach = radius_m + _DIST_EPS_M
doomed = [index for index in self.alive
if math.dist(self.cell_center(index), center) > reach]
self._discard(doomed)
def _kill_outside_wedge(self, station: Vec2, reading_deg: float) -> None:
"""删除角向偏离读数超出 δ_eff + 角向格裕量的幸存格。
角向格裕量 = asin(m / max(dist−m, 1)):格内任一点相对格心的方位
偏差不超过它(点距测点还可再近 m,故分母再内缩 m)。
"""
if not self.alive:
return
sx, sy = float(station[0]), float(station[1])
doomed = []
for index in self.alive:
gx, gy = self.cell_center(index)
dist = math.hypot(gx - sx, gy - sy)
if dist <= self.margin_m + 1.0:
continue # 格贴测点:格内方位可指向任意方向,不可删
angular_margin = math.degrees(
math.asin(min(0.9999, self.margin_m / (dist - self.margin_m))))
gap = angular_gap(math.degrees(math.atan2(gy - sy, gx - sx)), reading_deg)
if gap > self.half_angle_deg + angular_margin + _ANG_EPS_DEG:
doomed.append(index)
self._discard(doomed)
def apply_measure(self, station: Vec2, result: str,
svd_deg: float | None = None, *,
directional: bool = False) -> None:
"""按读数更新幸存集。真源保留性由 Voronoi 裕量保证(模块注释)。
``directional=True``(第四问):no_signal **不删空间候选**——定向源
的静默还可能来自朝向遮挡,单独删空间候选会误删真源;朝向维度由
:class:`~mcm26b.belief.joint.JointBelief` 维护。
"""
if result == NO_SIGNAL:
if directional:
return
# no_signal ⟹ dist > r_eff >= 1000:整格落入 B(s, 1000−m) 的删。
self._kill_disk(station, self.r_eff_min_m - self.margin_m)
elif result == NEAR:
self._kill_outside_disk(station, self.near_radius_m + self.margin_m)
elif result == DIRECTION:
if svd_deg is None:
raise ValueError("direction 读数必须携带 svd_deg")
# direction 隐含 dist <= r_max [题面]
self._kill_outside_disk(station, self.r_eff_max_m + self.margin_m)
self._kill_outside_wedge(station, svd_deg)
else:
raise ValueError(f"未知测向结果: {result!r}")
def alive_points(self, limit: int | None = None) -> list[Vec2]:
sorted_alive = sorted(self.alive)
if limit is not None and len(sorted_alive) > limit:
step = len(sorted_alive) / limit
indices = {sorted_alive[int(k * step) % len(sorted_alive)]
for k in range(limit)}
indices = sorted(indices)
else:
indices = sorted_alive
return [self.cell_center(index) for index in indices]
def nearest_alive_distance(self, point: Vec2, limit_m: float) -> float | None:
"""最近幸存格心距离;> limit_m 时返回 None(粗格环扫 + 下界剪枝)。
返回值语义:恰为最近幸存格心距离(若 <= limit_m);None 表示所有
幸存格心都超出 limit_m。粗格聚合半径已计入扫描范围,不漏判。
"""
if not self.alive:
return None
coarse_cell_m = self.spacing_m * _COARSE_FACTOR
coarse_reach = coarse_cell_m * math.sqrt(2) / 2
cw = self._coarse_width()
ci_base = int(math.floor((point[0] - self.x0) / coarse_cell_m))
cj_base = int(math.floor((point[1] - self.x0) / coarse_cell_m))
ring_max = int(math.ceil((limit_m + coarse_reach) / coarse_cell_m)) + 1
best: float | None = None
for ring in range(ring_max + 1):
lower_bound = max(0.0, ring - 1) * coarse_cell_m - coarse_reach
if best is not None and lower_bound >= best:
break # 更外环的所有格心都 >= 下界 >= 当前最优
for ci in range(max(0, ci_base - ring), min(cw - 1, ci_base + ring) + 1):
for cj in range(max(0, cj_base - ring), min(cw - 1, cj_base + ring) + 1):
if max(abs(ci - ci_base), abs(cj - cj_base)) != ring:
continue # 只扫当前环(Chebyshev 环)
if (ci * cw + cj) not in self._coarse:
continue
ccx = self.x0 + (ci + 0.5) * coarse_cell_m
ccy = self.x0 + (cj + 0.5) * coarse_cell_m
if math.hypot(ccx - point[0], ccy - point[1]) \
> limit_m + coarse_reach:
continue
i0, j0 = ci * _COARSE_FACTOR, cj * _COARSE_FACTOR
for i in range(i0, min(i0 + _COARSE_FACTOR, self.width)):
base = i * self.width
for j in range(j0, min(j0 + _COARSE_FACTOR, self.width)):
index = base + j
if index not in self.alive:
continue
dist = math.dist(self.cell_center(index), point)
if dist <= limit_m and (best is None or dist < best):
best = dist
return best
def guaranteed_no_signal(self, station: Vec2) -> bool:
"""真源与该测点距离必 > r_max(>= r_eff)⟹ 必 no_signal,测了也无信息。"""
return self.nearest_alive_distance(
station, self.r_eff_max_m + self.margin_m) is None
def mec(self) -> tuple[Vec2, float]:
"""幸存集最小包围圆(圆心, 半径);按 version 缓存。"""
if self._mec is None:
if not self.alive:
raise ValueError("幸存集为空(该频道已可证明缺席)")
points = [self.cell_center(index) for index in self.alive]
self._mec = smallest_enclosing_circle(points)
return self._mec
def subsample(self, size: int) -> list[Vec2]:
"""按索引序均匀抽取的子采样(索引序≈空间序,代表点空间铺开)。"""
return self.alive_points(limit=size)
def expected_survival(self, point: Vec2, subsample,
representatives) -> float:
"""E[log 幸存占比]:越小该测点的期望信息增益越大。
r_eff 不确定性解析积分:dist∈(1000,1500] 时按 r_eff~U(1000,1500)
取 p_ns=(dist−1000)/500 为 no_signal 分支概率。
- no_signal 分支:删 B(point, 1000−m),幸存 = 子采样中盘外比例;
- direction 分支:以代表点方位为期望读数,角域 ∩ B(point,1500+m) 内幸存;
- near 分支(dist <= 5+m):可行集收缩到 5m 盘,幸存占比趋近 0。
"""
if not self.alive or len(subsample) == 0 or len(representatives) == 0:
return 0.0
if _np is not None and isinstance(subsample, _np.ndarray) \
and isinstance(representatives, _np.ndarray):
return self._expected_survival_vector(point, subsample, representatives)
eps = 1.0 / len(subsample)
keep_ns = max(1.0 - self.fraction_within(subsample, point,
self.r_eff_min_m - self.margin_m),
eps)
log_keep_ns = math.log(keep_ns)
total = 0.0
for rep in representatives:
dist = math.hypot(rep[0] - point[0], rep[1] - point[1])
if dist <= self.near_radius_m + self.margin_m:
total += math.log(eps)
continue
p_ns = min(1.0, max(0.0, (dist - self.r_eff_min_m)
/ (self.r_eff_max_m - self.r_eff_min_m)))
reading = math.atan2(rep[1] - point[1], rep[0] - point[0])
survivors = 0
for p in subsample:
if math.dist(p, point) > self.r_eff_max_m + self.margin_m:
continue
gap = angular_gap(math.degrees(
math.atan2(p[1] - point[1], p[0] - point[0])),
math.degrees(reading))
if gap <= self.half_angle_deg + _ANG_EPS_DEG:
survivors += 1
keep_dir = max(survivors / len(subsample), eps)
total += p_ns * log_keep_ns + (1.0 - p_ns) * math.log(keep_dir)
return total / max(len(representatives), 1)
def fraction_within(self, samples, point: Vec2, radius_m: float) -> float:
if len(samples) == 0:
return 0.0
reach = radius_m + _DIST_EPS_M
if _np is not None and isinstance(samples, _np.ndarray):
distance = _np.hypot(samples[:, 0] - point[0], samples[:, 1] - point[1])
return float(_np.count_nonzero(distance <= reach)) / len(samples)
count = sum(1 for p in samples if math.dist(p, point) <= reach)
return count / len(samples)
def expected_posterior_radius(self, point: Vec2, subsample,
representatives) -> float:
"""一次测量后期望后验半径(米):选点主指标,直指清除证书阈值。
对每个代表假设(真值候选)按 r_eff~U(1000,1500) 解析积分
no_signal/direction 分支,near 分支(dist <= 5+m)后验为 5m 盘;
各分支半径用幸存子采样的质心最大距离近似(偏差 << 格距,排序用)。
与几何版 worst-posterior 同族,但作用于含负观测的非凸幸存集,
且以期望代替最坏(代表点覆盖幸存集主轴)。
"""
if not self.alive or len(subsample) == 0 or len(representatives) == 0:
return 0.0
if _np is not None and isinstance(subsample, _np.ndarray) \
and isinstance(representatives, _np.ndarray):
return self._expected_posterior_radius_vector(
point, subsample, representatives)
eps_radius = self.spacing_m / 2.0
total = 0.0
for rep in representatives:
dist = math.hypot(rep[0] - point[0], rep[1] - point[1])
if dist <= self.near_radius_m + self.margin_m:
total += self.near_radius_m + self.margin_m
continue
p_ns = min(1.0, max(0.0, (dist - self.r_eff_min_m)
/ (self.r_eff_max_m - self.r_eff_min_m)))
members_ns = [p for p in subsample
if math.dist(p, point) > self.r_eff_min_m - self.margin_m]
reading = math.atan2(rep[1] - point[1], rep[0] - point[0])
members_dir = []
for p in subsample:
if math.dist(p, point) > self.r_eff_max_m + self.margin_m:
continue
gap = angular_gap(math.degrees(
math.atan2(p[1] - point[1], p[0] - point[0])),
math.degrees(reading))
if gap <= self.half_angle_deg + _ANG_EPS_DEG:
members_dir.append(p)
r_ns = self._branch_radius(members_ns) if members_ns else eps_radius
r_dir = self._branch_radius(members_dir) if members_dir else eps_radius
total += p_ns * r_ns + (1.0 - p_ns) * r_dir
return total / max(len(representatives), 1)
\end{lstlisting}
\subsection*{问题三:信息增益选点与滚动任务排序}
统一任务表由站点扫频、定位测量与清除三类任务组成,每执行一项即重排;定位点按期望信息增益加移动代价评分。
\noindent\small \texttt{BayesLoop}:问题三闭环:网格候选集的正负观测更新、信息增益选点与滚动任务排序。
\begin{lstlisting}[style=corecode]
class BayesLoop(OmnidirectionalLoop):
"""问题3贝叶斯闭环:负观测感知 + 信息增益选点 + 滚动任务排序。"""
def _build_tasks(self) -> list[Task]:
tasks: list[Task] = []
for index, point in enumerate(self.ledger.points):
pending = self._fused_channels_at(index)
if not pending:
continue
switches = len(pending) - (1 if self.client.state.channel in pending else 0)
tasks.append(Task("station", tuple(point),
self.config.cost.measure_s * len(pending)
+ self.config.cost.switch_s * max(switches, 0),
station_index=index, station_channels=tuple(pending)))
for channel in self.belief.channels_by_status(ChannelStatus.DETECTED):
cert = self._best_certificate(channel)
if cert is not None and cert.certified \
and channel not in self._clear_failed:
tasks.append(Task("clear", cert.center,
self.config.cost.clear_success_s,
channel=channel))
continue
if self._defer_to_station(channel):
continue # 下一条方位线可由"本就要去的站点"顺路提供
point = self._select_point(channel)
if point is None:
continue
tasks.append(Task("locate", point,
self.config.cost.measure_s + self.config.cost.switch_s,
channel=channel))
return tasks
def _execute_station(self, index: int) -> None:
"""单站批量扫频:只测待测频道(跳测规则已过滤),扫完即返回重排。
P0-1:本站同时测量"已发现未定证"频道,使第二/第三条方位线在巡回
途中顺路取得;证书一旦成立即就地清除,不额外绕行。
"""
point = self.ledger.points[index]
pending = self._fused_channels_at(index)
unknown = [c for c, b in self.belief.channels.items()
if b.status is ChannelStatus.UNKNOWN]
self.skipped_measures += len(unknown) - len(
[c for c in pending if self.belief.channels[c].status
is ChannelStatus.UNKNOWN])
# 当前频道优先,减少换频;同站 near 直接清除不打断扫频
pending.sort(key=lambda c: (c != self.client.state.channel, c))
for channel in pending:
if self.stop_reason is not None or self._count_limit_reached() \
or self._virtual_cap_reached():
return
detected_before = (self.belief.channels[channel].status
is ChannelStatus.DETECTED)
outcome = self.client.measure(point, channel)
self.ledger.mark_measured(channel, index)
self._record_measure(channel, point, outcome)
if outcome.result in (DIRECTION, NEAR):
self._may_clear.add(channel)
if outcome.result == NEAR:
region = self.belief.candidate_set(channel)
cert = clear_certificate(region, channel, self.config)
if cert.certified:
if self._attempt_clear(channel, point, cert):
self._clear_failed.discard(channel)
else:
self._clear_failed.add(channel)
self._may_clear.discard(channel)
else:
self._anomaly("near_without_certificate", channel, position=point)
continue
if not detected_before:
continue
if outcome.result == NO_SIGNAL:
# 该频道已被发现:保证接收盘内不可能 no_signal,属守恒性异常
region = self.belief.candidate_set(channel)
safe = (safe_disk_radius(region, self.config)
if region is not None else -1.0)
if safe >= 0 and math.dist(point, region.mec()[0]) <= safe:
self._anomaly("no_signal_in_safe_disk", channel,
position=point)
continue
cert = self._best_certificate(channel)
if cert is not None and cert.certified:
if self._attempt_clear(channel, cert.center, cert):
self._clear_failed.discard(channel)
else:
self._clear_failed.add(channel)
self._may_clear.discard(channel)
def _record_measure(self, channel: int, position: Vec2, outcome) -> None:
"""belief(凸证书)与 PF(幸存集)同步更新;PF 空即证缺席。"""
self.belief.update_measure(channel, position, outcome.result,
outcome.svd_deg)
# 同一 0.5m 测点 × 同一频道的重复检测:历史上一次官方演练有 98.84%
# 的虚拟时间耗在空复扫上,本计数是该问题的回归指标。
key = _point_key(position)
seen = self._measured_points.setdefault(channel, set())
if key in seen:
self.redundant_rescan_measures += 1
seen.add(key)
self._clear_failed.discard(channel) # 新证据落地:允许再次尝试清除
belief = self.belief.channels[channel]
if belief.status in (ChannelStatus.CLEARED, ChannelStatus.ABSENT):
return # 已证毕:清除/缺席后的 no_signal 与距离无关,不进 PF
pf = self.particles[channel]
pf.apply_measure(position, outcome.result, outcome.svd_deg)
if not pf:
if belief.status is ChannelStatus.UNKNOWN:
self.belief.mark_absent(channel)
self.absence_evidence[channel] = "pf_empty"
if channel not in self.pf_absent_channels:
self.pf_absent_channels.append(channel)
else:
self._anomaly("pf_empty_while_detected", channel, position=position)
self._apply_count_bound()
def _best_certificate(self, channel: int) -> ClearCertificate | None:
"""凸证书优先;负观测挖洞使 PF 幸存集更紧时用网格证书补位。"""
region = self.belief.candidate_set(channel)
cert = clear_certificate(region, channel, self.config) \
if region is not None else None
if cert is not None and cert.certified:
return self._near_clear_target(cert)
if self.config.bayes.grid_cert_enabled:
pf = self.particles[channel]
if len(pf):
center, radius = pf.mec()
limit = (self.config.field.clear_radius_m - pf.margin_m
- self.config.search.clear_margin_m)
if radius <= limit:
return self._near_clear_target(
ClearCertificate(channel, center, radius, limit, True))
return cert
def _select_point(self, channel: int) -> Vec2 | None:
"""定位测量点:期望 log 幸存占比 + λ·移动时间 最小的候选。"""
context = self._probe_context(channel)
if context is None:
return None
candidates, info_score, center, position = context
bayes = self.config.bayes
def score(q: Vec2) -> float:
# 一次定位后仍需接近源才能清除;仅计到检测点的移动会偏好绕行。
move = bayes.lambda_move_s * (move_time(position, q, self.config.cost)
+ move_time(q, center, self.config.cost))
return info_score(q) + move
return min(candidates, key=lambda q: (score(q), _point_key(q)))
def _localize_and_clear(self, channel: int) -> bool:
"""阻塞式定位子循环(resolve 阶段安全网):PF 选点,失败退凸选点。"""
search = self.config.search
self._may_clear.add(channel)
for _ in range(search.max_localize_steps):
if self._virtual_cap_reached():
return False
cert = self._best_certificate(channel)
if cert is not None and cert.certified and channel in self._may_clear:
if self._attempt_clear(channel, cert.center, cert):
self._clear_failed.discard(channel)
return True
self._clear_failed.add(channel)
self._may_clear.discard(channel)
continue
point = self._select_point(channel)
if point is None:
region = self.belief.candidate_set(channel)
if region is None:
return False
point = select_measurement_point(
region, self.client.state.position, self.config).position
outcome = self.client.measure(point, channel)
self._record_measure(channel, point, outcome)
self._may_clear.add(channel)
if outcome.result == NO_SIGNAL:
region = self.belief.candidate_set(channel)
if region is not None and safe_disk_radius(region, self.config) >= 0 \
and math.dist(point, region.mec()[0]) \
<= safe_disk_radius(region, self.config):
self._anomaly("no_signal_in_safe_disk", channel, position=point)
cert = self._best_certificate(channel)
if cert is not None and cert.certified:
if self._attempt_clear(channel, cert.center, cert):
self._clear_failed.discard(channel)
return True
self._clear_failed.add(channel)
self._may_clear.discard(channel)
return self.belief.channels[channel].status is ChannelStatus.CLEARED
\end{lstlisting}
\subsection*{问题四:朝向候选、可见性选点与网格遍历清除}
每个幸存格附朝向档掩码构成联合候选集,可见性权重用于探点排序;定位停滞时以 $20$ m 方格遍历清除兜底。
\noindent\small \texttt{JointBelief}:单频道的位置与朝向联合候选集;朝向分 $24$ 档,负观测只在整格必落在最小接收半径内时才约束朝向。
\begin{lstlisting}[style=corecode]
class JointBelief:
"""单频道 (位置 × 朝向) 联合可行集。"""
def is_absent(self) -> bool:
"""联合集空 ⟹ 该频道不存在源(真源保留定理的二维推广)。"""
return not self.position.alive
def sample(self, size: int = 256) -> list[int]:
"""均匀抽样幸存格(评分用;证书始终用完整集合)。"""
cached = self._sample_cache
if cached is not None and cached[0] == self.position.version:
return cached[1]
alive = sorted(self.position.alive)
if len(alive) > size:
step = len(alive) / size
alive = [alive[int(k * step) % len(alive)] for k in range(size)]
self._sample_cache = (self.position.version, alive)
return alive
def expected_posterior_radius(self, reader: Vec2, size: int | None = None) -> float:
"""一次测量后的期望后验 MEC 半径(米);D1 值函数。
与 :meth:`ChannelParticles.expected_posterior_radius` 同族,差别只在
no_signal 分支概率:定向源不仅会因距离 > r_eff 静默,还会因朝向遮挡
静默,因此
p_vis(p) = P(r_eff ≥ |p−q|) × cover_ratio(p, q)
其中 ``cover_ratio`` 是该格朝向掩码中"可能覆盖 q 方位"的档位占比。
分支半径语义:
- direction(可见,非 near):读数为 q 处的方位 ±δ,后验收缩到该方位
楔形与当前可行集之交,半径取幸存子采样的质心最大距离;
- no_signal:定向源的静默**不删空间候选**,位置后验半径保持当前值;
- near:后验为 5m 盘。
返回值只用于**选点排序**,不参与任何安全判定;证书仍由凸/网格口径给出。
"""
pf = self.position
if not pf.alive:
return 0.0
if size is None:
size = self.config.bayes.score_subsample
indices = self.sample(max(2, size))
if not indices:
return 0.0
stride = max(len(indices) // max(self.config.bayes.mc_representatives, 1), 1)
representatives = indices[::stride]
if _np is not None:
return self._epr_vector(reader, indices, representatives)
return self._epr_python(reader, indices, representatives) # pragma: no cover
def prob_detect(self, reader: Vec2, size: int = 256) -> float:
"""在 ``reader`` 处能看见任一可行 (位置, 朝向) 的概率质量(启发式)。
r_eff ~ U[1000, 1500] 解析积分;朝向按可行档位占比计权。
只用于选点排序,不参与任何安全判定。
"""
cells = self.sample(size)
if not cells:
return 0.0
total = 0.0
span = self.field.r_eff_max_m - self.field.r_eff_min_m
for index in cells:
gx, gy = self.position.cell_center(index)
distance = math.hypot(reader[0] - gx, reader[1] - gy)
if distance > self.field.r_eff_max_m:
continue
if distance <= self.field.r_eff_min_m:
p_range = 1.0
else:
p_range = (self.field.r_eff_max_m - distance) / span
if p_range <= 0.0:
continue
bearing = math.degrees(math.atan2(reader[1] - gy, reader[0] - gx))
cover = _COVER_TABLE[_table_index(bearing)]
bits = bin(self._mask(index) & cover).count("1")
if bits:
total += p_range * bits / HEADING_BINS
return total / len(cells)
def apply(self, station: Vec2, result: str,
svd_deg: float | None = None) -> None:
"""按读数更新联合可行集(保守,真源 (位置, 朝向) 永不被删)。"""
if result == DIRECTION:
self.position.apply_measure(station, DIRECTION, svd_deg,
directional=True)
self._apply_heading(station, keep=True)
elif result == NEAR:
# near 同样要求可见(评测机先判朝向再判距离)
self.position.apply_measure(station, NEAR, directional=True)
self._apply_heading(station, keep=True)
elif result == NO_SIGNAL:
# 距离阈值必须取**保守侧**:只有"整格都必然落在 r_eff_min 圆盘内"
# 时,no_signal 才一定由朝向遮挡解释。格心到格内任一点 ≤ margin,
# 故阈值取 r_eff_min − margin;取 + margin 会把 1000~1010m 的模糊带
# 也当成"必须不覆盖",从而可能删掉真源的朝向档(破坏真源保留性)。
self._apply_heading(station, keep=False,
within_m=max(
0.0, self.field.r_eff_min_m
- self.position.margin_m))
else:
raise ValueError(f"未知测向结果: {result!r}")
def summary(self) -> dict:
"""联合集规模快照(写入策略 summary,供审计与论文复算)。"""
cells = self.sample(256)
bins = [bin(self._mask(index)).count("1") for index in cells]
return {"cells": len(self.position),
"sampled_mean_heading_bins": (round(sum(bins) / len(bins), 2)
if bins else 0.0),
"heading_bins_total": HEADING_BINS}
\end{lstlisting}
\noindent\small \texttt{DirectionalExtension}:问题四闭环:定向覆盖计划、遮挡解释与网格遍历清除兜底。
\begin{lstlisting}[style=corecode]
class DirectionalExtension(OmnidirectionalLoop):
"""问题4混合源闭环(P5):全向基线 + 定向观测处理 + 后备清除。"""
def _localize_and_clear(self, channel: int) -> bool:
"""定位子循环:选点测量收敛候选集,拿到清除证书后清除。
与 P4 的差异:保证接收盘内的 no_signal 说明探点被朝向遮挡
[观测推理],走"对跖复核 -> 环探"升级梯子(不计入定位步数,
由 max_blocked_probes 限流),之后仍未收敛再转后备清除。
"""
search = self.config.search
self._may_clear.add(channel)
steps = 0
while steps < search.max_localize_steps:
if self._virtual_cap_reached():
return False
region = self.belief.candidate_set(channel)
if region is None:
return False
cert = clear_certificate(region, channel, self.config)
if cert.certified and channel in self._may_clear:
if self._attempt_clear(channel, cert.center, cert):
return True
self._may_clear.discard(channel) # 证书失败:补测后才可重试
continue
if len(self.blocked_probes.get(channel, [])) >= search.max_blocked_probes:
break # 梯子用尽仍被遮挡:定位停滞,转后备清除
point = self._next_probe(region, channel)
within_range = (math.dist(point, region.mec()[0])
<= safe_disk_radius(region, self.config))
outcome = self.client.measure(point, channel)
self._observe(channel, point, outcome.result, outcome.svd_deg)
self._may_clear.add(channel)
if outcome.result == NO_SIGNAL:
if within_range and channel not in self.directional_observed:
self.directional_observed.append(channel)
self.blocked_probes.setdefault(channel, []).append(point)
continue # 盘外无信号也可能是距离不足,不据此断言定向
self.blocked_probes.pop(channel, None)
steps += 1
if self.belief.channels[channel].status is ChannelStatus.CLEARED:
return True
region = self.belief.candidate_set(channel)
if region is None:
return False
cert = clear_certificate(region, channel, self.config)
if cert.certified and channel in self._may_clear:
if self._attempt_clear(channel, cert.center, cert):
return True
self._may_clear.discard(channel)
return self._fallback_clear(channel, region)
def _observe(self, channel: int, point: Vec2, result: str,
svd_deg: float | None) -> None:
"""同步更新凸候选集与 (位置 × 朝向) 联合可行集。
联合集只在频道首次被看见(direction/near)时建立:从未被看见的
频道由覆盖证明负责缺席,不需要朝向维度。联合集一旦为空,直接给出
该频道的缺席证明(严格,与覆盖证明同一条件),可提前停止对其检测。
"""
self.belief.update_measure(channel, point, result, svd_deg)
if not self.config.search.directional_joint_belief:
return
joint = self.joint.get(channel)
if joint is None:
if result not in (DIRECTION, NEAR):
return
joint = JointBelief(self.config)
self.joint[channel] = joint
joint.apply(point, result, svd_deg)
if joint.is_absent() and self.belief.channels[channel].status \
is not ChannelStatus.CLEARED:
self.belief.mark_absent(channel)
self.absence_evidence[channel] = "joint_empty"
if channel not in self.joint_absent_channels:
self.joint_absent_channels.append(channel)
def _next_probe(self, region: CandidateRegion, channel: int) -> Vec2:
"""定向遮挡下的探点升级梯子 [工程]。
- 未被遮挡过:worst-posterior 选点器(保证接收盘内);
- 第 1 次被遮挡:对 MEC 圆心取对跖方向,偏移裁进保证接收盘;
- 之后:绕 MEC 圆心、半径 = ρ 的环探,方位自首个被遮挡点起
每次 +45°。几何边界必须如实陈述:环上可见弧没有下界——
源可位于环上(|g-c| = ρ)且朝向沿外切线,可见集退化为单个
切点,固定步进环探可以全部漏过。梯子只提高提前命中概率,
不构成接收或命中保证;完备性完全由后备清除兜底。
"""
visibility = self._visibility_probe(region, channel)
if visibility is not None:
return visibility
blocked = self.blocked_probes.get(channel, [])
center, radius = region.mec()
if not blocked:
return select_measurement_point(region, self.client.state.position,
self.config).position
safe = safe_disk_radius(region, self.config)
first = blocked[0]
if len(blocked) == 1 and safe >= 0:
dx, dy = center[0] - first[0], center[1] - first[1]
norm = math.hypot(dx, dy)
if norm <= 1e-9:
# 首个被遮挡探点恰在圆心:方向信息为零,沿 +x 走满安全盘。
return (center[0] + safe, center[1])
scale = min(1.0, safe / norm)
# 对跖方向,偏移模长 = min(norm, safe):圆心对称点裁进保证接收盘。
return (center[0] + dx * scale, center[1] + dy * scale)
base = math.atan2(first[1] - center[1], first[0] - center[0])
azimuth = base + math.radians(45.0 * (len(blocked) - 1))
return (center[0] + radius * math.cos(azimuth),
center[1] + radius * math.sin(azimuth))
def _fallback_clear(self, channel: int, region: CandidateRegion) -> bool:
"""定位停滞的后备:候选集有限圆盘覆盖逐格清除(格数即请求上界)。"""
search, session = self.config.search, self.config.session
width = max(p[0] for p in region.polygon) - min(p[0] for p in region.polygon)
height = max(p[1] for p in region.polygon) - min(p[1] for p in region.polygon)
cells = fallback_clear_cells(region.polygon, self.config)
entry = {"channel": channel, "cells": len(cells),
"bound": fallback_cell_bound(width, height, self.config),
"posterior_radius_m": round(region.mec()[1], 3)}
self.fallback_log.append(entry)
if len(cells) > search.fallback_max_cells:
# 宽候选集盲扫违背预算纪律:留待定位收敛,不得硬扫 [计划§4]。
entry["outcome"] = "skipped_cell_cap"
self._anomaly("fallback_cell_cap_exceeded", channel,
cells=len(cells), cap=search.fallback_max_cells)
return False
# 现实预算核对:N 次请求 × 单请求高分位估计 + 退出余量,不足不启动。
remaining = self.client.real_remaining_s
if remaining is None or len(cells) * session.real_rtt_p95_s \
+ session.reserve_exit_s > remaining:
entry["outcome"] = "skipped_budget"
return False
attempts = 0
order = self._fallback_order(cells)
for cell in order:
if self._virtual_cap_reached():
entry["outcome"] = "skipped_virtual_budget"
entry["attempts"] = attempts
return False
attempts += 1
outcome = self.client.clear(cell, channel)
if outcome.result == CLEAR_SUCCESS:
self.belief.mark_cleared(channel)
self.cleared.append(channel)
entry["outcome"] = "cleared"
entry["attempts"] = attempts
return True
# 候选集保守包含真源,格心距格内任一点 <= 10*sqrt(2) < 20m:全部
# 失败说明守恒性被打破,按异常留痕,不得据此删除候选。
entry["outcome"] = "exhausted"
entry["attempts"] = attempts
self._anomaly("fallback_exhausted_without_success", channel,
cells=len(cells))
return False
\end{lstlisting}
\endgroup