#!/usr/bin/env python3
"""
generate_atlas.py — Terrain-aware sequential city placement and infrastructure
generation for The Settled Reach Atlas (Phase 3, D-191 §3, §8, §9, #832).

Pipeline per body:
  1. Load body definition from body index.md frontmatter
  2. Simulate terrain via planet_simulation.simulate()
  3. Analyse terrain: continents, habitability, river mouths, cost grid
  4. Place cities sequentially (capital first, corridor growth, quadrant spread)
  5. Generate infrastructure: A* roads + rail MST connecting all cities
  6. Place gate terminal POI at largest population centre (occasionally scatter)
  7. Write updated markers.json (preserves existing rivers/oceans/mountains)

City names are left as empty strings; gemma_naming.py (#833) fills them.

Usage:
    python3 tooling/planet-gen/generate_atlas.py
    python3 tooling/planet-gen/generate_atlas.py --body GJ380c
    python3 tooling/planet-gen/generate_atlas.py --force
    python3 tooling/planet-gen/generate_atlas.py --dry-run
    python3 tooling/planet-gen/generate_atlas.py --seed 12345

Decisions: D-191 (atlas pipeline), D-188 (planet_class field name)
"""

import argparse
import hashlib
import heapq
import json
import math
import os
import sys
import time
from pathlib import Path

# ---------------------------------------------------------------------------
# Venv bootstrap
# ---------------------------------------------------------------------------
TOOLING_DIR = Path(__file__).resolve().parent
REPO_ROOT = (TOOLING_DIR / ".." / "..").resolve()
_venv_python = REPO_ROOT / ".venv" / "bin" / "python"
if _venv_python.exists() and Path(sys.executable).resolve() != _venv_python.resolve():
    os.execv(str(_venv_python), [str(_venv_python)] + sys.argv)

import numpy as np
import sqlite3
import yaml

from planet_simulation import simulate

# ---------------------------------------------------------------------------
# Constants
# ---------------------------------------------------------------------------

GRID_W = 512
GRID_H = 256

# Terrain cost grid (D-191 §3: A* on terrain cost grid)
COST_WATER = 1e9    # impassable
COST_MOUNTAIN = 10.0  # expensive
COST_RIVER = 0.5    # cheap corridor
COST_FLAT = 1.0     # baseline land

# City placement
MOUNTAIN_ELEV_THRESHOLD = 0.55   # normalised elevation above sea_level → mountain

SETTLEMENT_PATTERN_MODIFIERS = {
    "urban_concentrated": -1,
    "dispersed_rural": +1,
    "orbital_only": None,   # no surface cities
    "domed": None,          # special: forced to 1
    "cave": None,           # special: forced to 1
}

# Quadrant distribution (D-191 §3: after 2 cities in same quadrant, prefer others)
QUADRANT_SATURATION = 2

DB_PATH = REPO_ROOT / "server" / "data" / "systems.db"
WIKI_SYSTEMS = REPO_ROOT / "wiki" / "star-systems"
SYSTEMS_SCHEMA_PATH = REPO_ROOT / "server" / "data" / "systems-schema.sql"

# Delimiters for the single canonical atlas_* DDL block in
# server/data/systems-schema.sql. `_load_atlas_schema()` extracts everything
# between these markers at runtime so this file does not have to duplicate
# the schema (and drift from it).
_ATLAS_SCHEMA_BEGIN_MARKER = "-- BEGIN ATLAS INDEX"
_ATLAS_SCHEMA_END_MARKER = "-- END ATLAS INDEX"


# ---------------------------------------------------------------------------
# Generator metadata stamp (#855, #856)
# ---------------------------------------------------------------------------

def _file_sha1(*paths: Path) -> str:
    """Return SHA-1 hex of the concatenated content of one or more files.

    Files are sorted by path for determinism.  Missing files raise
    FileNotFoundError rather than silently skip — a ghost hash (empty-bytes
    digest) can mask real breakage when stored and current SHAs converge
    (#136 review H2).
    """
    h = hashlib.sha1()
    for p in sorted(paths):
        if not p.exists():
            raise FileNotFoundError(f"generator source not found: {p}")
        h.update(p.read_bytes())
    return h.hexdigest()


def _write_stamp(conn: sqlite3.Connection) -> None:
    """Upsert a meta row for generate_atlas after a successful run.

    Idempotent: running twice on the same source files writes the same SHA
    with an updated timestamp.  The meta table is created by the atlas schema
    migration executed in ensure_atlas_schema(); this function assumes it
    exists.

    Transaction ownership stays with the caller (matches the import_economics
    pattern) — no inner commit here.  Review H1 flagged the prior behaviour as
    a double-commit with the atlas data write that precedes it.
    """
    schema_sha = _file_sha1(SYSTEMS_SCHEMA_PATH)
    generator_sha = _file_sha1(Path(__file__))
    conn.execute(
        """INSERT OR REPLACE INTO meta
               (generator_name, schema_version, generator_sha, generated_at)
           VALUES ('generate_atlas', ?, ?, datetime('now'))""",
        (schema_sha, generator_sha),
    )


def _load_atlas_schema() -> str:
    """Return the atlas_* DDL block from systems-schema.sql.

    systems-schema.sql is the single source of truth for the atlas index
    tables (see the BEGIN ATLAS INDEX / END ATLAS INDEX markers). We extract
    just that block and run it through `executescript` so the generator
    works against any DB state (fresh or partially-migrated) without
    requiring a prior `wiki_sync.ensure_schema()` call and without
    maintaining a second copy of the DDL here.
    """
    if not SYSTEMS_SCHEMA_PATH.exists():
        raise RuntimeError(
            f"systems-schema.sql not found at {SYSTEMS_SCHEMA_PATH} — "
            "atlas generator cannot proceed without the canonical schema."
        )
    text = SYSTEMS_SCHEMA_PATH.read_text()
    try:
        start = text.index(_ATLAS_SCHEMA_BEGIN_MARKER)
        end = text.index(_ATLAS_SCHEMA_END_MARKER, start)
    except ValueError as e:
        raise RuntimeError(
            f"systems-schema.sql is missing the {_ATLAS_SCHEMA_BEGIN_MARKER}/"
            f"{_ATLAS_SCHEMA_END_MARKER} block — has the schema been "
            "restructured?"
        ) from e
    return text[start:end]


