This commit is contained in:
Horldsence
2026-09-12 16:58:25 +08:00
parent 1fa2264a17
commit 892451dca9
36 changed files with 8767 additions and 0 deletions
+13
View File
@@ -0,0 +1,13 @@
"""问题一:交会定位多边形直径与直径圆覆盖。
- ``optimal.py`` :生产实现 ``mcm26b.geometry.localization.solve``(半平面交 +
顶点对精确直径 + Welzl 最小包围圆 + 覆盖反例见证);
- ``baseline.py`` :基线实现(网格采样 + 轴对齐包围盒代理直径/覆盖);
- ``reference.py``:独立精确参考解(大框半平面裁剪 + 暴力直径/最小包围圆);
- ``cases.py`` :固定示例(``examples/q1``)与随机算例生成;
- ``compare.py`` :同案配对对比并落盘。
"""
from . import baseline, cases, optimal, reference
__all__ = ["optimal", "baseline", "reference", "cases"]
+98
View File
@@ -0,0 +1,98 @@
"""问题一基线实现:网格采样 + 轴对齐包围盒(bbox)代理。
这是"先上手"的朴素做法,常见于直接对可行点打网格再读包围盒:
1. 在固定大框上按步长采样,保留满足全部角域约束的点;
2. 直径用包围盒对角线代理 —— 只要区域是斜的/细长的就会系统性高估;
3. 覆盖判断用"包围盒中心 + 代理直径"的圆 —— 因为包围盒自身对其中心对称,
该判据几乎恒为"能覆盖",从而丢掉问题一真正要考的反例。
本实现刻意保留这些结构性缺陷,作为``optimal.py``的对照;它不读取最优实现的
任何中间量,两条代码路径完全独立。
"""
from __future__ import annotations
import math
from time import perf_counter
import numpy as np
from mcm26b.geometry.localization import Measurement, halfplanes
from project.common import Method
METHOD = Method(
key="baseline",
name="网格采样 + bbox 对角线代理",
kind="baseline",
note="轴对齐包围盒近似直径与覆盖,O(盒面积/步长²),与参考解同分辨率",
)
BOX_M = 2600.0
STEP_M = 5.0
def solve_case(case: dict) -> dict:
start = perf_counter()
measurements = [Measurement(row["x"], row["y"], row["bearing_deg"],
row.get("half_angle_deg", 1.0))
for row in case["measurements"]]
planes = [plane for m in measurements for plane in halfplanes(m)]
xs = np.arange(-BOX_M, BOX_M + STEP_M, STEP_M)
grid_x, grid_y = np.meshgrid(xs, xs) # 行=y,列=x
mask = np.ones(grid_x.shape, dtype=bool)
for a, b, c in planes:
mask &= (a * grid_x + b * grid_y <= c)
result = _evaluate(mask, xs)
result["method"] = METHOD.key
result["runtime_ms"] = round((perf_counter() - start) * 1000.0, 4)
result["algorithm"] = METHOD.name
return result
def _evaluate(mask: np.ndarray, axis: np.ndarray) -> dict:
if not mask.any():
return _empty("empty")
# 可行点触到采样框边界 => 区域至少延伸到框外,判无界
if (mask[0, :].any() or mask[-1, :].any()
or mask[:, 0].any() or mask[:, -1].any()):
return _empty("unbounded")
rows, cols = np.nonzero(mask)
x_min, x_max = axis[cols].min(), axis[cols].max()
y_min, y_max = axis[rows].min(), axis[rows].max()
width, height = x_max - x_min, y_max - y_min
diameter = math.hypot(width, height) # 包围盒对角线代理
center = ((x_min + x_max) / 2.0, (y_min + y_max) / 2.0)
px = axis[cols]
py = axis[rows]
max_distance = float(np.max(np.hypot(px - center[0], py - center[1])))
return {
"status": "polygon" if mask.sum() >= 2 else "point",
"diameter_m": diameter,
"diameter_infinite": False,
"covers": bool(max_distance <= diameter / 2.0 + 1e-9),
"mec_radius_m": None,
"max_excess_m": round(max_distance - diameter / 2.0, 6),
"outside_vertex_count": None,
"vertex_count": int(mask.sum()),
"bbox_width_m": round(width, 3),
"bbox_height_m": round(height, 3),
"runtime_ms": None,
}
def _empty(status: str) -> dict:
return {
"status": status,
"diameter_m": None,
"diameter_infinite": status == "unbounded",
"covers": None,
"mec_radius_m": None,
"max_excess_m": None,
"outside_vertex_count": None,
"vertex_count": 0,
}
+67
View File
@@ -0,0 +1,67 @@
"""问题一算例:官方示例目录 ``examples/q1`` + 随机交会算例。"""
from __future__ import annotations
import json
import math
import random
from project.common import ROOT
EXAMPLE_DIR = ROOT / "examples" / "q1"
EXAMPLE_IDS = ("normal", "triangle", "poor", "segment", "empty")
def load_examples() -> list[dict]:
"""读取固定算例;缺文件时跳过,保证离线可跑。"""
cases = []
for key in EXAMPLE_IDS:
path = EXAMPLE_DIR / f"{key}.json"
if not path.exists():
continue
raw = json.loads(path.read_text(encoding="utf-8-sig"))
cases.append({
"id": raw.get("id", key),
"name": raw.get("name", key),
"source": "examples/q1",
"measurements": [
{"x": float(row["x"]), "y": float(row["y"]),
"bearing_deg": float(row["bearing_deg"]),
"half_angle_deg": float(row.get("half_angle_deg", 1.0))}
for row in raw["measurements"]
],
})
return cases
def random_case(seed: int) -> dict:
"""随机三点/四点交会:目标在 Ω 内,测站分布在不同象限。
仅生成"区域有界且非退化"的良态算例,避免参考网格因细长区域失效;
生成后由 ``compare`` 再用最优实现复核状态。
"""
rng = random.Random(seed)
target = (rng.uniform(-900, 900), rng.uniform(-900, 900))
count = rng.choice((3, 3, 4, 5))
stations = []
base_angle = rng.uniform(0, 360)
for index in range(count):
angle = base_angle + 360 * index / count + rng.uniform(-25, 25)
radius = rng.uniform(500, 1100)
stations.append((target[0] + radius * math.cos(math.radians(angle)),
target[1] + radius * math.sin(math.radians(angle))))
measurements = []
for (sx, sy) in stations:
true_bearing = math.degrees(math.atan2(target[1] - sy, target[0] - sx)) % 360.0
# 观测误差固定在同一地点、落在 ±1° 内;读数保留两位小数
error = rng.uniform(-0.9, 0.9)
measurements.append({"x": round(sx, 2), "y": round(sy, 2),
"bearing_deg": round((true_bearing + error) % 360.0, 2),
"half_angle_deg": 1.0})
return {"id": f"rand-{seed}", "name": "随机交会", "source": "synthetic",
"target": [round(target[0], 2), round(target[1], 2)],
"measurements": measurements}
def random_cases(count: int, base_seed: int = 7000) -> list[dict]:
return [random_case(base_seed + index) for index in range(count)]
+150
View File
@@ -0,0 +1,150 @@
"""问题一同案配对对比:最优 vs 基线 vs 独立参考解。"""
from __future__ import annotations
from project.common import (Method, RESULTS_DIR, markdown_table, mean,
ratio, write_json, write_text)
from project.q1 import baseline, cases, optimal, reference
METHODS: tuple[Method, ...] = (optimal.METHOD, baseline.METHOD, reference.METHOD)
FINITE = ("point", "segment", "polygon")
def _run(case: dict) -> dict:
opt = optimal.solve_case(case)
base = baseline.solve_case(case)
ref = reference.solve_case(case)
finite = ref["status"] in FINITE and ref["diameter_m"] is not None
return {
"id": case["id"],
"name": case.get("name", case["id"]),
"source": case.get("source", "synthetic"),
"measurements": len(case["measurements"]),
"status_ref": ref["status"],
"status_opt": opt["status"],
"status_base": base["status"],
"status_opt_ok": opt["status"] == ref["status"],
"status_base_ok": base["status"] == ref["status"],
"diameter_ref_m": ref["diameter_m"],
"diameter_opt_m": opt["diameter_m"],
"diameter_base_m": base["diameter_m"],
"abs_err_opt_m": (round(abs(opt["diameter_m"] - ref["diameter_m"]), 4)
if finite and opt["diameter_m"] is not None else None),
"abs_err_base_m": (round(abs(base["diameter_m"] - ref["diameter_m"]), 4)
if finite and base["diameter_m"] is not None else None),
"covers_ref": ref["covers"],
"covers_opt": opt["covers"],
"covers_base": base["covers"],
"covers_opt_ok": (opt["covers"] == ref["covers"]
if ref["covers"] is not None and opt["covers"] is not None
else None),
"covers_base_ok": (base["covers"] == ref["covers"]
if ref["covers"] is not None and base["covers"] is not None
else None),
"mec_radius_opt_m": opt["mec_radius_m"],
"mec_radius_ref_m": ref["mec_radius_m"],
"outside_vertices_opt": opt["outside_vertex_count"],
"runtime_opt_ms": opt["runtime_ms"],
"runtime_base_ms": base["runtime_ms"],
"runtime_ref_ms": ref["runtime_ms"],
"reference_resolution_m": ref["resolution_m"],
}
def _aggregate(rows: list[dict]) -> dict:
polygon_rows = [r for r in rows if r["status_ref"] == "polygon"]
opt_errs = [r["abs_err_opt_m"] for r in polygon_rows
if r["abs_err_opt_m"] is not None]
base_errs = [r["abs_err_base_m"] for r in polygon_rows
if r["abs_err_base_m"] is not None]
cov_opt = [r["covers_opt_ok"] for r in rows if r["covers_opt_ok"] is not None]
cov_base = [r["covers_base_ok"] for r in rows if r["covers_base_ok"] is not None]
return {
"cases": len(rows),
"diameter_compared": len(opt_errs),
"status_accuracy_optimal": ratio(sum(r["status_opt_ok"] for r in rows), len(rows)),
"status_accuracy_baseline": ratio(sum(r["status_base_ok"] for r in rows), len(rows)),
"diameter_mae_optimal_m": mean(opt_errs),
"diameter_mae_baseline_m": mean(base_errs),
"diameter_max_err_optimal_m": round(max(opt_errs), 4) if opt_errs else 0.0,
"diameter_max_err_baseline_m": round(max(base_errs), 4) if base_errs else 0.0,
"coverage_compared": len(cov_opt),
"coverage_accuracy_optimal": ratio(sum(cov_opt), len(cov_opt)),
"coverage_accuracy_baseline": ratio(sum(cov_base), len(cov_base)),
"runtime_optimal_mean_ms": mean(r["runtime_opt_ms"] for r in rows),
"runtime_baseline_mean_ms": mean(r["runtime_base_ms"] for r in rows),
"runtime_reference_mean_ms": mean(r["runtime_ref_ms"] for r in rows),
}
def _render(rows: list[dict], aggregate: dict, title: str) -> str:
headers = ["算例", "来源", "参考状态", "最优状态", "基线状态",
"D参考(m)", "D最优(m)", "D基线(m)", "|ΔD|最优", "|ΔD|基线",
"覆盖参考", "覆盖最优", "覆盖基线"]
body = []
for r in rows:
body.append([
r["name"], r["source"], r["status_ref"], r["status_opt"], r["status_base"],
_fmt(r["diameter_ref_m"]), _fmt(r["diameter_opt_m"]),
_fmt(r["diameter_base_m"]),
_fmt(r["abs_err_opt_m"]), _fmt(r["abs_err_base_m"]),
_bool(r["covers_ref"]), _bool(r["covers_opt"]), _bool(r["covers_base"]),
])
text = [f"# {title}", "",
"参考解为独立的大框半平面裁剪 + 精确直径/MEC;",
"固定示例见 `examples/q1`,解析期望由 `tests/test_geometry*.py` 固化。",
"", markdown_table(headers, body), "",
"## 汇总", ""]
summary_rows = [
["状态判定准确率", _pct(aggregate["status_accuracy_optimal"]),
_pct(aggregate["status_accuracy_baseline"])],
["直径平均绝对误差 (m)", _fmt(aggregate["diameter_mae_optimal_m"], 3),
_fmt(aggregate["diameter_mae_baseline_m"], 3)],
["直径最大误差 (m)", _fmt(aggregate["diameter_max_err_optimal_m"], 3),
_fmt(aggregate["diameter_max_err_baseline_m"], 3)],
["覆盖判定准确率", _pct(aggregate["coverage_accuracy_optimal"]),
_pct(aggregate["coverage_accuracy_baseline"])],
["平均耗时 (ms)", _fmt(aggregate["runtime_optimal_mean_ms"], 3),
_fmt(aggregate["runtime_baseline_mean_ms"], 3)],
]
text.append(markdown_table(["指标", "最优", "基线"], summary_rows))
return "\n".join(text) + "\n"
def _fmt(value, digits: int = 2) -> str:
if value is None:
return "—"
if isinstance(value, bool):
return "是" if value else "否"
if isinstance(value, float):
return f"{value:.{digits}f}"
return str(value)
def _bool(value) -> str:
return "—" if value is None else ("是" if value else "否")
def _pct(value: float) -> str:
return f"{value*100:.1f}%"
def run(quick: bool = False) -> dict:
example_rows = [_run(case) for case in cases.load_examples()]
random_rows = [_run(case) for case in
cases.random_cases(12 if quick else 40)]
payload = {
"question": 1,
"methods": [m.__dict__ for m in METHODS],
"examples": example_rows,
"random": random_rows,
"aggregate_examples": _aggregate(example_rows),
"aggregate_random": _aggregate(random_rows),
}
write_json(RESULTS_DIR / "q1.json", payload)
text = _render(example_rows, payload["aggregate_examples"],
"问题一:固定示例对比")
text += "\n" + _render(random_rows, payload["aggregate_random"],
"问题一:随机交会算例对比")
write_text(RESULTS_DIR / "q1.md", text)
return payload
+59
View File
@@ -0,0 +1,59 @@
"""问题一最优实现:复用生产几何内核 ``mcm26b.geometry.localization.solve``。
算法(已由 ``tests/test_geometry.py`` / ``tests/test_geometry_certificates.py``
固化,本模块只做输入输出适配,不重写):
1. 每条带 ±δ 误差的示向度写成前向角域的三条半平面约束(含 0°/360° 跨界);
2. 枚举全部边界直线交点,保留满足全部约束者,单调链求凸包;
3. 检查是否存在非零衰退方向判定无界;否则枚举凸包顶点对取精确直径;
4. 用直径两端点中点圆检验覆盖,并给出圆外顶点、最大超出距离等反例见证;
5. Welzl 增量法求最小包围圆,给出 20 m 清除余量证书。
"""
from __future__ import annotations
from time import perf_counter
from mcm26b.geometry.localization import Measurement, solve
from project.common import Method
METHOD = Method(
key="optimal",
name="精确半平面交 + 顶点对直径 + MEC",
kind="optimal",
note="生产内核 mcm26b.geometry.localization.solve",
)
def solve_case(case: dict) -> dict:
"""输入 case["measurements"],输出归一化结果。"""
measurements = [Measurement(row["x"], row["y"], row["bearing_deg"],
row.get("half_angle_deg", 1.0))
for row in case["measurements"]]
start = perf_counter()
result = solve(measurements)
elapsed_ms = (perf_counter() - start) * 1000.0
return normalize(result, elapsed_ms)
def normalize(result: dict, elapsed_ms: float) -> dict:
status = result["status"]
diameter = result.get("diameter_m")
circle = result.get("diameter_circle") or {}
mec = result.get("minimum_enclosing_circle") or {}
covers = circle.get("covers") if status in ("point", "segment", "polygon") else None
return {
"method": METHOD.key,
"status": status,
"diameter_m": diameter,
"diameter_infinite": result.get("diameter_infinite"),
"covers": covers,
"mec_radius_m": mec.get("radius_m"),
"max_excess_m": circle.get("max_excess_m"),
"outside_vertex_count": len(circle.get("outside_vertex_indices", []))
if circle else None,
"vertex_count": len(result.get("vertices", [])),
"runtime_ms": round(elapsed_ms, 4),
"algorithm": METHOD.name,
}
+178
View File
@@ -0,0 +1,178 @@
"""问题一独立参考解:大框 Sutherland–Hodgman 半平面裁剪 + 精确直径/MEC。
与生产实现 ``mcm26b.geometry.localization.solve`` 的独立性:
- 生产实现枚举"半平面边界直线两两交点"再筛可行点;
- 本参考解从一个大方框出发,逐条半平面做 Sutherland–Hodgman 裁剪,
是另一条构造路径,且这里的半平面系数、直径、最小包围圆都在本文件内
独立编码(不复用生产实现的 ``halfplanes`` / ``diameter`` / ``mec``)。
只要区域有界,大框(±1e7 m)裁剪得到的顶点即精确多边形;顶点触到大框
边界则判无界。最小包围圆用"枚举所有 1/2/3 点支撑圆"的暴力算法,直径用
顶点对暴力枚举,二者都是 O(m²)/O(m³) 的朴素精确算法,m 为多边形顶点数。
"""
from __future__ import annotations
import math
from time import perf_counter
from project.common import Method
METHOD = Method(
key="reference",
name="大框半平面裁剪 + 精确直径/MEC",
kind="reference",
note="独立于生产顶点枚举的精确参考解",
)
BOX_M = 1.0e7
_PARALLEL_EPS = 1e-12
def _planes(case: dict) -> list[tuple[float, float, float]]:
"""独立重写 a*x+b*y<=c 的半平面构造(与题面角域定义一致)。"""
planes = []
for row in case["measurements"]:
half = math.radians(row.get("half_angle_deg", 1.0))
theta = math.radians(row["bearing_deg"] % 360.0)
sx, sy = row["x"], row["y"]
for angle, sign in ((theta - half, 1.0), (theta + half, -1.0)):
a = sign * math.sin(angle)
b = -sign * math.cos(angle)
planes.append((a, b, a * sx + b * sy))
# 前向约束:不允许把角域取成反向无限双楔
a, b = -math.cos(theta), -math.sin(theta)
planes.append((a, b, a * sx + b * sy))
return planes
def _clip(polygon, a, b, c, tol=1e-7):
"""Sutherland–Hodgman:保留 a*x+b*y<=c 的一侧。"""
if not polygon:
return []
output = []
count = len(polygon)
for index in range(count):
p, q = polygon[index], polygon[(index + 1) % count]
fp = a * p[0] + b * p[1] - c
fq = a * q[0] + b * q[1] - c
p_in, q_in = fp <= tol, fq <= tol
if p_in:
output.append(p)
if p_in != q_in:
ratio = fp / (fp - fq)
output.append((p[0] + ratio * (q[0] - p[0]),
p[1] + ratio * (q[1] - p[1])))
deduped = []
for point in output:
if not deduped or math.dist(point, deduped[-1]) > 1e-9:
deduped.append(point)
if len(deduped) > 1 and math.dist(deduped[0], deduped[-1]) <= 1e-9:
deduped.pop()
return deduped
def _brute_force_diameter(points):
best = 0.0
pair = (points[0], points[0])
for i, pi in enumerate(points):
for pj in points[i + 1:]:
value = math.dist(pi, pj)
if value > best:
best, pair = value, (pi, pj)
return best, pair
def _circumcircle(a, b, c):
ux, uy = b[0] - a[0], b[1] - a[1]
vx, vy = c[0] - a[0], c[1] - a[1]
det = 2.0 * (ux * vy - uy * vx)
if abs(det) <= _PARALLEL_EPS * max(
math.hypot(ux, uy) * math.hypot(vx, vy), 1e-300):
return None
uu, vv = ux * ux + uy * uy, vx * vx + vy * vy
center = (a[0] + (uu * vy - vv * uy) / det,
a[1] + (ux * vv - vx * uu) / det)
return center, math.dist(center, a)
def _brute_force_mec(points):
"""枚举 1/2/3 点支撑圆的精确最小包围圆(朴素但无迭代)。"""
xs = [p[0] for p in points]
ys = [p[1] for p in points]
best_c = ((min(xs) + max(xs)) / 2.0, (min(ys) + max(ys)) / 2.0)
best_r = max(math.dist(best_c, p) for p in points)
candidates: list[tuple[tuple[float, float], float]] = []
for i, p in enumerate(points):
for q in points[i + 1:]:
candidates.append((((p[0] + q[0]) / 2.0, (p[1] + q[1]) / 2.0),
math.dist(p, q) / 2.0))
for i, p in enumerate(points):
for j in range(i + 1, len(points)):
for k in range(j + 1, len(points)):
circle = _circumcircle(p, points[j], points[k])
if circle is not None:
candidates.append(circle)
for center, radius in candidates:
if radius >= best_r:
continue
if all(math.dist(center, p) <= radius + 1e-7 for p in points):
best_c, best_r = center, radius
return best_c, best_r
def solve_case(case: dict) -> dict:
start = perf_counter()
planes = _planes(case)
polygon = [(-BOX_M, -BOX_M), (BOX_M, -BOX_M), (BOX_M, BOX_M), (-BOX_M, BOX_M)]
for a, b, c in planes:
polygon = _clip(polygon, a, b, c)
if not polygon:
return _finish(_degenerate("empty"), start)
boundary_tol = 1e-6 * BOX_M
if any(abs(point[0]) >= BOX_M - boundary_tol
or abs(point[1]) >= BOX_M - boundary_tol for point in polygon):
return _finish(_degenerate("unbounded"), start)
if len(polygon) == 1:
summary = {"status": "point", "diameter_m": 0.0, "covers": True,
"mec_radius_m": 0.0, "max_excess_m": 0.0,
"vertex_count": 1}
else:
value, (a, b) = _brute_force_diameter(polygon)
if len(polygon) == 2:
summary = {"status": "segment", "diameter_m": value, "covers": True,
"mec_radius_m": value / 2.0, "max_excess_m": 0.0,
"vertex_count": 2}
else:
_mec_center, radius = _brute_force_mec(polygon)
diameter_center = ((a[0] + b[0]) / 2.0, (a[1] + b[1]) / 2.0)
max_distance = max(math.dist(diameter_center, p) for p in polygon)
summary = {
"status": "polygon",
"diameter_m": value,
"covers": bool(radius <= value / 2.0 + 1e-7),
"mec_radius_m": radius,
"max_excess_m": round(max_distance - value / 2.0, 9),
"vertex_count": len(polygon),
}
return _finish(summary, start)
def _degenerate(status: str) -> dict:
return {"status": status, "diameter_m": None, "covers": None,
"mec_radius_m": None, "max_excess_m": None, "vertex_count": 0}
def _finish(summary: dict, start: float) -> dict:
result = {
"method": METHOD.key,
"algorithm": METHOD.name,
"resolution_m": 0.0,
"diameter_infinite": summary["status"] == "unbounded",
"outside_vertex_count": None,
"runtime_ms": round((perf_counter() - start) * 1000.0, 4),
}
result.update(summary)
return result