From 9d67c3eba6fed5dc8d646f4cb0892d290802dd74 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Tue, 27 Oct 2026 20:25:43 +0900 Subject: [PATCH] =?UTF-8?q?feat(567-2):=20E04=20=EC=9E=84=EA=B3=84=20?= =?UTF-8?q?=ED=99=9C=EB=8F=99=EC=B2=B4=20=EA=B8=B0=EB=91=A5=20=EC=88=98=20?= =?UTF-8?q?=C2=B7=20100=20=EB=AF=B8=EB=A7=8C=20=EA=B2=BD=EA=B3=A0=20=C2=B7?= =?UTF-8?q?=20=EA=B8=B0=EB=B3=B8=20=EC=84=A4=EC=A0=95=20Scoops3D=20?= =?UTF-8?q?=EC=9E=AC=ED=98=84=20=EC=8B=9C=ED=97=98=20=C2=B7=20E07=20?= =?UTF-8?q?=EC=A7=91=EC=88=98=EB=A9=B4=EC=A0=81=20D=E2=88=9E?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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 --- E04_Stab3D/E04_Stab3D_Engine_Search.py | 12 +++- E07_StabArea/E07_StabArea_Engine.py | 68 ++++++++++++++++++++--- resources/tester/test_533_3_e04_engine.py | 35 ++++++++++++ resources/tester/test_536_2_e07_area.py | 18 +++++- 4 files changed, 124 insertions(+), 9 deletions(-) diff --git a/E04_Stab3D/E04_Stab3D_Engine_Search.py b/E04_Stab3D/E04_Stab3D_Engine_Search.py index c62c0fd7f..46d9863d6 100644 --- a/E04_Stab3D/E04_Stab3D_Engine_Search.py +++ b/E04_Stab3D/E04_Stab3D_Engine_Search.py @@ -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"), + }, ) diff --git a/E07_StabArea/E07_StabArea_Engine.py b/E07_StabArea/E07_StabArea_Engine.py index 1ff42172b..e6ce7ee40 100644 --- a/E07_StabArea/E07_StabArea_Engine.py +++ b/E07_StabArea/E07_StabArea_Engine.py @@ -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) diff --git a/resources/tester/test_533_3_e04_engine.py b/resources/tester/test_533_3_e04_engine.py index 25d77abdc..0321582d6 100644 --- a/resources/tester/test_533_3_e04_engine.py +++ b/resources/tester/test_533_3_e04_engine.py @@ -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 diff --git a/resources/tester/test_536_2_e07_area.py b/resources/tester/test_536_2_e07_area.py index 6bfe08965..a19e9a59e 100644 --- a/resources/tester/test_536_2_e07_area.py +++ b/resources/tester/test_536_2_e07_area.py @@ -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])