def ensure_atlas_schema(conn: sqlite3.Connection) -> None:
    """Apply the canonical atlas_* DDL from systems-schema.sql.

    Idempotent: all statements inside the block use CREATE TABLE / INDEX
    IF NOT EXISTS, so running this on an already-migrated DB is a no-op.

    Also ensures the meta stamp table exists (#855, #856).
    """
    conn.executescript(_load_atlas_schema())
    conn.executescript(
        """
        CREATE TABLE IF NOT EXISTS meta (
            generator_name TEXT PRIMARY KEY,
            schema_version TEXT NOT NULL,
            generator_sha TEXT NOT NULL,
            generated_at TEXT NOT NULL DEFAULT (datetime('now'))
        );
        """
    )


def _first_int(values, default: int = 0) -> int:
    """Coerce the first numeric element of an iterable (e.g. [r, c]) to int."""
    try:
        return int(values[0])
    except (TypeError, ValueError, IndexError):
        return default


def sync_markers_to_db(
    conn: sqlite3.Connection,
    body_id: str,
    markers: dict,
) -> dict:
    """Refresh atlas_* rows for a single body from its markers.json contents.

    Deletes any existing rows for the body (cascading across all seven atlas
    tables) and re-inserts from the markers dict. Returns a count dict for
    reporting.
    """
    # Wipe existing rows for this body — fully deterministic rebuild.
    for table in (
        "atlas_cities",
        "atlas_roads",
        "atlas_railroads",
        "atlas_pois",
        "atlas_rivers",
        "atlas_oceans",
        "atlas_mountain_ranges",
    ):
        conn.execute(f"DELETE FROM {table} WHERE body_id = ?", (body_id,))

    grid = markers.get("grid") or {}
    grid_w = int(grid.get("w", GRID_W))
    grid_h = int(grid.get("h", GRID_H))
    conn.execute(
        """INSERT INTO atlas_body_grids (body_id, grid_w, grid_h, updated_at)
           VALUES (?, ?, ?, datetime('now'))
           ON CONFLICT(body_id) DO UPDATE SET
               grid_w = excluded.grid_w,
               grid_h = excluded.grid_h,
               updated_at = excluded.updated_at""",
        (body_id, grid_w, grid_h),
    )

    counts = {
        "cities": 0, "roads": 0, "railroads": 0, "pois": 0,
        "rivers": 0, "oceans": 0, "mountain_ranges": 0,
    }

    for city in markers.get("cities") or []:
        local = city.get("id")
        if not local:
            continue
        center = city.get("center") or [0, 0]
        conn.execute(
            """INSERT INTO atlas_cities
               (city_id, body_id, local_id, name, kind, center_row, center_col, population)
               VALUES (?, ?, ?, ?, ?, ?, ?, ?)""",
            (
                f"{body_id}/{local}",
                body_id,
                local,
                city.get("name") or "",
                city.get("kind") or "city",
                _first_int(center, 0),
                _first_int(center[1:] if len(center) > 1 else [0], 0),
                int(city.get("population") or 0),
            ),
        )
        counts["cities"] += 1

    for road in markers.get("roads") or []:
        local = road.get("id")
        if not local:
            continue
        path = road.get("path") or []
        conn.execute(
            """INSERT INTO atlas_roads
               (road_id, body_id, local_id, name, kind, point_count)
               VALUES (?, ?, ?, ?, ?, ?)""",
            (
                f"{body_id}/{local}",
                body_id,
                local,
                road.get("name") or "",
                road.get("kind") or "commercial",
                len(path),
            ),
        )
        counts["roads"] += 1

    for rail in markers.get("railroads") or []:
        local = rail.get("id")
        if not local:
            continue
        path = rail.get("path") or []
        conn.execute(
            """INSERT INTO atlas_railroads
               (railroad_id, body_id, local_id, name, kind, point_count)
               VALUES (?, ?, ?, ?, ?, ?)""",
            (
                f"{body_id}/{local}",
                body_id,
                local,
                rail.get("name") or "",
                rail.get("kind") or "passenger_freight",
                len(path),
            ),
        )
        counts["railroads"] += 1

    for poi in markers.get("pois") or []:
        local = poi.get("id")
        if not local:
            continue
        center = poi.get("center") or [0, 0]
        conn.execute(
            """INSERT INTO atlas_pois
               (poi_id, body_id, local_id, name, kind, center_row, center_col)
               VALUES (?, ?, ?, ?, ?, ?, ?)""",
            (
                f"{body_id}/{local}",
                body_id,
                local,
                poi.get("name") or "",
                poi.get("kind") or "transit",
                _first_int(center, 0),
                _first_int(center[1:] if len(center) > 1 else [0], 0),
            ),
        )
        counts["pois"] += 1

    for river in markers.get("rivers") or []:
        local = river.get("id")
        if not local:
            continue
        path = river.get("path") or []
        conn.execute(
            """INSERT INTO atlas_rivers
               (river_id, body_id, local_id, name, point_count)
               VALUES (?, ?, ?, ?, ?)""",
            (
                f"{body_id}/{local}",
                body_id,
                local,
                river.get("name") or "",
                len(path),
            ),
        )
        counts["rivers"] += 1

    for water in markers.get("oceans") or []:
        local = water.get("id")
        if not local:
            continue
        center = water.get("center") or [0, 0]
        conn.execute(
            """INSERT INTO atlas_oceans
               (water_id, body_id, local_id, name, kind, center_row, center_col, area_fraction)
               VALUES (?, ?, ?, ?, ?, ?, ?, ?)""",
            (
                f"{body_id}/{local}",
                body_id,
                local,
                water.get("name") or "",
                water.get("kind") or "ocean",
                _first_int(center, 0),
                _first_int(center[1:] if len(center) > 1 else [0], 0),
                float(water.get("area_fraction") or 0.0),
            ),
        )
        counts["oceans"] += 1

    for rng_feat in markers.get("mountain_ranges") or []:
        local = rng_feat.get("id")
        if not local:
            continue
        center = rng_feat.get("center") or [0, 0]
        peak = rng_feat.get("peak") or center
        conn.execute(
            """INSERT INTO atlas_mountain_ranges
               (range_id, body_id, local_id, name, center_row, center_col,
                peak_row, peak_col, area_cells)
               VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?)""",
            (
                f"{body_id}/{local}",
                body_id,
                local,
                rng_feat.get("name") or "",
                _first_int(center, 0),
                _first_int(center[1:] if len(center) > 1 else [0], 0),
                _first_int(peak, 0),
                _first_int(peak[1:] if len(peak) > 1 else [0], 0),
                int(rng_feat.get("area_cells") or 0),
            ),
        )
        counts["mountain_ranges"] += 1

    return counts


