Aislo/resources/tester/test_537_1_e03_prob.py

# -*- coding: utf-8 -*-
"""537-1 E03 2D 확률 해석 · 민감도 — 엔진(`E03_StabSection_Engine_Prob`).

- 무한사면 해석해: F = c/(γ z sinβ cosβ) + (1 − m γw/γ) tanφ/tanβ · c · φ 정규(MIT OCW 1.151 App. 17 ·
  z 5 · m 0.5 · β 35° · c N(25, 5) · φ N(30°, 7.5°)) — Pf · μ · σ 정확값 = φ 위 수치 적분 · LHS 1만 ±2 %
  · c · tanφ 정규면 F 도 정규 → Pf = Φ(−β) · 절편 경로(임계면 고정 · 건조 무한사면) = 닫힌 꼴.
- Slide2 Tutorial 08(Probabilistic Analysis · LHS 1000 · Global Minimum · Bishop): FS 평균 1.145 · σ 0.1206
  · PF 10.2 % · RI 정규 1.204 · 대수정규 1.238.
- Slide2 Verification #36(Li & Lumb 1987 · Hassan & Wolff 1999 · r_u): 결정론 1.340(H&W 1.334) · 평균 1.355 ·
  PF 1.5 % · RI 대수정규 2.482(H&W 2.336).
- 고정 시드 재현 · LHS 층 · 절단 · 대수정규 평균 · 민감도(평균 = 결정론 F · 토네이도 순) · 매회 탐색 갈래.

    ./venv/Scripts/python.exe -m pytest resources/tester/test_537_1_e03_prob.py -q
"""

from __future__ import annotations

import math
import sys
from pathlib import Path

import numpy as np
import pytest
from scipy import integrate
from scipy.stats import norm

ROOT = Path(__file__).resolve().parents[2]
sys.path.insert(0, str(ROOT))

from E03_StabSection import E03_StabSection_Engine_Prob as prob  # noqa: E402
from E03_StabSection.E03_StabSection_Engine_Search import analyze  # noqa: E402


def var(name, prop, mean, std, rel=None, mat="s", dist="normal"):
    return {
        "name": name, "material": mat, "prop": prop, "dist": dist,
        "mean": mean, "std": std, "lo": rel, "hi": rel,
    }  # fmt: skip


def near(v, ref, tol):
    assert v is not None and abs(v - ref) <= tol * abs(ref), (v, ref)


# ---- 무한사면 해석해 ---------------------------------------------------------------
Z, BETA, G, GW = 5.0, math.radians(35.0), 20.0, 9.81
C_PHI = [var("c", "c", 25.0, 5.0), var("phi", "phi", 30.0, 7.5)]


def _ab(m):
    return 1 / (G * Z * math.sin(BETA) * math.cos(BETA)), (1 - m * GW / G) / math.tan(BETA)


def _exact(m):
    """φ 위 수치 적분 — Pf = ∫ Φ((1 − b tanφ − a μc)/(a σc)) f(φ) dφ · 평균 · 표준편차."""
    a, b = _ab(m)
    f = lambda p: norm.pdf(p, 30.0, 7.5)  # noqa: E731
    F = lambda p: a * 25.0 + b * math.tan(math.radians(p))  # noqa: E731
    pf = integrate.quad(lambda p: norm.cdf((1 - F(p)) / (a * 5.0)) * f(p), 0, 89, limit=200)[0]
    mu = integrate.quad(lambda p: F(p) * f(p), 0, 89)[0]
    m2 = integrate.quad(lambda p: (F(p) ** 2 + (a * 5.0) ** 2) * f(p), 0, 89)[0]
    return pf, mu, math.sqrt(m2 - mu * mu)


def test_infinite_slope_closed_form_pf_beta():
    """MIT App. 17 무한사면 — LHS 1만 Pf · β ±2 % · 평균 · 표준편차 ±1 %(φ 정규 · 각도)."""
    a, b = _ab(0.5)
    X = prob.sample(C_PHI, 10000, seed=0)
    s = prob.statistics(a * X[:, 0] + b * np.tan(np.radians(X[:, 1])))
    pf, mu, sd = _exact(0.5)
    near(pf, 0.2252, 0.01)  # 같은 문헌 MC 2×10⁶
    near(s["pf"], pf, 0.02)
    near(s["beta_normal"], (mu - 1) / sd, 0.02)
    near(s["mean"], mu, 0.01)
    near(s["std"], sd, 0.01)
    V = sd / mu
    near(
        s["beta_lognormal"],
        math.log(mu / math.sqrt(1 + V * V)) / math.sqrt(math.log1p(V * V)),
        0.02,
    )


