diff --git a/B05_Profile/B05_Profile_Engine_RouteInitial.py b/B05_Profile/B05_Profile_Engine_RouteInitial.py index d6340207..491b5674 100644 --- a/B05_Profile/B05_Profile_Engine_RouteInitial.py +++ b/B05_Profile/B05_Profile_Engine_RouteInitial.py @@ -1,27 +1,31 @@ -"""초기 계획노선 갈래 — 예상노선 추종 · 토공 최소 · 기울기 순한 노선 · 비중(PLAN 20장). +"""초기 계획노선 갈래 — 예상노선 추종 · 토공 최소 · 기울기 순한 노선 · 비중(PLAN 20 · 23장). 계산 자리 = 서버 단독(③). 결과는 화면 캐시용 — 이 모듈은 파일을 쓰지 않는다. 적용은 화면이 받은 노드로 기존 `POST /route/replan`([확인])을 부른다. 갈래(`mode`): - · follow — 예상노선 이격 최소(예상노선을 따르되 공통 제약은 지킴). - · earthwork — 절 · 성토 합 최소(DEM 옆경사로 가늠한 단면 토량). - · gentle — 종단기울기 최소(지반 기울기 제곱 합). - · weighted — 예상노선 이격 · 토공 · 기울기 세 목표를 **정규화해 가중합**(합 100). - 정규화 = 목표마다 복도 안 링크 평균값으로 나눔 — 단위가 다른 세 값을 1 m 당 - 「보통 크기 1」로 맞춘 뒤 비중을 곱한다. -탐색 = 예상노선 둘레 복도(`CORRIDOR_M`) 안 `SEARCH_CELL_M` 격자 · 16 방향 Dijkstra(scipy csgraph). -상태 = (칸, 진행 방향, 꺾은 뒤 곧게 온 걸음 수). 탐색 경로의 꺾임점이 곧 노드이고, 노드마다 -최소 곡선반지름 원호를 끼운다(`build_planned_polyline` · 노드별 곡선 — [확인] 재구성과 같은 규칙). + · follow — 예상노선 추적 = 전처리 초기값(`planned_route_initial.csv`)과 같은 규칙 · 같은 선. + 제약 탐색은 하지 않고, 넘은 곳은 `metrics.violations` 에 뜬다. + · earthwork — 후보 중 절 · 성토 합(`cut_m3 + fill_m3`) 최소. + · gentle — 후보 중 최대 종단기울기 최소(같으면 평균 최소). + · weighted — 후보 중 예상노선 이격 · 절성토 · 기울기(최대 + 평균) 세 값을 후보 범위로 + 0~1 정규화해 비중(합 100)으로 더한 값 최소. -공통 제약(모든 갈래 · **넘지 못함**) — 종단기울기 상한 · 최소 곡선반지름 · 계류 이격: +후보 = 목표 비중(토공 · 기울기 · 이격 · 기본 비중)을 바꿔 가며 격자 탐색(`..._Search`)한 노선들. +갈래는 **비교값(화면에 보이는 값)으로 고른다** — 탐색 목표(4 m 링크 기울기 · 옆경사)는 비교값 +(측점 20 m 기울기 · 절성토 가늠)과 어긋나 탐색 한 벌로는 갈래 순서가 안 맞는다(PLAN 23-1). + +공통 제약(follow 밖 모든 갈래 · **넘지 못함**) — 종단기울기 상한 · 최소 곡선반지름 · 계류 이격: · 종단기울기 — 지반 기울기가 상한을 넘는 링크는 탐색에서 뺌. 곡선 맞춤 뒤 측점(20 m) 사이 지반이 상한을 넘으면 탐색 상한을 조금씩 낮춰 다시 찾음(`GRADE_CAP_STEPS`). · 최소 곡선반지름 — 한 번에 이웃 방향(약 22.5°)으로만 꺾고, 꺾은 뒤 원호 두 개의 접선 길이만큼 곧게 가야 다시 꺾음 → 노드마다 최소 반지름 원호가 들어감. · 계류 이격 — 버퍼 안에서는 계류 쪽 · 반대쪽 ±45° 로만 움직임(건너기 · 빠져나가기). - · 그래도 길이 없을 때만(지형상 불가능) 벌점 탐색으로 가장 적게 넘긴 길을 내고, 넘긴 곳을 - `metrics.violations`(측점 · 값 · 사유)에 적음. `over_limit` 수 = 그 목록 수. + · 그래도 길이 없을 때만(지형상 불가능) 벌점 탐색으로 넘긴 길을 내고, 넘긴 곳을 + `metrics.violations`(측점 · 값 · 사유)에 적음. `over_limit` 수 = 그 목록 수. 이때 후보는 + 벌점 크기를 바꿔 가며 만들고(`RELAXED_GRADE_COSTS`), 가장 센 벌점(`SOFT_GRADE_COST`)의 + 기준 길보다 못 지킨 곳이 많은 후보는 뺀다 — 센 벌점 하나로만 찾으면 벌점이 목표를 덮어 + 네 갈래가 한 길로 모인다(PLAN 23-1 원인). 기준값 확정(2026-09-29 사용자 결정, PLAN 20장) — 아래 임시값이 그대로 **기본값**이 됐다 (지식DB `resources/knowledge/technical_info/01_임도` 후보 표 그대로). 종단기울기 상한 · @@ -33,34 +37,32 @@ from __future__ import annotations +import hashlib import json -import math -from dataclasses import dataclass from pathlib import Path from typing import Any import numpy as np -from scipy import ndimage -from scipy.sparse import csr_matrix -from scipy.sparse.csgraph import connected_components, dijkstra -from B05_Profile.B05_Profile_Engine_RouteInitial_Metrics import ( - GRADE_SMOOTH_M, - route_metrics, - smoothed_ground, +from B05_Profile.B05_Profile_Engine_RouteInitial_Metrics import GRADE_SMOOTH_M, route_metrics +from B05_Profile.B05_Profile_Engine_RouteInitial_Search import ( + CORRIDOR_M, + ROAD_WIDTH_M, + NoRoute, + RouteGraph, + SearchGrid, + Terrain, ) from common_util.common_util_route_polyline import build_planned_polyline +__all__ = ["NoRoute", "Terrain", "generate_initial_route"] + # ── 사용자 확정 대기(지식DB 후보 · PLAN 20장) ───────────────────────────────────── # 계류 이격(m) — 법정 이격 수치 없음. 타당성평가 「산지계류 100 m 이내 구간 비율」(별표1 비고, # 1/25,000 주요 수계)을 후보로 씀. STREAM_OFFSET_M = 100.0 # 이격을 따질 계류 — 도엽 하천중심선 `구분`. 「주요 수계」 가늠으로 세류 · 연결선은 뺌. STREAM_CLASSES = ("소하천", "지방하천", "국가하천") -# 토공 가늠 폭(m) = 유효너비 3 + 길어깨 0.5 × 2(별표2 간선). -ROAD_WIDTH_M = 4.0 -# 탐색 복도 반폭(m) — 예상노선에서 이보다 멀리는 가지 않음(법령 값 아님 · 알고리즘 값). -CORRIDOR_M = 100.0 # 기준값 확정(2026-09-29 사용자 결정) — 예전엔 여기 여섯 칸이 「사용자 확정 대기」였다. # 이제 후보 표 값 그대로 기본값으로 굳었고(화면에서 프로젝트마다 바꿀 수 있음), 빈 목록이 곧 # 「확정됨」이다(`criteria.pending`). @@ -74,45 +76,17 @@ MODE_WEIGHTS = { "earthwork": {"expected": 0.0, "earthwork": 100.0, "grade": 0.0}, "gentle": {"expected": 0.0, "earthwork": 0.0, "grade": 100.0}, } -SEARCH_CELL_M = 4.0 # 탐색 격자 칸(m) — 비용면(2 m)을 솎아 씀 -LENGTH_COST = 0.05 # 1 m 당 거리 비용 — 목표가 평평한 곳에서 헤매지 않게 -TURN_COST = 1.0 # 꺾을 때 벌점 = 최소 곡선반지름 원호 길이(R·θ) × 이 값 — 곧게 달리기 선호 +# 후보를 만들 탐색 비중 — 목표 하나씩 +CANDIDATE_WEIGHTS = (MODE_WEIGHTS["earthwork"], MODE_WEIGHTS["gentle"], MODE_WEIGHTS["follow"]) # 탐색 상한 = 상한 × 이 값들 차례로 — 곡선 맞춤 뒤 측점 사이 지반이 넘으면 다음 값으로 다시 찾음 GRADE_CAP_STEPS = (1.0, 0.93, 0.85) -# 모든 제약을 지키는 길이 없을 때만 쓰는 벌점(1 m 당) — 넘긴 곳을 가장 적게 -SOFT_GRADE_COST = 500.0 # × (넘은 기울기 / 상한) — 50 · 500 · 5000 에서 넘는 곳 수 같음(랩탑_보조) -SOFT_STREAM_COST = 20.0 # 계류 버퍼 안 나란히 달리기 -BOUND_BISECT = 7 # 지형상 불가능할 때 막는 기울기를 찾는 이분 탐색 횟수(구간 1/128 까지) -STREAM_CROSS_COS = math.cos(math.radians(45.0)) # 버퍼 안 이동 방향 한도(계류 쪽 · 반대쪽 ±45°) -# 16 방향(행, 열) — +x 에서 반시계. 나이트 걸음(1, 2) 등을 넣어 22.5° 안팎 방향을 곧게 달림. -DIRECTIONS = ( - (0, 1), - (1, 2), - (1, 1), - (2, 1), - (1, 0), - (2, -1), - (1, -1), - (1, -2), - (0, -1), - (-1, -2), - (-1, -1), - (-2, -1), - (-1, 0), - (-2, 1), - (-1, 1), - (-1, 2), -) - - -@dataclass -class Terrain: - """지반 격자 — 탐색 비용면(`B05_Profile_Engine_Solver._load_or_build_cost_surface`) 그대로.""" - - x: np.ndarray # 열 좌표(오름차순) - y: np.ndarray # 행 좌표(오름차순) - z: np.ndarray # [행, 열] 표고 - valid: np.ndarray # [행, 열] 지형 있음 +# 모든 제약을 지키는 길이 없을 때만 쓰는 벌점(1 m 당 × 넘은 기울기 / 상한) — 기준 길(가장 적게 +# 넘긴 길)은 센 값, 갈래 후보는 목표와 겨룰 만한 값들(목표는 1 m 당 평균 1 로 정규화됨) +SOFT_GRADE_COST = 500.0 +RELAXED_GRADE_COSTS = (50.0, 5.0) +# 후보 묶음은 갈래와 무관 — 같은 입력이면 다시 짜지 않음(최근 몇 벌만 · 프로세스 메모리) +POOL_CACHE_SIZE = 4 +_POOLS: dict[bytes, tuple[list[tuple[Any, dict[str, Any]]], float | None]] = {} def load_terrain_and_streams( @@ -214,10 +188,6 @@ def stream_lines(features: list[dict[str, Any]]) -> list[list[tuple[float, float return lines -class NoRoute(ValueError): - """제약을 지키며 시점과 종점을 잇는 길이 없음.""" - - def generate_initial_route( expected: list[tuple[float, float]], terrain: Terrain, @@ -232,9 +202,7 @@ def generate_initial_route( ) -> dict[str, Any]: """갈래 하나로 초기 계획노선을 만들고 비교값을 붙여 돌려준다. `max_grade` 는 비율(0.14). - 공통 제약을 넘지 않는 길을 먼저 찾고(탐색 상한을 `GRADE_CAP_STEPS` 차례로 낮춰 곡선 맞춤 뒤 - 지반 차까지 없앰), 그런 길이 없을 때만 벌점 탐색으로 가장 적게 넘긴 길을 낸다. - + follow 밖 갈래는 후보(`_candidates`) 중 제 목표 비교값이 가장 좋은 것을 고른다. `stream_offset_m`·`grade_smooth_m` 을 안 주면(None) 이 머리의 기본값을 쓴다 — 화면에서 프로젝트마다 겹쳐 쓸 수 있다(PLAN 20장 확정). """ @@ -245,25 +213,9 @@ def generate_initial_route( used = resolve_weights(mode, weights) stream_offset = STREAM_OFFSET_M if stream_offset_m is None else float(stream_offset_m) grade_smooth = GRADE_SMOOTH_M if grade_smooth_m is None else float(grade_smooth_m) - grid = _SearchGrid( - expected, terrain, streams, stream_offset_m=stream_offset, grade_smooth_m=grade_smooth - ) - def attempt(hard_cap: float, bound: float | None) -> tuple[Any, dict[str, Any]]: - nodes = _search_route( - grid, - expected, - used, - limit=max_grade, - hard_cap=hard_cap, - strict_stream=bound is None, - min_radius_m=min_radius_m, - ) - # 노드마다 곡선(묶지 않음) — [확인] 재구성(`/route/replan`)과 같은 규칙 → 적용 뒤도 같은 선 - outline = build_planned_polyline( - nodes, min_radius_m=min_radius_m, simplify=False, curve_flags=[True] * len(nodes) - ) - metrics = route_metrics( + def measure(outline: Any, bound: float | None, follow: bool = False) -> dict[str, Any]: + return route_metrics( outline.vertices, expected, terrain, @@ -275,26 +227,33 @@ def generate_initial_route( road_width_m=ROAD_WIDTH_M, bound_pct=None if bound is None else bound * 100.0, grade_smooth_m=grade_smooth, + follow=follow, ) - return outline, metrics - best: tuple[Any, dict[str, Any]] | None = None - for step in GRADE_CAP_STEPS: - try: - tried = attempt(max_grade * step, None) - except NoRoute: - break - if best is None or tried[1]["over_limit"]["total"] < best[1]["over_limit"]["total"]: - best = tried - if best[1]["over_limit"]["total"] == 0: - break bound = None - if best is None: - # 지형상 불가능 — 반지름 · 계류 규칙을 지키는 어떤 길도 넘어야 하는 가장 작은 기울기를 - # 구해(`_terrain_bound`) 그 위는 막고, 상한과 그 사이는 벌점으로 가장 적게 넘긴다. - bound = _terrain_bound(grid, expected, used, max_grade, min_radius_m) - best = attempt(bound + 1e-9, bound) - outline, metrics = best + if mode == "follow": + # 예상노선 추적 = 전처리 초기값(`planned_route_initial.csv`)과 같은 규칙 · 같은 입력 → + # 같은 선(`_write_planned_polyline` 기본값). 제약은 따지지 않고 넘은 곳만 목록에 뜸. + outline = build_planned_polyline(expected, min_radius_m=min_radius_m) + metrics = measure(outline, None, follow=True) + else: + key = _pool_key( + expected, terrain, streams, max_grade, min_radius_m, stream_offset, grade_smooth + ) + if key not in _POOLS: + grid = SearchGrid( + expected, + terrain, + streams, + stream_offset_m=stream_offset, + grade_smooth_m=grade_smooth, + ) + graph = RouteGraph(grid, expected, limit=max_grade, min_radius_m=min_radius_m) + while len(_POOLS) >= POOL_CACHE_SIZE: + _POOLS.pop(next(iter(_POOLS))) + _POOLS[key] = _candidates(graph, max_grade, min_radius_m, measure) + candidates, bound = _POOLS[key] + outline, metrics = _pick(mode, used, candidates) return { "mode": mode, "logic": mode, @@ -319,286 +278,98 @@ def generate_initial_route( } -# ── 탐색 ───────────────────────────────────────────────────────────────────────── - - -def _cell_size(coords: np.ndarray) -> float: - return float((coords[-1] - coords[0]) / (len(coords) - 1)) if len(coords) > 1 else 1.0 - - -def _distance_grid( - lines: list[list[tuple[float, float]]], xs: np.ndarray, ys: np.ndarray, cell: float -) -> np.ndarray: - """격자 칸마다 선까지 거리(m). 선이 격자에 안 닿으면 모두 무한.""" - marked = np.zeros((len(ys), len(xs)), dtype=bool) - for line in lines: - for (x0, y0), (x1, y1) in zip(line, line[1:]): - count = max(2, int(math.hypot(x1 - x0, y1 - y0) / (cell * 0.5)) + 1) - px = np.linspace(x0, x1, count) - py = np.linspace(y0, y1, count) - cols = np.rint((px - xs[0]) / cell).astype(int) - rows = np.rint((py - ys[0]) / cell).astype(int) - inside = (cols >= 0) & (cols < len(xs)) & (rows >= 0) & (rows < len(ys)) - marked[rows[inside], cols[inside]] = True - if not marked.any(): - return np.full(marked.shape, np.inf) - return ndimage.distance_transform_edt(~marked) * cell - - -class _SearchGrid: - """탐색 격자 — 갈래 · 탐색 상한이 바뀌어도 같은 것(복도 · 거리 · 기울기)을 한 번만 만든다.""" - - def __init__( - self, - expected: list[tuple[float, float]], - terrain: Terrain, - streams: list[list[tuple[float, float]]], - *, - stream_offset_m: float = STREAM_OFFSET_M, - grade_smooth_m: float = GRADE_SMOOTH_M, - ) -> None: - step = max(1, int(round(SEARCH_CELL_M / _cell_size(terrain.x)))) - self.xs, self.ys = terrain.x[::step], terrain.y[::step] - # 기울기 · 옆경사는 비교값과 같은 고른 지반에서 잰다(`smoothed_ground`) - self.z = smoothed_ground(terrain, grade_smooth_m)[::step, ::step] - self.stream_offset_m = stream_offset_m - self.cell = _cell_size(self.xs) - self.d_expected = _distance_grid([expected], self.xs, self.ys, self.cell) - self.d_stream = _distance_grid(streams, self.xs, self.ys, self.cell) - valid = np.asarray(terrain.valid[::step, ::step], dtype=bool) - self.corridor = (self.d_expected <= CORRIDOR_M) & valid - self.grad_y, self.grad_x = np.gradient(self.z, self.cell) - finite = np.where(np.isfinite(self.d_stream), self.d_stream, 0.0) - self.away_y, self.away_x = np.gradient(finite, self.cell) # 계류에서 멀어지는 쪽 - self.index = np.full(self.z.shape, -1, dtype=np.int64) - self.cells = int(self.corridor.sum()) - self.index[self.corridor] = np.arange(self.cells) - rows, cols = np.nonzero(self.corridor) - self.centres = np.column_stack([self.xs[cols], self.ys[rows]]) - - def nearest(self, point: tuple[float, float]) -> int: - dx = self.centres[:, 0] - point[0] - return int(np.argmin(np.hypot(dx, self.centres[:, 1] - point[1]))) - - def _pairs(self, dr: int, dc: int) -> tuple: - """한 방향으로 이웃한 칸 짝(앞 · 뒤 조각, 앞 · 뒤 번호, 둘 다 복도 안).""" - rows, cols = self.z.shape - first = (slice(max(0, -dr), rows - max(0, dr)), slice(max(0, -dc), cols - max(0, dc))) - second = (slice(max(0, dr), rows - max(0, -dr)), slice(max(0, dc), cols - max(0, -dc))) - a, b = self.index[first], self.index[second] - return first, second, a, b, (a >= 0) & (b >= 0) - - def grades(self, dr: int, dc: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """복도 안 한 방향 링크의 (앞 번호, 뒤 번호, 지반 기울기).""" - first, second, a, b, ok = self._pairs(dr, dc) - grade = np.abs(self.z[second] - self.z[first]) / (self.cell * math.hypot(dr, dc)) - return a[ok], b[ok], grade[ok] - - def links( - self, dr: int, dc: int, *, limit: float, hard_cap: float, strict_stream: bool - ) -> tuple: - """한 방향 링크 · 세 목표 · 벌점(1 m 당). `hard_cap` 넘는 기울기 링크는 뺀다. - - `limit`(상한)과 `hard_cap` 사이는 벌점 — 지형상 불가능할 때만 그 사이가 생긴다. - `strict_stream` 이면 계류 버퍼 안 나란히 달리기도 뺀다(아니면 벌점). - """ - first, second, a, b, ok = self._pairs(dr, dc) - norm = math.hypot(dr, dc) - grade = np.abs(self.z[second] - self.z[first]) / (self.cell * norm) - # 계류 버퍼 안에서는 계류 쪽 · 반대쪽 ±45° 로만(건너기 · 빠져나가기) - in_buffer = (self.d_stream[first] < self.stream_offset_m) | ( - self.d_stream[second] < self.stream_offset_m - ) - ax = self.away_x[first] + self.away_x[second] - ay = self.away_y[first] + self.away_y[second] - size = np.hypot(ax, ay) - across = np.abs(ax * dc + ay * dr) >= STREAM_CROSS_COS * norm * size - along = in_buffer & ~across & (size > 1e-6) - ok &= grade <= hard_cap - if strict_stream: - ok &= ~along - grade, along = grade[ok], along[ok] - over = np.maximum(0.0, grade - limit) - # 옆경사 = 진행 방향에 수직인 지반 기울기(두 칸 평균) - cross = np.abs( - (self.grad_x[first][ok] + self.grad_x[second][ok]) * -dr - + (self.grad_y[first][ok] + self.grad_y[second][ok]) * dc - ) / (2.0 * norm) - # 토공 가늠(1 m 당 m³) — 옆경사 위 반절 · 반성 단면 W²·s/4 + 상한 넘는 기울기가 측점 - # 반 간격(10 m) 동안 쌓는 높이 차 단면(지형상 불가능할 때만 0 이 아님) - earth = ROAD_WIDTH_M**2 * cross / 4.0 + ROAD_WIDTH_M * over * 10.0 - offset = 0.5 * (self.d_expected[first][ok] + self.d_expected[second][ok]) - penalty = SOFT_GRADE_COST * over / limit + SOFT_STREAM_COST * along - terms = np.stack([offset, earth, (grade / limit) ** 2]) - return a[ok], b[ok], self.cell * norm, terms, penalty - - -def _minimax_grade(grid: _SearchGrid, expected: list[tuple[float, float]]) -> float: - """복도 안에서 시점과 종점을 잇는 데 **어떤 길도 넘어야 하는** 가장 작은 기울기 상한(비율). - - 꺾기 규칙 없이 칸끼리만 이은 그래프에서 이분 탐색 — 그러니 실제 필요한 값의 아래 한계다. - """ - pieces = [grid.grades(dr, dc) for dr, dc in DIRECTIONS] - heads = np.concatenate([p[0] for p in pieces]) - tails = np.concatenate([p[1] for p in pieces]) - grades = np.concatenate([p[2] for p in pieces]) - start, end = grid.nearest(expected[0]), grid.nearest(expected[-1]) - - def joined(cap: float) -> bool: - keep = grades <= cap - graph = csr_matrix( - (np.ones(int(keep.sum())), (heads[keep], tails[keep])), - shape=(grid.cells, grid.cells), - ) - _, labels = connected_components(graph, directed=False) - return bool(labels[start] == labels[end]) - - low, high = 0.0, float(grades.max()) if grades.size else 0.0 - if not joined(high): - raise NoRoute("복도 안에서 시점과 종점이 이어지지 않습니다.") - for _ in range(30): - middle = 0.5 * (low + high) - low, high = (low, middle) if joined(middle) else (middle, high) - return high - - -def _terrain_bound( - grid: _SearchGrid, - expected: list[tuple[float, float]], - weights: dict[str, float], +def _candidates( + graph: RouteGraph, max_grade: float, min_radius_m: float, -) -> float: - """반지름 규칙까지 지키며 시점 · 종점을 잇는 데 **어떤 길도 넘어야 하는** 기울기(비율). + measure: Any, +) -> tuple[list[tuple[Any, dict[str, Any]]], float | None]: + """갈래가 고를 후보 노선들(폴리라인 · 비교값)과 지형상 불가능 때 막은 기울기(아니면 None). - 아래 한계 = 꺾기 규칙 없는 최소최대(`_minimax_grade`). 이어지는 값을 1.5 배씩 찾은 뒤 - 그 사이를 탐색 그래프 그대로 이분 탐색한다(계류는 벌점 — 기울기만 따짐). + 제약을 지키는 길이 있으면 비중마다 탐색 상한을 `GRADE_CAP_STEPS` 차례로 낮춰 찾은 노선들 · + 없으면 막은 기울기 안에서 벌점 크기 · 비중을 바꿔 찾은 노선들. 못 지킨 곳 수가 기준 + (엄격: 후보 최소 · 풂: 센 벌점 기준 길)보다 많은 후보는 뺀다. """ - def joined(cap: float) -> bool: - try: - _search_route( - grid, - expected, - weights, - limit=max_grade, - hard_cap=cap, - strict_stream=False, - min_radius_m=min_radius_m, - ) - except NoRoute: - return False - return True - - ceiling = max(float(grid.grades(dr, dc)[2].max(initial=0.0)) for dr, dc in DIRECTIONS) - low = max(_minimax_grade(grid, expected), max_grade) - high = low - while not joined(high): - if high >= ceiling: - raise NoRoute("복도 안에서 반지름 규칙을 지키며 시점과 종점을 잇는 길이 없습니다.") - low, high = high, min(high * 1.5, ceiling) - for _ in range(BOUND_BISECT): - if high - low <= 1e-4: - break - middle = 0.5 * (low + high) - low, high = (low, middle) if joined(middle) else (middle, high) - return high - - -def _turn_steps(min_radius_m: float, cell: float) -> int: - """꺾은 뒤 곧게 가야 할 칸 수 — 앞뒤 원호 접선이 한 직선에 들어가게. - - 원호를 끼우는 쪽(`build_planned_polyline`)은 직선의 절반까지만 접선으로 쓰므로 직선은 - 2 · R · tan(θ/2) 이상이어야 한다. θ = 이웃 방향 사이 가장 큰 각(나이트 걸음 26.57°). - 걸음은 칸 길이 단위로 센다 — 곧은 · 대각 걸음 1, 나이트 걸음 2. - """ - widest = math.atan2(1.0, 2.0) - return max(1, math.ceil(2.0 * min_radius_m * math.tan(widest / 2.0) / cell)) - - -def _search_route( - grid: _SearchGrid, - expected: list[tuple[float, float]], - weights: dict[str, float], - *, - limit: float, - hard_cap: float, - strict_stream: bool, - min_radius_m: float, -) -> list[tuple[float, float]]: - """노드(꺾임점) 목록 — 첫 · 끝은 예상노선 양 끝 그대로. - - 상태 = (칸, 방향, 층). 층 = 꺾은 뒤 곧게 온 칸 수 + 1(K + 1 에서 멈춤) · 층 0 = 시점. - 이웃 방향으로 꺾기는 K 칸 곧게 온 뒤에만(층 K + 1). 시점은 한 칸 더(K + 1 칸) — 기점을 - 격자 칸 가운데에서 예상노선 끝점으로 옮기며 첫 직선이 짧아지는 몫. 종점은 K 칸 곧게 - 들어온 상태를 먼저, 없으면 아무 상태(마지막 곡선 자리가 모자라면 목록에 뜸). - """ - count = len(DIRECTIONS) - steps = _turn_steps(min_radius_m, grid.cell) - full = steps + 1 # 꺾을 수 있는 층 - layers = full + 1 - links = [ - grid.links(dr, dc, limit=limit, hard_cap=hard_cap, strict_stream=strict_stream) - for dr, dc in DIRECTIONS - ] - if not any(link[0].size for link in links): - raise NoRoute("예상노선 둘레에 제약을 지키며 지날 수 있는 링크가 없습니다.") - - # 세 목표 정규화 — 복도 안 링크 평균이 1 이 되게(길이 가중). 평균이 0 이면 그 목표는 0. - lengths = np.concatenate([np.full(link[0].size, link[2]) for link in links]) - stacked = np.hstack([link[3] for link in links]) - scale = (stacked * lengths).sum(axis=1) / max(lengths.sum(), 1e-9) - order = ("expected", "earthwork", "grade") - share = np.array([weights[key] for key in order]) / max(sum(weights.values()), 1e-9) - factor = np.divide(share, scale, out=np.zeros_like(share), where=scale > 1e-12) - - def state(cell: np.ndarray | int, direction: int, walked: int) -> np.ndarray | int: - return (cell * count + direction) * layers + walked - - angles = [math.atan2(dr, dc) for dr, dc in DIRECTIONS] - heads, tails, costs = [], [], [] - for out, (a, b, length, terms, penalty) in enumerate(links): - move = length * (LENGTH_COST + factor @ terms + penalty) - units = max(1, int(length / grid.cell + 1e-9)) # 곧은 · 대각 1, 나이트 2 - for walked in range(layers): # 곧게 - heads.append(state(a, out, walked)) - tails.append(state(b, out, min(walked + units, full))) - costs.append(move) - for came in ((out - 1) % count, (out + 1) % count): # 이웃 방향에서 꺾어 들어옴 - bend = abs((angles[out] - angles[came] + math.pi) % (2.0 * math.pi) - math.pi) - heads.append(state(a, came, full)) - tails.append(state(b, out, min(units + 1, full))) - costs.append(move + TURN_COST * min_radius_m * bend) - size = grid.cells * count * layers - graph = csr_matrix( - (np.concatenate(costs), (np.concatenate(heads), np.concatenate(tails))), - shape=(size, size), - ) - start, end = grid.nearest(expected[0]), grid.nearest(expected[-1]) - distance, predecessors, _ = dijkstra( - graph, - indices=[state(start, d, 0) for d in range(count)], - return_predecessors=True, - min_only=True, - ) - straight = [state(end, d, full) for d in range(count)] - final = min(straight, key=lambda s: distance[s]) - if not np.isfinite(distance[final]): - final = min( - (state(end, d, w) for d in range(count) for w in range(layers)), - key=lambda s: distance[s], + def trace(nodes: list[tuple[float, float]], bound: float | None) -> tuple[Any, dict]: + # 노드마다 곡선(묶지 않음) — [확인] 재구성(`/route/replan`)과 같은 규칙 → 적용 뒤도 같은 선 + outline = build_planned_polyline( + nodes, min_radius_m=min_radius_m, simplify=False, curve_flags=[True] * len(nodes) ) - if not np.isfinite(distance[final]): - raise NoRoute("복도 안에서 제약을 지키며 시점과 종점을 잇는 길이 없습니다.") - chain = [int(final)] - while predecessors[chain[-1]] >= 0: - chain.append(int(predecessors[chain[-1]])) - chain.reverse() - # 꺾임점 = 방향이 바뀌는 칸. 첫 · 끝은 예상노선 끝점 그대로(기점 · 종점은 옮기지 않음). - nodes = [tuple(expected[0])] - for previous, current in zip(chain[1:], chain[2:]): - if (previous // layers) % count != (current // layers) % count: - cell = previous // layers // count - nodes.append((float(grid.centres[cell, 0]), float(grid.centres[cell, 1]))) - nodes.append(tuple(expected[-1])) - return nodes + return outline, measure(outline, bound) + + found: list[tuple[Any, dict[str, Any]]] = [] + for weights in CANDIDATE_WEIGHTS: + for step in GRADE_CAP_STEPS: + try: + nodes = graph.solve(weights, cap=max_grade * step, strict_stream=True) + except NoRoute: + break + found.append(trace(nodes, None)) + if found[-1][1]["over_limit"]["total"] == 0: + break + if not found: # 제약을 지키는 길 없음 — 비중과 무관하니 더 볼 것 없음 + break + if found: + least = min(item[1]["over_limit"]["total"] for item in found) + return [item for item in found if item[1]["over_limit"]["total"] == least], None + + # 지형상 불가능 — 반지름 · 계류 규칙을 지키는 어떤 길도 넘어야 하는 가장 작은 기울기를 + # 구해(`terrain_bound`) 그 위는 막고, 상한과 그 사이는 벌점으로 넘긴다. + bound = graph.terrain_bound() + cap = bound + 1e-9 + seed = trace( + graph.solve(DEFAULT_WEIGHTS, cap=cap, strict_stream=False, grade_cost=SOFT_GRADE_COST), + bound, + ) + found = [seed] + for weights in CANDIDATE_WEIGHTS: + for cost in RELAXED_GRADE_COSTS: + nodes = graph.solve(weights, cap=cap, strict_stream=False, grade_cost=cost) + found.append(trace(nodes, bound)) + most = seed[1]["over_limit"]["total"] + return [item for item in found if item[1]["over_limit"]["total"] <= most], bound + + +def _pool_key( + expected: list[tuple[float, float]], + terrain: Terrain, + streams: list[list[tuple[float, float]]], + *numbers: float, +) -> bytes: + """후보 묶음을 가를 입력 지문 — 지형 · 예상노선 · 계류 · 제약 값.""" + digest = hashlib.blake2b(digest_size=16) + for array in (terrain.x, terrain.y, terrain.z, terrain.valid, np.asarray(expected)): + digest.update(np.ascontiguousarray(array).tobytes()) + digest.update(repr((streams, numbers)).encode()) + return digest.digest() + + +def _pick( + mode: str, weights: dict[str, float], candidates: list[tuple[Any, dict[str, Any]]] +) -> tuple[Any, dict[str, Any]]: + """갈래 목표 비교값이 가장 좋은 후보. 같으면 앞 후보.""" + + def earth(metrics: dict[str, Any]) -> float: + return metrics["cut_m3"] + metrics["fill_m3"] + + def grade(metrics: dict[str, Any]) -> float: + return metrics["max_grade_pct"] + metrics["avg_grade_pct"] + + def offset(metrics: dict[str, Any]) -> float: + return metrics["offset_avg_m"] + + if mode == "earthwork": + return min(candidates, key=lambda item: earth(item[1])) + if mode == "gentle": + return min( + candidates, key=lambda item: (item[1]["max_grade_pct"], item[1]["avg_grade_pct"]) + ) + # 비중 — 세 값을 후보 범위로 0~1 로 맞춘 뒤 비중으로 더함(단위가 다른 값을 겨루게) + parts = [] + for key, value in (("expected", offset), ("earthwork", earth), ("grade", grade)): + values = [value(item[1]) for item in candidates] + low, span = min(values), max(values) - min(values) + parts.append([weights[key] * ((v - low) / span if span > 1e-9 else 0.0) for v in values]) + scores = [sum(column) for column in zip(*parts)] + return candidates[scores.index(min(scores))] diff --git a/B05_Profile/B05_Profile_Engine_RouteInitial_Metrics.py b/B05_Profile/B05_Profile_Engine_RouteInitial_Metrics.py index ddeb1fa2..bc5070b0 100644 --- a/B05_Profile/B05_Profile_Engine_RouteInitial_Metrics.py +++ b/B05_Profile/B05_Profile_Engine_RouteInitial_Metrics.py @@ -30,6 +30,8 @@ REASON_TERRAIN = "지형상 불가능 — 복도 안 어떤 길도 {bound:.1f}% REASON_GRADE_FIT = "곡선 맞춤 뒤 지반 차 — 측점 사이 지반이 상한을 넘음" REASON_RADIUS = "곡선 둘 자리 모자람 — 앞뒤 직선이 짧아 최소 반지름 원호가 안 들어감" REASON_STREAM = "계류 버퍼 안 — 건너기 · 기점 · 종점 봐줌을 넘음" +REASON_FOLLOW_GRADE = "예상노선 추적 — 예상노선이 지나는 지반이 상한을 넘음(제약 탐색 안 함)" +REASON_FOLLOW_STREAM = "예상노선 추적 — 예상노선이 계류 버퍼 안을 지남(제약 탐색 안 함)" REASON_STREAM_TERRAIN = "지형상 불가능 — 기울기 · 반지름을 지키며 계류 버퍼를 비켜 잇는 길 없음" @@ -110,10 +112,12 @@ def route_metrics( road_width_m: float, bound_pct: float | None, grade_smooth_m: float = GRADE_SMOOTH_M, + follow: bool = False, ) -> dict[str, Any]: """비교값과 제약 못 지킨 구간 목록. 종단기울기 = 측점(20 m) 사이 고른 지반 기울기. `bound_pct` 는 제약을 다 지키는 길이 없을 때만 — 복도 안 어떤 길도 넘어야 하는 기울기(%). + `follow` 는 예상노선 추적(제약 탐색 없이 예상노선 그대로) — 사유가 그 뜻으로 뜬다. """ grid = (terrain.y, terrain.x) sampler = RegularGridInterpolator(grid, terrain.z, bounds_error=False, fill_value=None) @@ -145,7 +149,7 @@ def route_metrics( float(ends[last]), float(grades[first : last + 1].max()) * 100.0, max_grade * 100.0, - terrain_reason if relaxed else REASON_GRADE_FIT, + REASON_FOLLOW_GRADE if follow else terrain_reason if relaxed else REASON_GRADE_FIT, ) ) @@ -181,7 +185,11 @@ def route_metrics( float(chain[last]), float(near[first : last + 1].min()), stream_offset_m, - REASON_STREAM_TERRAIN if relaxed else REASON_STREAM, + REASON_FOLLOW_STREAM + if follow + else REASON_STREAM_TERRAIN + if relaxed + else REASON_STREAM, ) ) diff --git a/B05_Profile/B05_Profile_Engine_RouteInitial_Search.py b/B05_Profile/B05_Profile_Engine_RouteInitial_Search.py new file mode 100644 index 00000000..53e3d625 --- /dev/null +++ b/B05_Profile/B05_Profile_Engine_RouteInitial_Search.py @@ -0,0 +1,361 @@ +"""초기 계획노선 갈래의 격자 탐색(PLAN 20 · 23장) — `B05_Profile_Engine_RouteInitial` 이 부름. + +탐색 = 예상노선 둘레 복도(`CORRIDOR_M`) 안 `SEARCH_CELL_M` 격자 · 16 방향 Dijkstra(scipy csgraph). +상태 = (칸, 진행 방향, 꺾은 뒤 곧게 온 걸음 수). 탐색 경로의 꺾임점이 곧 노드다. + +그래프(`RouteGraph`)는 요청마다 **한 번만** 짠다 — 복도 안 모든 링크를 상태 간선으로 펼쳐 두고, +탐색 상한 · 목표 비중 · 벌점 크기가 바뀌는 탐색은 간선 비용만 다시 계산한다(지형상 불가능 +상한 이분 탐색 · 후보 여러 벌을 한 요청 안에서). +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass + +import numpy as np +from scipy import ndimage +from scipy.sparse import csr_matrix +from scipy.sparse.csgraph import breadth_first_order, connected_components, dijkstra + +from B05_Profile.B05_Profile_Engine_RouteInitial_Metrics import GRADE_SMOOTH_M, smoothed_ground + +# 토공 가늠 폭(m) = 유효너비 3 + 길어깨 0.5 × 2(별표2 간선). +ROAD_WIDTH_M = 4.0 +# 탐색 복도 반폭(m) — 예상노선에서 이보다 멀리는 가지 않음(법령 값 아님 · 알고리즘 값). +CORRIDOR_M = 100.0 +SEARCH_CELL_M = 4.0 # 탐색 격자 칸(m) — 비용면(2 m)을 솎아 씀 +LENGTH_COST = 0.05 # 1 m 당 거리 비용 — 목표가 평평한 곳에서 헤매지 않게 +TURN_COST = 1.0 # 꺾을 때 벌점 = 최소 곡선반지름 원호 길이(R·θ) × 이 값 — 곧게 달리기 선호 +SOFT_STREAM_COST = 20.0 # 계류 버퍼 안 나란히 달리기(지형상 불가능할 때만 · 1 m 당) +BOUND_BISECT = 7 # 지형상 불가능할 때 막는 기울기를 찾는 이분 탐색 횟수(구간 1/128 까지) +STREAM_CROSS_COS = math.cos(math.radians(45.0)) # 버퍼 안 이동 방향 한도(계류 쪽 · 반대쪽 ±45°) +OBJECTIVES = ("expected", "earthwork", "grade") +# 16 방향(행, 열) — +x 에서 반시계. 나이트 걸음(1, 2) 등을 넣어 22.5° 안팎 방향을 곧게 달림. +DIRECTIONS = ( + (0, 1), + (1, 2), + (1, 1), + (2, 1), + (1, 0), + (2, -1), + (1, -1), + (1, -2), + (0, -1), + (-1, -2), + (-1, -1), + (-2, -1), + (-1, 0), + (-2, 1), + (-1, 1), + (-1, 2), +) + + +@dataclass +class Terrain: + """지반 격자 — 탐색 비용면(`B05_Profile_Engine_Solver._load_or_build_cost_surface`) 그대로.""" + + x: np.ndarray # 열 좌표(오름차순) + y: np.ndarray # 행 좌표(오름차순) + z: np.ndarray # [행, 열] 표고 + valid: np.ndarray # [행, 열] 지형 있음 + + +class NoRoute(ValueError): + """제약을 지키며 시점과 종점을 잇는 길이 없음.""" + + +def _cell_size(coords: np.ndarray) -> float: + return float((coords[-1] - coords[0]) / (len(coords) - 1)) if len(coords) > 1 else 1.0 + + +def _distance_grid( + lines: list[list[tuple[float, float]]], xs: np.ndarray, ys: np.ndarray, cell: float +) -> np.ndarray: + """격자 칸마다 선까지 거리(m). 선이 격자에 안 닿으면 모두 무한.""" + marked = np.zeros((len(ys), len(xs)), dtype=bool) + for line in lines: + for (x0, y0), (x1, y1) in zip(line, line[1:]): + count = max(2, int(math.hypot(x1 - x0, y1 - y0) / (cell * 0.5)) + 1) + px = np.linspace(x0, x1, count) + py = np.linspace(y0, y1, count) + cols = np.rint((px - xs[0]) / cell).astype(int) + rows = np.rint((py - ys[0]) / cell).astype(int) + inside = (cols >= 0) & (cols < len(xs)) & (rows >= 0) & (rows < len(ys)) + marked[rows[inside], cols[inside]] = True + if not marked.any(): + return np.full(marked.shape, np.inf) + return ndimage.distance_transform_edt(~marked) * cell + + +class SearchGrid: + """탐색 격자 — 복도 · 거리 · 기울기를 한 번만 만든다.""" + + def __init__( + self, + expected: list[tuple[float, float]], + terrain: Terrain, + streams: list[list[tuple[float, float]]], + *, + stream_offset_m: float, + grade_smooth_m: float = GRADE_SMOOTH_M, + ) -> None: + step = max(1, int(round(SEARCH_CELL_M / _cell_size(terrain.x)))) + self.xs, self.ys = terrain.x[::step], terrain.y[::step] + # 기울기 · 옆경사는 비교값과 같은 고른 지반에서 잰다(`smoothed_ground`) + self.z = smoothed_ground(terrain, grade_smooth_m)[::step, ::step] + self.stream_offset_m = stream_offset_m + self.cell = _cell_size(self.xs) + self.d_expected = _distance_grid([expected], self.xs, self.ys, self.cell) + self.d_stream = _distance_grid(streams, self.xs, self.ys, self.cell) + valid = np.asarray(terrain.valid[::step, ::step], dtype=bool) + self.corridor = (self.d_expected <= CORRIDOR_M) & valid + self.grad_y, self.grad_x = np.gradient(self.z, self.cell) + finite = np.where(np.isfinite(self.d_stream), self.d_stream, 0.0) + self.away_y, self.away_x = np.gradient(finite, self.cell) # 계류에서 멀어지는 쪽 + self.index = np.full(self.z.shape, -1, dtype=np.int64) + self.cells = int(self.corridor.sum()) + self.index[self.corridor] = np.arange(self.cells) + rows, cols = np.nonzero(self.corridor) + self.centres = np.column_stack([self.xs[cols], self.ys[rows]]) + + def nearest(self, point: tuple[float, float]) -> int: + dx = self.centres[:, 0] - point[0] + return int(np.argmin(np.hypot(dx, self.centres[:, 1] - point[1]))) + + def links(self, dr: int, dc: int, *, limit: float) -> tuple: + """복도 안 한 방향 링크 모두 — (앞 번호, 뒤 번호, 길이, 세 목표[3, n] · 1 m 당, + 넘은 기울기[n], 계류 버퍼 안 나란히[n], 지반 기울기[n]).""" + rows, cols = self.z.shape + first = (slice(max(0, -dr), rows - max(0, dr)), slice(max(0, -dc), cols - max(0, dc))) + second = (slice(max(0, dr), rows - max(0, -dr)), slice(max(0, dc), cols - max(0, -dc))) + a, b = self.index[first], self.index[second] + ok = (a >= 0) & (b >= 0) + norm = math.hypot(dr, dc) + grade = (np.abs(self.z[second] - self.z[first]) / (self.cell * norm))[ok] + # 계류 버퍼 안에서는 계류 쪽 · 반대쪽 ±45° 로만(건너기 · 빠져나가기) + in_buffer = (self.d_stream[first] < self.stream_offset_m) | ( + self.d_stream[second] < self.stream_offset_m + ) + ax = self.away_x[first] + self.away_x[second] + ay = self.away_y[first] + self.away_y[second] + size = np.hypot(ax, ay) + across = np.abs(ax * dc + ay * dr) >= STREAM_CROSS_COS * norm * size + along = (in_buffer & ~across & (size > 1e-6))[ok] + over = np.maximum(0.0, grade - limit) + # 옆경사 = 진행 방향에 수직인 지반 기울기(두 칸 평균) + cross = np.abs( + (self.grad_x[first][ok] + self.grad_x[second][ok]) * -dr + + (self.grad_y[first][ok] + self.grad_y[second][ok]) * dc + ) / (2.0 * norm) + # 토공 가늠(1 m 당 m³) — 옆경사 위 반절 · 반성 단면 W²·s/4 + 상한 넘는 기울기가 측점 + # 반 간격(10 m) 동안 쌓는 높이 차 단면(지형상 불가능할 때만 0 이 아님) + earth = ROAD_WIDTH_M**2 * cross / 4.0 + ROAD_WIDTH_M * over * 10.0 + offset = 0.5 * (self.d_expected[first][ok] + self.d_expected[second][ok]) + terms = np.stack([offset, earth, (grade / limit) ** 2]) + return a[ok], b[ok], self.cell * norm, terms, over, along, grade + + +def _turn_steps(min_radius_m: float, cell: float) -> int: + """꺾은 뒤 곧게 가야 할 칸 수 — 앞뒤 원호 접선이 한 직선에 들어가게. + + 원호를 끼우는 쪽(`build_planned_polyline`)은 직선의 절반까지만 접선으로 쓰므로 직선은 + 2 · R · tan(θ/2) 이상이어야 한다. θ = 이웃 방향 사이 가장 큰 각(나이트 걸음 26.57°). + 걸음은 칸 길이 단위로 센다 — 곧은 · 대각 걸음 1, 나이트 걸음 2. + """ + widest = math.atan2(1.0, 2.0) + return max(1, math.ceil(2.0 * min_radius_m * math.tan(widest / 2.0) / cell)) + + +class RouteGraph: + """상태 그래프 — 복도 안 모든 링크를 한 번 펼쳐 두고 탐색마다 비용 · 가릴 간선만 바꾼다. + + 상태 = (칸, 방향, 층). 층 = 꺾은 뒤 곧게 온 칸 수 + 1(K + 1 에서 멈춤) · 층 0 = 시점. + 이웃 방향으로 꺾기는 K 칸 곧게 온 뒤에만(층 K + 1) → 노드마다 최소 반지름 원호가 들어감. + 시점은 한 칸 더(K + 1 칸) — 기점을 격자 칸 가운데에서 예상노선 끝점으로 옮기며 첫 직선이 + 짧아지는 몫. 종점은 K 칸 곧게 들어온 상태를 먼저, 없으면 아무 상태(마지막 곡선 자리가 + 모자라면 목록에 뜸). + """ + + def __init__( + self, + grid: SearchGrid, + expected: list[tuple[float, float]], + *, + limit: float, + min_radius_m: float, + ) -> None: + self.grid, self.expected, self.limit = grid, expected, limit + count = len(DIRECTIONS) + self.layers = _turn_steps(min_radius_m, grid.cell) + 2 # 꺾을 수 있는 층 = layers - 1 + full = self.layers - 1 + angles = [math.atan2(dr, dc) for dr, dc in DIRECTIONS] + heads, tails, link_of, turn = [], [], [], [] + per_link = [] + offset = 0 + for out, (dr, dc) in enumerate(DIRECTIONS): + a, b, length, terms, over, along, grade = grid.links(dr, dc, limit=limit) + ids = np.arange(offset, offset + a.size) + offset += a.size + per_link.append((np.full(a.size, length), terms, over, along, grade)) + units = max(1, int(length / grid.cell + 1e-9)) # 곧은 · 대각 1, 나이트 2 + for walked in range(self.layers): # 곧게 + heads.append(self._state(a, out, walked)) + tails.append(self._state(b, out, min(walked + units, full))) + link_of.append(ids) + turn.append(np.zeros(a.size)) + for came in ((out - 1) % count, (out + 1) % count): # 이웃 방향에서 꺾어 들어옴 + bend = abs((angles[out] - angles[came] + math.pi) % (2.0 * math.pi) - math.pi) + heads.append(self._state(a, came, full)) + tails.append(self._state(b, out, min(units + 1, full))) + link_of.append(ids) + turn.append(np.full(a.size, TURN_COST * min_radius_m * bend)) + if offset == 0: + raise NoRoute("예상노선 둘레에 지날 수 있는 링크가 없습니다.") + self.length = np.concatenate([p[0] for p in per_link]) + self.terms = np.hstack([p[1] for p in per_link]) + self.over = np.concatenate([p[2] for p in per_link]) + self.along = np.concatenate([p[3] for p in per_link]) + self.grade = np.concatenate([p[4] for p in per_link]) + self.size = grid.cells * count * self.layers + # 간선 틀(CSR)은 한 번만 — 탐색마다 못 쓰는 간선을 걸러 내고 비용만 갈아 끼움 + link_of = np.concatenate(link_of) + order = csr_matrix( + (np.arange(1.0, link_of.size + 1.0), (np.concatenate(heads), np.concatenate(tails))), + shape=(self.size, self.size), + ) + self.indices, self.indptr = order.indices, order.indptr + edge = (order.data - 1.0).astype(np.int64) + self.edge_link = link_of[edge] + self.edge_turn = np.concatenate(turn)[edge] + start, end = grid.nearest(expected[0]), grid.nearest(expected[-1]) + self.starts = [self._state(start, d, 0) for d in range(count)] + self.straight = [self._state(end, d, full) for d in range(count)] + self.ends = [self._state(end, d, w) for d in range(count) for w in range(self.layers)] + + def _state(self, cell, direction: int, walked: int): + return (cell * len(DIRECTIONS) + direction) * self.layers + walked + + def _usable(self, cap: float, strict_stream: bool) -> np.ndarray: + """링크별 쓸 수 있음 — `cap` 넘는 기울기 · (엄격하면) 계류 버퍼 안 나란히 빼고.""" + usable = self.grade <= cap + return usable & ~self.along if strict_stream else usable + + def _graph(self, usable: np.ndarray, link_cost: np.ndarray | None) -> csr_matrix: + """쓸 수 있는 간선만 남긴 그래프 + 시점 상태 모두로 가는 원점(마지막 상태 `size`).""" + keep = usable[self.edge_link] + indptr = np.concatenate([[0], np.cumsum(keep)])[self.indptr] + indptr = np.append(indptr, indptr[-1] + len(self.starts)) + indices = np.concatenate([self.indices[keep], self.starts]) + if link_cost is None: + data = np.ones(indices.size) + else: + cost = link_cost[self.edge_link[keep]] + self.edge_turn[keep] + data = np.concatenate([cost, np.full(len(self.starts), 1e-9)]) + return csr_matrix((data, indices, indptr), shape=(self.size + 1, self.size + 1)) + + def joined(self, cap: float) -> bool: + """`cap` 까지 기울기를 풀면 시점 · 종점이 이어지나(계류는 벌점 — 기울기만 따짐).""" + graph = self._graph(self._usable(cap, False), None) + reached = breadth_first_order(graph, self.size, return_predecessors=False) + return bool(np.isin(self.ends, reached).any()) + + def solve( + self, + weights: dict[str, float], + *, + cap: float, + strict_stream: bool, + grade_cost: float = 0.0, + ) -> list[tuple[float, float]]: + """노드(꺾임점) 목록 — 첫 · 끝은 예상노선 양 끝 그대로. + + 세 목표는 쓸 수 있는 링크 평균이 1 이 되게 정규화(길이 가중)한 뒤 비중을 곱한다. + `grade_cost` = 상한을 넘은 기울기 벌점(1 m 당 × 넘은 비율 / 상한) — `cap` 이 상한보다 + 클 때(지형상 불가능)만 뜻이 있다. + """ + usable = self._usable(cap, strict_stream) + if not usable.any(): + raise NoRoute("예상노선 둘레에 제약을 지키며 지날 수 있는 링크가 없습니다.") + length = self.length[usable] + scale = (self.terms[:, usable] * length).sum(axis=1) / max(length.sum(), 1e-9) + share = np.array([weights[key] for key in OBJECTIVES]) / max(sum(weights.values()), 1e-9) + factor = np.divide(share, scale, out=np.zeros_like(share), where=scale > 1e-12) + penalty = grade_cost * self.over / self.limit + if not strict_stream: + penalty = penalty + SOFT_STREAM_COST * self.along + link_cost = self.length * (LENGTH_COST + factor @ self.terms + penalty) + distance, predecessors = dijkstra( + self._graph(usable, link_cost), indices=self.size, return_predecessors=True + ) + final = min(self.straight, key=lambda s: distance[s]) + if not np.isfinite(distance[final]): + final = min(self.ends, key=lambda s: distance[s]) + if not np.isfinite(distance[final]): + raise NoRoute("복도 안에서 제약을 지키며 시점과 종점을 잇는 길이 없습니다.") + chain = [int(final)] + while predecessors[chain[-1]] >= 0: + chain.append(int(predecessors[chain[-1]])) + chain = chain[::-1][1:] # 원점 뺌 + # 꺾임점 = 방향이 바뀌는 칸. 첫 · 끝은 예상노선 끝점 그대로(기점 · 종점은 옮기지 않음). + count, layers = len(DIRECTIONS), self.layers + nodes = [tuple(self.expected[0])] + for previous, current in zip(chain[1:], chain[2:]): + if (previous // layers) % count != (current // layers) % count: + cell = previous // layers // count + nodes.append((float(self.grid.centres[cell, 0]), float(self.grid.centres[cell, 1]))) + nodes.append(tuple(self.expected[-1])) + return nodes + + def minimax_grade(self) -> float: + """복도 안에서 시점과 종점을 잇는 데 **어떤 길도 넘어야 하는** 가장 작은 기울기 상한(비율). + + 꺾기 규칙 없이 칸끼리만 이은 그래프에서 이분 탐색 — 그러니 실제 필요한 값의 아래 한계다. + """ + grid = self.grid + pieces = [] + for dr, dc in DIRECTIONS: + a, b, *_, grade = grid.links(dr, dc, limit=self.limit) + pieces.append((a, b, grade)) + heads = np.concatenate([p[0] for p in pieces]) + tails = np.concatenate([p[1] for p in pieces]) + grades = np.concatenate([p[2] for p in pieces]) + start, end = grid.nearest(self.expected[0]), grid.nearest(self.expected[-1]) + + def joined(cap: float) -> bool: + keep = grades <= cap + graph = csr_matrix( + (np.ones(int(keep.sum())), (heads[keep], tails[keep])), + shape=(grid.cells, grid.cells), + ) + _, labels = connected_components(graph, directed=False) + return bool(labels[start] == labels[end]) + + low, high = 0.0, float(grades.max()) if grades.size else 0.0 + if not joined(high): + raise NoRoute("복도 안에서 시점과 종점이 이어지지 않습니다.") + for _ in range(12): # 아래 한계일 뿐 — 구간 1/4096 이면 넉넉함 + middle = 0.5 * (low + high) + low, high = (low, middle) if joined(middle) else (middle, high) + return high + + def terrain_bound(self) -> float: + """반지름 규칙까지 지키며 시점 · 종점을 잇는 데 **어떤 길도 넘어야 하는** 기울기(비율). + + 아래 한계 = 꺾기 규칙 없는 최소최대(`minimax_grade`). 이어지는 값을 1.5 배씩 찾은 뒤 + 그 사이를 상태 그래프 그대로 이분 탐색한다(계류는 벌점 — 기울기만 따짐). + """ + ceiling = float(self.grade.max(initial=0.0)) + low = max(self.minimax_grade(), self.limit) + high = low + while not self.joined(high): + if high >= ceiling: + raise NoRoute("복도 안에서 반지름 규칙을 지키며 시점과 종점을 잇는 길이 없습니다.") + low, high = high, min(high * 1.5, ceiling) + for _ in range(BOUND_BISECT): + if high - low <= 1e-4: + break + middle = 0.5 * (low + high) + low, high = (low, middle) if self.joined(middle) else (middle, high) + return high diff --git a/B05_Profile/B05_Profile_Router_Replan.py b/B05_Profile/B05_Profile_Router_Replan.py index 42edf351..07ffca4b 100644 --- a/B05_Profile/B05_Profile_Router_Replan.py +++ b/B05_Profile/B05_Profile_Router_Replan.py @@ -626,87 +626,5 @@ async def reset_route_plan(project_id: UUID) -> dict[str, Any] | JSONResponse: return {"status": "success", "project_id": str(project_id), **result} -@router.post("/{project_id}/route/initial", response_model=None) -async def initial_route( - project_id: UUID, request: RouteInitialRequest -) -> dict[str, Any] | JSONResponse: - """초기 계획노선 갈래 하나를 만들어 비교값과 함께 돌려준다(PLAN 20장). - - **파일을 쓰지 않는다** — 화면 캐시용. 적용은 받은 `nodes` 로 `POST /route/replan`([확인]). - """ - from B05_Profile.B05_Profile_Engine_Grade import legal_grade_limit_pct - from B05_Profile.B05_Profile_Engine_RouteInitial import ( - generate_initial_route, - load_terrain_and_streams, - ) - from common_util.common_util_surface_confirmation import get_surface_confirmation_params - - weights = request.weights.model_dump() if request.weights else None - if request.mode == "weighted" and weights and abs(sum(weights.values()) - 100.0) > 0.5: - return JSONResponse( - status_code=400, - content={"status": "error", "message": "비중 세 값의 합은 100 이어야 합니다."}, - ) - paths = await _project_paths(project_id) - if paths is None: - return JSONResponse(status_code=404, content=_PROJECT_PATH_MISSING) - project_root, _ = paths - expected = await asyncio.to_thread(_vertices_of, expected_route_csv_path(project_root)) - if not expected: - expected = await asyncio.to_thread(_vertices_of, design_route_csv_path(project_root)) - if len(expected) < 2: - return JSONResponse( - status_code=400, - content={"status": "error", "message": "예상노선이 없어 초기노선을 만들 수 없습니다."}, - ) - radius_m, _, _ = await _plan_criteria(project_id) - grade_class, design_speed, terrain_type = await _road_settings(project_id) - # 비포장 기준 — 포장 예외(18 %)는 포장 여부가 정해진 뒤 종단 설계가 따진다(사용자 확정 대기). - max_grade = legal_grade_limit_pct(grade_class, terrain_type, False, design_speed) / 100.0 - # 화면 제약 칸(PLAN 20장 확정) — 비운 칸은 위 프로젝트 기본값 그대로, 채운 칸만 겹쳐 쓴다. - criteria = request.criteria - if criteria and criteria.max_grade_pct is not None: - max_grade = criteria.max_grade_pct / 100.0 - if criteria and criteria.min_radius_m is not None: - radius_m = criteria.min_radius_m - stream_offset_m = criteria.stream_offset_m if criteria else None - grade_smooth_m = criteria.grade_smooth_m if criteria else None - pool = get_db_pool() - async with pool.acquire() as connection: - selection = await get_surface_confirmation_params(connection, str(project_id)) - started = time.perf_counter() - try: - terrain, streams = await asyncio.to_thread( - load_terrain_and_streams, project_root, selection - ) - result = await asyncio.to_thread( - functools.partial( - generate_initial_route, - [(x, y) for x, y in expected], - terrain, - streams, - mode=request.mode, - weights=weights, - max_grade=max_grade, - min_radius_m=radius_m, - stream_offset_m=stream_offset_m, - grade_smooth_m=grade_smooth_m, - ) - ) - except (FileNotFoundError, ValueError) as error: - return JSONResponse(status_code=409, content={"status": "error", "message": str(error)}) - logger.info( - "초기 계획노선 %s: project_id=%s %.1fs 연장 %.0fm 이격 평균 %.1fm", - request.mode, - project_id, - time.perf_counter() - started, - result["metrics"]["length_m"], - result["metrics"]["offset_avg_m"], - ) - current = await asyncio.to_thread(read_initial_choice, project_root) - return { - "status": "success", - "project_id": str(project_id), - **result, - "current_initial": current, - } +# `POST /route/initial` 은 나눈 파일에서 같은 `router` 에 붙는다 — 끝에서 불러 등록. +from B05_Profile import B05_Profile_Router_RouteInitial # noqa: E402, F401 diff --git a/B05_Profile/B05_Profile_Router_RouteInitial.py b/B05_Profile/B05_Profile_Router_RouteInitial.py new file mode 100644 index 00000000..b0639479 --- /dev/null +++ b/B05_Profile/B05_Profile_Router_RouteInitial.py @@ -0,0 +1,111 @@ +"""초기 계획노선 갈래 하나 — `POST /api/projects/{id}/route/initial`(PLAN 20 · 23장). + +`B05_Profile_Router_Replan`(700줄)에서 나눈 길. 같은 `router` 에 붙고 그 모듈 끝에서 불려 +등록된다 — main · 시험이 `Router_Replan.router` 하나만 실어도 이 길이 함께 실림. +""" + +import asyncio +import functools +import logging +import time +from typing import Any +from uuid import UUID + +from fastapi.responses import JSONResponse + +from B05_Profile import B05_Profile_Router_Replan as replan +from B05_Profile.B05_Profile_Engine_RouteInitial import read_initial_choice +from common_util.common_util_initial_snapshot import design_route_csv_path +from common_util.common_util_route_geometry import expected_route_csv_path +from config.config_db import get_db_pool + +logger = logging.getLogger(__name__) + + +@replan.router.post("/{project_id}/route/initial", response_model=None) +async def initial_route( + project_id: UUID, request: replan.RouteInitialRequest +) -> dict[str, Any] | JSONResponse: + """초기 계획노선 갈래 하나를 만들어 비교값과 함께 돌려준다(PLAN 20장). + + **파일을 쓰지 않는다** — 화면 캐시용. 적용은 받은 `nodes` 로 `POST /route/replan`([확인]). + """ + from B05_Profile.B05_Profile_Engine_Grade import legal_grade_limit_pct + from B05_Profile.B05_Profile_Engine_RouteInitial import ( + generate_initial_route, + load_terrain_and_streams, + ) + from common_util.common_util_surface_confirmation import get_surface_confirmation_params + + weights = request.weights.model_dump() if request.weights else None + if request.mode == "weighted" and weights and abs(sum(weights.values()) - 100.0) > 0.5: + return JSONResponse( + status_code=400, + content={"status": "error", "message": "비중 세 값의 합은 100 이어야 합니다."}, + ) + try: + paths = await replan._project_paths(project_id) + except LookupError: # 없는 프로젝트 — 저장소 조회가 None 대신 던짐 + paths = None + if paths is None: + return JSONResponse(status_code=404, content=replan._PROJECT_PATH_MISSING) + project_root, _ = paths + expected = await asyncio.to_thread(replan._vertices_of, expected_route_csv_path(project_root)) + if not expected: + expected = await asyncio.to_thread(replan._vertices_of, design_route_csv_path(project_root)) + if len(expected) < 2: + return JSONResponse( + status_code=400, + content={"status": "error", "message": "예상노선이 없어 초기노선을 만들 수 없습니다."}, + ) + radius_m, _, _ = await replan._plan_criteria(project_id) + grade_class, design_speed, terrain_type = await replan._road_settings(project_id) + # 비포장 기준 — 포장 예외(18 %)는 포장 여부가 정해진 뒤 종단 설계가 따진다(사용자 확정 대기). + max_grade = legal_grade_limit_pct(grade_class, terrain_type, False, design_speed) / 100.0 + # 화면 제약 칸(PLAN 20장 확정) — 비운 칸은 위 프로젝트 기본값 그대로, 채운 칸만 겹쳐 쓴다. + criteria = request.criteria + if criteria and criteria.max_grade_pct is not None: + max_grade = criteria.max_grade_pct / 100.0 + if criteria and criteria.min_radius_m is not None: + radius_m = criteria.min_radius_m + stream_offset_m = criteria.stream_offset_m if criteria else None + grade_smooth_m = criteria.grade_smooth_m if criteria else None + pool = get_db_pool() + async with pool.acquire() as connection: + selection = await get_surface_confirmation_params(connection, str(project_id)) + started = time.perf_counter() + try: + terrain, streams = await asyncio.to_thread( + load_terrain_and_streams, project_root, selection + ) + result = await asyncio.to_thread( + functools.partial( + generate_initial_route, + [(x, y) for x, y in expected], + terrain, + streams, + mode=request.mode, + weights=weights, + max_grade=max_grade, + min_radius_m=radius_m, + stream_offset_m=stream_offset_m, + grade_smooth_m=grade_smooth_m, + ) + ) + except (FileNotFoundError, ValueError) as error: + return JSONResponse(status_code=409, content={"status": "error", "message": str(error)}) + logger.info( + "초기 계획노선 %s: project_id=%s %.1fs 연장 %.0fm 이격 평균 %.1fm", + request.mode, + project_id, + time.perf_counter() - started, + result["metrics"]["length_m"], + result["metrics"]["offset_avg_m"], + ) + current = await asyncio.to_thread(read_initial_choice, project_root) + return { + "status": "success", + "project_id": str(project_id), + **result, + "current_initial": current, + } diff --git a/resources/tester/test_route_initial.py b/resources/tester/test_route_initial.py index 7fd4f980..b9bfe580 100644 --- a/resources/tester/test_route_initial.py +++ b/resources/tester/test_route_initial.py @@ -4,10 +4,15 @@ 공통 제약(종단기울기 상한 · 최소 곡선반지름 · 계류 이격)은 넘지 못함 — 지형상 불가능할 때만 `metrics.violations` 에 측점 · 값 · 사유로 뜨고 `over_limit.total` 은 그 수와 같다. +예상노선 추적(follow)은 전처리 초기값과 같은 선 — 제약 탐색 없이 넘은 곳만 목록(PLAN 23-1). +갈래는 제 목표 비교값이 가장 좋은 후보 — 실프로젝트(랩탑_보조 · 읽기만)에서 갈래가 갈림. """ +from pathlib import Path + import numpy as np import pytest +from shapely.geometry import LineString from B05_Profile.B05_Profile_Engine_RouteInitial import ( MODES, @@ -21,6 +26,9 @@ from B05_Profile.B05_Profile_Engine_RouteInitial import ( write_initial_choice, ) from B05_Profile.B05_Profile_Engine_RouteInitial_Metrics import grade_line, station_label +from common_util.common_util_route_polyline import build_planned_polyline + +SEARCHED = [mode for mode in MODES if mode != "follow"] GRADE = 0.14 RADIUS = 12.0 @@ -68,7 +76,7 @@ def hill() -> Terrain: ) -@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("mode", SEARCHED) def test_every_mode_keeps_grade_limit_on_hill(hill: Terrain, mode: str) -> None: result = _run(hill, mode) _kept(result) @@ -81,6 +89,30 @@ def test_every_mode_keeps_grade_limit_on_hill(hill: Terrain, mode: str) -> None: assert metrics["offset_max_m"] <= result["criteria"]["corridor_m"] + 4.0 +def test_follow_is_initial_polyline_and_lists_violations(hill: Terrain) -> None: + """추종 = 전처리 초기값 규칙(`build_planned_polyline` 기본값) — 언덕을 넘는 곳은 목록에.""" + result = _run(hill, "follow") + initial = build_planned_polyline(_line(200.0), min_radius_m=RADIUS) + assert result["planned"] == [[round(x, 4), round(y, 4)] for x, y in initial.vertices] + assert result["nodes"] == [node.as_dict() for node in initial.nodes] + assert result["relaxed"] is False and result["terrain_bound_pct"] is None + grades = [v for v in result["metrics"]["violations"] if v["kind"] == "grade"] + assert grades and all("예상노선 추적" in v["reason"] for v in grades) + assert result["metrics"]["over_limit"]["total"] == len(result["metrics"]["violations"]) + + +def _earth(result: dict) -> float: + return result["metrics"]["cut_m3"] + result["metrics"]["fill_m3"] + + +def test_each_branch_best_at_own_goal(hill: Terrain) -> None: + """토공 최소 = 절성토 가장 적음 · 기울기 순한 = 최대 기울기 가장 낮음(같은 후보 묶음).""" + runs = {mode: _run(hill, mode) for mode in SEARCHED} + assert all(_earth(runs["earthwork"]) <= _earth(r) + 1e-6 for r in runs.values()) + top = runs["gentle"]["metrics"]["max_grade_pct"] + assert all(top <= r["metrics"]["max_grade_pct"] + 1e-6 for r in runs.values()) + + def test_follow_stays_on_expected_where_allowed() -> None: """고른 경사(2 %) — 추종은 예상노선을 그대로(격자 칸 반 안쪽).""" ramp = _terrain(lambda xx, yy: 0.02 * xx) @@ -99,7 +131,7 @@ def test_stream_offset_is_hard() -> None: """평지 · 예상노선 20 m 옆을 나란히 흐르는 소하천 — 모든 갈래가 버퍼를 비켜 감.""" flat = _terrain(lambda xx, yy: np.zeros_like(xx)) stream = [[(0.0, 220.0), (800.0, 220.0)]] - for mode in MODES: + for mode in SEARCHED: result = _run(flat, mode, streams=stream) _kept(result) assert result["metrics"]["over_limit"]["stream"] == 0 @@ -237,3 +269,78 @@ def test_request_accepts_criteria_overrides() -> None: ) assert replan.initial is not None and replan.initial.criteria is not None assert replan.initial.criteria.max_grade_pct == 12 + + +# ── 실프로젝트(랩탑_보조) — 읽기만 · storage 가 없으면 건너뜀 ────────────────────────── +PROJECT = Path(__file__).resolve().parents[2] / "storage/1/3/936be972-11bc-46c2-8bf3-b15d8de7df0d" +SURFACE = PROJECT / "B05_Profile/route/cost_surface_classification_dtm_smooth.npz" + + +@pytest.fixture(scope="module") +def project_runs() -> dict: + """네 갈래 결과 — 비용면 캐시 · 계류 · 예상노선을 읽기만(파일 안 씀). 기준 14 % · R12.""" + if not SURFACE.is_file(): + pytest.skip("랩탑_보조 storage 없음") + from pyproj import Transformer + + from B04_PreProcess.B04_PreProcess_Router_Watershed import ( + STREAM_FILE, + _load_features, + _reproject_features, + ) + from common_util.common_util_crs import resolve_project_crs + from common_util.common_util_route_geometry import ( + expected_route_csv_path, + read_planned_route_csv, + ) + + cached = np.load(SURFACE, allow_pickle=False) + terrain = Terrain(cached["x"], cached["y"], cached["z"], cached["valid_mask"]) + to_metric = Transformer.from_crs("EPSG:4326", resolve_project_crs(PROJECT), always_xy=True) + features = _load_features(PROJECT / "B04_PreProcess" / "processed", STREAM_FILE) + streams = stream_lines(_reproject_features(features, to_metric)) + route = read_planned_route_csv(expected_route_csv_path(PROJECT)) + expected = [(float(v.x), float(v.y)) for v in route.vertices] + return { + mode: generate_initial_route( + expected, + terrain, + streams, + mode=mode, + weights=None, + max_grade=GRADE, + min_radius_m=RADIUS, + ) + for mode in MODES + } + + +def test_project_follow_is_preprocess_initial(project_runs: dict) -> None: + """추종 = 전처리 초기값 `planned_route_initial.csv`(1078 m) 그대로.""" + from common_util.common_util_route_geometry import ( + planned_route_initial_path, + read_planned_route_csv, + ) + + initial = read_planned_route_csv(planned_route_initial_path(PROJECT)) + stored = LineString([(v.x, v.y) for v in initial.vertices]) + follow = LineString(project_runs["follow"]["planned"]) + assert follow.hausdorff_distance(stored) < 0.01 + assert abs(follow.length - stored.length) < 0.01 + + +def test_project_branches_differ_by_goal(project_runs: dict) -> None: + """갈래끼리 선이 수 m 이상 갈리고 · 각 갈래가 제 목표 비교값에서 가장 좋음.""" + searched = {mode: project_runs[mode] for mode in SEARCHED} + for first in SEARCHED: + for second in SEARCHED: + if first < second: + gap = LineString(searched[first]["planned"]).hausdorff_distance( + LineString(searched[second]["planned"]) + ) + assert gap >= 3.0, (first, second, gap) + assert all(_earth(searched["earthwork"]) <= _earth(r) for r in searched.values()) + top = searched["gentle"]["metrics"]["max_grade_pct"] + assert all(top <= r["metrics"]["max_grade_pct"] for r in searched.values()) + for result in searched.values(): + assert result["metrics"]["over_limit"]["total"] == len(result["metrics"]["violations"])