"""노선 추천 후보의 균형 종단(PLAN 63-4 · 63-2 안 A).

후보 노선 하나에 종단을 세워 절 · 성토를 **횡단 ㎥** 로 맞춘 계획선과 그 토량을 낸다.
편집 창 [노선 추천]의 점수 · 미리 보기에만 쓴다 — 파일 · 캐시 · initial_snapshot · 확정 종단을
바꾸지 않는다. 확정 체인의 종단 규칙(1차 계획고 = 지반고 · 2026-09-02 사용자 확정)은 그대로다.

계산 자리 = 서버 단독(③) — 화면은 받은 측점 · 지반고 · 계획고를 그리기만 한다.

방법:
  1. 측점(20 m)마다 지반 횡단(±20 m · 0.5 m)을 DEM 에서 뜬다 — B05 측점 생성과 같은 꼴
     (좌 = +offset · 상단측 = 좌우 평균 표고가 높은 쪽 · 그쪽 절토).
  2. 측점마다 계획고 ±Δ 의 절 · 성 단면적 곡선 — B06 `compute_cross_design` 그대로(지반유형
     리핑암 = B06 기본값 · 곡선부 확폭 = B05 측점 생성과 같은 `_plan_radii` · `_curve_widenings`
     — 빼면 확정 뒤 성토가 1.7배로 벌어짐 · PLAN 77-7) · `AREA_STEP_M` 간격으로 재고 사이는 보간.
  3. 동적계획(앞 두 측점 상태)으로 측점 계획고를 고른다 — 목표 = 절토(다짐 환산) + 성토 ㎥ +
     λ × (성토 − 절토). λ 를 이분해 |절 − 성| / 큰 쪽이 허용(10 %) 안인 가장 작은 λ.
     · 종단기울기 상한 — 넘지 못함.
     · 종단곡선 — 측점마다 변화점(20 m)이라 곡선 반쪽이 직선의 0.45 를 못 넘어 L ≤ 18 m <
       법정 최소 L(20 m 이상). 그래서 대수차가 곡선 생략 한도(비포장 5 %) 안이어야 함 — 넘으면
       벌점(지형상 못 지킬 때만 넘음) · 넘은 변화점 수를 낸다.
  4. 배관 덮개 — 제약은 끈 채(확정 체인과 같음) 관 자리 계획고 − 지반고 < 덮개 여유인 관 수만.
"""

from __future__ import annotations

import time
from pathlib import Path
from typing import TYPE_CHECKING, Any

import numpy as np
from scipy.interpolate import PchipInterpolator, RegularGridInterpolator
from shapely.geometry import LineString, Point

from B05_Profile.B05_Profile_Engine_RouteInitial_Metrics import COST_EARTH_M3
from config.config_system_design import (
    EARTHWORK_CONVERSION_FACTORS,
    FOREST_ROAD_PROFILE_ALIGNMENT,
    FOREST_ROAD_PROFILE_CRITERIA,
)

if TYPE_CHECKING:
    from B05_Profile.B05_Profile_Engine_RouteInitial_Search import Terrain

STATION_M = 20.0  # 측점 간격 — 추천 비교값(`..._Metrics.STATION_M`)과 같음
TAIL_MERGE_M = 5.0  # 종점과 이보다 가까운 정규 측점은 뺌(짧은 끝 구간의 기울기 튐 방지)
CROSS_HALF_M = 20.0  # 지반 횡단 반폭 · 간격 — B06 횡단 파일과 같음
CROSS_STEP_M = 0.5
LEVEL_RANGE_M = 6.0  # 계획고 후보 = 상한 안 기준선 ± 이 폭
LEVEL_STEP_M = 0.25
AREA_STEP_M = 2.0  # 단면적 곡선을 재는 계획고 간격(사이는 PCHIP — 1 m 간격과 3 % 안)
GROUND_TYPE = "ripping_rock"  # B06 기본 지반유형(`attach_default_designs`)
CURVE_PENALTY_M3 = 500.0  # 대수차가 생략 한도를 1 %p 넘을 때 벌점(㎥)
LAMBDA_MAX = 0.95
LAMBDA_STEPS = 14
PIPE_SNAP_M = 5.0  # 관 자리가 노선에서 이 안이면 그 노선의 관으로 봄
GRADE_EPS = 1e-6


