Aislo/resources/tester/test_536_2_e07_area.py

"""536-2 — E07 광역 자연사면 위험 `E07_StabArea_Engine` · `_Router`.

무한사면 해석해(건조 · 포화 · 점착) · SINMAP 안정 지수(닫힌 꼴 = 표본 · 등급 경계) ·
D8 집수면적 견본(기운 평면 · V 골짜기) · 습윤지수 · SI 하한 / 상한 해석해 · 시나리오 둘(표준 · 최악) ·
가짜 복사본으로 띠 · 측점표 · 저장본 · 지문.
"""

import base64
import json
import math

import numpy as np
import pytest

from E07_StabArea import E07_StabArea_Engine as engine
from E07_StabArea import E07_StabArea_Router as router_mod

SOIL = dict(gamma_t=18.0, gamma_sat=20.0, c=0.0, c_hi=0.0, phi=35.0, phi_hi=35.0)
SAT = dict(scenario="saturated")


def _fs(beta_deg: float, m: float, **kw) -> float:
    p = engine.AreaParams(**{**SOIL, **kw})
    return float(engine.infinite_slope_fs(np.array([math.radians(beta_deg)]), np.array([m]), p)[0])


def test_무한사면_해석해_건조_포화_점착():
    t = math.tan(math.radians(35)) / math.tan(math.radians(30))
    assert _fs(30, 0.0) == pytest.approx(t, rel=1e-12)  # 건조 c = 0 → tanφ / tanβ
    assert _fs(30, 1.0) == pytest.approx((20 - 9.81) / 20 * t, rel=1e-12)  # 포화 → γ' / γsat 배
    # 손계산 — c 5 · 뿌리 2 · z 1.5 · β 35° · m 0.5
    b, z = math.radians(35), 1.5
    w = 18 * z * 0.5 + 20 * z * 0.5
    n = (18 * z * 0.5 + (20 - 9.81) * z * 0.5) * math.cos(b) ** 2
    hand = (7 + n * math.tan(math.radians(35))) / (w * math.sin(b) * math.cos(b))
    assert _fs(35, 0.5, c=5.0, c_root=2.0) == pytest.approx(hand, rel=1e-12)
    assert math.isinf(_fs(0, 1.0))  # 평평 = 활동 없음


def test_SINMAP_안정지수_닫힌꼴_표본():
    p = engine.AreaParams(
        gamma_t=17.5, gamma_sat=18.5, c=0.0, c_hi=8.0, phi=25.0, phi_hi=38.0, **SAT
    )
    beta = np.radians(np.array([10.0, 30.0, 35.0, 40.0, 60.0]))
    m = np.ones_like(beta)
    si = engine.stability_index(beta, m, p)
    a, b = engine._terms(beta, m, p)
    rng = np.random.default_rng(1)
    c = rng.uniform(0, 8, 400_000)
    t = rng.uniform(math.tan(math.radians(25)), math.tan(math.radians(38)), 400_000)
    for k in range(len(beta)):
        lo = 0 * a[k] + math.tan(math.radians(25)) * b[k]
        fs = c * a[k] + t * b[k]
        if lo >= 1:
            assert si[k] == pytest.approx(lo)  # Fs_min ≥ 1 → SI = Fs_min
        elif fs.max() < 1:
            assert si[k] == 0.0
        else:
            assert si[k] == pytest.approx((fs > 1).mean(), abs=3e-3)  # P(Fs > 1)
    assert si[0] >= 1 and 0 < si[1] < 1 and si[4] == 0  # 세 갈래 모두 지남
    # 폭 0(범위 없음) — 결정 Fs 와 같거나 0
    q = engine.AreaParams(**SOIL, **SAT)
    fs = engine.infinite_slope_fs(beta, m, q)
    si0 = engine.stability_index(beta, m, q)
    assert np.allclose(si0, np.where(fs >= 1, fs, 0.0))


def test_SINMAP_등급_경계():
    si = np.array([np.inf, 1.5, 1.4999, 1.25, 1.2499, 1.0, 0.9999, 0.5, 0.4999, 1e-9, 0.0, np.nan])
    assert engine.classify(si).tolist() == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 0]


