748 lines
34 KiB
Python
748 lines
34 KiB
Python
"""为 MATLAB 出图准备数据，并更新 docs/figures/manifest.json。
|
||
|
||
uv run python scripts/build_figure_data.py --prepare # 只写 .mat 数据
|
||
uv run python scripts/build_figure_data.py # 写数据 + 调 MATLAB 出图
|
||
uv run python scripts/build_figure_data.py --only geo-q3-seven-point-coverage
|
||
uv run python scripts/build_figure_data.py --list
|
||
|
||
分工：
|
||
|
||
- **本脚本（Python）**：解析行为日志与覆盖计划、调用生产求解器算几何、把绘图数据写成
|
||
``docs/figures/.data/*.mat``、维护 manifest（来源路径 + sha256 + 行数）；
|
||
- **MATLAB**：``scripts/matlab/*.m`` 负责排版与出图，只读 .mat，不做数据解析。
|
||
手动出图：MATLAB 中 ``cd scripts/matlab`` 后执行 ``run_all_figures``。
|
||
|
||
图内不放标题、图注与解释性文字，只保留坐标轴、图例与数据标注；
|
||
每张图的含义由文件名与 ``docs/figures/README.md`` 承担。
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import argparse
|
||
import hashlib
|
||
import json
|
||
import math
|
||
import subprocess
|
||
import sys
|
||
from pathlib import Path
|
||
from typing import Any, Callable, Sequence
|
||
|
||
import numpy as np
|
||
from scipy.io import savemat
|
||
|
||
ROOT = Path(__file__).resolve().parents[1]
|
||
OUT_DIR = ROOT / "docs" / "figures"
|
||
DATA_DIR = OUT_DIR / ".data"
|
||
MATLAB_DIR = ROOT / "scripts" / "matlab"
|
||
SESSIONS_DIR = ROOT / "referee-data" / "sessions"
|
||
ARCHIVE_DIR = ROOT / "examples" / "replay"
|
||
MATLAB_BIN = Path("/Applications/MATLAB_R2025a.app/bin/matlab")
|
||
|
||
# 题面常量（与 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
|
||
|
||
# 回放图指定的 q4 演练局（缺失时回退到自动挑选）
|
||
Q4_REPLAY_CASE = "GHSN-M38E-74GD-YLUD"
|
||
SEVEN_R = 1800.0
|
||
|
||
# 图名 → (MATLAB 函数, .mat 数据名)
|
||
FIGURES: tuple[tuple[str, str, str], ...] = (
|
||
("flow-overall-architecture", "fig_flow_overall", "flow_overall"),
|
||
("steps-state-machine", "fig_steps_state_machine", "steps_state_machine"),
|
||
("timing-attachment-example-105-111-194-199", "fig_timing_attachment",
|
||
"timing_attachment"),
|
||
("timing-measured-q3-10-sources", "fig_timing_measured", "timing_measured"),
|
||
("geo-q1-convex-intersection-zoom", "fig_geo_q1_meetup", "geo_q1_meetup"),
|
||
("geo-q1-diameter-circle-counterexample", "fig_geo_q1_counterexample",
|
||
"geo_q1_counterexample"),
|
||
("geo-q2-guaranteed-reception-disk", "fig_geo_q2_safe_disk", "geo_q2"),
|
||
("geo-q3-seven-point-coverage", "fig_geo_q3_coverage", "geo_q3"),
|
||
("geo-q4-minimum-cover-31-points", "fig_geo_q4_mincover", "geo_q4"),
|
||
("replay-webui-q3-omnidirectional", "fig_replay", "replay_q3"),
|
||
("replay-webui-q4-directional", "fig_replay", "replay_q4"),
|
||
)
|
||
|
||
SOURCES: dict[str, list[str]] = {
|
||
"flow-overall-architecture": ["src/mcm26b/config.py", "docs/B题建模分析与实施计划.md"],
|
||
"steps-state-machine": ["docs/协议规则表与状态机.md", "src/mcm26b/protocol/client.py"],
|
||
"timing-attachment-example-105-111-194-199": ["docs/协议规则表与状态机.md"],
|
||
"geo-q1-diameter-circle-counterexample": ["examples/q1/triangle.json"],
|
||
"geo-q2-guaranteed-reception-disk": ["src/mcm26b/config.py", "docs/algo-2.md"],
|
||
"geo-q3-seven-point-coverage": ["src/mcm26b/planner/seven.py", "docs/algo-3.md"],
|
||
"geo-q4-minimum-cover-31-points": [
|
||
"src/mcm26b/guarantee/plans/q4_mincover_v2.json", "docs/algo-4.md"],
|
||
}
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 数据来源
|
||
# ---------------------------------------------------------------------------
|
||
def source(path: str | Path, **extra: Any) -> dict[str, Any]:
|
||
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)
|
||
rel = str(p.relative_to(ROOT)) if p.is_relative_to(ROOT) else str(p)
|
||
return {"path": rel, "sha256_12": digest, "lines": lines, **extra}
|
||
|
||
|
||
def parse_jsonl(path: Path) -> dict[str, Any]:
|
||
"""兼容评测机 ``event=action`` 与客户端 ``event=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[str(ev.get("request_id", ""))] = ev
|
||
elif kind == "step" and str(ev.get("request_id", "")) not in fallback:
|
||
fallback[str(ev.get("request_id", ""))] = ev
|
||
if not actions:
|
||
for ev in fallback.values():
|
||
if not ev.get("path"):
|
||
continue
|
||
seq += 1
|
||
stamp = ev.get("logged_at_unix_s")
|
||
actions.append(dict(seq=seq, path=ev["path"], request=ev.get("request") or {},
|
||
response=ev.get("response") or {},
|
||
wall_ts_ms=stamp * 1000 if stamp else None))
|
||
if not actions:
|
||
raise ValueError(f"{path} 没有可回放的动作事件")
|
||
return dict(case_code=meta.get("case_code") or path.parent.name,
|
||
module=meta.get("module") or "?", robot_id=str(meta.get("robot_id") or "?"),
|
||
truth=list(meta.get("truth") or []), actions=actions, source=path)
|
||
|
||
|
||
def discover_runs() -> list[dict[str, Any]]:
|
||
runs: list[dict[str, Any]] = []
|
||
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_case(runs: Sequence[dict[str, Any]], case_code: str) -> dict[str, Any] | None:
|
||
"""按案例编码取局；不在本机留档里就返回 None，由调用方回退。"""
|
||
return next((r for r in runs if r["case_code"] == case_code), None)
|
||
|
||
|
||
def pick_run(runs: Sequence[dict[str, Any]], module: str, truth: int) -> dict[str, Any]:
|
||
cands = [r for r in runs if r["module"].startswith(module) and len(r["truth"]) >= truth
|
||
and any(a["path"] == "/measure" for a in r["actions"])]
|
||
if not cands:
|
||
cands = [r for r in runs if r["module"].startswith(module)]
|
||
return sorted(cands, key=lambda r: (abs(len(r["truth"]) - truth), len(r["actions"])))[0]
|
||
|
||
|
||
def decompose(run: dict[str, Any]) -> dict[str, Any]:
|
||
"""按附件规则逐动作复算移动 / 换频 / 动作耗时。"""
|
||
pos = (0.0, 0.0)
|
||
channel: int | None = None
|
||
rows: list[dict[str, Any]] = []
|
||
for act in run["actions"]:
|
||
req, res = act["request"], act["response"]
|
||
move_m = switch = action_s = 0.0
|
||
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":
|
||
action_s = CLEAR_OK_S if str(res.get("clear_result")) == "success" else CLEAR_NO_S
|
||
rows.append(dict(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),
|
||
hit=str(res.get("measure_result")) == "direction",
|
||
cleared=str(res.get("clear_result")) == "success",
|
||
channel=req.get("channel"),
|
||
x=req.get("position", {}).get("x"),
|
||
y=req.get("position", {}).get("y")))
|
||
last = run["actions"][-1]["response"]
|
||
return dict(rows=rows, virtual_s=float(last.get("virtual_time_s") or 0.0),
|
||
move_m=sum(r["move_m"] for r in rows),
|
||
measures=sum(1 for r in rows if r["path"] == "/measure"),
|
||
clears=sum(1 for r in rows if r["path"] == "/clear"))
|
||
|
||
|
||
def solved_channels(run: dict[str, Any]) -> list[int]:
|
||
return sorted({int(a["request"]["channel"]) for a in run["actions"]
|
||
if a["path"] == "/clear" and "channel" in a["request"]
|
||
and str(a["response"].get("clear_result")) == "success"})
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 几何：半平面裁剪与问题一求解
|
||
# ---------------------------------------------------------------------------
|
||
def clip_halfplane(poly: np.ndarray, a: float, b: float, c: float) -> np.ndarray:
|
||
"""Sutherland–Hodgman：半平面 a·x + b·y <= c 裁剪凸多边形。"""
|
||
if len(poly) == 0:
|
||
return poly
|
||
out: list[list[float]] = []
|
||
n = len(poly)
|
||
for i in range(n):
|
||
cur, nxt = poly[i], poly[(i + 1) % n]
|
||
dc = a * cur[0] + b * cur[1] - c
|
||
dn = a * nxt[0] + b * nxt[1] - c
|
||
if dc <= 1e-9:
|
||
out.append(list(cur))
|
||
if dn > 1e-9:
|
||
t = dc / (dc - dn)
|
||
out.append([cur[0] + t * (nxt[0] - cur[0]),
|
||
cur[1] + t * (nxt[1] - cur[1])])
|
||
elif dn <= 1e-9:
|
||
t = dc / (dc - dn)
|
||
out.append([cur[0] + t * (nxt[0] - cur[0]),
|
||
cur[1] + t * (nxt[1] - cur[1])])
|
||
return np.array(out) if out else np.empty((0, 2))
|
||
|
||
|
||
def angle_error(poly: np.ndarray, station: Sequence[float], bearing: float) -> float:
|
||
centre = poly.mean(axis=0)
|
||
dx, dy = centre[0] - station[0], centre[1] - station[1]
|
||
if abs(dx) < 1e-12 and abs(dy) < 1e-12:
|
||
return 0.0
|
||
return abs((math.degrees(math.atan2(dy, dx)) - bearing + 180.0) % 360.0 - 180.0)
|
||
|
||
|
||
def clip_wedge(region: np.ndarray, station: Sequence[float], bearing: float,
|
||
half_deg: float) -> np.ndarray:
|
||
"""用角域 W(s, θ ± δ) 裁剪区域：两个半平面各试一次，取方向一致的那一侧。"""
|
||
region = np.asarray(region, dtype=float)
|
||
if len(region) == 0:
|
||
return region
|
||
for edge in (bearing - half_deg, bearing + half_deg):
|
||
ang = math.radians(edge)
|
||
a, b = -math.sin(ang), math.cos(ang)
|
||
c = a * station[0] + b * station[1]
|
||
cands = [cand for cand in (clip_halfplane(region, a, b, c),
|
||
clip_halfplane(region, -a, -b, -c)) if len(cand) >= 3]
|
||
if not cands:
|
||
return np.empty((0, 2))
|
||
region = cands[0] if len(cands) == 1 else min(
|
||
cands, key=lambda cand: angle_error(cand, station, bearing))
|
||
return region
|
||
|
||
|
||
def wedge_polygon(station: Sequence[float], bearing: float, half_deg: float,
|
||
reach: float) -> np.ndarray:
|
||
a0, a1 = math.radians(bearing - half_deg), math.radians(bearing + half_deg)
|
||
ang = np.linspace(a0, a1, 24)
|
||
return np.vstack([[station[0], station[1]],
|
||
np.column_stack([station[0] + reach * np.cos(ang),
|
||
station[1] + reach * np.sin(ang)])])
|
||
|
||
|
||
def overview_window(stations: np.ndarray, verts: np.ndarray,
|
||
dc: dict[str, Any], mec: dict[str, Any]) -> tuple[float, np.ndarray]:
|
||
"""总览窗口：包含全部站点、定位区域与两个圆，中心取包围盒中心。"""
|
||
pts = [np.asarray(stations, dtype=float), np.asarray(verts, dtype=float)]
|
||
for circ in (dc, mec):
|
||
c = np.asarray(circ["center"], dtype=float)
|
||
r = float(circ["radius_m"])
|
||
pts.append(np.array([c + r, c - r]))
|
||
all_pts = np.vstack(pts)
|
||
lo, hi = all_pts.min(axis=0), all_pts.max(axis=0)
|
||
center = (lo + hi) / 2
|
||
half = float(np.max(hi - lo) / 2 * 1.12 + 40.0)
|
||
return half, center
|
||
|
||
|
||
def q1_solve(stations: Sequence[Sequence[float]]) -> dict[str, Any]:
|
||
"""调用生产求解器，保证图中数字与代码一致。"""
|
||
sys.path.insert(0, str(ROOT / "src"))
|
||
from mcm26b import Measurement, solve # noqa: PLC0415
|
||
res = solve([Measurement(x, y, b, h) for x, y, b, h in stations])
|
||
dc = res["diameter_circle"]
|
||
mec = res["minimum_enclosing_circle"]
|
||
return dict(vertices=np.array(res["vertices"]),
|
||
dc=dict(center=np.array(dc["center"]), radius_m=float(dc["radius_m"])),
|
||
mec=dict(center=np.array(mec["center"]), radius_m=float(mec["radius_m"])))
|
||
|
||
|
||
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 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 seven_bound(a: float = RING_RADIUS_M, radius: float = SEVEN_R) -> float:
|
||
return max(a / math.sqrt(3),
|
||
math.sqrt(radius ** 2 + a ** 2 - math.sqrt(3) * radius * a))
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 各图数据
|
||
# ---------------------------------------------------------------------------
|
||
def data_flow_overall() -> dict[str, Any]:
|
||
"""四个层次带：层标签、每层的框内短标签与语义色名。"""
|
||
bands = [
|
||
("L0", [("区域 Ω：半径 1800 m", 8.0, 26.0),
|
||
("移动 5 m/s · 检测 5 s", 36.5, 26.0),
|
||
("换频 1 s · 清除 5/3 s", 65.0, 26.0),
|
||
("源 10~16 个 · 频道 1~20", 93.5, 24.0)]),
|
||
("L1", [("证据层：候选集与粒子幸存集", 8.0, 34.0),
|
||
("决策层：选点与任务表", 45.0, 34.0),
|
||
("执行层：串行客户端", 82.0, 34.0)]),
|
||
("L2", [("问题 1 定位解算", 8.0, 26.0),
|
||
("问题 2 第二点选择", 36.5, 26.0),
|
||
("问题 3 全向搜索", 65.0, 26.0),
|
||
("问题 4 混合定向", 93.5, 24.0)]),
|
||
("L3", [("发现闭环", 8.0, 34.0),
|
||
("定位闭环", 45.0, 34.0),
|
||
("清除闭环", 82.0, 34.0)]),
|
||
]
|
||
colors = {"L0": "primary", "L1": "primary", "L2": "accent", "L3": "success"}
|
||
return dict(tags=[b[0] for b in bands],
|
||
labels=[[i[0] for i in b[1]] for b in bands],
|
||
colors=[colors[b[0]] for b in bands])
|
||
|
||
|
||
def data_steps_state_machine() -> dict[str, Any]:
|
||
return dict(labels=["IDLE", "倒计时 5 s", "OPEN 接口开放", "CLOSED 断连"],
|
||
edges=["muted", "warn", "success", "accent"])
|
||
|
||
|
||
EXAMPLE_STEPS = (("/enter", 0.0, 0.0, 0.0),
|
||
("/measure ch1", 100.0, 0.0, 5.0),
|
||
("/measure ch2", 0.0, 1.0, 5.0),
|
||
("/clear ch3", 80.0, 0.0, 3.0),
|
||
("/measure ch2", 0.0, 0.0, 5.0),
|
||
("/exit", 0.0, 0.0, 0.0))
|
||
|
||
|
||
def data_timing_attachment() -> dict[str, Any]:
|
||
return dict(labels=[s[0] for s in EXAMPLE_STEPS],
|
||
move=np.array([s[1] for s in EXAMPLE_STEPS]),
|
||
switch_s=np.array([s[2] for s in EXAMPLE_STEPS]),
|
||
action=np.array([s[3] for s in EXAMPLE_STEPS]),
|
||
xlabel="虚拟耗时 / s",
|
||
legend=["移动 L/5", "换频 1 s", "检测/清除动作"])
|
||
|
||
|
||
def data_timing_measured(run: dict[str, Any]) -> dict[str, Any]:
|
||
tl = decompose(run)
|
||
rows = tl["rows"]
|
||
seq = np.arange(1.0, len(rows) + 1.0)
|
||
steps = np.array([[r["move_s"], r["action_s"], r["switch_s"]] for r in rows],
|
||
dtype=float)
|
||
cum = np.cumsum(np.cumsum(steps, axis=1), axis=0) # 移动 / +动作 / +换频
|
||
return dict(seq=seq, cum=cum,
|
||
xlabel="动作序号", ylabel="累计虚拟时间 / s",
|
||
legend=["移动", "检测/清除动作", "换频"])
|
||
|
||
|
||
def data_geo_q1_meetup() -> tuple[dict[str, Any], Path]:
|
||
case = q1_practice_case()
|
||
if case:
|
||
stations = [(m["x"], m["y"], m["bearing_deg"], m["half_angle_deg"])
|
||
for m in case["measurements"]]
|
||
res = case["result"]
|
||
verts = np.array(res["vertices"])
|
||
dc = dict(center=np.array(res["diameter_circle"]["center"]),
|
||
radius_m=float(res["diameter_circle"]["radius_m"]))
|
||
mec = dict(center=np.array(res["minimum_enclosing_circle"]["center"]),
|
||
radius_m=float(res["minimum_enclosing_circle"]["radius_m"]))
|
||
src = case["_path"]
|
||
else:
|
||
src = ROOT / "examples" / "q1" / "triangle.json"
|
||
doc = json.loads(src.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 = q1_solve(stations)
|
||
verts, dc, mec = sol["vertices"], sol["dc"], sol["mec"]
|
||
|
||
zoom_c = verts.mean(axis=0)
|
||
span = float(np.abs(verts - zoom_c).max())
|
||
zoom = max(1.8 * span, 0.8 * mec["radius_m"] + span, 12.0)
|
||
xs = [s[0] for s in stations]
|
||
ys = [s[1] for s in stations]
|
||
half, view_c = overview_window(np.array([[s[0], s[1]] for s in stations]),
|
||
verts, dc, mec)
|
||
reach = 2.4 * half
|
||
wedges = [wedge_polygon((sx, sy), b, hd, reach) for sx, sy, b, hd in stations]
|
||
region = np.array([[zoom_c[0] - zoom, zoom_c[1] - zoom],
|
||
[zoom_c[0] + zoom, zoom_c[1] - zoom],
|
||
[zoom_c[0] + zoom, zoom_c[1] + zoom],
|
||
[zoom_c[0] - zoom, zoom_c[1] + zoom]], dtype=float)
|
||
for sx, sy, b, hd in stations:
|
||
region = clip_wedge(region, (sx, sy), b, hd)
|
||
data = dict(
|
||
stations=np.array([[s[0], s[1]] for s in stations]),
|
||
station_labels=[f"s{i + 1}" for i in range(len(stations))],
|
||
wedges=np.stack(wedges), zoom_wedges=np.stack(wedges),
|
||
zoom_region=region,
|
||
dc_center=dc["center"], dc_radius=dc["radius_m"],
|
||
mec_center=mec["center"], mec_radius=mec["radius_m"],
|
||
view_center=view_c, half_overview=float(half),
|
||
zoom_center=zoom_c, zoom_half=float(zoom),
|
||
legend1=["前向角域 W(s, θ ± 1°)", "定位区域 P"],
|
||
legend2=["直径圆", "最小包围圆 MEC"])
|
||
return data, Path(src)
|
||
|
||
|
||
def data_geo_q1_counterexample() -> dict[str, Any]:
|
||
src = ROOT / "examples" / "q1" / "triangle.json"
|
||
doc = json.loads(src.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 = q1_solve(stations)
|
||
verts, dc, mec = sol["vertices"], sol["dc"], sol["mec"]
|
||
zoom_c = verts.mean(axis=0)
|
||
zoom = max(1.9 * float(np.abs(verts - zoom_c).max()) + mec["radius_m"], 45.0)
|
||
half, view_c = overview_window(np.array([[s[0], s[1]] for s in stations]),
|
||
verts, dc, mec)
|
||
wedges = [wedge_polygon((sx, sy), b, hd, 3600.0) for sx, sy, b, hd in stations]
|
||
region = np.array([[zoom_c[0] - zoom, zoom_c[1] - zoom],
|
||
[zoom_c[0] + zoom, zoom_c[1] - zoom],
|
||
[zoom_c[0] + zoom, zoom_c[1] + zoom],
|
||
[zoom_c[0] - zoom, zoom_c[1] + zoom]], dtype=float)
|
||
for sx, sy, b, hd in stations:
|
||
region = clip_wedge(region, (sx, sy), b, hd)
|
||
return dict(
|
||
stations=np.array([[s[0], s[1]] for s in stations]),
|
||
station_labels=[f"s{i + 1}" for i in range(len(stations))],
|
||
wedges=np.stack(wedges), zoom_wedges=np.stack(wedges),
|
||
zoom_region=region,
|
||
dc_center=dc["center"], dc_radius=dc["radius_m"],
|
||
mec_center=mec["center"], mec_radius=mec["radius_m"],
|
||
view_center=view_c, half_overview=float(half),
|
||
zoom_center=zoom_c, zoom_half=float(zoom),
|
||
legend1=["前向角域 W(s, θ ± 1°)", "定位区域 P"],
|
||
legend2=["直径圆", "最小包围圆 MEC"])
|
||
|
||
|
||
def data_geo_q2() -> dict[str, Any]:
|
||
a = 1500.0 / (2 * math.cos(math.radians(1.0)))
|
||
safe = R_EFF - a
|
||
c = np.array([a * math.cos(math.radians(30.0)), a * math.sin(math.radians(30.0))])
|
||
return dict(omega_r=OMEGA_R, wedge=wedge_polygon((0.0, 0.0), 30.0, 1.0, 1500.0),
|
||
a=float(a), safe=float(safe), c=c, half=1900.0,
|
||
legend=["前向角域 W(s1)", "截断角域外接圆 B(c, a)",
|
||
"安全候选盘 B(c, 1000 − a)"])
|
||
|
||
|
||
def data_geo_q3() -> dict[str, Any]:
|
||
a = RING_RADIUS_M
|
||
stations = np.array([[0.0, 0.0]] + [[a * math.cos(k * math.pi / 3),
|
||
a * math.sin(k * math.pi / 3)] for k in range(6)])
|
||
return dict(omega_r=OMEGA_R, r_eff=R_EFF,
|
||
stations=stations,
|
||
worst=np.array([0.0, -SEVEN_R]),
|
||
half=float(RING_RADIUS_M + R_EFF + 90.0),
|
||
legend=["区域 Ω", "1000 m 接收盘", "覆盖站点", "最差点"])
|
||
|
||
|
||
def data_geo_q4() -> dict[str, Any]:
|
||
plan = load_plan("q4_mincover_v2")
|
||
pts = np.array(plan["points"], dtype=float)
|
||
theta = np.linspace(-math.pi / 2, math.pi / 2, 60)
|
||
half_disk = np.vstack([[0.0, 0.5],
|
||
np.column_stack([0.85 * np.cos(theta),
|
||
0.5 + 0.85 * np.sin(theta)])])
|
||
half = float(max(OMEGA_R, np.abs(pts).max()) * 1.06)
|
||
return dict(omega_r=OMEGA_R, points=pts, half=half, half_disk=half_disk,
|
||
radius_label="覆盖半径 1000 m",
|
||
heading_label="前向半盘 ±90°（朝向未知）",
|
||
legend=["区域 Ω", "31 点最小命中集"])
|
||
|
||
|
||
def data_replay(run: dict[str, Any], *, directional: bool) -> dict[str, Any]:
|
||
tl = decompose(run)
|
||
rows = tl["rows"]
|
||
truth = run["truth"]
|
||
solved = set(solved_channels(run))
|
||
truth_pts, truth_solved, truth_labels = [], [], []
|
||
for src in truth:
|
||
ch = int(src["channel"])
|
||
truth_pts.append([src["x"], src["y"]])
|
||
truth_solved.append(1.0 if ch in solved else 0.0)
|
||
truth_labels.append(f"ch{ch}")
|
||
# 定向源：朝向箭头（长度取该源有效接收半径）+ ±90° 前向半盘
|
||
headings, half_disks = [], []
|
||
for src in truth:
|
||
if not (directional and src.get("directional")):
|
||
continue
|
||
r = float(src.get("r_eff_m") or R_EFF)
|
||
ang = math.radians(float(src.get("heading_deg") or 0.0))
|
||
headings.append([src["x"], src["y"], r * math.cos(ang), r * math.sin(ang)])
|
||
arc = ang + np.linspace(-math.pi / 2, math.pi / 2, 49)
|
||
half_disks.append(np.vstack([[src["x"], src["y"]],
|
||
np.column_stack([src["x"] + r * np.cos(arc),
|
||
src["y"] + r * np.sin(arc)])]))
|
||
half_disk_arr = (np.stack(half_disks) if half_disks
|
||
else np.zeros((0, 50, 2), dtype=float))
|
||
path = [(0.0, 0.0)]
|
||
rays, clears, marks, events = [], [], [], []
|
||
for act, row in zip(run["actions"], rows):
|
||
if act["path"] == "/measure" and row["hit"]:
|
||
bearing = math.radians(float(act["response"].get("svd_deg") or 0.0))
|
||
rays.append([row["x"], row["y"], 220 * math.cos(bearing),
|
||
220 * math.sin(bearing)])
|
||
if act["path"] in ("/measure", "/clear") and row["x"] is not None:
|
||
if (row["x"], row["y"]) != path[-1]:
|
||
path.append((row["x"], row["y"]))
|
||
if act["path"] == "/clear":
|
||
clears.append([row["x"], row["y"], 1.0 if row["cleared"] else 0.0])
|
||
marks.append([row["virtual_s"], 1.0 if row["hit"] else 0.0,
|
||
1.0 if act["path"] == "/clear" and row["cleared"] else 0.0])
|
||
if len(events) < 10:
|
||
if act["path"] == "/measure":
|
||
detail = f"ch{row['channel']} ({row['x']:.0f},{row['y']:.0f})"
|
||
elif act["path"] == "/clear":
|
||
detail = f"ch{row['channel']} {'成功' if row['cleared'] else '未发现'}"
|
||
else:
|
||
detail = "启动会话" if act["path"] == "/enter" else "退出"
|
||
events.append([act["path"], detail, f"{row['virtual_s']:.0f}"])
|
||
legend = ["移动轨迹", "检测示向度", "已清除频道", "未清除频道"]
|
||
if directional:
|
||
legend = ["移动轨迹", "检测示向度", "定向源前向半盘 (±90°)",
|
||
"已清除频道", "未清除频道"]
|
||
# 世界视图窗口：容纳 Ω、轨迹、射线与定向半盘，避免半盘被画到框外
|
||
geom = [np.asarray(truth_pts, dtype=float), np.asarray(path, dtype=float)]
|
||
if rays:
|
||
rays_arr = np.asarray(rays, dtype=float)
|
||
geom.append(rays_arr[:, :2] + rays_arr[:, 2:])
|
||
if half_disk_arr.size:
|
||
geom.append(half_disk_arr.reshape(-1, 2))
|
||
half = float(max(OMEGA_R, np.abs(np.vstack(geom)).max()) * 1.06)
|
||
return dict(
|
||
omega_r=OMEGA_R, half=half,
|
||
out_png=("replay-webui-q4-directional.png" if directional
|
||
else "replay-webui-q3-omnidirectional.png"),
|
||
truth=np.array(truth_pts, dtype=float),
|
||
truth_solved=np.array(truth_solved, dtype=float).reshape(-1, 1),
|
||
truth_labels=truth_labels,
|
||
headings=np.array(headings, dtype=float).reshape(-1, 4),
|
||
half_disks=half_disk_arr,
|
||
path=np.array(path, dtype=float),
|
||
rays=np.array(rays, dtype=float).reshape(-1, 4),
|
||
clears=np.array(clears, dtype=float).reshape(-1, 3),
|
||
marks=np.array(marks, dtype=float).reshape(-1, 3),
|
||
events=events, total_s=float(tl["virtual_s"]), legend=legend)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 产出
|
||
# ---------------------------------------------------------------------------
|
||
def _as_cell(values: Any) -> np.ndarray:
|
||
"""字符串列表 → object 数组（MATLAB cell）；纯数值嵌套列表 → 数值矩阵。"""
|
||
if isinstance(values, str):
|
||
return np.array([values], dtype=object)
|
||
if isinstance(values, (list, tuple)):
|
||
if values and all(isinstance(v, (list, tuple)) and
|
||
all(isinstance(x, (int, float)) for x in v)
|
||
for v in values):
|
||
return np.asarray(values, dtype=float)
|
||
return np.array([_as_cell(v) if isinstance(v, (list, tuple)) else v
|
||
for v in values], dtype=object)
|
||
return np.asarray(values)
|
||
|
||
|
||
def write_mat(name: str, data: dict[str, Any]) -> None:
|
||
"""写 .mat；字符串字段转 object 数组，MATLAB 侧按 cell 读取。"""
|
||
DATA_DIR.mkdir(parents=True, exist_ok=True)
|
||
|
||
def wrap(value: Any, key: str) -> Any: # noqa: ARG001
|
||
if isinstance(value, (str, list, tuple)):
|
||
return _as_cell(value)
|
||
return value
|
||
|
||
payload = {key: wrap(value, key) for key, value in data.items()}
|
||
savemat(DATA_DIR / f"{name}.mat", payload, do_compression=True,
|
||
long_field_names=True)
|
||
|
||
|
||
def prepare(only: Sequence[str] | None) -> dict[str, list[dict[str, Any]]]:
|
||
"""写全部 .mat 数据，返回每张图的数据来源。"""
|
||
want = set(only) if only else None
|
||
runs = discover_runs()
|
||
q3 = pick_run(runs, "q3", 10) if runs else None
|
||
q4 = (pick_case(runs, Q4_REPLAY_CASE) if runs else None) or \
|
||
(pick_run(runs, "q4", 12) if runs else None)
|
||
sources: dict[str, list[dict[str, Any]]] = {}
|
||
|
||
def need(stem: str) -> bool:
|
||
return want is None or stem in want
|
||
|
||
build: list[tuple[str, str, Callable[[], dict[str, Any]]]] = [
|
||
("flow-overall-architecture", "flow_overall", data_flow_overall),
|
||
("steps-state-machine", "steps_state_machine", data_steps_state_machine),
|
||
("timing-attachment-example-105-111-194-199", "timing_attachment",
|
||
data_timing_attachment),
|
||
("geo-q1-diameter-circle-counterexample", "geo_q1_counterexample",
|
||
data_geo_q1_counterexample),
|
||
("geo-q2-guaranteed-reception-disk", "geo_q2", data_geo_q2),
|
||
("geo-q3-seven-point-coverage", "geo_q3", data_geo_q3),
|
||
("geo-q4-minimum-cover-31-points", "geo_q4", data_geo_q4),
|
||
]
|
||
for stem, name, maker in build:
|
||
if not need(stem):
|
||
continue
|
||
write_mat(name, maker())
|
||
sources[stem] = [source(p) for p in SOURCES.get(stem, [])]
|
||
|
||
if need("geo-q1-convex-intersection-zoom"):
|
||
data, src = data_geo_q1_meetup()
|
||
write_mat("geo_q1_meetup", data)
|
||
sources["geo-q1-convex-intersection-zoom"] = [source(src)]
|
||
if q3:
|
||
if need("timing-measured-q3-10-sources"):
|
||
write_mat("timing_measured", data_timing_measured(q3))
|
||
sources["timing-measured-q3-10-sources"] = [
|
||
source(q3["source"], case_code=q3["case_code"], module=q3["module"])]
|
||
if need("replay-webui-q3-omnidirectional"):
|
||
write_mat("replay_q3", data_replay(q3, directional=False))
|
||
sources["replay-webui-q3-omnidirectional"] = [
|
||
source(q3["source"], case_code=q3["case_code"])]
|
||
if q4 and need("replay-webui-q4-directional"):
|
||
write_mat("replay_q4", data_replay(q4, directional=True))
|
||
sources["replay-webui-q4-directional"] = [
|
||
source(q4["source"], case_code=q4["case_code"])]
|
||
return sources
|
||
|
||
|
||
def run_matlab(only: Sequence[str] | None, matlab_bin: Path) -> None:
|
||
"""逐张调用 MATLAB 出图。
|
||
|
||
本机实测：同一 MATLAB 批处理会话里渲染第二张图时图形子系统会崩溃，
|
||
因此每张图单独起一个会话（每次约数秒），并显式 exit。
|
||
"""
|
||
if not matlab_bin.is_file():
|
||
raise SystemExit(f"未找到 MATLAB：{matlab_bin}")
|
||
want = set(only) if only else None
|
||
log = MATLAB_DIR / "matlab_run.log"
|
||
lines: list[str] = []
|
||
for stem, fn, name in FIGURES:
|
||
if want and stem not in want:
|
||
continue
|
||
call = f"{fn}('{name}')" if fn == "fig_replay" else f"{fn}()"
|
||
script = (f"cd('{MATLAB_DIR}'); addpath('{MATLAB_DIR}'); "
|
||
f"{call}; exit")
|
||
proc = subprocess.run([str(matlab_bin), "-batch", script],
|
||
capture_output=True, text=True, check=False)
|
||
ok = (OUT_DIR / f"{stem}.png").is_file()
|
||
lines.append(f"{'OK ' if ok else 'FAIL'} {stem}")
|
||
if not ok:
|
||
tail = [ln for ln in (proc.stdout + proc.stderr).splitlines()
|
||
if ln.strip() and "Keyboard" not in ln and "TISFile" not in ln]
|
||
lines.extend(f" {ln}" for ln in tail[-6:])
|
||
log.write_text("\n".join(lines) + "\n", encoding="utf-8")
|
||
|
||
|
||
def finish(only: Sequence[str] | None, sources: dict[str, list[dict[str, Any]]]) -> None:
|
||
"""删掉历史遗留的 SVG，写 manifest。"""
|
||
for stale in OUT_DIR.glob("*.svg"):
|
||
stale.unlink()
|
||
want = set(only) if only else None
|
||
entries = []
|
||
for stem, fn, _name in FIGURES:
|
||
if want and stem not in want:
|
||
continue
|
||
png = OUT_DIR / f"{stem}.png"
|
||
entry = {"stem": stem, "png": str(png.relative_to(ROOT)),
|
||
"matlab": f"scripts/matlab/{fn}.m",
|
||
"sources": sources.get(stem, [])}
|
||
if not png.is_file():
|
||
entry["missing"] = True
|
||
entries.append(entry)
|
||
manifest = {
|
||
"generator": "scripts/build_figure_data.py（数据）+ scripts/matlab/*.m（MATLAB 出图）",
|
||
"command": "uv run python scripts/build_figure_data.py"
|
||
+ (f" --only {','.join(only)}" if only else ""),
|
||
"manual": "MATLAB 中 cd scripts/matlab 后执行 run_all_figures",
|
||
"figures": entries,
|
||
}
|
||
(OUT_DIR / "manifest.json").write_text(
|
||
json.dumps(manifest, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
|
||
|
||
|
||
def main() -> int:
|
||
parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
|
||
parser.add_argument("--only", default="", help="逗号分隔的图名子集（见 --list）")
|
||
parser.add_argument("--prepare", action="store_true", help="只写 .mat 数据，不调 MATLAB")
|
||
parser.add_argument("--list", action="store_true", help="列出全部图名")
|
||
parser.add_argument("--matlab", default=str(MATLAB_BIN), help="matlab 可执行文件路径")
|
||
args = parser.parse_args()
|
||
names = [stem for stem, _fn, _name in FIGURES]
|
||
if args.list:
|
||
print("\n".join(names))
|
||
return 0
|
||
only = [s.strip() for s in args.only.split(",") if s.strip()]
|
||
unknown = [s for s in only if s not in set(names)]
|
||
if unknown:
|
||
parser.error(f"未知图名 {unknown}")
|
||
|
||
sources = prepare(only or None)
|
||
if args.prepare:
|
||
print(f"[fig] 已写入 {len(list(DATA_DIR.glob('*.mat')))} 个 .mat → "
|
||
f"{DATA_DIR.relative_to(ROOT)}")
|
||
return 0
|
||
targets = only or names
|
||
for stem in targets: # 先清掉旧文件：只认本次生成的
|
||
(OUT_DIR / f"{stem}.png").unlink(missing_ok=True)
|
||
run_matlab(only or None, Path(args.matlab))
|
||
finish(only or None, sources)
|
||
missing = []
|
||
for stem in targets:
|
||
if (OUT_DIR / f"{stem}.png").is_file():
|
||
print(f"[fig] OK {stem}")
|
||
else:
|
||
missing.append(stem)
|
||
print(f"[fig] 缺失 {stem}")
|
||
print(f"[fig] manifest → {(OUT_DIR / 'manifest.json').relative_to(ROOT)}")
|
||
if missing:
|
||
print(f"[fig] {len(missing)} 张未生成：{', '.join(missing)}")
|
||
return 0
|
||
|
||
|
||
if __name__ == "__main__":
|
||
raise SystemExit(main())
|