#!/usr/bin/env python3
"""Audit cave primitive bounds against passage envelopes, allowing chamber joints only."""
from __future__ import annotations

import math
from typing import Any


INFLUENCE_REACH_M = 1.053334
SHELL_THICKNESS_M = 0.75


def _point(control: dict[str, Any], scale: float) -> tuple[float, float, float, float]:
    x, y, z = control["position_q"]
    return (x * scale, y * scale, z * scale, control["radius_q"] * scale)


def _knot_span(a, b) -> float:
    return max(1.0e-6, math.sqrt(math.sqrt((b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2)))


def _evaluate(controls, span: int, t: float):
    p1, p2 = controls[span], controls[span + 1]
    p0 = controls[span - 1] if span else controls[0]
    p3 = controls[min(span + 2, len(controls) - 1)]
    dt0, dt1, dt2 = _knot_span(p0, p1), _knot_span(p1, p2), _knot_span(p2, p3)
    t2, t3 = t * t, t * t * t
    result = []
    for axis in (0, 1, 3):
        m1 = dt1 * ((p1[axis] - p0[axis]) / dt0
                    - (p2[axis] - p0[axis]) / (dt0 + dt1)
                    + (p2[axis] - p1[axis]) / dt1)
        m2 = dt1 * ((p2[axis] - p1[axis]) / dt1
                    - (p3[axis] - p1[axis]) / (dt1 + dt2)
                    + (p3[axis] - p2[axis]) / dt2)
        value = ((2 * t3 - 3 * t2 + 1) * p1[axis]
                 + (-2 * t3 + 3 * t2) * p2[axis]
                 + (t3 - 2 * t2 + t) * m1 + (t3 - t2) * m2)
        result.append(value)
    radius = min(max(result[2], min(p1[3], p2[3])), max(p1[3], p2[3]))
    return result[0], result[1], p1[2] + (p2[2] - p1[2]) * t, radius


def _segments(passage: dict[str, Any], scale: float):
    controls = [_point(row, scale) for row in passage["control_points"]]
    route_end = min(max(1, int(passage.get("lower_mouth_control_index", len(controls) - 1))),
                    len(controls) - 1)
    along = 0.0
    for span in range(route_end):
        a, b = controls[span], controls[span + 1]
        chord = math.hypot(b[0] - a[0], b[1] - a[1])
        steps = max(1, math.ceil(chord / 2.0))
        prior = _evaluate(controls, span, 0.0)
        for step in range(1, steps + 1):
            current = _evaluate(controls, span, step / steps)
            length = math.dist(prior[:3], current[:3])
            yield prior, current, along, max(prior[3], current[3]) + INFLUENCE_REACH_M
            along += length
            prior = current


def _distance_squared(a, b, bounds, t: float) -> float:
    total = 0.0
    for axis in range(3):
        pos = a[axis] + (b[axis] - a[axis]) * t
        lo, hi = bounds[axis]
        gap = lo - pos if pos < lo else pos - hi if pos > hi else 0.0
        total += gap * gap
    return total


def _overlap_interval(a, b, bounds, radius: float):
    cuts = {0.0, 1.0}
    for axis in range(3):
        delta = b[axis] - a[axis]
        if abs(delta) < 1.0e-15:
            continue
        for edge in bounds[axis]:
            t = (edge - a[axis]) / delta
            if 0.0 < t < 1.0:
                cuts.add(t)
    cuts = sorted(cuts)
    best_t, best_d2 = 0.0, math.inf
    for lo, hi in zip(cuts, cuts[1:]):
        mid = (lo + hi) * 0.5
        quadratic = linear = 0.0
        for axis in range(3):
            delta = b[axis] - a[axis]
            pos = a[axis] + delta * mid
            low, high = bounds[axis]
            edge = low if pos < low else high if pos > high else None
            if edge is not None:
                quadratic += delta * delta
                linear += delta * (a[axis] - edge)
        candidate = min(hi, max(lo, -linear / quadratic)) if quadratic else lo
        d2 = _distance_squared(a, b, bounds, candidate)
        if d2 < best_d2:
            best_t, best_d2 = candidate, d2
    for t in (0.0, 1.0):
        d2 = _distance_squared(a, b, bounds, t)
        if d2 < best_d2:
            best_t, best_d2 = t, d2
    if best_d2 > radius * radius:
        return None

    radius2 = radius * radius
    if _distance_squared(a, b, bounds, 0.0) <= radius2:
        left = 0.0
    else:
        low, high = 0.0, best_t
        for _ in range(42):
            mid = (low + high) * 0.5
            if _distance_squared(a, b, bounds, mid) <= radius2:
                high = mid
            else:
                low = mid
        left = high
    if _distance_squared(a, b, bounds, 1.0) <= radius2:
        right = 1.0
    else:
        low, high = best_t, 1.0
        for _ in range(42):
            mid = (low + high) * 0.5
            if _distance_squared(a, b, bounds, mid) <= radius2:
                low = mid
            else:
                high = mid
        right = low
    return left, right, best_t


def _inside_chamber(point, chamber: dict[str, Any], scale: float) -> bool:
    cx, cy, cz = (value * scale for value in chamber["floor_center_q"])
    radius = chamber["radius_q"] * scale
    clear = chamber["clear_height_q"] * scale
    gx, gy = (value / 65536.0 for value in chamber["floor_gradient_q16"])
    floor = cz + gx * (point[0] - cx) + gy * (point[1] - cy)
    chamber_distance = max(math.hypot(point[0] - cx, point[1] - cy) - radius,
                           floor - point[2], point[2] - (floor + clear))
    return chamber_distance <= SHELL_THICKNESS_M


def _primitive_bounds(row: dict[str, Any]):
    return tuple((row["bounds_min_um"][i] / 1.0e6,
                  row["bounds_max_um"][i] / 1.0e6) for i in range(3))


def _bezier(points, t: float):
    u = 1.0 - t
    return tuple(u**3 * points[0][axis] + 3*u*u*t * points[1][axis]
                 + 3*u*t*t * points[2][axis] + t**3 * points[3][axis]
                 for axis in range(3))


def _distance_to_ellipse(radial: float, vertical: float,
                         radius: float, height: float) -> float:
    """Euclidean distance to the edge bore's horizontal/vertical ellipse."""
    radial = abs(radial)
    vertical = abs(vertical)
    if (radial / radius) ** 2 + (vertical / height) ** 2 <= 1.0:
        return 0.0

    lo, hi = 0.0, math.pi / 2.0
    ratio = (math.sqrt(5.0) - 1.0) / 2.0

    def squared_distance(angle: float) -> float:
        return ((radial - radius * math.cos(angle)) ** 2
                + (vertical - height * math.sin(angle)) ** 2)

    left = hi - ratio * (hi - lo)
    right = lo + ratio * (hi - lo)
    left_value, right_value = squared_distance(left), squared_distance(right)
    for _ in range(48):
        if left_value < right_value:
            hi, right, right_value = right, left, left_value
            left = hi - ratio * (hi - lo)
            left_value = squared_distance(left)
        else:
            lo, left, left_value = left, right, right_value
            right = lo + ratio * (hi - lo)
            right_value = squared_distance(right)
    return math.sqrt(min(left_value, right_value))


def _actual_edge_clearance(passage, edge, scale: float,
                           chamber: dict[str, Any]) -> float:
    """Return minimum true tube clearance outside the landing chamber.

    The route edge is a swept ellipse, with its recorded polyline as the axis.
    Comparing the tunnel tube to that ellipse avoids the false clashes caused by
    intersecting two enclosing boxes.
    """
    controls = [tuple(value / 1.0e6 for value in point)
                for point in edge["control_points_um"]]
    axis = [tuple(value / 1.0e6 for value in point)
            for point in edge.get("polyline_um", [])]
    if len(axis) < 2:
        axis = [_bezier(controls, index / 64.0) for index in range(65)]
    edge_radius = (edge["bore_radius_um"]
                   + edge["irregularity_amplitude_um"]) / 1.0e6
    edge_height = (edge["bore_height_um"]
                   + edge["irregularity_amplitude_um"]) / 1.0e6
    passage_segments = passage.get("_audit_segments") or tuple(_segments(passage, scale))
    minimum_clearance = math.inf
    for a, b, _, passage_radius in passage_segments:
        length = math.dist(a[:3], b[:3])
        steps = max(1, math.ceil(length / 0.05))
        for step in range(steps + 1):
            t = step / steps
            point = tuple(a[i] + (b[i] - a[i]) * t for i in range(3))
            if _inside_chamber(point, chamber, scale):
                continue
            tunnel_radius = (a[3] + (b[3] - a[3]) * t
                             + INFLUENCE_REACH_M)
            for start, end in zip(axis, axis[1:]):
                delta = tuple(end[i] - start[i] for i in range(3))
                length_squared = sum(value * value for value in delta)
                projection = (sum((point[i] - start[i]) * delta[i]
                                  for i in range(3)) / length_squared
                              if length_squared else 0.0)
                projection = min(1.0, max(0.0, projection))
                center = tuple(start[i] + projection * delta[i] for i in range(3))
                radial = math.hypot(point[0] - center[0], point[1] - center[1])
                vertical = abs(point[2] - center[2])
                clearance = (_distance_to_ellipse(radial, vertical,
                                                  edge_radius, edge_height)
                             - tunnel_radius)
                minimum_clearance = min(minimum_clearance, clearance)
    return minimum_clearance


def _test_bounds(bounds, strate_index: int, passages, scale: float,
                 *, allow_chamber: bool = True):
    for passage in passages:
        if strate_index not in (passage["upper_strate_index"], passage["lower_strate_index"]):
            continue
        chamber_key = ("upper_landing_chamber" if strate_index == passage["upper_strate_index"]
                       else "lower_landing_chamber")
        chamber = passage[chamber_key]
        segments = passage.get("_audit_segments")
        if segments is None:
            segments = tuple(_segments(passage, scale))
        for a, b, along, radius in segments:
            expanded = [(a[i] - radius, b[i] + radius) if a[i] <= b[i]
                        else (b[i] - radius, a[i] + radius) for i in range(3)]
            if any(expanded[i][1] < bounds[i][0] or expanded[i][0] > bounds[i][1]
                   for i in range(3)):
                continue
            overlap = _overlap_interval(a, b, bounds, radius)
            if overlap is None:
                continue
            first, last, closest = overlap
            if allow_chamber:
                points = [tuple(a[i] + (b[i] - a[i]) * t for i in range(3))
                          for t in (first, closest, last)]
                if all(_inside_chamber(point, chamber, scale) for point in points):
                    continue
            point = tuple(a[i] + (b[i] - a[i]) * closest for i in range(3))
            return passage["stable_id"], along + closest * math.dist(a[:3], b[:3]), point
    return None


def audit_layout(layout: dict[str, Any]) -> dict[str, Any]:
    scale = layout["core_voxel_size_um"] / (256.0 * 1.0e6)
    passages = [dict(row, _audit_segments=tuple(_segments(row, scale)))
                for row in layout.get("passages", [])]
    primitives = [*layout.get("cavern_rooms", []), *layout.get("cavern_edges", [])]
    checked = conflicts = chamber_joints = bounds_candidates = refined_false_positives = 0
    first_conflict = None
    for primitive in primitives:
        bounds = _primitive_bounds(primitive)
        checked += 1
        primitive_conflicted = False
        for passage in passages:
            found = _test_bounds(bounds, primitive["strate_index"], [passage], scale)
            if not found:
                continue
            bounds_candidates += 1
            passage_id, station, point = found
            if "control_points_um" in primitive:
                side = ("upper_landing_chamber"
                        if primitive["strate_index"] == passage["upper_strate_index"]
                        else "lower_landing_chamber")
                clearance = _actual_edge_clearance(passage, primitive, scale,
                                                   passage[side])
                if clearance >= 0.0:
                    refined_false_positives += 1
                    continue
            primitive_conflicted = True
            conflicts += 1
            if first_conflict is None:
                first_conflict = {"primitive": primitive["stable_id"],
                                  "passage": passage_id, "station_m": round(station, 3),
                                  "center_m": [round(value, 3) for value in point]}
            break
        if not primitive_conflicted:
            # Count the edge's named landing association as a confirmed joint.
            anchor = primitive.get("start_anchor_stable_id", "")
            if anchor.endswith(("/upper-landing-chamber", "/lower-landing-chamber")):
                chamber_joints += 1

    negative_controls = 0
    long_passage = max(passages, key=lambda row: len(row["control_points"]), default=None)
    if long_passage:
        controls = [_point(row, scale) for row in long_passage["control_points"]]
        index = min(max(1, len(controls) // 2), len(controls) - 2)
        center = controls[index]
        bounds = tuple((center[axis] - 0.25, center[axis] + 0.25) for axis in range(3))
        # Select the side whose chamber is farther from this forced tunnel overlap.
        strate = long_passage["lower_strate_index"]
        if _test_bounds(bounds, strate, [long_passage], scale, allow_chamber=False):
            negative_controls = 1
    return {"primitives_checked": checked, "bounds_candidates": bounds_candidates,
            "refined_edge_false_positives": refined_false_positives,
            "outside_chamber_conflicts": conflicts,
            "chamber_joints": chamber_joints, "negative_controls": negative_controls,
            "negative_controls_expected": 1, "first_conflict": first_conflict,
            "passed": conflicts == 0 and negative_controls == 1}