feat(567-2): E04 임계 활동체 기둥 수 · 100 미만 경고 · 기본 설정 Scoops3D 재현 시험 · E07 집수면적 D∞

- E04 best · stats 에 columns · few_columns(< 100 · Scoops3D limcol) - 시험 533_3: 기본 설정 그대로 · 상자 없이 제방 30° 1.3203 · D&G 1a 1.0385 · 제방 60° 0.7488(±2 %) · 기둥 경고 - E07 집수면적 = D∞(Tarboton 1997 · SINMAP 원본) · 내리막 없는 채운 평지만 D8 - 시험 536_2: D∞ 면 몫(22.5° = 반씩) · 30° 돌린 평면 줄무늬 표준편차 0.05(D8 0.16) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-27 20:25:43 +09:00
parent a843aa3ad3
commit 9d67c3eba6
4 changed files with 124 additions and 9 deletions
@@ -40,6 +40,8 @@ from scipy import ndimage
 from E04_Stab3D.E04_Stab3D_Engine_Column import Columns, Strata, solve

 TRUNCATION_JANBU = 0.10
+#: 임계 활동체 기둥(셀) 수 경고 문턱 — Scoops3D limcol 기본(≥ 200 개면 F 1 % 안 · Reid 외 2015 §5.3.4)
+MIN_COLUMNS = 100
 _STRUCT = np.zeros((3, 3, 3), bool)
 _STRUCT[1] = True  # 같은 후보 안에서만 8 이웃 · 후보끼리는 안 이음
 BATCH_CELLS = 3_000_000  # 한 번에 만드는 (후보 × 창) 칸 수 — 메모리 약 0.5 GB 아래
@@ -541,6 +543,9 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
        if pool:  # 끊기(InterruptedError) · 오류에도 일꾼을 닫음
            pool.shutdown(cancel_futures=True)
    f_out = np.where(np.isfinite(f_min), f_min, np.nan)
+    if "cells" in best:  # 기둥 수 · 경고(셀이 굵으면 F 오차가 큼 — 셀 줄임 권함)
+        best["columns"] = len(best["cells"])
+        best["few_columns"] = best["columns"] < MIN_COLUMNS
    return SearchResult(
        f_min=f_out,
        volume=volume,
@@ -549,7 +554,12 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
        trials=counters["trials"],
        evaluated=counters["evaluated"],
        seconds=time.perf_counter() - t0,
-        stats={"nodes": len(seen), "cells": int(np.isfinite(f_out).sum())},
+        stats={
+            "nodes": len(seen),
+            "cells": int(np.isfinite(f_out).sum()),
+            "columns": best.get("columns"),
+            "few_columns": best.get("few_columns"),
+        },
    )


@@ -3,8 +3,9 @@
 순수 계산(storage 모름) — D 예정노선 재해 안전성도 같은 엔진을 씀(한 벌).
 근거 = 지식DB `original/예정노선/02` §1.1 · §2 ① · `original/안정성평가/04`(붕적토 범위).

-  지형   DEM → 고른 격자(cell_m) · 경사 β(Horn 3×3) · D8 흐름(웅덩이 채움 · 수신 셀 = B04 배수유역
-         함수 재사용) · 단위 집수면적 a = 상류 셀 수 × cell²/cell(등고선 폭 = cell)
+  지형   DEM → 고른 격자(cell_m) · 경사 β(Horn 3×3) · D∞ 흐름(Tarboton 1997 · SINMAP 원본 · 웅덩이 채움
+         = B04 배수유역 함수 재사용 · 내리막 없는 평지만 D8) · 단위 집수면적 a = 상류 셀 몫 합 × cell²/cell
+         (등고선 폭 = cell)
  포화비 m  시나리오 둘
         sinmap(기본 · SINMAP 표준) = 습윤지수 W = min(x·a / sinβ, 1) · x = R/T ∈ [1/(T/R)_hi, 1/(T/R)_lo]
         saturated(최악 · 강우 극한) = 지표 포화 m = 1
@@ -98,8 +99,54 @@ def slope_rad(z: np.ndarray, cell_m: float) -> np.ndarray:
    return np.arctan(np.hypot(gx, gy))


