Merge remote-tracking branch 'origin/sub_laptop_5' into sub_laptop_7

E04 탐색 병렬(834a86e76) + 진행률 · 끊기 고리 함께 · 끊기 · 오류 때 일꾼 풀 닫음(try/finally) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-26 03:23:57 +09:00
parent 6c486da9dc 834a86e76e
commit 5d6b1b45e1
3 changed files with 228 additions and 99 deletions
@@ -244,6 +244,27 @@ def _finish(res: np.ndarray, drive: np.ndarray) -> np.ndarray:
    return out


+def _fixed_point(st: dict, f: np.ndarray, coef: dict, step) -> tuple[np.ndarray, np.ndarray]:
+    """F = step(기둥 묶음, F) 고정점 — 수렴한 활동체의 기둥은 빼며 셈(값은 빼지 않을 때와 같음).
+
+    coef = 기둥별 상수 배열(group 포함) · 반환 (F, 수렴 여부).
+    """
+    n = st["n"]
+    sub = coef
+    conv = np.zeros(n, bool)
+    for _ in range(MAX_ITER):
+        f_new = step(sub, f)
+        done = np.abs(f_new - f) < TOL
+        f = np.where(conv, f, f_new)
+        conv |= done
+        if conv.all():
+            break
+        live = ~conv[sub["group"]]
+        if live.mean() < 0.6:  # 남은 기둥만 — 매번 줄이면 줄이는 값이 더 듦
+            sub = {k: v[live] for k, v in sub.items()}
+    return f, conv
+
+
 def _bishop(st: dict, f0: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    grp, n = st["group"], st["n"]
    wq = st["w"] + st["q"]
@@ -252,19 +273,29 @@ def _bishop(st: dict, f0: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    lift = np.zeros(n)
    np.maximum.at(lift, grp, neg)
    f = np.where(np.isfinite(f0), np.maximum(f0, 0.05), 1.0) + lift
-    conv = np.zeros(n, bool)
-    for _ in range(MAX_ITER):
-        s_res, _ = _s_res(st, f)
-        s = s_res / np.maximum(f[grp], 1e-3)
-        nn = (wq - s * st["sin_a"]) / st["cos_psi"]  # 연직 평형에서 N
-        res = _gsum(grp, st["arm_s"] * s_res, n)
-        drive = _gsum(grp, -st["p"] * wq + st["e"] * st["kh"] * st["w"] + nn * st["m_n"], n)
-        f_new = _finish(res, drive)
-        done = np.abs(f_new - f) < TOL
-        f = np.where(conv, f, f_new)
-        conv |= done
-        if conv.all():
-            break
+    # S_res = num / mα · N = (W+Q − S sinα)/cosψ → 구동 = d1 − Σ S_res·k_n / F (기둥별 상수 미리)
+    num = st["c"] * st["a_h"] + np.maximum(wq - st["u_ah"], 0.0) * st["tan_phi"]
+    d1 = _gsum(
+        grp, -st["p"] * wq + st["e"] * st["kh"] * st["w"] + wq * st["m_n"] / st["cos_psi"], n
+    )
+    coef = {
+        "group": grp,
+        "num": num,
+        "cos_psi": st["cos_psi"],
+        "k": st["sin_a"] * st["tan_phi"],
+        "arm_s": st["arm_s"],
+        "k_n": st["sin_a"] * st["m_n"] / st["cos_psi"],
+    }
+
+    def step(c: dict, f: np.ndarray) -> np.ndarray:
+        g = c["group"]
+        fg = np.maximum(f, 1e-3)
+        s_res = c["num"] / (c["cos_psi"] + c["k"] / fg[g])
+        res = _gsum(g, c["arm_s"] * s_res, n)
+        drive = d1 - _gsum(g, s_res * c["k_n"], n) / fg
+        return _finish(res, drive)
+
+    f, conv = _fixed_point(st, f, coef, step)
    _, m_a = _s_res(st, f)
    ma_min = _gmin(grp, m_a, n)
    # 못 쓰는 면 — 수렴 안 함 · mα ≤ 0(바닥 수직력이 음) · F ≤ 0 → nan(후보에서 빠짐)
@@ -280,17 +311,23 @@ def _gmin(grp: np.ndarray, values: np.ndarray, n: int) -> np.ndarray:

 def _janbu(st: dict, f0: np.ndarray) -> np.ndarray:
    grp, n = st["group"], st["n"]
-    drive = _gsum(grp, (st["w"] + st["q"]) * st["tan_a"] + st["kh"] * st["w"], n)
+    wq = st["w"] + st["q"]
+    drive = _gsum(grp, wq * st["tan_a"] + st["kh"] * st["w"], n)
    f = np.where(np.isfinite(f0), np.maximum(f0, 0.05), 1.0)
-    conv = np.zeros(n, bool)
-    for _ in range(MAX_ITER):
-        s_res, _ = _s_res(st, f)
-        f_new = _finish(_gsum(grp, s_res / st["cos_a"], n), drive)
-        done = np.abs(f_new - f) < TOL
-        f = np.where(conv, f, f_new)
-        conv |= done
-        if conv.all():
-            break
+    num = st["c"] * st["a_h"] + np.maximum(wq - st["u_ah"], 0.0) * st["tan_phi"]
+    coef = {
+        "group": grp,
+        "num": num / st["cos_a"],
+        "cos_psi": st["cos_psi"],
+        "k": st["sin_a"] * st["tan_phi"],
+    }
+
+    def step(c: dict, f: np.ndarray) -> np.ndarray:
+        g = c["group"]
+        s_res = c["num"] / (c["cos_psi"] + c["k"] / np.maximum(f, 1e-3)[g])
+        return _finish(_gsum(g, s_res, n), drive)
+
+    f, conv = _fixed_point(st, f, coef, step)
    f[~conv | (f <= 0)] = np.nan
    return f

@@ -9,6 +9,7 @@
  구는 활동체 지표 평균 낙하 방향(Scoops3D 2.1).
 - 한 구가 지표를 여러 덩이로 자르면 덩이마다 셈(ndimage.label · 8 이웃).
 - 격자 가장자리 · DEM 빈칸에 닿은 덩이는 버림(잘린 질량 막음).
+- 후보 덩이는 프로세스로 나눠 셈 · 낸 차례대로 합침 — 일꾼 수와 무관하게 같은 값.

 다스리는 안전율(governing):
 - 못 뚫는 층(기본 발파암 = 연암)에 잘린 바닥 면적 비 ≥ 10 % → Janbu 보정으로 전환
@@ -21,9 +22,15 @@
 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
@@ -57,6 +64,7 @@ class SearchParams:
    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
@@ -107,13 +115,23 @@ def trial_columns(
    ixg = ix0[:, None, None] + off[None, None, :]
    iyg = iy0[:, None, None] + off[None, :, None]
    inb = (ixg >= 0) & (ixg < nx) & (iyg >= 0) & (iyg < ny)
-    ixc, iyc, inb = np.broadcast_arrays(np.clip(ixg, 0, nx - 1), np.clip(iyg, 0, ny - 1), inb)
-    shape = ixc.shape
+    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). 원소별(납작 · 창 모두)."""
+        """후보 t 의 활동면 — (rem, z_s, ∂z/∂x, ∂z/∂y). 원소별(납작)."""
        dxp = xs - centers[t, 0]
        dyp = ys - centers[t, 1]
        if d is None:
@@ -130,13 +148,16 @@ def trial_columns(
        return rem, centers[t, 2] - root, gx, gy

    # 1) 셀 중심으로 거르기 — 걸친 셀과 그 이웃만 셀 안 표본을 봄(창 대부분은 허공 · 땅속 밖)
-    t_grid = np.broadcast_to(np.arange(n_t)[:, None, None], shape)
-    xs0 = strata.x[0] + ixg * dx
-    ys0 = strata.y[0] + iyg * dx
-    rem0, zs0, _, _ = surface(t_grid, xs0, ys0)
+    # 창 칸 번호 → (후보, 셀) — 큰 창 배열을 펴지 않고 셈
+    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))
-    yy0, xx0 = iyc.ravel()[in0], ixc.ravel()[in0]
-    zb0 = np.maximum(zs0.ravel()[in0], strata.hard_top[yy0, xx0])
+    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
@@ -146,9 +167,7 @@ def trial_columns(
    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 = t_grid.ravel()[ci]
-    cy = iyc.ravel()[ci]
-    cx = ixc.ravel()[ci]
+    tt, cy, cx = cell_of(ci)
    hard = strata.hard_top[cy, cx]

    # 2) 부분 기둥 — sub×sub 점
@@ -211,6 +230,42 @@ def trial_columns(
    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 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)
@@ -287,6 +342,14 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
        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()
@@ -310,27 +373,20 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
                if not len(live):
                    continue
                w_now = (2 * math.ceil(r * max(ratio, 1) / strata.dx) + 3) ** 2
-                chunk = max(1, BATCH_CELLS // max(w_now, 1))
-                for s in range(0, len(live), chunk):
-                    part = live[s : s + chunk]
-                    d_part = dirs[part] if ratio != 1.0 else None
-                    cols, trial, big = trial_columns(
-                        strata,
-                        xyz[part],
-                        r,
-                        ratio,
-                        d_part,
-                        p.v_min,
-                        p.v_max,
-                        p.a_min,
-                        p.a_max,
-                        p.sub,
-                    )
+                # 메모리 = 일꾼 수 합쳐 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
-                    res = solve(strata, cols, p.groundwater, p.seismic_kh)
                    counters["evaluated"] += len(trial)
                    f, meth = governing(res)
                    improved |= _update(
@@ -344,7 +400,8 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
                        volume,
                        method,
                        node_map,
-                        max(p.refine_below, 1.2 * best["governing"]),
+                        thr,
+                        ref,
                    )
                    k = (
                        int(np.nanargmin(np.where(np.isfinite(f), f, np.inf)))
@@ -368,50 +425,57 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
                        )
        return improved

-    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:
+    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:
-                break
-            continue
-        seen.update(cand)
-        seeds = run(np.fromiter(cand, np.int64))
+                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)
+        # 타원체 — 구 탐색이 찾은 임계 노드(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)

+    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,
@@ -426,19 +490,35 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:


 def _update(
-    f, meth, res, cols, trial, node_ids, f_min, volume, method, node_map, refine_below=np.inf
+    f,
+    meth,
+    res,
+    cols,
+    trial,
+    node_ids,
+    f_min,
+    volume,
+    method,
+    node_map,
+    refine_below=np.inf,
+    ref=None,
 ) -> set[int]:
-    """셀마다 더 낮은 F 면 바꿈 — 같은 셀 여럿이면 가장 낮은 것(내림차순으로 써 마지막이 이김)."""
+    """셀마다 더 낮은 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]
-    nid = node_ids[trial[grp]]
-    node_map[iy, ix] = nid
-    return set(np.unique(nid[fc < refine_below]).tolist())
+    node_map[iy, ix] = node_ids[trial[grp]]
+    return seeds
@@ -85,3 +85,15 @@ def test_wet_cone_pore_pressure_ratio():
        groundwater={"mode": "ru", "ru": 0.4},
    )
    assert search(st, p).best["governing"] == pytest.approx(1.626, rel=0.03)
+
+
+def test_workers_same_result():
+    """프로세스 나눔 — 일꾼 수와 무관하게 셀별 최소 F · 전역 최소가 같음."""
+    st = _embankment(10 / np.tan(np.radians(30)), 19, 20, 8.64)
+    kw = dict(box=(-2, 10, -4, 4), h_min=4, h_max=20, spacing=2.0, coarse=2, max_depth=16)
+    kw.update(v_min=50, v_max=2e4, ratios=(1.0, 2.0))
+    one = search(st, SearchParams(workers=1, **kw))
+    two = search(st, SearchParams(workers=2, **kw))
+    np.testing.assert_array_equal(one.f_min, two.f_min)
+    assert one.best["governing"] == two.best["governing"]
+    assert one.evaluated == two.evaluated