def test_D8_집수면적_기운평면_V골짜기():
    cell = 2.0
    ny, nx = 5, 7
    xx = np.arange(nx)[None, :] * cell
    plane = np.repeat(100.0 - 0.5 * xx, ny, axis=0)  # 동쪽으로 기움 → 행마다 동쪽으로 흐름
    a = engine.catchment_area(plane, cell)
    assert np.allclose(a, (np.arange(nx)[None, :] + 1) * cell)
    yy = np.arange(9)[:, None] * cell
    valley = 100.0 + np.abs(yy - 8.0) * 1.0 - 0.1 * np.arange(nx)[None, :] * cell
    a = engine.catchment_area(valley, cell)
    assert a[4, -1] == pytest.approx(9 * nx * cell)  # 전부가 골 출구로
    assert np.all(a[4, :-1] < a[4, 1:])  # 골 따라 늘어남


def test_습윤지수():
    beta = np.radians(np.array([30.0, 30.0, 0.0]))
    a = np.array([100.0, 1e6, 50.0])
    w = engine.wetness(beta, a, 1 / 2000)
    assert w[0] == pytest.approx(100 / (2000 * 0.5))  # R a / (T sinβ)
    assert w[1] == 1.0 and w[2] == 1.0  # 넘침 · 평평 = 포화
    p = engine.AreaParams(**SOIL, tr_lo=2000, tr_hi=3000)
    assert engine.x_range(p) == (1 / 3000, 1 / 2000)
    assert engine.wet_m(beta, a, p).tolist() == w.tolist()  # 기본 = 가장 젖은 쪽
    assert engine.wet_m(beta, a, engine.AreaParams(**SOIL, **SAT)).tolist() == [1.0] * 3


SINMAP = dict(gamma_t=18.0, gamma_sat=18.0, c=2.0, c_hi=6.0, phi=30.0, phi_hi=40.0)


def _sinmap_fs(beta, a, c, tanphi, x, z=1.5, g=18.0):
    """Pack 1998 식 — FS = [C + cosθ(1 − W r) tanφ] / sinθ · C = c / (h ρs g) · h = z cosθ."""
    w = min(x * a / math.sin(beta), 1.0)
    h = z * math.cos(beta)
    return (c / (h * g) + math.cos(beta) * (1 - w * 9.81 / g) * tanphi) / math.sin(beta)


def test_SI_하한_상한_해석해():
    """γt = γsat 이면 엔진 = SINMAP 원식 · 하한 = (c, φ 아래 · T/R 아래) · 상한 = (위 · T/R 위)."""
    p = engine.AreaParams(**SINMAP, tr_lo=1000, tr_hi=3000)
    t1, t2 = math.tan(math.radians(30)), math.tan(math.radians(40))
    beta = np.radians(np.array([15.0, 33.0, 38.0, 45.0, 70.0]))
    area = np.array([20.0, 400.0, 300.0, 800.0, 5000.0])
    si = engine.stability_index(beta, area, p)
    rng = np.random.default_rng(7)
    n = 300_000
    c, t, x = rng.uniform(2, 6, n), rng.uniform(t1, t2, n), rng.uniform(1 / 3000, 1 / 1000, n)
    seen = set()
    for k, (b, ak) in enumerate(zip(beta, area)):
        lo = _sinmap_fs(b, ak, 2.0, t1, 1 / 1000)
        hi = _sinmap_fs(b, ak, 6.0, t2, 1 / 3000)
        fs = engine.infinite_slope_fs(
            beta[k : k + 1], engine.wet_m(beta[k : k + 1], area[k : k + 1], p), p
        )
        assert fs[0] == pytest.approx(lo, rel=1e-12)  # 결정 Fs = 하한
        if lo >= 1:
            assert si[k] == pytest.approx(lo, rel=1e-12)
            seen.add("lo")
        elif hi < 1:
            assert si[k] == 0.0
            seen.add("zero")
        else:
            w = np.minimum(x * ak / math.sin(b), 1)
            h = 1.5 * math.cos(b)
            mc = (c / (h * 18) + math.cos(b) * (1 - w * 9.81 / 18) * t) / math.sin(b)
            assert si[k] == pytest.approx((mc > 1).mean(), abs=3e-3)  # P(Fs > 1) · T/R 도 범위
            seen.add("prob")
    assert seen == {"lo", "zero", "prob"}


def test_시나리오_둘_최악이_더_나쁨():
    beta = np.radians(np.linspace(5, 60, 56))
    area = np.full_like(beta, 50.0)
    std = engine.stability_index(beta, area, engine.AreaParams(**SINMAP))
    worst = engine.stability_index(beta, area, engine.AreaParams(**SINMAP, **SAT))
    assert np.all(worst <= std + 1e-12)  # 완전 포화 = 아래 경계
    assert (engine.classify(worst) >= engine.RISK_CLASS).sum() > (
        engine.classify(std) >= engine.RISK_CLASS
    ).sum()
    # 다 젖는 집수(x·a / sinβ ≥ 1) 면 두 시나리오가 같음
    big = np.full_like(beta, 1e7)
    assert np.allclose(
        engine.stability_index(beta, big, engine.AreaParams(**SINMAP)), worst, rtol=1e-12
    )