def test_infinite_slope_normal_tanphi_exact():
    """c · tanφ 정규 · 독립 → F 정규 → Pf = Φ(−β) 정확(LHS 2만 ±2 %)."""
    a, b = _ab(0.5)
    vs = [var("c", "c", 25.0, 5.0), {**var("t", "c", 0.57735, 0.11547), "prop": "x"}]
    X = prob.sample(vs, 20000, seed=0)
    s = prob.statistics(a * X[:, 0] + b * X[:, 1])
    mu = a * 25.0 + b * 0.57735
    sd = math.hypot(a * 5.0, b * 0.11547)
    near(mu, 1.154412, 1e-5)
    near(sd, 0.163756, 1e-4)
    near(s["beta_normal"], (mu - 1) / sd, 0.02)
    near(s["pf"], norm.cdf(-(mu - 1) / sd), 0.02)


def _infinite_model():
    """건조 무한사면 β 35° · 지표 아래 5 m 나란한 면 · 모형 양 끝까지(끝 효과 없음)."""
    t, L = math.tan(BETA), 100.0
    model = {
        "materials": [{"id": "s", "name": "s", "gamma_t": G, "c": 25.0, "phi": 30.0}],
        "regions": [{"material_id": "s", "polygon": [[0, 0], [L, -t * L], [L, -t * L - 20], [0, -20]]}],
    }  # fmt: skip
    return model, {"kind": "polyline", "points": [[0, -Z], [L, -t * L - Z]]}


@pytest.mark.parametrize("method", ["ordinary", "bishop", "janbu"])
def test_fixed_surface_slices_equal_closed_form(method):
    """임계면 고정 절편 경로 = 무한사면 닫힌 꼴(물성 바꾼 줄마다) · LHS 1만 Pf ±2 %."""
    model, surface = _infinite_model()
    fixed = prob.Fixed(model, surface, 1)
    a, b = _ab(0.0)
    X = np.array([[25.0, 30.0], [10.0, 20.0], [40.0, 35.0]])
    got = fixed.F(C_PHI, X, method)
    assert got == pytest.approx(a * X[:, 0] + b * np.tan(np.radians(X[:, 1])), rel=1e-9)
    out = prob.run(model, surface=surface, direction=1, variables=C_PHI, n=10000, method=method)
    near(out["stats"]["pf"], _exact(0.0)[0], 0.02)


# ---- Slide2 대조 -------------------------------------------------------------------
T08 = {
    "materials": [{"id": "s", "name": "s", "gamma_t": 19.0, "c": 5.0, "phi": 30.0}],
    "regions": [{"material_id": "s", "polygon": [
        [0, 0], [130, 0], [130, 50], [80, 50], [50, 30], [0, 30]]}],
}  # fmt: skip
T08_VARS = [var("c", "c", 5, 1, 3), var("phi", "phi", 30, 3, 9), var("g", "gamma", 19, 0.5, 1.5)]


def _surface(model, method="bishop"):
    return analyze(model, direction=-1, noncircular=False, methods=(method,))["methods"][method][
        "min"
    ]


def test_slide2_tutorial08_global_minimum():
    """Tutorial 08 — LHS 1000 · 고정 시드 셋 평균: FS 평균 ±1 % · σ ±3 % · PF ±2 %p · RI ±3 %."""
    surface = _surface(T08)
    near(surface["F"], 1.142, 0.01)
    runs = [
        prob.run(T08, surface=surface, direction=-1, variables=T08_VARS, n=1000, seed=k)["stats"]
        for k in range(3)
    ]
    avg = {
        k: float(np.mean([r[k] for r in runs]))
        for k in ("mean", "std", "pf", "beta_normal", "beta_lognormal")
    }
    near(avg["mean"], 1.145, 0.01)
    near(avg["std"], 0.1206, 0.03)
    assert abs(avg["pf"] - 0.102) <= 0.02
    near(avg["beta_normal"], 1.204, 0.03)
    near(avg["beta_lognormal"], 1.238, 0.03)