# ---------------------------------------------------------------------------
# Grid helpers
# ---------------------------------------------------------------------------

# markers.json is stored in pixel space per D-191 §8 — the canonical format is
# `center: [row, col]` and `path: [[row, col], ...]`. Lat/lon is a display-time
# derivation in the atlas UI, not a storage format.


def quadrant(row: int, col: int) -> int:
    """Return quadrant index 0-3: NW=0, NE=1, SW=2, SE=3."""
    h = 0 if row < GRID_H // 2 else 1
    w = 0 if col < GRID_W // 2 else 1
    return h * 2 + w


def body_id_salt(body_id: str) -> int:
    """Deterministic per-body salt from body_id."""
    return int(hashlib.sha256(body_id.encode()).hexdigest()[:8], 16)


# ---------------------------------------------------------------------------
# Body definition loading
# ---------------------------------------------------------------------------

def load_body_def(body_dir: Path) -> dict | None:
    """Load body definition from body index.md frontmatter."""
    index_md = body_dir / "index.md"
    if not index_md.exists():
        return None
    content = index_md.read_text()
    if not content.startswith("---"):
        return None
    try:
        end = content.index("---", 3)
        bd = yaml.safe_load(content[3:end])
    except (ValueError, yaml.YAMLError):
        return None
    if not bd or "id" not in bd or "planet_class" not in bd:
        return None
    return bd


# ---------------------------------------------------------------------------
# Terrain analysis
# ---------------------------------------------------------------------------

def _analyse_terrain(terrain: dict) -> dict:
    """Analyse terrain for atlas city placement.

    Deterministic and RNG-free: continent flood-fill, habitability scoring,
    river-mouth detection, and the terrain cost grid are all pure functions
    of the terrain dict. Per-body variation comes from the RNG used later
    in `_score_capital_sites` and `place_cities`, not from this function.

    Returns:
        continents:       list of continent dicts {cells, area, id}
        habitability:     float32 (H, W) — 0=uninhabitable, 1=ideal
        river_mouths:     list of (row, col) — land cells where rivers meet water
        cost_grid:        float32 (H, W) — A* traversal cost per cell
        continent_map:    int32 (H, W) — continent label per land cell (-1=water)
        land_mask:        bool (H, W)
    """
    from scipy.ndimage import label

    elevation = terrain["elevation"]
    temperature = terrain["temperature"]
    moisture = terrain["moisture"]
    surface_water = terrain["surface_water"]
    river_grid = terrain.get("river_grid", np.zeros((GRID_H, GRID_W), bool))
    sea_level = terrain["sea_level"]

    land = ~surface_water

    # Normalised elevation above sea level
    elev_norm = np.where(land,
                         (elevation - sea_level) / (1.0 - sea_level + 1e-9),
                         0.0)
    mountains = land & (elev_norm > MOUNTAIN_ELEV_THRESHOLD)

    # --- Continent detection (flood-fill connected land) ---
    continent_labels, n_continents = label(land)
    continents = []
    for lbl in range(1, n_continents + 1):
        mask = continent_labels == lbl
        area = int(mask.sum())
        if area < 50:
            continue  # ignore tiny rock outcroppings
        ys, xs = np.where(mask)
        continents.append({
            "id": lbl,
            "area": area,
            "cells": list(zip(ys.tolist(), xs.tolist())),
            "_center": (int(ys.mean()), int(xs.mean())),
        })
    continents.sort(key=lambda c: c["area"], reverse=True)

    continent_map = np.where(land, continent_labels, -1).astype(np.int32)

    # --- Habitability scoring ---
    # temperature: 0.3–0.7 is ideal (mid-range)
    temp_score = 1.0 - np.abs(temperature - 0.5) * 2.0
    temp_score = np.clip(temp_score, 0.0, 1.0)

    # moisture: 0.3–0.6 is ideal
    moist_score = np.where(
        moisture < 0.3, moisture / 0.3,
        np.where(moisture < 0.6, 1.0, 1.0 - (moisture - 0.6) / 0.4)
    )
    moist_score = np.clip(moist_score, 0.0, 1.0)

    # slope: low slope = ideal (flat land)
    dy = np.gradient(elevation, axis=0)
    dx = np.gradient(elevation, axis=1)
    slope = np.sqrt(dy**2 + dx**2)
    slope_max = float(slope[land].max()) if land.any() else 1.0
    slope_score = 1.0 - np.clip(slope / (slope_max + 1e-9), 0.0, 1.0)

    habitability = (temp_score * 0.4 + moist_score * 0.3 + slope_score * 0.3)
    habitability = np.where(land, habitability, 0.0).astype(np.float32)

    # Coastal access bonus: land cells adjacent to water
    water_f = surface_water.astype(np.float32)
    from scipy.ndimage import uniform_filter
    water_proximity = uniform_filter(water_f, size=7)
    coastal_bonus = np.clip(water_proximity * 2.0, 0.0, 0.3)
    habitability = np.where(land, habitability + coastal_bonus.astype(np.float32), 0.0)
    habitability = np.clip(habitability, 0.0, 1.0).astype(np.float32)

    # --- River mouth detection ---
    # A river mouth is a river-grid cell that is land but adjacent to water
    river_mouths = []
    if river_grid.any() and surface_water.any():
        # Dilate surface_water by 1 cell
        from scipy.ndimage import binary_dilation
        water_dilated = binary_dilation(surface_water)
        # River cells on land that touch water
        mouth_mask = river_grid & land & water_dilated
        rys, rxs = np.where(mouth_mask)
        # Deduplicate by proximity (cluster to strongest mouth)
        if len(rys) > 0:
            used = np.zeros(len(rys), bool)
            for i in range(len(rys)):
                if used[i]:
                    continue
                used[i] = True
                r, c = int(rys[i]), int(rxs[i])
                # Suppress nearby duplicates within radius 10
                dist = np.sqrt((rys - r)**2 + (rxs - c)**2)
                used[dist < 10] = True
                river_mouths.append((r, c))

    # --- Terrain cost grid (D-191 §3) ---
    cost_grid = np.full((GRID_H, GRID_W), COST_FLAT, dtype=np.float32)
    cost_grid[surface_water] = COST_WATER
    cost_grid[mountains] = COST_MOUNTAIN
    cost_grid[river_grid & land] = COST_RIVER

    return {
        "continents": continents,
        "habitability": habitability,
        "river_mouths": river_mouths,
        "cost_grid": cost_grid,
        "continent_map": continent_map,
        "land_mask": land,
    }