+#: D∞ 삼각 면 여덟(Tarboton 1997 표 1) — (e1 = 곁 이웃 행·열, e2 = 대각 이웃 행·열)
+_FACETS = (
+    ((0, 1), (-1, 1)),
+    ((-1, 0), (-1, 1)),
+    ((-1, 0), (-1, -1)),
+    ((0, -1), (-1, -1)),
+    ((0, -1), (1, -1)),
+    ((1, 0), (1, -1)),
+    ((1, 0), (1, 1)),
+    ((0, 1), (1, 1)),
+)
+
+
+def dinf_receivers(z: np.ndarray, cell_m: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """D∞ — 셀마다 가장 가파른 면의 받는 셀 둘(평면 번호 · 없으면 −1) · e1 몫 = 1 − r/(π/4).
+
+    면 기울기 s1 = (e0 − e1)/d · s2 = (e1 − e2)/d · r = atan(s2/s1) 를 [0, π/4] 로 자름(끝이면 한 셀로).
+    NaN(밖) 이웃이 낀 면은 안 봄."""
+    ny, nx = z.shape
+    p = np.pad(z.astype(float), 1, constant_values=np.nan)
+    iy, ix = np.indices((ny, nx))
+    best = np.zeros((ny, nx))
+    r1 = np.full((ny, nx), -1, np.int64)
+    r2 = np.full((ny, nx), -1, np.int64)
+    f1 = np.ones((ny, nx))
+    quarter = math.pi / 4
+    with np.errstate(invalid="ignore"):
+        for (a, b), (c, d) in _FACETS:
+            e1 = p[1 + a : 1 + a + ny, 1 + b : 1 + b + nx]
+            e2 = p[1 + c : 1 + c + ny, 1 + d : 1 + d + nx]
+            s1 = (z - e1) / cell_m
+            s2 = (e1 - e2) / cell_m
+            r = np.arctan2(s2, s1)
+            s = np.where(
+                r < 0,
+                s1,
+                np.where(r > quarter, (z - e2) / (cell_m * math.sqrt(2)), np.hypot(s1, s2)),
+            )
+            up = np.isfinite(s) & (s > best)
+            best[up] = s[up]
+            r1[up] = ((iy + a) * nx + ix + b)[up]
+            r2[up] = ((iy + c) * nx + ix + d)[up]
+            f1[up] = 1.0 - np.clip(r, 0.0, quarter)[up] / quarter
+    return r1, r2, f1
+
+
 def catchment_area(z: np.ndarray, cell_m: float) -> np.ndarray:
-    """D8 단위 집수면적 a(m) = 자기 포함 상류 셀 수 × cell. 밖 = NaN.
+    """D∞ 단위 집수면적 a(m) = 자기 포함 상류 셀 몫 합 × cell. 밖 = NaN · 내리막 없는 셀(채운 평지) = D8.

    ponytail: 지형 격자 끝 밖 상류는 모름(복사 DTM = 노선 띠) — 띠 끝 셀은 a 가 작게 나옴 ·
    넓은 DEM(D 예정노선 등)이면 그대로 맞음.
@@ -107,15 +154,22 @@ def catchment_area(z: np.ndarray, cell_m: float) -> np.ndarray:
    filled, valid = condition_surface(z.astype(np.float32))
    spec = SimpleNamespace(n_rows=z.shape[0], n_cols=z.shape[1], cell_m=cell_m)
    receiver, _ = compute_receivers(spec, filled, valid)
+    r1, r2, f1 = dinf_receivers(np.where(valid, filled, np.nan), cell_m)
+    d8 = r1 < 0
+    r1[d8] = receiver.reshape(z.shape)[d8]
+    r2[d8] = -1
+    f1[d8] = 1.0
    flat = filled.ravel()
    order = np.flatnonzero(valid.ravel())
    order = order[np.argsort(-flat[order], kind="stable")].tolist()  # 높은 데서 낮은 데로
    acc = valid.ravel().astype(float).tolist()
-    recv = receiver.tolist()
+    one, two, share = r1.ravel().tolist(), r2.ravel().tolist(), f1.ravel().tolist()
    for i in order:
-        r = recv[i]
-        if r != i:
-            acc[r] += acc[i]
+        j, f = one[i], share[i]
+        if j != i:
+            acc[j] += acc[i] * f
+        if f < 1.0:
+            acc[two[i]] += acc[i] * (1.0 - f)
    out = np.asarray(acc).reshape(z.shape) * cell_m
    return np.where(valid, out, np.nan)

@@ -217,3 +217,38 @@ def test_ground_order_and_label_map(tmp_path):
    assert ids == ["fill", "weathered_soil", "weathered_rock", "soft_rock"]
    mapped = grid_engine.strata_from_grid(grid, label_map={"soil": "colluvium", "fill": "fill"})
    assert [m["id"] for m in mapped.materials][:2] == ["fill", "colluvium"]
+
+
+def _embankment(run_m, gamma, phi, c, half=20.0):
+    """Scoops3D 검증 제방(H 10 m) — 탐색 상자 없이 격자 전체 · y ±20 m."""
+    x = np.arange(-25, run_m + 35 + 1e-9, 1.0)
+    y = np.arange(-half, half + 1e-9, 1.0)
+    xx, _ = np.meshgrid(x, y)
+    z = np.clip(xx / run_m * 10.0, 0, 10.0)
+    mat = {"id": "m", "name": "m", "gamma_t": gamma, "gamma_sat": gamma, "c": c, "phi": phi}
+    return Strata(x, y, z, np.zeros((0,) + z.shape), [mat])
+
+
+@pytest.mark.parametrize(
+    ("run_m", "gamma", "phi", "c", "expected"),
+    [
+        (10 / math.tan(math.radians(30)), 19, 20, 8.64, 1.32),  # 표 6.8 제방 30°
+        (20.0, 20, 19.6, 3.0, 1.04),  # Donald & Giam 1a
+        (10 / math.tan(math.radians(60)), 19, 20, 8.64, 0.75),  # 표 6.11 제방 60°
+    ],
+)
+def test_default_settings_reproduce_scoops3d(run_m, gamma, phi, c, expected):
+    """567-2 — 실사용 기본 설정 그대로(노드 2 m · 거친 4 배 · 높이 2~30 m · dr 1 m · 부피 10~10,000 m³ ·
+    구 + 타원체) · 좁힌 상자 없이 → 문헌값 ±2 %(실측 +0.0 · −0.1 · −0.2 %)."""
+    best = search(_embankment(run_m, gamma, phi, c), SearchParams(workers=8)).best
+    assert best["governing"] == pytest.approx(expected, rel=0.02)
+    assert best["columns"] == len(best["cells"]) >= 100 and best["few_columns"] is False
+
+
+def test_few_columns_warning():
+    """567-2 — 임계 활동체 기둥 < 100(Scoops3D limcol) → few_columns 경고."""
+    st = _embankment(10 / math.tan(math.radians(60)), 19, 20, 8.64, half=8.0)
+    p = SearchParams(box=(-9, -3, -1, 1), h_min=4, h_max=10, v_min=10, v_max=60, ratios=(1.0,))
+    r = search(st, p)
+    assert r.best["columns"] < 100 and r.best["few_columns"] is True
+    assert r.stats["columns"] == r.best["columns"] and r.stats["few_columns"] is True
@@ -1,7 +1,7 @@
 """536-2 — E07 광역 자연사면 위험 `E07_StabArea_Engine` · `_Router`.

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

@@ -84,6 +84,22 @@ def test_D8_집수면적_기운평면_V골짜기():
    assert np.all(a[4, :-1] < a[4, 1:])  # 골 따라 늘어남


+def test_Dinf_면_몫_돌린평면():
+    """567-2 D∞(Tarboton 1997 · SINMAP 원본) — 면 안 22.5° 흐름 = 곁 · 대각 반씩 · 30° 돌린 평면
+    집수면적 = 흐름길 길이(D8 줄무늬 표준편차 0.16 → D∞ 0.05)."""
+    cell = 2.0
+    ii, jj = np.indices((5, 5))
+    r1, r2, f1 = engine.dinf_receivers(-(jj - math.tan(math.pi / 8) * ii) * cell, cell)
+    assert (r1[2, 2], r2[2, 2]) == (2 * 5 + 3, 1 * 5 + 3) and f1[2, 2] == pytest.approx(0.5)
+    n, th = 60, math.radians(30)
+    i, j = np.indices((n, n))
+    z = 500 - (j * math.cos(th) + i * math.sin(th)) * cell * 0.5
+    length = np.minimum((j + 0.5) / math.cos(th), (i + 0.5) / math.sin(th)) * cell
+    inner = (i > 5) & (j > 5) & (i < n - 5) & (j < n - 5)
+    ratio = engine.catchment_area(z, cell)[inner] / length[inner]
+    assert np.std(ratio) < 0.08 and np.median(np.abs(ratio - 1)) < 0.2
+
+
 def test_습윤지수():
    beta = np.radians(np.array([30.0, 30.0, 0.0]))
    a = np.array([100.0, 1e6, 50.0])