"""问题一独立参考解：大框 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