def route_sections(
    vertices: list[tuple[float, float]], terrain: Terrain
) -> tuple[np.ndarray, np.ndarray, list[list[dict[str, Any]]], list[str], list[dict[str, Any]]]:
    """측점 누가거리 · 중심 지반고 · 측점별 지반 횡단 샘플 · 단면유형(상단측 절토) · 곡선부 확폭."""
    xy = np.asarray(vertices, dtype=np.float64)
    keep = np.r_[True, np.hypot(*np.diff(xy, axis=0).T) > 1e-9]
    xy = xy[keep]
    seg = np.hypot(*np.diff(xy, axis=0).T)
    chain = np.concatenate([[0.0], np.cumsum(seg)])
    total = float(chain[-1])
    at = np.arange(0.0, total, STATION_M)
    if len(at) > 1 and total - at[-1] < TAIL_MERGE_M:
        at = at[:-1]
    at = np.append(at, total)
    index = np.clip(np.searchsorted(chain, at, side="right") - 1, 0, len(seg) - 1)
    tangent = np.diff(xy, axis=0)[index] / seg[index][:, None]
    left = np.column_stack([-tangent[:, 1], tangent[:, 0]])
    center = np.column_stack([np.interp(at, chain, xy[:, 0]), np.interp(at, chain, xy[:, 1])])
    offsets = np.arange(-CROSS_HALF_M, CROSS_HALF_M + 1e-9, CROSS_STEP_M)
    points = center[:, None, :] + left[:, None, :] * offsets[None, :, None]
    lookup = points.reshape(-1, 2)[:, ::-1]
    grid = (terrain.y, terrain.x)
    z = RegularGridInterpolator(grid, terrain.z, bounds_error=False, fill_value=np.nan)(lookup)
    valid = RegularGridInterpolator(
        grid, terrain.valid.astype(np.float64), method="nearest", bounds_error=False, fill_value=0
    )(lookup)
    z = z.reshape(len(at), len(offsets))
    ok = (valid.reshape(z.shape) > 0.5) & np.isfinite(z)
    middle = int(np.argmin(np.abs(offsets)))
    ground = np.where(ok[:, middle], z[:, middle], np.nan)
    if np.isnan(ground).all():
        raise ValueError("노선 위 지반고가 없어 종단을 세울 수 없습니다.")
    fine = ~np.isnan(ground)
    ground = np.interp(at, at[fine], ground[fine])
    sections, modes = [], []
    for row, flags in zip(z, ok):
        sections.append(
            [
                {"offset_m": float(o), "elevation_m": float(e) if f else None, "valid": bool(f)}
                for o, e, f in zip(offsets, row, flags)
            ]
        )
        lefts, rights = row[(offsets > 0) & flags], row[(offsets < 0) & flags]
        up = "right" if lefts.size and rights.size and rights.mean() > lefts.mean() else "left"
        modes.append(f"{up}_cut")
    return at, ground, sections, modes, _widenings(xy, chain, at)


def _widenings(xy: np.ndarray, chain: np.ndarray, at: np.ndarray) -> list[dict[str, Any]]:
    """측점마다 `compute_cross_design` 확폭 인자 — B05 측점 생성(곡선표 없을 때)과 같은 셈."""
    from B05_Profile.B05_Profile_Engine_Sections_Core import _curve_widenings, _plan_radii

    radii, outer = _plan_radii(xy, chain, at, float(chain[-1]))
    widenings, sides = _curve_widenings(at, radii, outer)
    return [
        {"plan_radius_m": r, "curve_outer_side": s, "curve_widening_m": w}
        for r, s, w in zip(radii, sides, widenings)
    ]


