Aislo/resources/tester/test_364_5_crowded_curves.py

"""364-5 한 측점을 다투는 가까운 두 곡선 — `B05_Profile_Engine_RouteInitial_Station._crowded`.

    ./venv/Scripts/python.exe -m pytest resources/tester/test_364_5_crowded_curves.py -q -s

자료 = 브레인 검증 판(요구 239) 다섯 갈래의 맞추기 전 선(`fixtures/b05_364_5_unsnapped_req239.json`)
+ 검증_sub1 지형. 측점 밖 곡선 수 · 넘은 곳 · 최대 기울기를 맞추기 전 · 옛 맞춤 · 새 맞춤으로 찍고
이웃 MC 20 m 안 측점 밖 쌍 0 · 측점 밖 합계가 옛 맞춤의 1/4 이하 · 이탈 한도 안을 봄.
"""

import json
from pathlib import Path

import numpy as np
import pytest
from scipy.interpolate import RegularGridInterpolator
from shapely.geometry import LineString

from B05_Profile.B05_Profile_Engine_RouteInitial import Terrain
from B05_Profile.B05_Profile_Engine_RouteInitial_Metrics import (
    SAMPLE_M,
    _resample,
    route_metrics,
    smoothed_ground,
)
from B05_Profile.B05_Profile_Engine_RouteInitial_Station import (
    MC_TOL_M,
    STATION_M,
    curve_mc_chainages,
    off_station_count,
    snap_curves_to_stations,
)
from common_util.common_util_route_geometry import read_planned_route_csv
from common_util.common_util_route_polyline import build_planned_polyline
from config.config_system_design import ROUTE_STATION_SNAP_MAX_GAP_M

HERE = Path(__file__).resolve().parent
PROJECT = HERE.parents[1] / "storage/1/3/7e420ed6-f940-4d12-a1a5-39cdcf70a391"
SURFACE = PROJECT / "B05_Profile/route/cost_surface_classification_dtm_smooth.npz"
CASES = json.loads((HERE / "fixtures/b05_364_5_unsnapped_req239.json").read_text(encoding="utf-8"))


def _ls(vertices: list) -> LineString:
    return LineString([xy for i, xy in enumerate(vertices) if not i or xy != vertices[i - 1]])


def _outline(case: dict):
    nodes = case["nodes"]
    return build_planned_polyline(
        [(n["x"], n["y"]) for n in nodes],
        min_radius_m=case["criteria"]["min_radius_m"],
        simplify=False,
        curve_flags=[n["radius_m"] is not None for n in nodes],
        radii=[n["radius_m"] for n in nodes],
    )


def _excess(terrain: Terrain, max_grade: float, smooth: float):
    """`generate_initial_route.grade_excess` 와 같은 잣대(20 m 고른 지반 기울기 − 상한)."""
    surface = RegularGridInterpolator(
        (terrain.y, terrain.x),
        smoothed_ground(terrain, smooth),
        bounds_error=False,
        fill_value=None,
    )

    def excess(points):
        xy, chain = _resample(points, SAMPLE_M)
        z = surface(xy[:, ::-1])
        span = min(int(round(STATION_M / SAMPLE_M)), len(z) - 1)
        if span < 1:
            return 0.0
        run = np.maximum(chain[span:] - chain[:-span], 1e-9)
        return max(0.0, float((np.abs(z[span:] - z[:-span]) / run).max()) - max_grade)

    return excess


def _off_mcs(outline) -> list[float]:
    return [
        mc
        for mc in curve_mc_chainages(outline)
        if abs(mc - STATION_M * round(mc / STATION_M)) > MC_TOL_M
    ]


def test_req239_five_modes_before_after() -> None:
    if not SURFACE.is_file():
        pytest.skip("검증_sub1 storage 없음")
    cached = np.load(SURFACE, allow_pickle=False)
    terrain = Terrain(cached["x"], cached["y"], cached["z"], cached["valid_mask"])
    expected = [
        (v.x, v.y)
        for v in read_planned_route_csv(PROJECT / "B05_Profile/route/expected_route.csv").vertices
    ]
    old_total = new_total = 0
    for case in CASES:
        cr = case["criteria"]
        before = _outline(case)
        follow = case["mode"] == "follow"
        excess = (
            None if follow else _excess(terrain, cr["max_grade_pct"] / 100, cr["grade_smooth_m"])
        )
        after = snap_curves_to_stations(before, cr["min_radius_m"], excess=excess)
        assert after is not None

        def grade(outline):
            m = route_metrics(
                outline.vertices,
                expected,
                terrain,
                [],
                outline.curves,
                max_grade=cr["max_grade_pct"] / 100,
                min_radius_m=cr["min_radius_m"],
                stream_offset_m=cr["stream_offset_m"],
                road_width_m=cr["road_width_m"],
                bound_pct=None,
                grade_smooth_m=cr["grade_smooth_m"],
            )
            return m["over_limit"]["total"], m["max_grade_pct"]

        (o0, g0), (o1, g1) = grade(before), grade(after)
        off = off_station_count(after)
        print(
            f"{case['mode']}: 측점 밖 {case['old']['off']}/{case['old']['curves']}(옛)"
            f" → {off}/{len(after.curves)} · 넘은 곳 {o0} → {o1}"
            f" · 최대 기울기 {g0:.1f} → {g1:.1f} %(맞추기 전 → 새)"
        )
        old_total += case["old"]["off"]
        new_total += off
        # 다투는 쌍 없음 — 이웃 측점 밖 곡선 MC 가 20 m 안이면 한 측점을 다툼
        mcs = _off_mcs(after)
        assert all(b - a >= STATION_M for a, b in zip(mcs, mcs[1:])), mcs
        gap = _ls(after.vertices).hausdorff_distance(_ls(before.vertices))
        assert gap <= ROUTE_STATION_SNAP_MAX_GAP_M + 0.05
        assert after.vertices[0] == pytest.approx(before.vertices[0], abs=1e-3)
        assert after.vertices[-1] == pytest.approx(before.vertices[-1], abs=1e-3)
    print(f"측점 밖 곡선 합 {old_total} → {new_total}")
    assert new_total * 4 <= old_total


def test_crowded_pair_drops_one() -> None:
    """11 m 사이 두 곡선 · 이탈 한도 1 m(측점 옮김 막힘) — 옛 맞춤은 둘 다 그 자리(44.8 · 55.8)."""
    points = [(0.0, 0.0), (45.0, 0.0), (54.0, 7.0), (58.0, 70.0)]
    before = build_planned_polyline(
        points, min_radius_m=12.0, simplify=False, curve_flags=[True] * 4
    )
    after = snap_curves_to_stations(before, 12.0, max_gap_m=1.0)
    assert after is not None and len(after.curves) == 1
    assert _ls(after.vertices).hausdorff_distance(_ls(before.vertices)) <= 1.0 + 0.05