eomsangdonandClaude Opus 5.5
496a32e205
E04 Search 충돌 풂 — 점 단위 덩이 나눔 · 전역 최소 다듬기(_set_best) 유지 + sub7 best.cells 를 _set_best 에 넣어 탐색 · 다듬기 두 길 모두 셀 전달 (요구 351) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
680 lines
29 KiB
Python
680 lines
29 KiB
Python
"""E04 3D 활동면 탐색 — Scoops3D 방식 구 격자 + 타원체 · 셀별 최소 안전율 · 임계 체적.
|
||
|
||
근거: Reid 외(2015) USGS TM 14-A1 3.2(탐색 격자 · 반지름 증분 · 부피/면적 제약 ·
|
||
거친→고운 탐색) — 방식만 · 코드 자체.
|
||
다른 점(정하고 보고):
|
||
- 노드 높이 = 노드 아래 지표 + h — 고저차 큰 산지에서 땅속 · 너무 높은 노드를 덜어 냄.
|
||
- 타원체 = 활동 단면(d−z)은 원 · 회전축 방향 반축만 ratio 배(Hungr 1989 꼴)
|
||
— 수직력이 축을 지나 Bishop 그대로. 방향 = 노드 지표 낙하 방향(5×5 평활).
|
||
구는 활동체 지표 평균 낙하 방향(Scoops3D 2.1).
|
||
- 한 구가 지표를 여러 덩이로 자르면 덩이마다 셈 — 셀 안 점 단위로 가름(_split_pieces ·
|
||
발끝을 살짝 비켜 지나는 구의 사면 몸통 · 발끝 앞 접시가 따로).
|
||
- 전역 최소 면은 마지막에 중심 · 반지름을 눈금 사이로 다듬음(_polish).
|
||
- 격자 가장자리 · DEM 빈칸에 닿은 덩이는 버림(잘린 질량 막음).
|
||
- 후보 덩이는 프로세스로 나눠 셈 · 낸 차례대로 합침 — 일꾼 수와 무관하게 같은 값.
|
||
|
||
다스리는 안전율(governing):
|
||
- 못 뚫는 층(기본 발파암 = 연암)에 잘린 바닥 면적 비 ≥ 10 % → Janbu 보정으로 전환
|
||
· 단 Bishop 이 더 낮으면 Bishop(작은 쪽). Hungr 1989 = 약층 복합면에서 Bishop 도
|
||
쓸 만함 · Zhang 1988 약층 경우 Janbu 보정 1.70 > Bishop 1.62 라 전환만은 덜 보수적.
|
||
- 10 % 문턱 = 기둥 몇 개 닿아 방법이 바뀌는 지도 이음매를 막음(Zhang 약층 잘림 비 33 %).
|
||
- 그 밖 → Hungr 3D Bishop · Bishop 이 수렴 못 하면 그 후보는 버림(Scoops3D 2.3.2).
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import math
|
||
import os
|
||
import pickle
|
||
import tempfile
|
||
import time
|
||
import weakref
|
||
from collections.abc import Callable
|
||
from concurrent.futures import ProcessPoolExecutor
|
||
from dataclasses import dataclass, field
|
||
from types import SimpleNamespace
|
||
|
||
import numpy as np
|
||
from scipy import ndimage
|
||
|
||
from E04_Stab3D.E04_Stab3D_Engine_Column import Columns, Strata, solve
|
||
|
||
TRUNCATION_JANBU = 0.10
|
||
_STRUCT = np.zeros((3, 3, 3), bool)
|
||
_STRUCT[1] = True # 같은 후보 안에서만 8 이웃 · 후보끼리는 안 이음
|
||
BATCH_CELLS = 3_000_000 # 한 번에 만드는 (후보 × 창) 칸 수 — 메모리 약 0.5 GB 아래
|
||
|
||
|
||
@dataclass
|
||
class SearchParams:
|
||
"""탐색 조건 — 길이 m · 부피 m³ · 면적 m²."""
|
||
|
||
box: tuple[float, float, float, float] | None = None # (x0, x1, y0, y1) · None = 격자 전체
|
||
node_mask: np.ndarray | None = (
|
||
None # (ny, nx) — 노드 아래 셀이 False 면 노드 뺌(노선 띠만 등 · Scoops3D 파일 탐색)
|
||
)
|
||
spacing: float = 2.0 # 가장 고운 노드 수평 간격
|
||
h_min: float = 2.0 # 노드 높이(지표 위) 범위 · 간격
|
||
h_max: float = 30.0
|
||
h_step: float = 2.0
|
||
coarse: int = 4 # 거친 격자 배수(Scoops3D multres)
|
||
dr: float = 1.0 # 반지름 증분
|
||
max_depth: float = 15.0 # 반지름 상한 = h + max_depth
|
||
v_min: float = 10.0
|
||
v_max: float = 10_000.0
|
||
a_min: float = 0.0
|
||
a_max: float = math.inf
|
||
ratios: tuple[float, ...] = (1.0, 2.0) # 구(늘 봄) + 타원체 축 방향 배수(구 임계 노드에서)
|
||
sub: int = 2 # 셀 안 표본 sub×sub(부분 기둥)
|
||
workers: int = 0 # 프로세스 수 · 0 = 자동(격자 5만 칸 이상이면 CPU − 2 · 12 까지)
|
||
refine_rounds: int = 4 # 가장 고운 간격에서 더 내려가 보는 횟수 상한
|
||
# 고운 탐색 · 타원체는 F < max(이 값, 전역 최소 × 1.2) 로 고친 노드만(높은 곳 지도는 거친 값)
|
||
refine_below: float = 2.0
|
||
groundwater: dict | None = None
|
||
seismic_kh: float = 0.0
|
||
# 화면 고리(533-3) — 결과 불변 · on_progress(percent, message) · should_stop() 참이면 InterruptedError
|
||
on_progress: Callable[[int, str], None] | None = None
|
||
should_stop: Callable[[], bool] | None = None
|
||
|
||
|
||
@dataclass
|
||
class SearchResult:
|
||
f_min: np.ndarray # 셀별 최소 안전율(다스리는 값) · 닿지 않은 셀 = nan
|
||
volume: np.ndarray # 그 임계 활동체 부피
|
||
method: np.ndarray # 0 = Bishop · 1 = Janbu 보정
|
||
best: dict # 전체 최소 활동체(중심 · 반지름 · 비 · 방향 · 방법별 F · 부피 · 셀)
|
||
trials: int = 0
|
||
evaluated: int = 0
|
||
seconds: float = 0.0
|
||
stats: dict = field(default_factory=dict)
|
||
|
||
|
||
def trial_columns(
|
||
strata: Strata,
|
||
centers: np.ndarray,
|
||
r: float,
|
||
ratio: float = 1.0,
|
||
directions: np.ndarray | None = None,
|
||
v_min: float = 0.0,
|
||
v_max: float = math.inf,
|
||
a_min: float = 0.0,
|
||
a_max: float = math.inf,
|
||
sub: int = 2,
|
||
) -> tuple[Columns, np.ndarray, np.ndarray]:
|
||
"""후보 T 개(같은 r · ratio)의 기둥 — (Columns, 덩이의 후보 번호, 후보별 가장 큰 덩이 부피).
|
||
|
||
부분 기둥 — 셀마다 sub×sub 점에서 활동면 · 지표(쌍선형)를 보고 걸친 비율 · 평균 깊이로
|
||
기둥을 만듦(Scoops3D 2.1 부분 기둥과 같은 뜻 · 가장자리 적분 오차 줄임). sub=1 = 셀 중심.
|
||
"""
|
||
centers = np.atleast_2d(np.asarray(centers, float))
|
||
n_t = len(centers)
|
||
dx = strata.dx
|
||
ny, nx = strata.z_top.shape
|
||
w = int(math.ceil(r * max(ratio, 1.0) / dx)) + 1
|
||
off = np.arange(-w, w + 1)
|
||
ix0 = np.rint((centers[:, 0] - strata.x[0]) / dx).astype(int)
|
||
iy0 = np.rint((centers[:, 1] - strata.y[0]) / dx).astype(int)
|
||
ixg = ix0[:, None, None] + off[None, None, :]
|
||
iyg = iy0[:, None, None] + off[None, :, None]
|
||
inb = (ixg >= 0) & (ixg < nx) & (iyg >= 0) & (iyg < ny)
|
||
nw = 2 * w + 1
|
||
shape = (n_t, nw, nw)
|
||
d = None if ratio == 1.0 else np.asarray(directions, float)
|
||
k2 = 1.0 / (ratio * ratio)
|
||
|
||
def remainder(t, xs, ys):
|
||
"""r² − (중심에서 수평 거리 꼴)² — 0 보다 크면 활동면이 그 점 위를 지남."""
|
||
dxp = xs - centers[t, 0]
|
||
dyp = ys - centers[t, 1]
|
||
if d is None:
|
||
return r * r - dxp * dxp - dyp * dyp
|
||
pp = dxp * d[t, 0] + dyp * d[t, 1]
|
||
qq = -dxp * d[t, 1] + dyp * d[t, 0]
|
||
return r * r - pp * pp - k2 * qq * qq
|
||
|
||
def surface(t, xs, ys):
|
||
"""후보 t 의 활동면 — (rem, z_s, ∂z/∂x, ∂z/∂y). 원소별(납작)."""
|
||
dxp = xs - centers[t, 0]
|
||
dyp = ys - centers[t, 1]
|
||
if d is None:
|
||
rem = r * r - dxp * dxp - dyp * dyp
|
||
root = np.sqrt(np.maximum(rem, 1e-12))
|
||
return rem, centers[t, 2] - root, dxp / root, dyp / root
|
||
dxd, dyd = d[t, 0], d[t, 1]
|
||
pp = dxp * dxd + dyp * dyd
|
||
qq = -dxp * dyd + dyp * dxd # 축 방향 a = (−d_y, d_x)
|
||
rem = r * r - pp * pp - k2 * qq * qq
|
||
root = np.sqrt(np.maximum(rem, 1e-12))
|
||
gx = (pp * dxd - k2 * qq * dyd) / root
|
||
gy = (pp * dyd + k2 * qq * dxd) / root
|
||
return rem, centers[t, 2] - root, gx, gy
|
||
|
||
# 1) 셀 중심으로 거르기 — 걸친 셀과 그 이웃만 셀 안 표본을 봄(창 대부분은 허공 · 땅속 밖)
|
||
# 창 칸 번호 → (후보, 셀) — 큰 창 배열을 펴지 않고 셈
|
||
def cell_of(flat):
|
||
t = flat // (nw * nw)
|
||
return t, iy0[t] + off[(flat // nw) % nw], ix0[t] + off[flat % nw]
|
||
|
||
rem0 = remainder(np.arange(n_t)[:, None, None], strata.x[0] + ixg * dx, strata.y[0] + iyg * dx)
|
||
in0 = np.flatnonzero(inb & (rem0 > 0))
|
||
t0, yy0, xx0 = cell_of(in0)
|
||
zs0 = centers[t0, 2] - np.sqrt(rem0.ravel()[in0])
|
||
zb0 = np.maximum(zs0, strata.hard_top[yy0, xx0])
|
||
act0 = np.zeros(shape, bool)
|
||
act0.ravel()[in0[zb0 < strata.z_top[yy0, xx0]]] = True
|
||
cand = act0 if sub == 1 else ndimage.binary_dilation(act0, structure=_STRUCT) & inb
|
||
ci = np.flatnonzero(cand)
|
||
big = np.zeros(n_t)
|
||
e = np.zeros(0, int)
|
||
empty = (Columns(e, e, e, np.zeros(0), np.zeros(0), np.zeros(0), np.zeros((0, 3))), e, big)
|
||
if not len(ci):
|
||
return empty
|
||
tt, cy, cx = cell_of(ci)
|
||
hard = strata.hard_top[cy, cx]
|
||
|
||
# 2) 부분 기둥 — sub×sub 점
|
||
fr = ((np.arange(sub) + 0.5) / sub - 0.5) if sub > 1 else np.zeros(1)
|
||
pts: dict[str, list] = {k: [] for k in ("c", "depth", "zs", "gx", "gy", "px", "py")}
|
||
slot = np.full((sub, sub, len(ci)), -1, np.int32) # (ka, kb, 셀) → 점 번호
|
||
n_pt = 0
|
||
for ka, fy in enumerate(fr):
|
||
for kb, fx in enumerate(fr):
|
||
xs = strata.x[0] + (cx + fx) * dx
|
||
ys = strata.y[0] + (cy + fy) * dx
|
||
rem, z_s, gx, gy = surface(tt, xs, ys)
|
||
zg = strata.ground_at(cy, cx, fy, fx)
|
||
z_b = np.maximum(z_s, hard)
|
||
k = np.flatnonzero((rem > 0) & (z_b < zg))
|
||
pts["c"].append(k)
|
||
slot[ka, kb, k] = n_pt + np.arange(len(k))
|
||
n_pt += len(k)
|
||
for key, val in (
|
||
("depth", zg - z_b),
|
||
("zs", z_s),
|
||
("gx", gx),
|
||
("gy", gy),
|
||
("px", xs),
|
||
("py", ys),
|
||
):
|
||
pts[key].append(val[k])
|
||
pt = {key: np.concatenate(v) for key, v in pts.items()}
|
||
if not n_pt:
|
||
return empty
|
||
lab_pt, n_lab = _split_pieces(
|
||
strata, slot, ci, n_t * nw * nw, nw, tt, cy, cx, hard, sub, surface, n_pt
|
||
)
|
||
# 기둥 = (셀, 덩이) 짝 — 두 덩이가 한 셀을 나누면 기둥 둘
|
||
key, inv = np.unique(pt["c"] * n_lab + lab_pt, return_inverse=True)
|
||
cell_k, lab = key // n_lab, key % n_lab + 1 # 덩이 번호 1.. (0 = 없음)
|
||
n_act = np.bincount(inv, minlength=len(key)).astype(float)
|
||
acc = {
|
||
k: np.bincount(inv, weights=pt[k], minlength=len(key))
|
||
for k in ("depth", "zs", "gx", "gy", "px", "py")
|
||
}
|
||
tt, cy, cx = tt[cell_k], cy[cell_k], cx[cell_k]
|
||
a_h = dx * dx
|
||
depth_cell = acc["depth"] / (sub * sub)
|
||
vol = np.bincount(lab, weights=depth_cell, minlength=n_lab + 1) * a_h
|
||
area = np.bincount(lab, weights=n_act / (sub * sub), minlength=n_lab + 1) * a_h
|
||
edge = np.bincount(lab, weights=strata.edge[cy, cx].astype(float), minlength=n_lab + 1)
|
||
trial_of = np.zeros(n_lab + 1, int)
|
||
trial_of[lab] = tt
|
||
np.maximum.at(big, trial_of[1:], vol[1:])
|
||
ok = (vol >= v_min) & (vol <= v_max) & (area >= a_min) & (area <= a_max) & (edge == 0)
|
||
ok[0] = False
|
||
keep_lab = np.flatnonzero(ok)
|
||
remap = np.full(n_lab + 1, -1)
|
||
remap[keep_lab] = np.arange(len(keep_lab))
|
||
sel = remap[lab] >= 0
|
||
na_s = n_act[sel]
|
||
mean = {k: acc[k][sel] / na_s for k in ("zs", "gx", "gy", "px", "py")}
|
||
frac = na_s / (sub * sub)
|
||
trial = trial_of[keep_lab]
|
||
cols = Columns(
|
||
iy=cy[sel],
|
||
ix=cx[sel],
|
||
group=remap[lab][sel],
|
||
z_slip=mean["zs"],
|
||
gx=mean["gx"],
|
||
gy=mean["gy"],
|
||
center=centers[trial],
|
||
direction=None if d is None else d[trial],
|
||
px=mean["px"],
|
||
py=mean["py"],
|
||
frac=frac,
|
||
depth=depth_cell[sel] / frac,
|
||
)
|
||
return cols, trial, big
|
||
|
||
|
||
_PART_OPTS = ("v_min", "v_max", "a_min", "a_max", "sub", "groundwater", "seismic_kh")
|
||
_PART: dict = {} # 일꾼 프로세스의 지층 · 조건
|
||
|
||
|
||
def _part_init(strata: Strata | str, opts: dict) -> None:
|
||
if isinstance(strata, str):
|
||
with open(strata, "rb") as fh:
|
||
strata = pickle.load(fh)
|
||
_PART.update(strata=strata, **opts)
|
||
|
||
|
||
def _open_pool(strata: Strata, opts: dict, n_workers: int) -> ProcessPoolExecutor:
|
||
"""일꾼 풀 — 지층은 임시 파일로 넘김(spawn 관으로 넘기면 일꾼이 차례로 뜸 · 12 개 6 s)."""
|
||
fd, path = tempfile.mkstemp(suffix=".e04.pkl")
|
||
with os.fdopen(fd, "wb") as fh:
|
||
pickle.dump(strata, fh, protocol=pickle.HIGHEST_PROTOCOL)
|
||
pool = ProcessPoolExecutor(n_workers, initializer=_part_init, initargs=(path, opts))
|
||
weakref.finalize(pool, os.remove, path) # 멈춤 · 오류로 빠져나가도 지움
|
||
return pool
|
||
|
||
|
||
def _part(job: tuple, k: dict | None = None) -> tuple:
|
||
"""후보 덩이 하나 — 기둥 + 솔버. 셀 합치기에 쓰는 칸(iy · ix · group)만 돌려줌(옮길 양 줄임)."""
|
||
xyz, dirs, r, ratio = job
|
||
k = k or _PART
|
||
st = k["strata"]
|
||
cols, trial, big = trial_columns(
|
||
st, xyz, r, ratio, dirs, k["v_min"], k["v_max"], k["a_min"], k["a_max"], k["sub"]
|
||
)
|
||
res = solve(st, cols, k["groundwater"], k["seismic_kh"]) if len(trial) else None
|
||
light = SimpleNamespace(
|
||
iy=cols.iy.astype(np.int32), ix=cols.ix.astype(np.int32), group=cols.group.astype(np.int32)
|
||
)
|
||
return light, trial, big, res
|
||
|
||
|
||
def _split_pieces(strata, slot, ci, n_win, nw, tt, cy, cx, hard, sub, surface, n_pt):
|
||
"""활동면이 지표를 자른 점들을 덩이로 — 이웃 점(4 방향)이 둘 다 땅속이고 사이도 땅속이면 이음.
|
||
|
||
한 셀 안에서도 덩이가 갈림(발끝을 살짝 비켜 지나는 구 = 사면 몸통 · 발끝 앞 접시 따로).
|
||
지표는 셀 중심 사이 쌍선형(축 방향 직선) · 활동면은 아래로 볼록 → 두 점 사이 깊이는
|
||
지표 꺾임(셀 중심 줄)에서만 끝점보다 낮아짐 → 셀 중심 줄을 사이에 둔 짝만 그 줄에서 봄.
|
||
slot = (ka, kb, 셀) → 점 번호(−1 없음) · ci = 창 칸 번호 · n_win = 창 칸 수.
|
||
반환 (점별 덩이 번호, 수).
|
||
"""
|
||
from scipy.sparse import coo_matrix
|
||
from scipy.sparse.csgraph import connected_components
|
||
|
||
fr = ((np.arange(sub) + 0.5) / sub - 0.5) if sub > 1 else np.zeros(1)
|
||
rows, cols_ = [], []
|
||
|
||
def link(a, b):
|
||
both = (a >= 0) & (b >= 0)
|
||
rows.append(a[both])
|
||
cols_.append(b[both])
|
||
|
||
# 셀 안 — 꺾임(f = 0)을 사이에 둔 짝은 꺾임 점에서 땅속인지 봄
|
||
for k in range(sub - 1):
|
||
for m in range(sub):
|
||
for axis in ("x", "y"):
|
||
a = slot[m, k] if axis == "x" else slot[k, m]
|
||
b = slot[m, k + 1] if axis == "x" else slot[k + 1, m]
|
||
c = np.flatnonzero((a >= 0) & (b >= 0))
|
||
if fr[k] < 0 < fr[k + 1]:
|
||
fy, fx = (fr[m], 0.0) if axis == "x" else (0.0, fr[m])
|
||
xs = strata.x[0] + (cx[c] + fx) * strata.dx
|
||
ys = strata.y[0] + (cy[c] + fy) * strata.dx
|
||
rem, z_s, _, _ = surface(tt[c], xs, ys)
|
||
zg = strata.ground_at(cy[c], cx[c], fy, fx)
|
||
c = c[(rem > 0) & (np.maximum(z_s, hard[c]) < zg)]
|
||
rows.append(a[c])
|
||
cols_.append(b[c])
|
||
# 셀 사이 — 이웃 셀(창 안 · 줄 넘김 막음)의 맞닿은 점끼리
|
||
cell_at = np.full(n_win + nw, -1, np.int32) # 창 칸 → 셀 번호(끝 넘침 칸 −1)
|
||
cell_at[ci] = np.arange(len(ci), dtype=np.int32)
|
||
for step, pos in ((1, ci % nw), (nw, (ci // nw) % nw)):
|
||
j = cell_at[ci + step]
|
||
nb = np.flatnonzero((j >= 0) & (pos < nw - 1))
|
||
j = j[nb]
|
||
for m in range(sub):
|
||
if step == 1:
|
||
link(slot[m, sub - 1, nb], slot[m, 0, j])
|
||
else:
|
||
link(slot[sub - 1, m, nb], slot[0, m, j])
|
||
r_ = np.concatenate(rows)
|
||
c_ = np.concatenate(cols_)
|
||
graph = coo_matrix((np.ones(len(r_), np.int8), (r_, c_)), shape=(n_pt, n_pt))
|
||
n_lab, lab = connected_components(graph, directed=False)
|
||
return lab, n_lab
|
||
|
||
|
||
def governing(res: dict) -> tuple[np.ndarray, np.ndarray]:
|
||
"""다스리는 안전율 · 방법(0 Bishop · 1 Janbu 보정) — 머리 주석 규칙."""
|
||
jc = np.where(np.isfinite(res["janbu_corrected"]), res["janbu_corrected"], np.inf)
|
||
use_j = (res["truncated_ratio"] >= TRUNCATION_JANBU) & (jc < res["bishop"])
|
||
f = np.where(use_j, jc, res["bishop"])
|
||
f = np.where(np.isfinite(res["bishop"]), f, np.nan)
|
||
return f, use_j.astype(np.int8)
|
||
|
||
|
||
def evaluate_surface(
|
||
strata: Strata,
|
||
center,
|
||
r: float,
|
||
ratio: float = 1.0,
|
||
direction=None,
|
||
groundwater: dict | None = None,
|
||
seismic_kh: float = 0.0,
|
||
sub: int = 2,
|
||
) -> dict:
|
||
"""활동면 하나(구 · 타원체) — 가장 큰 덩이의 방법별 안전율 · 부피."""
|
||
dirs = None if direction is None else np.asarray(direction, float)[None]
|
||
cols, _, _ = trial_columns(strata, np.asarray(center, float)[None], r, ratio, dirs, sub=sub)
|
||
if not len(cols.group):
|
||
raise ValueError("활동면이 지표를 자르지 않음")
|
||
res = solve(strata, cols, groundwater, seismic_kh)
|
||
k = int(np.argmax(res["volume"]))
|
||
out = {key: (v[k] if isinstance(v, np.ndarray) else v) for key, v in res.items()}
|
||
out["governing"], out["method"] = (float(x[k]) for x in governing(res))
|
||
return out
|
||
|
||
|
||
def _local_fall(strata: Strata) -> tuple[np.ndarray, np.ndarray]:
|
||
zs = ndimage.uniform_filter(np.nan_to_num(strata.z_top, nan=float(np.nanmean(strata.z_top))), 5)
|
||
gy, gx = np.gradient(zs, strata.dx)
|
||
norm = np.hypot(gx, gy)
|
||
norm[norm == 0] = 1.0
|
||
return -gx / norm, -gy / norm
|
||
|
||
|
||
def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
|
||
"""거친 격자 → 고운 격자(Scoops3D 3.2.2) 로 셀별 최소 안전율 · 임계 체적."""
|
||
p = params or SearchParams()
|
||
t0 = time.perf_counter()
|
||
ny, nx = strata.z_top.shape
|
||
x0, x1, y0, y1 = p.box or (strata.x[0], strata.x[-1], strata.y[0], strata.y[-1])
|
||
nxs = int((x1 - x0) // p.spacing) + 1
|
||
nys = int((y1 - y0) // p.spacing) + 1
|
||
nhs = int(round((p.h_max - p.h_min) / p.h_step)) + 1
|
||
fall_x, fall_y = _local_fall(strata)
|
||
|
||
f_min = np.full((ny, nx), np.inf)
|
||
volume = np.full((ny, nx), np.nan)
|
||
method = np.full((ny, nx), -1, np.int8)
|
||
node_map = np.full((ny, nx), -1, np.int64)
|
||
best = {"governing": np.inf}
|
||
at_best: dict = {} # 전역 최소의 노드 번호 · 방향(다듬기용 · 결과에 안 냄)
|
||
seen: set[int] = set()
|
||
counters = {"trials": 0, "evaluated": 0}
|
||
stage = [0, 40, "거친 탐색"] # 진행률 — 시작 % · 폭 % · 이름
|
||
|
||
def tick(done: float) -> None:
|
||
if p.should_stop and p.should_stop():
|
||
raise InterruptedError
|
||
if p.on_progress:
|
||
p.on_progress(int(stage[0] + stage[1] * done), str(stage[2]))
|
||
|
||
def node_xyz(ids: np.ndarray) -> np.ndarray:
|
||
ix = ids % nxs
|
||
iy = (ids // nxs) % nys
|
||
ih = ids // (nxs * nys)
|
||
x = x0 + ix * p.spacing
|
||
y = y0 + iy * p.spacing
|
||
cx = np.clip(np.rint((x - strata.x[0]) / strata.dx).astype(int), 0, nx - 1)
|
||
cy = np.clip(np.rint((y - strata.y[0]) / strata.dx).astype(int), 0, ny - 1)
|
||
z = strata.z_top[cy, cx] + p.h_min + ih * p.h_step
|
||
return np.column_stack([x, y, z]), cy, cx, p.h_min + ih * p.h_step
|
||
|
||
# 후보 덩이 = 프로세스로 나눠 셈(_part) · 결과는 낸 차례대로 합침 → 일꾼 수와 무관한 같은 값
|
||
n_workers = p.workers or (
|
||
min(12, max(1, (os.cpu_count() or 2) - 2)) if ny * nx >= 50_000 else 1
|
||
)
|
||
opts = {k: getattr(p, k) for k in _PART_OPTS}
|
||
pool = _open_pool(strata, opts, n_workers) if n_workers > 1 else None
|
||
local = {"strata": strata, **opts}
|
||
|
||
def run(ids: np.ndarray, ratios: tuple[float, ...] = (1.0,)) -> set[int]:
|
||
"""노드들의 모든 반지름 · 비 — 셀 최소를 고친 노드 번호를 돌려줌."""
|
||
improved: set[int] = set()
|
||
xyz0, cy0, cx0, _ = node_xyz(ids)
|
||
keep = np.isfinite(xyz0[:, 2])
|
||
if p.node_mask is not None:
|
||
keep &= p.node_mask[cy0, cx0]
|
||
ids = ids[keep]
|
||
if not len(ids):
|
||
return improved
|
||
xyz, cy, cx, h = node_xyz(ids)
|
||
dirs = np.column_stack([fall_x[cy, cx], fall_y[cy, cx]])
|
||
r_lo = np.maximum(p.dr, np.floor(0.5 * h / p.dr) * p.dr)
|
||
r_hi = h + p.max_depth
|
||
for k_ratio, ratio in enumerate(ratios):
|
||
done = np.zeros(len(ids), bool)
|
||
radii = np.arange(r_lo.min(), r_hi.max() + 1e-9, p.dr)
|
||
for k_r, r in enumerate(radii):
|
||
tick((k_ratio + k_r / len(radii)) / len(ratios))
|
||
live = np.flatnonzero(~done & (r_lo <= r + 1e-9) & (r <= r_hi + 1e-9))
|
||
if not len(live):
|
||
continue
|
||
w_now = (2 * math.ceil(r * max(ratio, 1) / strata.dx) + 3) ** 2
|
||
# 메모리 = 일꾼 수 합쳐 BATCH_CELLS · 일꾼마다 한 덩이 이상(작은 고운 탐색도 나눔)
|
||
chunk = max(1, BATCH_CELLS // (w_now * n_workers))
|
||
chunk = min(chunk, -(-len(live) // n_workers))
|
||
parts = [live[s : s + chunk] for s in range(0, len(live), chunk)]
|
||
jobs = [(xyz[q], dirs[q] if ratio != 1.0 else None, r, ratio) for q in parts]
|
||
outs = pool.map(_part, jobs) if pool else (_part(j, local) for j in jobs)
|
||
# 고친 노드 판정 기준 = 이 반지름 시작 값(덩이 나눔과 무관)
|
||
ref = f_min.copy()
|
||
thr = max(p.refine_below, 1.2 * best["governing"])
|
||
for part, (cols, trial, big, res) in zip(parts, outs):
|
||
done[part[big > p.v_max]] = True
|
||
counters["trials"] += len(part)
|
||
if not len(trial):
|
||
continue
|
||
counters["evaluated"] += len(trial)
|
||
f, meth = governing(res)
|
||
improved |= _update(
|
||
f,
|
||
meth,
|
||
res,
|
||
cols,
|
||
trial,
|
||
ids[part],
|
||
f_min,
|
||
volume,
|
||
method,
|
||
node_map,
|
||
thr,
|
||
ref,
|
||
)
|
||
k = _set_best(best, res, f, meth, xyz[part][trial], r, ratio, cols, nx)
|
||
if k >= 0:
|
||
at_best.update(
|
||
node=int(ids[part][trial[k]]),
|
||
direction=dirs[part][trial[k]] if ratio != 1.0 else None,
|
||
)
|
||
return improved
|
||
|
||
try:
|
||
step = max(1, int(p.coarse))
|
||
grid = np.stack(
|
||
np.meshgrid(
|
||
np.arange(0, nhs, step),
|
||
np.arange(0, nys, step),
|
||
np.arange(0, nxs, step),
|
||
indexing="ij",
|
||
),
|
||
-1,
|
||
).reshape(-1, 3)
|
||
ids = (grid[:, 0] * nys + grid[:, 1]) * nxs + grid[:, 2]
|
||
seen.update(ids.tolist())
|
||
seeds = run(ids)
|
||
rounds = passes = 0
|
||
while seeds and (step > 1 or rounds < p.refine_rounds):
|
||
if step == 1:
|
||
rounds += 1
|
||
step = max(1, step // 2)
|
||
passes += 1
|
||
stage[:] = [min(40 + 8 * (passes - 1), 77), 8, f"고운 탐색 {passes}"]
|
||
cand = set()
|
||
for nid in seeds:
|
||
ih, rem = divmod(nid, nys * nxs)
|
||
iy, ix = divmod(rem, nxs)
|
||
for dh in (-step, 0, step):
|
||
for dy in (-step, 0, step):
|
||
for dxs in (-step, 0, step):
|
||
h2, y2, x2 = ih + dh, iy + dy, ix + dxs
|
||
if 0 <= h2 < nhs and 0 <= y2 < nys and 0 <= x2 < nxs:
|
||
cid = (h2 * nys + y2) * nxs + x2
|
||
if cid not in seen:
|
||
cand.add(cid)
|
||
if not cand:
|
||
if step == 1:
|
||
break
|
||
continue
|
||
seen.update(cand)
|
||
seeds = run(np.fromiter(cand, np.int64))
|
||
|
||
# 타원체 — 구 탐색이 찾은 임계 노드(F < refine_below)에서만(시간 1/3 · 노선 1 km 실측)
|
||
ell = tuple(r for r in p.ratios if r != 1.0)
|
||
limit = max(p.refine_below, 1.2 * best["governing"])
|
||
crit = np.unique(node_map[(node_map >= 0) & (f_min < limit)])
|
||
if ell and len(crit):
|
||
stage[:] = [85, 14, "타원체"]
|
||
run(crit, ell)
|
||
|
||
if at_best:
|
||
stage[:] = [99, 1, "다듬기"]
|
||
tick(0.0)
|
||
_polish(strata, p, best, at_best, f_min, volume, method, node_map)
|
||
|
||
finally:
|
||
if pool: # 끊기(InterruptedError) · 오류에도 일꾼을 닫음
|
||
pool.shutdown(cancel_futures=True)
|
||
f_out = np.where(np.isfinite(f_min), f_min, np.nan)
|
||
return SearchResult(
|
||
f_min=f_out,
|
||
volume=volume,
|
||
method=method,
|
||
best=best,
|
||
trials=counters["trials"],
|
||
evaluated=counters["evaluated"],
|
||
seconds=time.perf_counter() - t0,
|
||
stats={"nodes": len(seen), "cells": int(np.isfinite(f_out).sum())},
|
||
)
|
||
|
||
|
||
def _set_best(
|
||
best: dict, res: dict, f, meth, centers, r: float, ratio: float, cols, nx: int
|
||
) -> int:
|
||
"""덩이들 가운데 best 보다 낮은 것이 있으면 best 를 바꾸고 그 번호 · 없으면 −1."""
|
||
if not np.isfinite(f).any():
|
||
return -1
|
||
k = int(np.nanargmin(np.where(np.isfinite(f), f, np.inf)))
|
||
if not f[k] < best["governing"]:
|
||
return -1
|
||
best.clear()
|
||
best.update(
|
||
{key: (v[k].tolist() if isinstance(v, np.ndarray) else v) for key, v in res.items()}
|
||
)
|
||
best.update(
|
||
governing=float(f[k]),
|
||
method=int(meth[k]),
|
||
center=np.asarray(centers[k], float).tolist(),
|
||
radius=float(r),
|
||
ratio=float(ratio),
|
||
# 활동체 셀(iy·nx+ix) — 한 구가 자른 여러 덩이 중 이 덩이만
|
||
cells=np.unique(cols.iy[cols.group == k] * nx + cols.ix[cols.group == k]).tolist(),
|
||
)
|
||
return k
|
||
|
||
|
||
def _polish(strata, p, best, at_best, f_min, volume, method, node_map) -> None:
|
||
"""전역 최소 다듬기 — 그 면의 중심 · 반지름을 노드 · 반지름 눈금 사이로 옮겨 봄(Nelder-Mead).
|
||
|
||
급경사 전역 최소는 발끝에 거의 닿는 면 — 눈금(노드 간격 · dr)에 걸린 값이 연속 최소보다
|
||
1~2 % 높음(60° 제방 1 m 격자 0.7595 → 0.7519 · Scoops3D 0.75). 비 · 방향은 그대로.
|
||
"""
|
||
from scipy.optimize import minimize
|
||
|
||
ratio = best["ratio"]
|
||
d = at_best["direction"]
|
||
d = None if d is None else np.asarray(d, float)[None]
|
||
|
||
def trial(v):
|
||
cols, tr, _ = trial_columns(
|
||
strata, v[None, :3], v[3], ratio, d, p.v_min, p.v_max, p.a_min, p.a_max, p.sub
|
||
)
|
||
if not len(tr):
|
||
return None
|
||
res = solve(strata, cols, p.groundwater, p.seismic_kh)
|
||
return (cols, tr, res, *governing(res))
|
||
|
||
def objective(v):
|
||
if p.node_mask is not None: # 둘레 네 셀이 모두 노드 칸이어야(가장자리 밖으로 안 샘)
|
||
fy = (v[1] - strata.y[0]) / strata.dx
|
||
fx = (v[0] - strata.x[0]) / strata.dx
|
||
iy = np.clip([math.floor(fy), math.ceil(fy)], 0, len(strata.y) - 1)
|
||
ix = np.clip([math.floor(fx), math.ceil(fx)], 0, len(strata.x) - 1)
|
||
if not p.node_mask[np.ix_(iy, ix)].all():
|
||
return np.inf
|
||
out = trial(v)
|
||
if out is None or not np.isfinite(out[3]).any():
|
||
return np.inf
|
||
return float(np.nanmin(out[3]))
|
||
|
||
v0 = np.r_[best["center"], best["radius"]]
|
||
step = np.array([p.spacing, p.spacing, p.h_step, p.dr]) / 2
|
||
simplex = np.vstack([v0, v0 + np.diag(step)])
|
||
# 찾은 노드 · 반지름에서 한 눈금 안(그 사이를 놓친 것만 고침) · box 안
|
||
lo, hi = v0 - 2 * step, v0 + 2 * step
|
||
if p.box is not None:
|
||
lo[:2] = np.maximum(lo[:2], (p.box[0], p.box[2]))
|
||
hi[:2] = np.minimum(hi[:2], (p.box[1], p.box[3]))
|
||
o = minimize(
|
||
objective,
|
||
v0,
|
||
method="Nelder-Mead",
|
||
bounds=list(zip(lo, hi)),
|
||
options={"initial_simplex": simplex, "xatol": 1e-3, "fatol": 1e-5, "maxiter": 400},
|
||
)
|
||
if not o.fun < best["governing"] * (1 - 1e-6):
|
||
return
|
||
cols, tr, res, f, meth = trial(o.x)
|
||
_update(f, meth, res, cols, tr, np.array([at_best["node"]]), f_min, volume, method, node_map)
|
||
_set_best(
|
||
best,
|
||
res,
|
||
f,
|
||
meth,
|
||
np.repeat(o.x[None, :3], len(tr), 0),
|
||
float(o.x[3]),
|
||
ratio,
|
||
cols,
|
||
strata.z_top.shape[1],
|
||
)
|
||
|
||
|
||
def _update(
|
||
f,
|
||
meth,
|
||
res,
|
||
cols,
|
||
trial,
|
||
node_ids,
|
||
f_min,
|
||
volume,
|
||
method,
|
||
node_map,
|
||
refine_below=np.inf,
|
||
ref=None,
|
||
) -> set[int]:
|
||
"""셀마다 더 낮은 F 면 바꿈 — 같은 셀 여럿이면 가장 낮은 것(내림차순으로 써 마지막이 이김).
|
||
|
||
돌려주는 노드 = ref(없으면 f_min) 보다 낮춘 후보의 노드 중 F < refine_below.
|
||
"""
|
||
fc = f[cols.group]
|
||
ok = np.isfinite(fc)
|
||
iy, ix, grp, fc = cols.iy[ok], cols.ix[ok], cols.group[ok], fc[ok]
|
||
order = np.argsort(-fc, kind="stable")
|
||
iy, ix, grp, fc = iy[order], ix[order], grp[order], fc[order]
|
||
ref = f_min if ref is None else ref
|
||
hit = (fc < ref[iy, ix] * (1 - 1e-6)) & (fc < refine_below)
|
||
seeds = set(np.unique(node_ids[trial[grp[hit]]]).tolist())
|
||
better = fc < f_min[iy, ix] * (1 - 1e-6)
|
||
iy, ix, grp, fc = iy[better], ix[better], grp[better], fc[better]
|
||
f_min[iy, ix] = fc
|
||
volume[iy, ix] = res["volume"][grp]
|
||
method[iy, ix] = meth[grp]
|
||
node_map[iy, ix] = node_ids[trial[grp]]
|
||
return seeds
|