def _area_curves(
    sections: list[list[dict[str, Any]]],
    modes: list[str],
    ground: np.ndarray,
    levels: np.ndarray,
    curves: list[dict[str, Any]] | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """측점 · 계획고 후보마다 (절토 자연 ㎡, 절토 다짐 환산 ㎡, 성토 ㎡)."""
    from B06_Section.B06_Section_Engine_Design import compute_cross_design

    factors = EARTHWORK_CONVERSION_FACTORS
    shape = levels.shape
    cut, cut_c, fill = np.zeros(shape), np.zeros(shape), np.zeros(shape)
    for i, (samples, mode) in enumerate(zip(sections, modes)):
        widen = curves[i] if curves else {}
        delta = levels[i] - ground[i]
        lo = np.floor(delta.min() / AREA_STEP_M) * AREA_STEP_M
        at = np.arange(lo, delta.max() + AREA_STEP_M, AREA_STEP_M)
        rows = []
        for d in at:
            try:
                r = compute_cross_design(
                    samples,
                    float(ground[i] + d),
                    ground_type=GROUND_TYPE,
                    section_mode=mode,
                    **widen,
                )
            except ValueError:  # 지반 샘플이 모자란 측점 — 토량 모름(0)
                rows.append((0.0, 0.0, 0.0))
                continue
            rock = factors.get(r.get("cut_rock_kind") or GROUND_TYPE, factors[GROUND_TYPE])
            rows.append(
                (
                    r["cut_area_m2"],
                    r["cut_soil_area_m2"] * factors["soil"]["compacted"]
                    + r["cut_rock_area_m2"] * rock["compacted"],
                    r["fill_area_m2"],
                )
            )
        table = np.asarray(rows)
        # 단면적은 계획고에 대해 볼록(대략 2차)이라 직선 보간은 크게 잡힌다 — 단조 3차(PCHIP)
        if len(at) > 1:
            curve = PchipInterpolator(at, table, axis=0)(delta)
        else:
            curve = np.repeat(table, len(delta), axis=0)
        cut[i], cut_c[i], fill[i] = curve[:, 0], curve[:, 1], curve[:, 2]
    return cut, cut_c, fill


def _reference(ground: np.ndarray, spans: np.ndarray, max_grade: float) -> np.ndarray:
    """상한 안 기준선 — 앞 · 뒤로 기울기를 누른 두 선의 평균(`..._Metrics.grade_line` 꼴)."""
    forward, backward = ground.copy(), ground.copy()
    for i in range(1, len(ground)):
        rise = max_grade * spans[i - 1]
        forward[i] = min(max(ground[i], forward[i - 1] - rise), forward[i - 1] + rise)
    for i in range(len(ground) - 2, -1, -1):
        rise = max_grade * spans[i]
        backward[i] = min(max(ground[i], backward[i + 1] - rise), backward[i + 1] + rise)
    return 0.5 * (forward + backward)


def _dp(
    levels: np.ndarray, cost: np.ndarray, spans: np.ndarray, max_grade: float, delta_max: float
) -> np.ndarray | None:
    """측점별 계획고 후보 번호 — 기울기 상한은 막고 대수차 넘침은 벌점. 길이 없으면 None."""
    n = len(levels)
    if n == 2:
        steep = np.abs(levels[1][None, :] - levels[0][:, None]) / spans[0] > max_grade + GRADE_EPS
        table = cost[0][:, None] + cost[1][None, :] + np.where(steep, np.inf, 0.0)
        if not np.isfinite(table.min()):
            return None
        return np.array(np.unravel_index(int(table.argmin()), table.shape))
    steep = np.abs(levels[1][None, :] - levels[0][:, None]) / spans[0] > max_grade + GRADE_EPS
    table = cost[0][:, None] + cost[1][None, :] + np.where(steep, np.inf, 0.0)  # [앞앞, 앞]
    back = []
    for i in range(1, n - 1):
        g_in = (levels[i][None, :] - levels[i - 1][:, None]) / spans[i - 1]  # [a, b]
        g_out = (levels[i + 1][None, :] - levels[i][:, None]) / spans[i]  # [b, c]
        over = np.maximum(np.abs(g_out[None, :, :] - g_in[:, :, None]) - delta_max, 0.0)
        total = table[:, :, None] + CURVE_PENALTY_M3 * 100.0 * over  # [a, b, c]
        arg = total.argmin(axis=0)
        table = np.take_along_axis(total, arg[None], axis=0)[0]
        table += np.where(np.abs(g_out) > max_grade + GRADE_EPS, np.inf, 0.0) + cost[i + 1][None, :]
        back.append(arg)
    if not np.isfinite(table.min()):
        return None
    b, c = np.unravel_index(int(table.argmin()), table.shape)
    path = [int(c), int(b)]
    for arg in reversed(back):
        path.append(int(arg[path[-1], path[-2]]))
    return np.array(path[::-1])


def balanced_profile(
    chainage: np.ndarray,
    ground: np.ndarray,
    sections: list[list[dict[str, Any]]],
    modes: list[str],
    *,
    max_grade: float,
    design_speed: int = 20,
    curves: list[dict[str, Any]] | None = None,
) -> dict[str, Any]:
    """상한(종단기울기 · 종단곡선)을 지키며 절 · 성토를 횡단 ㎥ 로 맞춘 계획선과 토량.

    시점 · 종점은 지반고(확정 체인 기본 오프셋 0 과 같음) — 그 길이 없을 때만 풀어 다시 찾는다.
    """
    started = time.perf_counter()
    spans = np.diff(chainage)
    weight = np.zeros(len(chainage))
    weight[:-1] += spans / 2.0
    weight[1:] += spans / 2.0
    offsets = np.arange(-LEVEL_RANGE_M, LEVEL_RANGE_M + 1e-9, LEVEL_STEP_M)
    levels = _reference(ground, spans, max_grade)[:, None] + offsets[None, :]
    cut, cut_c, fill = (
        weight[:, None] * a for a in _area_curves(sections, modes, ground, levels, curves)
    )
    criteria = (
        FOREST_ROAD_PROFILE_CRITERIA["design_speed"].get(design_speed)
        or (FOREST_ROAD_PROFILE_CRITERIA["design_speed"][20])
    )
    skip = FOREST_ROAD_PROFILE_CRITERIA["vertical_curve_skip_delta_pct"] / 100.0
    room = 2.0 * FOREST_ROAD_PROFILE_ALIGNMENT["curve_tangent_max_ratio"] * float(spans.min())
    fits = criteria["min_curve_length_m"] <= room
    delta_max = max(skip, room / criteria["min_vertical_radius_m"]) if fits else skip
    tolerance = FOREST_ROAD_PROFILE_ALIGNMENT["balance_tolerance_percent"] / 100.0
    pinned = np.full(levels.shape, np.inf)
    for end in (0, -1):
        pinned[end, int(np.argmin(np.abs(levels[end] - ground[end])))] = 0.0
    pinned[1:-1] = 0.0

    def solve(lam: float) -> dict[str, Any] | None:
        base = cut_c + fill + lam * (fill - cut_c)
        path = _dp(levels, base + pinned, spans, max_grade, delta_max)
        if path is None:
            path = _dp(levels, base, spans, max_grade, delta_max)
        if path is None:
            return None
        pick = np.arange(len(path))
        c, f = float(cut_c[pick, path].sum()), float(fill[pick, path].sum())
        return {
            "lambda": lam,
            "path": path,
            "cut": float(cut[pick, path].sum()),
            "cut_c": c,
            "fill": f,
            "imbalance": (f - c) / max(c, f, 1e-9),
        }

    tried = [solve(0.0)]
    if tried[0] is None:
        raise ValueError("종단기울기 상한 안에서 계획고를 고를 수 없습니다.")
    first = tried[0]
    if abs(first["imbalance"]) > tolerance:
        sign = 1.0 if first["imbalance"] > 0 else -1.0
        low, high = 0.0, LAMBDA_MAX
        for _ in range(LAMBDA_STEPS):
            mid = 0.5 * (low + high)
            result = solve(sign * mid)
            if result is None:
                break
            tried.append(result)
            if abs(result["imbalance"]) <= tolerance or result["imbalance"] * sign < 0:
                high = mid
            else:
                low = mid
    within = [r for r in tried if abs(r["imbalance"]) <= tolerance]
    best = (
        min(within, key=lambda r: abs(r["lambda"]))
        if within
        else min(tried, key=lambda r: abs(r["imbalance"]))
    )
    plan = levels[np.arange(len(chainage)), best["path"]]
    grades = np.diff(plan) / spans
    changes = np.abs(np.diff(grades))
    return {
        "stations_m": [round(float(v), 2) for v in chainage],
        "ground_m": [round(float(v), 3) for v in ground],
        "plan_m": [round(float(v), 3) for v in plan],
        "cut_m3": round(best["cut"], 1),  # 자연 상태
        "cut_compacted_m3": round(best["cut_c"], 1),
        "fill_m3": round(best["fill"], 1),
        "borrow_m3": round(max(best["fill"] - best["cut_c"], 0.0), 1),
        "spoil_m3": round(max(best["cut_c"] - best["fill"], 0.0), 1),
        "imbalance_pct": round(best["imbalance"] * 100.0, 1),
        "balanced": abs(best["imbalance"]) <= tolerance,
        "tolerance_pct": round(tolerance * 100.0, 1),
        "max_grade_pct": round(float(np.abs(grades).max()) * 100.0 if grades.size else 0.0, 2),
        "grade_limit_pct": round(max_grade * 100.0, 2),
        "curve_over": int((changes > delta_max + GRADE_EPS).sum()),
        "curve_delta_max_pct": round(delta_max * 100.0, 2),
        "max_cut_depth_m": round(float(np.max(ground - plan)), 2),
        "max_fill_height_m": round(float(np.max(plan - ground)), 2),
        "seconds": round(time.perf_counter() - started, 2),
    }


def route_profile(
    vertices: list[tuple[float, float]], terrain: Terrain, *, max_grade: float, design_speed: int
) -> dict[str, Any]:
    """후보 노선 하나의 균형 종단 — 측점 횡단을 뜨고 `balanced_profile`."""
    chainage, ground, sections, modes, curves = route_sections(vertices, terrain)
    if len(chainage) < 2:
        raise ValueError("노선이 짧아 종단을 세울 수 없습니다.")
    return balanced_profile(
        chainage,
        ground,
        sections,
        modes,
        max_grade=max_grade,
        design_speed=design_speed,
        curves=curves,
    )


def cover_shortfall(
    profile: dict[str, Any],
    vertices: list[tuple[float, float]],
    pipes: list[tuple[float, float, float]],
) -> dict[str, int]:
    """관(x, y, 덮개 여유 m) 중 이 노선 위(`PIPE_SNAP_M` 안)인 것 · 그중 덮개가 모자란 수.

    덮개 = 관 자리 계획고 − 지반고(측점 사이 직선 보간). 제약으로 쓰지 않고 세기만 한다."""
    line = LineString(vertices)
    stations = np.asarray(profile["stations_m"])
    rise = np.asarray(profile["plan_m"]) - np.asarray(profile["ground_m"])
    on, short = 0, 0
    for x, y, clearance in pipes:
        point = Point(x, y)
        if line.distance(point) > PIPE_SNAP_M:
            continue
        on += 1
        if float(np.interp(line.project(point), stations, rise)) < clearance:
            short += 1
    return {"pipes_on_route": on, "cover_short": short}


def project_pipes(project_root: Path) -> list[tuple[float, float, float]]:
    """사업지에 놓인 관(x, y, 덮개 여유 m) — 배수유역 편집분을 읽기만 한다(좌표 없는 관은 뺌)."""
    from common_util.common_util_drainage_pipes import facility_clearance_m, read_pipe_points_file
    from config.config_system import (
        DRAINAGE_CACHE_DIRNAME,
        DRAINAGE_EDITS_DIRNAME,
        DRAINAGE_PIPE_POINTS_FILENAME,
    )

    path = (
        project_root
        / "B04_PreProcess"
        / DRAINAGE_CACHE_DIRNAME
        / DRAINAGE_EDITS_DIRNAME
        / DRAINAGE_PIPE_POINTS_FILENAME
    )
    return [
        (float(p.x), float(p.y), facility_clearance_m(p.facility, p.options))
        for p in read_pipe_points_file(path)
        if p.x is not None and p.y is not None
    ]


def attach_profile(
    metrics: dict[str, Any],
    vertices: list[tuple[float, float]],
    terrain: Terrain,
    *,
    max_grade: float,
    design_speed: int,
) -> None:
    """후보 비교값에 균형 종단을 붙이고 토공 점수를 그것으로 바꾼다(한 번만 · 묶음 캐시에 남음).

    점수 = 절토(자연) + 성토 + 반입 + 사토 ㎥ — 가늠값(폭 4 m 띠 · 비탈면 뺌)은
    `earth_estimate_m3` 로 남기고, 공사비 어림의 토공 몫도 이 점수로 바꾼다."""
    if "profile" in metrics:
        return
    try:
        profile = route_profile(vertices, terrain, max_grade=max_grade, design_speed=design_speed)
    except ValueError as error:
        metrics["profile"] = {"error": str(error)}
        return
    estimate = metrics["cut_m3"] + metrics["fill_m3"]
    metrics["cost_index_base"] = metrics["cost_index"]  # 「종단 반영」 끔이 되돌아갈 값
    metrics["profile"] = profile
    metrics["earth_estimate_m3"] = round(estimate, 1)
    metrics["cost_index"] = round(
        metrics["cost_index"] + COST_EARTH_M3 * (profile_earth(profile) - estimate)
    )


def profile_earth(profile: dict[str, Any]) -> float:
    """균형 종단 토공 점수 — 절토(자연) + 성토 + 반입 + 사토 ㎥."""
    return profile["cut_m3"] + profile["fill_m3"] + profile["borrow_m3"] + profile["spoil_m3"]


def plain_metrics(metrics: dict[str, Any]) -> dict[str, Any]:
    """균형 종단을 붙이기 전 비교값(복사본) — 「종단 반영」 끔(원지반 가늠 점수).

    묶음 캐시의 비교값은 종단을 붙인 채 남으므로, 끌 때는 붙인 것을 걷어 낸다."""
    if "profile" not in metrics and "cost_index_base" not in metrics:
        return metrics
    plain = {
        key: value
        for key, value in metrics.items()
        if key not in ("profile", "earth_estimate_m3", "cost_index_base")
    }
    if "cost_index_base" in metrics:
        plain["cost_index"] = metrics["cost_index_base"]
    return plain