# ---------------------------------------------------------------------------
# City count computation (D-191 §8)
# ---------------------------------------------------------------------------

def distribute_population(total_population: int, n_cities: int) -> list[int]:
    """Split a body's population across n_cities using a simple geometric decay.

    Capital gets ~50%, each subsequent city gets half of what the previous one
    did, and the remaining tail is absorbed by the capital so totals still match.
    Returns a list of int populations (rounded to the nearest 10,000 for realism).
    """
    if n_cities <= 0 or total_population <= 0:
        return []
    if n_cities == 1:
        return [int(total_population)]

    weights = [0.5 ** i for i in range(n_cities)]
    weight_sum = sum(weights)
    shares = [w / weight_sum for w in weights]

    populations = []
    remaining = total_population
    for i, share in enumerate(shares):
        if i == len(shares) - 1:
            pop = remaining
        else:
            pop = int(round(total_population * share / 10_000) * 10_000)
            pop = min(pop, remaining - (len(shares) - i - 1) * 10_000)
            pop = max(pop, 10_000)
        populations.append(pop)
        remaining -= pop
    return populations


def compute_city_count(population: int, settlement_pattern: str | None) -> int:
    """Compute number of cities for a body.

    Base count: floor(log10(pop / 1_000_000)), minimum 1 if inhabited.
    Settlement pattern modifiers from D-191 §8.
    """
    if population <= 0:
        return 0

    if settlement_pattern in ("orbital_only",):
        return 0

    if settlement_pattern in ("domed", "cave"):
        return 1

    base = int(math.floor(math.log10(max(population, 1_000_001) / 1_000_000)))
    base = max(1, base)  # at least 1 city if inhabited

    mod = SETTLEMENT_PATTERN_MODIFIERS.get(settlement_pattern, 0) or 0
    result = max(1, base + mod)
    return min(result, 12)  # cap at 12 to avoid overcrowding small maps


# ---------------------------------------------------------------------------
# City placement (D-191 §3)
# ---------------------------------------------------------------------------

def _score_capital_sites(
    analysis: dict,
    rng: np.random.Generator,
    noise_factor: float = 0.25,
) -> np.ndarray:
    """Score all land cells for capital placement.

    Capital scoring factors:
    - Habitability score (temperature, moisture, slope)
    - River mouth bonus (~50% of capitals at river mouths per D-191 §3)
    - ±25% noise variation per seed
    """
    habitability = analysis["habitability"]
    land = analysis["land_mask"]
    river_mouths = analysis["river_mouths"]

    # Base score
    score = habitability.copy()

    # River mouth bonus (makes ~50% of capital choices go there).
    # Build a single sparse accumulator with all mouth points set, then
    # run gaussian_filter once — O(1) filter calls rather than O(n_mouths).
    if river_mouths:
        from scipy.ndimage import gaussian_filter
        mouth_field = np.zeros((GRID_H, GRID_W), dtype=np.float32)
        for r, c in river_mouths:
            mouth_field[r, c] = 1.0
        river_bonus = gaussian_filter(mouth_field, sigma=8.0)
        river_bonus = river_bonus / (river_bonus.max() + 1e-9)
        score = score * 0.5 + river_bonus * 0.5

    # ±25% noise for seed variation (D-191 §3)
    noise = rng.uniform(1.0 - noise_factor, 1.0 + noise_factor,
                        size=(GRID_H, GRID_W)).astype(np.float32)
    score = score * noise
    score = np.where(land, score, 0.0)

    return score


def _quadrant_counts(placed: list[tuple[int, int]]) -> np.ndarray:
    """Return count of placed cities per quadrant (4 quadrants)."""
    counts = np.zeros(4, dtype=int)
    for r, c in placed:
        counts[quadrant(r, c)] += 1
    return counts


def _dijkstra_from_sources(
    cost_grid: np.ndarray,
    sources: list[tuple[int, int]],
) -> np.ndarray:
    """Multi-source Dijkstra on cost_grid from all source cells.

    Returns distance array (H, W). Uses 8-connectivity.
    """
    dist = np.full((GRID_H, GRID_W), np.inf, dtype=np.float64)
    heap = []
    for r, c in sources:
        if cost_grid[r, c] < 1e8:
            dist[r, c] = 0.0
            heapq.heappush(heap, (0.0, r, c))

    while heap:
        d, r, c = heapq.heappop(heap)
        if d > dist[r, c]:
            continue
        for dr in (-1, 0, 1):
            for dc in (-1, 0, 1):
                if dr == 0 and dc == 0:
                    continue
                nr, nc = r + dr, c + dc
                if nr < 0 or nr >= GRID_H or nc < 0 or nc >= GRID_W:
                    continue
                edge = cost_grid[nr, nc]
                if edge >= 1e8:
                    continue
                step = math.sqrt(2) if dr != 0 and dc != 0 else 1.0
                nd = d + edge * step
                if nd < dist[nr, nc]:
                    dist[nr, nc] = nd
                    heapq.heappush(heap, (nd, nr, nc))

    return dist


