diff --git a/README.md b/README.md index 355d51e..64b6981 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,10 @@ - [四问建模分析与 P0–P7 实施计划](docs/B题建模分析与实施计划.md) - [协议规则表、字段表与状态机（P0 交付）](docs/协议规则表与状态机.md) - [问题一完整解答](docs/问题一解答.md) -- [问题三搜索与总耗时改进方案（含七点覆盖证明和自建消融）](docs/algo-3.md) +- [问题一算法说明（当前实现口径）](docs/algo-1.md) +- [问题二算法说明（保证接收的第二点选择）](docs/algo-2.md) +- [问题三算法说明（七点覆盖证明、贝叶斯闭环与滚动调度）](docs/algo-3.md) +- [问题四算法说明（定向覆盖最小命中集、联合信念与后备清除）](docs/algo-4.md) - [第三、四问系统性架构优化与运行提速方案（运行历史复核、最小命中集、联合信念）](docs/第三四问架构优化与运行提速方案.md) - [第三、四问图论与动态规划优化方案（定点实测、路由 gap、定向定位 DP、证书硬化）](docs/第三四问图论与动态规划优化方案.md) - [1-4 问最优与基线实现对比框架](project/README.md)（`uv run python -m project.run`，结果在 `project/results/`） @@ -140,8 +143,8 @@ client.exit() | P2 几何内核（角域交、直径、MEC、证书） | ✅ 完成（迁入 `mcm26b.geometry`） | | P3 问题二选点（K1、候选区域 Q、采样评估） | ✅ 完成（belief + 选点器，测试固化） | | P4 问题3全向闭环（49点覆盖、定位队列、清除、停止） | ✅ 官方问题3演练单案例通过（10/10全清，[验证报告](output/q3/practice/20260911-port2026/report.md)） | -| P5 问题4定向扩展（联合约束、后备清除、预算） | ✅ 自建完成（官方演练待执行） | -| P6 加速与消融 | ✅ 七点批量（`--strategy seven`，官方演练 12/12）+ 贝叶斯闭环（`--strategy bayes` 默认，自建消融 −30.0%，[算法说明](docs/algo-3-bayes.md)） | +| P5 问题4定向扩展（联合约束、后备清除、预算） | ✅ 自建完成（官方演练待执行；[算法说明](docs/algo-4.md)） | +| P6 加速与消融 | ✅ 七点批量（`--strategy seven`，官方演练 12/12）+ 贝叶斯闭环（`--strategy bayes` 默认，自建消融 −30.0%，[算法说明](docs/algo-3.md)、[设计细节](docs/algo-3-bayes.md)） | | P7 图论/动态规划优化（S1–S5） | ✅ 已落地（[实测回填](docs/第三四问图论与动态规划优化方案.md#8-本轮实测回填2026-09-12-实施结果)）：路由度量入库、覆盖证书硬化（15 m/72 缺口 566→0）、Q4 后备盲清 125.7→4.3 次/局；同时修复联合信念 `no_signal` 朝向规则反向缺陷 | ### 图论/DP 轮（2026-09-12，S1–S5）要点 diff --git a/docs/algo-1.md b/docs/algo-1.md index 81cbe69..e073d28 100644 --- a/docs/algo-1.md +++ b/docs/algo-1.md @@ -1,185 +1,266 @@ -# 问题一：有界测向误差下定位区域的直径与圆覆盖 +# 问题一算法说明：有界测向误差下的角域交会、直径与最小包围圆 -## 摘要 +> 本文按**当前仓库实现**重写，替代早期只描述理论方案的版本。 +> 问题一生产入口为 `src/mcm26b/geometry/localization.py::solve`，CLI 为 +> `python -m mcm26b.geometry`；问题二至四只复用其中的 +> `smallest_enclosing_circle` 与同一套角域半平面语义，并不直接调用 +> `solve`（复用与差异见 §7）。 -针对同一干扰源在多个检测点处的带误差方位观测，建立前向角域交集模型，将定位区域表示为二维线性不等式组的可行域。首先推导角域的半平面表达，给出兼顾空集、点、线段和无界区域的几何构造算法；随后证明，有界定位区域的直径必由一对顶点取得，且判断是否存在直径等于区域直径的覆盖圆，只需检验任意一组最远顶点的中点。通过符合题设测向条件的等边三角形反例，说明直径圆一般不能覆盖定位区域。进一步以最小包围圆连接定位精度和后续清除条件，并给出算法复杂度、数值处理原则及既有演练观测的复算结果。全部结论以有界误差模型为前提，不引入误差独立、正态或均匀分布假设。 +## 1. 问题与实现入口 -**关键词：** 交会定位；有界误差；半平面交；凸多边形直径；最小包围圆。 +同一频道的同一静止干扰源，在若干检测点 $s_i=(s_{ix},s_{iy})$ 处得到示向度 +$\theta_i$。题设测向误差不超过 $\pm1^\circ$，同一地点重复测量误差不变； +示向度从正东逆时针计量。问题一要求定位区域的直径，并判断“以该直径为 +直径的圆”能否覆盖定位区域。 -## 1. 问题界定与思路来源 - -题目要求由若干检测点坐标及同一干扰源的示向度，求交会定位多边形的直径，并判断以该直径为直径的圆能否覆盖定位区域。题面规定测向误差在 $[-1^\circ,1^\circ]$ 内，同一地点重复测量的误差不变；示向度从正东沿逆时针方向计量。[题面第1页、附录2第(1)、(3)项](C:/Users/HelixCopex/code/MATH/issue/B题.pdf) - -由此得到三步建模思路。第一，单次观测不能确定一条精确射线，却可以确定一个包含真实方向的角域，因此采用集合相交而非直线最小二乘。第二，角域可以写成半平面交，问题随之转化为凸几何问题。第三，“两点最远距离”和“覆盖全部位置所需半径”是不同的优化量，应分别计算直径和最小包围圆。 - -半平面、凸组合和凸集交运算属于标准凸分析工具，其背景见 Boyd 与 Vandenberghe 的《Convex Optimization》第2章；本文针对测向模型给出独立推导，不将一般理论本身视为新的研究贡献。[文献2](https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf) - -本文严格区分两个集合：问题一的 $P$ 仅表示测向角域之交；第二问及在线搜索中的 $K$ 可进一步结合源位置圆域和接收距离先验。向 $P$ 加入圆盘后，边界可能出现圆弧，不能继续把它当成原问题中的纯多边形。 - -## 2. 符号及基本假设 - -| 符号 | 含义 | -| --- | --- | -| $s_i=(s_{ix},s_{iy})$ | 第 $i$ 个检测点，坐标单位为米 | -| $\theta_i$ | 第 $i$ 次有效示向度 | -| $\delta_i$ | 误差半角，题设取 $1^\circ$ | -| $g=(x,y)$ | 未知干扰源位置 | -| $W_i$ | 第 $i$ 次测向对应的闭前向角域 | -| $P$ | 所有角域的交集 | -| $V=\{v_1,\ldots,v_m\}$ | 非空有界 $P$ 的顶点集 | -| $D(P)$ | 区域内任意两点距离的上确界 | -| $B(c,r)$ | 中心为 $c$、半径为 $r$ 的闭圆盘 | -| $\rho(P)$ | 覆盖 $P$ 的最小圆盘半径 | - -假设各次有效示向度对应同一频道的同一个静止源。仅把 `direction` 观测代入角域模型，不能给 `near` 或 `no_signal` 人为赋予角度。理论推导允许 $0\le\delta_i<90^\circ$，用于讨论退化情况，题目主模型始终采用 $\delta_i=1^\circ$。 - -程序可额外计入两位小数舍入对应的 $0.005^\circ$ 裕量，但必须显式声明。题设的 $1^\circ$ 误差界与实现中的浮点容差是两个概念，均不能在结果解释时省略。 - -## 3. 由角域推导半平面约束 - -定义方向向量 $u(\alpha)=(\cos\alpha,\sin\alpha)$，二维叉积为 $[a,b]=a_xb_y-a_yb_x$。三角函数计算前将角度换成弧度。 - -令 $u_i^-=u(\theta_i-\delta_i)$、$u_i^+=u(\theta_i+\delta_i)$。当 $0<\delta_i<90^\circ$ 时，真实位置位于两条边界射线之间，等价于 $[u_i^-,g-s_i]\ge0$ 和 $[u_i^+,g-s_i]\le0$。第一条要求 $g-s_i$ 在下边界的左侧，第二条要求其在上边界的右侧。 - -将第一条展开，得到 $\sin(\theta_i-\delta_i)x-\cos(\theta_i-\delta_i)y\le\sin(\theta_i-\delta_i)s_{ix}-\cos(\theta_i-\delta_i)s_{iy}$。 - -将第二条展开，得到 $-\sin(\theta_i+\delta_i)x+\cos(\theta_i+\delta_i)y\le-\sin(\theta_i+\delta_i)s_{ix}+\cos(\theta_i+\delta_i)s_{iy}$。 - -同时保留前向约束 $u(\theta_i)\cdot(g-s_i)\ge0$。该约束在正半角情形下冗余，但当 $\delta_i=0$ 时，前两条只限定一条直线，加入它才得到正确的射线。角域顶点本身也包含在闭集合内，这是保守几何处理，并不表示该点仍能获得示向度。 - -因此，$P=\bigcap_{i=1}^{n}W_i=\{g:Ag\le b\}$。采用向量不等式后，$359^\circ$ 与 $0^\circ$ 跨界无需分段比较角度。每个 $W_i$ 是闭凸集，故 $P$ 也是闭凸集；它可能为空、为单点、为线段、为有界凸多边形，也可能无界。 - -## 4. 定位区域的构造与分类 - -### 4.1 边界交点与可行性 - -把约束统一写成 $a_jx+b_jy\le c_j$。两条边界的行列式为 $\Delta_{jk}=a_jb_k-a_kb_j$。当 $\Delta_{jk}\ne0$ 时，交点坐标为 $x=(c_jb_k-c_kb_j)/\Delta_{jk}$、$y=(a_jc_k-a_kc_j)/\Delta_{jk}$。保留满足全部约束的交点，再去重并构造凸包。 - -**命题1：** 至少有一次观测时，非空的角域交集一定存在极点，因此可通过可行边界交点枚举判空。 - -**证明：** 取第一个角域的轴向量 $u_1$，其中任意 $g$ 满足 $u_1\cdot(g-s_1)\ge\cos\delta_1\,\lVert g-s_1\rVert$，且 $\cos\delta_1>0$。在非空 $P$ 上最小化线性函数 $u_1\cdot(g-s_1)$；以任一可行点的函数值截取下水平集，得到非空紧集，故最小值能够取得。最小化点集位于一条支撑直线上，且有界，所以为点或线段。取该点或线段端点即得到 $P$ 的极点。二维线性约束的极点至少有两条法向量线性无关的有效边界，否则沿所有有效边界的公共切向可作双向小扰动，与极点定义矛盾。因此极点包含在上述交点枚举中。证毕。 - -无观测时，$P=\mathbb R^2$，应直接判为无界，不能用“没有交点”判断为空集。 - -### 4.2 无界性判定 - -对非空 $P$ 定义衰退锥 $\mathcal R=\{d:Ad\le0\}$。若存在 $d\ne0$，则任意 $g_0\in P$ 与 $t\ge0$ 都满足 $g_0+td\in P$，所以 $P$ 无界。 - -反之，若 $P$ 无界，可取 $g_k\in P$ 使 $\lVert g_k-g_0\rVert\to\infty$。单位向量 $(g_k-g_0)/\lVert g_k-g_0\rVert$ 存在收敛子列，极限 $d$ 满足 $\lVert d\rVert=1$。将 $Ag_k\le b$ 减去 $Ag_0$，再除以距离并取极限，得到 $Ad\le0$。因此 $P$ 有界当且仅当 $\mathcal R=\{0\}$。 - -本模型的衰退锥被一个小于 $180^\circ$ 的闭角域约束，若包含非零向量，必含一条极射线，该射线平行于某条有效约束边界。故只需检查每个法向量 $(a_j,b_j)$ 的两个切向方向 $(b_j,-a_j)$、$(-b_j,a_j)$ 是否满足所有齐次不等式。 - -### 4.3 分类约定 - -可行交点为空时返回空集；存在非零衰退方向时返回无界。否则求可行交点的凸包，按顶点数和仿射维数区分单点、线段和二维多边形。空集的直径在本文约定为未定义；无界集合的直径为 $+\infty$。程序使用空值及单独的无穷标志表示二者，避免写入非标准 JSON 数值。 - -显示窗口只用于绘图裁剪。任意大矩形都不能取代衰退锥判定，更不能把截断图形的直径称为原无界区域的直径。 - -## 5. 直径的顶点定理 - -**定理1：** 若 $P=\operatorname{conv}(V)$ 非空且有界，则 $D(P)=\max_{i,j}\lVert v_i-v_j\rVert$。 - -**证明：** 任取 $x,y\in P$，可写成 $x=\sum_i\lambda_i v_i$、$y=\sum_j\mu_jv_j$，其中权重非负且分别求和为1。于是 $x-y=\sum_{i,j}\lambda_i\mu_j(v_i-v_j)$。由三角不等式，$\lVert x-y\rVert\le\sum_{i,j}\lambda_i\mu_j\lVert v_i-v_j\rVert\le\max_{i,j}\lVert v_i-v_j\rVert$。最远顶点对又属于 $P$，上界可以达到，结论成立。证毕。 - -因此，小规模问题中直接枚举所有顶点对即可，计算复杂度为 $O(m^2)$。单点的直径为0，线段的直径为端点距离。这里的直接枚举不是近似算法；其缺点仅在于大规模时计算量较高。 - -## 6. 直径圆覆盖的充要条件 - -**定理2：** 取任意一组最远顶点 $a,b$，令 $D=\lVert a-b\rVert$、$o=(a+b)/2$。存在直径为 $D$ 的覆盖圆，当且仅当所有顶点均满足 $\lVert v_i-o\rVert\le D/2$。 - -**证明：** 若半径为 $D/2$、中心为 $z$ 的圆同时包含 $a,b$，则 $D\le\lVert a-z\rVert+\lVert z-b\rVert\le D$。两个不等式必须同时取等号，所以 $z$ 在两点连线上，且到两点距离相等，只能是中点 $o$。当 $D=0$ 时结论直接成立。另一方面，圆盘是凸集，包含全部顶点等价于包含它们的凸包。因此只需在这个唯一可能的圆心处检查全部顶点。证毕。 - -该定理不仅回答“以某条直径线段为直径的圆是否覆盖”，也回答“能否移动一个同样大小的圆使其覆盖”。若检查失败，移动圆心或更换最远顶点对均不能解决。 - -定义最大超出量 $E=\max_i\lVert v_i-o\rVert-D/2$。在精确算术下，$E=0$ 表示覆盖成立，$E>0$ 表示失败；不会出现严格负值，因为最远点对就在圆周上。取达到最大值的顶点，即得到一份可以独立复核的覆盖失败证据。 - -## 7. 符合测向条件的反例 - -取误差半角均为 $1^\circ$，三个检测点及读数如下。 - -| 检测点 | 坐标（米） | 示向度 | +| 层次 | 文件 | 内容 | | --- | --- | --- | -| $s_1$ | $(-1000,0)$ | $1^\circ$ | -| $s_2$ | $(540,-500\sqrt3)$ | $121^\circ$ | -| $s_3$ | $(520,520\sqrt3)$ | $241^\circ$ | +| 生产内核 | `src/mcm26b/geometry/localization.py` | `Measurement`、`halfplanes`、`convex_hull`、`diameter`、`smallest_enclosing_circle`、`solve` | +| CLI | `src/mcm26b/geometry/__main__.py` | 读 JSON、可加两位小数舍入裕量、输出结果 | +| 统一入口 | `src/mcm26b/cli.py` | `mcm26b q1` | +| 对比框架 | `project/q1/{optimal,baseline,reference,cases,compare}.py` | 生产实现 / 朴素基线 / 独立大框参考解 | +| 回归 | `tests/test_geometry.py`、`tests/test_geometry_certificates.py` | 解析形状、真值保留、随机对拍、证书边界 | +| 结果 | `examples/q1/*.json`、`project/results/q1.md` | 固定解析算例与 40 个随机交会算例的配对结果 | -每个角域的下边界方向分别为 $0^\circ,120^\circ,240^\circ$。对应的三个左半平面是 $y\ge0$、$\sqrt3x+y\le40\sqrt3$、$y\le\sqrt3x$，交集为顶点 $A=(0,0)$、$B=(40,0)$、$C=(20,20\sqrt3)$ 的等边三角形。 +```bash +.venv/bin/python -m mcm26b.geometry examples/q1/triangle.json +.venv/bin/python -m mcm26b.geometry examples/q1/normal.json --rounding-margin-deg 0.005 +``` -还需验证每个角域的另一条边不会切去这个三角形。以 $s_1$ 为例，到三个顶点的方向分别为 $0^\circ,0^\circ,\arctan(20\sqrt3/1020)$，最后一项约为 $1.945^\circ<2^\circ$，均在 $[0^\circ,2^\circ]$ 内。其余两站由绕三角形中心旋转 $120^\circ$ 的对称性成立。因此三个完整前向角域之交恰为上述三角形。 +`--rounding-margin-deg` 默认 $0$；库内在线求解默认把服务器两位小数舍入的 +$\le0.005^\circ$ 偏差计入有效半角，这是**实现口径**，不是题面新增误差。 +问题一 CLI 只有显式传参才附加该裕量，以便同时报告纯 $1^\circ$ 与 +$1.005^\circ$ 两套敏感性结果。 -将真实源取为三角形中心 $g_*=(20,20\sqrt3/3)$。以第一站为例，真实方向约为 $0.649^\circ$，与 $1^\circ$ 的读数偏差小于 $1^\circ$；检测距离约为1020米，小于1500米。另两站同理。该反例符合题设方位误差和可能的有效接收半径。 +## 2. 模型：把角域写成半平面交 -此时 $D=40$ 米，而以 $AB$ 为直径的圆中心为 $(20,0)$，到 $C$ 的距离为 $20\sqrt3>20$ 米，不能覆盖。三角形外接圆即最小包围圆，其中心为 $(20,20\sqrt3/3)$，半径为 $40/\sqrt3\approx23.094011$ 米。这也证明了“直径不超过40米”不足以保证一次清除。 +令方向向量 $u(\alpha)=(\cos\alpha,\sin\alpha)$，二维叉积 +$[a,b]=a_xb_y-a_yb_x$。设误差半角为 $\delta$（题设主模型 $\delta=1^\circ$）， +真实源 $g$ 属于闭前向角域 -## 8. 最小包围圆及其清除含义 +$$ +W(s,\theta,\delta)=\{g:\ [u(\theta-\delta),g-s]\ge0,\ +[u(\theta+\delta),g-s]\le0,\ u(\theta)\cdot(g-s)\ge0\}. +$$ -### 8.1 存在性、唯一性与支撑点 +对应实现中的三条半平面 $a x+b y\le c$ 为 -定义 $f(z)=\max_{v\in V}\lVert v-z\rVert$，则 $\rho(P)=\min_z f(z)$。由于 $V$ 有限，$f$ 连续且随 $\lVert z\rVert\to\infty$ 而趋于无穷，最小值存在。 +| 约束 | $a$ | $b$ | $c$ | +| --- | --- | --- | --- | +| 下边界左侧 | $\sin(\theta-\delta)$ | $-\cos(\theta-\delta)$ | $a s_x+b s_y$ | +| 上边界右侧 | $-\sin(\theta+\delta)$ | $\cos(\theta+\delta)$ | $a s_x+b s_y$ | +| 前向 | $-\cos\theta$ | $-\sin\theta$ | $a s_x+b s_y$ | -若两个不同中心 $c_1,c_2$ 都以半径 $r$ 覆盖 $V$，中点 $c$ 对每个顶点满足 $\lVert v-c\rVert^2=(\lVert v-c_1\rVert^2+\lVert v-c_2\rVert^2)/2-\lVert c_1-c_2\rVert^2/420$。 + +## 4. 伪代码 + +```text +solve(measurements): + if measurements 为空: return unbounded + 平移/缩放，生成每条读数的三条半平面 + V = [] + for 每一对边界直线: + 求交点 p + if p 满足全部半平面: V.append(p) + if V 为空: return empty + if 存在非零衰退切向 d: return unbounded + H = convex_hull(V) + D, (a,b) = argmax_{v_i,v_j∈H} |v_i-v_j| + o = (a+b)/2 + covers = max_{v∈H} |v-o| <= D/2 + 输出 diameter、diameter_circle、witness + (m,rho) = Welzl_MEC(H) + 输出 minimum_enclosing_circle 与 rho+eps<=20 清除证书 +``` + +## 5. 复杂度与数值实现 + +| 步骤 | 复杂度 | +| --- | --- | +| 生成半平面 | $O(k)$，$k=3n$ | +| 边界交点枚举与可行性判定 | $O(k^3)$ | +| 凸包 | $O(k^2\log k)$ 最坏 | +| 直径 | $O(m^2)$ | +| 直径圆检查 | $O(m)$ | +| Welzl 最小圆 | 随机排列下期望 $O(m)$；固定排列不声称最坏线性 | + +主要数值约定： + +- 坐标归一化后计算，输出换回原坐标并重新复核包含性； +- 仅对**精确平行**直线对跳过求交；近平行对保留； +- 位置可行性容差随坐标幅值缩放；衰退锥判定使用独立紧阈值； +- `solve` 主动限制最多 100 个检测点，超过时要求改用半平面交算法； +- 所有输出用 `null`/显式标志表示空集和无界，不写非法 JSON 浮点数。 + +## 6. 已验证算例与结果 + +固定解析算例（`examples/q1`，当前代码实测）： + +| 算例 | 状态 | 直径 (m) | MEC 半径 (m) | 直径圆覆盖 | 最大超出 (m) | MEC 可直接清除 | +| --- | --- | ---: | ---: | --- | ---: | --- | +| `normal` 三点交会 | polygon | 26.986 | 13.943 | 否（2 个圆外顶点） | 3.964 | 是 | +| `triangle` 等边三角形反例 | polygon | 40.000 | 23.094 | 否（1 个圆外顶点） | 14.641 | 否 | +| `segment` 零误差退化线段 | segment | 100.000 | 50.000 | 是 | 0.000 | 否 | +| `poor` 近共线检测点 | unbounded | — | — | — | — | — | +| `empty` 不相容观测 | empty | — | — | — | — | — | + +`triangle` 的三个顶点是 $(0,0)$、$(40,0)$、$(20,20\sqrt3)$：以直径 +$AB$ 为直径的圆圆心 $(20,0)$、半径 20，到 $C$ 的距离为 $20\sqrt3>20$， +故不能覆盖；最小包围圆半径 $40/\sqrt3\approx23.094$，也不能一次清除。 +其中 $A=(0,0)$ 坐标是浮点下的 $-6.4\times10^{-13}$，上表按解析值取整。 + +既有官方演练观测（案例 `U5ZW-P9KN-TCAG-R7ST`，不公开源坐标、未执行清除） +的复算结果： + +| 误差半角 | 直径 (m) | MEC 半径 (m) | 直径圆最大超出 (m) | | --- | ---: | ---: | ---: | | $1^\circ$ | 23.996723 | 12.068583 | 0.604002 | | $1.005^\circ$ | 24.119088 | 12.130377 | 0.608486 | -两种口径下直径圆均不覆盖区域，最小包围半径均小于20米。该演练未公开源坐标，且没有执行清除，故只用于观测交会和计算一致性复核，不能据此声称实测定位误差或全清成功率。数据见[原始三点观测](C:/Users/HelixCopex/code/MATH/output/q1/practice/U5ZW-P9KN-TCAG-R7ST_observations.json)，敏感性数值见[已有复算报告](C:/Users/HelixCopex/code/MATH/output/q1/optimization_report.json)。 +两种口径下直径圆都不覆盖定位区域，MEC 半径都小于 20 m，但该演练没有 +真实清除动作，不能据此声称实测定位误差或全清成功率。 -## 11. 结论与实现对应 +`project/q1` 的 40 个随机交会算例（冻结算例与独立参考解对拍）： -问题一的准确答案为：将测向角域表示为半平面交，先判断非空性及有界性，再由最远顶点对得到直径；以该直径为直径的圆一般不能覆盖定位区域，覆盖的充要条件是任意一组最远顶点的中点到全部顶点的距离均不超过半径。最小包围圆提供独立且适合后续清除的精度指标。 +| 方法 | 状态准确率 | 直径平均绝对误差 | 直径圆覆盖准确率 | 平均耗时 | +| --- | ---: | ---: | ---: | ---: | +| 生产 `solve` | 100% | 0.000 m | 100% | 0.086 ms | +| 朴素基线（网格 + 包围盒对角线） | 100% | 4.583 m | 52.5% | 15.932 ms | -对应实现为 [localization.py](C:/Users/HelixCopex/code/MATH/geometry/localization.py)。在项目根目录可运行 `python -m geometry examples/q1/triangle.json`；追加 `--rounding-margin-deg 0.005` 可复算额外舍入裕量。本文整理理论与已有计算结果，不改变原始观测或既有程序。 +## 7. 与问题二/三/四当前实现的复用关系和差异（核对结果） + +问题一给出的几何内核是后续所有几何决策的基础，但**复用方式是“语义 + 最小圆函数”， +不是直接调用 `solve`**。代码层面的核对结果如下。 + +| 后续问题 | 实际复用了什么 | 实际没有复用 / 做了什么替换 | +| --- | --- | --- | +| 问题二（`planner.select_second_point`） | 同一套角域半平面语义；`smallest_enclosing_circle` 求候选集 MEC；用 MEC 半径构造保证接收盘 | 不调用 `solve`；候选集不是纯角域交，而是 $K_1=\Omega\cap W\cap B(s,1500)$；用 192 边形外包圆盘 + Sutherland–Hodgman 裁剪，而不是枚举直线交点；不计算直径，也不使用“直径中点圆”判据 | +| 问题三（`belief.CandidateRegion`、`planner.bayes_loop`） | 复用 `smallest_enclosing_circle`；角域边界公式与问题一相同；MEC 作为清除证书 | 在线信念是有界凸外包多边形，圆盘边界是外切多边形而非精确曲线；不再判无界；no_signal 的负观测由网格粒子滤波处理（问题一模型里没有负观测）；不用直径，而用“MEC 半径 + 外包余量 + 工程余量 $\le20$” | +| 问题四（`planner.DirectionalExtension`、`belief.joint`） | 继续复用问题一的 MEC 与角域语义；（位置 × 朝向）联合信念也建立在同一网格上 | 定向源的 no_signal 不能排除位置，不能沿用问题一/三的纯角度或全向距离解释；覆盖语义改为 `forward_half_disk`（见 `algo-4.md`） | + +代码证据： + +- `src/mcm26b/belief/__init__.py` 只从 `geometry.localization` 导入 + `smallest_enclosing_circle`； +- `src/mcm26b/particles.py` 同样只导入 `smallest_enclosing_circle`； +- `生产/运行代码中 `geometry.localization.solve` 的调用者只有 + `geometry/__main__.py` 和 `project/q1/optimal.py`，即问题一 CLI 与其 + 对比框架；测试文件会直接调用它做解析与随机对拍，不在此列； +- 在线 `clip_wedge` 的半平面系数与问题一 `halfplanes` 一致，但作用在 + 凸多边形裁剪流程上；`solve` 的“纯角域 + 精确顶点枚举 + 无界判定” + 是针对问题一的一站式分析工具，不适合直接嵌入带圆盘、负观测和状态更新 + 的在线闭环。 + +**差异一句话总结：** 问题一回答“固定角域观测下源可能在哪里、直径圆能否覆盖”， +是纯几何精确解；问题二至四把其中的角域半平面和最小圆升级为“带距离先验的 +保守信念 + 主动动作 + 证据证明”的闭环，问题三再加非凸负观测，问题四再加 +朝向维度。 + +## 8. 结论与适用边界 + +问题一的当前实现是：三条半平面 → 可行顶点枚举/凸包 → 衰退锥判有界 → +顶点对精确直径 → 任意直径中点唯一候选的覆盖检验 → Welzl 最小包围圆与 +20 m 证书。其结论只依赖题面有界误差，不引入误差分布假设。空集、无界、 +点、线段和多边形状态都有明确输出。 + +未解决问题一的数值敏感性时，不应把在线闭环里的外包多边形半径直接当成 +问题一的精确直径；两者口径不同。反过来，在线闭环也不能把 `solve` 直接套在 +$K_1$ 上，因为 $K_1$ 含圆盘和距离上界，边界可能含圆弧。 ## 参考文献 -1. 全国大学生数学建模竞赛组委会. 2026年高教社杯全国大学生数学建模竞赛B题：无线电干扰源的快速自动定位与清除[Z]. 2026. [本地题面](C:/Users/HelixCopex/code/MATH/issue/B题.pdf)，第1页、附录2。 -2. BOYD S, VANDENBERGHE L. Convex Optimization[M]. Cambridge: Cambridge University Press, 2004. [作者公开版本](https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf)，第2章。用于凸几何背景，本文的题目化证明在正文给出。 -3. WELZL E. Smallest enclosing disks (balls and ellipsoids)[C]//MAURER H. New Results and New Trends in Computer Science. Lecture Notes in Computer Science, vol. 555. Berlin: Springer, 1991: 359–370. [作者公开论文](https://people.inf.ethz.ch/emo/PublFiles/SmallEnclDisk_LNCS555_91.pdf)。用于随机增量最小包围圆算法来源。 -4. CGAL Project. CGAL::Min_circle_2: Bounding Volumes[EB/OL]. [在线文档](https://doc.cgal.org/latest/Bounding_volumes/classCGAL_1_1Min__circle__2.html)，访问日期：2026-09-10。用于最小圆及支撑点定义的交叉核对。 +1. 全国大学生数学建模竞赛组委会. 2026 年高教社杯全国大学生数学建模竞赛 B 题：无线电干扰源的快速自动定位与清除[Z]. 2026. 本地题面 `issue/B题.pdf`。 +2. BOYD S, VANDENBERGHE L. Convex Optimization[M]. Cambridge University Press, 2004. 凸集、半平面与凸组合背景。 +3. WELZL E. Smallest enclosing disks (balls and ellipsoids)[C]//New Results and New Trends in Computer Science, LNCS 555, 1991: 359–370. 随机增量最小包围圆。 +4. 本项目. `docs/algo-2.md`（问题二当前实现）. 问题二对本文半平面语义与最小圆的复用说明。 +5. 本项目. `docs/algo-3.md`、`docs/algo-4.md`. 在线闭环中的保守信念、负观测与联合信念差异。 diff --git a/docs/algo-2.md b/docs/algo-2.md index bddb5c5..9258a8b 100644 --- a/docs/algo-2.md +++ b/docs/algo-2.md @@ -1,248 +1,312 @@ -# 问题二：具有接收保证的第二检测点选择 +# 问题二算法说明：保证接收的第二检测点选择（当前实现） -## 摘要 +> 本文按**当前仓库实际运行路径**重写。生产入口是 +> `src/mcm26b/planner.select_second_point`，其核心选点器是 +> `select_measurement_point`；`src/mcm26b/geometry/planning.py` 中的 +> `first_bearing_candidate_disk` 是解析充分证书与独立校核工具， +> 并不在运行时被 `select_second_point` 调用（见 §6 核对结果）。 -针对只获得一次全向干扰源示向度的情形，提出以接收保证为约束、以后验最小包围半径为主要指标的主动选点方法。首先将源位置圆域、首次角域与有效接收距离上界相交，构造包含真实源的位置候选集；随后把第二检测点的安全条件写成一族闭圆盘的交，并证明其非空性等价于候选集最小包围半径不超过最小有效接收半径。利用首次角域狭窄的题设特征，进一步推导截断角域的最小包围圆，得到半径约250米、能够直接计算的安全候选圆盘。为比较候选点的定位效果，建立兼顾 `direction` 和 `near` 的保守后验模型，推导交会角的局部几何意义，并给出采用外包多边形及观测区间划分的有限候选评估算法。严格接收保证与启发式性能优化分别陈述，不将有限采样结果解释为连续问题的全局最优解。 +## 1. 问题、实现入口与交付物 -**关键词：** 主动测向；鲁棒选点；安全候选区域；最坏后验半径；保守几何近似。 +问题二给定一个全向干扰源的首次有效示向度 $(s,\theta)$，要求给出第二个 +检测点 $q$ 的选择策略、定位效果与候选区域。两个约束必须同时满足： -## 1. 问题界定与思路来源 +1. **接收保证**：$q$ 到所有仍可能为真源的位置都不超过有效接收半径下界 + $R_-=1000$ m，否则第二次测量可能返回 `no_signal`，几何收益落空； +2. **定位效果**：在保证接收的候选点中，使第二次测量后的保守位置集合 + 尽量小。当前实现以“最坏后验最小包围半径”为主指标，以移动耗时为次指标。 -题目第二问给定某全向干扰源在一个检测点处的示向度，要求提出第二检测点的选择策略、说明其定位效果并给出候选区域。这里的源类型明确为全向，不能混入第四问的朝向可见性条件。[题面第1页问题2](C:/Users/HelixCopex/code/MATH/issue/B题.pdf) +| 层次 | 文件 | 内容 | +| --- | --- | --- | +| 在线入口 | `src/mcm26b/planner/__init__.py::select_second_point` | 由首次 `direction` 构造 $K_1$，返回 `SecondPointPlan` | +| 核心选点 | 同文件 `select_measurement_point` | 任意候选集的保证接收盘内选点 | +| 几何证书 | `src/mcm26b/geometry/planning.py::first_bearing_candidate_disk` | 截断角域解析外包圆与安全候选盘 | +| 信念构造 | `src/mcm26b/belief/__init__.py::BeliefState`、`CandidateRegion` | $K_1$ 的保守凸外包与更新 | +| 对比框架 | `project/q2/{optimal,baseline,reference,region,cases,compare}.py` | 生产选点器 / 两种朴素基线 / 保证盘内网格参考 | +| 回归 | `tests/test_planner.py`、`tests/test_geometry_certificates.py`、`tests/test_belief.py` | 接收保证、后验收缩、near 分支、解析证书 | +| 结果 | `project/results/q2.{md,json}` | 30 个随机首次示向度算例配对结果 | -本问有两个先后不同的问题：第二点是否仍能收到信号；在能收到信号的点中，如何使新测量更有价值。如果只追求较大交会角而不控制到真实源的距离，第二次测量可能返回 `no_signal`，预期几何收益无法实现。因此先构造保证接收的可行域，再优化测向后的不确定区域。 +## 2. 首次观测后的候选集 $K_1$ -方法来源分为三部分。候选位置集合来自题设有界误差与接收距离，沿用问题一的角域约束；安全选点来自“对每个可能源位置都满足距离条件”的鲁棒可行性；定位收益来自主动测量设计，以最坏后验半径代替未经题面支持的平均误差。凸函数与鲁棒建模的一般背景参见 Boyd 与 Vandenberghe 的《Convex Optimization》第3章、第6.4节；本题的安全圆盘公式、扇形推导及区间评估方法在下文逐步建立。[文献3](https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf) +### 2.1 精确定义 -需要强调：本题没有给出源位置概率分布，也没有给出可直接使用的独立高斯测向噪声模型。因此不能把卡尔曼滤波、协方差或费舍尔信息计算直接作为严格保证。交会角仍有明确几何意义，但应作为局部解释和候选筛选依据。 +记首次检测点为 $s$、首次示向度为 $\theta$、误差半角 +$\delta_{\rm eff}=1^\circ+0.005^\circ$（题设误差 + 两位小数舍入裕量）， +有效接收半径上界 $R_+=1500$ m，源位置圆域 $\Omega=B(0,1800)$。则 -## 2. 符号、题设事实与建模约定 +$$ +K_1=\Omega\cap W(s,\theta,\delta_{\rm eff})\cap B(s,R_+). +$$ -| 符号 | 含义 | -| --- | --- | -| $s$、$\theta$ | 首次检测位置与示向度 | -| $\delta$ | 误差半角，默认 $1^\circ$ | -| $g$ | 未知的全向干扰源位置 | -| $\Omega=B(0,1800)$ | 干扰源位置圆域 | -| $R_g\in[R_-,R_+]$ | 未知有效接收半径，$R_-=1000$、$R_+=1500$ 米 | -| $K_1$ | 首次观测后的保守位置候选集 | -| $q$ | 第二检测点 | -| $\mathcal Q(K_1)$ | 对 $K_1$ 内所有位置均保证接收的第二点区域 | -| $\rho(K)$ | 非空有界集合 $K$ 的最小包围半径 | -| $J(q)$ | 第二点的最坏保守后验半径 | +真实源在合法观测下必属于 $K_1$。注意 $B(s,R_+)$ 的半径是**接收半径上界**： +`direction` 读数本身就隐含测得距离不超过 $R_+$；而 $B(s,R_-)$ 不在 $K_1$ +之中，因为第二次是否收到信号正是本问要保证的事。 -采用以下题设事实：有效接收半径在1000至1500米之间；检测点到全向源距离不超过5米时无法得到示向度，可直接光学定位并清除；清除距离不超过20米；机器狗移动速度为5米/秒，检测需要5秒，同一频道连续测量不发生额外频道切换。[题面附录2第(2)、(4)至(9)项](C:/Users/HelixCopex/code/MATH/issue/B题.pdf) +### 2.2 实际保守表示 -1800米圆域限制的是源位置，不是机器狗的搜索位置。检测点可以位于该圆域外；接口坐标分量的绝对值上限属于数值防护限制，不能误解成要求 $q\in\Omega$。[题面第4页](C:/Users/HelixCopex/code/MATH/issue/B题.pdf) +在线实现不直接处理圆弧边界，而用以下外切多边形表示 $K_1$： -如果把读数舍入额外计入误差预算，可显式采用 $\delta=1.005^\circ$。下文解析推导对 $0\le\delta\le45^\circ$ 成立，但这一扩展区间不是题目误差范围。 +- 先取 $\Omega$ 的 192 边外切正多边形 $P_M(0,1800)$； +- 再与 $B(s,1500)$ 的 192 边外切正多边形求交； +- 最后用 Sutherland–Hodgman 逐半平面裁剪两条测向边界。 -## 3. 首次观测对应的源位置集合 +设圆盘外切多边形相对圆盘的最大径向超出为 $R(1-\cos(\pi/M))$。实现把每次 +圆盘裁剪的该量累加到 `CandidateRegion.slack_m`，因此维护的多边形满足 +$K_1\subseteq P_{\rm out}$，且 `slack_m` 是可审计的外包半径余量。 -设 $u=u(\theta)=(\cos\theta,\sin\theta)$。由第一次有效示向度，真实位置属于前向角域 $W(s,\theta,\delta)$；由已经收到信号，必有 $\lVert g-s\rVert\le R_g\le R_+$；由源位置限制，又有 $g\in\Omega$。故定义 $K_1=\Omega\cap W(s,\theta,\delta)\cap B(s,R_+)$。 +这条路径与问题一的区别：问题一在纯角域上枚举边界直线交点，精确处理 +空集/无界；问题二候选集一开始就被 $\Omega$ 和 $B(s,R_+)$ 限成有界集， +因此可以只用凸多边形裁剪。角域系数与问题一 `halfplanes` 完全一致： +$\sin(\theta-\delta)x-\cos(\theta-\delta)y\le c$ 与 +$-\sin(\theta+\delta)x+\cos(\theta+\delta)y\le c$。 -该集合为闭凸紧集。真实观测一致时 $g\in K_1$，所以它非空。首次返回 `direction` 还意味着 $\lVert g-s\rVert>5$，但若把这个开圆盘排除，会破坏凸性。本文暂不使用该排除信息，保留包含真源的凸外包集合；由此得到的保证可能更保守，不会误删真实位置。 +## 3. 保证接收区域 -这一区别决定了后续结论的口径：本文给出的保证对较大的 $K_1$ 中所有位置成立，因而对实际源成立；但在 $K_1$ 上定义的最坏后验指标可能比利用全部观测细节的精确模型更悲观。 +### 3.1 一般集合形式 -若计算得到 $K_1=\varnothing$，说明输入、频道、误差界或数值处理与题设不一致。此时应返回不一致状态，不能继续给空集计算圆心。问题一的纯角域求解器也不能直接把上述圆弧边界当作直线边界处理。 +对非空紧候选集 $K$，定义 -## 4. 安全候选区域的精确定义及基本性质 +$$ +F_K(q)=\max_{g\in K}\lVert q-g\rVert,\qquad +\mathcal Q(K)=\{q:F_K(q)\le R_-\}=\bigcap_{g\in K}B(g,R_-). +$$ -### 4.1 接收距离条件 +$F_K$ 是凸函数且 1-Lipschitz，$\mathcal Q(K)$ 是闭凸集。 -若第二检测点满足 $\lVert q-g\rVert\le R_-$，由于 $R_g\ge R_-$，即可保证全向源信号能够接收。因此定义 $F_K(q)=\max_{g\in K}\lVert q-g\rVert$，并令 $\mathcal Q(K)=\{q:F_K(q)\le R_-\}=\bigcap_{g\in K}B(g,R_-)$。 +**定理 1.** $\mathcal Q(K)\ne\varnothing$ 当且仅当 +$\rho(K)=\min_q F_K(q)\le R_-$，其中 $\rho(K)$ 是 $K$ 的最小包围半径。 -这里限制的是第二点到所有可能源位置的距离，而不是到第一次检测点的距离。保证接收也不等于保证返回示向度，因为距离不超过5米时可能返回 `near`；该分支的定位效果更好，不能视为策略失败。 +证明：$\min_q F_K(q)$ 就是最小包围半径；存在 $q$ 使 $F_K(q)\le R_-$ 与 +其最小值不超过 $R_-$ 等价。证毕。 -### 4.2 凸性、稳定性和非空条件 +### 3.2 多边形外包给出的充分安全盘 -**命题1：** 对任意非空紧集 $K$，函数 $F_K$ 为凸函数，且满足 $\lvert F_K(q_1)-F_K(q_2)\rvert\le\lVert q_1-q_2\rVert$。 +设在线维护的多边形 $P_{\rm out}\supseteq K$，其最小包围圆为 +$B(m,\rho_{\rm out})$。任取 -**证明：** 对每个固定的 $g$，距离函数关于 $q$ 凸；逐点取最大值保持凸性。又由三角不等式，$\lVert q_1-g\rVert\le\lVert q_2-g\rVert+\lVert q_1-q_2\rVert$，对 $g$ 取最大值得到一个方向的不等式，再交换两点即得结论。证毕。 +$$ +q\in B(m,\ R_- - \rho_{\rm out}-\eta),\qquad \eta=10^{-3}\ {\rm m}, +$$ -因此 $\mathcal Q(K)$ 是闭凸集。非空时，它包含在任一 $B(g_0,R_-)$ 内，所以也有界。 +则对任意真源 $g\in K$： -**定理1：** $\mathcal Q(K)\ne\varnothing$ 当且仅当 $\rho(K)\le R_-$。 +$$ +\lVert q-g\rVert\le\lVert q-m\rVert+\lVert m-g\rVert +\le (R_--\rho_{\rm out}-\eta)+\rho_{\rm out} +=R_--\eta<1000. +$$ -**证明：** 最小包围半径的定义为 $\rho(K)=\min_q F_K(q)$。存在 $q$ 使 $F_K(q)\le R_-$，与其最小值不超过 $R_-$ 完全等价。最小值存在，因为 $F_K$ 连续且随 $\lVert q\rVert\to\infty$ 趋于无穷。证毕。 +实现中 `safe_disk_radius(region)` 直接返回 +$R_--\rho_{\rm out}-10^{-3}$；因为 $\rho_{\rm out}$ 已是包含真实候选集 +的多边形最小包围半径，这个盘内的点到真实源的三角不等式上界严格小于 +1000 m。1 mm 的数值余量避免浮点贴边被 $r_{\rm eff}=1000$ 的源判成 +`no_signal`。若该半径 $<0$，说明当前凸外包的 MEC 半径已超过 1000 m： -当 $\rho(K)=R_-$ 时，安全区域仅包含唯一最小包围圆中心；唯一性可由两个候选中心的中点距离恒等式证明，见[问题一第8节](C:/Users/HelixCopex/code/MATH/docs/algo-1.md)。因此应使用不严格不等式，不能把等号情形误判为空。 +- 这**不能**证明真实 $\mathcal Q(K)$ 为空，只说明这个充分证书不够紧； +- 实现退化为探索 MEC 圆心，并在 `reason` 中明确标注“可能 `no_signal`”。 -### 4.3 外包圆给出的充分区域 +### 3.3 首次截断角域的解析候选盘 -若已知 $K\subseteq B(c,r)$ 且 $r\le R_-$，则任意满足 $\lVert q-c\rVert\le R_--r$ 的点，都有 $\lVert q-g\rVert\le\lVert q-c\rVert+\lVert c-g\rVert\le R_-$。因此 $B(c,R_--r)\subseteq\mathcal Q(K)$。 +完全去掉 $\Omega$ 后，完整的截断角域 +$S=W(s,\theta,\delta)\cap B(s,R_+)$ 有一个精确最小包围圆： -这给出了容易计算的充分候选区域，但一般不等于整个安全区域。一个宽松外包圆若满足 $r>R_-$，只能说明该外包圆不能提供证书，不能证明真实的 $\mathcal Q(K)$ 为空。 +$$ +c=s+a\,u(\theta),\qquad a=\frac{R_+}{2\cos\delta},\qquad +S\subseteq B(c,a), +$$ -## 5. 首次截断角域的解析最小包围圆 +对 $0\le\delta\le45^\circ$ 成立。简证：设 +$g=s+t\,u(\theta+\varphi)$，$0\le t\le R_+$， +$|\varphi|\le\delta$，则 +$\lVert g-c\rVert^2=t^2+a^2-2at\cos\varphi$。固定 $\varphi$ 时该式 +关于 $t$ 是凸二次函数，最大值在 $t=0$ 或 $t=R_+$ 处：前者为 $a^2$， +后者不超过 $R_+^2+a^2-2aR_+\cos\delta=a^2$（代入 +$a=R_+/(2\cos\delta)$），故 $S\subseteq B(c,a)$。最小性由三点 +$s$、$s+R_+u(\theta\pm\delta)$ 的加权距离恒等式给出：权重 +$\lambda_\pm=1/(4\cos^2\delta)$、 +$\lambda_0=1-1/(2\cos^2\delta)$ 非负且和为 1，且 +$c=\lambda_0s+\lambda_+p_++\lambda_-p_-$，故任何覆盖这三点的圆半径 +至少为 $a$。于是 -### 5.1 坐标变换 +$$ +\mathcal Q_0=B(c,\ R_- - a) +$$ -暂时去掉只会缩小集合的 $\Omega$ 约束，考虑完整截断角域 $S=W(s,\theta,\delta)\cap B(s,R_+)$。平移并旋转，使 $s$ 成为原点、$u$ 成为横轴。其点可写成 $g=s+t\,u(\theta+\varphi)$，其中 $0\le t\le R_+$、$-\delta\le\varphi\le\delta$。 +是保证接收的充分候选区域。题设 $\delta_{\rm eff}\le1.005^\circ$ 时： -令 $a=R_+/(2\cos\delta)$，并取 $c=s+au$。 - -### 5.2 包含性证明 - -**定理2：** 对 $0\le\delta\le45^\circ$，截断角域 $S$ 的最小包围圆为 $B(c,a)$。 - -**证明的第一部分：** 对任意 $g\in S$，有 $\lVert g-c\rVert^2=t^2+a^2-2at\cos\varphi$。固定 $\varphi$ 时，右侧关于 $t$ 为凸二次函数，区间最大值在端点取得。当 $t=0$ 时等于 $a^2$；当 $t=R_+$ 时，由 $\cos\varphi\ge\cos\delta$ 得其不超过 $R_+^2+a^2-2aR_+\cos\delta=a^2$。故整个 $S$ 被该圆覆盖。 - -### 5.3 最小性证明 - -取角域中的三个点 $p_0=s$、$p_+=s+R_+u(\theta+\delta)$、$p_-=s+R_+u(\theta-\delta)$。它们均在 $B(c,a)$ 的圆周上。 - -令 $\lambda_+=\lambda_-=1/(4\cos^2\delta)$，$\lambda_0=1-1/(2\cos^2\delta)$。由于 $\delta\le45^\circ$，三权重非负且和为1，直接代入可得 $c=\lambda_0p_0+\lambda_+p_++\lambda_-p_-$。 - -对任意中心 $z$，利用加权距离恒等式，有 $\sum_j\lambda_j\lVert p_j-z\rVert^2=a^2+\lVert c-z\rVert^2\ge a^2$。若半径为 $r$ 的圆覆盖这三个点，则左侧不超过 $r^2$，故 $r\ge a$。任何覆盖 $S$ 的圆都必须覆盖这三个点，结合第一部分的包含性，得到最小半径恰为 $a$。当 $\delta=0$ 时两个远端点重合，结论退化为长度 $R_+$ 的线段直径圆，仍成立。证毕。 - -该推导不依赖源位置先验，也不依赖位置概率分布。它给出的是完整截断角域的精确最小圆；与 $\Omega$ 相交后，真实 $K_1$ 的最小圆可能更小。 - -## 6. 本题可直接采用的第二点候选区域 - -由 $K_1\subseteq S\subseteq B(c,a)$ 和第4.3节，可取 $\mathcal Q_0=B(c,b)$，其中 $c=s+au$、$a=1500/(2\cos\delta)$、$b=1000-a$。 - -| 误差半角 | 轴向中心距离 $a$ / 米 | 安全候选半径 $b$ / 米 | +| 误差半角 | 轴向中心距离 $a$ / m | 安全候选半径 $R_- - a$ / m | | --- | ---: | ---: | | $1^\circ$ | 750.114246 | 249.885754 | | $1.005^\circ$ | 750.115391 | 249.884609 | -两种口径下 $b>0$，故本题首次有效测向后的安全候选区域可以直接证明非空。无需将“若安全区域为空就去最小圆圆心”作为正常流程的保底步骤。 +该解析盘由 `first_bearing_candidate_disk` 实现并测试。它给出的是 +$S$ 的最小圆；由于 $K_1\subseteq S$，$B(c,R_- - a)\subseteq\mathcal Q(K_1)$ +是严格充分的。实际在线路径使用比 $S$ 更紧的 $K_1$ 凸外包，因此其安全盘 +通常大于该解析下界（例：`q2-4200` 的 $\rho_{\rm out}=674.59$ m， +安全盘半径 325.41 m，而解析下界约 249.89 m）。 -几何上，$\mathcal Q_0$ 位于首测示向方向前方约750米，以约250米为半径。位于圆盘轴向前后端的点、两侧横向点以及圆心都保证接收。横向点只是可供比较的候选点，并未被证明总是最优。 +## 4. 实际选点算法 -为抵消坐标计算和动作提交中的工程舍入，可取缩小区域 $\mathcal Q_\eta=B(c,b-\eta)$，其中 $0<\eta5$ 的候选位置，令真实方位为 $h_q(g)=\operatorname{atan2}(g_y-q_y,g_x-q_x)$。可能出现的读数集合定义为 $\mathcal B(q)=\{\beta:\exists g\in K_1,\ \lVert g-q\rVert>5,\ d_{\mathbb S^1}(\beta,h_q(g))\le\delta\}$。 +1. 在候选集边界上以 `witness_spacing_m=60` m 采样见证点 + $g$（`region.witnesses`，含全部顶点）； +2. 对每个 $g$，真实方位 $\phi=\operatorname{bearing}(q,g)$，枚举允许的 + 读数区间端点 $\phi\pm\delta_{\rm eff}$； +3. 对每个端点，复制候选多边形并用 `clip_wedge(q, 读数, δ_eff)` 得到后验 + 凸外包，记录其 MEC 半径； +4. 额外处理 `near` 分支：以 $B(q,5)$ 裁剪候选集，若交集非空，记录其 MEC + 半径；它代表“源离第二点不超过 5 m，直接光学定位清除”的分支； +5. 以上全部有效半径取最大值，就是该候选点的评估值。 -这里用所有允许误差描述可能读数，不是假设各地点的误差可以随机独立重抽样。同一个位置的重复检测不产生新的独立信息；所构造的 $\mathcal Q_0$ 与首测位置相距至少约500米，天然避免原地重复测量。 +该评估是**启发式上界近似**：它只对目标边界采样和读数区间端点取值， +没有覆盖连续集合上所有可能读数；但它可在候选点之间做一致比较。 +1/4 间距（15 m）的加密复核只提高输出数值的可信度，不改变“采样启发式” +这一事实。定理 1 保证接收，选点最优性没有全局证明。 -### 7.2 最坏后验半径指标 +### 4.3 清除与 `near` -定义 $J(q)$ 为所有可能 `direction` 分支的 $\rho(K_d(q,\beta))$，以及非空 `near` 分支的 $\rho(K_n(q))$ 的上确界。等价地，它是未来所有允许观测下最坏的保守最小包围半径；不存在的观测分支不参与取最大值。 +若第二次测量返回 `near`，源在 $q$ 的 5 m 内。由 `clear_certificate`： +候选集 $B(q,5)$ 的 MEC 半径至多略大于 5 m（外包多边形余量很小）， +工程余量后仍远小于 20 m，因此可以在当前位置直接 `/clear`，不需要再取 +第三个测向读数。若返回 `direction`，则用新角域更新 $K$ 并继续按 +问题三/四的定位流程处理；后验半径若仍大于 20 m，**不能**声称两次定位 +已经完成。 -选点的主要目标为在安全候选区域中降低 $J(q)$。该指标直接服务于一次清除要求：若 $J(q)\le20$，则每个可能分支都有半径不超过20米的保守定位圆。它比仅比较面积更符合任务，因为面积较小的集合仍可能十分狭长。 +## 5. 参数与复杂度 -对于任意有效分支，后验都是 $K_1$ 的子集。集合包含关系给出 $\rho(K_2)\le\rho(K_1)$，故 $J(q)\le\rho(K_1)$。但不一定严格小于，因此不能仅由“增加一次观测”推断定位精度必然严格改善，更不能保证只需两次测向便达到20米阈值。 +| 参数 | 默认值 | 作用 | +| --- | ---: | --- | +| `polygon_sides` | 192 | 圆盘外切多边形边数；$\Omega$ 处径向超出约 0.24 m | +| `bearing_rounding_margin_deg` | 0.005° | 两位小数舍入的显式裕量 | +| `witness_spacing_m` | 60 m | 粗评边界见证间距 | +| `candidate_rings` | (0, 0.4, 0.75, 1.0) | 安全盘内半径比例 | +| `candidate_angles` | 12 | 每环方向数；连同圆心共 37 个候选 | +| `near_radius_m` | 5 m | `near` 分支判定 | +| `clear_margin_m` | 0.5 m | 清除证书的额外工程余量 | -### 7.3 移动成本的引入 +若候选集顶点数为 $m$、候选点数为 $N_q$、边界见证数为 $N_g$，单次选点的 +粗评估约为 $O(N_qN_g\cdot{\rm clip}+N_q{\rm MEC})$，`clip` 为 +Sutherland–Hodgman 裁剪的 $O(m)$ 代价；再加一次 1/4 间距加密复核。 +这是明显可并行/缓存的小规模在线计算，不追求全局最优。 -连续测量同一频道时，本次动作的虚拟时间增量为 $C(q)=\lVert q-s\rVert/5+5$ 秒。推荐先按后验半径择优，再在几何效果接近的候选点中比较移动成本，以避免任意混合米和秒。 +## 6. 实现核对结果：文档旧说法与当前运行路径的差异 -用 $\widehat J(q)$ 表示实际计算的评估量，在第10节的保守方案中取为上界 $\overline J(q)$。设有限候选点中的最小评估值为 $J_{\min}$。可先保留满足 $\widehat J(q)\le J_{\min}+\tau$ 的点，再选择其中 $C(q)$ 最小者，其中 $\tau\ge0$ 是显式给定的定位容忍量，单位为米；$\tau=0$ 对应严格的字典序选择。本文不预设一个未经实验支持的最优 $\tau$。 +早期说明只写到 `geometry/planning.py` 的解析候选盘，容易让人误以为 +`select_second_point` 直接调用它。当前代码的实际情况是： -## 8. 交会角的局部几何解释 +- `select_second_point` 构造 `CandidateRegion`（$\Omega\cap B(s,1500)$ 的 + 192 边形外包）并用 `clip_wedge` 交入首次角域； +- 真正的选点发生在 `select_measurement_point`：它在上述凸外包的 MEC + 保证盘内做 37 点候选评估； +- `first_bearing_candidate_disk` 是解析证书函数，用于证明 + $B(s+au,1000-a)$ 的充分性，并被 + `tests/test_geometry_certificates.py` 以稠密角域采样校核； +- 生产路径不调用 `first_bearing_candidate_disk`，但两者给出同一性质的 + 保证：选出的点必能收到 `direction`/`near`。生产路径通常比解析盘更紧， + 因为它利用了 $K_1$ 完整形状而不是仅用其外包扇区。 -在某个假定真实位置 $g$ 附近，令两站距离为 $r_1,r_2$，真实视线单位向量为 $u_1,u_2$，对应法向量为 $n_1,n_2$。当角误差很小且位置扰动相对站距较小时，方位函数的一阶展开给出 $n_i\cdot\Delta g\approx r_i e_i$，其中 $e_i$ 用弧度表示，满足 $\lvert e_i\rvert\le\delta_{\rm rad}$。 +因此，问题二的**严格保证**来自“安全盘点到候选集所有点的距离 $\le R_-$”， +解析盘是充分但不唯一的一种实现方式；**最优性**来自有限候选评估， +两者必须在说明中分开陈述。 -于是局部不确定区域近似由两条带状约束限定：$\lvert n_1\cdot\Delta g\rvert\le r_1\delta_{\rm rad}$、$\lvert n_2\cdot\Delta g\rvert\le r_2\delta_{\rm rad}$。若视线交会角为 $\phi$，法向量矩阵的行列式绝对值为 $\lvert\sin\phi\rvert$。对法向坐标中的矩形作线性逆变换，得到近似面积 $A_{\rm loc}\approx4r_1r_2\delta_{\rm rad}^2/\lvert\sin\phi\rvert$。 +## 7. 已验证结果 -这解释了两点：近于平行或反向的视线使误差显著放大；在距离不变的局部比较下，接近直角的交会较有利。然而真实源位置未知，$r_2$ 和 $\phi$ 随 $g$ 改变，且有限角域并非无限带状区域。因此不能由该近似推出“总应横移固定距离”或“在某个假定源位置形成90度就是全局最优”。本文只将其用于优先考虑轴线两侧候选点，最终仍以后验集合指标评价。 +`project/q2` 在 30 个随机首次示向度算例上，用同一 $K_1$、同一评估器 +对四种方法配对比较（初始 $\rho$ 均值 681.47 m）： -## 9. 圆弧先验的保守离散 +| 方法 | 保证接收率 | 最坏后验均值 | 终值比 $\rho'/\rho$ | 平均移动 | 结论 | +| --- | ---: | ---: | ---: | ---: | --- | +| 生产选点器（保证盘内最小最坏后验） | **100%** | **73.89 m** | **0.104** | 170.4 s | 保证接收且后验最小 | +| 基线① 垂直 1000 m 最大交会角 | **0%** | 80.11 m | 0.116 | 200.0 s | 忽略接收距离，真源在远端必 `no_signal` | +| 基线② 直接取 MEC 圆心 | 100% | 340.75 m | 0.500 | 136.3 s | 保证接收但几乎不收缩 | +| 独立参考：保证盘内 9×9 网格枚举 | 100% | 75.63 m | 0.106 | 170.3 s | 生产采样接近网格最优 | -### 9.1 外包切线多边形 +代表性单例 `q2-4200`： -直接将圆周采样点连接得到的是内接多边形，会错误排除圆盘边缘的真实位置。为了复用第一问的半平面工具，选择 $M\ge3$ 个均匀法向量 $v_j=u(2\pi j/M)$，定义 $P_M(o,R)=\{x:v_j\cdot(x-o)\le R,\ j=0,\ldots,M-1\}$。 +| 指标 | 数值 | +| --- | ---: | +| $K_1$ MEC 圆心 | (1172.85, −225.68) | +| $K_1$ MEC 半径 $\rho$ | 674.59 m | +| 保证接收盘半径 | 325.41 m | +| 选中第二点 | (1257.07, 88.64) | +| 最坏后验半径 | 54.45 m | +| 移动耗时 | 181.42 s | -由柯西不等式，任意 $x\in B(o,R)$ 都满足所有约束，因此 $B(o,R)\subseteq P_M(o,R)$。相邻切线交点到圆心的距离为 $R\sec(\pi/M)$，故单个圆盘外包的最大径向超出不超过 $R[\sec(\pi/M)-1]$。 +## 8. 与问题一及问题三/四的关系（核对结果） -令 $\widetilde K_1=P_M(0,1800)\cap W(s,\theta,\delta)\cap P_M(s,1500)$，便有 $K_1\subseteq\widetilde K_1$。对任何第二读数再求半平面交，仍然保留真实源。上述单圆径向误差不能未经论证直接当成最终定位半径的误差界；靠近切触或退化时尤其需要加密检查。 +| 对比 | 问题一 | 问题二（本文） | 问题三/四中的去向 | +| --- | --- | --- | --- | +| 信息模型 | 纯角域 $W$ | 角域 $W$ + 距离先验 $K_1=\Omega\cap W\cap B(s,R_+)$ | 问题三在任意后验 $K_c$ 上重复使用；问题四再叠加朝向维度 | +| 求解方式 | 枚举半平面边界交点、凸包、精确直径 | 圆盘外切多边形 + Sutherland–Hodgman 裁剪 | 问题三/四沿用 `CandidateRegion` 凸外包 | +| 决策类型 | 对固定观测做几何分析 | 主动选择下一个测量点 | 问题三是多源闭环；问题四是混合源闭环 | +| 主要输出 | 直径与直径圆覆盖 | `SecondPointPlan`（位置、最坏后验、移动耗时） | 问题三把它推广到定位子过程；问题四把它作为联合信念选点的基础候选之一 | +| 保证口径 | 几何精确（在输入模型内） | 接收严格保证；最坏后验是采样启发式 | 清除证书严格；选点与调度是启发式 | -安全选点仍使用第6节的解析 $\mathcal Q_\eta$。即使粗糙外包多边形的最小圆偏大，也不能据此推翻解析安全区域的非空结论。 +关键区别： -### 9.2 多边形安全区域的可选扩展 +1. 问题二不直接复用问题一的 `solve`，而是复用其角域语义和 + `smallest_enclosing_circle`；候选集多出圆盘和距离上界，不能当作 + “纯角域问题一”处理。 +2. 问题一的输出是直径和直径圆，问题二的清除判断改用最小包围半径 + $\rho\le20$；直径圆本身不参与在线清除。 +3. 问题三对同一个 `select_measurement_point` 做任意后验集上的复用， + 但同时引入网格粒子滤波处理 $K$ 无法表示的负观测（`no_signal`）； + 问题二没有负观测更新，因为它假设第二次测量点已在保证接收盘内。 +4. 问题四不能直接使用“安全盘内必收到信号”的结论：定向源在安全盘内 + 也可能因朝向遮挡返回 `no_signal`。问题四保留问题二的安全盘作为 + 基础候选，但新增（位置 × 朝向）联合信念来判断可见性（见 + `algo-4.md`）。 -若采用任意有界凸外包多边形 $\widetilde K=\operatorname{conv}\{v_j\}$，则 $\max_{g\in\widetilde K}\lVert q-g\rVert=\max_j\lVert q-v_j\rVert$，证明与第一问的凸组合论证相同。因此其保证接收区域由有限个约束 $\lVert q-v_j\rVert\le1000$ 给出，是圆盘交。 +## 9. 结论 -也可先解最小包围圆问题，再构造充分安全圆盘。由于 $\widetilde K$ 是外包集合，得到的区域可能更小，但其中每个点仍安全。若该离散区域为空，优先加密或回到解析候选盘，不把离散失败误报为真实问题无解。 +问题二的当前实现可以概括为： -## 10. 未来读数的保守区间评估 +- **候选集**：$K_1=\Omega\cap W(s,\theta,\delta_{\rm eff})\cap B(s,R_+)$， + 用 192 边形外包圆盘 + 半平面裁剪保守维护； +- **保证接收**：$B(m,R_--\rho_{\rm out}-10^{-3})$ 内的点必不返回 + `no_signal`；解析扇区最小圆给出可独立复核的充分候选盘 + $B(s+au,1000-a)$，$a=1500/(2\cos\delta)$； +- **选点**：在保证盘内取 37 个候选，以边界见证点 × 读数端点的最坏后验 + 包围半径为第一指标、移动时间为第二指标，选出后再做 1/4 间距加密复核； +- **清除**：若第二次返回 `near`，直接就地清除；否则把新角域交入后验， + 未达到 20 m 证书前继续测量。 -直接枚举少量假定源位置和误差端点，只得到被枚举情形的最大值，可能漏掉真正最坏读数。为了使评估输出具有上界意义，可以划分全部可能方位角，而非仅点采样。 - -将整周按圆周意义划分为若干区间 $I_j$，中心为 $\beta_j$，宽度为 $h_j$，要求 $\delta+h_j/2<90^\circ$。跨越 $0^\circ$ 的区间按圆周距离处理。 - -**命题2：** 对任意 $\beta\in I_j$，有 $W(q,\beta,\delta)\subseteq W(q,\beta_j,\delta+h_j/2)$。 - -**证明：** 对左侧任意非顶点位置，其方向到 $\beta$ 的圆周距离不超过 $\delta$；而 $\beta$ 到区间中心的距离不超过 $h_j/2$。由圆周距离的三角不等式，其方向到中心的距离不超过两者之和。顶点在两个闭角域中均被保留。证毕。 - -因此定义 $U_j(q)=\widetilde K_1\cap W(q,\beta_j,\delta+h_j/2)$，则任何落在该区间的真实读数所对应的保守后验，都包含于 $U_j(q)$，从而 $\rho(K_d(q,\beta))\le\rho(U_j(q))$。 - -对于 `near` 分支，可以取统一上界 $\min\{5,\rho(\widetilde K_1)\}$。该上界不需要精确计算圆弧交集；即使某点实际不会触发 `near`，纳入此上界只增加保守性。 - -令 $\overline J(q)$ 等于所有非空 $U_j(q)$ 的最小包围半径，以及上述 `near` 上界的最大值，即得到 $J(q)\le\overline J(q)$。该结论在精确几何计算下成立。浮点实现还需给出包含性和数值误差处理，不能把无误差预算的计算值称为机器可验证的严格上界。 - -通过将某个角区间二分，子区间扩张后的角域包含于父区间的扩张角域，所以相应方向上界不会增大。通过保留旧切线并增加新切线，也可使外包多边形嵌套缩小。这给出稳定的自适应加密方向，但没有在本文中证明连续选点的统一收敛速率。 - -该评估可能较慢，实际可先用少量读数采样粗筛，再对入围点使用完整角区间上界。前一阶段仅决定计算资源分配；最后报告的保守上界必须来自完整覆盖的区间，不能来自粗筛样本最大值。 - -## 11. 可执行的第二点选择算法 - -**算法2：安全区域内的有限候选主动测向。** - -1. 输入 $s,\theta,\delta$，核验首次观测为同一全向源的有效示向度，构造 $K_1$ 或其保守表示。 -2. 计算 $a=1500/(2\cos\delta)$、$c=s+au(\theta)$、$b=1000-a$，选择满足误差预算的缩小半径 $b-\eta$。 -3. 构造有限候选集：包含圆心 $c$，并取 $q_{\ell j}=c+r_\ell[\cos\psi_j\,u+\sin\psi_j\,u_\perp]$，其中 $01000$ | 不存在对该集合全部位置保证接收的点；圆心仅是最小化最坏距离的探索选择 | -| 外包圆半径大于1000 | 只能判定该证书不足，不能直接判真实安全区域为空 | -| 测量轴接近 $0^\circ/360^\circ$ | 使用向量约束和圆周区间，避免角度断点错误 | -| 两次测向交会较差 | 保留完整后验并继续测量，不将小面积当成小直径 | -| 源类型为定向 | 距离证书仍成立，但接收还受朝向限制，不能直接沿用本问的接收保证 | - -对于一般集合无安全点的情形，最小包围圆中心确实解 $\min_q\max_{g\in K}\lVert q-g\rVert$，但其最坏距离仍可能大于1000米。改善交会角若要求偏离该中心，还会牺牲最坏距离最优性，必须明确是另一个多目标探索策略，不能把两项互相冲突的要求同时称为已保证。 - -## 13. 验证设计及当前实现边界 - -本问的验证分为数学证书与策略性能两层。数学层检查：截断角域被解析圆覆盖；安全候选到所有允许位置的距离不超过1000米；外包离散不删除真实点；角区间扩张包含所有区间内读数。数值抽样只能辅助检查这些证明，不能取代它们。 - -性能层应在相同自建源位置与确定的空间误差场上，比较轴向点、两个横向点、有限候选最坏后验策略。误差场在同一位置保持固定；不得让每次调用重新独立生成误差。记录第二次测量类型、后验直径和包围半径、额外移动、计算耗时及是否达到20米清除阈值。至少包含误差边界、源域边界、接收半径1000米、近距离和近共线案例。没有这些实验前，不报告未经验证的性能提升比例。 - -目前 [planning.py](C:/Users/HelixCopex/code/MATH/geometry/planning.py) 中的 `first_bearing_candidate_disk` 已实现第5至6节的解析证书；问题一内核可用于纯半平面区域的顶点和最小圆计算。**本仓库第二问交付的定位是"安全选点的启发式基线"**：解析安全候选盘证书已实现并证明（选中点必能收到信号），而选点评分采用边界采样点 × 读数区间端点的有限枚举近似（含 near 分支）。第9至11节描述的完整选点器数学方案，特别是第10节覆盖全部未来读数的保守区间上界，尚未在现有仓库中实现；采样最大值只用于候选点之间的相对比较与粗筛，不得当作严格最坏情形上界报告，也没有相应正式成绩。该实现边界不影响解析候选区域的证明。 - -对后续实现而言，应新增接收任意半平面约束的交集入口，或使用独立的多边形裁剪函数，再复用第一问的最小包围圆。不能把圆盘先验伪装成额外测向观测传入现有 `solve(measurements)` 接口。 - -## 14. 结论 - -第二问的候选区域可以明确取为前方约750米处、半径约250米的圆盘，其解析式为 $\mathcal Q_0=B(s+au(\theta),1000-a)$，其中 $a=1500/(2\cos\delta)$。这一候选区域对题设全向源提供接收保证，且在首次有效测向条件下严格非空。 - -在该区域内，以最坏保守后验最小包围半径为主指标、移动代价为次指标，采用有限候选和逐步加密进行选择。交会角解释局部几何优劣，外包多边形和读数区间扩张保证不低估相关后验集合；全局最优和“两次即可清除”均不作为未经证明的结论。 +严格保证与启发式性能优化在文档和代码中始终分开：接收保证有证明， +全局最优性没有证明；后者依赖自建的有限候选与固定评估口径。 ## 参考文献 -1. 全国大学生数学建模竞赛组委会. 2026年高教社杯全国大学生数学建模竞赛B题：无线电干扰源的快速自动定位与清除[Z]. 2026. [本地题面](C:/Users/HelixCopex/code/MATH/issue/B题.pdf)，问题2、附录2及第4页坐标约束。 -2. 模拟器配套材料. 附件2：模拟器通信接口说明及编程指南[Z]. 2026. [本地附件](C:/Users/HelixCopex/code/MATH/issue/附件/附件2.docx)。用于 `direction`、`near`、`no_signal` 等程序观测字段与动作接口的对应。 -3. BOYD S, VANDENBERGHE L. Convex Optimization[M]. Cambridge: Cambridge University Press, 2004. [作者公开版本](https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf)，第3章、第6.4节。用于凸性及鲁棒建模的一般背景，不作为本题候选圆盘公式的出处。 -4. 本项目. 问题一：有界测向误差下定位区域的直径与圆覆盖[Z]. 2026. [algo-1.md](C:/Users/HelixCopex/code/MATH/docs/algo-1.md)。用于角域半平面表示、顶点判定和最小包围圆的衔接。 - -**来源说明：** 第5节的截断角域最小圆及第6节的安全圆盘，是结合题设参数在本文中推导的结果；第8节的交会角关系由方位函数一阶展开推得；第9至11节是据集合包含关系构造的保守计算方案。以上陈述说明本项目的推导路径，不宣称这些一般几何思想在既有文献中首次出现。 +1. 全国大学生数学建模竞赛组委会. 2026 年高教社杯全国大学生数学建模竞赛 B 题[Z]. 2026. 本地题面 `issue/B题.pdf`。 +2. 模拟器配套材料. 附件 2：模拟器通信接口说明及编程指南[Z]. 2026. 本地附件 `issue/附件/附件2.docx`。 +3. BOYD S, VANDENBERGHE L. Convex Optimization[M]. Cambridge University Press, 2004. 凸性、最坏情形与鲁棒可行性背景。 +4. 本项目. `docs/algo-1.md`（问题一当前实现与复用差异）。 +5. 本项目. `docs/algo-3.md`、`docs/algo-4.md`（选点器在后续问题中的推广与边界）。 diff --git a/docs/algo-3-bayes.md b/docs/algo-3-bayes.md index 0db2fd7..28b1606 100644 --- a/docs/algo-3-bayes.md +++ b/docs/algo-3-bayes.md @@ -1,5 +1,11 @@ # 算法说明：问题3 贝叶斯闭环（`bayes` 策略） +> **历史补充文档.** 当前规范算法、最新默认参数与现行成绩口径见 +> [`algo-3.md`](algo-3.md)；本文保留 2026-09-11 的设计与消融细节， +> 其中 `lambda_move_s=0.0105`、站点环半径 1200 m 等数字可能已被 +> 后续回填（现行默认见 `config.BayesConfig`，λ=0.05、环半径 1159.7 m）。 +> 两者冲突时以代码与 `algo-3.md` 为准。 + 对应实现：`src/mcm26b/particles.py`（感知层）、`src/mcm26b/planner/bayes_loop.py` （决策层）、`scripts/benchmark_q3_bayes.py`（消融）、`scripts/tune_q3_bayes.py` （Optuna 调参）。基线对比：`docs/algo-3.md` 的七点批量（`seven`，A2）。 diff --git a/docs/algo-3.md b/docs/algo-3.md index c55cb29..6475546 100644 --- a/docs/algo-3.md +++ b/docs/algo-3.md @@ -1,414 +1,456 @@ -# 问题三：具有覆盖证明的多源搜索与总耗时优化 +# 问题三算法说明：全向多源搜索的覆盖证明、贝叶斯闭环与任务调度（当前实现） -## 摘要 +> 本文汇总当前仓库中问题三的全部实际算法，生产默认策略为 +> `src/mcm26b/planner/bayes_loop.py::BayesLoop`（CLI 策略名 `bayes`）。 +> `grid`（`OmnidirectionalLoop`，49 点）与 `seven` +> （`SevenBatch`，七点解析覆盖 + 本站批量）保留为可复现对照。 +> 问题三当前生产不包含任何机器学习策略；历史 HRL/RL 代码已移除， +> 自建消融中手写基线稳定优于学习策略（见 `docs/第三四问架构优化与运行提速方案.md`）。 -针对源数为10至16、频道互异、位置未知的全向无线电干扰源，建立“发现—定位—清除—完整性证明”的在线决策模型，以全部清除为约束、以移动、检测、换频道和清除的虚拟总耗时为优化目标。沿用问题一的角域交与最小包围圆、问题二的安全接收区域，将改进重点放在全域覆盖规模与多源任务调度。首先从官方演练日志识别移动与无源频道验证的主要开销；随后证明中心加六个环形检测点即可覆盖题设圆域，并提出同站批量扫频、发现任务入队、就近定位清除的可执行基线。进一步规划利用无信号观测维护剩余区域、复用顺路检测点、优化安全清除位置及滚动比较剩余任务耗时的方法。给出停止条件、几何保守性、算法复杂度和分阶段验收方案。12个自建场景的配对消融表明，七点批量策略全部清除成功，平均逐局虚拟时间降幅为58.55%；该结果用于选择下一步实现方向，不作为官方成绩或全局最优性证明。 +## 1. 问题定位与策略族 -**关键词：** 全向源搜索；覆盖证书；集合定位；任务调度；虚拟总耗时；滚动规划。 +问题三中每局有 $10\sim16$ 个全向干扰源，位置未知、频道互异（1..20）， +源位于 $\Omega=B(0,1800)$ 内；有效接收半径 $r_{\rm eff}$ 未知，但满足 +$1000\le r_{\rm eff}\le1500$ m。策略端只能看到协议响应 +（`direction` / `near` / `no_signal`），目标是发现、定位并清除全部源， +使虚拟总耗时尽量小，同时给出“不存在漏源”的证明或明确未完成。 -> **延伸实现（2026-09-11）**：本文规划的无信号负观测更新、开放路径 2-opt -> 任务排序等已落地为贝叶斯闭环策略 `bayes`（网格幸存集滤波 + 信息增益 -> 选点 + 滚动任务调度，参数经 Optuna 回填），自建消融相对本文七点批量 -> 策略平均虚拟时间再降 30.0%（均值 3726s）。设计、严格性证明与消融数据 -> 见 [algo-3-bayes.md](algo-3-bayes.md)。 +| 策略 | 类 | 发现保证 | 定位/决策 | 角色 | +| --- | --- | --- | --- | --- | +| `grid` | `OmnidirectionalLoop` | 700 m 方格 49 点，单元对角 $700\sqrt2\approx989.95\le1000$ | 发现后立即就近定位，无信息增益 | 历史基线 | +| `seven` | `SevenBatch` | 原点 + 半径 1200 m 六等分环，解析覆盖界 968.90 m | 本站先扫频，再按候选集 MEC 圆心最近服务队列 | A2 中间对照 | +| `bayes` | `BayesLoop` | 七站环（默认半径 1159.7 m，界 984.54 m）+ PF 负观测 | 网格粒子滤波 + 信息增益选点 + 滚动任务排序 | **生产默认** | -## 1. 问题界定与思路来源 +七站族还有 `SevenImmediate`（A1，发现即离站）与 `SevenSweep` +（A3，先扫完全域再统一服务）。它们的差异只在发现/定位的调度顺序， +七点覆盖证明与距离语义相同。 -第三问要求在不知道干扰源位置和具体数量的条件下，自动发现并清除全部全向源，尽量缩短总时间。每个频道最多一个源，因此本题的数据关联由频道确定；机器狗无需解决无法区分多个信号身份的问题。源位于半径1800米的圆域内，最小有效接收半径为1000米。这两个尺度之间的关系使稀疏、有证明的覆盖成为可能。[题面](C:/Users/HelixCopex/code/MATH/issue/B题.pdf)、[接口规则](C:/Users/HelixCopex/code/MATH/docs/协议规则表与状态机.md) +### 1.1 统一计时模型 -本文与前两问分工如下：问题一回答“当前观测允许源在哪里、能否一次清除”；问题二回答“下一次去哪里测才安全且有效”；问题三回答“先搜哪些位置、扫哪些频道、何时中断搜索去清除，以及何时能够确信任务结束”。全局调度会调用前两问的几何内核，但不能只把第二问的局部精度目标逐个重复就视为总时间优化。 +忽略微秒量化时，一次动作回合的虚拟耗时为 -参考本项目[问题3—4论文列表](C:/Users/HelixCopex/code/MATH/docs/B题问题3-4相关论文列表.md)，选取并核对以下相关思想。 +$$ +T=\frac{L}{5}+5N_m+N_s+5N_c^++3N_c^-, +$$ -| 来源 | 本文借鉴的内容 | 与本题的区别 | +其中 $L$ 为移动距离（速度 5 m/s），$N_m$ 为检测次数，$N_s$ 为换频次数， +$N_c^+$ 为成功清除次数，$N_c^-$ 为失败清除次数。`/clear` 本身不切换 +测向频道，`/exit` 不增加虚拟时间，访问顺序是起点固定、终点自由的**开放 +路径问题**，不是必须回到原点的 TSP 回路。官方内部按整数微秒累计，项目用 +`evaluation.recompute_virtual_time` 逐动作复算并校验。 + +## 2. 三层在线架构 + +`BayesLoop` 在 `OmnidirectionalLoop` 的会话/预算/收尾骨架上重构为三层： + +1. **保证层**：七点覆盖账本 `CoverageLedger` + 解析覆盖定理；负责 + “全测无信号 $\Rightarrow$ 该频道不存在源”。 +2. **信念层**： + - 凸候选集 `CandidateRegion`（正观测方向/近距离约束，用于清除证书）； + - 网格幸存集 `ChannelParticles`（正/负观测统一更新，用于缺席证明、 + 跳测与信息增益）； + - 两者同读一套响应流，任何一方证毕即结束该频道。 +3. **决策层**：统一任务表（站点扫频 / 定位测量 / 清除），开放路径 + 最近邻 + 2-opt 滚动重排；每执行一个任务后用真实响应重规划。 + +问题三不复用问题四的（位置 × 朝向）联合信念；`no_signal` 在全向源上 +可以排除 $B(s,1000)$。 + +## 3. 证据层：逐频道候选集与负观测 + +### 3.1 凸候选集 `CandidateRegion` + +对每个已收到过有效信号的频道 $c$ 维护凸多边形 $P_c^{\rm out}$： + +- 初始为 $\Omega$ 的 192 边外切正多边形； +- 每次 `direction`： + 1. 与 $B(s,1500)$（192 边外切多边形）求交； + 2. 与角域 $W(s,\theta,\delta_{\rm eff})$ 求交，$\delta_{\rm eff}=1^\circ+0.005^\circ$； +- `near`：与 $B(s,5)$ 外切多边形求交； +- `no_signal`：**不收缩凸候选集**。对全向源，负观测信息由 `ChannelParticles` + 处理；让凸表示保留所有正观测候选，可以避免因模型误差误删真源。 + +圆盘外切多边形的径向超出逐次累加到 `slack_m`，所有裁剪容差只向 +“保留更多候选”的方向倾斜，因此真源始终包含在 $P_c^{\rm out}$ 内。 + +清除证书使用 $P_c^{\rm out}$ 的最小包围圆 $B(m_c,\rho_c)$： + +$$ +\rho_c+\text{slack}_c+0.5\le20 +\quad\Longrightarrow\quad +\text{在 }m_c\text{ 处一次 `/clear` 必命中真源}. +$$ + +其中 $0.5$ m 是工程余量，slack 是外包多边形的半径余量；该判据是保守的 +充分条件，而不是“必须达到的最优半径”。 + +### 3.2 网格幸存集 `ChannelParticles` + +每个频道维护 $\Omega$ 外扩裕量上的方格网格，间距 +$h=14.65$ m，约 4.7 万格；Voronoi 裕量 + +$$ +m=\frac{\sqrt2}{2}h\approx10.359\ {\rm m} +$$ + +是格心到格内任一点的最大距离。所有删除判定都内缩 $m$：**只有整个格子与 +可行集不相交时才删除该格**。 + +| 观测 | 位置/距离更新 | 实现 | | --- | --- | --- | -| Bayram、Vander Hook、Isler，2016，列表第10项 | 联合规划测向位置和访问顺序，同时计入移动时间与测量时间 | 原文以GDOP描述不确定度；其近似比不能转移为本文有界角域和清除任务的近似保证 | -| Olofsson等，2020，列表第22项 | 根据已观测区域分配后续搜索资源，在有限预算下规划路径 | 原文使用PHD和期望发现数；本文用确定性的逐频道剩余区域，不用概率很小替代无源证明 | -| Perera等，2023，列表第4项 | 将稀疏覆盖路线与无线电源定位结合，并分别评价路径与采样代价 | 原文依赖RSSI；本题接口没有RSS/RSSI数值，不能直接使用其估计器 | +| `no_signal` @ $s$ | 距离 $>r_{\rm eff}\ge1000$，删除格心距 $s\le1000-m$ 的格 | `_kill_disk(s,1000-m)` | +| `direction` @ $s$，读数 $\theta$ | 距离 $\le1500$；角向偏差不超过 $\delta_{\rm eff}$ + 格内角裕量 | 先删距 $s>1500+m$ 的格，再删角向超出的格 | +| `near` @ $s$ | 距离 $\le5$ | 删格心距 $s>5+m$ 的格 | -上述内容分别见[大学研究记录及论文摘要](https://experts.umn.edu/en/publications/gathering-bearing-data-for-target-localization-6/)、[作者公开论文](https://torarnj.folk.ntnu.no/AURO2020a.pdf)、[arXiv原文入口](https://arxiv.org/abs/2312.03493)。这里只迁移任务组织方法，以下七点几何证明和题目化调度公式在本文独立推导。 +**定理（真源保留）.** 真源 $P$ 所在格的格心 $g$ 满足 +$|g-P|\le m$；与真值一致的观测不会让该格与可行集完全不相交，因此 +“整格不相交”判定永不删除该格。于是： -论文列表中建议的RSS梯度逼近不能直接用于本题：`/measure` 只返回无信号、近距离或示向度，没有连续场强值。Atanasov等的源搜索也依赖场信号、随机梯度或测量模型，其收敛结论不证明本题仅测向闭环收敛。[原论文摘要](https://arxiv.org/abs/1402.0051) 本文也不假定源位置均匀分布、每次误差独立重抽样或某一交会角在所有场景下最优。 +- **缺席证明**：幸存集为空 $\Longrightarrow$ 可行集为空 + $\Longrightarrow$ 该频道无源； +- **跳测规则**：若某测点到所有幸存格心的最小距离 $>1500+m$，则在该点 + 测量必为 `no_signal` 且不产生任何新排除，可直接跳过，省 5 秒检测和 + 1 秒换频； +- **PF 清除证书**：若幸存集 MEC 半径 $\rho_{\rm pf}$ 满足 + $\rho_{\rm pf}+m+0.5\le20$，则以该 MEC 圆心清除必命中真源。 -## 2. 现有方案的实测瓶颈 +`BayesLoop` 在每次测量后同步更新凸候选集与 PF，PF 只在频道尚未 +CLEARED/ABSENT 时接收新证据。 -官方问题3演练 `DHFK-PHD2-SR2P-FXK9` 已证明现有实现能在该局全清10个源，完整数据见[演练验证报告](C:/Users/HelixCopex/code/MATH/output/q3/practice/20260911-port2026/report.md)。该局在策略运行结束后才读取官方源数，源真值没有进入决策。 +### 3.3 数量上界证明 -| 分项 | 本局数值 | 对改进的含义 | -| --- | ---: | --- | -| 总虚拟时间 | 16217.454814秒 | 作为官方基线，不与自建场景混算 | -| 移动距离 / 耗时 | 64737.274米 / 12947.454814秒 | 移动约占79.84%，应优先减少站点和往返 | -| 检测次数 / 耗时 | 540次 / 2700秒 | 无源频道各扫49点，共490次，占全部检测90.74% | -| 换频道耗时 | 520秒 | 可以优化，但节省一次换频仅相当于少走5米 | -| 清除耗时 | 50秒，10次全部成功 | 清除可靠性已得到单局支持，无需先追求更小定位半径 | -| 最后一个源清除时刻 | 8454.332399秒 | 其后7763.122415秒用于剩余任务，约占总时间47.87% | +题面给定每局源总数 $10\sim16$ 且每频道至多一个源。若 +`DETECTED`/`CLEARED` 的不同频道数达到 16，则所有仍未出现有效信号的 +`UNKNOWN` 频道都不可能有源，可直接标记为 `ABSENT`，证据来源记为 +`count_bound`。这条规则与位置无关，是严格计数证明，不需要等七站扫完。 +`summary.completion_basis` 会区分 `count_cap` 与 `per_channel_proof`， +避免把“达到数量上界”误读成“逐频道缺席证明”。 -现有 `OmnidirectionalLoop._coverage_phase` 的结构还引入额外往返：在一个覆盖站按频道扫描，遇到 `direction` 就立即定位清除，再回原覆盖站继续下一频道。由原始动作日志统计，从非网格点返回网格检测点共10次，累计约9996.714米、1999.343秒。这是实际行程分类，不能全部当成可直接删除的收益，因为改序后仍需访问下一任务；但它明确说明同站任务应先批量处理。 +## 4. 发现保证：七点解析覆盖与覆盖账本 -当前49点网格同时服务第三问和第四问。它的定向覆盖证明要求更强的几何条件；第三问无需为未知发射朝向付出同样的采样代价。改进第三问时，应显式区分全向策略和第四问策略，不能全局替换公用网格。 +### 4.1 `full_disk` 语义 -## 3. 状态、信息与优化目标 +全向源在测点 $s$ 返回 `no_signal`，意味着真实距离严格大于 +$r_{\rm eff}\ge1000$。因此若集合 $S$ 满足 -### 3.1 符号与可用信息 +$$ +\Omega\subseteq\bigcup_{s\in S}B(s,1000), +$$ -| 符号 | 含义 | +则任何合法源必在某个 $s\in S$ 的 1000 m 内，在该点不可能返回 +`no_signal`。于是“某频道在全部 $S$ 点都测得 `no_signal`”构成严格的 +缺席证明。 + +`CoverageLedger` 按 `(channel, point_index)` 二维记录，只有所有点都已测 +（或显式跳过）才允许 `coverage_complete(channel)=True`；“路过某点”不算 +“测过某频道”。 + +### 4.2 七点构造与解析界 + +取原点加半径 $a$ 的六等分环： + +$$ +S(a)=\{(0,0)\}\cup\left\{a(\cos(k\pi/3),\sin(k\pi/3)):\ k=0,\dots,5\right\}. +$$ + +**定理 1.** 当 $01000$ 时退化为探索 MEC 圆心，可能 `no_signal`）。盘内取 +37 个候选点（圆心 + 3 环 × 12 方向），以“边界见证点 × 读数区间端点 + +near 分支”的最坏后验 MEC 半径为主指标、移动时间为次指标；选出后再按 +1/4 间距加密复核。该函数最初为问题二实现，在问题三中服务于 +`SevenBatch`、`OmnidirectionalLoop` 的定位子循环，以及 `BayesLoop` +的阻塞式安全后备 `_localize_and_clear`。 + +### 5.2 `bayes` 的信息增益选点 + +`BayesLoop._select_point` 在 PF 幸存集上选择定位测点，候选集分两种： + +- 凸候选集安全盘非空：取 `select_measurement_point` 的环形候选，并额外 + 加入以当前位置到 MEC 圆心方向为轴的对称候选，以及零移动的当前位置； +- $\rho>1000$、无保证安全盘：取 PF MEC 的 $1.0/1.5/2.0$ 倍探索环和 + PF MEC 圆心。 + +对候选 $q$，默认评分模式 `entropy` 是期望对数幸存占比： + +$$ +\text{score}(q)= +\mathbb E[\log(\text{后验幸存占比})] ++\lambda\left(\frac{\lVert p-q\rVert}{5} ++\frac{\lVert q-m_{\rm MEC}\rVert}{5}\right), +\qquad \lambda=0.05. +$$ + +其中期望对代表假设（从幸存格均匀抽样）和 $r_{\rm eff}\sim U(1000,1500)$ +解析积分： + +- `no_signal` 分支（概率 $(d-1000)/500$）：幸存格收缩到 + $B(q,1000-m)$ 之外； +- `direction` 分支：以代表点方位为期望读数，用 $\pm\delta_{\rm eff}$ + 角域裁剪后的子采样估计幸存占比； +- `near` 分支（$d\le5+m$）：后验为 5 m 盘，占比趋近零。 + +可选 `score_mode=radius` 直接最小化期望后验 MEC 半径。 +`score_subsample=800`、`mc_representatives=7` 是工程参数；均匀源分布、 +零误差代表读数和 $r_{\rm eff}$ 均匀分布都只是**选点评分假设**，不参与 +任何安全证明。安全证明只使用题面给定的 $\pm1.005^\circ$ 角域与距离界。 + +当前配置中以下“更强的耦合”均已实现但实测为净负收益，默认关闭： +`fuse_localize_into_tour=False`（站点顺带定位）、`exact_order_limit=0` +（Held-Karp 精确任务序）、`route_aware_probe=False`（路线感知探点）。 +代码与开关保留，便于复现与后续研究。 + +### 5.3 清除动作 + +`BayesLoop._best_certificate` 的优先级： + +1. **凸证书**：`clear_certificate(region)`； +2. **PF 证书**（`grid_cert_enabled=True`）：负观测挖洞使 PF 幸存集更紧时， + 若 PF MEC 半径 $+m+0.5\le20$，用 PF MEC 圆心作为清除点。 + +得到证书后，`_near_clear_target` 还会把清除点从证书圆心向当前位置方向 +平移一小段：设当前证书半径 $r$、限值 $limit$，最多平移 +$\Delta=\min(\text{当前距离},limit-r-10^{-6})$，新目标点仍保证所有候选 +源在 $r+\Delta\le limit$ 内。该平移可省去最后一段接近路程，不放松证书。 + +定位子循环每频道最多 `max_localize_steps=8` 次追加测量；若在安全接收盘内 +出现 `no_signal`，记录异常 `no_signal_in_safe_disk`，不静默忽略。 + +## 6. 决策层：统一任务表与滚动重排 + +### 6.1 任务类型 + +`BayesLoop._build_tasks` 在每次决策点构建： + +| 任务 | 位置 | 服务时间 | 触发条件 | +| --- | --- | --- | --- | +| `station` | 七站之一 | $5\times$待测频道数 + 换频 | 该站仍有 `UNKNOWN` 频道未测，且跳测规则未排除 | +| `locate` | 信息增益选点 | 5 s + 可能的换频 | 频道 `DETECTED` 且暂无清除证书，且未委托给站点 | +| `clear` | 证书中心（可平移） | 5 s 成功（失败 3 s） | 频道有清除证书且未处于“清除失败后待补测”状态 | + +站点任务的服务时间计算当前频道优先：若机器人当前测向频道包含在待测集合， +则少一次换频。清除动作不改变测向频道。 + +### 6.2 跳测规则 + +`guaranteed_no_signal(station)`：若该频道所有幸存格心到站点的距离都 +大于 $1500+m$，则该次测量必为 `no_signal` 且不会产生新排除，任务构造时 +直接跳过。跳测缓存按 PF `version` 生效；幸存集单调收缩，因此判定单调。 +这条规则是负观测给问题三带来的主要检测节省之一。 + +### 6.3 开放路径排序 + +任务表规模通常为 6～12 项，代价包含移动、检测、换频，终点自由。默认： + +1. 最近邻构造初始序； +2. 开放路径 2-opt 改进，最多 `two_opt_passes=4` 轮； +3. 换频代价按实际频道状态沿路径逐步传递，而不是固定服务费； +4. 执行排序后的**第一个任务**，拿到响应后立即回到决策点重新构建/排序 + （滚动时域）。 + +可选的 Held-Karp 精确 DP（含“已访问集合 × 结束时频道”状态）与 +`route_aware_probe` 路线感知精修均已实现；在当前冻结基准集上， +精确 DP 与 2-opt 差距在噪声内（甚至略差），默认关闭。 + +### 6.4 已停靠位置的机会测量 + +`_measure_at_stop` 在每次决策点评估当前位置是否值得顺带补测已发现频道： +只对尚未在该 0.5 m 点测过、无证书、到凸候选集所有顶点距离均不超过 +$1000-10^{-3}$ m、且预测后验半径至少收缩 +`opportunistic_radius_ratio=0.5` 的频道追加一次测量。每次仍支付 5 秒与 +实际换频成本，只是省去移动；它不是“免费检测”。 + +## 7. 停止、预算与异常收尾 + +### 7.1 完成证明 + +循环停止后，`summary.completed` 的依据分两类： + +- `count_cap`：清除数达到 16（题面源数上界）； +- `per_channel_proof`：每个频道状态都是 `CLEARED` 或 `ABSENT`。 + +缺席证据在 `summary.absence_evidence` 中逐频道记录来源： + +| 证据 | 条件 | 严格性 | +| --- | --- | --- | +| `pf_empty` | 某频道 PF 幸存集为空 | 严格（Voronoi 真源保留） | +| `coverage` | 七点/49 点全部测得 `no_signal`（或计划上显式跳过） | 严格（解析/校核覆盖） | +| `count_bound` | 已有 16 个不同频道的有效存在证据 | 严格（题面数量上界） | + +`persistent_search` **默认 False**：一旦“每频道已清除或证明不存在”就立即 +`/exit`，因为虚拟时间是成绩口径。显式 `--persistent` 会继续复扫未清除 +频道直到清除数达 16、实机预算耗尽或虚拟安全余量不足；历史官方演练证明 +无条件复扫会浪费大量时间，因此该项目只作为可审计的保险开关保留。 + +### 7.2 预算与中断 + +- 虚拟安全余量 `virtual_reserve_s=10000` s：剩余虚拟时间不足以完成 + 一个保守单动作上界时停止发起新动作，避免触及 100 小时上限； +- 现实预算以服务器 `/enter` 返回的 `remaining_real_duration_s` 为准； + `DeadlineExceeded` 时 stop_reason 记为 `real_budget_exhausted`； +- 传输、协议或守恒性异常都会尽力 `/exit` 并输出未完成 summary， + 不静默丢会话；`completed=False` 与具体 `stop_reason` 同时留档。 + +## 8. 基线与自建消融 + +### 8.1 调度变体 + +| 版本 | 发现/定位关系 | 说明 | +| --- | --- | --- | +| A0 `grid` | 49 点，发现后立即离站定位 | 历史生产基线 | +| A1 `SevenImmediate` | 七点，发现后立即离站定位 | 只替换覆盖点 | +| A2 `SevenBatch` | 七点，本站先扫完再按最近候选圆心服务队列 | 中间对照，曾用于官方演练 | +| A3 `SevenSweep` | 七点，先扫完全部必要站点再统一服务 | 批量调度的另一极端 | + +A2 相对 A0 的收益主要来自“减少站间往返 + 每站批量扫频 + 无源频道只需 +7 次而不是 49 次检测”。在早期 12 个自建场景上，A0/A1/A2/A3 平均虚拟 +时间分别为 13788.04 / 7476.09 / 5496.20 / 6094.73 秒，A2 逐局平均降幅 +58.55%；更晚的 14 场景消融中，`bayes` 相对 A2 再降 30.0% +（3726.2 s vs 5352.7 s）。这些数字来自自建模拟器与固定误差场，不是官方成绩。 + +### 8.2 冻结基准集当前结果 + +`project/q3` 在 18 个冻结 core 场景（seeds 100–105 × 10/13/16 源， +`harness/suites.py` 版本化清单）上的配对结果： + +| 方法 | 有效局 | 全清率 | 虚拟时间均值 | 中位 | 最大 | 检测均值 | 请求均值 | 位移均值 | +| --- | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | +| `bayes` | 18/18 | 100% | **3188.7 s** | 3111.3 s | 3663.7 s | 111.1 | 124.1 | 12338 m | +| `seven` | 18/18 | 100% | 5361.7 s | 5367.5 s | 6462.4 s | 100.8 | 113.8 | 23508 m | +| `grid` | 18/18 | 100% | 14935.7 s | 15275.1 s | 16691.4 s | 430.6 | 443.6 | 61550 m | + +配对变化：`bayes` 相对 `grid` 虚拟时间 **−78.4%**，18/18 场景更快； +所有 54 局零缺席误判、零异常、正常退出且计时复算一致。另有边界场景 +（源在 $\Omega$ 边界、$r_{\rm eff}=1000$）验证七点覆盖的扇区接缝。 + +### 8.3 官方演练记录（历史口径） + +- 旧 49 点实现通过官方问题三演练案例 `DHFK-PHD2-SR2P-FXK9`： + 10/10 全清、540 次检测、10 次清除，虚拟总时间 16217.45 s， + 现实 7.629 s； +- A2 七点批量通过官方案例 `J4ER-Y9EP-Q379-UPFA`：12/12 全清、 + 109 次检测、虚拟 5232.13 s、现实 3.406 s； +- 后续一次官方问题三演练（`3BW5-9R54-NZNA-MYB9`）记录 13/13 全清、 + 虚拟 4058.71 s。 + +以上案例的源场景、实现版本、误差场均不同，不能直接横向比较，也不能把 +历史官方数字当成当前 `bayes` 同案例重跑结果。当前默认 `bayes` 的官方向 +成绩仍需新的演练验证。 + +## 9. 复杂度、参数与保守性边界 + +| 模块 | 计算量（近似） | | --- | --- | -| $\Omega=B(0,R)$，$R=1800$ | 源位置区域，单位米 | -| $R_-=1000$，$R_+=1500$ | 有效接收半径下界、上界 | -| $p,f,t$ | 当前机器人位置、测向频道、官方虚拟时刻 | -| $K_c$ | 已发现频道 $c$ 的保守位置候选集 | -| $E_c$ | 尚未发现频道 $c$ 仍可能存在源的位置集合 | -| $\mathcal U,\mathcal D,\mathcal L,\mathcal C,\mathcal A$ | 未发现、待定位、可保证清除、已清除、已证明无源频道集 | -| $\mathcal P$ | 曾获得有效存在证据的不同频道集，包含已清除频道 | -| $\delta=1.005^\circ$ | 题设1°误差加两位小数舍入裕量 | -| $S$ | 具有全域覆盖证明的检测站集合 | -| $L,N_m,N_s,N_c^+,N_c^-$ | 移动长度、检测数、换频数、成功与失败清除数 | - -源静止、同位置同频道的测向误差固定。只有该频道实际收到的观测才能更新其状态；经过一个位置不代表在那里扫过所有频道。不同频道的无信号证据不能互相替代。 - -### 3.2 统一的时间口径 - -忽略微秒量化时，整局虚拟耗时为 - -$$T=\frac{L}{5}+5N_m+N_s+5N_c^++3N_c^-.$$ - -`/clear` 的目标频道不改变测向机频道，因而不引入换频代价；`/exit` 不要求回到原点，也不增加虚拟时间。因此访问顺序是从当前位置出发、终点自由的**开放路径**问题，不能默认使用需要返回起点的旅行商回路。 - -官方内部按微秒累计，验收时按每次移动耗时量化到整数微秒，再累加检测、换频和清除的整数秒。连续实数公式用于规划，逐动作微秒公式用于复算。原演练汇总的541处偏差已查明为累计浮点与微秒精度口径不同，详见[独立审计](C:/Users/HelixCopex/code/MATH/output/q3/practice/20260911-port2026/verification.json)。 - -优化目标按优先级定义为：保持真源候选和有效完成条件；在现实1200秒、虚拟360000秒上限内全清；在可完成策略中降低 $T$。现实计算和网络耗时单独约束，不加进虚拟时间公式。由于源场景和误差场未知，理想模型可写为对允许场景集合 $\Xi$ 的 $\inf_\pi\sup_{\xi\in\Xi}T(\pi,\xi)$，但本文不声称解出了这个连续在线极小极大问题，也不声称已证明任意合法场景都能在两个预算内完成。 - -## 4. 逐频道证据与安全停止 - -### 4.1 观测更新 - -首次 `direction` 在 $s$ 返回 $\theta$ 时，建立 - -$$K_c=\Omega\cap B(s,R_+)\cap W(s,\theta,\delta).$$ - -以后每次有效测向再交入相应距离上界和角域。`near` 时交入 $B(s,5)$，当前检测点本身就满足20米清除条件。角域的半平面表达、保守圆盘外包和最小包围圆计算沿用[问题一](C:/Users/HelixCopex/code/MATH/docs/algo-1.md)、[问题二](C:/Users/HelixCopex/code/MATH/docs/algo-2.md)。 - -全向源在 $s$ 返回 `no_signal`，意味着若该频道仍存在未清除源，则其距离严格大于实际接收半径，进而大于1000米。因此对未发现频道可维护 - -$$E_c=\Omega\setminus\bigcup_{s\in\mathcal N_c}B(s,R_-),$$ - -其中 $\mathcal N_c$ 只含该频道实际测得无信号的位置。有限精度实现可只排除 $B(s,R_- -\eta)$，$\eta>0$ 是有依据的数值裕量；少删不会漏源。该负观测更新仅用于第三问，不能直接用于第四问的定向源。 - -圆盘差通常非凸。不能把扣除圆盘后的集合强行表示成单个凸多边形，更不能丢弃狭窄或小面积分量。第一阶段保留当前凸定位内核，仅在未知频道的覆盖账本中使用负观测；后续如把无信号用于已发现源的精定位，应使用多片区域或再取保守凸外包。 - -### 4.2 完成与计数推理 - -安全停止条件是 $\mathcal C\cup\mathcal A=\{1,\ldots,20\}$，且本局 `/exit` 已被接受。无源证明可来自完整覆盖、剩余集合的可靠判空，或下述数量上界。 - -由于每频道至多一个源且总数不超过16，当 $|\mathcal P|=16$ 时，所有未出现过有效信号的频道即可由数量上界判为无源，**不必等待16个都清除才停止扫未知频道**。随后仍须清除 $\mathcal P\setminus\mathcal C$，不能直接退出。成功清除16个则立即满足全清条件。 - -只发现或清除10个不够，因为10是下界。反过来，若无源证据已排除10个频道，则其余10个频道都必须有源，但仍需获取其定位观测；这种计数推理不能替代位置搜索。若已知存在与无源证据互相矛盾，或可存在频道数不足10，应报告协议、模型或数值不一致，不能把矛盾解释为成功。 - -## 5. 用七个检测站替代49点网格 - -### 5.1 布点 - -取原点及半径 $a$ 圆周上的六个等角点： - -$$S(a)=\{(0,0)\}\cup\{a(\cos(k\pi/3),\sin(k\pi/3)):k=0,\ldots,5\}.$$ - -推荐首个工程版本固定 $a=1200$ 米。站点坐标为 $(0,0)$、$(1200,0)$、$(600,\pm600\sqrt3)$、$(-600,\pm600\sqrt3)$、$(-1200,0)$。这些点全部位于源圆域内；允许域外检测这一规则仍然保留，但本方案不需要依赖它。 - -### 5.2 连续覆盖证明 - -**定理1：** 若 $01800$ 的单元才可直接判为域外；边界相交不确定时继续保留或细分。 - -一次该频道无信号观测位于 $s$，若整个单元满足 - -$$\|z-s\|+h\le1000-\eta,$$ - -则单元全部落入保证可收盘，可以安全删除。不得只检查单元中心，也不得用近似覆盖面积达到100%作为判空证据。若单元仍与真实残余相交，就保留它，直至可给出整单元删除证书。 - -**命题2：** 以上规则维持 $E_c\subseteq\widetilde E_c$，故 $\widetilde E_c=\varnothing$ 足以证明无源。反过来，$\widetilde E_c\ne\varnothing$ 可能只是离散保守性，不能据此断言有源。 - -七点有31.098米几何余量；可使最终方格边长不超过20米，其直径不超过28.285米。对任意与 $\Omega$ 相交的单元，取其中一个域内点，该点距某站不超过968.902米；单元全部点距该站不超过997.186米。因此在工程裕量小于约2.81米且正确处理边界相交的情况下，七站无信号足以最终删除这些单元。这里用整单元的最大距离证书；若只用 $\|z-s\|+h$ 的较松证书仍未通过，继续细分或直接使用七点账本证明，不错误地停在粗格残片上。 - -### 7.2 候选点与频道选择 - -新增候选点来自：未完成的七个锚点、当前位置、已规划的定位点/清除点、残余区域边界附近点。对 $q$ 和频道集合 $F$，定义 - -$$G(q,F)=\sum_{c\in F}\operatorname{area}\bigl(\widetilde E_c\cap B(q,1000-\eta)\bigr),$$ - -并以 $G(q,F)/(\|p-q\|/5+5|F|+N_{\rm switch}(F))$ 作为探索排序指标。面积是搜索潜力代理，不是发现概率或保证的信息增益；面积计算可近似，但只影响候选排序。无源判定仍由整区域覆盖证书给出。对极小、线状或退化残余，还要保留单元数、最远未覆盖距离及锚点完成优先级，不能因面积接近零而永久忽略。 - -在当前位置按每个候选频道的 $G(q,\{c\})$ 与5秒检测、可能1秒换频排序，只有优于当前预演方案时才添加。未来站点若要删除，必须重新证明该频道的剩余 $\widetilde E_c$ 被保留站点的保证接收盘覆盖；一次扫描所得“没有发现”只删该次盘内区域，不能整片按经验跳过。 - -### 7.3 完备后备 - -七点中未完成的“站点—频道”记录始终保留。自适应探索若连续不产生新的可靠覆盖、计算超时或动态证书不确定，则恢复固定七点访问。设置有限的额外探测额度或要求每次替换都完成至少一个新的覆盖记录，使启发式不能无限拖延覆盖完成。 - -这给出确定性的搜索后备：自适应部分只允许节省已证明冗余的动作，不允许删除尚未证明冗余的覆盖责任。该设计尚未进入当前生产实现，也未包含在第12节的七点消融中。 - -## 8. 定位和清除也应优化剩余时间 - -### 8.1 定位到“够清除”就停止 - -设保守多边形 $\widetilde K_c$ 的最小包围圆为 $B(m_c,\rho_c)$。清除阈值使用 $r_{\rm safe}=20-\varepsilon$，其中 $\varepsilon$ 覆盖实际计算和坐标提交误差。若 $\rho_c\le r_{\rm safe}$，就具有一次清除证书；没有必要为了获得更小半径继续增加检测或横移。 - -外包多边形本身已经包含真集，它的多余面积使决策更保守；不能把“外包多边形与原圆盘的径向差”未经传播分析当作浮点误差保证。现有实现还扣除了 `slack_m` 和工程清除余量，本次消融完整保留其行为；后续统一外包与数值误差口径时另行验收。 - -### 8.2 在所有保证清除点中选择短路径位置 - -若 $V_c$ 为 $\widetilde K_c$ 的顶点，则 - -$$\mathcal C_c=\bigcap_{v\in V_c}B(v,r_{\rm safe})$$ - -是凸集，非空当且仅当 $\rho(\widetilde K_c)\le r_{\rm safe}$。圆盘凸性保证覆盖全部顶点等价于覆盖整个多边形，故任意 $q\in\mathcal C_c$ 都能清除真源。 - -无后续任务时，选择 $\arg\min_{q\in\mathcal C_c}\|q-p\|$；已计划下一站 $b$ 时选择 $\arg\min_{q\in\mathcal C_c}(\|p-q\|+\|q-b\|)$。两者都是凸优化问题，可用二阶锥形式表达，也可在圆盘交的候选边界上求解后逐顶点复核。 - -轻量实现可先使用内部安全圆盘 $B(m_c,r_{\rm safe}-\rho_c)\subseteq\mathcal C_c$：将 $p$ 投影到该圆盘即可得到一个不远于前往 $m_c$ 的保证清除点。这只是对该充分圆盘的最近点，并非对整个 $\mathcal C_c$ 最优。若圆盘半径为0就退化到 $m_c$。 - -### 8.3 安全测点的目标由半径改为任务耗时 - -问题二提供的安全区域为 $\mathcal Q(K_c)=\bigcap_{g\in K_c}B(g,1000)$，可用解析盘或保守多边形构造充分区域。选点仍必须先满足接收保证，不能为走近一个猜测坐标而接受未经评估的失联风险。 - -在可安全测量的候选点中，使用以秒计量的短程代价 - -$$J_H(q)=\frac{\|p-q\|}{5}+5+\mathbf1_{f\ne c}+\max_{z\in\mathcal Z(q)}\widehat V_{H-1}(K_c\mid q,z;q,c),$$ - -其中 $z$ 包含 `near` 及允许的示向度，$\widehat V$ 估计完成该源及返回后续任务的时间。若某分支已能保证清除，其终端代价包括移动至 $\mathcal C_c$ 的时间和5秒清除；若未能清除，可继续一层安全测向，或赋予经过标定的剩余定位代价。未来测点、清除位置与访问顺序共同影响 $\widehat V$。 - -这个形式解释了为什么“后验半径更小”不一定“完成更快”：走远很久将半径降到5米，可能劣于近处测量后已经达到19米。所有实际比较项统一为秒，避免随意相加米、面积和时间。有限候选和有限视野只给出预演策略的启发式排序。 - -若要把 $\max_z$ 报告为严格上界，必须采用问题二第10节覆盖全部读数的角区间扩张，并对区域误差与终端代价也作保守上界。当前“边界位置采样×误差端点”并无该保证，不能仅因取了样本最大值就声称最坏总时间已受控。推荐先保留现有局部选点，验证搜索调度收益，再单独启用这一成本模型。 - -## 9. 保守性、有限结束与预算边界 - -**命题3（不漏发现）：** 固定七点对每个仍未知频道完成扫描，或由可靠的残余覆盖/数量证书替代相应扫描，必能发现所有合法全向源。 - -**证明：** 固定点由定理1保证；动态替代保持该频道覆盖域包含所有仍可能位置；数量证书只在已有16个不同源的存在证据时排除其他频道。因此真源不可能被当作无源跳过。证毕。 - -**命题4（证书清除）：** 若维护的外包集合包含真源，实际提交点满足对该外包的20米清除约束，则在题设规则下 `/clear` 必成功。几何证书通过而官方失败时必须记录异常并检查模型、数值及协议状态，不能在同一证书上循环重试。 - -上述命题不证明任意安全测点序列都会有限步收缩到20米。现有每轮8次测量是工程限额，不能作为收敛定理。拟在后续增加有限后备：对已发现频道的有界外包区域，用边长20米方格的中心覆盖所有相交格，任意格内点至格心不超过 $10\sqrt2<20$ 米；遍历该有限清除候选即可在模型成立时清除源。单元相交计算必须向保留更多格的方向保守。 - -后备仅在定位停滞时使用，并须先检查实际枚举格数、路径和预算。设当前到首格及格间的计划路长为 $L_f$，共 $M_f$ 格，最迟第 $M_f$ 格成功的虚拟耗时上界为 - -$$T_f\le L_f/5+3(M_f-1)+5.$$ - -使用相交格数硬顶、现实高分位单请求估计 $M_f\,t_{\rm RTT}^{95}$ 加退出预留，以及虚拟剩余预算共同筛选。现有第四问中的后备格算法可供复用，但不能跳过第三问的验证。硬顶或预算不允许时，保留“未完成”，不将耗尽预算说成证明无源。 - -因此本文区分三层结论：覆盖与证书的数学安全性；预算足够且正确遍历有限后备时的有限完成性；现实/虚拟双预算内的性能，目前由实验评价，尚无所有合法场景的统一完成上界。 - -## 10. 推荐算法与伪代码 - -**算法3A：七点覆盖＋同站批量扫频＋就近定位清除（优先落地版本）。** - -```text -读取第三问参数，生成七点 S(1200)，初始化逐频道观测和覆盖账本 -调用 /enter，记录官方时间预算 -while 尚有未知频道的站点记录未完成： - 若预算不足：标记未完成，转退出 - 选择仍有未知频道待测的最近站点 s - 在 s 批量测量待查频道（当前测向频道若在其中，优先测） - no_signal：记录该频道在 s 的覆盖证据 - direction：更新该频道 K，加入待定位队列 - near：在 s 清除并依据响应更新状态 - 对本站结束后的待定位队列： - 按当前位置到候选圆心距离选最近频道 - 调用问题二安全选点与问题一清除证书完成其定位子循环 - 清除成功后重新排序；未完成者保留，单轮不无限重试 - 若已清除16个：结束任务 -根据逐频道七点账本标记无源；继续处理尚未清除的已发现频道 -仅在全频道均已清除/有无源证明时标记完成 -调用 /exit；退出响应失败单独记录，不能宣称正常结束 -``` - -**算法3B：在3A之上逐项增加的完整改进规划。** 每次 `direction/near` 更新存在集合 $\mathcal P$，达到16时提前结束未知频道搜索；在各实际测量点更新负观测区域并进行有证书的站点替代；用插入耗时比较“先搜下一站”与“先做某个源”；对候选路线重排；以预计剩余时间选择安全测点，在可清除区域中选择省路位置；定位停滞时按预算调用有限后备。每执行一个动作就用真实响应更新状态，禁止整段盲执行预测的后续动作。 - -```mermaid -flowchart TD - A[进入并读取预算] --> B[更新逐频道证据与完成条件] - B -->|已全清且无遗漏证明完整| X[退出并保存日志] - B -->|预算不足| Y[按未完成退出] - B -->|仍有任务| C[候选覆盖站与定位清除任务] - C --> D[比较短程总耗时并保持覆盖责任] - D --> E[同站批量扫描或执行一次定位清除动作] - E --> F[接收官方响应并验证候选及证书] - F --> B -``` - -第一版和后续规划分别命名，防止把3B尚未实现的负观测跳站、时间型选点和后备算法算进3A的测试成绩。 - -## 11. 计算量与工程实现 - -七点固定账本最多 $7\times20=140$ 个布尔记录。最近站选择检查7个点，每次按至多20个频道统计待测项；16个源的贪心任务队列重排规模很小。固定站路线的开放路径最优排序如需精确解，可使用子集动态规划，规模为 $O(k^2 2^k)$、$k\le7$；但加入观测依赖的定位任务后不再是同一个静态问题。 - -定位计算仍主要消耗在多个候选点、候选源见证点、读数分支的几何裁剪和最小圆计算。假设候选测点 $Q$ 个、读数区间 $B$ 个、外包顶点 $V$ 个，单层保守评估包含 $QB$ 次裁剪和最小圆计算，可表示为 $O(QB[V+T_{\rm MEC}(V)])$；两层分支不能简单称为线性。复用候选集缓存、先粗筛少量入围点，再完整区间验证，能够控制现实计算预算。 - -残余区域若使用 $H$ 个单元，20频道掩码空间为 $O(20H)$，一次频道负观测朴素检查 $O(H)$；可以预计算七锚点的单元覆盖掩码，并通过空间索引限制动态点的更新范围。分辨率只控制证书的保守程度，不得因内存不足删除尚未证实的残余格。 - -| 模块 | 当前事实 | 建议改动 | -| --- | --- | --- | -| `coverage` | 共用49点与逐频道账本 | 新增第三问七点构造；保留第四问网格及其独立证明 | -| `planner.OmnidirectionalLoop` | 发现后立即离站定位 | 把“本站扫频”与“队列服务”拆开，引入第三问策略选项 | -| `belief` | 无信号不裁剪候选 | 增加第三问独立残余覆盖层，避免破坏第四问朝向语义 | -| `geometry/planning`、`planner` | 有安全选点与最小圆证书 | 后续增加最近安全清除点和时间型候选评分 | -| `evaluation` | 浮点累计、固定微秒容差 | 单独兼容官方逐步微秒时钟和自建连续时钟，保留两种原始口径 | -| `protocol` | 串行、幂等、预算退出、JSONL日志 | 延用；增加的动作仍走同一客户端，不能并发测多个频道 | - -实际路径分别位于 [coverage](C:/Users/HelixCopex/code/MATH/src/mcm26b/coverage/__init__.py)、[planner](C:/Users/HelixCopex/code/MATH/src/mcm26b/planner/__init__.py)、[belief](C:/Users/HelixCopex/code/MATH/src/mcm26b/belief/__init__.py)、[evaluation](C:/Users/HelixCopex/code/MATH/src/mcm26b/evaluation/__init__.py)。前两问文档中部分旧路径仍写为根目录 `geometry`，复用时以当前 `src/mcm26b/geometry` 为准。 - -## 12. 已完成的自建配对消融 - -### 12.1 实验口径 - -为验证优先级，新增独立脚本 [benchmark_q3_search.py](C:/Users/HelixCopex/code/MATH/scripts/benchmark_q3_search.py)，只连接进程内 `LocalSimulator`，不连接官方端口、不修改生产策略。数据保存在 [search-ablation](C:/Users/HelixCopex/code/MATH/output/q3/search-ablation/benchmark.json)，同时逐局保存源场景、固定误差场种子、接受动作日志与哈希。 - -共12个场景：源数10和16各一组，每组包含 `random_omni_scene` 种子0、1、2、3、4的5个场景，以及1个源位于1800米边界、接收半径均为1000米的场景。随机生成器径向采样不是圆域面积均匀分布，并特意含近边界源，不将其称为真实场景分布。各变体使用同一局源坐标、接收半径、频道和固定空间误差场；不同位置的合成误差机制不代表官方误差场。边界例包含扇区接缝附近的源。 - -比较四种版本，均沿用现有安全选点与清除内核、舍入裕量0.005°和8次追加检测限额： - -| 版本 | 唯一涉及的搜索与调度变化 | +| PF 更新 | 与删除格/角向扫描的幸存格数相关；粗格索引与向量化（numpy 可选）控制常数 | +| 七点覆盖 | 最多 $7\times20=140$ 个 `(站, 频道)` 记录 | +| 任务排序 | 最近邻 $O(n^2)$ + 2-opt $O(n^2)$/轮；精确 DP 为 $O(2^nn\cdot|C|)$，默认关闭 | +| 凸候选集 | 192 边形外切 + $O(m)$ Sutherland–Hodgman 裁剪；MEC 用 Welzl | +| 选点评分 | $O(N_q\times(\text{子采样}))$；候选约 37 个，子采样 800 | + +保守性边界： + +- 所有“证毕”结论依赖题面观测语义和保守外包/内缩；若接口语义、 + 接收半径范围或误差界改变，必须重新校核； +- 性能结论依赖自建场景分布、固定误差场和排序启发式，不构成对全部合法 + 场景的最坏时间保证； +- PF 的均匀先验与信息增益评分是选点启发式，不是官方概率模型； +- `persistent_search` 关闭时，如果官方场景违反题面数量上界或存在模型外 + 行为，不能由本算法保证发现；官方接口语义变化时必须记录异常而不是 + 强行把未完成解释为成功。 + +## 10. 与问题一、二、四当前实现的复用关系和差异（核对结果） + +| 复用/差异 | 具体情况 | | --- | --- | -| A0：49点即时定位 | 当前生产基线 | -| A1：七点即时定位 | 只把覆盖点换成七点，仍发现即离站 | -| A2：七点本站批量 | 本站扫频完成后，按最近候选圆心服务队列；当前频道优先 | -| A3：七点全域先扫 | 扫完所有必要覆盖站后，再按最近候选圆心定位；near仍原地清除 | +| 问题一几何内核 | `belief/__init__.py` 与 `particles.py` 只从 `geometry.localization` 导入 `smallest_enclosing_circle`；在线 `clip_wedge` 复刻同一半平面公式，但工作在凸多边形裁剪流程，不调用 `solve` | +| 问题二选点器 | `select_measurement_point` 直接从问题二推广而来，问题三用它做基线定位与安全后备；`safe_disk_radius` 同样复用 | +| 问题一没有的 | 负观测 `no_signal`、PF 非凸挖洞、覆盖缺席证明、数量上界、任务排序、滚动时域、预算/异常治理 | +| 问题二没有的 | 多频道数据关联、发现覆盖、清除动作、在线顺序决策；问题二只做一次第二点选择 | +| 问题四的差异 | 问题四不能用七点 `full_disk` 缺席证明；`no_signal` 可能由朝向遮挡造成，因此改用 `forward_half_disk` 31 点覆盖 +（位置 × 朝向）联合信念 + 预算化后备清除（见 `algo-4.md`） | -A2/A3都包含当前频道优先这一小优化，因此A1与A2之差不是“纯粹只改变一个排序变量”的严格单因素实验。后续应再拆分消融换频顺序。3B的动态残余区域、16个已发现源提前停止搜索、时间型选点、最近安全清除点及后备均未在这组实验启用。 +一句话概括：**问题三 = 问题一的保守几何 + 问题二的保证接收选点 + 逐频道非凸负观测 + 覆盖/数量缺席证明 + 开放路径滚动调度**。 -### 12.2 全部12个场景结果 +## 11. 复现与参考 -| 版本 | 全清且正常完成 | 平均虚拟时间 / 秒 | 平均检测数 | 平均移动距离 / 米 | 平均逐局时间降幅 | -| --- | ---: | ---: | ---: | ---: | ---: | -| A0 | 12/12 | 13788.04 | 418.00 | 56183.95 | — | -| A1 | 12/12 | 7476.09 | 100.42 | 34146.70 | 43.25% | -| A2 | 12/12 | 5496.20 | 100.92 | 24177.66 | **58.55%** | -| A3 | 12/12 | 6094.73 | 101.67 | 27164.04 | 53.90% | +```bash +# 冻结基准集：bayes vs seven_batch（runner 内置策略名） +uv run python scripts/run_suite.py --problem 3 --suite core --jobs 6 \ + --variants bayes,seven_batch --label q3-core -每个版本共清除156个源，48局合计624个源全部清除，无失败清除、无错误缺席判定、无定位异常。A2相对A0的逐局降幅范围为37.50%至70.84%。降幅按 $\frac1n\sum_i(1-T_{A2,i}/T_{A0,i})$ 计算，不能与两个均值之比混称同一个指标。 +# 项目对比框架（含 grid 基线，生成 project/results/q3.{md,json}） +uv run python -m project.run --question 3 -另以保存的自建场景真值独立复核48局动作日志，逐步检查正观测后的候选保留真源、清除证书及实际清除距离、全清集合、计时、检测数、退出响应与日志哈希，全部通过，见[证据复核结果](C:/Users/HelixCopex/code/MATH/output/q3/search-ablation/evidence-audit.json)。这些真值只用于离线验收，不进入策略。 - -| 源数（每组6场景） | A0平均时间 / 秒 | A1平均时间 / 秒 | A2平均时间 / 秒 | A3平均时间 / 秒 | -| --- | ---: | ---: | ---: | ---: | -| 10 | 15516.58 | 6313.47 | **4836.11** | 5467.38 | -| 16 | 12059.50 | 8638.71 | **6156.28** | 6722.08 | - -10源时，未知频道的无源证明开销较大，所以七点替换收益明显；16源时，可通过清除数量上界提前停止，原网格未必扫完，但批量调度仍显著减少移动。A1与A2检测数相近，A2平均移动距离进一步下降，表明减少往返比单纯少测一次更有价值。这组自建证据支持优先实现A2，而非先投入复杂概率图、强化学习或全局非线性优化。 - -四版本平均本地计算墙钟约1.7至1.8秒/局，记录在数据中；它不含官方网络、真实程序截止或官方UI开销，不能当作新方案官方现实耗时。原官方动作日志只能验证旧路径上的响应，不能用它推断新测点的读数。只有在官方新演练执行新策略后，才能报告官方提速。 - -### 12.3 复现 - -```powershell -.venv/Scripts/python.exe scripts/benchmark_q3_search.py --seeds 0,1,2,3,4 +# 回归测试 +uv run pytest tests/test_particles.py tests/test_bayes_loop_q3.py \ + tests/test_planner.py tests/test_coverage_plan.py tests/test_safety_valves.py -q ``` -也可使用 `uv run python` 调用同一脚本。输出包含原始逐局数据，失败不会从汇总中剔除。当前实验是小规模方案筛选，不能证明所有合法场景、所有固定误差场或不同环半径下的性能。 +主要文件： -## 13. 分阶段改进与验收计划 +- `src/mcm26b/planner/bayes_loop.py`：`bayes` 决策层； +- `src/mcm26b/planner/__init__.py`：`OmnidirectionalLoop`、问题二选点器、清除证书； +- `src/mcm26b/planner/seven.py`：七点策略族与覆盖校验； +- `src/mcm26b/particles.py`：网格幸存集与信息增益； +- `src/mcm26b/belief/__init__.py`：凸候选集； +- `src/mcm26b/coverage/__init__.py`、`src/mcm26b/guarantee/cover.py`：覆盖账本与校核； +- `src/mcm26b/routing.py`：开放路径 DP/2-opt/MST 下界； +- `src/mcm26b/planner/bayes_loop.py::_exact_order`：可选 Held-Karp。 -| 阶段 | 交付 | 必须通过的验收 | -| --- | --- | --- | -| I：七点＋本站批量 | 将A2接入第三问CLI选项，保留A0可切换；原始日志格式兼容 | 七点连续覆盖证明与坐标误差预算；10/16源及边界例全清；第四问覆盖策略未被改动 | -| II：存在计数与扫频次序 | $|\mathcal P|=16$ 提前停止未知频道扫描；同站频道优先独立消融 | 已发现数、已清除数不混淆；发现10个不能停；clear不改变测向频道 | -| III：残余覆盖与顺路观测 | 保守单元账本、负观测复用、有证书跳站 | 边界及单元接缝不漏；未知频道各有独立覆盖责任；启发式停滞自动恢复七点 | -| IV：时间型局部动作 | 最近安全清除点、短程任务预演、插入排序 | 实际提交点满足全部外包顶点约束；near和跨0°读数正确；采样评分不冒充严格上界 | -| V：定位停滞与预算后备 | 有界区域有限清除覆盖，明确上限与退出 | 清除格子覆盖全部候选；格数/两种时间预算不足均记未完成；故障与失败可回放 | - -每个阶段都在相同场景和同一固定误差场上进行配对比较。建议将种子0—4保留为本次开发集；最终验证另生成独立种子100—199，并分别覆盖10至16个源、接收半径1000和1500米、边界源、团簇源、同站near、接近共线读数、误差恰在±1°及空间相关误差场。所有超时、矛盾和未全清样例进入结果表，不允许只报成功局。 - -主要指标为全清且正常退出比例、平均/中位/P95/最大虚拟总时间、移动长度、检测数、换频数、清除失败数、证明无源的成本；同时报告最后清除时刻、现实计算时间与总请求数。小样本不确定性用逐局配对差分及其范围表达，增加样本后可在声明抽样假设的前提下给出区间估计；不从当前12例推出总体成功率为100%。 - -官方验收继续使用 socket 接口，逐条保存请求响应，最后保存官方结果及加密日志；优先使用官方支持的同案例复现功能进行配对，如果无法复现同一案例，就分别报告案例差异，不将两次随机局的时间差当成纯算法收益。本轮只完成方案研究与自建实验，没有再次启动官方演练。 - -## 14. 结论与实现边界 - -第三问首先应把覆盖责任从49点缩减为已证明足够的七点，并将覆盖站内扫频和离站定位拆开。该改进直接针对官方日志中的主要开销，在保持几何定位方法不变时，已经得到12个自建场景的配对支持。推荐先落地A2，再独立验证存在数量上界、残余覆盖和时间型测点/清除点优化。 - -初稿形成时，七点与三种调度变体存在于离线消融脚本，官方仅验证A0；随后A2已通过下述一次官方演练。生产 `src/mcm26b` 默认策略仍未修改，A2由独立运行脚本调用。3B给出了可实现的进一步规划，但没有将其预期收益写成已测结果。固定覆盖、观测集合和清除证书负责可靠性；候选评分、任务排序和半径调参负责效率，二者的保证范围分别陈述。 - -## 15. 补充：A2官方演练已通过 - -在方案研究之后，用户启动新的问题3演练，于2026-09-11通过2026端口实际执行A2。案例 `J4ER-Y9EP-Q379-UPFA` 的官方源数为12，12个全部清除，109次检测、12次清除无失败，现实耗时3.406秒，总虚拟耗时5232.132109秒（1小时27分12秒），移动22645.661米。其余8个频道各有完整七点无信号记录，最终正常退出。 - -123个请求响应已严格离线回放，清除证书和官方微秒计时全部通过复核。原始数据、官方产物、运行代码哈希及复现入口见[本局验证报告](C:/Users/HelixCopex/code/MATH/output/q3/practice/20260911-seven-batch-01/report.md)。因此A2现已有自建配对实验及一个官方案例的有效性证据。 - -相较之前A0官方10源案例的16217.454814秒，本局虚拟时间少67.74%，但案例和源数不同，不能把该差异称为严格算法提速率；第12节的58.55%仍仅指12个自建场景的配对平均结果。3B的进一步改进尚未启用。 +延伸阅读：`docs/algo-3-bayes.md`（早期贝叶斯闭环设计与消融细节）、 +`docs/第三四问架构优化与运行提速方案.md`（运行历史与瓶颈）、 +`docs/第三四问图论与动态规划优化方案.md`（图论/DP 轮实测与负面结论）。 ## 参考文献 -1. 全国大学生数学建模竞赛组委会. 2026年高教社杯全国大学生数学建模竞赛B题：无线电干扰源的快速自动定位与清除[Z]. 2026. [题面](C:/Users/HelixCopex/code/MATH/issue/B题.pdf)；模拟器配套材料《附件2：模拟器通信接口说明及编程指南》[本地文件](C:/Users/HelixCopex/code/MATH/issue/附件/附件2.docx)。 -2. 本项目. 问题一：有界测向误差下定位区域的直径与圆覆盖；问题二：具有接收保证的第二检测点选择[Z]. 2026. [algo-1](C:/Users/HelixCopex/code/MATH/docs/algo-1.md)、[algo-2](C:/Users/HelixCopex/code/MATH/docs/algo-2.md)。 -3. BAYRAM H, VANDER HOOK J, ISLER V. Gathering Bearing Data for Target Localization[J]. IEEE Robotics and Automation Letters, 2016, 1(1): 369–374. DOI: 10.1109/LRA.2016.2521387. [作者所在大学的研究记录](https://experts.umn.edu/en/publications/gathering-bearing-data-for-target-localization-6/)。核对题目、作者、页码及移动与测量联合代价；未套用其GDOP近似比。 -4. OLOFSSON J, HENDEBY G, LAUKNES T R, JOHANSEN T A. Multi-Agent Informed Path Planning Using the Probability Hypothesis Density[J]. Autonomous Robots, 2020. [作者公开稿](https://torarnj.folk.ntnu.no/AURO2020a.pdf)，第1—2节。用于未观测区域与预算路径的建模对照。 -5. PERERA A, TRAN V P, ANAVATTI S, KASMARIK K, GARRATT M. Radio Source Localization using Sparse Signal Measurements from Uncrewed Ground Vehicles[EB/OL]. arXiv:2312.03493, 2023. [原文入口](https://arxiv.org/abs/2312.03493)。用于稀疏覆盖与实验设计的对照，RSSI估计器未迁移。 -6. ATANASOV N A, LE NY J, PAPPAS G J. Distributed Algorithms for Stochastic Source Seeking with Mobile Robot Networks: Technical Report[EB/OL]. arXiv:1402.0051v2, 2014. [原文入口](https://arxiv.org/abs/1402.0051)。用于说明场强梯度和测量模型假设的适用边界，不引用其结论证明本算法收敛。 - -在线来源核对日期：2026-09-11。本文的七点构造、覆盖界和具体调度规划来自对题设参数的推导与本项目实验，不宣称一般几何和路径规划思想在文献中首次出现。 +1. 全国大学生数学建模竞赛组委会. 2026 年高教社杯全国大学生数学建模竞赛 B 题[Z]. 2026. 本地题面 `issue/B题.pdf`。 +2. 模拟器配套材料. 附件 1《模拟器使用说明》、附件 2《通信接口说明及编程指南》[Z]. 2026. 本地 `issue/附件/`。 +3. 本项目. `docs/algo-1.md`、`docs/algo-2.md`、`docs/algo-4.md`。 +4. 本项目. `docs/算法审核与优化-20260911.md`、`docs/第三四问架构优化与运行提速方案.md`、`docs/第三四问图论与动态规划优化方案.md`。 diff --git a/docs/algo-4.md b/docs/algo-4.md new file mode 100644 index 0000000..5c8da37 --- /dev/null +++ b/docs/algo-4.md @@ -0,0 +1,485 @@ +# 问题四算法说明：混合全向/定向源搜索与清除闭环（当前实现） + +> 本文汇总当前仓库中问题四的全部实际算法。生产策略是 +> `src/mcm26b/planner/__init__.py::DirectionalExtension`，CLI 策略名 +> `directional`，入口为 `mcm26b run --problem 4`。 +> 与问题三不同，问题四的源可能包含定向源；`no_signal` 不再等价于 +> “距离超过接收半径”，因此问题三的七站全向缺席证明**不能直接搬用**。 +> 当前仓库尚无官方问题四演练记录；本文性能数字均来自自建 +> `LocalSimulator` 冻结基准集，明确标注为测试假设。 + +## 1. 问题定位与核心难点 + +问题四每局有 $10\sim16$ 个源，其中一部分为定向源，其余为全向源； +位置在 $\Omega=B(0,1800)$ 内，频道互异。定向源有一个朝向 +$h$，只有当接收点位于其**前向半平面**内、且距离不超过 $r_{\rm eff}$ +时才可能返回有效读数： + +$$ +r=\lVert q-P\rVert\le r_{\rm eff}\quad\text{且}\quad +\operatorname{angular\_gap}\big(\operatorname{bearing}(P\to q),h\big) +\le90^\circ . +$$ + +其中 $r_{\rm eff}$ 对每个源未知，满足 $1000\le r_{\rm eff}\le1500$ m。 +若 $r>r_{\rm eff}$ 或朝向不覆盖接收点，结果为 `no_signal`；若可见且 +$r\le5$ m，结果为 `near`；若可见且 $r>5$ m，结果为 `direction`， +数值为接收点到源的方位加有界误差。 + +因此问题四的困难不是一般的“覆盖搜索”，而是两个额外的证明/决策问题： + +1. **发现**：如何用有限的访问点证明“某个频道没有源”，即使该源是 + 朝外发射的定向源； +2. **定位与清除**：定向源的 `no_signal` 可能只是朝向遮挡，不能像全向源 + 那样删除位置候选；需要推理“源可能朝哪边”，并在定位停滞时仍能可靠 + 清除。 + +当前实现的办法是： + +- 发现层用 `forward_half_disk` 语义的 31 点最小命中集覆盖计划 + （`q4_mincover_v2`）； +- 感知层在原有位置网格上叠加 24 档朝向掩码，得到（位置 × 朝向）联合 + 可行集 `JointBelief`； +- 决策层用联合信念的可见性概率筛选探点，失败时按预算进入 ≤20 m 的 + 有限格后备清除； +- 调度层用“本站批量扫频 + 按候选 MEC 圆心最近服务队列”的滚动顺序， + 避免发现后立即远行、再返回同站扫描的往返。 + +| 层次 | 文件/类 | 作用 | +| --- | --- | --- | +| 覆盖计划 | `src/mcm26b/guarantee/cover.py`、`guarantee/plans/q4_mincover_v2.json` | `forward_half_disk` 最小命中集与运行期校核 | +| 位置网格 | `src/mcm26b/particles.py::ChannelParticles` | 位置幸存集与 Voronoi 保守裕量 | +| 朝向联合信念 | `src/mcm26b/belief/joint.py::JointBelief` | 每格 24 档朝向掩码、可见性/期望后验半径 | +| 决策闭环 | `src/mcm26b/planner/__init__.py::DirectionalExtension` | 批量扫频、可见性探点、后备清除、收尾 | +| 对比框架 | `project/q4/{optimal,baseline,compare}.py` | `directional` vs `bayes`/`grid` 直接搬用 | +| 回归 | `tests/test_loop_q4.py`、`tests/test_joint_belief.py`、`tests/test_coverage_plan.py` | 真源保留、朝外边界源、后备纪律、证书校核 | +| 结果 | `project/results/q4.{md,json}` | 11 个 core+boundary 场景配对结果 | + +## 2. 发现覆盖：`forward_half_disk` 与最小命中集 + +### 2.1 覆盖语义 + +设访问点集合为 $S$。对全向源，`no_signal` 在 $s$ 处排除 $B(s,1000)$， +因此完整覆盖只需 + +$$ +\Omega\subseteq\bigcup_{s\in S}B(s,1000). +$$ + +对定向源，还需要朝向条件。严格的缺席证明要求： + +$$ +\forall P\in\Omega,\ \forall h\in\mathbb S^1,\quad +\exists s\in S:\quad +\lVert s-P\rVert\le1000 +\ \wedge\ (s-P)\cdot u(h)\ge0 . +$$ + +也就是说，对任意位置和任意朝向，至少有一个被访问点同时落在该源的 +1000 m 距离内和 90° 前向半平面内。这个语义称为 `forward_half_disk`。 +问题三的七站环不满足该语义：`tests/test_coverage_plan.py` 中有专门用例 +验证“七站点集在 `forward_half_disk` 下有缺口”。 + +### 2.2 49 点 legacy 充分条件 + +旧实现使用 $S=\{(700i,700j):i,j\in[-3,3]\}$ 的 49 个点： + +- 任意 $P\in[-2100,2100]^2$ 到某个网格点的欧氏距离不超过 + $700\sqrt2\approx989.95\le1000$； +- 因此对该网格点，无论源朝向如何（前向半平面总是 $180^\circ$）， + 距离条件都成立。 + +这是 49 点充分条件，保证发现但代价很高：最近邻开放巡回约 36.4 km； +对无源频道还要在 49 个点各测一次。它现在保留为 legacy 对照与回退计划 +`square49`。 + +### 2.3 31 点最小命中集 `q4_mincover_v2` + +当前默认计划由离线贪心集合覆盖搜索得到，并带有可运行期校验证书。 +生成参数：搜索分辨率 50 m / 48 档，设计距离裕量 +$r_{\rm design}=965$ m，前向角裕量 $5^\circ$；随后用证书分辨率 +15 m / 72 档（$3\,256\,200$ 个 $(P,h)$ 采样）做精细剪枝，每删一个点都 +用同一证书口径复核。 + +`q4_mincover_v2` 的证书摘要： + +| 字段 | 值 | +| --- | ---: | +| 点数 | 31 | +| 覆盖语义 | `forward_half_disk` | +| 校核分辨率 | 15 m / 72 朝向档 | +| 采样目标数 | 3 256 200 | +| 缺口数 | 0 | +| 最坏命中距离 | 990.051 m | +| 最近邻开放巡回 | 23 883.8 m（约 4776.8 s 虚拟移动） | +| 生成器 SHA-256 | `adf4fe6e…434c0e1d`（完整值见 JSON） | + +旧 `q4_mincover_v1`（31 点）在粗采样 75 m/24 档下看似零缺口，但在 +25 m/48 档有 125 个缺口、15 m/72 档有 566 个缺口，属于采样过拟合； +`q4_mincover_v2` 修正了这一点，点数不变而细网格缺口归零。 + +**严格连续空间证书的已知残余风险.** 15 m/72 采样格上的零缺口加上设计 +裕量可以覆盖相邻未采样 $(P,h)$，但默认 v2 的最坏命中距离 990.051 m +略高于严格连续界要求的约 989.4 m（需满足 +$r_{\rm verify}+15\sqrt2/2\le1000$，即 $r_{\rm verify}\le989.4$ m）。 +仓库另有 36 点备选计划（最坏命中 973.7 m、巡回约 25.4 km，见 +`output/q4/cover-certificate-v3.json`）可关闭该连续性缺口；当前默认的 +v2 在 15 m/72 离散证书上零缺口，连续空间残余约 0.65 m，已作为已知 +风险登记。若论文要求严格连续证明，应改用 36 点方案并重新验收性能。 + +### 2.4 运行期校验 + +`load_plan(name, verify=True)` 会读取计划 JSON 中的证书字段，并用 +**证书自带的采样分辨率**重新校核；缺口不为 0 直接拒绝启动。默认 +`search.verify_coverage_plan=True`。校核成本在进程内按计划名缓存一次， +15 m/72 档约数秒（见 `scripts/run_suite.py` 的 `plan_verify_wall_s`）。 +运行期不会用比证书更粗的采样去“自证”。 + +## 3. （位置 × 朝向）联合信念 `JointBelief` + +### 3.1 网格与朝向档 + +位置维度复用 `ChannelParticles`：间距 $h=14.65$ m、Voronoi 裕量 +$m=\sqrt2/2\,h\approx10.359$ m；朝向维度把整周分为 +24 档（每档 $15^\circ$，半档 $7.5^\circ$）。每个位置格附带一个 24 位 +掩码，表示该格仍可能具有的朝向；缺省掩码为全档。 + +联合集合的保守性目标是：**真源 (位置, 朝向) 永不被删除**（在该离散 +表示意义下，即真源所在格的朝向至少保留一档，且位置格不被删）。因此 +所有判定都向“保留更多”的一侧取整，并计入格内角裕量 + +$$ +\alpha_{\rm cell}=\arcsin\!\left(\frac{m}{\max(d-m,1)}\right), +$$ + +即格内任一点相对格心到测点方位的最大偏差上界。 + +### 3.2 观测更新规则 + +| 观测 | 位置维度 | 朝向维度 | +| --- | --- | --- | +| `direction` @ $s$，读数 $\beta$ | 删除距 $s>1500+m$ 的格；删除角向偏离读数超过 $\delta_{\rm eff}+\alpha_{\rm cell}$ 的格 | 朝向必须**可能覆盖** $s$：保留与方位 $\operatorname{bearing}(g\to s)$ 夹角不超过 $90^\circ+7.5^\circ+1.5^\circ+\alpha_{\rm cell}$ 的档 | +| `near` @ $s$ | 删除距 $s>5+m$ 的格 | 同 `direction`：朝向必须可能覆盖 $s$ | +| `no_signal` @ $s$ | **位置不删**（可能是距离不足或朝向遮挡） | 仅对“整格必在 $r_{\rm eff,min}=1000$ 圆盘内”的格，保留与方位夹角大于 $90^\circ-7.5^\circ-1.5^\circ-\alpha_{\rm cell}$ 的档，即只保留**确定看不见** $s$ 的朝向档；掩码变空则删除该位置格 | + +其中 `1.5°` 是朝向查表分辨率带来的保守余量，$7.5^\circ$ 是半档。 +`no_signal` 的距离阈值使用 $1000-m$（保守侧）：只有在确信整格都在 +1000 m 内时，才把静默解释为朝向遮挡；否则保留全部朝向档，避免删掉 +真源。 + +> 2026-09-12 修复记录：早期实现把 `no_signal` 分支的朝向掩码方向写反， +> 会删除真源朝向档；现已改为“只保留确定看不见的档”，并由 +> `tests/test_joint_belief.py::test_true_source_survives_consistent_observations` +> 以随机观测序列固化。 + +联合集一旦为空（所有位置格都被删除），就给出该频道不存在的严格缺席 +证明，`absence_evidence` 记为 `joint_empty`。这与覆盖证明是两类独立证据。 + +### 3.3 可见性与值函数 + +`JointBelief` 提供两个只用于选点排序的启发式量： + +1. `prob_detect(q)`：在测点 $q$ 能看见该频道任一可行 + (位置, 朝向) 的概率质量。对抽样格心 $P$： + - 距离概率：$d\le1000$ 时为 1；$d\ge1500$ 时为 0；中间线性插值 + （对应 $r_{\rm eff}\sim U[1000,1500]$）； + - 朝向覆盖比例：该格掩码中与“$P\to q$ 方位”夹角不超过 + $90^\circ$ 的档位占比； + - 两者相乘后对格取平均。 +2. `expected_posterior_radius(q)`：一次测量后的期望后验 MEC 半径（D1 + 值函数）。对每个代表格： + +$$ +p_{\rm vis}=(1-p_{\rm silent\_range})\cdot\text{cover\_ratio}, +\quad +p_{\rm ns}=1-p_{\rm vis}, +$$ + + - `direction` 分支：后验半径用“以 $P\to q$ 方位为读数、$\pm\delta_{\rm eff}$ + 角域裁剪后的子采样质心最大距离”近似； + - `no_signal` 分支：位置候选**不缩小**，半径保持当前 PF MEC 半径； + - `near` 分支：后验为 5 m 盘。 + 最后对代表假设取平均。 + +这两个量的先验（均匀格点、$r_{\rm eff}$ 均匀分布）是工程假设，不承担 +任何安全责任；安全责任由覆盖证书、清除证书和联合真源保留定理承担。 + +## 4. 决策闭环 `DirectionalExtension` + +### 4.1 覆盖阶段：本站批量扫频 + 最近队列服务 + +```text +while 未达到数量上界且未停止: + index = 尚有 UNKNOWN 待测的最近覆盖点（31 点计划） + if 不存在: break + pending = 该点尚未测过的 UNKNOWN 频道 + pending 按“当前测向频道优先”排序 + for channel in pending: + outcome = /measure(point, channel) + observe(channel, point, outcome) # 更新凸信念 + JointBelief + if outcome == near 且候选集有清除证书: + 就地 /clear + queue = 所有 DETECTED 且无证书的频道 + while queue: + channel = 距其候选集 MEC 圆心最近的频道 + 从 queue 移除并调用 _localize_and_clear(channel) +``` + +与问题三 `SevenBatch` 的结构类似，但覆盖点集、缺席语义和定位逻辑完全不同。 +`observe` 对 `direction`/`near` 更新凸候选集和联合信念；对 `no_signal` +只更新联合信念的朝向维度，绝不用于删除凸候选集的位置。 + +覆盖阶段结束后（所有计划点上的 `UNKNOWN` 频道都测过）： +`_mark_absent_after_coverage` 把覆盖账本完整的 `UNKNOWN` 频道标为 +`ABSENT`；然后 `_resolve_detected` 再对剩余 `DETECTED` 频道做第二轮 +定位。 + +### 4.2 定位子循环 + +`DirectionalExtension._localize_and_clear(channel)`： + +```text +步数 steps = 0，被遮挡探点表 blocked = [] +while steps < max_localize_steps=8: + if 凸/PF 清除证书已成立: 尝试 /clear；成功则返回 + if len(blocked) >= max_blocked_probes=4: break # 定位停滞 -> 后备 + q = _next_probe(region, channel) + within_range = (|q - region.MEC 圆心| <= safe_disk_radius(region)) + outcome = /measure(q, channel) + observe(...) + if outcome == no_signal: + if within_range: 标记 directional_observed + blocked.append(q); continue # 不消耗 steps + else: # direction 或 near + blocked.clear(); steps += 1 +循环结束后再检查证书；若仍不可清除 -> _fallback_clear(channel) +``` + +要点： + +- 只有**在保证接收盘内**出现 `no_signal` 才标记“观测到定向”，避免把 + 距离不足误诊为朝向遮挡； +- 被遮挡探点不重复选取；最多 4 次，随后转入后备清除； +- `direction`/`near` 读数才计入定位步数； +- 每一步后都重新评估证书，滚动推进。 + +### 4.3 `_next_probe`：联合信念选点与后备梯子 + +候选点族 `_probe_candidates` 包含： + +1. `select_measurement_point` 给出的保证接收盘内最坏后验最优候选 + （复用问题二/问题三的选点器）； +2. 若已有被遮挡点，取关于 MEC 圆心的对跖方向点，并裁进保证接收盘； +3. 抵近点序列 `center + α*(current-center)`（代码保留，默认 + `approach_fractions=()` 未启用）； +4. 以 MEC 圆心为心、半径 $0.5\rho,\ 1.0\rho,\ 1.5\rho$ 的 8 方向环探。 + +候选去重到 0.5 m，并可排除已测过且被遮挡的点。 + +若联合信念可用，优先调用 `_visibility_probe`（默认口径 +`probe_score_mode=nearest_visible`）： + +1. 对每个候选计算 `prob_detect(q)`； +2. 只保留概率不低于最大值 0.7 倍的候选 + （`probe_visibility_ratio=0.7`），避免选中“期望上有点信息、实际必被 + 朝向遮挡”的点； +3. 在保留候选中选择 + $\text{move}(p\to q)+\text{move}(q\to m_{\rm MEC})+6$ 秒最小者 + （6 = 检测 5 s + 换频 1 s）。 + +若联合信念不可用/为空或候选为空，则退回确定性梯子： +未被遮挡过用最坏后验选点器；第一次被遮挡用对跖点；之后绕 MEC 圆心 +半径 $\rho$ 的环上每次 +45° 试探。几何边界必须如实说明：环上可见弧 +没有正的下界（源可在环上且朝向沿外切线，可见集退化为切点），因此 +梯子只是提高命中概率，**完备性由后备清除兜底**。 + +### 4.4 后备清除：有限格盲清 + +定位停滞时，`_fallback_clear` 对凸候选多边形生成 20 m 格心的有限覆盖： + +1. 用全局 20 m 半开网格与候选多边形求交（带正容差，宁可多保留格）； +2. 格内任一点到格心不超过 $10\sqrt2\approx14.14\le20$ m，因此若真源 + 仍在候选多边形内，逐格 `/clear` 至少命中一次； +3. 格数上界 + $(\lceil w/20\rceil+2)(\lceil h/20\rceil+2)$； +4. 若格数超过 `fallback_max_cells=400`，记录 + `fallback_cell_cap_exceeded` 并跳过，不硬扫； +5. 若现实剩余预算不足以覆盖 + `格数 × real_rtt_p95_s + reserve_exit_s`，记录 `skipped_budget` + 并跳过； +6. 通过检查时，按最近邻开放路径逐格 `/clear`；命中成功即标记清除； + 全部失败则记录 `fallback_exhausted_without_success`，**不删除候选集**。 + +后备清除只依赖“清除与朝向无关”这一题面事实，是定向源定位失败时唯一 +不依赖可见性的完备性底线；但它受格数上限和现实预算约束，预算不足时 +会诚实地保持未完成。`fallback_optimize_tour=False` 为默认：实测把 +格间路径排得更短反而会让“命中即停”的真源顺序变差，故不启用。 + +## 5. 调度、完成条件与预算 + +### 5.1 与问题三调度层的关系 + +`DirectionalExtension` 继承 `OmnidirectionalLoop` 的会话、预算、异常 +和 summary 骨架，但不使用 `BayesLoop` 的统一任务表/信息增益调度； +它采用“每站批量扫频 + 本站服务队列”的滚动策略（类似 A2 的推广）。 +因此： + +- 检测/换频/清除的计时口径与问题三完全一致； +- 覆盖点、缺席证据、定位方法、清除证书不同； +- 默认不启用 Held-Karp/2-opt 任务表；后备清除走最近邻。 + +### 5.2 完成条件 + +- 清除数达到 16：直接停止并按 `count_cap` 完成； +- 否则要求每个频道 `CLEARED` 或 `ABSENT`，按 `per_channel_proof` 完成； +- 缺席证据来源： + - `coverage`：该频道在 31 点计划上全部测得 `no_signal`； + - `joint_empty`：联合可行集为空； +- 定位/清除尚未完成或预算耗尽：`completed=False`， + `stop_reason=unresolved_channels` / `real_budget_exhausted` / + `virtual_time_cap` 等，不把未完成说成成功。 + +## 6. 关键参数 + +| 参数 | 默认值 | 说明 | +| --- | ---: | --- | +| `q4_coverage_plan` | `q4_mincover_v2` | 31 点定向最小命中集 | +| `verify_coverage_plan` | `True` | 运行期按证书分辨率校核 | +| `directional_joint_belief` | `True` | 启用（位置 × 朝向）联合信念 | +| `probe_score_mode` | `nearest_visible` | 可见性过滤 + 移动最省 | +| `probe_visibility_ratio` | 0.7 | 可见概率保留比例 | +| `probe_exclude_blocked` | `True` | 不重复测已遮挡点 | +| `approach_fractions` | `()` | 抵近点默认关闭 | +| `max_localize_steps` | 8 | 每轮定位成功读数上限 | +| `max_blocked_probes` | 4 | 被遮挡探点上限，超限转后备 | +| `fallback_clear_grid_m` | 20 m | 后备清除格边长；$10\sqrt2<20$ | +| `fallback_max_cells` | 400 | 后备格数硬顶 | +| `fallback_optimize_tour` | `False` | 后备走格 NP+2-opt 默认关闭 | +| `persistent_search` | `False` | 找全证明后立即退出 | + +## 7. 已验证结果与失效机理 + +### 7.1 冻结基准集 + +`project/q4` 在 11 个冻结 core+boundary 场景上的配对结果： + +| 方法 | 有效局 | 全清率 | 缺席误判 | 虚拟时间均值 | 中位 | 最大 | 后备清除均值 | +| --- | ---: | ---: | ---: | ---: | ---: | ---: | ---: | +| `directional` | 11/11 | **100%** | **0** | 9069.9 s | 8803.3 s | 12081.1 s | 24.8 | +| `bayes` 直接搬用 | 0/11 | 77% | **24** | 3890.8 s | 3646.2 s | 5054.4 s | 0 | +| `grid` 直接搬用 | 0/11 | 67% | 0 | 18179.2 s | 18345.2 s | 20457.9 s | 0 | + +两条搬用基线的失效机理不同： + +- `bayes` 使用七站**全向**缺席证明；在朝外发射、$r_{\rm eff}=1000$ m 的 + 边界定向源场景中会把真源直接判成 `ABSENT`。它“快而错”，累计 24 次 + 缺席误判，虚拟时间不能作为有效优势； +- `grid` 的 49 点覆盖对定向源发现是完备的，但没有朝向推理与后备清除， + 定位停滞时只能停在 `DETECTED`，属于“稳而不全”。 + +`directional` 在 11/11 场景全清、零缺席误判；代价是时间更长，平均 +24.8 次后备清除。纯 core 9 场景（不含朝外边界压力）的均值约 8602.0 s、 +最大 9630.4 s、后备 4.3 次/局；boundary 场景（朝外发射源）会把后备 +推到数十至上百次，单局最大 12081.1 s（16 源 8 定向边界场景）。 +因此第四问的时间对“源是否背离圆心发射”非常敏感，不存在对所有合法 +场景的均匀性能保证。 + +### 7.2 覆盖计划消融 + +| 计划 | 点数 | 15 m/72 缺口 | 最坏命中 | NN 巡回 | +| --- | ---: | ---: | ---: | ---: | +| `square49` | 49 | 0 | ≤1000 m（构造） | 36.4 km | +| `q4_mincover_v1` | 31 | 566 | 1000.0 m | 22.9 km | +| `q4_mincover_v2`（默认） | 31 | **0** | 990.051 m | 23.9 km | +| 36 点严格连续备选 | 36 | 0 | 973.7 m | 25.4 km | + +改用最小命中集后，覆盖点减少了 37%，巡回长度从 36.4 km 降到 23.9 km； +但第四问时间的大头仍是 23.9 km 巡回（约 4776.8 s）和逐频道缺席检测。 +把七站式的“有源提前停止”规则用于 31 点计划可在一定程度上减少检测， +但最小命中集的巡回与检测开销仍是该问题的成本地板（详见 +`docs/第三四问图论与动态规划优化方案.md` §2、§8）。 + +### 7.3 回归测试 + +问题四相关测试覆盖： + +- 联合信念真源 (位置, 朝向) 在随机观测序列后不丢； +- 朝外发射的边界定向源不漏、不被误判缺席； +- 保证接收盘内 `no_signal` 不删空间候选、判向不越界； +- 后备清除格覆盖候选多边形每个见证点；格数上限/预算检查跳过而非硬扫； +- 全清、零异常、逐微秒计时复算一致。 + +当前全仓 `pytest` **268 项通过**（2026-09-12 本轮复核）。 + +## 8. 与问题一、二、三当前实现的复用关系和差异（核对结果） + +| 对比 | 复用 | 差异 | +| --- | --- | --- | +| 问题一 | 角域半平面语义、Welzl 最小包围圆 `smallest_enclosing_circle`、MEC 清除证书思路 | 问题四不计算直径/直径圆；`no_signal` 不能解释为距离排除；覆盖语义是 `forward_half_disk` 而非纯角度交集 | +| 问题二 | `select_measurement_point`、`safe_disk_radius` 作为基础候选与“是否在保证接收盘内”的判据；清除点可向当前位置平移 | 定向源在安全盘内仍可能因朝向 `no_signal`，因此问题二的“盘内必收到信号”保证不成立；新增联合信念和可见性筛选 | +| 问题三 | PF 网格位置维度、覆盖账本、会话/预算/异常 summary、清除动作 | 问题三的负观测可删位置候选（`full_disk` 七点缺席证明）；问题四的负观测只约束朝向，覆盖计划改用 31 点 `forward_half_disk`，不使用七站全向缺席证明；问题四也不使用 `BayesLoop` 的信息增益任务表 | +| 问题四独有 | — | 24 档朝向掩码、可见性概率选点、对跖/环探后备梯子、20 m 格后备清除及其格数/预算纪律 | + +关键结论：**问题四不是“问题三加一点朝向”，而是发现证明、负观测语义和 +清除兜底三处都必须重做**。直接套用问题三的 `bayes`（七站缺席证明）会在 +朝外定向源上产生缺席误判；直接套用 `grid`（49 点全向）可以发现但无法 +可靠定位清除。当前 `directional` 的路线是：最小命中集保证发现、 +联合信念推理可见性、预算化后备清除保证清除底线。 + +## 9. 适用边界与风险声明 + +- **覆盖证书**：默认 v2 是 15 m/72 档离散证书加设计裕量，零缺口； + 严格连续空间证书还存在约 0.65 m 的已知残余，正式论文若要求严格连续 + 证明应改用 36 点备选计划； +- **联合信念**：真源保留证明依赖评测机观测语义与题面误差界一致； + 先验分布、朝向均匀抽样和可见性概率都只是决策启发式； +- **后备清除**：完备性只在格数不超过 400、现实预算足够且候选集确实 + 包含真源时成立；任一条件不满足就诚实地保持未完成； +- **性能**：11/11 全清来自自建模拟器与冻结场景，不是官方成绩； + 当前仓库尚无官方问题四演练。第四问对定向源朝向分布敏感，不能把 + core 均值 8602 s 或边界均值 9069.9 s 外推为官方期望； +- **调度**：站点批量与最近队列是贪心，不是全局最优；图论/DP 方案 + 中的插入、Or-opt、精确任务序等在本场景实测为负收益或未推广到问题四， + 默认关闭。 + +## 10. 复现 + +```bash +# 冻结基准集：directional vs 问题三策略直接搬用 +uv run python scripts/run_suite.py --problem 4 --suite all --jobs 6 \ + --variants directional --label q4-current + +# 项目对比框架（生成 project/results/q4.{md,json}） +uv run python -m project.run --question 4 + +# 覆盖证书生成与细网格校核（默认计划的来源） +uv run python scripts/verify_q4_cover.py --spacing 100 --radius 2850 \ + --target-step 50 --headings 48 --design-radius 965 --heading-margin 5 \ + --verify-step 15 --verify-headings 72 --fine-prune \ + --write-plan q4_mincover_v2 --out output/q4/cover-certificate-v2.json + +# 回归 +uv run pytest tests/test_joint_belief.py tests/test_loop_q4.py \ + tests/test_coverage_plan.py tests/test_particles.py -q +``` + +主要文件： + +- `src/mcm26b/planner/__init__.py::DirectionalExtension`、`fallback_clear_cells`、 + `fallback_cell_bound`：闭环、探点与后备清除； +- `src/mcm26b/belief/joint.py`：联合信念； +- `src/mcm26b/particles.py`：位置幸存集； +- `src/mcm26b/guarantee/cover.py`、`guarantee/plans/q4_mincover_v2.json`： + 覆盖计划与证书； +- `src/mcm26b/routing.py`：后备/离线巡回度量。 + +延伸阅读：`docs/第三四问架构优化与运行提速方案.md`（第四问失败机理与 +P0-2/P0-3 实测）、`docs/第三四问图论与动态规划优化方案.md`（联合信念 +修复记录与 D1 探点口径消融）、`docs/algo-3.md`（问题三全向闭环）。 + +## 参考文献 + +1. 全国大学生数学建模竞赛组委会. 2026 年高教社杯全国大学生数学建模竞赛 B 题[Z]. 2026. 本地题面 `issue/B题.pdf`。 +2. 模拟器配套材料. 附件 2《通信接口说明及编程指南》[Z]. 2026. 本地附件 `issue/附件/附件2.docx`。 +3. 本项目. `docs/algo-1.md`、`docs/algo-2.md`、`docs/algo-3.md`。 +4. 本项目. `docs/第三四问架构优化与运行提速方案.md`、`docs/第三四问图论与动态规划优化方案.md`、`docs/算法审核与优化-20260911.md`。 diff --git a/docs/paper/main.pdf b/docs/paper/main.pdf index 2515c4e..d2b1d31 100644 Binary files a/docs/paper/main.pdf and b/docs/paper/main.pdf differ diff --git a/docs/第三四问图论与动态规划优化方案.md b/docs/第三四问图论与动态规划优化方案.md index 09adaa4..ad62d79 100644 --- a/docs/第三四问图论与动态规划优化方案.md +++ b/docs/第三四问图论与动态规划优化方案.md @@ -2,7 +2,7 @@ - 版本：2026-09-12 - 前置文档：[第三四问架构优化与运行提速方案](第三四问架构优化与运行提速方案.md)（架构与 P0–P2 项）、 - [algo-3-bayes.md](algo-3-bayes.md)（现役第三问策略）、[algo-3.md](algo-3.md)（七点覆盖证明） + [algo-3-bayes.md](algo-3-bayes.md)（现役第三问策略）、[algo-3.md](algo-3.md)（问题三当前算法说明） - 范围：只讨论"用图论 / 动态规划类算法进一步压低**虚拟耗时**"这一件事； 题面约束、协议语义、保守性不变量全部不动 - 本文所有数字均为**本轮新增实测**（自建内存回环、冻结场景、现役代码）， diff --git a/scripts/build_docs_figures.py b/scripts/build_docs_figures.py new file mode 100644 index 0000000..004d83c --- /dev/null +++ b/scripts/build_docs_figures.py @@ -0,0 +1,1309 @@ +"""从行为日志 JSONL 复现 WebUI"回放"视图，并出论文级 PNG/SVG。 + +本脚本是 docs/figures/ 全部插图的唯一生成入口： + + uv run python scripts/build_docs_figures.py # 全量 + uv run python scripts/build_docs_figures.py --only flow,replay + uv run python scripts/build_docs_figures.py --list + +三类数据来源： + +1. **题目/规则常量**（``mcm26b.config`` 与 ``docs/`` 中的解析结论）——流程图、 + 计时示例、七点/31 点覆盖示意图，全部由公式现算，不抄数字； +2. **实测行为日志 JSONL**——执行时序分解图与回放复现图。优先读取仓库内 + ``referee-data/sessions/*/log.jsonl``（评测机日志，含真值），其次读取 + ``examples/replay/*.jsonl``（去重后的归档，供无本机评测机数据时复现）； +3. **官方演练留档**（``output/q1/practice/*_observations.json``，若存在）—— + 问题一交会定位区域图；缺失时回落到 ``examples/q1/triangle.json`` 反例算例。 + +所有图同时输出 SVG（WebUI/文档内联预览）与 PNG（论文插图，300+ dpi 等效）， +并把数据来源与 sha256 写入 ``docs/figures/manifest.json``。 +""" + +from __future__ import annotations + +import argparse +import json +import math +import shutil +import sys +from dataclasses import dataclass +from pathlib import Path +from typing import Any, Iterable, Sequence + +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +from figkit import (DEFAULT_RASTERIZER, PALETTE, ROOT, Figure, # noqa: E402 + write_manifest) + +OUT_DIR = ROOT / "docs" / "figures" +ARCHIVE_DIR = ROOT / "examples" / "replay" +SESSIONS_DIR = ROOT / "referee-data" / "sessions" + +# 题面常量（与 src/mcm26b/config.py 一致；此处仅用于画图，不参与求解） +OMEGA_R = 1800.0 +R_EFF = 1000.0 +SPEED_MPS = 5.0 +MEASURE_S = 5.0 +SWITCH_S = 1.0 +CLEAR_OK_S = 5.0 +CLEAR_NO_S = 3.0 +RING_RADIUS_M = 1159.7 +SEVEN_R = 1800.0 +GRID_STEP_M = 700.0 +Q4_DESIGN_R = 965.0 +Q4_HEADING_MARGIN_DEG = 5.0 +Q4_WORST_HIT_M = 990.051 + + +# =========================================================================== +# 行为日志解析：评测机日志与参赛队客户端会话日志统一成同一事件流 +# =========================================================================== +@dataclass +class LogRun: + """一局可回放的行为日志（与 referee `replay.js` 的归一化口径一致）。""" + + label: str + case_code: str + module: str + robot_id: str + truth: list[dict[str, Any]] + actions: list[dict[str, Any]] # 保序：{seq, path, request, response} + source: Path + + @property + def clears(self) -> int: + return sum(1 for a in self.actions if a["path"] == "/clear") + + @property + def measures(self) -> int: + return sum(1 for a in self.actions if a["path"] == "/measure") + + @property + def virtual_s(self) -> float: + last = self.actions[-1]["response"] if self.actions else {} + return float(last.get("virtual_time_s") or 0.0) + + def solution_channels(self) -> list[int]: + return sorted({int(a["request"]["channel"]) for a in self.actions + if a["path"] == "/clear" and "channel" in a["request"]}) + + +def parse_jsonl(path: Path) -> LogRun: + """兼容两种日志：评测机 ``event=action/session_start`` 与客户端 + ``event=accepted/attempt/step``（取最终 accepted，重试不重复计）。""" + meta: dict[str, Any] = {} + actions: list[dict[str, Any]] = [] + fallback: dict[str, dict[str, Any]] = {} + seq = 0 + for line in path.read_text(encoding="utf-8").splitlines(): + line = line.strip() + if not line: + continue + try: + ev = json.loads(line) + except json.JSONDecodeError: + continue + kind = ev.get("event") + if kind == "session_start": + meta = ev + elif kind == "action": + seq += 1 + actions.append(dict(seq=ev.get("seq", seq), path=ev.get("path", ""), + request=ev.get("request") or {}, + response=ev.get("response") or {}, + wall_ts_ms=ev.get("wall_ts_ms"))) + elif kind == "accepted": + fallback[ev.get("request_id", "")] = ev + elif kind == "step" and ev.get("request_id") not in fallback: + fallback[ev.get("request_id", "")] = ev + if not actions: # 客户端日志：attempt → accepted → step + for ev in fallback.values(): + if not ev.get("path"): + continue + seq += 1 + actions.append(dict(seq=seq, path=ev["path"], request=ev.get("request") or {}, + response=ev.get("response") or {}, + wall_ts_ms=ev.get("logged_at_unix_s", 0) * 1000 + if ev.get("logged_at_unix_s") else None)) + if not actions: + raise ValueError(f"{path} 中没有可回放的动作事件") + return LogRun(label=path.stem, case_code=meta.get("case_code") or path.parent.name, + module=meta.get("module", "?"), robot_id=str(meta.get("robot_id", "?")), + truth=list(meta.get("truth") or []), actions=actions, source=path) + + +@dataclass +class Timeline: + """逐动作的计时分解：与评测机引擎/附件规则同源。""" + + rows: list[dict[str, Any]] + total_virtual_s: float + total_real_s: float + total_move_m: float + + +def decompose(run: LogRun) -> Timeline: + """按附件计时规则逐动作复算：move/switch/measure/clear 四类耗时。""" + pos = (0.0, 0.0) + channel: int | None = None + rows: list[dict[str, Any]] = [] + total_move = 0.0 + for act in run.actions: + req, res = act["request"], act["response"] + move_m = switch = action_s = 0.0 + detail = "" + if act["path"] in ("/measure", "/clear") and "position" in req: + tgt = (float(req["position"]["x"]), float(req["position"]["y"])) + move_m = math.dist(pos, tgt) + if move_m > 0: + pos = tgt + if act["path"] == "/measure": + ch = int(req.get("channel", 0)) + switch = 0.0 if channel in (None, ch) else SWITCH_S + channel = ch + action_s = MEASURE_S + elif act["path"] == "/clear": + ok = bool(res.get("clear_result") == "cleared") + action_s = CLEAR_OK_S if ok else CLEAR_NO_S + detail = "清除成功" if ok else "未发现" + rows.append(dict( + seq=act["seq"], path=act["path"], move_s=move_m / SPEED_MPS, + switch_s=switch, action_s=action_s, move_m=move_m, + virtual_s=float(res.get("virtual_time_s") or 0.0), detail=detail, + channel=req.get("channel"), x=req.get("position", {}).get("x"), + y=req.get("position", {}).get("y"))) + total_move += move_m + wall = [r for r in run.actions if r.get("wall_ts_ms")] + total_real = (wall[-1]["wall_ts_ms"] - wall[0]["wall_ts_ms"]) / 1000.0 if len(wall) > 1 else 0.0 + return Timeline(rows=rows, total_virtual_s=run.virtual_s, + total_real_s=total_real, total_move_m=total_move) + + +# =========================================================================== +# 数据源选择 +# =========================================================================== +def discover_runs() -> list[LogRun]: + """枚举可回放的日志：优先评测机留档，其次 examples/replay 归档。""" + runs: list[LogRun] = [] + if SESSIONS_DIR.is_dir(): + for log in sorted(SESSIONS_DIR.glob("*/log.jsonl")): + try: + runs.append(parse_jsonl(log)) + except ValueError: + continue + if not runs and ARCHIVE_DIR.is_dir(): + for log in sorted(ARCHIVE_DIR.glob("*.jsonl")): + try: + runs.append(parse_jsonl(log)) + except ValueError: + continue + return runs + + +def pick_run(runs: Sequence[LogRun], *, module: str, truth: int) -> LogRun | None: + """挑一局"完整、干净、点数匹配"的日志：优先全清且动作数适中的那局。""" + cands = [r for r in runs + if r.module.startswith(module) and len(r.truth) >= truth + and r.measures > 0 and r.solution_channels()] + if not cands: + cands = [r for r in runs if r.module.startswith(module)] + if not cands: + return None + return sorted(cands, key=lambda r: (abs(len(r.truth) - truth), len(r.actions)))[0] + + +def q1_practice_case() -> dict[str, Any] | None: + """官方演练问题一留档（输出的定位多边形/圆都由当时的求解器算出）。""" + for path in sorted((ROOT / "output" / "q1" / "practice").glob("*_observations.json")): + doc = json.loads(path.read_text(encoding="utf-8")) + if doc.get("result"): + doc["_path"] = path + return doc + return None + + +# =========================================================================== +# 通用画法 +# =========================================================================== +def layer_band(fig: Figure, y: float, h: float, tag: str, label: str, + color: str, *, x: float = 40.0, w: float = 120.0) -> None: + """左侧"层"标签带 + 右侧标题文字。""" + fig.rect(x, y, w, h, fill=color, stroke=None, radius=6.0) + fig.text(x + w / 2, y + h / 2, tag, size=13, font="hei_bold", + color=PALETTE["paper"], anchor="middle", align="middle") + fig.text(x + w + 12, y + 16, label, size=13.5, font="hei_bold", color=color) + + +def box(fig: Figure, x: float, y: float, w: float, h: float, title: str, body: str = "", + *, fill: str | None = None, stroke: str = PALETTE["primary"], + title_font: str = "hei_bold", title_size: float = 12.5, + body_size: float = 10.5, body_lines: int | None = None, + text_color: str | None = None, radius: float = 6.0) -> None: + """圆角框 + 居中标题（+ 可选多行说明）。""" + fig.rect(x, y, w, h, fill=fill if fill is not None else PALETTE["primary_soft"], + stroke=stroke, sw=1.2, radius=radius) + lines = body.split("\n") if body else [] + n = len(lines) + if body: + block = title_size + 4.0 + n * (body_size + 3.0) + ty = y + (h - block) / 2 + title_size * 0.5 + else: + ty = y + h / 2 + fig.text(x + w / 2, ty, title, size=title_size, font=title_font, + color=text_color or PALETTE["ink"], anchor="middle", align="middle") + for i, line in enumerate(lines): + fig.text(x + w / 2, ty + title_size * 0.5 + 4.0 + body_size * 0.5 + i * (body_size + 3.0), + line, size=body_size, color=PALETTE["muted"], anchor="middle", align="middle") + if body_lines is not None: + assert n <= body_lines, f"说明行数 {n} 超过预留 {body_lines}：{title}" + + +def panel_title(fig: Figure, x: float, y: float, text: str, *, size: float = 12.5, + color: str | None = None) -> None: + fig.text(x, y, text, size=size, font="hei_bold", color=color or PALETTE["primary"]) + + +def footer_note(fig: Figure, text: str, *, y: float | None = None) -> None: + fig.text(24, y if y is not None else fig.height - 22, text, size=10, + color=PALETTE["muted"]) + + +def axis_x(fig: Figure, x0: float, x1: float, y: float, label: str, *, + ticks: Sequence[float] | None = None, unit: str = "", size: float = 10.0, + color: str | None = None, fmt: str = "{:.0f}") -> None: + color = color or PALETTE["muted"] + fig.line(x0, y, x1, y, stroke=color, sw=1.0) + fig.text((x0 + x1) / 2, y + 30, label, size=size + 0.5, color=PALETTE["ink"], + anchor="middle") + if ticks: + for t in ticks: + px = x0 + (x1 - x0) * t + fig.line(px, y - 3, px, y + 3, stroke=color, sw=1.0) + fig.text(px, y + 14, (fmt.format(t) if not unit else fmt.format(t) + unit), + size=size, color=color, anchor="middle") + + +# =========================================================================== +# 图 1：总体思路与算法主线 +# =========================================================================== +def fig_flow() -> Figure: + fig = Figure(1180, 700, title="B 题总体思路与算法主线", + caption="四问共享同一套"协议计时 + 有界误差几何 + 信念更新 + 覆盖计划"主线") + fig.source("src/mcm26b/config.py") + fig.source("docs/B题建模分析与实施计划.md") + + # L0 题面抽象 + layer_band(fig, 66, 46, "L0", "题面抽象：把文字约束变成可计算的集合", PALETTE["primary"], w=110) + l0 = [("区域 Ω：半径 1800 m 的圆盘", PALETTE["faint"]), + ("动作耗时：5 m/s · 检测 5 s · 换频 1 s · 清除 5/3 s", PALETTE["faint"]), + ("观测语义：direction / no_signal / near", PALETTE["faint"]), + ("源：10~16 个，频道 1~20，定向源覆盖 ±90°", PALETTE["faint"])] + for i, (t, fill) in enumerate(l0): + box(fig, 168 + i * 250, 62, 236, 54, t, fill=fill, stroke=PALETTE["faint"], + title_font="song", title_size=10.5, text_color=PALETTE["ink"]) + + # L1 三层在线架构 + layer_band(fig, 140, 192, "L1", "在线三层：证据 → 决策 → 执行（依赖只能向下）", + PALETTE["primary"], w=110) + layers = [ + ("证据层", "候选集 K_c：角域交 + 负观测裁剪\n粒子/网格幸存集：位置 × 朝向", PALETTE["primary_soft"]), + ("决策层", "选点：信息增益 / 最坏后验半径\n任务表：检测 · 定位 · 清除 · 滚动 2-opt", PALETTE["accent_soft"]), + ("执行层", "SimulatorClient：串行 + 幂等重试\n本地截止时钟 + 预留 /exit", PALETTE["success_soft"]), + ] + for i, (t, b, fill) in enumerate(layers): + box(fig, 168 + i * 336, 168, 320, 140, t, b, fill=fill, + stroke=PALETTE["primary"]) + if i < 2: + fig.arrow(168 + i * 336 + 322, 238, 168 + (i + 1) * 336 - 4, 238, + stroke=PALETTE["ink"], sw=1.6) + + # 四问主线 + layer_band(fig, 356, 130, "L2", "四问主线：同一内核逐问扩展（复用而不是重写）", + PALETTE["accent"], w=110) + qs = [ + ("问题 1\n定位解算", "半平面交 P=⋂Wᵢ\nD(P) 顶点对枚举\nMEC 半径 ≤ 20 m 即保证清除", + PALETTE["primary_soft"]), + ("问题 2\n第二点选择", "候选集 K₁ = Ω∩W(s₁)∩B(s₁,1500)\n安全候选盘半径 R/(2cosδ) − 1000\n最坏后验半径最小化", + PALETTE["accent_soft"]), + ("问题 3\n全向搜索", "七点覆盖 a=1159.7 m ⇒ 984.5 m < 1000 m\n贝叶斯闭环 + 任务表调度\n全清即证明，立即 /exit", + PALETTE["success_soft"]), + ("问题 4\n混合定向", "forward_half_disk 31 点最小命中集\n位置 × 朝向联合信念\n有限格后备盲清", + PALETTE["warn_soft"]), + ] + for i, (t, b, fill) in enumerate(qs): + box(fig, 168 + i * 250, 384, 236, 96, t, b, fill=fill, + stroke=PALETTE["primary"], title_size=12.0, body_size=9.5) + if i < 3: + fig.arrow(168 + i * 250 + 238, 432, 168 + (i + 1) * 250 - 2, 432, + stroke=PALETTE["muted"], sw=1.2) + + # 三条闭环 + layer_band(fig, 508, 92, "L3", "三条闭环：发现 → 定位 → 清除 → 覆盖记账", + PALETTE["success"], w=110) + loops = [ + ("发现闭环", "覆盖点巡回扫描，任一频道亮灯即登记新源"), + ("定位闭环", "同频道 ≥2 点交会 → 直径/MEC 判可清除"), + ("清除闭环", "clear 成功即销账；未发现转后备盲清"), + ] + for i, (t, b) in enumerate(loops): + box(fig, 168 + i * 336, 530, 320, 52, t + "：" + b, fill=PALETTE["paper"], + stroke=PALETTE["success"], title_font="song", title_size=10.5) + + fig.arrow(600, 486, 600, 528, stroke=PALETTE["muted"], sw=1.2) + footer_note(fig, "口径：计时/几何/覆盖常数均来自 src/mcm26b/config.py；" + "各问实现的核对结果见 docs/algo-1..4.md。") + return fig + + +# =========================================================================== +# 图 2：执行步骤状态机（协议回合） +# =========================================================================== +def fig_steps() -> Figure: + fig = Figure(1100, 560, title="执行步骤：接口状态机与一次动作回合", + caption="题目没有显式协议；实际系统是"窗口驱动 + 严格串行"的动作回合") + fig.source("docs/协议规则表与状态机.md") + fig.source("src/mcm26b/protocol/client.py") + + # 上：会话状态机 + panel_title(fig, 40, 44, "(a) 会话状态机：窗口开放与截止") + states = [("IDLE", 90, PALETTE["faint"]), ("COUNTDOWN 5 s", 300, PALETTE["warn_soft"]), + ("OPEN：接口开放", 560, PALETTE["success_soft"]), + ("CLOSED：断连", 830, PALETTE["accent_soft"])] + for i, (t, x, fill) in enumerate(states): + box(fig, x, 72, 180, 52, t, fill=fill, stroke=PALETTE["muted"], title_size=11.5) + if i < 3: + fig.arrow(x + 182, 98, states[i + 1][1] - 4, 98, stroke=PALETTE["ink"], sw=1.4) + fig.text(560, 138, "/enter 成功 ⇒ 现实预算 1200 s 起算；" + "虚拟世界限时 360000 s（两者独立）", + size=10.5, color=PALETTE["muted"], anchor="middle") + fig.text(830, 156, "剩 60 s 预留 /exit", size=10, color=PALETTE["accent"], anchor="middle") + + # 下：单次动作回合 + panel_title(fig, 40, 206, "(b) 单次动作回合：一次请求的完整生命周期") + steps = [ + ("本地预校验", "channel ∈ [1,20]\n坐标有限且 ≤2e6\nrequest_id 格式", PALETTE["faint"]), + ("串行发送", "同一时刻只有\n一个未决请求\n（不并发）", PALETTE["primary_soft"]), + ("等服务端计时", "accepted=true 才更新\n虚拟时刻以响应为准\naccepted=false 不改状态", PALETTE["warn_soft"]), + ("更新本地状态", "位置 / 当前频道\n剩余现实预算\nJSONL 会话日志", PALETTE["success_soft"]), + ("幂等重试（异常路径）", "网络错 / 429 / 500：\n同 ID 同内容原样重发\n400/409 立即失败", PALETTE["accent_soft"]), + ] + x0, w, gap = 44, 186, 26 + for i, (t, b, fill) in enumerate(steps): + x = x0 + i * (w + gap) + box(fig, x, 240, w, 104, t, b, fill=fill, stroke=PALETTE["muted"], + title_size=11.0, body_size=9.5) + if i < len(steps) - 1: + fig.arrow(x + w + 2, 292, x + w + gap - 4, 292, stroke=PALETTE["ink"], sw=1.3) + + # 计时规则表 + panel_title(fig, 40, 382, "(c) 每个动作的虚拟耗时（附件1 表1 / 附件2 第4节）") + rows = [ + ("/enter · /exit", "0 s", "不改变虚拟时刻，/exit 不切频道"), + ("移动（/measure、/clear 携带 position）", "直线距离 / 5 m·s⁻¹", "同点不动则 0 s"), + ("切换测向频道（仅 /measure 触发）", "1 s", "/clear 不切换频道"), + ("检测 /measure", "5 s", "无论返回 direction 还是 no_signal"), + ("清除 /clear 成功 / 未发现", "5 s / 3 s", "5 s = 光学定位 3 s + 激光清除 2 s"), + ] + top, rh = 412, 26 + fig.rect(44, top, 1012, rh, fill=PALETTE["primary_soft"], stroke=None) + for j, head in enumerate(("指令 / 动作", "虚拟耗时", "说明")): + fig.text([64, 560, 700][j], top + rh / 2, head, size=10.5, font="hei_bold", + align="middle") + for i, (a, b, c) in enumerate(rows): + y = top + rh * (i + 1) + if i % 2 == 0: + fig.rect(44, y, 1012, rh, fill="#f7f9fb", stroke=None) + fig.text(64, y + rh / 2, a, size=10.5, align="middle") + fig.text(560, y + rh / 2, b, size=10.5, font="times_bold", color=PALETTE["accent"], + align="middle") + fig.text(700, y + rh / 2, c, size=10, color=PALETTE["muted"], align="middle") + footer_note(fig, "总分：T = L/5 + 5N_m + N_s + 5N_c⁺ + 3N_c⁻（虚拟时间即成绩口径）") + return fig + + +# =========================================================================== +# 图 3：附件计时示例逐步推演 +# =========================================================================== +EXAMPLE_STEPS = [ + ("/enter", 0.0, 0.0, 0.0, "初始频道 1，虚拟时刻 0"), + ("/measure (300,400) ch1", 100.0, 0.0, 5.0, "移动 500 m；频道未变"), + ("/measure (300,400) ch2", 0.0, 1.0, 5.0, "换频 1 s，位置不变"), + ("/clear (300,0) ch3", 80.0, 0.0, 3.0, "未发现目标；/clear 不切频道"), + ("/measure (300,0) ch2", 0.0, 0.0, 5.0, "频道仍是 2（第 3 步之后未变）"), + ("/exit", 0.0, 0.0, 0.0, "退出，虚拟时刻停在 199 s"), +] + + +def fig_timing_rule() -> Figure: + fig = Figure(1160, 470, title="附件计时示例逐步推演（105 / 111 / 194 / 199）", + caption="协议计时的唯一权威例子：每一步的分项与累计都必须逐项复现") + fig.source("docs/协议规则表与状态机.md") + + # 左：堆叠条 + panel_title(fig, 40, 44, "(a) 每一步的耗时构成与累计虚拟时刻") + x0, w, gap, base = 60, 120, 34, 300 + scale = 0.62 # px per second + acc = 0.0 + for i, (cmd, mv, sw, ac, note) in enumerate(EXAMPLE_STEPS): + x = x0 + i * (w + gap) + total = mv + sw + ac + y = base + for value, color in ((mv, PALETTE["series"][0]), (sw, PALETTE["series"][3]), + (ac, PALETTE["series"][2])): + h = value * scale + if h > 0: + fig.rect(x, y - h, w, h, fill=color, stroke=None) + fig.text(x + w / 2, y - h / 2, f"{value:.0f}", size=10.5, + font="times_bold", color=PALETTE["paper"], anchor="middle", + align="middle") + y -= h + acc += total + fig.text(x + w / 2, base + 20, cmd.split(" ")[0], size=10.5, font="times_bold", + anchor="middle") + fig.text(x + w / 2, base + 36, (cmd.split(" ", 1)[1] if " " in cmd else ""), + size=9, color=PALETTE["muted"], anchor="middle") + fig.text(x + w / 2, base - total * scale - 18, f"Σ={acc:.0f}", size=11, + font="times_bold", color=PALETTE["accent"], anchor="middle") + if i: + fig.arrow(x - gap + 4, base + 8, x - 6, base + 8, stroke=PALETTE["faint"], + sw=1.0, head=6) + fig.line(x0, base + 2, x0 + 6 * (w + gap) - gap, base + 2, stroke=PALETTE["muted"], sw=1.2) + fig.legend_chips(60, 70, [("移动 L/5", PALETTE["series"][0]), + ("换频 1 s", PALETTE["series"][3]), + ("动作 5 s / 3 s", PALETTE["series"][2])], size=10.5) + + # 右：核对表 + panel_title(fig, 700, 44, "(b) 与附件原文逐步核对") + top, rh = 66, 30 + heads = ("步", "指令", "本步", "累计", "关键点") + widths = (34, 250, 60, 60, 316) + xs = [700] + for wd in widths[:-1]: + xs.append(xs[-1] + wd) + fig.rect(700, top, sum(widths), rh, fill=PALETTE["primary_soft"], stroke=None) + for x, head in zip(xs, heads): + fig.text(x + 8, top + rh / 2, head, size=10.5, font="hei_bold", align="middle") + acc = 0.0 + for i, (cmd, mv, sw, ac, note) in enumerate(EXAMPLE_STEPS): + y = top + rh * (i + 1) + acc += mv + sw + ac + if i == 3: + fig.rect(700, y, sum(widths), rh, fill=PALETTE["warn_soft"], stroke=None) + elif i % 2 == 0: + fig.rect(700, y, sum(widths), rh, fill="#f7f9fb", stroke=None) + fig.text(xs[0] + 8, y + rh / 2, str(i + 1), size=10.5, align="middle") + fig.text(xs[1] + 8, y + rh / 2, cmd, size=10, align="middle") + fig.text(xs[2] + 8, y + rh / 2, f"{mv + sw + ac:.0f}", size=10.5, font="times_bold", + align="middle") + fig.text(xs[3] + 8, y + rh / 2, f"{acc:.0f}", size=10.5, font="times_bold", + color=PALETTE["accent"], align="middle") + fig.text(xs[4] + 8, y + rh / 2, note, size=9.5, color=PALETTE["muted"], + align="middle") + fig.text(700, top + rh * 7 + 18, "验收：uv run mcm26b sim-check 复现 [105, 111, 194, 199]", + size=10.5, color=PALETTE["success"]) + fig.text(700, top + rh * 7 + 38, + "tests/test_protocol_timing.py 逐步断言移动 / 换频 / 动作三个分项", + size=10, color=PALETTE["muted"]) + footer_note(fig, "注：第 4 步 /clear 未发现目标，只计入 3 s 光学定位；" + "第 5 步频道仍为 2，因此不产生换频耗时。") + return fig + + +# =========================================================================== +# 图 4：实测日志的虚拟时间分解 +# =========================================================================== +def fig_timing_run(run: LogRun) -> Figure: + tl = decompose(run) + fig = Figure(1160, 600, title=f"实测执行时序分解：{run.case_code}（{run.module}，" + f"{len(run.truth)} 源）", + caption="同一局的虚拟时间构成、累计曲线与现实耗时对照") + fig.source(run.source, case_code=run.case_code, module=run.module, + actions=len(run.actions)) + + left, right = 70, 1120 + mid = 560 + # (a) 累计虚拟时间 + panel_title(fig, 40, 44, "(a) 累计虚拟时间（按动作类型堆叠）") + top, bottom = 76, 300 + tv = max(tl.total_virtual_s, 1.0) + acc = {"move": 0.0, "switch": 0.0, "action": 0.0} + xs, series = [], [] + for row in tl.rows: + acc["move"] += row["move_s"] + acc["switch"] += row["switch_s"] + acc["action"] += row["action_s"] + xs.append(row["seq"]) + series.append(dict(acc)) + n = max(len(xs), 1) + + def px(seq: float) -> float: + return left + (right - left) * (seq - 1) / max(n - 1, 1) + + def py(value: float) -> float: + return bottom - (bottom - top) * (value / (tv * 1.08)) + + for f in (0.25, 0.5, 0.75, 1.0): + y = py(tv * f) + fig.line(left, y, right, y, stroke=PALETTE["faint"], sw=0.8) + fig.text(left - 8, y, f"{tv * f:.0f}", size=10, font="times", anchor="end", + align="middle") + fig.text(left - 52, (top + bottom) / 2, "虚拟时间 / s", size=10.5, anchor="middle", + rotate=-90) + layers = ("move", "switch", "action") + colors = {"move": PALETTE["series"][0], "switch": PALETTE["series"][3], + "action": PALETTE["series"][2]} + for y_prev in (0.0,): + pass + for i, key in enumerate(layers): + upper = [s[key] + sum(s[k] for k in layers[i + 1:]) for s in series] + lower = [sum(s[k] for k in layers[i + 1:]) for s in series] + pts = [(px(x), py(u)) for x, u in zip(xs, upper)] + \ + [(px(x), py(l)) for x, l in zip(reversed(xs), reversed(lower))] + fig.polygon(pts, fill=colors[key], stroke=None, opacity=0.9) + fig.line(left, bottom, right, bottom, stroke=PALETTE["muted"], sw=1.2) + fig.legend_chips(72, 62, [("移动 L/5", colors["move"]), ("换频 1 s", colors["switch"]), + ("检测/清除动作", colors["action"])], size=10.5) + fig.text(right, bottom + 22, f"动作序号 1 … {n}", size=10, anchor="end") + + # (b) 分项占比 + panel_title(fig, 40, 334, "(b) 虚拟时间构成占比") + totals = {k: sum(s[k] for s in series) for k in layers} + bar_x, bar_y, bar_w, bar_h = 70, 366, 420, 26 + cur = bar_x + for key in ("move", "action", "switch"): + seg = bar_w * (totals[key] / tv) + fig.rect(cur, bar_y, seg, bar_h, fill=colors[key], stroke=None) + if seg > 34: + fig.text(cur + seg / 2, bar_y + bar_h / 2, f"{totals[key] / tv * 100:.0f}%", + size=10.5, font="times_bold", color=PALETTE["paper"], + anchor="middle", align="middle") + cur += seg + labels = [("移动", "move"), ("检测/清除", "action"), ("换频", "switch")] + ly = bar_y + bar_h + 22 + for label, key in labels: + fig.rect(bar_x, ly - 5, 10, 10, fill=colors[key], stroke=None, radius=2) + fig.text(bar_x + 16, ly, f"{label} {totals[key]:,.0f} s", size=10.5, + align="middle") + bar_x += 140 + fig.text(70, bar_y + bar_h + 52, + f"合计 {tv:,.1f} s ＝ 移动 {totals['move']:,.0f} + 换频 {totals['switch']:,.0f} " + f"+ 动作 {totals['action']:,.0f}", size=10.5, color=PALETTE["ink"]) + fig.text(70, bar_y + bar_h + 74, + f"检测 {run.measures} 次 · 清除 {run.clears} 次 · " + f"移动总里程 {tl.total_move_m / 1000:,.1f} km · " + f"清除频道 {len(run.solution_channels())}/{len(run.truth)}", + size=10.5, color=PALETTE["muted"]) + + # (c) 现实耗时 + panel_title(fig, 40, 470, "(c) 现实耗时对照") + panel_title(fig, 620, 470, "(d) 逐动作耗时（前 40 个动作）") + rw = tl.total_real_s + real_x, real_y, real_w = 70, 506, 420 + ratio = tv / rw if rw > 0 else 0.0 + fig.rect(real_x, real_y, real_w, 20, fill=PALETTE["primary_soft"], stroke=None) + fig.rect(real_x, real_y, real_w * min(1.0, 1.0), 20, fill=PALETTE["primary"], + stroke=None) + fig.text(real_x, real_y + 40, + f"虚拟 {tv:,.1f} s ÷ 现实 {rw:.2f} s ≈ {ratio:,.0f} 倍加速", + size=11, font="times_bold", color=PALETTE["accent"]) + fig.text(real_x, real_y + 60, + "现实耗时只受网络与服务端开销影响：本地回环约毫秒级，" + "因此成绩口径恒为虚拟时间。", size=9.5, color=PALETTE["muted"]) + + bx, by, bw, bh = 620, 494, 500, 76 + subset = tl.rows[:40] + mx = max((r["move_s"] + r["switch_s"] + r["action_s"]) for r in subset) or 1.0 + for i, row in enumerate(subset): + x = bx + bw * i / max(len(subset), 1) + wpx = bw / max(len(subset), 1) * 0.72 + y = by + bh + for value, key in ((row["move_s"], "move"), (row["switch_s"], "switch"), + (row["action_s"], "action")): + h = bh * value / mx + if h > 0: + fig.rect(x, y - h, wpx, h, fill=colors[key], stroke=None) + y -= h + fig.line(bx, by + bh, bx + bw, by + bh, stroke=PALETTE["muted"], sw=1.0) + fig.text(bx + bw / 2, by + bh + 24, "动作序号（1…40）", size=10, anchor="middle") + fig.text(bx + bw + 6, by + 6, f"峰值 {mx:.1f} s", size=9.5, color=PALETTE["muted"]) + footer_note(fig, "复算口径与 referee.simulator 同源：移动=直线距离/5 m·s⁻¹、" + "换频 1 s、检测 5 s、清除成功 5 s / 未发现 3 s。") + return fig + + +# =========================================================================== +# 图 5/6：问题一几何 +# =========================================================================== +def _world(fig: Figure, x: float, y: float, w: float, h: float, half: float, *, + grid: float = 250.0) -> tuple[Any, Any]: + """建立世界坐标 → 画布坐标的映射（等比例、y 轴向上）。""" + def sx(wx: float) -> float: + return x + w / 2 + wx / half * (w / 2) + + def sy(wy: float) -> float: + return y + h / 2 - wy / half * (h / 2) + + def poly(points: Iterable[Sequence[float]]) -> list[list[float]]: + return [[sx(px), sy(py)] for px, py in points] + + fig.rect(x, y, w, h, fill="#fbfcfd", stroke=PALETTE["faint"], sw=1.0) + step = grid + while step < half: + fig.line(sx(-half), sy(step), sx(half), sy(step), stroke="#eef2f6", sw=0.8) + fig.line(sx(-half), sy(-step), sx(half), sy(-step), stroke="#eef2f6", sw=0.8) + fig.line(sx(step), sy(-half), sx(step), sy(half), stroke="#eef2f6", sw=0.8) + fig.line(sx(-step), sy(-half), sx(-step), sy(half), stroke="#eef2f6", sw=0.8) + step *= 2 + fig.line(sx(-half), sy(0), sx(half), sy(0), stroke=PALETTE["faint"], sw=1.0) + fig.line(sx(0), sy(-half), sx(0), sy(half), stroke=PALETTE["faint"], sw=1.0) + return sx, sy + + +def _wedge(fig: Figure, sx: Any, sy: Any, station: Sequence[float], bearing_deg: float, + half_deg: float, reach: float, *, color: str, opacity: float = 0.28, + sw: float = 1.2, dash: str | None = None) -> None: + """前向角域扇形：从检测点沿示向度 ±half 度、长度 reach。""" + a0 = math.radians(bearing_deg - half_deg) + a1 = math.radians(bearing_deg + half_deg) + pts = [station] + steps = 12 + for i in range(steps + 1): + a = a0 + (a1 - a0) * i / steps + pts.append([station[0] + reach * math.cos(a), station[1] + reach * math.sin(a)]) + fig.polygon([[sx(px), sy(py)] for px, py in pts], fill=color, stroke=None, + opacity=opacity) + for a in (a0, a1): + fig.line(sx(station[0]), sy(station[1]), + sx(station[0] + reach * math.cos(a)), sy(station[1] + reach * math.sin(a)), + stroke=color, sw=sw, dash=dash) + fig.line(sx(station[0]), sy(station[1]), + sx(station[0] + reach * math.cos(math.radians(bearing_deg))), + sy(station[1] + reach * math.sin(math.radians(bearing_deg))), + stroke=color, sw=sw + 0.4) + + +def fig_q1_meetup() -> Figure: + case = q1_practice_case() + fig = Figure(1160, 540, title="问题一：角域交会、直径与最小包围圆", + caption="半平面交给出定位区域 P；直径与最小包围圆决定"能否保证清除"") + if case: + title = f"官方演练案例 {case['case_code']}（频道 {case['channel']}，三点示向度 ±1°）" + source = case["_path"] + else: + title = "反例算例 examples/q1/triangle.json（三点示向度 ±1°）" + source = ROOT / "examples" / "q1" / "triangle.json" + panel_title(fig, 40, 44, "(a) " + title) + + half = 60.0 + gx, gy, gw, gh = 60, 70, 560, 400 + sx, sy = _world(fig, gx, gy, gw, gh, half, grid=20.0) + if case: + stations = [(m["x"], m["y"], m["bearing_deg"], m["half_angle_deg"]) + for m in case["measurements"]] + verts = case["result"]["vertices"] + d_pair = case["result"]["diameter_pair"] + dc = case["result"]["diameter_circle"] + mec = case["result"]["minimum_enclosing_circle"] + mec_r = mec["radius_m"] + E = dc["max_vertex_distance_m"] - dc["radius_m"] + D = case["result"]["diameter_m"] + else: + doc = json.loads(Path(source).read_text(encoding="utf-8")) + stations = [(m["x"], m["y"], m["bearing_deg"], m.get("half_angle_deg", 1.0)) + for m in doc["measurements"]] + sol = _solve_q1(stations) + verts, d_pair = sol["vertices"], sol["diameter_pair"] + dc, mec = sol["diameter_circle"], sol["mec"] + mec_r = mec["radius_m"] + E = dc["max_vertex_distance_m"] - dc["radius_m"] + D = math.dist(*d_pair) + for s_x, s_y, bearing, hd in stations: + _wedge(fig, sx, sy, (s_x, s_y), bearing, hd, 150.0, + color=PALETTE["series"][0], opacity=0.14) + fig.circle(sx(s_x), sy(s_y), 4.0, fill=PALETTE["primary"], stroke=None) + fig.text(sx(s_x) + 7, sy(s_y) - 8, f"({s_x:.0f}, {s_y:.0f})", size=9, + color=PALETTE["muted"], font="times") + poly = [[sx(px), sy(py)] for px, py in verts] + fig.polygon(poly, fill=PALETTE["accent_soft"], stroke=PALETTE["accent"], sw=1.4, + opacity=0.95) + fig.circle(sx(dc["center"][0]), sy(dc["center"][1]), + dc["radius_m"] * gw / (2 * half), fill=None, + stroke=PALETTE["series"][4], sw=1.2, dash="5 4") + fig.circle(sx(mec["center"][0]), sy(mec["center"][1]), + mec["radius_m"] * gw / (2 * half), fill=None, + stroke=PALETTE["success"], sw=1.6) + fig.circle(sx(mec["center"][0]), sy(mec["center"][1]), 3.5, fill=PALETTE["success"], + stroke=None) + fig.legend_chips(74, 486, [("前向角域 Wᵢ", PALETTE["series"][0]), + ("定位区域 P", PALETTE["accent"]), + ("直径圆（不覆盖）", PALETTE["series"][4]), + ("最小包围圆 MEC", PALETTE["success"])], size=10) + fig.text(gx + gw - 6, gy + gh - 14, "单位：米", size=9.5, color=PALETTE["muted"], + anchor="end") + + # 右栏：结论 + rx = 660 + panel_title(fig, rx, 44, "(b) 关键量与判定") + rows = [ + ("定位区域状态", "多边形（有界）"), + ("直径 D(P)", f"{ctx['D']:.4f} m"), + ("最小包围圆半径 ρ", f"{ctx['rho']:.4f} m"), + ("直径圆最大超出 E", f"{ctx['E']:+.4f} m"), + ("直径圆是否覆盖 P", "否" if ctx["E"] > 0 else "是"), + ("ρ ≤ 20 m ⇒ 可保证清除", "是" if ctx["rho"] <= 20 else "否"), + ] + top, rh = 76, 34 + for i, (k, v) in enumerate(rows): + y = top + i * rh + fill = "#f7f9fb" if i % 2 == 0 else None + if fill: + fig.rect(rx, y, 460, rh, fill=fill, stroke=None) + fig.text(rx + 10, y + rh / 2, k, size=11, align="middle") + fig.text(rx + 450, y + rh / 2, v, size=11.5, font="times_bold", + color=PALETTE["accent"] if i in (3, 4, 5) else PALETTE["ink"], + anchor="end", align="middle") + box(fig, rx, top + 6 * rh + 16, 460, 96, + "两个证明要点", + "① 直径只需枚举顶点对：凸集上最远点必为顶点；\n" + "② 直径圆唯一可能圆心是直径中点：若 a、b 同在该圆内，\n" + " 三角不等式两处同时取等 ⇒ 圆心只能是中点。", + fill=PALETTE["primary_soft"], stroke=PALETTE["primary"], title_size=11.5, + body_size=10) + box(fig, rx, top + 6 * rh + 124, 460, 86, "结论", + "直径圆不覆盖定位区域（E > 0），题目要求的是\n" + "\"存在一点可保证清除\"：取 MEC 圆心即可，\n" + "因为 ρ ≤ 20 m 落在清除半径内。", + fill=PALETTE["success_soft"], stroke=PALETTE["success"], title_size=11.5, + body_size=10) + fig.source(source) + footer_note(fig, "求解入口：uv run mcm26b q1 .json；" + "算法与退化情形见 docs/algo-1.md。") + return fig + + +def _solve_q1(stations: Sequence[Sequence[float]]) -> dict[str, Any]: + """复用生产实现求解（避免图上出现与代码不一致的手算数字）。 + + 返回 ``{"vertices", "diameter_pair", "diameter_circle", "mec"}``， + 其中两个圆都是 ``{"center": [x, y], "radius_m": r}``。 + """ + sys.path.insert(0, str(ROOT / "src")) + from mcm26b import Measurement, solve # noqa: PLC0415 + ms = [Measurement(x, y, b, h) for x, y, b, h in stations] + res = solve(ms) + dc, mec = res["diameter_circle"], res["minimum_enclosing_circle"] + return dict(vertices=res["vertices"], diameter_pair=res["diameter_pair"], + diameter_circle=dict(center=dc["center"], radius_m=dc["radius_m"], + max_vertex_distance_m=dc["max_vertex_distance_m"]), + mec=dict(center=mec["center"], radius_m=mec["radius_m"])) + + +def fig_q1_counterexample() -> Figure: + """题目要的"直径圆不覆盖"反例：等边三角形三点观测。""" + doc = json.loads((ROOT / "examples" / "q1" / "triangle.json").read_text(encoding="utf-8")) + stations = [(m["x"], m["y"], m["bearing_deg"], m.get("half_angle_deg", 1.0)) + for m in doc["measurements"]] + verts, pair, dc, mec = _solve_q1(stations) + fig = Figure(880, 520, title="问题一：直径圆可以覆盖吗？——符合题设的反例", + caption="三点等边布站、示向度精度 ±1°，定位区域仍是"芒果形"四边形") + fig.source("examples/q1/triangle.json") + half = 1200.0 + gx, gy, gw, gh = 60, 60, 760, 400 + sx, sy = _world(fig, gx, gy, gw, gh, half, grid=400.0) + for s_x, s_y, bearing, hd in stations: + _wedge(fig, sx, sy, (s_x, s_y), bearing, hd, 1500.0, + color=PALETTE["series"][0], opacity=0.10) + fig.circle(sx(s_x), sy(s_y), 4.5, fill=PALETTE["primary"], stroke=None) + fig.text(sx(s_x) + 8, sy(s_y) - 6, f"({s_x:.0f}, {s_y:.0f})", size=9.5, + color=PALETTE["primary"], font="times") + fig.polygon([[sx(px), sy(py)] for px, py in verts], fill=PALETTE["accent_soft"], + stroke=PALETTE["accent"], sw=1.6) + fig.circle(sx(dc[0][0]), sy(dc[0][1]), dc[1] * gw / (2 * half), fill=None, + stroke=PALETTE["series"][4], sw=1.4, dash="6 4") + fig.circle(sx(mec[0][0]), sy(mec[0][1]), mec[1] * gw / (2 * half), fill=None, + stroke=PALETTE["success"], sw=1.6) + for px, py in pair: + fig.circle(sx(px), sy(py), 4.0, fill=PALETTE["accent"], stroke=None) + fig.arrow(sx(dc[0][0]), sy(dc[0][1]), sx(verts[1][0]), sy(verts[1][1]), + stroke=PALETTE["accent"], sw=1.0) + fig.text(sx(verts[1][0]) + 6, sy(verts[1][1]) + 14, "最大超出顶点", size=9.5, + color=PALETTE["accent"]) + fig.legend_chips(74, 486, [("定位区域 P", PALETTE["accent"]), + ("直径圆（不覆盖）", PALETTE["series"][4]), + ("最小包围圆 MEC", PALETTE["success"])], size=10) + footer_note(fig, f"D(P) = {math.dist(*pair):.4f} m，MEC 半径 {mec[1]:.4f} m；" + "反例构造与证明见 docs/问题一解答.md §六。") + return fig + + +# =========================================================================== +# 图 7：问题二保证接收盘 +# =========================================================================== +def _q2_params() -> dict[str, float]: + sys.path.insert(0, str(ROOT / "src")) + from mcm26b.config import ProblemConfig # noqa: PLC0415 + cfg = ProblemConfig() + delta = math.radians(cfg.field.bearing_error_half_deg) + r_max = cfg.field.r_eff_max_m + a = r_max / (2 * math.cos(delta)) # 截断角域最小包围圆半径 + safe_r = R_EFF - a # 保证接收盘半径 + d0 = (R_EFF + r_max) / 2 # 最坏后验半径的直观位置 + return dict(delta_deg=cfg.field.bearing_error_half_deg, a=a, safe_r=safe_r, + center_d=a, d0=d0 - a, reach=r_max) + + +def fig_q2_region() -> Figure: + p = _q2_params() + fig = Figure(1120, 540, title="问题二：保证接收的候选盘与第二点选择", + caption="候选集 K₁ 由三个闭集相交给出；解析安全盘中心在 s₁ + a·u") + fig.source("src/mcm26b/config.py") + half = 1700.0 + gx, gy, gw, gh = 60, 70, 620, 420 + sx, sy = _world(fig, gx, gy, gw, gh, half, grid=500.0) + # Ω 边界 + fig.circle(sx(0), sy(0), OMEGA_R * gw / (2 * half), fill=None, + stroke=PALETTE["primary"], sw=1.4) + fig.text(sx(0), sy(OMEGA_R) - 14, "Ω：R = 1800 m", size=10, + color=PALETTE["primary"], anchor="middle") + # 首次检测点与角域 + s1 = (0.0, 0.0) + bearing = 30.0 + _wedge(fig, sx, sy, s1, bearing, 1.0, 1500.0, color=PALETTE["accent"], opacity=0.25, + sw=1.4) + fig.circle(sx(*s1), sy(*s1), 5.0, fill=PALETTE["accent"], stroke=None) + fig.text(sx(*s1) + 8, sy(*s1) - 10, "s₁（首次检测点）", size=10.5, + color=PALETTE["accent"]) + a = p["a"] + cx = a * math.cos(math.radians(bearing)) + cy = a * math.sin(math.radians(bearing)) + fig.circle(sx(cx), sy(cy), a * gw / (2 * half), fill=None, stroke=PALETTE["series"][4], + sw=1.2, dash="6 4") + fig.circle(sx(cx), sy(cy), R_EFF * gw / (2 * half), fill=PALETTE["success_soft"], + stroke=PALETTE["success"], sw=1.6, opacity=0.45) + fig.circle(sx(cx), sy(cy), 4.0, fill=PALETTE["success"], stroke=None) + fig.arrow(sx(*s1), sy(*s1), sx(cx), sy(cy), stroke=PALETTE["series"][4], sw=1.3) + fig.text((sx(*s1) + sx(cx)) / 2, (sy(*s1) + sy(cy)) / 2 - 10, + f"a = R_max/(2cos δ) = {a:.2f} m", size=10, color=PALETTE["series"][4]) + fig.legend_chips(74, 508, [("首次角域 W(s₁)", PALETTE["accent"]), + ("截断角域的外接圆 B(c,a)", PALETTE["series"][4]), + ("安全候选盘 B(c, 1000−a)", PALETTE["success"])], size=10) + + # 右栏：解析结论 + rx = 720 + panel_title(fig, rx, 44, "解析结论（δ = ±1°）") + rows = [ + ("首次角域截断半径 R_max", "1500 m"), + ("a = R_max / (2cos δ)", f"{a:.4f} m"), + ("安全盘半径 1000 − a", f"{p['safe_r']:.4f} m"), + ("安全盘圆心 c = s₁ + a·u", f"({cx:.2f}, {cy:.2f})"), + ] + top, rh = 76, 32 + for i, (k, v) in enumerate(rows): + y = top + i * rh + if i % 2 == 0: + fig.rect(rx, y, 360, rh, fill="#f7f9fb", stroke=None) + fig.text(rx + 10, y + rh / 2, k, size=10.5, align="middle") + fig.text(rx + 350, y + rh / 2, v, size=11, font="times_bold", + color=PALETTE["accent"], anchor="end", align="middle") + box(fig, rx, top + 4 * rh + 14, 360, 118, "为什么它是"安全"的", + "截断角域（扇形）的最小包围圆恰为 B(c, a)，\n" + "故扇形 ⊂ B(c, a)。取 B(c, 1000−a) 内任一点，\n" + "与扇形内任一点的距离 ≤ a + (1000−a) = 1000 m，\n" + "即第二点必然落在源的 1000 m 内 ⇒ 必然收到信号。", + fill=PALETTE["primary_soft"], stroke=PALETTE["primary"], title_size=11.5, + body_size=9.5) + box(fig, rx, top + 4 * rh + 142, 360, 96, "选点目标", + "在安全盘内最小化"最坏保守后验包围半径"：\n" + "主目标 = 后验区域直径，次目标 = 移动 + 检测成本；\n" + "交会角只作同分筛选与可解释性说明。", + fill=PALETTE["success_soft"], stroke=PALETTE["success"], title_size=11.5, + body_size=9.5) + footer_note(fig, f"δ = ±1° 时安全盘半径约 {p['safe_r']:.2f} m（严格非空）；" + "完整推导见 docs/algo-2.md §3。") + return fig + + +# =========================================================================== +# 图 8：问题三覆盖计划 +# =========================================================================== +def fig_q3_coverage() -> Figure: + fig = Figure(1120, 600, title="问题三：七点解析覆盖与"发现即证明"", + caption="七点计划用解析上界换掉网格：任意源都在某站点 1000 m 内") + fig.source("docs/algo-3.md") + fig.source("src/mcm26b/planner/seven.py") + half = 2200.0 + gx, gy, gw, gh = 60, 70, 620, 440 + sx, sy = _world(fig, gx, gy, gw, gh, half, grid=600.0) + fig.circle(sx(0), sy(0), OMEGA_R * gw / (2 * half), fill="#f6f9fc", + stroke=PALETTE["primary"], sw=1.4) + ring = [(RING_RADIUS_M * math.cos(k * math.pi / 3), RING_RADIUS_M * math.sin(k * math.pi / 3)) + for k in range(6)] + fig.circle(sx(0), sy(0), RING_RADIUS_M * gw / (2 * half), fill=None, + stroke=PALETTE["faint"], sw=1.0, dash="4 4") + for px, py in ring: + fig.circle(sx(px), sy(py), 5.0, fill=PALETTE["accent"], stroke=None) + fig.circle(sx(0), sy(0), 5.0, fill=PALETTE["accent"], stroke=None) + fig.text(sx(0), sy(0) - 12, "S₀ 原点", size=10, color=PALETTE["accent"], + anchor="middle") + fig.text(sx(0), sy(RING_RADIUS_M) - 16, f"环半径 a = {RING_RADIUS_M:.1f} m", + size=10, color=PALETTE["accent"], anchor="middle") + fig.circle(sx(0), sy(0), R_EFF * gw / (2 * half), fill=PALETTE["success_soft"], + stroke=PALETTE["success"], sw=1.2, opacity=0.35) + fig.text(sx(0) + 6, sy(R_EFF) - 16, "1000 m 接收半径", size=10, + color=PALETTE["success"]) + d_bound = max(RING_RADIUS_M / math.sqrt(3), math.sqrt( + SEVEN_R ** 2 + RING_RADIUS_M ** 2 - math.sqrt(3) * SEVEN_R * RING_RADIUS_M)) + fig.legend_chips(74, 528, [("站点（原点 + 六等分环）", PALETTE["accent"]), + ("1000 m 接收盘", PALETTE["success"])], size=10) + fig.text(74, 556, f"解析上界 D(a) = max{{a/√3, √(R²+a²−√3Ra)}} = {d_bound:.4f} m < 1000 m " + f"⇒ 15.46 m 几何余量", size=10.5, color=PALETTE["success"]) + + # 右栏：对比 + rx = 720 + panel_title(fig, rx, 44, "(a) 三套发现计划的代价对比") + heads = ("计划", "点数", "保证", "每轮检测") + widths = (120, 56, 116, 88) + top, rh = 76, 30 + x = rx + fig.rect(rx, top, sum(widths), rh, fill=PALETTE["primary_soft"], stroke=None) + for wd, head in zip(widths, heads): + fig.text(x + wd / 2, top + rh / 2, head, size=10.5, font="hei_bold", + anchor="middle", align="middle") + x += wd + rows = [ + ("七点解析（默认）", "7", "a=1159.7 m ⇒ 984.54 m", "7 × 频道数"), + ("49 点网格（基线）", "49", "700√2 = 989.95 m", "49 × 频道数"), + ("最小命中集（题4）", "31", "990.05 m（定向）", "31 × 频道数"), + ] + for i, row in enumerate(rows): + y = top + rh * (i + 1) + x = rx + if i == 0: + fig.rect(rx, y, sum(widths), rh, fill=PALETTE["success_soft"], stroke=None) + elif i % 2 == 0: + fig.rect(rx, y, sum(widths), rh, fill="#f7f9fb", stroke=None) + for wd, cell in zip(widths, row): + fig.text(x + wd / 2, y + rh / 2, cell, size=10, anchor="middle", align="middle") + x += wd + + panel_title(fig, rx, 208, "(b) 为什么可以"发现即证明"") + box(fig, rx, 236, 380, 120, + "覆盖账本 + 负观测", + "每个频道维护"是否已被某站点在 1000 m 内测过"：\n" + "· 收到 direction ⇒ 该源已发现，转入定位/清除；\n" + "· 收到 no_signal ⇒ 只能排除以该站点为心、\n" + " 半径 1000 m 的圆盘（r_eff ≥ 1000 m）；\n" + "· 当某频道在全部站点都测过且全为 no_signal，\n" + " 由覆盖证明知该频道必无源 ⇒ 可安全销账。", + fill=PALETTE["primary_soft"], stroke=PALETTE["primary"], title_size=11.5, + body_size=9.5) + box(fig, rx, 372, 380, 138, "闭环与停止", + "发现 → 同频道第二点交会 → ρ ≤ 20 m 判定 → clear；\n" + "未发现则换点重测，必要时转后备盲清（问题四）。\n" + "所有频道销账 ⇒ 完成证明成立，立即 /exit——\n" + "虚拟时间就是成绩，不做无意义的复扫\n" + "（--persistent 可显式开启"扫到限时"保险模式）。", + fill=PALETTE["success_soft"], stroke=PALETTE["success"], title_size=11.5, + body_size=9.5) + footer_note(fig, "七个站点的巡回长度 6a ≈ 6.96 km；a 的可行下界 " + "a_min ≈ 1122.96 m（零余量），默认留 15 m 以上余量。") + return fig + + +# =========================================================================== +# 图 9：问题四最小命中集与联合信念 +# =========================================================================== +def _load_plan(name: str) -> dict[str, Any]: + path = ROOT / "src" / "mcm26b" / "guarantee" / "plans" / f"{name}.json" + return json.loads(path.read_text(encoding="utf-8")) + + +def fig_q4_cover() -> Figure: + plan = _load_plan("q4_mincover_v2") + pts = [(float(x), float(y)) for x, y in plan["points"]] + cert = plan["certificate"] + fig = Figure(1120, 600, title="问题四：定向源的 31 点最小命中集与联合信念", + caption="定向源只有 ±90° 前向半球能被发现：覆盖语义从"圆盘"改成"半盘 + 朝向"") + fig.source("src/mcm26b/guarantee/plans/q4_mincover_v2.json", + plan=plan["name"], certificate=cert["uncovered"]) + fig.source("src/mcm26b/planner/directional.py") + half = 2300.0 + gx, gy, gw, gh = 60, 70, 620, 440 + sx, sy = _world(fig, gx, gy, gw, gh, half, grid=600.0) + fig.circle(sx(0), sy(0), OMEGA_R * gw / (2 * half), fill="#f6f9fc", + stroke=PALETTE["primary"], sw=1.4) + fig.text(sx(0), sy(OMEGA_R) - 14, "Ω：R = 1800 m", size=10, + color=PALETTE["primary"], anchor="middle") + for i, (px, py) in enumerate(pts): + fig.circle(sx(px), sy(py), 4.0, fill=PALETTE["accent"], stroke=None) + if i < 6: + fig.text(sx(px) + 6, sy(py) - 6, f"P{i + 1}", size=8.5, + color=PALETTE["muted"], font="times") + fig.text(74, 500, f"31 个覆盖点（{plan['name']}）：离线贪心最小命中集，" + f"设计半径 {Q4_DESIGN_R:.0f} m + 前向角裕量 {Q4_HEADING_MARGIN_DEG:.0f}°", + size=10.5, color=PALETTE["ink"]) + fig.text(74, 522, f"运行期证书：{cert['targets']:,} 个 (位置, 朝向) 采样格全部被命中，" + f"缺口 {cert['uncovered']}，最坏命中距离 {cert['worst_hit_distance_m']:.1f} m", + size=10.5, color=PALETTE["success"]) + fig.text(74, 544, f"NN 巡回 {cert['nn_tour_m'] / 1000:.1f} km " + f"≈ {cert['nn_tour_virtual_s']:.0f} s 虚拟时间；" + "对照旧 q4_mincover_v1：粗采样过拟合，15 m/72 档下有 566 个缺口", + size=10.5, color=PALETTE["muted"]) + + rx = 720 + panel_title(fig, rx, 44, "(a) 覆盖语义：forward_half_disk") + _half_disk(fig, rx + 10, 76, 190, 150, PALETTE["accent"]) + box(fig, rx + 210, 76, 170, 150, "与问题三的差别", + "问题三 full_disk：\n任何位置负观测都\n能排除 1000 m 圆盘。\n\n" + "问题四：负观测只能\n排除"背对它"的朝向，\n位置信息几乎不更新。", + fill=PALETTE["warn_soft"], stroke=PALETTE["warn"], title_size=11, + body_size=9.5) + panel_title(fig, rx, 250, "(b) 位置 × 朝向联合信念") + box(fig, rx, 278, 380, 128, "JointBelief：网格 × 朝向档", + "每个覆盖点观测后只更新"哪些朝向还活着"：\n" + "· direction 命中 ⇒ 保留该朝向档，位置信息由后续交会补；\n" + "· no_signal ⇒ 删去"朝向正对本站点"的档（含格内角裕量）；\n" + "· 可见性必须与引擎一致：距离 ≤ r_eff 且夹角 ≤ 90°；\n" + "· 后备盲清：联合信念坍缩到有限格时，按可疑格就近清。", + fill=PALETTE["primary_soft"], stroke=PALETTE["primary"], title_size=11.5, + body_size=9.5) + box(fig, rx, 420, 380, 118, "本轮修复的正确性缺陷", + "旧版 no_signal 朝向裁剪规则反向、且漏计格内角裕量，\n" + "会删掉真源所在朝向档。现已修复，并由\n" + "tests/test_joint_belief.py 的随机观测序列保留性测试固化。", + fill=PALETTE["accent_soft"], stroke=PALETTE["accent"], title_size=11.5, + body_size=9.5) + footer_note(fig, "后备盲清 125.7 → 4.3 次/局，Q4 最坏虚拟时间 10891 → 9630 s；" + "负面结论与开关见 docs/第三四问图论与动态规划优化方案.md §8。") + return fig + + +def _half_disk(fig: Figure, x: float, y: float, w: float, h: float, color: str) -> None: + """画一个 forward_half_disk：圆盘 + 前向半平面高亮 + 朝向箭头。""" + cx, cy, r = x + w / 2, y + h / 2, min(w, h) / 2 - 8 + fig.circle(cx, cy, r, fill="#f7f9fb", stroke=PALETTE["muted"], sw=1.0) + pts = [[cx, cy]] + for i in range(25): + a = math.radians(-90 + 180 * i / 24) + pts.append([cx + r * math.cos(a), cy + r * math.sin(a)]) + fig.polygon(pts, fill=color, stroke=None, opacity=0.35) + fig.arrow(cx, cy, cx + r * 0.86, cy, stroke=color, sw=1.4, head=7) + fig.line(cx, cy - r, cx, cy + r, stroke=color, sw=1.0, dash="4 4") + fig.text(cx, cy + r + 14, "朝向未知 ⇒ ±90°", size=9.5, color=color, + anchor="middle") + fig.text(cx, cy - r - 8, "覆盖半径 1000 m", size=9.5, color=PALETTE["muted"], + anchor="middle") + + +# =========================================================================== +# 图 10/11：WebUI 回放预览复现 +# =========================================================================== +def replay_figure(run: LogRun, *, width: float = 1200, height: float = 760, + title: str | None = None) -> Figure: + """复现 referee WebUI"回放"Tab 的画布语义：世界 Ω、真值源、运动轨迹、 + 检测/清除事件与虚拟时刻进度条。""" + tl = decompose(run) + fig = Figure(width, height, title=title or f"WebUI 回放复现：{run.case_code}", + caption="与 referee 静态回放器同一语义：移动 5 m/s 直线展开、" + "检测附实测示向度射线、清除先光学定位再清除") + fig.source(run.source, case_code=run.case_code, module=run.module, + truth=len(run.truth), actions=len(run.actions)) + + map_x, map_y, map_w = 40.0, 64.0, 700.0 + map_h = height - map_y - 132 + half = 2100.0 + sx, sy = _world(fig, map_x, map_y, map_w, map_h, half, grid=1000.0) + fig.circle(sx(0), sy(0), OMEGA_R * map_w / (2 * half), fill="#f7fafc", + stroke=PALETTE["primary"], sw=1.6) + fig.text(sx(0), sy(OMEGA_R) - 16, "Ω（半径 1800 m）", size=10.5, + color=PALETTE["primary"], anchor="middle") + + # 真值源 + solved = set(run.solution_channels()) + for src in run.truth: + ch = int(src["channel"]) + directional = bool(src.get("directional")) + color = PALETTE["accent"] if ch in solved else PALETTE["warn"] + px, py = sx(float(src["x"])), sy(float(src["y"])) + fig.circle(px, py, 6.0, fill=color, stroke=PALETTE["paper"], sw=1.2) + fig.text(px + 9, py - 4, f"ch{ch}", size=9.5, font="times_bold", color=color) + if directional: + heading = math.radians(float(src.get("heading_deg") or 0.0)) + fig.arrow(px, py, px + 26 * math.cos(heading), py - 26 * math.sin(heading), + stroke=PALETTE["accent"], sw=1.2, head=6) + + # 运动轨迹（按动作序列连接 position） + path: list[list[float]] = [[sx(0), sy(0)]] + pos = (0.0, 0.0) + for act in run.actions: + req = act["request"] + if act["path"] in ("/measure", "/clear") and "position" in req: + tgt = (float(req["position"]["x"]), float(req["position"]["y"])) + if tgt != pos: + path.append([sx(*tgt)]) + pos = tgt + if len(path) > 1: + fig.polyline(path, fill=None, stroke=PALETTE["primary"], sw=1.4, opacity=0.85) + fig.circle(sx(0), sy(0), 5.0, fill=PALETTE["primary"], stroke=PALETTE["paper"], sw=1.2) + fig.text(sx(0) + 8, sy(0) + 12, "起点 (0,0)", size=9.5, color=PALETTE["primary"]) + + # 检测射线（direction）与清除点 + for act, row in zip(run.actions, tl.rows): + if act["path"] == "/measure" and str(act["response"].get("measure_result")) == "direction": + bearing = float(act["response"].get("svd_deg") or 0.0) + if row["x"] is None: + continue + px, py = sx(float(row["x"])), sy(float(row["y"])) + ex = px + 150 * math.cos(math.radians(bearing)) + ey = py - 150 * math.sin(math.radians(bearing)) + fig.line(px, py, ex, ey, stroke=PALETTE["success"], sw=0.9, opacity=0.5) + fig.circle(px, py, 2.6, fill=PALETTE["success"], stroke=None) + elif act["path"] == "/clear": + if row["x"] is None: + continue + px, py = sx(float(row["x"])), sy(float(row["y"])) + ok = act["response"].get("clear_result") == "cleared" + fig.circle(px, py, 8.0, fill=None, + stroke=PALETTE["success"] if ok else PALETTE["muted"], sw=1.6) + if ok: + fig.line(px - 5, py - 5, px + 5, py + 5, stroke=PALETTE["success"], sw=1.8) + fig.line(px - 5, py + 5, px + 5, py - 5, stroke=PALETTE["success"], sw=1.8) + + # 右栏：状态与事件 + rx = map_x + map_w + 28 + rw = width - rx - 32 + panel_title(fig, rx, 44, "回放状态") + chips = [ + (f"案例 {run.case_code}", PALETTE["primary"]), + (f"{run.module}", PALETTE["muted"]), + (f"参赛队号 {run.robot_id}", PALETTE["muted"]), + (f"真值源 {len(run.truth)}（定向 {sum(1 for s in run.truth if s.get('directional'))}）", + PALETTE["accent"]), + (f"全清 {len(solved)}/{len(run.truth)}", PALETTE["success"]), + (f"动作 {len(run.actions)} · 检测 {run.measures} · 清除 {run.clears}", + PALETTE["primary"]), + ] + y = 74 + for label, color in chips: + fig.rect(rx, y, rw, 24, fill=PALETTE["paper"], stroke=color, sw=1.1, radius=12) + fig.text(rx + 12, y + 12, label, size=10, color=color, align="middle") + y += 30 + + panel_title(fig, rx, y + 12, "事件流节选（前 12 条动作）") + y += 40 + fig.rect(rx, y, rw, 22, fill=PALETTE["primary_soft"], stroke=None) + for x, head in ((rx + 8, "#"), (rx + 40, "指令"), (rx + 150, "参数"), (rx + 268, "虚拟 t/s")): + fig.text(x, y + 11, head, size=9.5, font="hei_bold", align="middle") + y += 22 + for i, (act, row) in enumerate(zip(run.actions, tl.rows)): + if i >= 12: + break + fill = "#f7f9fb" if i % 2 == 0 else None + if fill: + fig.rect(rx, y, rw, 20, fill=fill, stroke=None) + detail = "" + if act["path"] == "/measure": + detail = f"ch{row['channel']} ({row['x']:.0f},{row['y']:.0f})" + elif act["path"] == "/clear": + detail = f"ch{row['channel']} {row['detail']}" + elif act["path"] == "/enter": + detail = "启动" + elif act["path"] == "/exit": + detail = "退出" + fig.text(rx + 8, y + 10, str(i + 1), size=9, font="times", align="middle") + fig.text(rx + 40, y + 10, act["path"], size=9, font="times_bold", align="middle") + fig.text(rx + 150, y + 10, detail, size=9, font="times", align="middle") + fig.text(rx + 268, y + 10, f"{row['virtual_s']:.0f}", size=9, font="times", + align="middle") + y += 20 + fig.text(rx, y + 16, f"… 共 {len(run.actions)} 条动作；" + f"终点虚拟时刻 {tl.total_virtual_s:,.1f} s", + size=9.5, color=PALETTE["muted"]) + + # 底部：进度条（对齐回放器的可拖动时间轴） + bar_y = height - 92 + fig.rect(map_x, bar_y, map_w, 12, fill=PALETTE["faint"], stroke=None, radius=6) + fig.rect(map_x, bar_y, map_w, 12, fill=PALETTE["primary"], stroke=None, radius=6) + for act, row in zip(run.actions, tl.rows): + frac = row["virtual_s"] / max(tl.total_virtual_s, 1e-9) + px = map_x + map_w * frac + if act["path"] == "/clear": + fig.line(px, bar_y - 5, px, bar_y + 17, stroke=PALETTE["success"], sw=1.6) + elif act["path"] == "/measure" and str(act["response"].get("measure_result")) == "direction": + fig.line(px, bar_y - 3, px, bar_y + 15, stroke=PALETTE["accent"], sw=1.2) + fig.text(map_x, bar_y + 30, "0 s", size=9.5, color=PALETTE["muted"], font="times") + fig.text(map_x + map_w / 2, bar_y + 30, "时间轴刻度：竖线 = 检测到信号（红）/ 清除（绿）", + size=9.5, color=PALETTE["muted"], anchor="middle") + fig.text(map_x + map_w, bar_y + 30, f"{tl.total_virtual_s:,.0f} s", size=9.5, + color=PALETTE["muted"], anchor="end", font="times") + fig.legend_chips(map_x, bar_y + 54, + [("已清除频道", PALETTE["success"]), ("未清除频道（真值）", PALETTE["warn"]), + ("定向源朝向", PALETTE["accent"]), ("示向度射线", PALETTE["success"])], + size=9.5) + return fig + + +# =========================================================================== +# 入口 +# =========================================================================== +def build(only: Sequence[str] | None, scale: float, rasterizer: str) -> list[dict[str, Any]]: + want = set(only) if only else None + runs = discover_runs() + q3 = pick_run(runs, module="q3", truth=10) + q4 = pick_run(runs, module="q4", truth=12) + + jobs: list[tuple[str, Any]] = [ + ("flow", lambda: fig_flow()), + ("steps", lambda: fig_steps()), + ("timing-rule", lambda: fig_timing_rule()), + ("geo-q1", lambda: fig_q1_meetup()), + ("geo-q1-counterexample", lambda: fig_q1_counterexample()), + ("geo-q2", lambda: fig_q2_region()), + ("geo-q3", lambda: fig_q3_coverage()), + ("geo-q4", lambda: fig_q4_cover()), + ] + if q3: + jobs.append(("timing-run", lambda q3=q3: fig_timing_run(q3))) + jobs.append(("replay-q3", lambda q3=q3: replay_figure(q3))) + if q4: + jobs.append(("replay-q4", lambda q4=q4: replay_figure( + q4, title=f"WebUI 回放复现（定向）：{q4.case_code}"))) + + out: list[dict[str, Any]] = [] + for stem, make in jobs: + if want and stem not in want: + continue + fig = make() + meta = fig.write(OUT_DIR, stem, scale=scale, rasterizer=rasterizer) + out.append(meta) + print(f"[fig] {stem:<24} {meta['size'][0]:.0f}×{meta['size'][1]:.0f} " + f"→ {meta.get('png', meta['svg'])}") + return out + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--only", default="", help="逗号分隔的图名子集（见 --list）") + parser.add_argument("--scale", type=float, default=2.0, + help="PNG 相对 SVG 的像素倍率（论文插图建议 2.0~3.0）") + parser.add_argument("--rasterizer", default=DEFAULT_RASTERIZER) + parser.add_argument("--list", action="store_true", help="列出全部图名后退出") + args = parser.parse_args() + names = ["flow", "steps", "timing-rule", "timing-run", "geo-q1", + "geo-q1-counterexample", "geo-q2", "geo-q3", "geo-q4", + "replay-q3", "replay-q4"] + if args.list: + print("\n".join(names)) + return 0 + only = [s.strip() for s in args.only.split(",") if s.strip()] + entries = build(only or None, args.scale, args.rasterizer) + command = ("uv run python scripts/build_docs_figures.py" + + (f" --only {args.only}" if args.only else "")) + path = write_manifest(OUT_DIR, entries, command) + print(f"[fig] manifest → {path.relative_to(ROOT)}（{len(entries)} 张图）") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/figkit.py b/scripts/figkit.py new file mode 100644 index 0000000..40eed26 --- /dev/null +++ b/scripts/figkit.py @@ -0,0 +1,372 @@ +"""论文级图表工具：统一配色、场景图 DSL、SVG/PNG 双出与来源清单。 + +设计约束（对应 docs/figures/README.md 的"可复现"承诺）： + +- **零新增运行时依赖**：几何/连线由本模块直接写 SVG，PNG 由 ``rsvg-convert`` + 光栅化（与项目既有可视化链路一致）；中文与西文文字层用 Pillow 以论文自带 + 的 ``docs/paper/fonts/`` 字体单独合成，保证与 ``main.tex`` 同字体、同字号； +- **一份场景图两处消费**：SVG 用于 WebUI/文档内联预览（走系统字体栈）， + PNG 用于论文插图（走 bundled 字体）；两者坐标完全一致； +- **来源可追溯**：每张图强制携带 ``sources``（文件路径 + sha256 前 12 位 + + 行数）与生成命令，写入 ``manifest.json``，避免"图与数据对不上"。 + +文字排版约定：**主风格只提供 SVG 精度的文字**，Pillow 侧复用同一套 +``(x, y, anchor, font, size, dy)`` 参数重画一遍；``dy`` 表示"文字视觉中心相对 +基线的偏移比例"，由 font-size 缩放，用来让框内文字纵向居中。 +""" + +from __future__ import annotations + +import hashlib +import json +import math +import shutil +import subprocess +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any, Iterable, Sequence + +ROOT = Path(__file__).resolve().parents[1] +FONT_DIR = ROOT / "docs" / "paper" / "fonts" +DEFAULT_RASTERIZER = "rsvg-convert" + +# --------------------------------------------------------------------------- +# 论文配色：低饱和、单色可打印，正文/强调/真值/成功/失败五类语义色 + 顺序色阶 +# --------------------------------------------------------------------------- +PALETTE: dict[str, str] = { + "ink": "#1c2430", # 主文字/坐标轴 + "muted": "#6b7684", # 次要文字与网格 + "faint": "#dfe4ea", # 分隔线 + "paper": "#ffffff", + "primary": "#1f4e79", # 主线：协议/流程 + "primary_soft": "#dce6f1", + "accent": "#b03a2e", # 强调：真值源/关键结论 + "accent_soft": "#f6dcd8", + "success": "#1e7a5f", # 成功：清除/全清 + "success_soft": "#d8ece5", + "warn": "#b8860b", # 警示：后备/未达标 + "warn_soft": "#f6ecd2", + "series": ("#1f4e79", "#b03a2e", "#1e7a5f", "#b8860b", "#5b4b8a", "#3d7ea6"), +} + +FONTS: dict[str, str] = { + "song": "SongtiSC-Regular.ttf", # 中文正文：宋体 + "song_bold": "SongtiSC-Bold.ttf", + "hei": "HeitiSC-Light.ttf", # 中文标题：黑体 + "hei_bold": "HeitiSC-Medium.ttf", + "times": "Times New Roman.ttf", # 西文与数字 + "times_bold": "Times New Roman Bold.ttf", +} + +_SVG_FAMILY = { + "song": "'Songti SC','SimSun',serif", + "song_bold": "'Songti SC','SimSun',serif", + "hei": "'Heiti SC','Microsoft YaHei',sans-serif", + "hei_bold": "'Heiti SC','Microsoft YaHei',sans-serif", + "times": "'Times New Roman',Times,serif", + "times_bold": "'Times New Roman',Times,serif", +} + + +# --------------------------------------------------------------------------- +# 场景图 DSL +# --------------------------------------------------------------------------- +@dataclass +class Figure: + """一张图的场景描述：宽高以 CSS 像素（≈1/96 inch）计。""" + + width: float + height: float + title: str = "" + caption: str = "" + background: str = PALETTE["paper"] + elements: list[dict[str, Any]] = field(default_factory=list) + sources: list[dict[str, Any]] = field(default_factory=list) + notes: list[str] = field(default_factory=list) + + # ---------------------------------------------------------------- 图元 + def rect(self, x: float, y: float, w: float, h: float, **kw: Any) -> "Figure": + self.elements.append(dict(kind="rect", x=x, y=y, w=w, h=h, **kw)) + return self + + def circle(self, cx: float, cy: float, r: float, **kw: Any) -> "Figure": + self.elements.append(dict(kind="circle", cx=cx, cy=cy, r=r, **kw)) + return self + + def line(self, x1: float, y1: float, x2: float, y2: float, **kw: Any) -> "Figure": + self.elements.append(dict(kind="line", x1=x1, y1=y1, x2=x2, y2=y2, **kw)) + return self + + def polyline(self, points: Sequence[Sequence[float]], **kw: Any) -> "Figure": + self.elements.append(dict(kind="polyline", points=[list(p) for p in points], **kw)) + return self + + def polygon(self, points: Sequence[Sequence[float]], **kw: Any) -> "Figure": + self.elements.append(dict(kind="polygon", points=[list(p) for p in points], **kw)) + return self + + def path(self, d: str, **kw: Any) -> "Figure": + self.elements.append(dict(kind="path", d=d, **kw)) + return self + + def arrow(self, x1: float, y1: float, x2: float, y2: float, **kw: Any) -> "Figure": + """直线箭头：``head`` 为箭头边长（像素），``dash`` 支持虚线。""" + self.elements.append(dict(kind="arrow", x1=x1, y1=y1, x2=x2, y2=y2, **kw)) + return self + + def text(self, x: float, y: float, s: str, *, size: float = 12.0, + font: str = "song", color: str | None = None, anchor: str = "start", + align: str = "top", bold: bool = False, rotate: float | None = None, + dy: float | None = None, opacity: float = 1.0) -> "Figure": + """文字。``(x, y)`` 为参考点；``align`` ∈ {top, middle, bottom} 决定 + 参考点落在文字块的哪条水平线上（``anchor`` 决定左右）。 + """ + self.elements.append(dict( + kind="text", x=x, y=y, s=s, size=size, font=font, + color=color or PALETTE["ink"], anchor=anchor, align=align, + bold=bold, rotate=rotate, dy=dy, opacity=opacity)) + return self + + def legend_chips(self, x: float, y: float, items: Iterable[tuple[str, str]], + *, size: float = 11.0, gap: float = 16.0, box: float = 10.0, + font: str = "song") -> float: + """一行"色块 + 文字"图例，返回结束时的 x 坐标。""" + cur = x + for label, color in items: + self.rect(cur, y - box * 0.75, box, box, fill=color, stroke=None, radius=2.0) + self.text(cur + box + 5.0, y + box * 0.35, label, size=size, font=font, + align="middle") + cur += box + 5.0 + _estimate_width(label, size, font) + gap + return cur + + def source(self, path: str | Path, **extra: Any) -> "Figure": + """登记本图的数据来源（相对仓库根），含 sha256 前 12 位与行数。""" + p = Path(path) + p = p if p.is_absolute() else ROOT / p + digest, lines = "-", 0 + if p.is_file(): + data = p.read_bytes() + digest = hashlib.sha256(data).hexdigest()[:12] + lines = data.count(b"\n") + (0 if data.endswith(b"\n") else 1) + self.sources.append(dict(path=str(p.relative_to(ROOT)) if p.is_relative_to(ROOT) else str(p), + sha256_12=digest, lines=lines, **extra)) + return self + + # ---------------------------------------------------------------- 输出 + def to_svg(self) -> str: + body: list[str] = [_svg_element(el) for el in self.elements] + head = (f'') + meta = f"{_esc(self.title)}\n" + if self.caption: + meta += f"{_esc(self.caption)}\n" + return "\n".join([head, meta, f'', *body, "", ""]) + + def write(self, out_dir: str | Path, stem: str, *, scale: float = 2.0, + rasterizer: str = DEFAULT_RASTERIZER, png: bool = True) -> dict[str, Any]: + """写出 ``.svg``（+ ``.png``），返回产出清单。""" + out = Path(out_dir) + out.mkdir(parents=True, exist_ok=True) + svg_path = out / f"{stem}.svg" + svg_path.write_text(self.to_svg(), encoding="utf-8") + produced: dict[str, Any] = {"stem": stem, "title": self.title, + "caption": self.caption, + "svg": str(svg_path.relative_to(ROOT)), + "size": [self.width, self.height], + "sources": self.sources, "notes": self.notes} + if png: + png_path = out / f"{stem}.png" + _rasterize(svg_path, png_path, scale=scale, rasterizer=rasterizer) + _compose_text(self, png_path, scale=scale) + produced["png"] = str(png_path.relative_to(ROOT)) + produced["png_scale"] = scale + return produced + + +# --------------------------------------------------------------------------- +# SVG 生成 +# --------------------------------------------------------------------------- +def _n(value: float) -> str: + """紧凑数字：去掉多余小数位，避免 SVG 体积膨胀。""" + if abs(value - round(value)) < 1e-9: + return str(int(round(value))) + return f"{value:.3f}".rstrip("0").rstrip(".") + + +def _esc(text: str) -> str: + return (str(text).replace("&", "&").replace("<", "<").replace(">", ">")) + + +def _paint(kw: dict[str, Any]) -> str: + parts: list[str] = [] + fill = kw.get("fill", "none") + if fill is None: + fill = "none" + parts.append(f'fill="{fill}"') + stroke = kw.get("stroke") + if stroke: + parts.append(f'stroke="{stroke}" stroke-width="{_n(kw.get("sw", 1.0))}"') + if kw.get("dash"): + dash = kw["dash"] if isinstance(kw["dash"], str) else \ + " ".join(_n(v) for v in kw["dash"]) + parts.append(f'stroke-dasharray="{dash}"') + parts.append(f'stroke-linejoin="{kw.get("linejoin", "miter")}"') + if kw.get("opacity") is not None: + parts.append(f'opacity="{_n(kw["opacity"])}"') + if kw.get("radius"): + parts.append(f'rx="{_n(kw["radius"])}" ry="{_n(kw["radius"])}"') + return " ".join(parts) + + +def _text_svg(el: dict[str, Any]) -> str: + size = el["size"] + dy = _baseline_offset(el) + anchor = {"start": "start", "middle": "middle", "end": "end"}[el["anchor"]] + family = _SVG_FAMILY[el["font"]] + weight = ' font-weight="bold"' if el["bold"] or el["font"].endswith("_bold") else "" + style = ' font-style="italic"' if el.get("italic") else "" + extra = f' opacity="{_n(el["opacity"])}"' if el.get("opacity", 1.0) < 1 else "" + rotate = f' transform="rotate({_n(el["rotate"])} {_n(el["x"])} {_n(el["y"])})"' \ + if el.get("rotate") else "" + return (f'{_esc(el["s"])}') + + +def _baseline_offset(el: dict[str, Any]) -> float: + """参考点 → SVG 基线（alphabetic）的 y 偏移。""" + if el.get("dy") is not None: + return float(el["dy"]) + if el["align"] == "middle": + return el["size"] * 0.35 + if el["align"] == "bottom": + return 0.0 + return el["size"] * 0.82 + + +def _svg_element(el: dict[str, Any]) -> str: + kind = el["kind"] + if kind == "text": + return _text_svg(el) + if kind == "rect": + return (f'') + if kind == "circle": + return f'' + if kind == "line": + return (f'') + if kind in ("polyline", "polygon"): + pts = " ".join(f"{_n(px)},{_n(py)}" for px, py in el["points"]) + return f'<{kind} points="{pts}" {_paint(el)}/>' + if kind == "path": + return f'' + if kind == "arrow": + return _arrow_svg(el) + raise ValueError(f"未知图元：{kind}") + + +def _arrow_svg(el: dict[str, Any]) -> str: + x1, y1, x2, y2 = el["x1"], el["y1"], el["x2"], el["y2"] + color = el.get("stroke") or PALETTE["ink"] + sw = el.get("sw", 1.4) + head = el.get("head", 8.0) + ang = math.atan2(y2 - y1, x2 - x1) + bx, by = x2 - head * math.cos(ang), y2 - head * math.sin(ang) + left = (bx - head * 0.5 * math.sin(ang), by + head * 0.5 * math.cos(ang)) + right = (bx + head * 0.5 * math.sin(ang), by - head * 0.5 * math.cos(ang)) + dash = f' stroke-dasharray="{el["dash"]}"' if el.get("dash") else "" + return "\n".join([ + f'', + f'', + ]) + + +# --------------------------------------------------------------------------- +# PNG：rsvg-convert 出几何层 + Pillow 出文字层（论文同款字体） +# --------------------------------------------------------------------------- +def _rasterize(svg_path: Path, png_path: Path, *, scale: float, + rasterizer: str) -> None: + exe = shutil.which(rasterizer) + if not exe: + raise RuntimeError(f"未找到 {rasterizer}：PNG 输出需要它（brew install librsvg）") + cmd = [exe, "-f", "png", "-o", str(png_path), "--zoom", str(scale), str(svg_path)] + proc = subprocess.run(cmd, capture_output=True, text=True, check=False) + if proc.returncode != 0: + raise RuntimeError(f"{rasterizer} 失败：{proc.stderr.strip()}") + + +def _compose_text(fig: Figure, png_path: Path, *, scale: float) -> None: + """把文字层用 Pillow 重画在光栅图上：与基础几何同色同坐标，字体来自 + ``docs/paper/fonts/``（缺失时回退系统字体）。""" + from PIL import Image, ImageDraw, ImageFont # 可选依赖：仅生成图时导入 + + texts = [el for el in fig.elements if el["kind"] == "text"] + if not texts: + return + img = Image.open(png_path).convert("RGBA") + layer = Image.new("RGBA", img.size, (0, 0, 0, 0)) + draw = ImageDraw.Draw(layer) + cache: dict[tuple[str, float], Any] = {} + for el in texts: + key = (el["font"], round(el["size"] * scale, 2)) + if key not in cache: + cache[key] = _load_font(el["font"], el["size"] * scale, ImageFont) + font = cache[key] + x, y = el["x"] * scale, (el["y"] + _baseline_offset(el)) * scale + anchor = {"start": "ls", "middle": "ms", "end": "rs"}[el["anchor"]] + draw.text((x, y), el["s"], font=font, fill=_rgba(el["color"], el.get("opacity", 1.0)), + anchor=anchor) + Image.alpha_composite(img, layer).convert("RGB").save(png_path, optimize=True) + + +def _load_font(name: str, size: float, image_font: Any) -> Any: + path = FONT_DIR / FONTS[name] + if path.is_file(): + return image_font.truetype(str(path), max(6, int(round(size)))) + fallback = { + "song": "/System/Library/Fonts/Supplemental/Songti.ttc", + "song_bold": "/System/Library/Fonts/Supplemental/Songti.ttc", + "hei": "/System/Library/Fonts/Supplemental/Songti.ttc", + "hei_bold": "/System/Library/Fonts/Supplemental/Songti.ttc", + "times": "/System/Library/Fonts/Supplemental/Times New Roman.ttf", + "times_bold": "/System/Library/Fonts/Supplemental/Times New Roman Bold.ttf", + }[name] + return image_font.truetype(fallback, max(6, int(round(size)))) + + +def _rgba(color: str, opacity: float) -> tuple[int, int, int, int]: + value = color.lstrip("#") + return (int(value[0:2], 16), int(value[2:4], 16), int(value[4:6], 16), + max(0, min(255, int(round(opacity * 255))))) + + +def _estimate_width(text: str, size: float, font: str) -> float: + """粗略宽度估计（仅用于图例排版留白）：CJK 按 1.0 em，西文按 0.5 em。""" + wide = sum(1 for ch in text if ord(ch) > 0x2E80) + narrow = len(text) - wide + factor = 1.0 if font.startswith(("song", "hei")) else 0.98 + return (wide * size + narrow * size * 0.52) * factor + + +def write_manifest(out_dir: str | Path, entries: Sequence[dict[str, Any]], + command: str) -> Path: + """落盘 manifest.json：图清单 + 数据来源 + 复现命令。""" + out = Path(out_dir) + out.mkdir(parents=True, exist_ok=True) + doc = { + "generator": "scripts/build_docs_figures.py", + "command": command, + "palette": {k: v for k, v in PALETTE.items() if k != "series"}, + "fonts": {"dir": "docs/paper/fonts", "map": FONTS}, + "figures": list(entries), + } + path = out / "manifest.json" + path.write_text(json.dumps(doc, ensure_ascii=False, indent=2) + "\n", encoding="utf-8") + return path