def _decode(text: str, dtype: str, grid: dict) -> np.ndarray:
    return np.frombuffer(base64.b64decode(text), dtype=dtype).reshape(grid["ny"], grid["nx"])


@pytest.fixture()
def project(tmp_path, monkeypatch):
    """가짜 E 복사본 — 서쪽이 낮은 30° 비탈(동서 200 m · 남북 100 m) · 노선 = 남북 중심선."""
    x = np.arange(0.0, 200.5, 0.5)
    y = np.arange(0.0, 100.5, 0.5)
    z = np.repeat((x * math.tan(math.radians(30)))[None, :], len(y), axis=0)
    valid = np.ones_like(z, bool)
    valid[:, x > 180] = False  # 동쪽 끝 = 지형 없음
    dem = tmp_path / "dtm.npz"
    np.savez(
        dem,
        x=x.astype(np.float32),
        y=y.astype(np.float32),
        z=z.astype(np.float32),
        valid_mask=valid,
    )
    stations = [
        {"chainage_m": float(c), "label": f"No.{i}", "center_x": 100.0, "center_y": float(c)}
        for i, c in enumerate(range(0, 101, 20))
    ]
    route = [
        {"sequence_num": i, "model_x": 100.0, "model_y": float(c), "chainage_m": float(c)}
        for i, c in enumerate((0, 100))
    ]
    src = router_mod.source
    monkeypatch.setattr(src, "terrain_grid_path", lambda _r: dem)
    monkeypatch.setattr(src, "stations", lambda _r: stations)
    monkeypatch.setattr(src, "db_snapshot", lambda _r: {"route_points": route})
    return tmp_path


def test_길_띠_측점표_저장본(project):
    base = router_mod.defaults(project)
    assert base["scenario"] == "sinmap" and (base["tr_lo"], base["tr_hi"]) == (2000.0, 3000.0)
    settings = {**base, "band_m": 30.0, "cell_m": 2.0, "scenario": "saturated", "c_root": 0.0}
    assert settings["material_id"] == "colluvium" and settings["c_hi"] == settings["c"] + 5
    steps = []
    summary = router_mod.run(project, project, settings, "fp", lambda p, m: steps.append(p))
    doc = router_mod.load(project)
    assert steps == sorted(steps) and doc["fingerprint"] == "fp"
    g = doc["grid"]
    cls = _decode(doc["cls"], "u1", g)
    fs = _decode(doc["fs"], "<u2", g)
    # 붕적토 포화 30° — 손계산 Fs
    b = math.radians(30)
    n = (18.5 - 9.81) * 1.5 * math.cos(b) ** 2
    hand = (5 + n * math.tan(math.radians(27))) / (18.5 * 1.5 * math.sin(b) * math.cos(b))
    inner = fs[cls > 0]
    assert np.median(inner) / 1000 == pytest.approx(hand, abs=2e-3)
    assert summary["fs_min"] == pytest.approx(hand, abs=5e-3)
    # 띠 = 중심선 ± 30 m(격자 상자도 띠만큼) · 등급 면적 합 = 덮인 면적
    assert g["nx"] * g["dx"] <= 62 + 1e-9 and g["ny"] * g["dx"] <= 102 + 1e-9
    assert sum(summary["areas_m2"].values()) == pytest.approx(summary["covered_area_m2"])
    assert summary["covered_area_m2"] == pytest.approx((cls > 0).sum() * 4.0)
    rows = doc["stations"]
    assert [r["chainage_m"] for r in rows] == [0, 20, 40, 60, 80, 100]
    assert all(r["label"] == f"No.{i}" and r["worst_class"] >= 4 for i, r in enumerate(rows))
    assert summary["risk_stations"] == [0, 20, 40, 60, 80, 100]
    assert sum(r["cells"] for r in rows) == (cls > 0).sum()
    # 지문 = 설정에 따라 바뀜
    assert router_mod.fingerprint("s", settings) != router_mod.fingerprint(
        "s", {**settings, "c": 6}
    )
    json.dumps(doc)  # 그대로 JSON