def test_slide2_36_li_lumb_ru():
    """#36 r_u 확률 변수 · MC 1만 — 결정론 ±1 % · 평균 ±1.5 % · PF ±1 %p · RI 대수정규 = Hassan & Wolff ±3 %."""
    model = {
        "materials": [{"id": "s", "name": "s", "gamma_t": 18.0, "c": 18.0, "phi": 30.0}],
        "regions": [{"material_id": "s", "polygon": [[0, 0], [20, 0], [20, 15], [15, 15], [5, 5], [0, 5]]}],
        "groundwater": {"mode": "ru", "ru": 0.2},
    }  # fmt: skip
    surface = _surface(model)
    near(surface["F"], 1.340, 0.01)
    vs = [
        var("c", "c", 18, 3.6), var("phi", "phi", 30, 3), var("g", "gamma", 18, 0.9),
        var("ru", "ru", 0.2, 0.02, mat=None),
    ]  # fmt: skip
    s = prob.run(model, surface=surface, direction=-1, variables=vs, n=10000, sampling="mc")[
        "stats"
    ]
    near(s["mean"], 1.355, 0.015)
    assert abs(s["pf"] - 0.015) <= 0.01
    near(s["beta_lognormal"], 2.336, 0.03)


# ---- 표본 · 재현 · 민감도 · 갈래 ---------------------------------------------------------
def test_seed_reproducible_and_lhs_strata():
    """같은 시드 = 같은 표본 · 다른 시드 = 다름 · LHS = 변수마다 n 칸에 하나씩 · 절단 안 · 상관."""
    a = prob.sample(T08_VARS, 500, seed=7)
    assert np.array_equal(a, prob.sample(T08_VARS, 500, seed=7))
    assert not np.array_equal(a, prob.sample(T08_VARS, 500, seed=8))
    for j, v in enumerate(T08_VARS):
        lo, hi = v["mean"] - v["lo"], v["mean"] + v["hi"]
        assert a[:, j].min() >= lo and a[:, j].max() <= hi
        pa, pb = norm.cdf(-3), norm.cdf(3)
        u = (norm.cdf((a[:, j] - v["mean"]) / v["std"]) - pa) / (pb - pa)
        assert sorted(np.floor(u * 500).astype(int)) == list(range(500))
    rho = prob.sample(C_PHI, 5000, seed=1, correlations=[(0, 1, -0.5)])
    assert np.corrcoef(rho.T)[0, 1] == pytest.approx(-0.5, abs=0.03)


def test_truncation_and_lognormal_moments():
    """대수정규 = 평균 · σ 그대로 · 음수 없음 · 정규 c 는 물리 한계 0 에서 절단."""
    ln = prob.sample([var("c", "c", 10.0, 5.0, dist="lognormal")], 20000, seed=2)[:, 0]
    near(ln.mean(), 10.0, 0.01)
    near(ln.std(), 5.0, 0.03)
    assert ln.min() > 0
    nc = prob.sample([var("c", "c", 5.0, 5.0)], 5000, seed=2)[:, 0]
    assert nc.min() >= 0


def test_sensitivity_and_tornado():
    """민감도 — 평균 줄 F = 결정론 F · 곡선 51 점 · c ↑ = F ↑ · 토네이도 = 폭 큰 순."""
    surface = _surface(T08)
    out = prob.run(T08, surface=surface, direction=-1, variables=T08_VARS, n=200)
    sens = out["sensitivity"]
    assert sens["F_mean"] == pytest.approx(surface["F"], rel=1e-9)
    curves = {c["name"]: c for c in sens["curves"]}
    assert len(curves["c"]["F"]) == 51 and curves["c"]["x"][0] == pytest.approx(2.0)
    assert np.all(np.diff(curves["c"]["F"]) > 0) and np.all(np.diff(curves["phi"]["F"]) > 0)
    swings = [t["swing"] for t in sens["tornado"]]
    assert swings == sorted(swings, reverse=True) and sens["tornado"][0]["name"] in ("c", "phi")
    assert out["stats"]["n"] == 200 and len(out["F"]) == 200


def test_search_mode_matches_fixed_at_low_variance():
    """매회 탐색 — 변동이 아주 작으면 고정면과 같은 F(±1 %) · 표본 수 · 시간 칸."""
    surface = _surface(T08)
    tiny = [var("c", "c", 5, 0.01), var("phi", "phi", 30, 0.01)]
    out = prob.run(T08, surface=surface, direction=-1, variables=tiny, n=4, mode="search")
    for f in out["F"]:
        near(f, surface["F"], 0.01)
    assert out["mode"] == "search" and out["elapsed_s"]["samples"] > 0


if __name__ == "__main__":
    sys.exit(pytest.main([__file__, "-q"]))