def place_cities(
    analysis: dict,
    n_cities: int,
    rng: np.random.Generator,
    noise_factor: float = 0.25,
    body_population: int = 0,
) -> list[dict]:
    """Place cities sequentially per D-191 §3.

    Returns a list of city dicts in the canonical markers.json schema
    (D-191 §8 — pixel space):
        {id, name, kind, center: [row, col], population}
    with internal `_row`, `_col`, `_continent_id` fields used by downstream
    infrastructure generation. The underscore-prefixed fields are stripped
    before writing to disk.
    """
    if n_cities <= 0:
        return []

    continents = analysis["continents"]
    habitability = analysis["habitability"]
    cost_grid = analysis["cost_grid"]
    continent_map = analysis["continent_map"]
    land = analysis["land_mask"]

    if not continents:
        return []

    placed_coords: list[tuple[int, int]] = []
    placed_continents: set[int] = set()
    cities = []

    # ── 1. Capital ──────────────────────────────────────────────────────────
    capital_score = _score_capital_sites(analysis, rng, noise_factor)

    # Suppress edge pixels (poles tend to be degenerate)
    capital_score[:5, :] = 0
    capital_score[-5:, :] = 0

    best = int(np.argmax(capital_score))
    cap_r, cap_c = divmod(best, GRID_W)

    continent_id = int(continent_map[cap_r, cap_c]) if continent_map[cap_r, cap_c] > 0 else 1
    placed_coords.append((cap_r, cap_c))
    placed_continents.add(continent_id)

    cities.append({
        "id": "city_0",
        "name": "",
        "kind": "capital",
        "center": [cap_r, cap_c],
        "population": 0,  # filled after all cities are placed
        "_row": cap_r,
        "_col": cap_c,
        "_continent_id": continent_id,
    })

    if n_cities == 1:
        _assign_populations(cities, body_population)
        return cities

    # Minimum exclusion radius: varies by city count (more cities → tighter spacing)
    excl_radius = max(20, GRID_H // (n_cities + 1))

    # ── 2. Subsequent cities (corridor growth + quadrant spread) ────────────
    for i in range(1, n_cities):
        q_counts = _quadrant_counts(placed_coords)

        # Multi-source Dijkstra from all placed cities along cost grid
        corridor_dist = _dijkstra_from_sources(cost_grid, placed_coords)

        # Score for new city: habitability weighted by corridor distance
        # Prefer connected but not too-close (distance 50–300 from existing cities)
        ideal_min = excl_radius
        ideal_max = min(300, excl_radius * 6)
        dist_score = np.where(
            (corridor_dist >= ideal_min) & (corridor_dist < ideal_max),
            1.0 - (corridor_dist - ideal_min) / (ideal_max - ideal_min + 1e-9),
            np.where(corridor_dist < ideal_min, 0.0, 0.1),
        ).astype(np.float32)

        site_score = habitability * 0.6 + dist_score * 0.4

        # Quadrant penalty: if a quadrant is already at QUADRANT_SATURATION,
        # reduce its score so subsequent cities prefer emptier quadrants.
        q_penalty = np.ones((GRID_H, GRID_W), np.float32)
        for qr in (0, 1):
            for qc in (0, 1):
                q_idx = qr * 2 + qc
                if q_counts[q_idx] >= QUADRANT_SATURATION:
                    rlo, rhi = (0, GRID_H // 2) if qr == 0 else (GRID_H // 2, GRID_H)
                    clo, chi = (0, GRID_W // 2) if qc == 0 else (GRID_W // 2, GRID_W)
                    q_penalty[rlo:rhi, clo:chi] = 0.3

        site_score = site_score * q_penalty

        # New-continent bonus at cities 3-4 (D-191 §3: port on new continent)
        if i in (2, 3) and len(continents) > 1:
            unexplored_conts = [
                c for c in continents
                if c["id"] not in placed_continents and c["area"] > 200
            ]
            if unexplored_conts:
                cont = unexplored_conts[0]
                # Find coastline cells on that continent
                coast_bonus = np.zeros((GRID_H, GRID_W), np.float32)
                cont_mask = continent_map == cont["id"]
                # Coastal cells: continent land cells adjacent to water
                from scipy.ndimage import binary_dilation
                water_dilated = binary_dilation(~analysis["land_mask"])
                coast_mask = cont_mask & water_dilated
                coast_bonus[coast_mask] = 0.4
                site_score = site_score + coast_bonus

        # ±25% noise
        noise = rng.uniform(1.0 - noise_factor, 1.0 + noise_factor,
                            size=(GRID_H, GRID_W)).astype(np.float32)
        site_score = site_score * noise
        site_score = np.where(land, site_score, 0.0)

        # Suppress already-placed city zones
        for pr, pc in placed_coords:
            rlo = max(0, pr - excl_radius)
            rhi = min(GRID_H, pr + excl_radius)
            clo = max(0, pc - excl_radius)
            chi = min(GRID_W, pc + excl_radius)
            site_score[rlo:rhi, clo:chi] = 0.0

        if site_score.max() < 1e-6:
            break  # no more habitable land

        best = int(np.argmax(site_score))
        r, c = divmod(best, GRID_W)
        cont_id = int(continent_map[r, c]) if continent_map[r, c] > 0 else continent_id
        placed_coords.append((r, c))
        placed_continents.add(cont_id)

        cities.append({
            "id": f"city_{i}",
            "name": "",
            "kind": "city",
            "center": [r, c],
            "population": 0,
            "_row": r,
            "_col": c,
            "_continent_id": cont_id,
        })

    _enforce_unique_city_coords(cities, analysis["land_mask"])
    _assign_populations(cities, body_population)
    return cities


def _assign_populations(cities: list[dict], body_population: int) -> None:
    """Fill the `population` field on each city from the body's total."""
    pops = distribute_population(body_population, len(cities))
    for city, pop in zip(cities, pops):
        city["population"] = int(pop)


def _enforce_unique_city_coords(cities: list[dict], land: np.ndarray) -> None:
    """Guarantee no two cities share the same (row, col).

    Rare on real heightmaps but possible when grids are small and the
    quadrant-saturation penalty pushes candidates into tight corners. If
    two cities land on identical pixels the MST treats them as zero-
    distance nodes and A* produces an empty path, silently skipping the
    edge. We deterministically perturb duplicates by walking outward in
    a fixed spiral until a free, walkable land cell is found; the search
    order is fully determined by `city_index` so the operation stays
    byte-equivalent across runs.
    """
    if len(cities) <= 1:
        return

    # Fixed spiral offsets — small radius first, then widen.
    spiral: list[tuple[int, int]] = []
    for radius in range(1, 12):
        for dr in range(-radius, radius + 1):
            for dc in range(-radius, radius + 1):
                if abs(dr) == radius or abs(dc) == radius:
                    spiral.append((dr, dc))

    occupied: set[tuple[int, int]] = set()
    for idx, city in enumerate(cities):
        center = (city["_row"], city["_col"])
        if center not in occupied:
            occupied.add(center)
            continue
        # Duplicate — walk the spiral for the first free land cell.
        for dr, dc in spiral:
            nr, nc = center[0] + dr, center[1] + dc
            if 0 <= nr < GRID_H and 0 <= nc < GRID_W and land[nr, nc] and (nr, nc) not in occupied:
                print(
                    f"    warning: city {idx} at {center} collided with an "
                    f"earlier placement — deterministically perturbed to ({nr}, {nc})"
                )
                city["_row"] = nr
                city["_col"] = nc
                city["center"] = [nr, nc]
                occupied.add((nr, nc))
                break
        else:
            # Fallback: land is saturated — just keep the duplicate; the
            # MST edge will collapse but the body is a degenerate case.
            print(
                f"    warning: city {idx} at {center} has no free neighbor "
                f"— duplicate allowed (degenerate body)"
            )
            occupied.add(center)


# ---------------------------------------------------------------------------
# Infrastructure generation (D-191 §3)
# ---------------------------------------------------------------------------

def _astar_path(
    cost_grid: np.ndarray,
    start: tuple[int, int],
    goal: tuple[int, int],
) -> list[tuple[int, int]]:
    """A* pathfinding on cost_grid (8-connectivity). Returns path or []."""
    sr, sc = start
    gr, gc = goal

    def h(r, c):
        return math.sqrt((r - gr)**2 + (c - gc)**2)

    dist = {(sr, sc): 0.0}
    prev = {}
    heap = [(h(sr, sc), 0.0, sr, sc)]

    while heap:
        _, d, r, c = heapq.heappop(heap)
        if r == gr and c == gc:
            # Reconstruct path
            path = []
            cur = (gr, gc)
            while cur in prev:
                path.append(cur)
                cur = prev[cur]
            path.append((sr, sc))
            path.reverse()
            return path
        if d > dist.get((r, c), float("inf")) + 1e-9:
            continue
        for dr in (-1, 0, 1):
            for dc in (-1, 0, 1):
                if dr == 0 and dc == 0:
                    continue
                nr, nc = r + dr, c + dc
                if nr < 0 or nr >= GRID_H or nc < 0 or nc >= GRID_W:
                    continue
                edge = cost_grid[nr, nc]
                if edge >= 1e8:
                    continue
                step = math.sqrt(2) if dr != 0 and dc != 0 else 1.0
                nd = d + edge * step
                if nd < dist.get((nr, nc), float("inf")):
                    dist[(nr, nc)] = nd
                    prev[(nr, nc)] = (r, c)
                    heapq.heappush(heap, (nd + h(nr, nc), nd, nr, nc))

    return []


def _mst_edges(n_cities: int, city_coords: list[tuple[int, int]]) -> list[tuple[int, int]]:
    """Prim's MST on Euclidean distances between cities. Returns list of (i, j) edges."""
    if n_cities <= 1:
        return []

    in_tree = {0}
    edges = []

    while len(in_tree) < n_cities:
        best_cost = float("inf")
        best_edge = None
        for i in in_tree:
            for j in range(n_cities):
                if j in in_tree:
                    continue
                r1, c1 = city_coords[i]
                r2, c2 = city_coords[j]
                cost = math.sqrt((r1 - r2)**2 + (c1 - c2)**2)
                if cost < best_cost:
                    best_cost = cost
                    best_edge = (i, j)
        if best_edge is None:
            break
        edges.append(best_edge)
        in_tree.add(best_edge[1])

    return edges


def _subsample_path(path: list[tuple[int, int]], max_points: int = 64) -> list[list[int]]:
    """Subsample a path to at most max_points for compact JSON storage."""
    if len(path) <= max_points:
        return [[r, c] for r, c in path]
    step = len(path) / max_points
    result = []
    for i in range(max_points):
        idx = int(i * step)
        r, c = path[idx]
        result.append([r, c])
    # Always include endpoint
    r, c = path[-1]
    if result[-1] != [r, c]:
        result.append([r, c])
    return result


def generate_infrastructure(
    analysis: dict,
    cities: list[dict],
) -> tuple[list[dict], list[dict]]:
    """Generate roads and railroads connecting cities.

    Roads: A* paths on terrain cost grid connecting each MST edge.
    Railroads: same MST edges but via a cost grid favouring corridors.

    Returns (roads, railroads) — each is a list of path dicts.
    """
    if len(cities) < 2:
        return [], []

    cost_grid = analysis["cost_grid"]
    city_coords = [((c["_row"], c["_col"])) for c in cities]
    n = len(city_coords)

    mst = _mst_edges(n, city_coords)

    roads = []
    railroads = []
    road_cost = cost_grid.copy()
    # Rail cost: slightly prefer following roads (lower cost after first pass)
    rail_cost = cost_grid.copy()

    for edge_idx, (i, j) in enumerate(mst):
        start = city_coords[i]
        goal = city_coords[j]

        road_path = _astar_path(road_cost, start, goal)
        if road_path:
            roads.append({
                "id": f"road_{edge_idx}",
                "name": "",
                "kind": "commercial",
                "path": _subsample_path(road_path),
            })
            # After placing a road, reduce cost along it for rail (rail follows roads)
            for r, c in road_path:
                if rail_cost[r, c] < 1e8:
                    rail_cost[r, c] = max(0.3, rail_cost[r, c] * 0.5)

        rail_path = _astar_path(rail_cost, start, goal)
        if rail_path:
            railroads.append({
                "id": f"railroad_{edge_idx}",
                "name": "",
                "kind": "passenger_freight",
                "path": _subsample_path(rail_path),
            })

    return roads, railroads


# ---------------------------------------------------------------------------
# Gate terminal POI placement (D-191 §8)
# ---------------------------------------------------------------------------

def place_gate_terminal(
    cities: list[dict],
    rng: np.random.Generator,
) -> list[dict]:
    """Place a gate terminal POI at the largest population centre (sometimes scatter).

    Primary placement: city 0 (capital, largest population center).
    Occasional scatter: 15% chance to a secondary city (per D-191 §8).

    Returns a list with a single POI dict in the canonical template schema:
    `{id, name, kind: "transit", center: [row, col]}`. The name is left empty
    for gemma_naming.py (#833) to fill.
    """
    if not cities:
        return []

    scatter_prob = 0.15
    gate_city = cities[0]
    if len(cities) > 1 and rng.random() < scatter_prob:
        gate_city = cities[int(rng.integers(1, len(cities)))]

    return [{
        "id": "poi_0",
        "name": "",
        "kind": "transit",
        "center": [gate_city["_row"], gate_city["_col"]],
    }]


# ---------------------------------------------------------------------------
# Markers.json update
# ---------------------------------------------------------------------------

class AtlasGridMismatch(Exception):
    """Raised when a markers.json grid header does not match generator constants.

    If a hand-authored template ships with, say, `{"w": 1024, "h": 512}`
    and the generator overlays new cities computed against the 512 × 256
    cost grid, every coordinate is half-scale and every marker is broken.
    Fail loud here rather than silently produce corrupt output.
    """


def load_markers(body_dir: Path) -> dict:
    """Load existing markers.json, return empty structure if missing.

    Validates the grid header against the generator constants so the
    caller can trust that overlaid city/road/poi coordinates are in
    the same pixel space as the loaded geographic features. A mismatch
    raises `AtlasGridMismatch` — regenerating a markers.json against
    a different grid would corrupt every coordinate in it.
    """
    markers_path = body_dir / "markers.json"
    if markers_path.exists():
        with open(markers_path) as f:
            markers = json.load(f)
        grid = markers.get("grid") or {}
        g_w = grid.get("w")
        g_h = grid.get("h")
        if g_w != GRID_W or g_h != GRID_H:
            raise AtlasGridMismatch(
                f"{markers_path} has grid {{w: {g_w}, h: {g_h}}} but the "
                f"atlas generator runs against {{w: {GRID_W}, h: {GRID_H}}}. "
                "Either regenerate the heightmap pipeline at the generator "
                "resolution, or update GRID_W / GRID_H to match the source."
            )
        return markers
    return {
        "grid": {"w": GRID_W, "h": GRID_H},
        "rivers": [],
        "oceans": [],
        "mountain_ranges": [],
        "roads": [],
        "cities": [],
        "railroads": [],
        "pois": [],
    }


# ---------------------------------------------------------------------------
# DB queries
# ---------------------------------------------------------------------------

def query_inhabited_bodies(conn: sqlite3.Connection) -> list[dict]:
    """Query systems.db for all inhabited bodies with terrain_reference set."""
    rows = conn.execute("""
        SELECT
            b.body_id,
            b.system_id,
            b.terrain_reference,
            b.population,
            b.settlement_pattern,
            b.planet_class,
            b.economic_role,
            COALESCE(b.cultural_corridor, s.cultural_corridor) AS cultural_corridor
        FROM bodies b
        JOIN star_systems s ON b.system_id = s.system_id
        WHERE b.inhabited = 1
          AND b.terrain_reference IS NOT NULL
          AND b.population > 0
        ORDER BY b.system_id, b.body_id
    """).fetchall()

    return [
        {
            "body_id": row[0],
            "system_id": row[1],
            "terrain_reference": row[2],
            "population": row[3] or 0,
            "settlement_pattern": row[4],
            "planet_class": row[5],
            "economic_role": row[6],
            "cultural_corridor": row[7],
        }
        for row in rows
    ]


# ---------------------------------------------------------------------------
# Main pipeline
# ---------------------------------------------------------------------------

def process_body(
    body_info: dict,
    seed: int,
    noise_factor: float,
    dry_run: bool,
    force: bool,
    verbose: bool,
) -> dict:
    """Process one body.

    Returns a dict with:
        status:  'generated' | 'already_populated' | 'gas_giant' | 'error'
        message: human-readable detail (on error)
        markers: the markers.json dict if one was produced or loaded (for DB sync)
    """
    body_id = body_info["body_id"]
    terrain_ref = body_info["terrain_reference"]
    population = body_info["population"]
    settlement_pattern = body_info["settlement_pattern"]

    # Locate body directory from terrain_reference (repo-root relative).
    body_dir = REPO_ROOT / Path(terrain_ref).parent

    if not body_dir.exists():
        return {"status": "error", "message": f"body_dir not found: {body_dir}"}

    # Incremental check: if the markers file already has cities, we keep it as
    # the source of truth and return it for DB sync instead of regenerating.
    # `load_markers` validates the grid header; a mismatch on a hand-authored
    # template is reported as an error so downstream DB sync doesn't silently
    # index a file that disagrees with the generator's pixel space.
    if not force:
        markers_path = body_dir / "markers.json"
        if markers_path.exists():
            try:
                existing = load_markers(body_dir)
            except AtlasGridMismatch as e:
                return {"status": "error", "message": str(e)}
            except json.JSONDecodeError:
                existing = None  # corrupted markers.json → regenerate below
            if existing and existing.get("cities"):
                return {
                    "status": "already_populated",
                    "markers": existing,
                }

    bd = load_body_def(body_dir)
    if not bd:
        return {"status": "error", "message": f"no body definition for {body_id}"}

    # Seeded RNG: deterministic per (world_seed, body_id)
    body_salt = body_id_salt(body_id)
    rng = np.random.default_rng(seed ^ body_salt)

    try:
        terrain = simulate(bd)
    except Exception as e:
        return {"status": "error", "message": f"simulate failed: {e}"}

    if not terrain:
        return {"status": "gas_giant"}

    analysis = _analyse_terrain(terrain)

    n_cities = compute_city_count(population, settlement_pattern)

    if verbose:
        print(f"    {body_id}: pop={population:,} pattern={settlement_pattern} "
              f"→ {n_cities} cities, {len(analysis['river_mouths'])} river mouths, "
              f"{len(analysis['continents'])} continents")

    cities = place_cities(
        analysis, n_cities, rng, noise_factor, body_population=population
    )
    roads, railroads = generate_infrastructure(analysis, cities)
    pois = place_gate_terminal(cities, rng)

    # Load existing markers (preserves rivers/oceans/mountain_ranges) and
    # overlay the new generator output. `load_markers` asserts the loaded
    # grid header matches the generator constants — a mismatch would
    # silently corrupt every coordinate in the file.
    try:
        markers = load_markers(body_dir)
    except AtlasGridMismatch as e:
        return {"status": "error", "message": str(e)}

    output_cities = [
        {k: v for k, v in city.items() if not k.startswith("_")}
        for city in cities
    ]
    markers["cities"] = output_cities
    markers["roads"] = roads
    markers["railroads"] = railroads
    markers["pois"] = pois

    if not dry_run:
        with open(body_dir / "markers.json", "w") as f:
            json.dump(markers, f, indent=2)

    return {"status": "generated", "markers": markers}


def main():
    parser = argparse.ArgumentParser(
        description="Atlas generation — terrain-aware city placement and infrastructure (#832)"
    )
    parser.add_argument("--db", default=str(DB_PATH), help="Path to systems.db")
    parser.add_argument("--body", help="Process only this body_id")
    parser.add_argument("--force", action="store_true",
                        help="Regenerate even if cities already populated")
    parser.add_argument("--dry-run", action="store_true",
                        help="Run analysis but do not write markers.json")
    parser.add_argument("--seed", type=int, default=42,
                        help="World seed for deterministic placement (default: 42)")
    parser.add_argument("--noise", type=float, default=0.25,
                        help="City placement noise factor ±N (default: 0.25 = ±25%%)")
    parser.add_argument("--verbose", action="store_true",
                        help="Print per-body detail")
    args = parser.parse_args()

    db_path = Path(args.db)
    if not db_path.exists():
        print(f"error: {db_path} not found", file=sys.stderr)
        sys.exit(1)

    print(f"\n  Atlas Generation Pipeline (#832)")
    print(f"  DB:   {db_path}")
    print(f"  Seed: {args.seed}")
    if args.dry_run:
        print(f"  Mode: DRY RUN (no markers.json written, no DB sync)")
    if args.force:
        print(f"  Force: enabled (will overwrite existing city placements)")
    print()

    conn = sqlite3.connect(str(db_path))
    conn.execute("PRAGMA foreign_keys=ON")
    ensure_atlas_schema(conn)

    bodies = query_inhabited_bodies(conn)
    if args.body:
        bodies = [b for b in bodies if b["body_id"] == args.body]
        if not bodies:
            print(f"error: body '{args.body}' not found or has no terrain_reference",
                  file=sys.stderr)
            conn.close()
            sys.exit(1)

    print(f"  {len(bodies)} inhabited bodies with terrain_reference\n")

    t_total = time.time()
    n_generated = 0
    n_already = 0
    n_gas = 0
    n_errors = 0
    total_counts = {
        "cities": 0, "roads": 0, "railroads": 0, "pois": 0,
        "rivers": 0, "oceans": 0, "mountain_ranges": 0,
    }

    for i, body_info in enumerate(bodies):
        body_id = body_info["body_id"]
        t0 = time.time()

        result = process_body(
            body_info,
            seed=args.seed,
            noise_factor=args.noise,
            dry_run=args.dry_run,
            force=args.force,
            verbose=args.verbose,
        )

        elapsed = time.time() - t0
        status = result["status"]

        if status == "generated":
            n_generated += 1
            if not args.dry_run:
                counts = sync_markers_to_db(conn, body_id, result["markers"])
                for k, v in counts.items():
                    total_counts[k] += v
            print(f"  [{i+1}/{len(bodies)}] {body_id:20s} generated ({elapsed:.1f}s)")
        elif status == "already_populated":
            n_already += 1
            if not args.dry_run:
                counts = sync_markers_to_db(conn, body_id, result["markers"])
                for k, v in counts.items():
                    total_counts[k] += v
            if args.verbose:
                print(f"  [{i+1}/{len(bodies)}] {body_id:20s} already populated — DB sync only")
        elif status == "gas_giant":
            n_gas += 1
            if args.verbose:
                print(f"  [{i+1}/{len(bodies)}] {body_id:20s} gas giant — no surface")
        elif status == "error":
            n_errors += 1
            print(f"  [{i+1}/{len(bodies)}] {body_id:20s} ERROR: {result.get('message', '')}")

    if not args.dry_run:
        # Stamp generator metadata (#855, #856) together with the atlas data
        # in a single commit — atlas data + stamp land atomically, and the
        # stamp function itself no longer commits (H1).  Failure of the stamp
        # write rolls back the atlas data too rather than leaving a stamped-
        # but-missing-data intermediate state.
        try:
            _write_stamp(conn)
            conn.commit()
            print("  Stamped: generate_atlas")
        except Exception as exc:  # noqa: BLE001
            conn.rollback()
            print(f"  WARNING: failed to write generator stamp: {exc}", file=sys.stderr)
            print("  Atlas data NOT committed — regen required.", file=sys.stderr)
    conn.close()

    elapsed_total = time.time() - t_total

    print(f"\n  Done: {elapsed_total:.0f}s")
    print(f"    generated:         {n_generated}")
    print(f"    already populated: {n_already}")
    print(f"    gas giants:        {n_gas}")
    print(f"    errors:            {n_errors}")

    if args.dry_run:
        print(f"\n  Dry run — no markers.json files were written and no DB rows were touched.")
    else:
        print(f"\n  Atlas index refreshed:")
        for k in ("cities", "roads", "railroads", "pois", "rivers", "oceans", "mountain_ranges"):
            print(f"    atlas_{k:16s} {total_counts[k]:6d} rows")
        if n_generated > 0:
            print(f"\n  {n_generated} markers.json files updated with cities, roads, railroads, pois.")
            print(f"  Run gemma_naming.py (#833) to fill city and feature names.")

    print()

    if n_errors > 0:
        sys.exit(1)


if __name__ == "__main__":
    main()