"""P4b 网格粒子滤波验收：真源保留定理、证书正确性与查询原语。 核心不变式（docs/algo-3-bayes.md §2）： - 真源保留：任意观测序列（|误差| <= 1°、r_eff ∈ [1000,1500]、题面语义） 后，真源所在格必在幸存集中——由此 MEC 证书与缺席证明才严格成立； - 缺席证明：无源频道在排除盘覆盖 Ω 后幸存集必须清空； - 查询：nearest_alive_distance 与暴力扫描一致；version 单调递增。 """ import math import unittest from dataclasses import replace from random import Random from mcm26b.config import ProblemConfig from mcm26b.geometry.localization import smallest_enclosing_circle from mcm26b.particles import ChannelParticles, angular_gap from mcm26b.protocol.rules import DIRECTION, NEAR, NO_SIGNAL from mcm26b.planner.seven import seven_points CONFIG = ProblemConfig() rng_shared = Random(99) def _make_particle(spacing_m=25.0): return ChannelParticles(CONFIG, spacing_m=spacing_m) def _cell_index(particle, point): i = int(math.floor((point[0] - particle.x0) / particle.spacing_m)) j = int(math.floor((point[1] - particle.x0) / particle.spacing_m)) return i * particle.width + j def _observe(source, station, rng, cleared=False): """按题面语义产生一次读数（误差在 [-1°,1°] 内随机，含两位小数舍入）。""" if cleared: return NO_SIGNAL, None dist = math.dist(station, (source[0], source[1])) if dist > source[2]: return NO_SIGNAL, None if dist <= CONFIG.field.near_radius_m: return NEAR, None true_bearing = math.degrees(math.atan2(source[1] - station[1], source[0] - station[0])) error = (2.0 * rng.random() - 1.0) * CONFIG.field.bearing_error_half_deg return DIRECTION, round(true_bearing + error, 2) % 360.0 class TruthPreservationTests(unittest.TestCase): def test_truth_cell_survives_arbitrary_sequences(self): """随机源 × 随机测点 × 随机合法误差：真源所在格永不被删。""" rng = Random(20260911) for trial in range(8): particle = _make_particle() radius = CONFIG.field.omega_radius_m * math.sqrt(rng.random()) angle = rng.uniform(0.0, 2.0 * math.pi) source = (radius * math.cos(angle), radius * math.sin(angle), rng.uniform(1000.0, 1500.0)) truth_cell = _cell_index(particle, source) self.assertIn(truth_cell, particle.alive) stations = [(rng.uniform(-1800, 1800), rng.uniform(-1800, 1800)) for _ in range(25)] + seven_points() for station in stations: result, svd = _observe(source, station, rng) particle.apply_measure(station, result, svd) self.assertIn(truth_cell, particle.alive, f"trial {trial} station {station}: 真源格被删") self.assertTrue(particle.alive, "含源频道幸存集不允许为空") # 真源距幸存集 <= Voronoi 裕量（清除证书的几何依据） nearest = min(math.dist(p, source[:2]) for p in (particle.cell_center(i) for i in particle.alive)) self.assertLessEqual(nearest, particle.margin_m + 1e-6) def test_absence_proof_when_exclusions_cover_domain(self): """无源频道：七站 no_signal 的排除盘覆盖 Ω，幸存集必须清空。""" particle = _make_particle() for station in seven_points(): particle.apply_measure(station, NO_SIGNAL, None) self.assertEqual(len(particle), 0, "七站覆盖半径 968.9m < 1000−m，缺席必须可证") def test_direction_result_never_empties_feasible_set(self): """direction 分支（角域 ∩ B(s,1500)）不允许把含源频道剪空。""" rng = Random(7) particle = _make_particle() source = (900.0, -700.0, 1200.0) for _ in range(40): station = (rng.uniform(-1700, 1700), rng.uniform(-1700, 1700)) result, svd = _observe(source, station, rng) particle.apply_measure(station, result, svd) if result == DIRECTION: self.assertTrue(particle.alive) class CertificateTests(unittest.TestCase): def test_clear_certificate_implies_clear_success(self): """证书蕴含关系：MEC(幸存)+裕量 <= 20−0.5 ⟹ 真源距清除点 <= 20。 用小测试场（omega=200m）+ 细网格（h=5m）使证书可达：两次方向 交会后 near 读数把可行集收缩到 5m 盘，网格证书触发，然后按 /clear 语义验证真源到证书圆心的距离不超过 20m。 """ small = replace(CONFIG, field=replace(CONFIG.field, omega_radius_m=200.0)) particle = ChannelParticles(small, spacing_m=5.0) limit = (CONFIG.field.clear_radius_m - particle.margin_m - CONFIG.search.clear_margin_m) source = (60.0, 40.0, 1500.0) # r_eff 取题面上界；小场内必收得到信号 # 两个远站方向交会（读数带 ±1° 合法误差），可行集收缩为小菱形 for station in ((-150.0, -40.0), (40.0, -150.0)): result, svd = _observe(source, station, rng_shared) particle.apply_measure(station, result, svd) self.assertTrue(particle.alive) # near：源距站 3m（<=5m），可行集收缩到 B(station, 5+裕量) station = (source[0] + 3.0, source[1]) particle.apply_measure(station, NEAR, None) center, radius = particle.mec() self.assertLessEqual(radius, limit, "5m 盘收缩后网格证书必须可触发") self.assertLessEqual(radius + particle.margin_m + CONFIG.search.clear_margin_m, CONFIG.field.clear_radius_m) self.assertLessEqual(math.dist(center, source[:2]), CONFIG.field.clear_radius_m, "证书触发时真源到清除点必须 <= 20m") def test_mec_contains_all_alive_cells(self): particle = _make_particle() particle.apply_measure((400.0, 300.0), DIRECTION, 25.0) center, radius = particle.mec() for point in particle.alive_points(): self.assertLessEqual(math.dist(point, center), radius + 1e-6) class QueryPrimitiveTests(unittest.TestCase): def test_nearest_alive_distance_matches_bruteforce(self): rng = Random(3) particle = _make_particle() for _ in range(6): station = (rng.uniform(-1500, 1500), rng.uniform(-1500, 1500)) particle.apply_measure(station, NO_SIGNAL, None) for _ in range(30): probe = (rng.uniform(-1900, 1900), rng.uniform(-1900, 1900)) limit = rng.uniform(100.0, 2500.0) got = particle.nearest_alive_distance(probe, limit) brute = min((math.dist(p, probe) for p in particle.alive_points()), default=None) if brute is not None and brute <= limit: self.assertIsNotNone(got) self.assertAlmostEqual(got, brute, delta=1e-6) else: self.assertIsNone(got, f"limit={limit} brute={brute} 应返回 None") def test_version_and_mec_cache(self): particle = _make_particle() before = particle.version particle.apply_measure((0.0, 0.0), NO_SIGNAL, None) self.assertGreater(particle.version, before) mec_once = particle.mec() self.assertIs(particle.mec(), mec_once, "version 不变时 MEC 应命中缓存") def test_angular_gap(self): self.assertAlmostEqual(angular_gap(350.0, 10.0), 20.0) self.assertAlmostEqual(angular_gap(10.0, 350.0), 20.0) self.assertAlmostEqual(angular_gap(45.0, 45.0), 0.0) if __name__ == "__main__": unittest.main()