import { INF, MAP_H, MAP_W, SIZE, MinHeap, clamp, fbm, hash2, indexOf, inside, nearMapEdge, pickEntities, rand, valueNoise, xyOf } from "./mapUtils.js";
import {
  aStar,
  averagePathField,
  compactPathArray,
  distanceToNearest,
  getDegree,
  incrementDegree,
  influenceFromPaths,
  influenceFromPoints,
  makeTransportCost,
  nearestConnectable,
  neighbors8,
  pathCompactness,
  pathEndpointDistance,
  pathLength,
  pathOverlapRatio,
  samplePath,
  smoothPathByLineOfSight,
} from "./mapGeneratorHelpers.js";
import { LANDUSE, isBuiltLanduse, isUrbanResidentialLanduse } from "./landuseCodes.js";

export function generateMapFeatures(seed, terrain) {
  const {
    elevation,
    moisture,
    slope,
    sea,
    river,
    floodplain,
    plain,
    agriculture,
    ridgeField,
    valleyField,
    basinField,
    coastalLowland,
    flowAccum,
    arcSpineField,
    branchRidgeField,
    depositionalLowland,
    alluvialFanField,
    deltaField,
    portSuitability,
    crossingSuitability,
    passSuitability,
    prefectureMask,
    prefectureRegionId,
  } = terrain;

  function pickPoints(scoreArray, { threshold, max, minDistance, seedOffset = 0, predicate = () => true }) {
    const candidates = [];
    for (let y = 2; y < MAP_H - 2; y++) {
      for (let x = 2; x < MAP_W - 2; x++) {
        const i = indexOf(x, y);
        if (!predicate(x, y, i)) continue;
        const score = scoreArray[i] + hash2(x, y, seed + seedOffset) * 0.08;
        if (score >= threshold) candidates.push({ x, y, score });
      }
    }
    return pickEntities(candidates, { max, minDistance, threshold, seed: seed + seedOffset });
  }

  let ports = pickPoints(portSuitability, {
    threshold: 0.3 + rand(seed, 1001) * 0.08,
    max: 3 + Math.floor(rand(seed, 1002) * 7),
    minDistance: 10,
    seedOffset: 1000,
    predicate: (x, y, i) => !sea[i],
  }).map((p) => {
    const i = indexOf(p.x, p.y);
    let seaEdge = 0;
    for (let dy = -3; dy <= 3; dy++) for (let dx = -3; dx <= 3; dx++) {
      const nx = p.x + dx;
      const ny = p.y + dy;
      if (inside(nx, ny) && sea[indexOf(nx, ny)]) seaEdge += 1 / (1 + Math.hypot(dx, dy));
    }
    const harborPotential = p.score + coastalLowland[i] * 0.28 + river[i] * 0.08 + seaEdge * 0.025 - slope[i] * 0.2;
    return { ...p, harborPotential, seaEdge, portClass: "fishing", kind: "Fishing Port" };
  }).sort((a, b) => b.harborPotential - a.harborPotential)
    .map((p, n) => {
      const isLakeLike = p.seaEdge < 0.25 && river[indexOf(p.x, p.y)] > 0.32;
      const portClass = isLakeLike ? "lake" : n === 0 ? "major" : n < 3 && p.harborPotential > 0.34 ? "regional" : "fishing";
      const kind = portClass === "major" ? "Major Port" : portClass === "regional" ? "Regional Port" : portClass === "lake" ? "Lake Port" : "Fishing Port";
      return { ...p, portClass, kind, score: p.harborPotential };
    });
  if (!ports.some((p) => p.portClass === "major")) {
    const fallbackMajor = ports.find((p) => p.portClass !== "lake") || ports[0];
    if (fallbackMajor) {
      fallbackMajor.portClass = "major";
      fallbackMajor.kind = "Major Port";
      fallbackMajor.score += 0.16;
    }
  }
  const majorPorts = ports.filter((p) => p.portClass === "major");
  const commercialPorts = ports.filter((p) => p.portClass === "major" || p.portClass === "regional");

  let crossings = pickPoints(crossingSuitability, {
    threshold: 0.28 + rand(seed, 1011) * 0.08,
    max: 8 + Math.floor(rand(seed, 1012) * 15),
    minDistance: 8,
    seedOffset: 1010,
    predicate: (x, y, i) => !sea[i],
  }).map((p) => ({ ...p, kind: "River Crossing" }));

  let passes = pickPoints(passSuitability, {
    threshold: 0.16 + rand(seed, 1021) * 0.08,
    max: 4 + Math.floor(rand(seed, 1022) * 10),
    minDistance: 9,
    seedOffset: 1020,
    predicate: (x, y, i) => !sea[i],
  }).map((p) => ({ ...p, kind: "Pass" }));

  const settlementCluster = new Float32Array(SIZE);
  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const depositional = (depositionalLowland?.[i] || 0) + (alluvialFanField?.[i] || 0) * 0.65 + (deltaField?.[i] || 0) * 0.85;
      const spineBarrier = (arcSpineField?.[i] || 0) * 0.62 + (branchRidgeField?.[i] || 0) * 0.42;
      const valleyCorridor = clamp(valleyField[i] * 0.62 + river[i] * 0.16);
      const lowlandCorridor = clamp(coastalLowland[i] * 0.38 + basinField[i] * 0.34 + plain[i] * 0.24 + agriculture[i] * 0.18 + depositional * 0.22);
      const terrainGate = clamp(1.0 - slope[i] * 1.18 - ridgeField[i] * 0.52 - spineBarrier * 0.34 - Math.max(0, elevation[i] - 0.62) * 1.35, 0.08, 1);
      const localPatch = valueNoise(x * 0.7, y * 0.7, seed + 1037, 10);
      const broadPatch = fbm(x * 0.32 + 71, y * 0.32 - 19, seed + 1038);
      settlementCluster[i] = clamp((valleyCorridor + lowlandCorridor) * terrainGate * (0.72 + broadPatch * 0.42 + localPatch * 0.18));
    }
  }

  const settlementScore = new Float32Array(SIZE);
  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      let nearFeature = 0;
      for (const p of [...ports, ...crossings, ...passes]) nearFeature = Math.max(nearFeature, 1 / (1 + Math.hypot(x - p.x, y - p.y) / 4));
      const depositional = (depositionalLowland?.[i] || 0) + (alluvialFanField?.[i] || 0) * 0.62 + (deltaField?.[i] || 0) * 0.95;
      const spineBarrier = (arcSpineField?.[i] || 0) * 0.48 + (branchRidgeField?.[i] || 0) * 0.34;
      const riverPull = Math.min(0.36, river[i] * 0.14 + valleyField[i] * 0.16 + depositional * 0.08);
      const mountainVillage = valleyField[i] * clamp(elevation[i] - 0.42, 0, 0.3) * 0.52;
      const remoteMountainPenalty = Math.max(0, elevation[i] - 0.58) * Math.max(0, ridgeField[i] + spineBarrier - 0.22) * (1 - valleyField[i]) * 0.75;
      const base = agriculture[i] * 0.50 + plain[i] * 0.14 + nearFeature * 0.23 + riverPull + basinField[i] * 0.13 + coastalLowland[i] * 0.08 + depositional * 0.13 + mountainVillage - slope[i] * 0.48 - ridgeField[i] * 0.24 - spineBarrier * 0.12 - floodplain[i] * 0.06 - remoteMountainPenalty;
      settlementScore[i] = clamp(base * (0.74 + settlementCluster[i] * 0.66) + settlementCluster[i] * 0.13);
    }
  }

  // Capacity-first settlement context.  These rasters are intentionally small
  // and reusable: all village, town, city-capacity, density and land-use passes
  // should read from the same human-geography interpretation of the terrain.
  const developable = new Float32Array(SIZE);
  const ruralSuitability = new Float32Array(SIZE);
  const townSuitability = new Float32Array(SIZE);
  const valleySettlement = new Float32Array(SIZE);
  const coastalSettlement = new Float32Array(SIZE);
  const confluenceField = new Float32Array(SIZE);

  function localConfluenceScore(x, y) {
    let arms = 0;
    let strong = 0;
    for (const [dx, dy] of [[1,0],[-1,0],[0,1],[0,-1],[1,1],[-1,1],[1,-1],[-1,-1]]) {
      const nx = x + dx;
      const ny = y + dy;
      if (!inside(nx, ny)) continue;
      const rv = river[indexOf(nx, ny)];
      if (rv > 0.22) arms++;
      if (rv > 0.38) strong++;
    }
    return clamp((arms >= 3 ? 0.16 : arms === 2 ? 0.08 : 0) + strong * 0.035);
  }

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const depositional = (depositionalLowland?.[i] || 0) + (alluvialFanField?.[i] || 0) * 0.62 + (deltaField?.[i] || 0) * 0.85;
      const highPenalty = Math.max(0, elevation[i] - 0.56);
      const lowSlope = clamp(1 - slope[i] * 2.25);
      const confluence = localConfluenceScore(x, y);
      confluenceField[i] = confluence;
      developable[i] = clamp(
        plain[i] * 0.34 +
        agriculture[i] * 0.22 +
        basinField[i] * 0.22 +
        valleyField[i] * 0.24 +
        coastalLowland[i] * 0.16 +
        depositional * 0.20 +
        lowSlope * 0.10 -
        slope[i] * 0.84 -
        ridgeField[i] * 0.54 -
        highPenalty * 1.12 -
        floodplain[i] * 0.04
      );
      valleySettlement[i] = clamp(
        valleyField[i] * 0.48 +
        river[i] * 0.16 +
        confluence * 0.72 +
        depositional * 0.14 +
        basinField[i] * 0.08 +
        lowSlope * 0.12 -
        slope[i] * 0.58 -
        ridgeField[i] * 0.32 -
        highPenalty * 0.74
      );
      coastalSettlement[i] = clamp(
        coastalLowland[i] * 0.48 +
        (portSuitability?.[i] || 0) * 0.28 +
        (deltaField?.[i] || 0) * 0.18 +
        plain[i] * 0.10 -
        slope[i] * 0.52 -
        ridgeField[i] * 0.22
      );
      ruralSuitability[i] = clamp(
        settlementScore[i] * 0.48 +
        agriculture[i] * 0.40 +
        developable[i] * 0.22 +
        valleySettlement[i] * 0.24 +
        coastalSettlement[i] * 0.16 -
        Math.max(0, elevation[i] - 0.64) * 0.54
      );
      townSuitability[i] = clamp(
        settlementScore[i] * 0.30 +
        developable[i] * 0.38 +
        valleySettlement[i] * 0.24 +
        coastalSettlement[i] * 0.20 +
        confluence * 0.34 +
        basinField[i] * 0.14 +
        plain[i] * 0.10 -
        slope[i] * 0.34 -
        ridgeField[i] * 0.16
      );
    }
  }

  function buildSettlementRegionStats() {
    const stats = new Map();
    for (let y = 1; y < MAP_H - 1; y++) {
      for (let x = 1; x < MAP_W - 1; x++) {
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const regionId = regionIdAt(x, y);
        if (regionId < 0) continue;
        let st = stats.get(regionId);
        if (!st) {
          st = { area: 0, developable: 0, valley: 0, coast: 0, basin: 0, plain: 0, town: 0 };
          stats.set(regionId, st);
        }
        st.area++;
        if (developable[i] > 0.18) st.developable++;
        if (valleySettlement[i] > 0.22) st.valley++;
        if (coastalSettlement[i] > 0.24) st.coast++;
        if (basinField[i] > 0.22) st.basin++;
        if (plain[i] > 0.24) st.plain++;
        if (townSuitability[i] > 0.28) st.town++;
      }
    }
    return stats;
  }

  const settlementRegionStats = buildSettlementRegionStats();

  function pickRegionalSettlementPoints(scoreArray, {
    totalMax,
    minDistance,
    seedOffset,
    threshold,
    quotaForRegion,
    kind,
    extraScore = () => 0,
    predicate = () => true,
  }) {
    const byRegion = new Map();
    for (let y = 2; y < MAP_H - 2; y++) {
      for (let x = 2; x < MAP_W - 2; x++) {
        const i = indexOf(x, y);
        if (sea[i] || !predicate(x, y, i)) continue;
        const regionId = regionIdAt(x, y);
        if (regionId < 0) continue;
        const score = scoreArray[i] + extraScore(x, y, i) + hash2(x, y, seed + seedOffset) * 0.055;
        if (score < threshold) continue;
        if (!byRegion.has(regionId)) byRegion.set(regionId, []);
        byRegion.get(regionId).push({ x, y, score, kind, regionId });
      }
    }
    const picked = [];
    for (const [regionId, candidates] of [...byRegion.entries()].sort((a, b) => a[0] - b[0])) {
      const quota = quotaForRegion(regionId, settlementRegionStats.get(regionId) || { area: 0 });
      if (quota <= 0) continue;
      picked.push(...pickEntities(candidates, {
        max: quota,
        minDistance,
        threshold,
        seed: seed + seedOffset + regionId * 997,
      }));
    }
    return picked
      .sort((a, b) => b.score - a.score)
      .slice(0, totalMax)
      .sort((a, b) => a.regionId - b.regionId || b.score - a.score);
  }

  const villageScore = new Float32Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    villageScore[i] = clamp(ruralSuitability[i] * 0.72 + valleySettlement[i] * 0.26 + coastalSettlement[i] * 0.18 + settlementCluster[i] * 0.12);
  }

  let villages = pickRegionalSettlementPoints(villageScore, {
    threshold: 0.26 + rand(seed, 1031) * 0.055,
    totalMax: 105 + Math.floor(rand(seed, 1032) * 35),
    minDistance: 4,
    seedOffset: 1030,
    kind: "Village",
    quotaForRegion: (regionId, st) => {
      if (!st || st.developable < 12) return 0;
      const base = 1.8 + st.developable / 58 + st.valley / 32 + st.coast / 42 + st.basin / 70;
      const selectedBonus = regionId === 0 ? 7.5 : 0;
      return Math.round(clamp(base + selectedBonus + rand(seed, 1035 + regionId * 13) * 2.4, regionId === 0 ? 10 : 2, regionId === 0 ? 26 : 12));
    },
    extraScore: (x, y, i) => confluenceField[i] * 0.10,
  }).map((p, n) => {
    const i = indexOf(p.x, p.y);
    const settlementType = valleySettlement[i] > 0.42 ? "Valley Village" : coastalSettlement[i] > 0.44 ? "Coastal Village" : "Village";
    const population = Math.round((350 + Math.pow(rand(seed, 18000 + n * 17 + p.x * 3 + p.y), 1.9) * 5200 + ruralSuitability[i] * 3200) / 100) * 100;
    return { ...p, kind: settlementType, population };
  });

  const marketScore = new Float32Array(SIZE);
  for (let y = 4; y < MAP_H - 4; y++) {
    for (let x = 4; x < MAP_W - 4; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;

      let villagePull = 0;
      let nearbyVillages = 0;
      for (const v of villages) {
        const d = Math.hypot(x - v.x, y - v.y);
        if (d < 24) {
          villagePull += 1 / (1 + d / 2.8);
          nearbyVillages++;
        }
      }

      let featurePull = 0;
      for (const p of [...ports, ...crossings]) featurePull = Math.max(featurePull, 1 / (1 + Math.hypot(x - p.x, y - p.y) / 3.2));
      const valleyMouth = valleyField[i] > 0.22 && (plain[i] > 0.22 || basinField[i] > 0.18 || coastalLowland[i] > 0.18) ? 0.12 : 0;
      marketScore[i] = clamp(
        townSuitability[i] * 0.66 +
        villagePull * 0.56 +
        featurePull * 0.26 +
        confluenceField[i] * 0.30 +
        valleyMouth +
        plain[i] * 0.12 +
        basinField[i] * 0.12 +
        (depositionalLowland?.[i] || 0) * 0.08 +
        (deltaField?.[i] || 0) * 0.07 +
        nearbyVillages * 0.006 -
        slope[i] * 0.22 -
        ridgeField[i] * 0.10 -
        (arcSpineField?.[i] || 0) * 0.06
      );
    }
  }

  let markets = pickRegionalSettlementPoints(marketScore, {
    threshold: 0.30 + rand(seed, 1041) * 0.055,
    totalMax: 34 + Math.floor(rand(seed, 1042) * 16),
    minDistance: 8,
    seedOffset: 1040,
    kind: "Market Town",
    quotaForRegion: (regionId, st) => {
      if (!st || st.town < 8) return 0;
      const base = 0.8 + st.developable / 230 + st.valley / 120 + st.coast / 140 + st.basin / 160;
      const selectedBonus = regionId === 0 ? 3.0 : 0;
      return Math.round(clamp(base + selectedBonus + rand(seed, 1045 + regionId * 17) * 1.3, regionId === 0 ? 4 : 1, regionId === 0 ? 10 : 4));
    },
    extraScore: (x, y, i) => (distanceToNearest(ports, x, y) < 7 ? 0.07 : 0) + (distanceToNearest(crossings, x, y) < 5 ? 0.06 : 0),
  }).map((p, n) => {
    const i = indexOf(p.x, p.y);
    const kind = coastalSettlement[i] > 0.46 && distanceToNearest(ports, p.x, p.y) < 9 ? "Port Town" : valleySettlement[i] > 0.42 ? "Valley Market Town" : "Market Town";
    const population = Math.round((4200 + Math.pow(rand(seed, 18100 + n * 19 + p.x * 5 + p.y), 1.55) * 26000 + marketScore[i] * 16000) / 1000) * 1000;
    return { ...p, kind, population };
  });

  const defenseScore = new Float32Array(SIZE);
  for (let y = 3; y < MAP_H - 3; y++) {
    for (let x = 3; x < MAP_W - 3; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const hillShoulder = clamp(1 - Math.abs(elevation[i] - 0.50) / 0.24);
      let riverArms = 0;
      for (const [nx, ny] of neighbors8(x, y)) if (river[indexOf(nx, ny)] > 0.32) riverArms++;
      const confluence = riverArms >= 3 ? 0.38 : riverArms === 2 ? 0.18 : 0;
      const roadJunctionProxy = (
        (distanceToNearest(markets, x, y) < 7 ? 1 : 0) +
        (distanceToNearest(crossings, x, y) < 6 ? 1 : 0) +
        (distanceToNearest(passes, x, y) < 7 ? 1 : 0) +
        (distanceToNearest(commercialPorts, x, y) < 8 ? 1 : 0)
      ) >= 2 ? 0.32 : 0;
      const hillEdge = plain[i] > 0.2 && elevation[i] > 0.36 && elevation[i] < 0.62 && (slope[i] > 0.12 || ridgeField[i] > 0.12) ? 0.3 : 0;
      const mountainRidgeCastle = elevation[i] > 0.56 && ridgeField[i] > 0.3 && valleyField[i] > 0.1 ? 0.28 : 0;
      const validCastleSite = confluence > 0 || roadJunctionProxy > 0 || hillEdge > 0 || mountainRidgeCastle > 0;
      defenseScore[i] = validCastleSite
        ? clamp(hillShoulder * 0.28 + confluence + roadJunctionProxy + hillEdge + mountainRidgeCastle + slope[i] * 0.05 - floodplain[i] * 0.42 - coastalLowland[i] * 0.12)
        : 0;
    }
  }

  let castles = pickPoints(defenseScore, {
    threshold: 0.34 + rand(seed, 1051) * 0.08,
    max: 2 + Math.floor(rand(seed, 1052) * 4),
    minDistance: 15,
    seedOffset: 1050,
    predicate: (x, y, i) => !sea[i] && defenseScore[i] > 0,
  }).map((p) => ({
    ...p,
    kind: elevation[indexOf(p.x, p.y)] > 0.55 ? "Mountain Castle" : elevation[indexOf(p.x, p.y)] > 0.38 ? "Hilltop Castle" : "Flatland Castle",
  }));

  function normalEdgePenalty(x, y) {
    if (nearMapEdge(x, y, 1)) return INF;
    if (nearMapEdge(x, y, 2)) return 7;
    if (nearMapEdge(x, y, 4)) return 2.8;
    return 0;
  }

  function premodernCost(x, y) {
    const i = indexOf(x, y);
    if (sea[i]) return INF;
    const crossingBonus = distanceToNearest(crossings, x, y) < 4 ? 0.65 : 0;
    const passBonus = distanceToNearest(passes, x, y) < 4 ? 0.45 : 0;
    const riverPenalty = river[i] > 0.28 ? (crossingBonus ? 0.45 : 2.4) : 0;
    const highMountain = elevation[i] > 0.72 ? 4.2 : elevation[i] > 0.58 ? 1.4 : 0;
    return Math.max(0.35, 1 + slope[i] * 5.8 + riverPenalty + highMountain + floodplain[i] * 0.62 - plain[i] * 0.32 - valleyField[i] * 0.42 - coastalLowland[i] * 0.12 - passBonus + normalEdgePenalty(x, y) + hash2(x, y, seed + 111) * 0.16);
  }

  const premodernRoads = [];
  function addPremodernRoad(a, b) {
    const path = aStar(a, b, premodernCost);
    if (path.length > 3) premodernRoads.push(path);
  }

  for (const castle of castles) {
    const near = pickEntities([...markets, ...ports, ...crossings, ...passes].map((p) => ({ ...p, score: 1 / (1 + Math.hypot(p.x - castle.x, p.y - castle.y)) })), { max: 2 + Math.floor(rand(seed, castle.x + castle.y) * 3), minDistance: 1, threshold: 0 });
    for (const p of near) addPremodernRoad(castle, p);
  }
  for (const market of markets) {
    const near = pickEntities([...markets.filter((p) => p !== market), ...ports, ...crossings].map((p) => ({ ...p, score: 1 / (1 + Math.hypot(p.x - market.x, p.y - market.y)) })), { max: 1 + Math.floor(rand(seed, market.x + market.y + 20) * 3), minDistance: 1, threshold: 0 });
    for (const p of near) addPremodernRoad(market, p);
  }

  function urbanSiteSuitability(p) {
    const i = indexOf(p.x, p.y);
    if (sea[i]) return 0;
    const portBonus = p.kind === "Port Town" || p.portClass === "major" || p.portClass === "regional" ? 0.18 : 0;
    const historicalBonus = p.kind === "Market City" || p.kind === "Castle Town" || p.kind === "Market Town" || p.kind === "Valley Market Town" ? 0.07 : 0;
    return clamp(
      (developable?.[i] || 0) * 0.56 +
      (townSuitability?.[i] || 0) * 0.24 +
      plain[i] * 0.22 +
      agriculture[i] * 0.10 +
      basinField[i] * 0.16 +
      coastalLowland[i] * 0.16 +
      valleyField[i] * 0.14 +
      (confluenceField?.[i] || 0) * 0.16 +
      portBonus + historicalBonus -
      slope[i] * 0.54 -
      ridgeField[i] * 0.32 -
      Math.max(0, elevation[i] - 0.55) * 1.25
    );
  }

  function estimateUrbanCapacity(p, radius = 22, densityBias = 1.0) {
    if (!p || !inside(p.x, p.y)) return 0;
    const centerRegion = regionIdAt(p.x, p.y);
    let capacity = 0;
    const r = Math.ceil(radius);
    for (let dy = -r; dy <= r; dy++) {
      for (let dx = -r; dx <= r; dx++) {
        const x = p.x + dx;
        const y = p.y + dy;
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (sea[i]) continue;
        if (centerRegion >= 0 && regionIdAt(x, y) !== centerRegion) continue;
        const d = Math.hypot(dx, dy);
        if (d > radius) continue;
        const dev = developable?.[i] || urbanSiteSuitability({ x, y });
        if (dev < 0.045) continue;
        const radial = clamp(1 - d / Math.max(1, radius));
        const coreDensity = 2600 + 6200 * Math.pow(radial, 1.55);
        const suburbanDensity = 650 + 2600 * Math.pow(radial, 0.85);
        const terrainMultiplier = clamp(0.58 + plain[i] * 0.28 + basinField[i] * 0.18 + coastalLowland[i] * 0.16 + valleyField[i] * 0.14 - slope[i] * 0.44 - ridgeField[i] * 0.20, 0.28, 1.22);
        capacity += dev * (coreDensity * 0.34 + suburbanDensity * 0.66) * terrainMultiplier * densityBias;
      }
    }
    const site = urbanSiteSuitability(p);
    const hardSitePenalty = site < 0.18 || elevation[indexOf(p.x, p.y)] > 0.68 || slope[indexOf(p.x, p.y)] > 0.78 || ridgeField[indexOf(p.x, p.y)] > 0.76;
    if (hardSitePenalty) capacity = Math.min(capacity, 95000);
    else if (site < 0.28) capacity = Math.min(capacity, 220000);
    return Math.max(26000, Math.round(capacity / 1000) * 1000);
  }

  function cityPopulationCap(p) {
    const pop = p?.population || 0;
    const radius = p?.isPrefecturalCapital ? 34 : p?.isRegionalCapital ? 29 : pop >= 650000 ? 30 : pop >= 250000 ? 24 : pop >= 90000 ? 18 : 13;
    const bias = p?.isPrefecturalCapital ? 1.22 : p?.isRegionalCapital ? 1.12 : 1.0;
    return estimateUrbanCapacity(p, radius, bias);
  }

  function regionIdAt(x, y) {
    const i = indexOf(x, y);
    if (prefectureMask[i]) return 0;
    const id = prefectureRegionId?.[i];
    return id !== undefined && id >= 0 ? id : -1;
  }

  function inFocusedPrefecture(point) {
    if (!point || !inside(point.x, point.y)) return false;
    const i = indexOf(point.x, point.y);
    return Boolean(prefectureMask[i] && !sea[i]);
  }

  function isHumanRegionCell(i) {
    return !sea[i] && (prefectureMask[i] || ((prefectureRegionId?.[i] ?? -1) >= 0));
  }

  function inHumanRegion(point) {
    return Boolean(point && inside(point.x, point.y) && isHumanRegionCell(indexOf(point.x, point.y)));
  }

  function sameGeneratedRegion(a, b) {
    if (!a || !b) return false;
    const ar = regionIdAt(a.x, a.y);
    const br = regionIdAt(b.x, b.y);
    return ar >= 0 && ar === br;
  }

  function pathTouchesHumanRegion(path) {
    return Boolean(path?.some(([x, y]) => inside(x, y) && isHumanRegionCell(indexOf(x, y))));
  }

  function populationDensityProxyForCapital(i) {
    return settlementScore[i] * 0.18 + marketScore[i] * 0.12;
  }

  function buildRegionLandStats() {
    const areaByRegion = new Map();
    for (let i = 0; i < SIZE; i++) {
      if (sea[i]) continue;
      const regionId = prefectureMask[i] ? 0 : ((prefectureRegionId?.[i] ?? -1) >= 0 ? prefectureRegionId[i] : -1);
      if (regionId < 0) continue;
      areaByRegion.set(regionId, (areaByRegion.get(regionId) || 0) + 1);
    }
    return areaByRegion;
  }

  function fallbackCapitalCandidate(regionId = 0) {
    const inRegion = (p) => p && regionIdAt(p.x, p.y) === regionId && !sea[indexOf(p.x, p.y)];
    const pools = [...markets, ...ports, ...villages].filter(inRegion);
    let best = null;
    let bestScore = -INF;
    for (const p of pools) {
      const i = indexOf(p.x, p.y);
      const score = urbanSiteSuitability(p) * 1.6 + plain[i] * 0.32 + populationDensityProxyForCapital(i) + (p.kind?.includes("Port") ? 0.18 : 0) + (p.score || 0);
      if (score > bestScore) { bestScore = score; best = p; }
    }
    if (best) return { ...best, kind: "Market City", score: bestScore, regionId };

    for (let y = 4; y < MAP_H - 4; y++) {
      for (let x = 4; x < MAP_W - 4; x++) {
        const i = indexOf(x, y);
        if (regionIdAt(x, y) !== regionId || sea[i]) continue;
        const score = plain[i] * 0.72 + agriculture[i] * 0.24 + basinField[i] * 0.18 + coastalLowland[i] * 0.14 - slope[i] * 0.72 - ridgeField[i] * 0.32;
        if (score > bestScore) { bestScore = score; best = { x, y, score, kind: "Market City", regionId }; }
      }
    }
    return best;
  }

  function buildCityCandidatePool(castleTowns) {
    const historicalCandidates = [
      ...castleTowns.map((p) => ({ ...p, score: p.score + 0.4 })),
      ...ports.map((p) => ({ ...p, kind: "Port Town", score: p.score + 0.28 })),
      ...markets.map((p) => ({ ...p, kind: "Market City", score: p.score + 0.12 })),
    ];
    const terrainCandidates = [];
    for (let y = 4; y < MAP_H - 4; y++) {
      for (let x = 4; x < MAP_W - 4; x++) {
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const regionId = regionIdAt(x, y);
        if (regionId < 0) continue;
        const localCapacity = developable[i] * 0.40 + townSuitability[i] * 0.34 + marketScore[i] * 0.16 + confluenceField[i] * 0.08;
        const score = localCapacity + settlementScore[i] * 0.22 + plain[i] * 0.12 + basinField[i] * 0.12 + coastalLowland[i] * 0.10 + valleyField[i] * 0.10 - slope[i] * 0.30 - ridgeField[i] * 0.16 + hash2(x, y, seed + 1666) * 0.035;
        if (score > 0.29) terrainCandidates.push({ x, y, score, kind: "City Site", regionId });
      }
    }
    return [...historicalCandidates, ...terrainCandidates].map((p) => {
      const i = indexOf(p.x, p.y);
      const suitability = urbanSiteSuitability(p);
      return {
        ...p,
        regionId: p.regionId ?? regionIdAt(p.x, p.y),
        urbanSuitability: suitability,
        score: p.score + suitability * 0.72 - slope[i] * 0.20 - ridgeField[i] * 0.16 - Math.max(0, elevation[i] - 0.58) * 0.78,
      };
    }).filter((p) => p.regionId >= 0 && (p.urbanSuitability >= 0.10 || p.kind === "Castle Town"));
  }

  function cityQuotaForRegion(regionId, area) {
    const st = settlementRegionStats.get(regionId) || { developable: area, town: 0, valley: 0, coast: 0 };
    const base = 0.7 + Math.sqrt(Math.max(1, st.developable || area)) / 64 + (st.town || 0) / 360 + rand(seed, 1650 + regionId * 19) * 0.85;
    const selectedBonus = regionId === 0 ? 2.5 : 0;
    const min = regionId === 0 ? 5 : (st.developable > 90 ? 1 : 0);
    const max = regionId === 0 ? 8 : (st.developable > 760 ? 3 : st.developable > 320 ? 2 : 1);
    return Math.round(clamp(base + selectedBonus, min, max));
  }

  function makeModernCity(candidate, localIndex, regionId, candidateCount) {
    const selectedRegion = regionId === 0;
    const isRegionCapital = localIndex === 0;
    const rank = isRegionCapital
      ? (selectedRegion ? "Prefectural Capital" : "Regional Capital")
      : localIndex < 4 ? "Regional Center" : "Small City";
    const r = rand(seed, 1600 + regionId * 97 + localIndex * 13 + candidate.x * 3 + candidate.y);
    const rawScale = Math.pow(1 - localIndex / Math.max(1, candidateCount + 1), 1.55) * 0.58 + Math.pow(r, 3.4) * 0.42;
    const rankBase = rank === "Prefectural Capital" ? 420000 : rank === "Regional Capital" ? 260000 : rank === "Regional Center" ? 115000 : 26000;
    const rankSpread = rank === "Prefectural Capital" ? 1450000 : rank === "Regional Capital" ? 820000 : rank === "Regional Center" ? 520000 : 185000;
    const pi = indexOf(candidate.x, candidate.y);
    const suitability = candidate.urbanSuitability ?? urbanSiteSuitability(candidate);
    const geographyBoost = clamp(plain[pi] * 0.34 + agriculture[pi] * 0.18 + basinField[pi] * 0.2 + coastalLowland[pi] * 0.18 + valleyField[pi] * 0.12 + suitability * 0.24 + (candidate.kind === "Port Town" ? 0.22 : 0));
    const rawPopulation = Math.round((rankBase + rankSpread * Math.pow(rawScale + geographyBoost * 0.18, 1.75)) / 1000) * 1000;
    const capacityRadius = rank === "Prefectural Capital" ? 34 : rank === "Regional Capital" ? 29 : rank === "Regional Center" ? 22 : 15;
    const capacityBias = rank === "Prefectural Capital" ? 1.22 : rank === "Regional Capital" ? 1.10 : 1.0;
    const capacityPop = estimateUrbanCapacity(candidate, capacityRadius, capacityBias);
    const floorPop = rank === "Prefectural Capital" ? 260000 : rank === "Regional Capital" ? 90000 : rank === "Regional Center" ? 48000 : 18000;
    const limitedPopulation = Math.min(rawPopulation, capacityPop * 1.08);
    const population = Math.round((capacityPop >= floorPop ? Math.max(floorPop, limitedPopulation) : limitedPopulation) / 1000) * 1000;
    const footprintCells = Math.max(10, population / (rank === "Small City" ? 2200 : rank === "Regional Center" ? 2800 : 3600));
    const urbanRadius = clamp(6.5 + Math.sqrt(footprintCells / Math.PI) * 1.75 + (isRegionCapital ? (selectedRegion ? 2.7 : 1.6) : rank === "Regional Center" ? 1.1 : 0), 7, selectedRegion ? 34 : 30);
    const coreRadius = clamp(2.2 + Math.sqrt(population) / 340, 2.6, selectedRegion ? 9 : 8);
    const urbanWeight = clamp(0.74 + Math.log10(Math.max(10000, population)) * 0.36 + Math.sqrt(Math.max(1, footprintCells)) / 130, 1.15, selectedRegion ? 3.15 : 2.8);
    return {
      ...candidate,
      population,
      urbanRadius,
      coreRadius,
      urbanWeight,
      rank,
      regionId,
      isRegionalCapital: isRegionCapital,
      isPrefecturalCapital: false,
      kind: selectedRegion && isRegionCapital ? "Prefectural Capital" : isRegionCapital ? "Regional Capital" : candidate.kind || "City",
    };
  }

  function generateUrbanCentersByRegion() {
    const castleTowns = castles.map((c) => ({ x: c.x, y: c.y, score: c.score + 0.45, kind: "Castle Town" }));
    const cityCandidates = buildCityCandidatePool(castleTowns);
    const areaByRegion = buildRegionLandStats();
    const regionIds = [...areaByRegion.keys()].sort((a, b) => a - b);
    const modernCities = [];
    const debug = { regionCount: regionIds.length, citiesByRegion: {}, candidateCountByRegion: {} };

    for (const regionId of regionIds) {
      const area = areaByRegion.get(regionId) || 0;
      if (area < 140) continue;
      const candidates = cityCandidates.filter((p) => p.regionId === regionId);
      const quota = cityQuotaForRegion(regionId, area);
      const threshold = (regionId === 0 ? 0.33 : 0.35) + rand(seed, 1062 + regionId * 31) * 0.12;
      let picked = pickEntities(candidates, {
        max: quota,
        minDistance: regionId === 0 ? 9 : 10,
        threshold,
        seed: seed + 1060 + regionId * 101,
      });

      if (picked.length === 0 || (regionId === 0 && !picked.some((city) => prefectureMask[indexOf(city.x, city.y)]))) {
        const fallback = fallbackCapitalCandidate(regionId);
        if (fallback) picked = [fallback, ...picked];
      }

      picked = picked
        .filter((city, index, arr) => arr.findIndex((other) => other.x === city.x && other.y === city.y) === index)
        .sort((a, b) => {
          const ai = indexOf(a.x, a.y);
          const bi = indexOf(b.x, b.y);
          const aScore = (a.score || 0) + urbanSiteSuitability(a) * 0.85 + plain[ai] * 0.16 - slope[ai] * 0.24;
          const bScore = (b.score || 0) + urbanSiteSuitability(b) * 0.85 + plain[bi] * 0.16 - slope[bi] * 0.24;
          return bScore - aScore;
        })
        .slice(0, quota);

      const regionCities = picked.map((city, localIndex) => makeModernCity(city, localIndex, regionId, candidates.length));
      debug.citiesByRegion[regionId] = regionCities.length;
      debug.candidateCountByRegion[regionId] = candidates.length;
      modernCities.push(...regionCities);
    }

    if (!modernCities.some((city) => prefectureMask[indexOf(city.x, city.y)])) {
      const fallback = fallbackCapitalCandidate(0);
      if (fallback) modernCities.unshift(makeModernCity(fallback, 0, 0, 1));
    }

    let selectedCapitalIndex = -1;
    let selectedCapitalScore = -INF;
    for (let i = 0; i < modernCities.length; i++) {
      const city = modernCities[i];
      const ci = indexOf(city.x, city.y);
      if (!prefectureMask[ci] || sea[ci]) continue;
      const suitability = urbanSiteSuitability(city);
      const score = suitability * 900000 + (city.population || 0) * 0.55 + (city.score || 0) * 120000 - slope[ci] * 180000 - Math.max(0, elevation[ci] - 0.58) * 360000;
      if (score > selectedCapitalScore) { selectedCapitalScore = score; selectedCapitalIndex = i; }
    }

    if (selectedCapitalIndex >= 0) {
      const cap = modernCities[selectedCapitalIndex];
      modernCities[selectedCapitalIndex] = {
        ...cap,
        rank: "Prefectural Capital",
        kind: "Prefectural Capital",
        isPrefecturalCapital: true,
        isRegionalCapital: true,
        population: Math.round(Math.min(Math.max(cap.population || 0, 620000), estimateUrbanCapacity(cap, 34, 1.25) * 1.08) / 1000) * 1000,
        urbanRadius: Math.max(cap.urbanRadius || 0, 18),
        coreRadius: Math.max(cap.coreRadius || 0, 5.5),
        urbanWeight: Math.max(cap.urbanWeight || 0, 2.15),
      };
    }

    for (let i = 0; i < modernCities.length; i++) {
      if (i !== selectedCapitalIndex) modernCities[i] = { ...modernCities[i], isPrefecturalCapital: false };
    }

    modernCities.sort((a, b) => {
      const aIn = prefectureMask[indexOf(a.x, a.y)] ? 1 : 0;
      const bIn = prefectureMask[indexOf(b.x, b.y)] ? 1 : 0;
      if (a.isPrefecturalCapital !== b.isPrefecturalCapital) return a.isPrefecturalCapital ? -1 : 1;
      if (aIn !== bIn) return bIn - aIn;
      if (a.isRegionalCapital !== b.isRegionalCapital) return a.isRegionalCapital ? -1 : 1;
      return ((b.population || 0) + (b.score || 0) * 90000) - ((a.population || 0) + (a.score || 0) * 90000);
    });

    return { castleTowns, modernCities, urbanHierarchyDebug: debug };
  }

  const { castleTowns, modernCities, urbanHierarchyDebug } = generateUrbanCentersByRegion();

  const capital = modernCities.find((city) => city.isPrefecturalCapital && prefectureMask[indexOf(city.x, city.y)]) || modernCities.find((city) => prefectureMask[indexOf(city.x, city.y)]) || markets.find((p) => prefectureMask[indexOf(p.x, p.y)]) || ports.find((p) => prefectureMask[indexOf(p.x, p.y)]) || { x: Math.floor(MAP_W / 2), y: Math.floor(MAP_H / 2), score: 1, population: 0, urbanRadius: 12, coreRadius: 4, urbanWeight: 1, isPrefecturalCapital: true };

  const generatedRegionIdsForTransport = [...new Set([
    ...modernCities,
    ...ports,
    ...markets,
    ...castles,
    ...villages,
  ].map((p) => regionIdAt(p.x, p.y)).filter((id) => id >= 0))].sort((a, b) => a - b);

  function regionNodes(nodes, regionId) {
    return nodes.filter((p) => p && regionIdAt(p.x, p.y) === regionId && !sea[indexOf(p.x, p.y)]);
  }

  function primaryNodeForRegion(regionId) {
    if (regionId === 0) return capital;
    return modernCities.find((city) => city.regionId === regionId && city.isRegionalCapital)
      || modernCities.find((city) => regionIdAt(city.x, city.y) === regionId)
      || markets.find((p) => regionIdAt(p.x, p.y) === regionId)
      || ports.find((p) => regionIdAt(p.x, p.y) === regionId)
      || villages.find((p) => regionIdAt(p.x, p.y) === regionId)
      || null;
  }

  function buildRegionalTransportLinks(nodes, makeConfig, seedSalt = 0) {
    const links = [];
    for (const regionId of generatedRegionIdsForTransport) {
      const primary = primaryNodeForRegion(regionId);
      const localNodes = uniqueByPosition([primary, ...regionNodes(nodes, regionId)]);
      if (localNodes.length < 2) continue;
      const config = makeConfig(regionId, localNodes.length);
      const localLinks = buildHierarchicalLinks(localNodes, {
        ...config,
        seedOffset: (config.seedOffset || 0) + seedSalt + regionId * 997,
      });
      links.push(...localLinks.map((link) => ({ ...link, regionId })));
    }
    return links;
  }

  function uniqueByPosition(nodes) {
    const seen = new Set();
    const out = [];
    for (const node of nodes.filter(Boolean)) {
      const key = `${node.x},${node.y}`;
      if (seen.has(key)) continue;
      seen.add(key);
      out.push(node);
    }
    return out;
  }


  const urbanFootprint = new Uint8Array(SIZE);
  const urbanCoreFootprint = new Uint8Array(SIZE);
  const oldUrbanFootprint = new Uint8Array(SIZE);
  const ruralSettlementFootprint = new Uint8Array(SIZE);

  function canUrbanizeCell(i, centerRegion) {
    if (sea[i] || !isHumanRegionCell(i)) return false;
    const [x, y] = xyOf(i);
    if (centerRegion >= 0 && regionIdAt(x, y) !== centerRegion) return false;
    if (elevation[i] > 0.78 || slope[i] > 0.68 || ridgeField[i] > 0.78) return false;
    return (developable[i] || 0) > 0.035 || (townSuitability[i] || 0) > 0.18;
  }

  function growSettlementFootprint(center, {
    targetCells,
    coreCells = 0,
    maxRadius,
    mask,
    coreMask = null,
    oldMask = null,
    seedOffset = 0,
    minSupport = -0.08,
    allowSmallMountainValleys = false,
  }) {
    if (!center || !inside(center.x, center.y)) return 0;
    const start = indexOf(center.x, center.y);
    const centerRegion = regionIdAt(center.x, center.y);
    if (sea[start] || !isHumanRegionCell(start)) return 0;
    const heap = new MinHeap();
    const queued = new Set([start]);
    const selected = [];
    heap.push({ i: start, f: -10 });
    const hardLimit = Math.ceil(maxRadius + 2);

    while (heap.length > 0 && selected.length < targetCells) {
      const cur = heap.pop();
      if (!cur) break;
      const i = cur.i;
      const [x, y] = xyOf(i);
      const d = Math.hypot(x - center.x, y - center.y);
      if (d > maxRadius) continue;
      if (!canUrbanizeCell(i, centerRegion)) {
        if (!(allowSmallMountainValleys && valleySettlement[i] > 0.34 && slope[i] < 0.58 && elevation[i] < 0.72)) continue;
      }
      const support =
        developable[i] * 1.12 +
        townSuitability[i] * 0.36 +
        plain[i] * 0.16 +
        basinField[i] * 0.14 +
        coastalLowland[i] * 0.13 +
        valleyField[i] * 0.12 -
        slope[i] * 0.62 -
        ridgeField[i] * 0.26 -
        Math.max(0, elevation[i] - 0.60) * 0.70 -
        d / Math.max(1, maxRadius) * 0.12;
      if (support < minSupport && selected.length > 0) continue;
      selected.push(i);
      mask[i] = 1;
      if (oldMask && selected.length <= Math.max(2, Math.round(targetCells * 0.18))) oldMask[i] = 1;
      if (coreMask && selected.length <= coreCells) coreMask[i] = 1;

      for (const [nx, ny] of neighbors8(x, y)) {
        if (Math.abs(nx - center.x) > hardLimit || Math.abs(ny - center.y) > hardLimit) continue;
        const ni = indexOf(nx, ny);
        if (queued.has(ni) || sea[ni]) continue;
        const nd = Math.hypot(nx - center.x, ny - center.y);
        if (nd > maxRadius + 1) continue;
        const sameRegion = centerRegion < 0 || regionIdAt(nx, ny) === centerRegion;
        if (!sameRegion) continue;
        const terrainCost =
          nd / Math.max(1, maxRadius) * 0.92 +
          slope[ni] * 1.35 +
          ridgeField[ni] * 0.70 +
          Math.max(0, elevation[ni] - 0.58) * 1.10 +
          floodplain[ni] * 0.08 -
          developable[ni] * 1.42 -
          townSuitability[ni] * 0.40 -
          plain[ni] * 0.18 -
          valleyField[ni] * 0.14 -
          coastalLowland[ni] * 0.12 +
          hash2(nx, ny, seed + seedOffset) * 0.055;
        queued.add(ni);
        heap.push({ i: ni, f: terrainCost });
      }
    }
    return selected.length;
  }

  for (let n = 0; n < modernCities.length; n++) {
    const city = modernCities[n];
    const pop = city.population || 40000;
    const densityPerCell = city.isPrefecturalCapital ? 3600 : city.isRegionalCapital ? 3300 : pop >= 180000 ? 2850 : 2200;
    const targetCells = Math.round(clamp(pop / densityPerCell, city.isRegionalCapital ? 22 : 8, city.isPrefecturalCapital ? 620 : city.isRegionalCapital ? 420 : 190));
    const coreCells = Math.round(clamp(pop / (city.isRegionalCapital ? 27000 : 36000), pop >= 120000 ? 3 : 1, city.isPrefecturalCapital ? 58 : city.isRegionalCapital ? 42 : 18));
    const maxRadius = clamp(Math.max(city.urbanRadius || 8, Math.sqrt(targetCells / Math.PI) * 2.25), 7, city.isPrefecturalCapital ? 35 : city.isRegionalCapital ? 30 : 22);
    const made = growSettlementFootprint(city, {
      targetCells,
      coreCells,
      maxRadius,
      mask: urbanFootprint,
      coreMask: urbanCoreFootprint,
      oldMask: oldUrbanFootprint,
      seedOffset: 6000 + n * 31,
      minSupport: city.isRegionalCapital ? -0.12 : -0.06,
    });
    city.urbanFootprintCells = made;
    city.coreFootprintCells = Math.min(coreCells, made);
    city.urbanRadius = Math.max(city.urbanRadius || 0, clamp(Math.sqrt(Math.max(1, made) / Math.PI) * 1.85, 6, city.isRegionalCapital ? 34 : 26));
    city.coreRadius = Math.max(city.coreRadius || 0, clamp(Math.sqrt(Math.max(1, city.coreFootprintCells) / Math.PI) * 1.25, 2, 9));
  }

  for (let n = 0; n < markets.length; n++) {
    const town = markets[n];
    const pop = town.population || 9000;
    const targetCells = Math.round(clamp(pop / 3600, 2, 13));
    const maxRadius = clamp(3.5 + Math.sqrt(targetCells) * 1.5, 4, 10);
    const made = growSettlementFootprint(town, {
      targetCells,
      coreCells: 0,
      maxRadius,
      mask: oldUrbanFootprint,
      oldMask: oldUrbanFootprint,
      seedOffset: 7000 + n * 23,
      minSupport: -0.14,
      allowSmallMountainValleys: true,
    });
    town.urbanFootprintCells = made;
  }

  for (const village of villages) {
    const i = indexOf(village.x, village.y);
    if (!sea[i]) ruralSettlementFootprint[i] = 1;
    if ((village.population || 0) > 2600 || valleySettlement[i] > 0.42) {
      for (const [nx, ny] of neighbors8(village.x, village.y)) {
        const ni = indexOf(nx, ny);
        if (!sea[ni] && regionIdAt(nx, ny) === regionIdAt(village.x, village.y) && (ruralSuitability[ni] > 0.22 || valleySettlement[ni] > 0.32)) ruralSettlementFootprint[ni] = 1;
      }
    }
  }

  const populationDensity = new Float32Array(SIZE);
  let maxPopulationDensity = 0;
  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      let density = 0;
      if (urbanCoreFootprint[i]) density += 1.35;
      else if (oldUrbanFootprint[i]) density += 0.72;
      else if (urbanFootprint[i]) density += 0.52;
      else if (ruralSettlementFootprint[i]) density += 0.12;
      for (const city of modernCities) {
        const populationScale = clamp((Math.log10(Math.max(10000, city.population || 10000)) - 4) / 2.25, 0.12, 1.60);
        const d = Math.hypot(city.x - x, city.y - y);
        const urbanR = Math.max(5, city.urbanRadius || 11);
        const coreR = Math.max(2.4, city.coreRadius || 4);
        const inFootprint = urbanFootprint[i] ? 1 : 0;
        density += populationScale * (inFootprint ? 0.92 : 0.82) / (1 + Math.pow(d / urbanR, 2.28));
        density += populationScale * 0.72 * Math.exp(-(d * d) / (coreR * coreR * 2.05));
      }
      for (const market of markets) {
        const d = Math.hypot(market.x - x, market.y - y);
        density += 0.18 / (1 + Math.pow(d / 6.2, 2.0));
      }
      for (const village of villages) {
        const d = Math.hypot(village.x - x, village.y - y);
        density += 0.038 / (1 + Math.pow(d / 3.5, 2));
      }
      const terrainFactor = urbanFootprint[i] || oldUrbanFootprint[i]
        ? clamp(0.68 + plain[i] * 0.32 + agriculture[i] * 0.08 + basinField[i] * 0.16 + coastalLowland[i] * 0.12 + valleyField[i] * 0.10 - slope[i] * 0.46 - ridgeField[i] * 0.18, 0.34, 1.20)
        : clamp(0.42 + plain[i] * 0.58 + agriculture[i] * 0.18 + basinField[i] * 0.22 + coastalLowland[i] * 0.18 + valleyField[i] * 0.14 - slope[i] * 0.92 - ridgeField[i] * 0.42 - Math.max(0, elevation[i] - 0.58) * 0.90, 0.018, 1.18);
      density *= terrainFactor;
      populationDensity[i] = density;
      if (density > maxPopulationDensity) maxPopulationDensity = density;
    }
  }
  if (maxPopulationDensity > 0) {
    for (let i = 0; i < SIZE; i++) populationDensity[i] = clamp(populationDensity[i] / maxPopulationDensity);
  }

  function densityValue(x, y) {
    return populationDensity[indexOf(x, y)] || 0;
  }

  function midDensityAffinity(x, y) {
    const d = densityValue(x, y);
    return clamp(1 - Math.abs(d - 0.38) / 0.38);
  }

  function transportTier(p) {
    const pop = p?.population || 0;
    if (p?.isPrefecturalCapital || p?.rank === "Prefectural Capital") return 0;
    if (pop >= 900000) return 1;
    if (pop >= 360000) return 2;
    if (pop >= 180000) return 3;
    if (pop >= 90000) return 4;
    if (p?.portClass === "major") return 2;
    if (p?.portClass === "regional") return 3;
    if (p?.kind === "Market Town") return 4;
    if (p?.kind?.includes("Castle")) return 5;
    return 6;
  }

  function nodeKey(p) {
    return `${p.x},${p.y}`;
  }

  function addUniqueNode(list, node) {
    if (!node) return;
    const key = nodeKey(node);
    if (!list.some((p) => nodeKey(p) === key)) list.push(node);
  }

  function pointLineDistanceXY(x, y, a, b) {
    const vx = b.x - a.x;
    const vy = b.y - a.y;
    const len2 = vx * vx + vy * vy;
    if (len2 <= 0.0001) return Math.hypot(x - a.x, y - a.y);
    const t = clamp(((x - a.x) * vx + (y - a.y) * vy) / len2, 0, 1);
    return Math.hypot(x - (a.x + vx * t), y - (a.y + vy * t));
  }

  function segmentProgressXY(x, y, a, b) {
    const vx = b.x - a.x;
    const vy = b.y - a.y;
    const len2 = vx * vx + vy * vy;
    if (len2 <= 0.0001) return 0;
    return clamp(((x - a.x) * vx + (y - a.y) * vy) / len2, 0, 1);
  }

  function cellTransportCorridorScore(x, y, mode = "road") {
    if (!inside(x, y)) return -INF;
    const i = indexOf(x, y);
    if (sea[i]) return -INF;
    const barrier = mountainBarrierPenalty(x, y, mode === "express" ? "express" : mode === "rail" ? "rail" : "road");
    if (barrier >= INF) return -INF;
    const density = densityValue(x, y);
    const midDensity = midDensityAffinity(x, y);
    const lowland = plain[i] * 0.42 + basinField[i] * 0.24 + coastalLowland[i] * 0.22 + valleyField[i] * 0.34 + agriculture[i] * 0.08;
    const terrainCost = slope[i] * (mode === "rail" ? 1.35 : mode === "express" ? 1.18 : 1.0)
      + ridgeField[i] * 0.62
      + Math.max(0, elevation[i] - (mode === "rail" ? 0.50 : 0.56)) * 1.18
      + barrier * (mode === "rail" ? 0.0048 : mode === "express" ? 0.0036 : 0.0038);
    if (mode === "rail") return density * 1.95 + lowland + coastalLowland[i] * 0.22 + valleyField[i] * 0.24 - terrainCost;
    if (mode === "express") return density * 0.82 + midDensity * 1.10 + lowland * 0.62 - Math.max(0, density - 0.88) * 1.35 - terrainCost;
    return density * 1.28 + midDensity * 0.28 + lowland + valleyField[i] * 0.18 - terrainCost;
  }

  function sampleCorridorValue(a, b, mode = "road", samples = 14) {
    let total = 0;
    let count = 0;
    for (let k = 1; k < samples; k++) {
      const t = k / samples;
      const x = Math.round(a.x + (b.x - a.x) * t);
      const y = Math.round(a.y + (b.y - a.y) * t);
      if (!inside(x, y)) continue;
      const score = cellTransportCorridorScore(x, y, mode);
      if (score <= -INF / 2) continue;
      total += score;
      count++;
    }
    return count ? total / count : -3;
  }

  function pairTransportScore(a, b, mode = "road") {
    const d = Math.hypot(a.x - b.x, a.y - b.y);
    if (d < 4) return -INF;
    const demand = Math.sqrt(Math.max(0.01, transportDemand(a)) * Math.max(0.01, transportDemand(b)));
    const hierarchyDelta = Math.max(0, transportTier(b) - transportTier(a));
    const corridor = sampleCorridorValue(a, b, mode);
    const sameCorridor = sameCorridorAffinity(a, b);
    const distancePenalty = mode === "express" ? d / 58 : mode === "rail" ? d / 48 : d / 42;
    const hierarchyBonus = hierarchyDelta * (mode === "express" ? 0.16 : 0.10);
    const portBonus = (a.portClass || b.portClass) ? (mode === "rail" ? 0.34 : mode === "express" ? 0.18 : 0.26) : 0;
    return demand * (mode === "express" ? 1.08 : mode === "rail" ? 1.18 : 1.0) + corridor * 0.72 + sameCorridor + hierarchyBonus + portBonus - distancePenalty;
  }

  function buildHierarchicalLinks(nodes, { mode = "road", maxLinks = 10, extraLinks = 3, minDistance = 10, maxDistance = 70, maxDegree = 3, seedOffset = 0 } = {}) {
    const unique = [];
    const seen = new Set();
    for (const node of nodes.filter(Boolean)) {
      const i = indexOf(node.x, node.y);
      if (!inside(node.x, node.y) || sea[i]) continue;
      const key = nodeKey(node);
      if (seen.has(key)) continue;
      seen.add(key);
      unique.push({ ...node, transportTier: transportTier(node), demand: transportDemand(node) });
    }
    const ranked = unique.sort((a, b) => a.transportTier - b.transportTier || b.demand - a.demand || b.score - a.score);
    const degree = new Map();
    const usedPairs = new Set();
    const links = [];

    function pairKey(a, b) {
      const ak = nodeKey(a);
      const bk = nodeKey(b);
      return ak < bk ? `${ak}|${bk}` : `${bk}|${ak}`;
    }
    function tryAdd(a, b, force = false) {
      if (!a || !b || nodeKey(a) === nodeKey(b)) return false;
      const d = Math.hypot(a.x - b.x, a.y - b.y);
      if (d < minDistance || d > maxDistance) return false;
      const key = pairKey(a, b);
      if (usedPairs.has(key)) return false;
      if (!force && (getDegree(degree, a) >= maxDegree || getDegree(degree, b) >= maxDegree)) return false;
      usedPairs.add(key);
      incrementDegree(degree, a);
      incrementDegree(degree, b);
      links.push({ a, b, score: pairTransportScore(a, b, mode), distance: d });
      return true;
    }

    for (let i = 1; i < ranked.length && links.length < maxLinks; i++) {
      const child = ranked[i];
      const parentCandidates = ranked.slice(0, i)
        .filter((parent) => transportTier(parent) <= transportTier(child) && Math.hypot(parent.x - child.x, parent.y - child.y) <= maxDistance)
        .map((parent) => ({ parent, score: pairTransportScore(parent, child, mode) - getDegree(degree, parent) * 0.16 - Math.max(0, getDegree(degree, child) - 1) * 0.22 }))
        .sort((a, b) => b.score - a.score);
      if (parentCandidates[0]) tryAdd(parentCandidates[0].parent, child, true);
    }

    const candidates = [];
    for (let i = 0; i < ranked.length; i++) {
      for (let j = i + 1; j < ranked.length; j++) {
        const a = ranked[i];
        const b = ranked[j];
        const d = Math.hypot(a.x - b.x, a.y - b.y);
        if (d < minDistance || d > maxDistance) continue;
        candidates.push({ a, b, score: pairTransportScore(a, b, mode) + hash2(a.x + b.x, a.y + b.y, seed + seedOffset + i * 31 + j * 37) * 0.05 });
      }
    }
    candidates.sort((a, b) => b.score - a.score);
    let addedExtra = 0;
    for (const c of candidates) {
      if (links.length >= maxLinks || addedExtra >= extraLinks) break;
      if (tryAdd(c.a, c.b)) addedExtra++;
    }
    return links;
  }

  function pickCorridorWaypoints(a, b, mode = "road", maxCount = 2) {
    const d = Math.hypot(a.x - b.x, a.y - b.y);
    if (d < 20) return [];
    const width = mode === "express" ? 10.5 : mode === "rail" ? 8.5 : 9.5;
    const minProgress = 0.18;
    const maxProgress = 0.82;
    const candidates = [];
    const minX = Math.max(1, Math.floor(Math.min(a.x, b.x) - width - 3));
    const maxX = Math.min(MAP_W - 2, Math.ceil(Math.max(a.x, b.x) + width + 3));
    const minY = Math.max(1, Math.floor(Math.min(a.y, b.y) - width - 3));
    const maxY = Math.min(MAP_H - 2, Math.ceil(Math.max(a.y, b.y) + width + 3));
    for (let y = minY; y <= maxY; y++) {
      for (let x = minX; x <= maxX; x++) {
        const progress = segmentProgressXY(x, y, a, b);
        if (progress < minProgress || progress > maxProgress) continue;
        const lineD = pointLineDistanceXY(x, y, a, b);
        if (lineD > width) continue;
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const cellScore = cellTransportCorridorScore(x, y, mode);
        if (cellScore <= -INF / 2) continue;
        const centerBias = -Math.abs(progress - 0.5) * 0.15;
        const linePenalty = lineD / width * (mode === "express" ? 0.42 : 0.32);
        const density = densityValue(x, y);
        const densityGate = mode === "express" ? midDensityAffinity(x, y) * 0.22 : density * 0.20;
        const score = cellScore + densityGate + centerBias - linePenalty + hash2(x, y, seed + 6400 + mode.length * 101) * 0.05;
        candidates.push({ x, y, score, progress, kind: `${mode} corridor waypoint` });
      }
    }
    if (!candidates.length) return [];
    const count = Math.min(maxCount, d > 62 ? 2 : 1);
    return pickEntities(candidates, {
      max: count,
      minDistance: Math.max(7, Math.floor(d / 4.2)),
      threshold: mode === "express" ? -0.42 : -0.30,
      seed: seed + 6500 + Math.round(a.x * 13 + a.y * 17 + b.x * 19 + b.y * 23),
    }).sort((p, q) => p.progress - q.progress);
  }

  function makeDensityAwareTransportCost(baseCost, mode, guidePoints = []) {
    return (x, y, cx, cy) => {
      const base = baseCost(x, y, cx, cy);
      if (base >= INF) return base;
      const density = densityValue(x, y);
      const midDensity = midDensityAffinity(x, y);
      const cityDistance = distanceToNearest(modernCities, x, y);
      let guidePull = 0;
      for (const p of guidePoints) guidePull = Math.max(guidePull, 1 / (1 + Math.hypot(x - p.x, y - p.y) / 5.8));
      if (mode === "rail") {
        const lowDemandPenalty = Math.max(0, 0.08 - density) * 2.4;
        return Math.max(0.30, base + lowDemandPenalty - density * 0.72 - guidePull * 0.34);
      }
      if (mode === "express") {
        const coreAvoid = cityDistance < 2.6 ? 5.0 : cityDistance < 5.8 ? 1.4 : 0;
        const lowDemandPenalty = Math.max(0, 0.10 - density) * 3.2;
        return Math.max(0.42, base + lowDemandPenalty + coreAvoid - midDensity * 0.36 - Math.min(density, 0.72) * 0.18 - guidePull * 0.20);
      }
      const lowDemandPenalty = Math.max(0, 0.06 - density) * 1.5;
      return Math.max(0.28, base + lowDemandPenalty - density * 0.44 - midDensity * 0.12 - guidePull * 0.26);
    };
  }

  function routeThroughTransportCorridor(a, b, mode, baseCost, existingPaths, hubs, avoidPoints, options = {}) {
    const start = routePoint(a, mode, a.x * 31 + a.y * 37 + (options.salt || 0));
    const goal = routePoint(b, mode, b.x * 31 + b.y * 37 + 17 + (options.salt || 0));
    const via = pickCorridorWaypoints(start, goal, mode, options.maxWaypoints ?? (mode === "express" ? 1 : 2));
    const terminals = [start, ...via, goal];
    const endpointSet = [start, goal, ...via];
    const guidedBaseCost = makeDensityAwareTransportCost(baseCost, mode, via);
    const path = [];
    for (let i = 0; i < terminals.length - 1; i++) {
      const from = terminals[i];
      const to = terminals[i + 1];
      const cost = makeTransportCost(
        guidedBaseCost,
        [...existingPaths, path],
        hubs,
        endpointSet,
        options.corridorRadius ?? (mode === "express" ? 5 : mode === "rail" ? 5 : 3),
        options.corridorStrength ?? (mode === "express" ? 9.4 : mode === "rail" ? 8.8 : 5.8),
        avoidPoints,
        options.avoidRadius ?? (mode === "express" ? 8.0 : mode === "rail" ? 2.5 : 3.2),
        options.avoidStrength ?? (mode === "express" ? 12.0 : mode === "rail" ? 4.2 : 5.4),
      );
      const segment = aStar(from, to, cost);
      if (segment.length < 2) return { path: [], via, start, goal };
      if (path.length) path.push(...segment.slice(1));
      else path.push(...segment);
    }
    return { path, via, start, goal };
  }

  function nearPassPoint(x, y, radius = 5) {
    return distanceToNearest(passes, x, y) <= radius;
  }

  function mountainBarrierPenalty(x, y, type = "rail") {
    const i = indexOf(x, y);
    const e = elevation[i];
    const s = slope[i];
    const pass = nearPassPoint(x, y, type === "express" ? 7 : type === "rail" ? 6 : 5);
    if (e > 0.84) return INF;
    if (pass && e > 0.80 && s > 0.16) return INF;
    if (!pass && e > 0.78) return INF;
    if (!pass && e > 0.70 && s > 0.16) return INF;
    if (!pass && e > 0.66 && s > 0.28) return INF;
    if (!pass && e > 0.72) return type === "express" ? 260 : type === "rail" ? 330 : type === "minor" ? 80 : 155;
    if (!pass && e > 0.64 && s > 0.20) return type === "express" ? 145 : type === "rail" ? 180 : type === "minor" ? 54 : 96;
    const passDiscount = pass ? (type === "minor" ? 0.35 : 0.22) : 1;
    const mountain = Math.max(0, e - 0.48);
    const steep = Math.max(0, s - 0.15);
    const typeFactor = type === "express" ? 360 : type === "rail" ? 430 : type === "minor" ? 115 : 210;
    return (mountain * mountain * typeFactor + steep * steep * 150 + ridgeField[i] * 9.5) * passDiscount;
  }

  function transportAccessPoint(node, mode = "road", salt = 0) {
    if (!node || nearMapEdge(node.x, node.y, 1) || node.kind === "External Gateway") return node;
    const minR = mode === "express" ? 6 : mode === "rail" ? 2 : 3;
    const maxR = mode === "express" ? 16 : mode === "rail" ? 6 : 8;
    let best = null;
    let bestScore = -INF;
    for (let dy = -maxR; dy <= maxR; dy++) {
      for (let dx = -maxR; dx <= maxR; dx++) {
        const d = Math.hypot(dx, dy);
        if (d < minR || d > maxR) continue;
        const x = node.x + dx;
        const y = node.y + dy;
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const barrier = mode === "express" || mode === "rail" ? mountainBarrierPenalty(x, y, mode) : mountainBarrierPenalty(x, y, "road");
        if (barrier >= INF) continue;
        const targetD = (minR + maxR) * 0.5;
        const flatness = plain[i] * 1.0 + agriculture[i] * 0.2 + valleyField[i] * 0.26 + coastalLowland[i] * 0.16 - slope[i] * 1.22 - ridgeField[i] * 0.72 - Math.max(0, elevation[i] - 0.58) * 2.35;
        const ring = -Math.abs(d - targetD) * 0.08;
        const riverPenalty = river[i] > 0.5 ? 0.45 : river[i] * 0.12;
        const density = densityValue(x, y);
        const densityAffinity = mode === "rail" ? density * 0.9 : mode === "express" ? density * 0.70 + midDensityAffinity(x, y) * 0.18 - Math.max(0, density - 0.92) * 0.35 : density * 0.24;
        const noise = hash2(x, y, seed + salt + (mode === "rail" ? 6000 : mode === "express" ? 7000 : 5000)) * 0.12;
        const score = flatness + densityAffinity + ring - riverPenalty - barrier * 0.012 + noise;
        if (score > bestScore) {
          bestScore = score;
          best = { x, y, score: node.score || 0.5, kind: `${mode} Access`, parent: node };
        }
      }
    }
    return best || node;
  }

  function routePoint(node, mode, salt = 0) {
    return transportAccessPoint(node, mode, salt);
  }

  const townAvoidNodes = [...modernCities, ...markets, ...ports];

  const urbanCenters = modernCities.map((city, n) => {
    let best = { x: city.x, y: city.y, score: city.score + 0.5 };
    let bestScore = -INF;
    const searchR = Math.max(2, Math.round(city.coreRadius));
    for (let dy = -searchR; dy <= searchR; dy++) {
      for (let dx = -searchR; dx <= searchR; dx++) {
        const x = city.x + dx;
        const y = city.y + dy;
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const d = Math.hypot(dx, dy);
        const score = plain[i] * 0.54 + agriculture[i] * 0.16 - slope[i] * 0.36 - d * 0.06 + hash2(x, y, seed + 1700 + n) * 0.07;
        if (score > bestScore) { bestScore = score; best = { x, y, score: city.score + 0.5, cityIndex: n, parent: city }; }
      }
    }
    return { ...best, kind: city.rank === "Prefectural Capital" ? "Central Business District" : "Urban Center", population: Math.round(city.population * (city.rank === "Prefectural Capital" ? 0.18 : 0.12)), insidePrefecture: Boolean(prefectureMask[indexOf(best.x, best.y)]) };
  });

  function railCost(x, y) {
    const i = indexOf(x, y);
    if (sea[i]) return INF;
    const barrier = mountainBarrierPenalty(x, y, "rail");
    if (barrier >= INF) return INF;
    const density = densityValue(x, y);
    const highPenalty = Math.max(0, elevation[i] - 0.52) * 14 + barrier;
    const riverPenalty = river[i] > 0.5 ? 1.6 : river[i] > 0.25 ? 0.7 : 0;
    return Math.max(0.42, 1 + slope[i] * 22 + highPenalty + riverPenalty + floodplain[i] * 0.28 - density * 0.88 - plain[i] * 0.28 - valleyField[i] * 0.62 - coastalLowland[i] * 0.48 + ridgeField[i] * 1.4 + normalEdgePenalty(x, y) + hash2(x, y, seed + 222) * 0.08);
  }

  const railways = [];
  const branchRailways = [];
  const railDegree = new Map();
  const railCore = [capital];
  const railHubs = [...modernCities, ...commercialPorts];

  function addRailRoute(a, b, bucket = railways) {
    const existingRails = [...railways, ...branchRailways];
    const { path } = routeThroughTransportCorridor(a, b, "rail", railCost, existingRails, railHubs, townAvoidNodes, {
      salt: a.x * 19 + a.y * 23 + b.x * 7 + b.y * 11,
      maxWaypoints: bucket === railways ? 2 : 1,
      corridorRadius: 5,
      corridorStrength: bucket === railways ? 10.8 : 8.4,
      avoidRadius: 2.4,
      avoidStrength: 4.2,
    });
    const length = pathLength(path);
    const direct = pathEndpointDistance(path);
    const overlap = pathOverlapRatio(path, existingRails, 2);
    const densityPurpose = averagePathField(path, populationDensity) * 1.22 + averagePathField(path, plain) * 0.24 + averagePathField(path, valleyField) * 0.22 + averagePathField(path, coastalLowland) * 0.16;
    const isMain = bucket === railways;
    if (path.length > 3 && direct >= (isMain ? 18 : 12) && length >= (isMain ? 22 : 14) && pathCompactness(path) < (isMain ? 3.25 : 3.5) && overlap < (isMain ? 0.32 : 0.22) && densityPurpose > (isMain ? 0.20 : 0.13)) {
      bucket.push(path);
      incrementDegree(railDegree, a);
      incrementDegree(railDegree, b);
      return true;
    }
    return false;
  }

  function addRailRouteRelaxed(a, b, bucket = railways) {
    if (!a || !b || Math.hypot(a.x - b.x, a.y - b.y) < 6) return false;
    const existingRails = [...railways, ...branchRailways, ...bucket].filter(Boolean);
    const start = routePoint(a, "rail", a.x * 83 + a.y * 89);
    const goal = routePoint(b, "rail", b.x * 97 + b.y * 101);
    const path = aStar(start, goal, makeTransportCost(railCost, existingRails, railHubs, [start, goal], 4, 7.4, townAvoidNodes, 2.0, 3.0));
    const direct = pathEndpointDistance(path);
    const densityPurpose = averagePathField(path, populationDensity) * 0.95 + averagePathField(path, plain) * 0.20 + averagePathField(path, valleyField) * 0.28 + averagePathField(path, coastalLowland) * 0.18;
    if (path.length > 4 && pathTouchesHumanRegion(path) && direct >= 8 && pathLength(path) >= 10 && pathLength(path) < 150 && pathCompactness(path) < 5.1 && pathOverlapRatio(path, existingRails, 2) < 0.58 && densityPurpose > 0.055) {
      bucket.push(path);
      incrementDegree(railDegree, a);
      incrementDegree(railDegree, b);
      return true;
    }
    return false;
  }

  function addRailRouteForced(a, b, bucket = railways) {
    if (!a || !b || Math.hypot(a.x - b.x, a.y - b.y) < 6) return false;
    const existingRails = [...railways, ...branchRailways, ...bucket].filter(Boolean);
    const start = routePoint(a, "rail", a.x * 109 + a.y * 113);
    const goal = routePoint(b, "rail", b.x * 127 + b.y * 131);
    const softRailCost = (x, y) => {
      const i = indexOf(x, y);
      if (sea[i]) return INF;
      const rawBarrier = mountainBarrierPenalty(x, y, "rail");
      if (rawBarrier >= INF && elevation[i] > 0.88) return INF;
      const barrier = rawBarrier >= INF ? 80 + Math.max(0, elevation[i] - 0.68) * 140 + slope[i] * 28 : rawBarrier * 0.42;
      return Math.max(0.38, 1 + slope[i] * 11.5 + barrier + Math.max(0, elevation[i] - 0.60) * 5.8 + ridgeField[i] * 0.55 + (river[i] > 0.5 ? 0.95 : river[i] * 0.28) - densityValue(x, y) * 0.55 - valleyField[i] * 0.34 - plain[i] * 0.16 - coastalLowland[i] * 0.16 + normalEdgePenalty(x, y) * 0.3);
    };
    const path = aStar(start, goal, makeTransportCost(softRailCost, existingRails, railHubs, [start, goal], 4, 6.2, townAvoidNodes, 1.8, 2.2));
    const direct = pathEndpointDistance(path);
    if (path.length > 4 && pathTouchesHumanRegion(path) && direct >= 7 && pathLength(path) < 180 && pathCompactness(path) < 7.2 && pathOverlapRatio(path, existingRails, 2) < 0.76) {
      bucket.push(path);
      incrementDegree(railDegree, a);
      incrementDegree(railDegree, b);
      return true;
    }
    return false;
  }

  const transportCities = modernCities.filter((city) => (city.population || 0) >= 120000);
  const railBackboneNodes = [capital, ...transportCities, ...commercialPorts.filter((p) => p.portClass !== "fishing")];
  const mainRailLinks = buildHierarchicalLinks(railBackboneNodes, {
    mode: "rail",
    maxLinks: 3 + Math.floor(rand(seed, 1070) * 3),
    extraLinks: 1,
    minDistance: 16,
    maxDistance: 72,
    maxDegree: 3,
    seedOffset: 1070,
  });
  for (const link of mainRailLinks) {
    if (addRailRoute(link.a, link.b, railways) || addRailRouteRelaxed(link.a, link.b, railways) || addRailRouteForced(link.a, link.b, railways)) {
      addUniqueNode(railCore, link.a);
      addUniqueNode(railCore, link.b);
    }
  }

  const branchRailNodes = [capital, ...railCore, ...modernCities.filter((city) => city !== capital && (city.population || 0) < 220000), ...majorPorts];
  const branchRailLinks = buildHierarchicalLinks(branchRailNodes, {
    mode: "rail",
    maxLinks: 5 + Math.floor(rand(seed, 1071) * 4),
    extraLinks: 1,
    minDistance: 12,
    maxDistance: 54,
    maxDegree: 2,
    seedOffset: 1071,
  });
  for (const link of branchRailLinks) {
    if (railways.length && (addRailRoute(link.a, link.b, branchRailways) || addRailRouteRelaxed(link.a, link.b, branchRailways) || addRailRouteForced(link.a, link.b, branchRailways))) {
      addUniqueNode(railCore, link.a);
      addUniqueNode(railCore, link.b);
    }
  }

  const regionalMainRailLinks = buildRegionalTransportLinks(
    [...modernCities.filter((city) => (city.population || 0) >= 105000), ...commercialPorts.filter((p) => p.portClass !== "fishing")],
    (regionId, count) => ({
      mode: "rail",
      maxLinks: Math.min(regionId === 0 ? 5 : 3, Math.max(1, count - 1)),
      extraLinks: regionId === 0 ? 1 : 0,
      minDistance: 14,
      maxDistance: 66,
      maxDegree: 3,
      seedOffset: 2070,
    }),
    2070,
  );
  for (const link of regionalMainRailLinks) {
    if (addRailRoute(link.a, link.b, railways) || addRailRouteRelaxed(link.a, link.b, railways) || addRailRouteForced(link.a, link.b, railways)) {
      addUniqueNode(railCore, link.a);
      addUniqueNode(railCore, link.b);
    }
  }

  const regionalBranchRailLinks = buildRegionalTransportLinks(
    [...modernCities.filter((city) => (city.population || 0) >= 45000), ...commercialPorts, ...markets],
    (regionId, count) => ({
      mode: "rail",
      maxLinks: Math.min(regionId === 0 ? 6 : 3, Math.max(1, count - 1)),
      extraLinks: 0,
      minDistance: 10,
      maxDistance: 48,
      maxDegree: 2,
      seedOffset: 2071,
    }),
    2071,
  );
  for (const link of regionalBranchRailLinks) {
    if (railways.length && (addRailRoute(link.a, link.b, branchRailways) || addRailRouteRelaxed(link.a, link.b, branchRailways) || addRailRouteForced(link.a, link.b, branchRailways))) {
      addUniqueNode(railCore, link.a);
      addUniqueNode(railCore, link.b);
    }
  }

  compactPathArray(railways, { minLength: 10, maxOverlap: 0.40, maxCount: Math.max(5, generatedRegionIdsForTransport.length * 3) });
  compactPathArray(branchRailways, { minLength: 8, maxOverlap: 0.30, maxCount: Math.max(9, generatedRegionIdsForTransport.length * 4) });

  if (railways.length === 0) {
    for (const regionId of generatedRegionIdsForTransport) {
      const primary = primaryNodeForRegion(regionId);
      const targets = regionNodes(modernCities, regionId)
        .filter((city) => primary && nodeKey(city) !== nodeKey(primary) && (city.population || 0) >= 60000)
        .sort((a, b) => (b.population || 0) - (a.population || 0))
        .slice(0, 2);
      for (const target of targets) {
        if (addRailRouteForced(primary, target, railways)) break;
      }
      if (railways.length >= Math.max(1, Math.min(3, generatedRegionIdsForTransport.length))) break;
    }
  }

  const railInfluence = influenceFromPaths([...railways, ...branchRailways], 5);
  const stationCandidates = [
    ...modernCities.map((p, i) => ({ ...routePoint(p, "rail", 1900 + i), score: p.score + 0.46, kind: "Major Station", population: p.population })),
    ...railways.flatMap((path) => samplePath(path, 18 + Math.floor(rand(seed, path.length) * 14))).map((p) => ({ ...p, kind: "Station", score: 0.52 + agriculture[indexOf(p.x, p.y)] * 0.2 })),
    ...branchRailways.flatMap((path) => samplePath(path, 16 + Math.floor(rand(seed, path.length + 99) * 14))).map((p) => ({ ...p, kind: "Station", score: 0.42 + agriculture[indexOf(p.x, p.y)] * 0.2 })),
  ];

  let stations = pickEntities(stationCandidates, { max: 14 + Math.floor(rand(seed, 1080) * 22), minDistance: 6, threshold: 0.38, seed: seed + 1080 });

  const industrialScore = new Float32Array(SIZE);
  for (let y = 4; y < MAP_H - 4; y++) {
    for (let x = 4; x < MAP_W - 4; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const nearPort = 1 / (1 + distanceToNearest(majorPorts.length ? majorPorts : commercialPorts, x, y) / 5);
      const nearCity = distanceToNearest(modernCities, x, y);
      const cityEdge = nearCity > 5 && nearCity < 20 ? 0.22 : nearCity <= 5 ? -0.25 : 0;
      industrialScore[i] = clamp(plain[i] * 0.24 + coastalLowland[i] * 0.24 + railInfluence[i] * 0.38 + nearPort * 0.58 + river[i] * 0.04 + cityEdge - slope[i] * 0.36 - ridgeField[i] * 0.18 - floodplain[i] * 0.03);
    }
  }

  let industrialZones = pickPoints(industrialScore, {
    threshold: 0.31 + rand(seed, 1091) * 0.09,
    max: 4 + Math.floor(rand(seed, 1092) * 13),
    minDistance: 10,
    seedOffset: 1090,
    predicate: (x, y, i) => !sea[i],
  }).map((p) => ({ ...p, kind: "Industrial Zone" }));

  function roadCost(x, y) {
    const i = indexOf(x, y);
    if (sea[i]) return INF;
    const barrier = mountainBarrierPenalty(x, y, "road");
    if (barrier >= INF) return INF;
    const density = densityValue(x, y);
    const nodeAvoid = distanceToNearest(townAvoidNodes, x, y) < 2.2 ? 2.0 : 0;
    return Math.max(0.35, 1 + slope[i] * 17.8 + barrier + Math.max(0, elevation[i] - 0.54) * 9.2 + nodeAvoid + (river[i] > 0.45 ? 0.85 : 0) + floodplain[i] * 0.22 - density * 0.50 - plain[i] * 0.22 - valleyField[i] * 0.28 - coastalLowland[i] * 0.20 + ridgeField[i] * 1.15 + normalEdgePenalty(x, y) + hash2(x, y, seed + 333) * 0.08);
  }

  function expresswayCost(x, y) {
    const i = indexOf(x, y);
    if (sea[i]) return INF;
    const barrier = mountainBarrierPenalty(x, y, "express");
    if (barrier >= INF) return INF;
    const density = densityValue(x, y);
    const midDensity = midDensityAffinity(x, y);
    const cityDistance = distanceToNearest(modernCities, x, y);
    const coreAvoid = cityDistance < 2.2 ? 18.0 : cityDistance < 4.5 ? 7.0 : cityDistance < 7.5 ? 2.0 : 0;
    const marketAvoid = distanceToNearest(markets, x, y) < 2.5 ? 1.8 : 0;
    const lowDensityPenalty = density < 0.10 ? (0.10 - density) * 4.8 : 0;
    const urbanCorridorBonus = density * 1.18 + midDensity * 0.28 + (cityDistance >= 4 && cityDistance <= 16 ? 0.40 : 0);
    const constructionCost = 0.58 + slope[i] * 28.0 + barrier * 1.10 + Math.max(0, elevation[i] - 0.60) * 16.0 + ridgeField[i] * 1.45 + (river[i] > 0.45 ? 1.15 : river[i] * 0.45);
    return Math.max(0.50, 1.18 + constructionCost + coreAvoid + marketAvoid + lowDensityPenalty - urbanCorridorBonus - plain[i] * 0.12 - valleyField[i] * 0.12 - coastalLowland[i] * 0.12 + normalEdgePenalty(x, y) + hash2(x, y, seed + 444) * 0.015);
  }

  const nationalRoads = [];
  const roadDegree = new Map();
  function transportDemand(p) {
    const pop = Math.sqrt(Math.max(0, p.population || 0)) / 700;
    const capitalBoost = p.isPrefecturalCapital || p.rank === "Prefectural Capital" ? 2.1 : 0;
    const portBoost = p.portClass === "major" ? 1.4 : p.portClass === "regional" ? 0.8 : p.portClass ? 0.35 : 0;
    const historyBoost = p.kind?.includes("Castle") ? 0.55 : p.kind === "Market Town" ? 0.42 : 0;
    const gatewayBoost = p.kind === "External Gateway" ? 1.1 : 0;
    return 0.35 + pop + capitalBoost + portBoost + historyBoost + gatewayBoost;
  }

  function sameCorridorAffinity(a, b) {
    const ai = indexOf(a.x, a.y);
    const bi = indexOf(b.x, b.y);
    return Math.min(0.6, (basinField[ai] + basinField[bi]) * 0.14 + (valleyField[ai] + valleyField[bi]) * 0.10 + (coastalLowland[ai] + coastalLowland[bi]) * 0.10);
  }

  const roadTargetCandidates = [
    ...modernCities.filter((p) => inHumanRegion(p) && (p.population || 0) >= 110000),
    ...ports.filter(inHumanRegion),
    ...markets.filter(inHumanRegion),
    ...castles.filter(inHumanRegion),
  ].map((p) => ({ ...p, demand: transportDemand(p), score: (p.score || 0.4) + transportDemand(p) * 0.34 + ((p.population || 0) >= 220000 ? 0.30 : 0.05) + densityValue(p.x, p.y) * 0.20 }));
  const pickedRoadTargets = pickEntities(roadTargetCandidates, {
    max: 8 + Math.floor(rand(seed, 1101) * 9),
    minDistance: 9,
    threshold: 0.1,
    seed: seed + 1100,
  });
  const roadTargets = [
    capital,
    ...pickedRoadTargets
      .filter((p) => Math.hypot(p.x - capital.x, p.y - capital.y) > 2)
      .sort((a, b) => transportDemand(b) - transportDemand(a)),
  ];
  const roadHubs = [...modernCities, ...ports, ...markets, ...stations];
  const roadCore = [capital];

  function addNationalRoad(a, b) {
    const existing = [...nationalRoads, ...railways, ...branchRailways];
    const { path } = routeThroughTransportCorridor(a, b, "road", roadCost, existing, roadHubs, townAvoidNodes, {
      salt: a.x * 31 + a.y * 37 + b.x * 13 + b.y * 17,
      maxWaypoints: 2,
      corridorRadius: 3,
      corridorStrength: 7.0,
      avoidRadius: 2.8,
      avoidStrength: 4.4,
    });
    const direct = pathEndpointDistance(path);
    const urbanPasses = modernCities.filter((city) => path.some(([x, y]) => Math.hypot(x - city.x, y - city.y) <= Math.max(6, Math.min(13, (city.urbanRadius || 8) * 0.78)))).length;
    const densityPurpose = averagePathField(path, populationDensity) * 1.28 + averagePathField(path, plain) * 0.18 + averagePathField(path, valleyField) * 0.16;
    const passBonusOk = urbanPasses >= 1 || densityPurpose > 0.16 || direct >= 20;
    if (path.length > 3 && direct >= 12 && pathLength(path) >= 14 && pathCompactness(path) < 4.15 && pathOverlapRatio(path, existing, 2) < 0.74 && passBonusOk) {
      nationalRoads.push(path);
      incrementDegree(roadDegree, a);
      incrementDegree(roadDegree, b);
      return true;
    }
    return false;
  }

  function addNationalRoadRelaxed(a, b, bucket = nationalRoads) {
    if (!a || !b || Math.hypot(a.x - b.x, a.y - b.y) < 5) return false;
    const existing = [...nationalRoads, ...bucket, ...railways, ...branchRailways];
    const { path } = routeThroughTransportCorridor(a, b, "road", roadCost, existing, roadHubs, townAvoidNodes, {
      salt: a.x * 47 + a.y * 53 + b.x * 59 + b.y * 61,
      maxWaypoints: 2,
      corridorRadius: 4,
      corridorStrength: 8.2,
      avoidRadius: 2.3,
      avoidStrength: 3.2,
    });
    const direct = pathEndpointDistance(path);
    const densityPurpose = averagePathField(path, populationDensity) * 1.06 + averagePathField(path, plain) * 0.14 + averagePathField(path, valleyField) * 0.18 + averagePathField(path, coastalLowland) * 0.10;
    const usefulInside = pathTouchesHumanRegion(path);
    if (path.length > 3 && usefulInside && direct >= 8 && pathLength(path) >= 10 && pathCompactness(path) < 5.05 && pathOverlapRatio(path, existing, 2) < 0.86 && densityPurpose > 0.070) {
      bucket.push(path);
      incrementDegree(roadDegree, a);
      incrementDegree(roadDegree, b);
      return true;
    }
    return false;
  }

  function coverageRoadCost(x, y) {
    const i = indexOf(x, y);
    if (sea[i]) return INF;
    const rawBarrier = mountainBarrierPenalty(x, y, "road");
    if (rawBarrier >= INF && elevation[i] > 0.84) return INF;
    const barrier = rawBarrier >= INF ? 90 + Math.max(0, elevation[i] - 0.66) * 160 + slope[i] * 34 : rawBarrier * 0.38;
    const density = densityValue(x, y);
    return Math.max(0.34, 1 + slope[i] * 9.8 + barrier + Math.max(0, elevation[i] - 0.58) * 4.8 + ridgeField[i] * 0.45 + (river[i] > 0.5 ? 0.75 : 0) - density * 0.86 - plain[i] * 0.20 - valleyField[i] * 0.30 - coastalLowland[i] * 0.16 + normalEdgePenalty(x, y) * 0.4 + hash2(x, y, seed + 338) * 0.04);
  }

  function nationalRoadCorridorScore(path) {
    if (!path || !path.length) return -1;
    const density = averagePathField(path, populationDensity);
    const lowland = averagePathField(path, plain);
    const valley = averagePathField(path, valleyField);
    const coast = averagePathField(path, coastalLowland);
    const avgSlope = averagePathField(path, slope);
    const avgElevation = averagePathField(path, elevation);
    // Low-density valley/coastal corridors are allowed.  The score is meant to
    // reject truly roadless mountain/ridge alignments, not rural national roads.
    return density * 0.46 + lowland * 0.24 + valley * 0.24 + coast * 0.16 - avgSlope * 0.18 - Math.max(0, avgElevation - 0.60) * 0.15;
  }

  function isBackcountryNationalRoad(path, a, b) {
    const corridor = nationalRoadCorridorScore(path);
    const endpointWeight = nationalRoadPopulationWeight(a) + nationalRoadPopulationWeight(b);
    const direct = pathEndpointDistance(path);
    const valley = averagePathField(path, valleyField);
    const coast = averagePathField(path, coastalLowland);
    const lowland = averagePathField(path, plain);
    const naturalCorridor = valley * 0.8 + coast * 0.65 + lowland * 0.55;
    const endpointDensity = Math.max(densityValue(a.x, a.y), densityValue(b.x, b.y));
    const remoteEndpoint = endpointDensity < 0.10 && endpointWeight < 52000;
    const longRemote = direct > 34 && corridor < 0.075 && naturalCorridor < 0.16;
    return (corridor < 0.045 && naturalCorridor < 0.13) || (remoteEndpoint && corridor < 0.070 && naturalCorridor < 0.18) || longRemote;
  }

  function addNationalRoadCoverageFallback(a, b, bucket = nationalRoads) {
    if (!a || !b || Math.hypot(a.x - b.x, a.y - b.y) < 5) return false;
    const existing = [...nationalRoads, ...bucket, ...railways, ...branchRailways].filter(Boolean);
    const path = aStar(a, b, makeTransportCost(coverageRoadCost, existing, roadHubs, [a, b], 4, 8.2, townAvoidNodes, 2.0, 3.0));
    const direct = pathEndpointDistance(path);
    const usefulInside = pathTouchesHumanRegion(path);
    if (path.length > 4 && usefulInside && direct >= 7 && pathLength(path) < 150 && pathCompactness(path) < 6.9 && !isBackcountryNationalRoad(path, a, b)) {
      bucket.push(path);
      incrementDegree(roadDegree, a);
      incrementDegree(roadDegree, b);
      return true;
    }
    return false;
  }

  const roadLinks = buildHierarchicalLinks(roadTargets, {
    mode: "road",
    maxLinks: 7 + Math.floor(rand(seed, 1102) * 5),
    extraLinks: 3,
    minDistance: 13,
    maxDistance: 62,
    maxDegree: 3,
    seedOffset: 1102,
  });
  for (const link of roadLinks) {
    if (addNationalRoad(link.a, link.b)) {
      addUniqueNode(roadCore, link.a);
      addUniqueNode(roadCore, link.b);
    }
  }

  // National roads should behave like long trunk corridors: they intentionally
  // pass near as many urbanized cells/cities as possible, unlike expressways.
  const trunkCities = modernCities
    .filter((city) => inHumanRegion(city) && (city.population || 0) >= 90000)
    .slice()
    .sort((a, b) => a.x - b.x || a.y - b.y);
  for (let i = 0; i < trunkCities.length - 1; i++) {
    const a = trunkCities[i];
    const b = trunkCities[i + 1];
    const d = a && b ? Math.hypot(a.x - b.x, a.y - b.y) : 0;
    if (a && b && sameGeneratedRegion(a, b) && d >= 13 && d <= 58 && getDegree(roadDegree, a) < 5 && getDegree(roadDegree, b) < 5) addNationalRoad(a, b);
  }
  // Add a second, sparse north-south / inland-coastal layer so towns are not
  // only chained left-to-right.  This helps yellow national roads pass through
  // multiple towns instead of ending as isolated spurs.
  const verticalTrunkCities = trunkCities.slice().sort((a, b) => a.y - b.y || a.x - b.x);
  for (let i = 0; i < verticalTrunkCities.length - 2; i += 3) {
    const a = verticalTrunkCities[i];
    const b = verticalTrunkCities[Math.min(verticalTrunkCities.length - 1, i + 2)];
    const d = a && b ? Math.hypot(a.x - b.x, a.y - b.y) : 0;
    if (a && b && sameGeneratedRegion(a, b) && d >= 20 && d <= 62 && (a.population || 0) >= 110000 && (b.population || 0) >= 110000 && getDegree(roadDegree, a) < 5 && getDegree(roadDegree, b) < 5) addNationalRoad(a, b);
  }

  function addRegionalNationalRoadBackbones() {
    let added = 0;
    const regionalRoadLinks = buildRegionalTransportLinks(
      [
        ...modernCities.filter((city) => (city.population || 0) >= 45000),
        ...ports.filter((p) => p.portClass !== "fishing"),
        ...markets,
        ...castleTowns,
      ],
      (regionId, count) => ({
        mode: "road",
        maxLinks: Math.min(regionId === 0 ? 8 : 5, Math.max(1, count - 1)),
        extraLinks: regionId === 0 ? 2 : 1,
        minDistance: 9,
        maxDistance: 56,
        maxDegree: regionId === 0 ? 4 : 3,
        seedOffset: 2102,
      }),
      2102,
    );
    for (const link of regionalRoadLinks) {
      if (!sameGeneratedRegion(link.a, link.b)) continue;
      if (addNationalRoad(link.a, link.b) || addNationalRoadRelaxed(link.a, link.b)) {
        addUniqueNode(roadCore, link.a);
        addUniqueNode(roadCore, link.b);
        added++;
      }
    }
    for (const regionId of generatedRegionIdsForTransport) {
      const localCities = regionNodes(modernCities, regionId)
        .filter((city) => (city.population || 0) >= 70000)
        .sort((a, b) => a.x - b.x || a.y - b.y);
      for (let i = 0; i < localCities.length - 1; i++) {
        const a = localCities[i];
        const b = localCities[i + 1];
        const d = Math.hypot(a.x - b.x, a.y - b.y);
        if (d >= 10 && d <= 52 && getDegree(roadDegree, a) < 6 && getDegree(roadDegree, b) < 6 && (addNationalRoad(a, b) || addNationalRoadRelaxed(a, b))) added++;
      }
    }
    return added;
  }

  const regionalNationalRoadsAdded = addRegionalNationalRoadBackbones();

  function uniqueByCell(nodes) {
    const seen = new Set();
    const out = [];
    for (const node of nodes.filter(Boolean)) {
      if (!inside(node.x, node.y) || sea[indexOf(node.x, node.y)]) continue;
      const key = nodeKey(node);
      if (seen.has(key)) continue;
      seen.add(key);
      out.push(node);
    }
    return out;
  }

  function internalNationalRoadCellCount(paths = nationalRoads) {
    const seen = new Set();
    for (const path of paths) {
      for (const [x, y] of path || []) {
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (prefectureMask[i] && !sea[i]) seen.add(`${x},${y}`);
      }
    }
    return seen.size;
  }


  function nationalRoadPopulationWeight(node) {
    if (!node) return 0;
    const pop = Math.max(0, node.population || 0);
    if (pop > 0) return pop;
    if (node.portClass === "major") return 180000;
    if (node.portClass === "regional") return 90000;
    if (node.portClass) return 35000;
    if (node.kind === "Market Town" || node.kind === "Market City") return 55000;
    if (node.kind?.includes("Castle")) return 45000;
    return 18000;
  }

  function nationalRoadPopulationCoverage(paths = nationalRoads, radius = 7.0) {
    const ruralNodes = villages
      .filter((v) => prefectureMask[indexOf(v.x, v.y)] && !sea[indexOf(v.x, v.y)] && (transportDemand(v) > 0.24 || settlementCluster[indexOf(v.x, v.y)] > 0.33))
      .sort((a, b) => transportDemand(b) - transportDemand(a))
      .slice(0, 18);
    const nodes = uniqueByCell([
      capital,
      ...modernCities.filter((city) => inFocusedPrefecture(city) && ((city.population || 0) >= 60000 || city.isPrefecturalCapital)),
      ...ports.filter((p) => inFocusedPrefecture(p) && p.portClass !== "fishing"),
      ...markets.filter(inFocusedPrefecture),
      ...castleTowns.filter(inFocusedPrefecture),
      ...ruralNodes,
    ]);
    let total = 0;
    let covered = 0;
    const uncovered = [];
    for (const node of nodes) {
      const weight = nationalRoadPopulationWeight(node);
      if (weight <= 0) continue;
      total += weight;
      const d = nearestPathCellDistance(node, paths);
      if (d <= radius) covered += weight;
      else uncovered.push({ node, weight, distance: d, score: weight * (1 + Math.min(2.8, d / 12)) + transportDemand(node) * 48000 });
    }
    uncovered.sort((a, b) => b.score - a.score);
    const uncoveredPopulation = Math.max(0, total - covered);
    return { ratio: total ? covered / total : 1, total, covered, uncoveredPopulation, uncovered };
  }

  // Metropolitan national roads are split by role: yellow radial roads connect
  // the large city to neighbouring cities/ports; white ring roads are generated
  // later as ordinary urban ring roads.
  const metroRoadHubs = uniqueByCell([capital, ...modernCities.filter((city) => city !== capital && (city.population || 0) >= 240000)])
    .filter((city) => inHumanRegion(city))
    .slice(0, Math.max(4, generatedRegionIdsForTransport.length + 2));

  function addMetroRadialNationalRoads() {
    let added = 0;
    for (const hub of metroRoadHubs) {
      const maxRadials = (hub.population || 0) >= 900000 || hub.isPrefecturalCapital ? 5 : 3;
      const bySector = new Map();
      const candidates = uniqueByCell([
        ...modernCities.filter((city) => city !== hub && (city.population || 0) >= 70000),
        ...ports.filter((p) => p.portClass !== "fishing"),
        ...markets,
      ]);
      for (const node of candidates) {
        if (!sameGeneratedRegion(hub, node)) continue;
        const d = Math.hypot(node.x - hub.x, node.y - hub.y);
        if (d < 9 || d > 46) continue;
        const angle = Math.atan2(node.y - hub.y, node.x - hub.x);
        const sector = Math.floor(((angle + Math.PI) / (Math.PI * 2)) * 8);
        const coveredPenalty = nearestPathCellDistance(node, nationalRoads) <= 6.2 ? 0.72 : 0;
        const score = pairTransportScore(hub, node, "road") + transportDemand(node) * 0.44 + densityValue(node.x, node.y) * 0.22 - getDegree(roadDegree, node) * 0.16 - coveredPenalty;
        const old = bySector.get(sector);
        if (!old || score > old.score) bySector.set(sector, { node, score, d });
      }
      const sectorTargets = [...bySector.values()].sort((a, b) => b.score - a.score);
      let made = 0;
      for (const { node } of sectorTargets) {
        if (made >= maxRadials) break;
        if (getDegree(roadDegree, hub) >= 10 || getDegree(roadDegree, node) >= 7) continue;
        if (addNationalRoad(hub, node) || addNationalRoadRelaxed(hub, node) || addNationalRoadCoverageFallback(hub, node)) {
          made++;
          added++;
        }
      }
    }
    return added;
  }

  function ensureInternalNationalRoadCoverage(maxAdded = 7) {
    let added = 0;
    const minimumCells = Math.max(72, Math.floor((MAP_W + MAP_H) * 0.58));
    const targetPopulationCoverage = 0.85;
    const maxUncoveredPopulation = 50000;
    const minimumNetworkPaths = 10;
    const hasEnoughCells = () => internalNationalRoadCellCount() >= minimumCells && nationalRoads.length >= minimumNetworkPaths;
    const coverageState = () => nationalRoadPopulationCoverage(nationalRoads, 8.5);
    const hasEnoughPopulationCoverage = () => {
      const state = coverageState();
      return state.ratio >= targetPopulationCoverage && state.uncoveredPopulation <= maxUncoveredPopulation;
    };
    if (hasEnoughPopulationCoverage() && nationalRoads.length >= minimumNetworkPaths) return added;

    const populationNodes = uniqueByCell([
      capital,
      ...modernCities.filter((city) => inFocusedPrefecture(city) && (city.population || 0) >= 85000),
      ...ports.filter((p) => inFocusedPrefecture(p) && (p.portClass === "major" || p.portClass === "regional")),
      ...markets.filter((p) => inFocusedPrefecture(p) && transportDemand(p) > 0.35),
      ...castleTowns.filter((p) => inFocusedPrefecture(p) && transportDemand(p) > 0.35),
    ]).sort((a, b) => nationalRoadPopulationWeight(b) - nationalRoadPopulationWeight(a) || transportDemand(b) - transportDemand(a));
    const internalNodes = uniqueByCell([
      capital,
      ...modernCities.filter((city) => prefectureMask[indexOf(city.x, city.y)] && (city.population || 0) >= 85000),
      ...ports.filter((p) => prefectureMask[indexOf(p.x, p.y)] && (p.portClass === "major" || p.portClass === "regional")),
      ...markets.filter((p) => prefectureMask[indexOf(p.x, p.y)] && transportDemand(p) > 0.35),
      ...castleTowns.filter((p) => prefectureMask[indexOf(p.x, p.y)] && transportDemand(p) > 0.35),
    ]).sort((a, b) => nationalRoadPopulationWeight(b) - nationalRoadPopulationWeight(a) || transportDemand(b) - transportDemand(a));
    if (populationNodes.length < 2 && internalNodes.length < 2) return added;

    function targetIsWorthNationalRoad(node) {
      if (!node || node === capital) return false;
      const w = nationalRoadPopulationWeight(node);
      if (w >= 135000) return true;
      if (node.portClass === "major" || node.portClass === "regional") return true;
      if (densityValue(node.x, node.y) >= 0.18 && w >= 75000) return true;
      if (w >= 18000 && transportDemand(node) >= 0.28 && (valleyField[indexOf(node.x, node.y)] > 0.20 || coastalLowland[indexOf(node.x, node.y)] > 0.18 || plain[indexOf(node.x, node.y)] > 0.36)) return true;
      return false;
    }

    function tryCoverageLink(a, b) {
      if (!a || !b || Math.hypot(a.x - b.x, a.y - b.y) < 6) return false;
      const before = coverageState().ratio;
      const oldCount = nationalRoads.length;
      if (!(addNationalRoadRelaxed(a, b) || addNationalRoadCoverageFallback(a, b))) return false;
      const path = nationalRoads[nationalRoads.length - 1];
      const after = coverageState().ratio;
      if (isBackcountryNationalRoad(path, a, b) && after - before < 0.025) {
        nationalRoads.splice(oldCount, nationalRoads.length - oldCount);
        return false;
      }
      return true;
    }

    const currentCoverage = coverageState();
    const uncoveredPopulationTargets = currentCoverage.uncovered.map((item) => item.node).filter(targetIsWorthNationalRoad);
    const internalSpanNodes = internalNodes.length >= 2 ? internalNodes : populationNodes;
    const byX = internalSpanNodes.slice().sort((a, b) => a.x - b.x);
    const byY = internalSpanNodes.slice().sort((a, b) => a.y - b.y);
    const edgeBackstops = uniqueByCell([byX[0], byX[byX.length - 1], byY[0], byY[byY.length - 1]])
      .filter((node) => targetIsWorthNationalRoad(node) && Math.hypot(node.x - capital.x, node.y - capital.y) >= 9 && densityValue(node.x, node.y) >= 0.12)
      .sort((a, b) => nationalRoadPopulationWeight(b) - nationalRoadPopulationWeight(a));
    const ruralBackboneTargets = villages
      .filter((v) => prefectureMask[indexOf(v.x, v.y)] && !sea[indexOf(v.x, v.y)] && nearestPathCellDistance(v, nationalRoads) > 7.0)
      .filter((v) => transportDemand(v) >= 0.28 && (valleyField[indexOf(v.x, v.y)] > 0.20 || coastalLowland[indexOf(v.x, v.y)] > 0.18 || plain[indexOf(v.x, v.y)] > 0.36))
      .sort((a, b) => transportDemand(b) - transportDemand(a))
      .slice(0, 5);
    const primaryTargets = uniqueByCell([
      ...uncoveredPopulationTargets.slice(0, 5),
      ...edgeBackstops.slice(0, 2),
      ...ruralBackboneTargets,
      ...populationNodes.filter(targetIsWorthNationalRoad).slice(0, 4),
    ]).filter((node) => node !== capital && Math.hypot(node.x - capital.x, node.y - capital.y) >= 7);

    for (const target of primaryTargets) {
      if (added >= maxAdded || hasEnoughPopulationCoverage()) break;
      if (nearestPathCellDistance(target, nationalRoads) <= 6.5) continue;
      if (tryCoverageLink(capital, target)) added++;
    }

    const orderedByPopulation = populationNodes
      .filter((node) => targetIsWorthNationalRoad(node) && nearestPathCellDistance(node, nationalRoads) > 7.0)
      .sort((a, b) => nationalRoadPopulationWeight(b) - nationalRoadPopulationWeight(a));
    for (const target of orderedByPopulation) {
      if (added >= maxAdded || hasEnoughPopulationCoverage()) break;
      const anchor = populationNodes
        .filter((node) => nodeKey(node) !== nodeKey(target) && nearestPathCellDistance(node, nationalRoads) <= 5.8)
        .sort((a, b) => Math.hypot(a.x - target.x, a.y - target.y) - Math.hypot(b.x - target.x, b.y - target.y))[0] || capital;
      const d = Math.hypot(anchor.x - target.x, anchor.y - target.y);
      if (d < 7 || d > 50) continue;
      if (getDegree(roadDegree, anchor) >= 7 || getDegree(roadDegree, target) >= 5) continue;
      if (tryCoverageLink(anchor, target)) added++;
    }

    // Cell-count backstop is deliberately weak: use it only when both network
    // shape and population coverage are poor.  This avoids forcing yellow roads
    // into sparsely inhabited mountain or peninsula tips just to hit 100% coverage.
    if (!hasEnoughCells() && coverageState().ratio < 0.80) {
      for (const chain of [byX, byY]) {
        if (added >= maxAdded) break;
        for (let i = 0; i < chain.length - 1; i += 3) {
          if (added >= maxAdded || hasEnoughCells() || hasEnoughPopulationCoverage()) break;
          const a = chain[i];
          const b = chain[i + 1];
          const d = Math.hypot(a.x - b.x, a.y - b.y);
          if (d < 10 || d > 42) continue;
          if (!targetIsWorthNationalRoad(a) && !targetIsWorthNationalRoad(b)) continue;
          if (getDegree(roadDegree, a) >= 6 || getDegree(roadDegree, b) >= 6) continue;
          if (tryCoverageLink(a, b)) added++;
        }
      }
    }
    return added;
  }

  function ensureRegionalNationalRoadCoverage(maxAddedPerRegion = 3) {
    let added = 0;
    for (const regionId of generatedRegionIdsForTransport) {
      const primary = primaryNodeForRegion(regionId);
      if (!primary) continue;
      const localNodes = uniqueByCell([
        primary,
        ...regionNodes(modernCities, regionId).filter((city) => (city.population || 0) >= 60000),
        ...regionNodes(ports, regionId).filter((p) => p.portClass !== "fishing"),
        ...regionNodes(markets, regionId),
        ...regionNodes(castleTowns, regionId),
        ...regionNodes(villages, regionId).filter((v) => transportDemand(v) >= 0.26).slice(0, 4),
      ]).sort((a, b) => nationalRoadPopulationWeight(b) - nationalRoadPopulationWeight(a) || transportDemand(b) - transportDemand(a));
      let made = 0;
      for (const target of localNodes) {
        if (made >= maxAddedPerRegion) break;
        if (nodeKey(target) === nodeKey(primary) || nearestPathCellDistance(target, nationalRoads) <= 6.5) continue;
        const anchor = localNodes
          .filter((node) => nodeKey(node) !== nodeKey(target) && nearestPathCellDistance(node, nationalRoads) <= 6.0)
          .sort((a, b) => Math.hypot(a.x - target.x, a.y - target.y) - Math.hypot(b.x - target.x, b.y - target.y))[0] || primary;
        const d = Math.hypot(anchor.x - target.x, anchor.y - target.y);
        if (d < 7 || d > 56) continue;
        if (addNationalRoadRelaxed(anchor, target) || addNationalRoadCoverageFallback(anchor, target)) {
          made++;
          added++;
        }
      }
    }
    return added;
  }

  const metroRadialNationalRoadsAdded = addMetroRadialNationalRoads();
  const internalNationalRoadFallbacks = ensureInternalNationalRoadCoverage(8);
  const regionalNationalRoadFallbacks = ensureRegionalNationalRoadCoverage(3);

  const expressways = [];
  const expressDegree = new Map();
  const expressCore = [capital];

  function snapPathToExistingExpressways(path, existingPaths, radius = 2.4) {
    if (!path?.length || !existingPaths?.length) return path || [];
    const snapped = [];
    const skipEnd = Math.min(5, Math.floor(path.length / 5));
    for (let pi = 0; pi < path.length; pi++) {
      const [x, y] = path[pi];
      let best = null;
      let bestD = radius;
      if (pi >= skipEnd && pi < path.length - skipEnd) {
        for (const existing of existingPaths) {
          for (const [ex, ey] of existing) {
            const d = Math.hypot(x - ex, y - ey);
            if (d < bestD) { bestD = d; best = [ex, ey]; }
          }
        }
      }
      const next = best || [x, y];
      const last = snapped[snapped.length - 1];
      if (!last || last[0] !== next[0] || last[1] !== next[1]) snapped.push(next);
    }
    return snapped;
  }


  function addExpressway(a, b, bucket = expressways) {
    const existing = [...expressways, ...nationalRoads, ...railways, ...branchRailways];
    let { path } = routeThroughTransportCorridor(a, b, "express", expresswayCost, existing, roadHubs, townAvoidNodes, {
      salt: a.x * 41 + a.y * 43 + b.x * 19 + b.y * 29,
      maxWaypoints: 1,
      corridorRadius: 5,
      corridorStrength: 10.8,
      avoidRadius: 8.5,
      avoidStrength: 14.0,
    });
    path = smoothPathByLineOfSight(path, (x, y) => expresswayCost(x, y, x, y) < INF && slope[indexOf(x, y)] < 0.54 && elevation[indexOf(x, y)] < 0.78, 12);
    path = snapPathToExistingExpressways(path, expressways, 2.4);
    const direct = pathEndpointDistance(path);
    const densityPurpose = averagePathField(path, populationDensity) * 0.8 + averagePathField(path, plain) * 0.16 + averagePathField(path, coastalLowland) * 0.12;
    const turnScore = pathTurnScore(path);
    const deviation = pathLateralDeviationRatio(path);
    const compact = pathCompactness(path);
    if (path.length > 8 && direct >= 26 && pathLength(path) >= 30 && compact < 2.28 && turnScore < 0.64 && deviation < 0.36 && pathOverlapRatio(path, existing, 2) < 0.34 && densityPurpose > 0.11) {
      bucket.push(path);
      incrementDegree(expressDegree, a);
      incrementDegree(expressDegree, b);
      return true;
    }
    return false;
  }

  function addExpresswayRelaxed(a, b, bucket = expressways) {
    if (!a || !b || Math.hypot(a.x - b.x, a.y - b.y) < 14) return false;
    const existing = [...expressways, ...nationalRoads, ...railways, ...branchRailways].filter(Boolean);
    let { path } = routeThroughTransportCorridor(a, b, "express", expresswayCost, existing, roadHubs, townAvoidNodes, {
      salt: a.x * 149 + a.y * 151 + b.x * 157 + b.y * 163,
      maxWaypoints: 1,
      corridorRadius: 6,
      corridorStrength: 8.4,
      avoidRadius: 7.5,
      avoidStrength: 10.0,
    });
    path = smoothPathByLineOfSight(path, (x, y) => expresswayCost(x, y, x, y) < INF && slope[indexOf(x, y)] < 0.62 && elevation[indexOf(x, y)] < 0.84, 10);
    path = snapPathToExistingExpressways(path, expressways, 2.6);
    const direct = pathEndpointDistance(path);
    const densityPurpose = averagePathField(path, populationDensity) * 0.62 + averagePathField(path, plain) * 0.18 + averagePathField(path, valleyField) * 0.12 + averagePathField(path, coastalLowland) * 0.12;
    if (path.length > 7 && pathTouchesHumanRegion(path) && direct >= 18 && pathLength(path) >= 20 && pathCompactness(path) < 3.05 && pathTurnScore(path) < 0.82 && pathLateralDeviationRatio(path) < 0.54 && pathOverlapRatio(path, existing, 2) < 0.48 && densityPurpose > 0.065) {
      bucket.push(path);
      incrementDegree(expressDegree, a);
      incrementDegree(expressDegree, b);
      return true;
    }
    return false;
  }

  const expressNodes = [capital, ...modernCities.filter((p) => inHumanRegion(p) && (p.population || 0) >= 220000), ...majorPorts.filter((p) => inHumanRegion(p) && p.portClass === "major")];
  // Expressways are intentionally light in this urban-model iteration.  Full
  // expressway routing is expensive and will be revisited with the transport
  // rewrite; for now, derive at most one express corridor from an existing trunk.
  const expressLinks = [];
  for (const link of expressLinks) {
    if (!sameGeneratedRegion(link.a, link.b)) continue;
    if (addExpressway(link.a, link.b) || addExpresswayRelaxed(link.a, link.b)) {
      addUniqueNode(expressCore, link.a);
      addUniqueNode(expressCore, link.b);
    }
  }

  // Skipped: relaxed all-region expressway fallback uses repeated A* searches.

  if (expressways.length === 0) {
    const trunkCandidate = nationalRoads
      .filter((path) => path && path.length >= 26 && pathEndpointDistance(path) >= 18)
      .map((path) => ({
        path,
        score: pathEndpointDistance(path) * 0.12 + averagePathField(path, populationDensity) * 3.0 + averagePathField(path, plain) * 0.7 + averagePathField(path, coastalLowland) * 0.45 - pathCompactness(path) * 0.25,
      }))
      .sort((a, b) => b.score - a.score)[0];
    if (trunkCandidate) expressways.push(trunkCandidate.path);
  }

  const ringRoads = [];
  const ringExpressways = [];
  const ringRailways = [];

  function ringAnchorCandidates(city, mode, targetRadius, sectors = 8) {
    const anchors = [];
    const minR = Math.max(5, targetRadius - 5);
    const maxR = targetRadius + 7;
    for (let s = 0; s < sectors; s++) {
      const angle0 = (s / sectors) * Math.PI * 2;
      let best = null;
      let bestScore = -INF;
      for (let dy = -Math.ceil(maxR); dy <= Math.ceil(maxR); dy++) {
        for (let dx = -Math.ceil(maxR); dx <= Math.ceil(maxR); dx++) {
          const d = Math.hypot(dx, dy);
          if (d < minR || d > maxR) continue;
          const angle = Math.atan2(dy, dx);
          let delta = Math.abs(Math.atan2(Math.sin(angle - angle0), Math.cos(angle - angle0)));
          if (delta > Math.PI / sectors * 0.95) continue;
          const x = city.x + dx;
          const y = city.y + dy;
          if (!inside(x, y)) continue;
          const i = indexOf(x, y);
          if (sea[i] || regionIdAt(x, y) !== regionIdAt(city.x, city.y)) continue;
          const barrier = mode === "road" ? mountainBarrierPenalty(x, y, "road") : mountainBarrierPenalty(x, y, mode === "express" ? "express" : "rail");
          if (barrier >= INF) continue;
          const density = densityValue(x, y);
          const densityTerm = mode === "rail" ? density * 0.75 : mode === "express" ? midDensityAffinity(x, y) * 0.72 : density * 0.28 + midDensityAffinity(x, y) * 0.22;
          const score = plain[i] * 0.72 + agriculture[i] * 0.12 + densityTerm - slope[i] * 1.25 - Math.max(0, elevation[i] - 0.58) * 1.3 - barrier * 0.01 - Math.abs(d - targetRadius) * 0.035 + hash2(x, y, seed + 4100 + s * 37 + mode.length * 101) * 0.08;
          if (score > bestScore) {
            bestScore = score;
            best = { x, y, score, kind: `${mode} ring anchor`, parent: city };
          }
        }
      }
      if (best) anchors.push(best);
    }
    return anchors;
  }

  function ringCost(baseCost, city, targetRadius, mode) {
    return (x, y, cx, cy) => {
      const base = baseCost(x, y, cx, cy);
      if (base >= INF) return base;
      const d = Math.hypot(x - city.x, y - city.y);
      const tooClose = Math.max(0, targetRadius * 0.46 - d);
      const tooFar = Math.max(0, d - targetRadius * 1.55);
      const bandPenalty = tooClose * 0.34 + tooFar * 0.16 + Math.abs(d - targetRadius) * 0.018;
      const density = densityValue(x, y);
      const densityBias = mode === "rail" ? -density * 0.42 : mode === "express" ? -midDensityAffinity(x, y) * 0.32 + Math.max(0, density - 0.82) * 0.8 : -density * 0.12;
      return Math.max(0.36, base + bandPenalty + densityBias);
    };
  }

  function softRingRailCost(x, y) {
    const i = indexOf(x, y);
    const barrier = mountainBarrierPenalty(x, y, "rail");
    if (sea[i] || barrier >= INF) return INF;
    const density = densityValue(x, y);
    return Math.max(0.38, 1 + slope[i] * 14 + barrier + Math.max(0, elevation[i] - 0.56) * 22 + (river[i] > 0.5 ? 1.3 : river[i] * 0.6) - density * 0.62 - plain[i] * 0.20 + normalEdgePenalty(x, y) + hash2(x, y, seed + 7222) * 0.05);
  }

  function softRingExpressCost(x, y) {
    const i = indexOf(x, y);
    const barrier = mountainBarrierPenalty(x, y, "express");
    if (sea[i] || barrier >= INF) return INF;
    return Math.max(0.38, 1 + slope[i] * 13 + barrier + Math.max(0, elevation[i] - 0.58) * 20 + (river[i] > 0.5 ? 1.0 : river[i] * 0.5) - midDensityAffinity(x, y) * 0.42 - plain[i] * 0.14 + normalEdgePenalty(x, y) + hash2(x, y, seed + 7444) * 0.05);
  }

  function addEnvironmentalRing(city, mode, bucket, baseCost, existingPaths, targetRadius) {
    const anchors = ringAnchorCandidates(city, mode, targetRadius, mode === "road" ? 7 : 8);
    if (anchors.length < 3) return 0;
    let made = 0;
    const cost = ringCost(baseCost, city, targetRadius, mode);
    for (let i = 0; i < anchors.length - (anchors.length < 4 ? 1 : 0); i++) {
      const a = anchors[i];
      const b = anchors[(i + 1) % anchors.length];
      if (Math.hypot(a.x - b.x, a.y - b.y) > targetRadius * 1.85) continue;
      const path = aStar(a, b, makeTransportCost(cost, [...existingPaths, ...bucket], roadHubs, [a, b], mode === "road" ? 3 : 4, mode === "road" ? 4.8 : 7.0, townAvoidNodes, mode === "express" ? 3.8 : 2.2, mode === "express" ? 4.8 : 2.8));
      if (path.length >= 5 && path.length <= targetRadius * 8.0) {
        bucket.push(path);
        made++;
      }
    }
    return made;
  }

  function flexibleRingAnchors(city, targetRadius, maxAnchors = 6) {
    const candidates = [];
    const maxR = targetRadius + 11;
    const minR = Math.max(5, targetRadius * 0.45);
    for (let dy = -Math.ceil(maxR); dy <= Math.ceil(maxR); dy++) {
      for (let dx = -Math.ceil(maxR); dx <= Math.ceil(maxR); dx++) {
        const d = Math.hypot(dx, dy);
        if (d < minR || d > maxR) continue;
        const x = city.x + dx;
        const y = city.y + dy;
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (sea[i] || !prefectureMask[i] || elevation[i] > 0.82) continue;
        const score = plain[i] * 0.7 + midDensityAffinity(x, y) * 0.32 + densityValue(x, y) * 0.2 - slope[i] * 1.15 - Math.max(0, elevation[i] - 0.58) * 0.88 - Math.abs(d - targetRadius) * 0.02 + hash2(x, y, seed + 7555) * 0.06;
        candidates.push({ x, y, score, angle: Math.atan2(dy, dx), kind: "flexible ring anchor", parent: city });
      }
    }
    return pickEntities(candidates, { max: maxAnchors, minDistance: 5, threshold: -1, seed: seed + city.x * 83 + city.y * 89 })
      .sort((a, b) => a.angle - b.angle);
  }

  function addLooseEnvironmentalRing(city, bucket, baseCost, targetRadius) {
    let anchors = ringAnchorCandidates(city, "road", targetRadius, 6);
    if (anchors.length < 3) anchors = flexibleRingAnchors(city, targetRadius, 6);
    if (anchors.length < 2) return 0;
    let made = 0;
    for (let i = 0; i < anchors.length; i++) {
      const a = anchors[i];
      const b = anchors[(i + 1) % anchors.length];
      const path = aStar(a, b, (x, y, cx, cy) => {
        const base = baseCost(x, y, cx, cy);
        if (base >= INF) return INF;
        const d = Math.hypot(x - city.x, y - city.y);
        const band = Math.max(0, targetRadius * 0.42 - d) * 0.22 + Math.max(0, d - targetRadius * 1.7) * 0.14 + Math.abs(d - targetRadius) * 0.012;
        return Math.max(0.3, base + band);
      });
      if (path.length >= 4 && path.length <= targetRadius * 9.0) {
        bucket.push(path);
        made++;
      }
    }
    return made;
  }

  const mediumRingCities = modernCities.filter((c) => inHumanRegion(c) && (c.population || 0) >= 160000).slice(0, Math.max(5, generatedRegionIdsForTransport.length * 2));
  let metroRingRoadSegmentsAdded = 0;
  for (const city of mediumRingCities) {
    const radius = clamp(8 + Math.sqrt(city.population || 100000) / 175, 10, 22);
    metroRingRoadSegmentsAdded += addEnvironmentalRing(city, "road", ringRoads, roadCost, [...nationalRoads, ...railways, ...branchRailways], radius);
  }
  const largeRingCities = uniqueByCell([...metroRoadHubs, ...modernCities.filter((c) => inHumanRegion(c) && (c.population || 0) >= 420000)]).slice(0, Math.max(3, generatedRegionIdsForTransport.length));
  for (const city of largeRingCities) {
    const roadRadius = clamp(10 + Math.sqrt(city.population || 400000) / 155, 13, 28);
    const railRadius = Math.max(8, roadRadius - 4);
    const roadRingSegments = addEnvironmentalRing(city, "road", ringRoads, roadCost, [...nationalRoads, ...expressways, ...railways, ...branchRailways], roadRadius);
    metroRingRoadSegmentsAdded += roadRingSegments || addLooseEnvironmentalRing(city, ringRoads, roadCost, roadRadius);
    // Expressway rings are intentionally disabled; expressways stay as sparse interurban corridors.
    const railRingSegments = addEnvironmentalRing(city, "rail", ringRailways, railCost, [...railways, ...branchRailways, ...nationalRoads, ...expressways], railRadius);
    if (railRingSegments === 0) addLooseEnvironmentalRing(city, ringRailways, softRingRailCost, railRadius);
  }
  ringExpressways.length = 0;
  compactPathArray(ringRoads, { minLength: 8, maxOverlap: 0.32, maxCount: Math.max(18, generatedRegionIdsForTransport.length * 5) });
  compactPathArray(ringRailways, { minLength: 8, maxOverlap: 0.26, maxCount: Math.max(8, generatedRegionIdsForTransport.length * 2) });

  const gatewayCandidates = [];
  for (let x = 0; x < MAP_W; x++) for (const y of [0, MAP_H - 1]) {
    const i = indexOf(x, y);
    if (!sea[i]) {
      const density = densityValue(x, y);
      gatewayCandidates.push({
        x, y,
        side: y === 0 ? "N" : "S",
        score: plain[i] * 0.9 + agriculture[i] * 0.35 + valleyField[i] * 0.42 + coastalLowland[i] * 0.28 + density * 0.55 + (1 - slope[i]) * 0.42 - Math.max(0, elevation[i] - 0.56) * 1.8 - ridgeField[i] * 0.42
      });
    }
  }
  for (let y = 0; y < MAP_H; y++) for (const x of [0, MAP_W - 1]) {
    const i = indexOf(x, y);
    if (!sea[i]) {
      const density = densityValue(x, y);
      gatewayCandidates.push({
        x, y,
        side: x === 0 ? "W" : "E",
        score: plain[i] * 0.9 + agriculture[i] * 0.35 + valleyField[i] * 0.42 + coastalLowland[i] * 0.28 + density * 0.55 + (1 - slope[i]) * 0.42 - Math.max(0, elevation[i] - 0.56) * 1.8 - ridgeField[i] * 0.42
      });
    }
  }

  const minExternalGatewayCount = Math.min(5, gatewayCandidates.length);
  const targetGatewayCount = Math.min(gatewayCandidates.length, 4 + Math.floor(rand(seed, 1201) * 3));
  let externalGateways = pickEntities(gatewayCandidates, {
    max: targetGatewayCount,
    minDistance: 20,
    threshold: 0.18,
    seed: seed + 1201,
  }).map((p) => ({ ...p, kind: "External Gateway" }));
  if (externalGateways.length < minExternalGatewayCount) {
    const fallbackGateways = gatewayCandidates
      .slice()
      .sort((a, b) => b.score - a.score);
    for (const gate of fallbackGateways) {
      if (externalGateways.some((p) => Math.hypot(p.x - gate.x, p.y - gate.y) < 18)) continue;
      externalGateways.push({ ...gate, kind: "External Gateway" });
      if (externalGateways.length >= minExternalGatewayCount) break;
    }
  }

  function externalRoadCost(goal) {
    return (x, y) => {
      const i = indexOf(x, y);
      if (sea[i]) return INF;
      const barrier = mountainBarrierPenalty(x, y, "road");
      if (barrier >= INF) return INF;
      const borderPenalty = nearMapEdge(x, y, 1) && !(Math.abs(x - goal.x) <= 2 && Math.abs(y - goal.y) <= 2) ? 7 : nearMapEdge(x, y, 3) ? 1.5 : 0;
      const density = densityValue(x, y);
      const nodeAvoid = distanceToNearest(townAvoidNodes, x, y) < 2.2 ? 2.0 : 0;
      return Math.max(0.35, 1 + slope[i] * 13 + barrier + Math.max(0, elevation[i] - 0.54) * 7.5 + nodeAvoid + (river[i] > 0.45 ? 0.9 : 0) + floodplain[i] * 0.24 - density * 0.3 - plain[i] * 0.24 + borderPenalty + hash2(x, y, seed + 333) * 0.06);
    };
  }
  function externalExpresswayCost(goal) {
    return (x, y) => {
      const i = indexOf(x, y);
      if (sea[i]) return INF;
      const barrier = mountainBarrierPenalty(x, y, "express");
      if (barrier >= INF) return INF;
      const borderPenalty = nearMapEdge(x, y, 1) && !(Math.abs(x - goal.x) <= 2 && Math.abs(y - goal.y) <= 2) ? 7 : nearMapEdge(x, y, 3) ? 1.5 : 0;
      const density = densityValue(x, y);
      const cityDistance = distanceToNearest(modernCities, x, y);
      const coreAvoid = cityDistance < 2.2 ? 16.0 : cityDistance < 4.5 ? 6.0 : cityDistance < 7.5 ? 1.8 : 0;
      const lowDensityPenalty = density < 0.10 ? (0.10 - density) * 4.2 : 0;
      const urbanCorridorBonus = density * 1.04 + midDensityAffinity(x, y) * 0.24 + (cityDistance >= 4 && cityDistance <= 16 ? 0.32 : 0);
      const constructionCost = 0.55 + slope[i] * 23.0 + barrier * 1.06 + Math.max(0, elevation[i] - 0.60) * 12.0 + ridgeField[i] * 1.20 + (river[i] > 0.45 ? 1 : river[i] * 0.38);
      return Math.max(0.50, 1.14 + constructionCost + coreAvoid + lowDensityPenalty + floodplain[i] * 0.18 + borderPenalty - urbanCorridorBonus - plain[i] * 0.10 + hash2(x, y, seed + 444) * 0.04);
    };
  }
  function externalRailCost(goal) {
    return (x, y) => {
      const i = indexOf(x, y);
      if (sea[i]) return INF;
      const barrier = mountainBarrierPenalty(x, y, "rail");
      if (barrier >= INF) return INF;
      const borderPenalty = nearMapEdge(x, y, 1) && !(Math.abs(x - goal.x) <= 2 && Math.abs(y - goal.y) <= 2) ? 9 : nearMapEdge(x, y, 3) ? 1.8 : 0;
      const density = densityValue(x, y);
      return Math.max(0.42, 1 + slope[i] * 22 + barrier + Math.max(0, elevation[i] - 0.52) * 14 + (river[i] > 0.45 ? 1.2 : 0) + borderPenalty - density * 1.0 - plain[i] * 0.28 + hash2(x, y, seed + 222) * 0.05);
    };
  }

  const externalRoads = [];
  const externalExpressways = [];
  const externalRailways = [];
  const nationalRoadBranchRoads = [];

  function selectExternalStart(pool, gate, degreeMap, maxDegree = 2) {
    const sorted = pool
      .filter(Boolean)
      .map((p) => ({ ...p, d: Math.hypot(p.x - gate.x, p.y - gate.y), degree: getDegree(degreeMap, p) }))
      .sort((a, b) => a.d + a.degree * 16 + (a.degree >= maxDegree ? 30 : 0) - (b.d + b.degree * 16 + (b.degree >= maxDegree ? 30 : 0)));
    return sorted.find((p) => p.degree < maxDegree) || sorted[0] || capital;
  }

  externalGateways.forEach((gate, idx) => {
    // Always lay a national-road class gateway link first.  Expressways are
    // additional sparse corridors; they should not replace the ordinary trunk
    // road connection to the neighbouring prefecture.
    const roadStartRaw = selectExternalStart([...roadCore, ...modernCities, ...ports, ...markets], gate, roadDegree, 4);
    const roadStart = routePoint(roadStartRaw, "road", gate.x * 53 + gate.y * 59);
    const roadExisting = [...nationalRoads, ...externalRoads, ...expressways, ...externalExpressways, ...railways, ...branchRailways];
    let roadPath = aStar(roadStart, gate, makeTransportCost(externalRoadCost(gate), roadExisting, roadHubs, [roadStart, gate], 3, 7.0, townAvoidNodes, 3.2, 5.6));
    if (roadPath.length > 6) {
      externalRoads.push(roadPath);
      incrementDegree(roadDegree, roadStartRaw);
      incrementDegree(roadDegree, gate);
      addUniqueNode(roadCore, gate);
    }

    const makeExpressLink = false;
    if (makeExpressLink) {
      const expressStartRaw = selectExternalStart([...expressCore, ...roadCore, ...modernCities, ...ports], gate, expressDegree, 3);
      const expressStart = routePoint(expressStartRaw, "express", gate.x * 71 + gate.y * 73);
      const expressExisting = [...expressways, ...externalExpressways, ...nationalRoads, ...externalRoads, ...railways, ...branchRailways];
      let expressPath = aStar(expressStart, gate, makeTransportCost(externalExpresswayCost(gate), expressExisting, roadHubs, [expressStart, gate], 4, 8.2, townAvoidNodes, 5.4, 7.8));
      expressPath = smoothPathByLineOfSight(expressPath, (x, y) => externalExpresswayCost(gate)(x, y, x, y) < INF && slope[indexOf(x, y)] < 0.56 && elevation[indexOf(x, y)] < 0.80, 12);
      expressPath = snapPathToExistingExpressways(expressPath, [...expressways, ...externalExpressways], 2.4);
      const direct = pathEndpointDistance(expressPath);
      const densityPurpose = averagePathField(expressPath, populationDensity) * 0.72 + averagePathField(expressPath, plain) * 0.14 + averagePathField(expressPath, coastalLowland) * 0.10;
      const turnScore = pathTurnScore(expressPath);
      const deviation = pathLateralDeviationRatio(expressPath);
      if (expressPath.length > 6 && direct >= 18 && pathLength(expressPath) >= 20 && pathCompactness(expressPath) < 2.40 && turnScore < 0.66 && deviation < 0.42 && densityPurpose > 0.08) {
        externalExpressways.push(expressPath);
        incrementDegree(expressDegree, expressStartRaw);
        incrementDegree(expressDegree, gate);
        expressCore.push(gate);
      }
    }

    if ((idx === 0 || rand(seed, 1220 + idx) > 0.5) && modernCities.length > 0) {
      const railStartRaw = selectExternalStart([...railCore, ...modernCities, ...ports], gate, railDegree, 2);
      const railStart = routePoint(railStartRaw, "rail", gate.x * 61 + gate.y * 67);
      const railExisting = [...railways, ...branchRailways, ...externalRailways, ...nationalRoads, ...externalRoads, ...expressways, ...externalExpressways];
      const railPath = aStar(railStart, gate, makeTransportCost(externalRailCost(gate), railExisting, railHubs, [railStart, gate], 4, 8.2, townAvoidNodes, 2.5, 4.4));
      if (railPath.length > 6) {
        externalRailways.push(railPath);
        incrementDegree(railDegree, railStartRaw);
        incrementDegree(railDegree, gate);
      }
    }
  });

  function nearestOtherTrunkCell(node, ownPath, paths, minDistance = 3.5) {
    let best = null;
    let bestD = INF;
    for (const path of paths) {
      if (!path || path === ownPath) continue;
      for (let k = 0; k < path.length; k += Math.max(1, Math.floor(path.length / 46))) {
        const [x, y] = path[k];
        const d = Math.hypot(node.x - x, node.y - y);
        if (d < bestD) {
          bestD = d;
          best = { x, y, d };
        }
      }
    }
    return best && bestD >= minDistance ? best : null;
  }

  function terminalIsConnectedToNationalRoad(node, ownPath = null, extraPaths = [], radius = 3.8) {
    if (nearMapEdge(node.x, node.y, 4) || distanceToNearest(externalGateways, node.x, node.y) <= 4.2) return true;
    const paths = [...nationalRoads, ...externalRoads, ...extraPaths];
    for (const path of paths) {
      if (!path || path === ownPath) continue;
      const step = Math.max(1, Math.floor(path.length / 64));
      for (let k = 0; k < path.length; k += step) {
        const [x, y] = path[k];
        if (Math.hypot(node.x - x, node.y - y) <= radius) return true;
      }
    }
    return false;
  }

  function terminalImportance(node) {
    let best = 0;
    for (const city of modernCities) {
      const d = Math.hypot(node.x - city.x, node.y - city.y);
      if (d > 6.5) continue;
      if (city.isPrefecturalCapital || city.rank === "Prefectural Capital") best = Math.max(best, 4);
      else if ((city.population || 0) >= 180000) best = Math.max(best, 3);
      else if ((city.population || 0) >= 90000) best = Math.max(best, 2);
      else best = Math.max(best, 1);
    }
    for (const port of ports) {
      const d = Math.hypot(node.x - port.x, node.y - port.y);
      if (d > 6.5) continue;
      if (port.portClass === "major") best = Math.max(best, 3);
      else if (port.portClass === "regional") best = Math.max(best, 2);
      else best = Math.max(best, 1);
    }
    for (const market of markets) if (Math.hypot(node.x - market.x, node.y - market.y) <= 5.5) best = Math.max(best, 1);
    for (const castle of castles) if (Math.hypot(node.x - castle.x, node.y - castle.y) <= 5.5) best = Math.max(best, 1);
    return best;
  }

  function repairNationalRoadDeadEnds() {
    const repairs = [];
    const trunkPaths = () => [...nationalRoads, ...externalRoads, ...repairs];
    let repairCount = 0;
    for (let pass = 0; pass < 3; pass++) {
      for (const path of nationalRoads) {
        if (!path || path.length < 8) continue;
        const terminals = [
          { x: path[0][0], y: path[0][1] },
          { x: path[path.length - 1][0], y: path[path.length - 1][1] },
        ];
        for (const terminal of terminals) {
          if (repairCount >= 28) break;
          if (terminalIsConnectedToNationalRoad(terminal, path, repairs, 3.8)) continue;
          const target = nearestOtherTrunkCell(terminal, path, trunkPaths(), 4.0);
          if (!target || target.d > 38) continue;
          const existing = [...nationalRoads, ...externalRoads, ...repairs];
          const repairPath = aStar(terminal, target, makeTransportCost(roadCost, existing, roadHubs, [terminal, target], 2, 7.4, townAvoidNodes, 2.8, 4.6));
          if (repairPath.length >= 4 && repairPath.length <= 68 && pathCompactness(repairPath) < 4.6 && pathOverlapRatio(repairPath, existing, 2) < 0.76) {
            repairs.push(repairPath);
            repairCount++;
          }
        }
      }
    }
    nationalRoads.push(...repairs);
    return repairs.length;
  }
  const nationalRoadDeadEndRepairs = 0;

  function demoteUnresolvedNationalRoadBranches() {
    let demoted = 0;
    for (let i = nationalRoads.length - 1; i >= 0; i--) {
      const path = nationalRoads[i];
      if (!path || path.length < 8) continue;
      const a = { x: path[0][0], y: path[0][1] };
      const b = { x: path[path.length - 1][0], y: path[path.length - 1][1] };
      const aConnected = terminalIsConnectedToNationalRoad(a, path, [], 3.8);
      const bConnected = terminalIsConnectedToNationalRoad(b, path, [], 3.8);
      const deadCount = (aConnected ? 0 : 1) + (bConnected ? 0 : 1);
      if (!deadCount) continue;
      const aImportance = terminalImportance(a);
      const bImportance = terminalImportance(b);
      const importantTrunk = Math.max(aImportance, bImportance) >= 3 || (aImportance >= 2 && bImportance >= 2 && pathEndpointDistance(path) >= 24);
      const looksLikeBranch = deadCount >= 2 || !importantTrunk || pathLength(path) < 34;
      if (!looksLikeBranch) continue;
      nationalRoads.splice(i, 1);
      nationalRoadBranchRoads.push(path);
      demoted++;
    }
    return demoted;
  }
  const nationalRoadBranchDemotions = demoteUnresolvedNationalRoadBranches();
  const postDemotionInternalNationalRoadFallbacks = ensureInternalNationalRoadCoverage(3);

  function terminalIsConnectedToRail(node, ownPath = null, extraPaths = [], radius = 3.8) {
    if (nearMapEdge(node.x, node.y, 4) || distanceToNearest(externalGateways, node.x, node.y) <= 4.2) return true;
    const paths = [...railways, ...branchRailways, ...externalRailways, ...ringRailways, ...extraPaths];
    for (const path of paths) {
      if (!path || path === ownPath) continue;
      const step = Math.max(1, Math.floor(path.length / 64));
      for (let k = 0; k < path.length; k += step) {
        const [x, y] = path[k];
        if (Math.hypot(node.x - x, node.y - y) <= radius) return true;
      }
    }
    return false;
  }

  function railTerminalImportance(node) {
    let best = 0;
    for (const city of modernCities) {
      const d = Math.hypot(node.x - city.x, node.y - city.y);
      if (d > 6.5) continue;
      if (city.isPrefecturalCapital || city.rank === "Prefectural Capital") best = Math.max(best, 4);
      else if ((city.population || 0) >= 220000) best = Math.max(best, 3);
      else if ((city.population || 0) >= 120000) best = Math.max(best, 2);
      else best = Math.max(best, 1);
    }
    for (const port of ports) {
      const d = Math.hypot(node.x - port.x, node.y - port.y);
      if (d > 6.5) continue;
      if (port.portClass === "major") best = Math.max(best, 2);
      else if (port.portClass === "regional") best = Math.max(best, 1);
    }
    for (const station of stations) if (Math.hypot(node.x - station.x, node.y - station.y) <= 5.0) best = Math.max(best, 1);
    return best;
  }

  function repairRailDeadEnds() {
    const repairs = [];
    const currentPaths = () => [...railways, ...branchRailways, ...externalRailways, ...ringRailways, ...repairs];
    for (let pass = 0; pass < 3; pass++) {
      for (const path of [...railways, ...branchRailways, ...externalRailways]) {
        if (!path || path.length < 8) continue;
        for (const terminal of [{ x: path[0][0], y: path[0][1] }, { x: path[path.length - 1][0], y: path[path.length - 1][1] }]) {
          if (terminalIsConnectedToRail(terminal, path, repairs, 3.8)) continue;
          let target = nearestOtherTrunkCell(terminal, path, currentPaths(), 4.0);
          if ((!target || target.d > 48) && railTerminalImportance(terminal) >= 2) {
            const candidates = [...modernCities, ...ports, ...stations]
              .map((p) => ({ p, d: Math.hypot(terminal.x - p.x, terminal.y - p.y) }))
              .filter(({ d }) => d >= 6 && d <= 42)
              .sort((a, b) => a.d - b.d);
            for (const { p } of candidates) {
              const access = routePoint(p, "rail", 9400 + p.x * 17 + p.y * 19);
              if (terminalIsConnectedToRail(access, path, repairs, 3.8)) { target = access; break; }
            }
          }
          if (!target || (target.d && target.d > 52)) continue;
          const existing = currentPaths();
          const repairPath = aStar(terminal, target, makeTransportCost(railCost, existing, railHubs, [terminal, target], 4, 8.4, townAvoidNodes, 2.0, 3.8));
          if (repairPath.length >= 4 && repairPath.length <= 60 && pathCompactness(repairPath) < 4.0 && pathTurnScore(repairPath) < 0.86 && pathOverlapRatio(repairPath, existing, 2) < 0.84) {
            repairs.push(repairPath);
          }
        }
      }
    }
    branchRailways.push(...repairs);
    return repairs.length;
  }
  const railDeadEndRepairs = 0;

  let throughExpresswayAdded = false;
  function throughExpresswayCost(a, b) {
    return (x, y) => {
      const i = indexOf(x, y);
      if (sea[i]) return INF;
      const barrier = mountainBarrierPenalty(x, y, "express");
      if (barrier >= INF) return INF;
      const nearEndpoint = Math.min(Math.hypot(x - a.x, y - a.y), Math.hypot(x - b.x, y - b.y)) <= 3;
      const borderPenalty = nearEndpoint ? 0 : nearMapEdge(x, y, 3) ? 2.2 : 0;
      const density = densityValue(x, y);
      const cityDistance = distanceToNearest(modernCities, x, y);
      const coreAvoid = cityDistance < 2.0 ? 12.0 : cityDistance < 4.5 ? 4.8 : cityDistance < 7.5 ? 1.4 : 0;
      const lowDensityPenalty = density < 0.08 ? (0.08 - density) * 4.0 : 0;
      const urbanCorridorBonus = density * 1.00 + midDensityAffinity(x, y) * 0.22 + (cityDistance >= 4 && cityDistance <= 18 ? 0.34 : 0);
      return Math.max(0.48, 1.16 + slope[i] * 21.0 + barrier * 1.08 + Math.max(0, elevation[i] - 0.62) * 13.0 + coreAvoid + lowDensityPenalty + (river[i] > 0.45 ? 1.0 : 0) + borderPenalty - urbanCorridorBonus - plain[i] * 0.14 - coastalLowland[i] * 0.12 + hash2(x, y, seed + 9101) * 0.03);
    };
  }

  function permissiveThroughExpresswayCost(a, b) {
    return (x, y) => {
      const i = indexOf(x, y);
      if (sea[i]) return 24 + nearMapEdge(x, y, 2) * 2 + hash2(x, y, seed + 9202) * 0.2;
      const rawBarrier = mountainBarrierPenalty(x, y, "express");
      const tunnelBarrier = rawBarrier >= INF ? 120 + Math.max(0, elevation[i] - 0.66) * 260 + slope[i] * 55 : rawBarrier;
      const nearEndpoint = Math.min(Math.hypot(x - a.x, y - a.y), Math.hypot(x - b.x, y - b.y)) <= 3;
      const borderPenalty = nearEndpoint ? 0 : nearMapEdge(x, y, 3) ? 2.0 : 0;
      const cityDistance = distanceToNearest(modernCities, x, y);
      const coreAvoid = cityDistance < 2.0 ? 9.0 : cityDistance < 4.5 ? 3.5 : 0;
      const density = densityValue(x, y);
      const urbanCorridorBonus = density * 0.85 + midDensityAffinity(x, y) * 0.18 + (cityDistance >= 4 && cityDistance <= 18 ? 0.24 : 0);
      return Math.max(0.52, 1.18 + slope[i] * 14.0 + tunnelBarrier * 0.42 + Math.max(0, elevation[i] - 0.66) * 18.0 + coreAvoid + borderPenalty - urbanCorridorBonus - plain[i] * 0.10 + hash2(x, y, seed + 9201) * 0.035);
    };
  }

  function pointToSegmentDistance(p, a, b) {
    const vx = b.x - a.x;
    const vy = b.y - a.y;
    const len2 = vx * vx + vy * vy;
    if (len2 <= 0.0001) return Math.hypot(p.x - a.x, p.y - a.y);
    const t = clamp(((p.x - a.x) * vx + (p.y - a.y) * vy) / len2, 0, 1);
    return Math.hypot(p.x - (a.x + vx * t), p.y - (a.y + vy * t));
  }

  function pathTurnScore(path) {
    if (!path || path.length < 3) return 0;
    let total = 0;
    let count = 0;
    for (let i = 1; i < path.length - 1; i++) {
      const [x0, y0] = path[i - 1];
      const [x1, y1] = path[i];
      const [x2, y2] = path[i + 1];
      const ax = x1 - x0;
      const ay = y1 - y0;
      const bx = x2 - x1;
      const by = y2 - y1;
      const al = Math.hypot(ax, ay);
      const bl = Math.hypot(bx, by);
      if (al < 0.01 || bl < 0.01) continue;
      const dot = clamp((ax * bx + ay * by) / (al * bl), -1, 1);
      total += Math.acos(dot);
      count++;
    }
    return count ? total / count : 0;
  }

  function pathLateralDeviationRatio(path) {
    if (!path || path.length < 3) return 0;
    const a = { x: path[0][0], y: path[0][1] };
    const b = { x: path[path.length - 1][0], y: path[path.length - 1][1] };
    const direct = Math.max(1, Math.hypot(b.x - a.x, b.y - a.y));
    let maxDeviation = 0;
    for (let i = 1; i < path.length - 1; i++) {
      const p = { x: path[i][0], y: path[i][1] };
      maxDeviation = Math.max(maxDeviation, pointToSegmentDistance(p, a, b));
    }
    return maxDeviation / direct;
  }

  function chooseThroughExpresswayVia(a, b) {
    const candidates = [capital, ...modernCities.filter((city) => (city.population || 0) >= 90000)];
    let best = null;
    let bestScore = -INF;
    for (const city of candidates) {
      if (!city || !inHumanRegion(city) || sea[indexOf(city.x, city.y)]) continue;
      const access = routePoint(city, "express", city.x * 73 + city.y * 79 + 9301);
      const lineD = pointToSegmentDistance(access, a, b);
      const density = densityValue(access.x, access.y);
      const popScore = Math.sqrt(Math.max(0, city.population || 0)) / 520;
      const score = density * 2.8 + popScore + (city.isPrefecturalCapital ? 0.9 : 0) - lineD / 38 - mountainBarrierPenalty(access.x, access.y, "express") * 0.004;
      if (score > bestScore) { bestScore = score; best = access; }
    }
    return best;
  }

  function addThroughExpressway() {
    if (externalGateways.length < 2) return false;
    let bestPair = null;
    let bestScore = -INF;
    for (let i = 0; i < externalGateways.length; i++) {
      for (let j = i + 1; j < externalGateways.length; j++) {
        const a = externalGateways[i];
        const b = externalGateways[j];
        const d = Math.hypot(a.x - b.x, a.y - b.y);
        const opposite = (a.side === "N" && b.side === "S") || (a.side === "S" && b.side === "N") || (a.side === "W" && b.side === "E") || (a.side === "E" && b.side === "W");
        const score = d + (opposite ? 42 : 0) - Math.abs((a.score || 0) - (b.score || 0)) * 3;
        if (score > bestScore) {
          bestScore = score;
          bestPair = [a, b];
        }
      }
    }
    if (!bestPair) return false;
    const [a, b] = bestPair;
    const existing = [...externalExpressways, ...expressways, ...nationalRoads, ...railways, ...branchRailways];
    const via = chooseThroughExpresswayVia(a, b);
    let path = [];
    let viaUsed = false;
    if (via) {
      let first = aStar(a, via, makeTransportCost(throughExpresswayCost(a, via), existing, roadHubs, [a, via], 5, 9.6, townAvoidNodes, 4.2, 6.0));
      first = smoothPathByLineOfSight(first, (x, y) => throughExpresswayCost(a, via)(x, y, x, y) < INF && elevation[indexOf(x, y)] < 0.88, 11);
      let second = aStar(via, b, makeTransportCost(throughExpresswayCost(via, b), [...existing, first], roadHubs, [via, b], 5, 9.6, townAvoidNodes, 4.2, 6.0));
      second = smoothPathByLineOfSight(second, (x, y) => throughExpresswayCost(via, b)(x, y, x, y) < INF && elevation[indexOf(x, y)] < 0.88, 11);
      if (first.length > 6 && second.length > 6) { path = [...first, ...second.slice(1)]; viaUsed = true; }
    }
    if (path.length < 12) {
      path = aStar(a, b, makeTransportCost(throughExpresswayCost(a, b), existing, roadHubs, [a, b], 5, 9.4, townAvoidNodes, 6.0, 9.0));
      path = smoothPathByLineOfSight(path, (x, y) => throughExpresswayCost(a, b)(x, y, x, y) < INF && elevation[indexOf(x, y)] < 0.88, 11);
    }
    path = snapPathToExistingExpressways(path, [...externalExpressways, ...expressways], 2.8);
    if (!viaUsed && (path.length < 12 || pathEndpointDistance(path) < Math.min(MAP_W, MAP_H) * 0.45 || pathCompactness(path) > 3.25)) {
      path = aStar(a, b, makeTransportCost(permissiveThroughExpresswayCost(a, b), existing, roadHubs, [a, b], 4, 7.2, townAvoidNodes, 4.0, 6.5));
      path = smoothPathByLineOfSight(path, (x, y) => permissiveThroughExpresswayCost(a, b)(x, y, x, y) < INF, 12);
      path = snapPathToExistingExpressways(path, [...externalExpressways, ...expressways], 2.8);
    }
    const turnScore = pathTurnScore(path);
    const deviation = pathLateralDeviationRatio(path);
    const compactness = pathCompactness(path);
    if (path.length < 12 || pathEndpointDistance(path) < Math.min(MAP_W, MAP_H) * 0.42 || compactness > (viaUsed ? 4.10 : 3.05) || turnScore > (viaUsed ? 0.72 : 0.60) || deviation > (viaUsed ? 0.52 : 0.38)) return false;
    externalExpressways.push(path);
    incrementDegree(expressDegree, a);
    incrementDegree(expressDegree, b);
    throughExpresswayAdded = true;
    return true;
  }
  // addThroughExpressway();

  function nearestPathCellDistance(node, paths) {
    let best = INF;
    for (const path of paths) {
      for (const [x, y] of path) best = Math.min(best, Math.hypot(node.x - x, node.y - y));
    }
    return best;
  }

  function pruneHighMountainTransport(paths, threshold = 0.82) {
    for (let i = paths.length - 1; i >= 0; i--) {
      if (paths[i].some(([x, y]) => elevation[indexOf(x, y)] > threshold)) paths.splice(i, 1);
    }
  }
  for (const paths of [railways, branchRailways, ringRailways, externalRailways, expressways]) pruneHighMountainTransport(paths, 0.82);

  const expressInfluence = influenceFromPaths([...expressways, ...externalExpressways], 6);
  const roadInfluence = influenceFromPaths([...nationalRoads, ...ringRoads, ...expressways, ...externalRoads, ...externalExpressways], 4);

  const icCandidates = [];
  for (const path of [...expressways, ...externalExpressways]) {
    icCandidates.push(...samplePath(path, 11 + Math.floor(rand(seed, path.length + 333) * 5)).map((p) => ({ ...p, score: 0.62 + plain[indexOf(p.x, p.y)] * 0.24 + midDensityAffinity(p.x, p.y) * 0.16, kind: "Interchange" })));
    for (const city of modernCities) {
      let best = null;
      let bestDistance = 999;
      for (const [x, y] of path) {
        const d = Math.hypot(x - city.x, y - city.y);
        if (d < bestDistance) { bestDistance = d; best = { x, y }; }
      }
      if (best && bestDistance > 4 && bestDistance < 18) icCandidates.push({ ...best, score: 0.8 + city.score * 0.1, kind: "Urban Interchange" });
    }
  }

  let interchanges = pickEntities(icCandidates, { max: 14 + Math.floor(rand(seed, 1130) * 18), minDistance: 7, threshold: 0.44, seed: seed + 1130 });

  const icAccessRoads = [];
  const nationalRoadAccessPoints = nationalRoads.flatMap((path) => samplePath(path, 8));
  for (const ic of interchanges) {
    const accessTargets = [
      ...industrialZones.map((p) => ({ ...p, score: 0.95 / (1 + Math.hypot(p.x - ic.x, p.y - ic.y) / 7) })),
      ...modernCities.map((p) => ({ ...routePoint(p, "road", 8200 + p.x * 7 + p.y), score: 0.72 / (1 + Math.hypot(p.x - ic.x, p.y - ic.y) / 10) })),
      ...nationalRoadAccessPoints.map((p) => ({ ...p, score: 0.62 / (1 + Math.hypot(p.x - ic.x, p.y - ic.y) / 6), kind: "National Road Access" })),
    ];
    const target = pickEntities(accessTargets, { max: 1, minDistance: 1, threshold: 0, seed: seed + 1134 + ic.x * 3 + ic.y })[0];
    if (!target || Math.hypot(target.x - ic.x, target.y - ic.y) > 22) continue;
    const path = aStar(ic, target, roadCost);
    if (path.length > 2 && path.length < 36) icAccessRoads.push(path);
  }

  const logisticsScore = new Float32Array(SIZE);
  for (let y = 4; y < MAP_H - 4; y++) {
    for (let x = 4; x < MAP_W - 4; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const nearIC = 1 / (1 + distanceToNearest(interchanges, x, y) / 3);
      const cityPenalty = distanceToNearest(modernCities, x, y) < 5 ? 0.28 : 0;
      logisticsScore[i] = clamp(nearIC * 0.56 + plain[i] * 0.24 + roadInfluence[i] * 0.22 + expressInfluence[i] * 0.16 - slope[i] * 0.32 - cityPenalty);
    }
  }

  let logisticsParks = pickPoints(logisticsScore, {
    threshold: 0.32 + rand(seed, 1141) * 0.1,
    max: 3 + Math.floor(rand(seed, 1142) * 13),
    minDistance: 9,
    seedOffset: 1140,
    predicate: (x, y, i) => !sea[i],
  }).map((p) => ({ ...p, kind: "Logistics Park" }));

  const cityInfluence = influenceFromPoints(modernCities, 34, (p) => p.urbanWeight || 1.2);
  const cityCoreInfluence = influenceFromPoints(urbanCenters, 11, (p) => p.parent?.coreRadius ? 1.35 + p.parent.coreRadius / 5 : 1.2);
  const stationInfluence = influenceFromPoints(stations, 10, () => 1);
  const railInfluence2 = influenceFromPaths([...railways, ...branchRailways, ...ringRailways, ...externalRailways], 6);
  const satelliteScore = new Float32Array(SIZE);
  const largeCitiesForSatellites = modernCities.filter((c) => (c.population || 0) >= 320000);
  for (let y = 4; y < MAP_H - 4; y++) {
    for (let x = 4; x < MAP_W - 4; x++) {
      const i = indexOf(x, y);
      if (sea[i] || !isHumanRegionCell(i)) continue;
      let ringPull = 0;
      let parent = null;
      for (const city of largeCitiesForSatellites) {
        const d = Math.hypot(city.x - x, city.y - y);
        const ideal = clamp(11 + Math.sqrt(city.population || 320000) / 150, 13, 27);
        const v = clamp(1 - Math.abs(d - ideal) / 9);
        if (v > ringPull) { ringPull = v; parent = city; }
      }
      if (!parent) continue;
      const railPull = Math.max(railInfluence2[i], stationInfluence[i] * 0.84);
      const separated = distanceToNearest(modernCities, x, y) > 7 ? 1 : 0;
      satelliteScore[i] = clamp(ringPull * 0.42 + railPull * 0.38 + populationDensity[i] * 0.14 + plain[i] * 0.2 + basinField[i] * 0.08 + agriculture[i] * 0.05 - slope[i] * 0.86 - ridgeField[i] * 0.34 - Math.max(0, elevation[i] - 0.56) * 0.72 + separated * 0.1 + hash2(x, y, seed + 1160) * 0.035);
    }
  }
  let satelliteCities = pickPoints(satelliteScore, {
    threshold: 0.43 + rand(seed, 1161) * 0.07,
    max: Math.min(14, 2 + largeCitiesForSatellites.length * 4 + Math.floor(rand(seed, 1162) * 4)),
    minDistance: 8,
    seedOffset: 1160,
    predicate: (x, y, i) => isHumanRegionCell(i),
  }).map((p, n) => {
    const parent = largeCitiesForSatellites.slice().sort((a, b) => Math.hypot(a.x - p.x, a.y - p.y) - Math.hypot(b.x - p.x, b.y - p.y))[0];
    const basePop = parent ? parent.population * (0.045 + rand(seed, 1165 + n) * 0.11) : 42000 + rand(seed, 1165 + n) * 90000;
    return { ...p, kind: "Satellite City", parentCityIndex: parent ? modernCities.indexOf(parent) : -1, population: Math.round(basePop / 1000) * 1000, urbanRadius: 5 + Math.sqrt(basePop) / 135, coreRadius: 1.5 + Math.sqrt(basePop) / 420, urbanWeight: 0.55 + Math.sqrt(basePop) / 720 };
  });
  const satelliteInfluence = influenceFromPoints(satelliteCities, 16, (p) => p.urbanWeight || 0.8);
  const oldCoreInfluence = influenceFromPoints([...castleTowns, ...markets, ...ports], 12, () => 1);
  const industrialInfluence = influenceFromPoints(industrialZones, 9, () => 1);
  const logisticsInfluence = influenceFromPoints(logisticsParks, 9, () => 1);
  const interchangeInfluence = influenceFromPoints(interchanges, 8, () => 1);
  const premodernInfluence = influenceFromPaths(premodernRoads, 4);
  const villageInfluence = influenceFromPoints(villages, 7, () => 1);

  const newTownScore = new Float32Array(SIZE);
  for (let y = 4; y < MAP_H - 4; y++) {
    for (let x = 4; x < MAP_W - 4; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const dCity = distanceToNearest(modernCities, x, y);
      const ring = dCity > 8 && dCity < 22 ? 1 : 0;
      const uplandTerrace = elevation[i] > 0.36 && elevation[i] < 0.58 && slope[i] < 0.34 && ridgeField[i] < 0.34 ? 0.24 : 0;
      newTownScore[i] = clamp(ring * 0.34 + stationInfluence[i] * 0.24 + roadInfluence[i] * 0.08 + railInfluence2[i] * 0.1 + plain[i] * 0.14 + uplandTerrace + agriculture[i] * 0.06 - slope[i] * 0.72 - ridgeField[i] * 0.22 - floodplain[i] * 0.22 - satelliteInfluence[i] * 0.18);
    }
  }

  let newTowns = pickPoints(newTownScore, {
    threshold: 0.32 + rand(seed, 1151) * 0.1,
    max: 2 + Math.floor(rand(seed, 1152) * 10),
    minDistance: 11,
    seedOffset: 1150,
    predicate: (x, y, i) => !sea[i],
  }).map((p) => ({ ...p, kind: "New Town" }));

  const minorRoads = [...nationalRoadBranchRoads];
  const trunkNodes = [...markets, ...modernCities, ...stations.slice(0, 24), ...crossings.slice(0, 16)];
  const roadNetInfluence = influenceFromPaths([...nationalRoads, ...ringRoads, ...expressways, ...ringExpressways, ...externalRoads, ...externalExpressways, ...premodernRoads], 3);

  function minorRoadCost(x, y) {
    const i = indexOf(x, y);
    if (sea[i] || elevation[i] > 0.72) return INF;
    const barrier = mountainBarrierPenalty(x, y, "minor");
    if (barrier >= INF) return INF;
    return Math.max(0.3, 1 + slope[i] * 8.4 + barrier * 0.55 + Math.max(0, elevation[i] - 0.58) * 4.4 + floodplain[i] * 0.18 + (river[i] > 0.5 ? 1.0 : 0.18 * river[i]) - plain[i] * 0.24 - valleyField[i] * 0.36 - coastalLowland[i] * 0.12 + ridgeField[i] * 0.58 - roadNetInfluence[i] * 0.35 + normalEdgePenalty(x, y) + hash2(x, y, seed + 555) * 0.15);
  }

  const connectedPairs = new Set();
  function addMinorRoad(a, b) {
    if (!a || !b || Math.hypot(a.x - b.x, a.y - b.y) > 38) return;
    const key = `${a.x},${a.y}|${b.x},${b.y}`;
    const reverseKey = `${b.x},${b.y}|${a.x},${a.y}`;
    if (connectedPairs.has(key) || connectedPairs.has(reverseKey)) return;
    connectedPairs.add(key);
    const path = aStar(a, b, minorRoadCost);
    if (path.length > 2 && path.length < 90) minorRoads.push(path);
  }

  function nearestRoadAccessNode(node, paths, sampleStep = 7) {
    let best = null;
    let bestD = INF;
    for (const path of paths) {
      if (!path || path.length === 0) continue;
      const step = Math.max(1, Math.floor(path.length / Math.max(8, Math.ceil(path.length / sampleStep))));
      for (let k = 0; k < path.length; k += step) {
        const [x, y] = path[k];
        const d = Math.hypot(node.x - x, node.y - y);
        if (d < bestD) { bestD = d; best = { x, y, kind: "Road access", d }; }
      }
    }
    return best;
  }

  // Branches from the yellow national-road network are drawn as ordinary white
  // roads.  This keeps the national-road layer as a through-network while still
  // connecting local towns, ports, castle towns, and suburban/new-town nodes.
  const nationalAccessPaths = [...nationalRoads, ...externalRoads, ...ringRoads, ...premodernRoads];
  const localTownNodes = [...modernCities, ...ports, ...markets, ...castleTowns, ...satelliteCities, ...newTowns]
    .filter((node, idx, arr) => idx === arr.findIndex((p) => Math.hypot(p.x - node.x, p.y - node.y) < 2.5));
  for (const town of localTownNodes) {
    const nearestTrunk = nearestRoadAccessNode(town, nationalAccessPaths, 6);
    if (nearestTrunk && nearestTrunk.d > 2.8 && nearestTrunk.d < 42) addMinorRoad(town, nearestTrunk);
  }
  const neighborTownNodes = localTownNodes.slice().sort((a, b) => (b.population || 0) - (a.population || 0));
  for (const town of neighborTownNodes.slice(0, 28)) {
    const neighbor = pickEntities(neighborTownNodes.filter((p) => p !== town).map((p) => ({ ...p, score: 1 / (1 + Math.hypot(p.x - town.x, p.y - town.y)) })), { max: 1, minDistance: 1, threshold: 0 })[0];
    if (neighbor && Math.hypot(neighbor.x - town.x, neighbor.y - town.y) < 22) addMinorRoad(town, neighbor);
  }

  for (const village of villages.slice().sort((a, b) => (b.population || 0) - (a.population || 0)).slice(0, 45)) {
    if (rand(seed, village.x * 13 + village.y * 17) < 0.78) {
      const target = pickEntities(trunkNodes.map((p) => ({ ...p, score: 1 / (1 + Math.hypot(p.x - village.x, p.y - village.y)) })), { max: 1, minDistance: 1, threshold: 0 })[0];
      if (target && Math.hypot(target.x - village.x, target.y - village.y) < 34) addMinorRoad(village, target);
    }
  }
  for (const market of markets.slice().sort((a, b) => (b.population || 0) - (a.population || 0)).slice(0, 32)) {
    const localVillages = pickEntities(villages.map((v) => ({ ...v, score: 1 / (1 + Math.hypot(v.x - market.x, v.y - market.y)) })), { max: 2, minDistance: 1, threshold: 0 });
    for (const v of localVillages) if (Math.hypot(v.x - market.x, v.y - market.y) < 24) addMinorRoad(market, v);
  }
  for (const pass of passes.slice(0, 8)) {
    const target = pickEntities([...markets, ...villages].map((p) => ({ ...p, score: 1 / (1 + Math.hypot(p.x - pass.x, p.y - pass.y)) })), { max: 1, minDistance: 1, threshold: 0 })[0];
    if (target) addMinorRoad(pass, target);
  }
  for (const port of ports) {
    const target = pickEntities([...markets, ...villages, ...stations.slice(0, 18)].map((p) => ({ ...p, score: 1 / (1 + Math.hypot(p.x - port.x, p.y - port.y)) })), { max: 1, minDistance: 1, threshold: 0 })[0];
    if (target && Math.hypot(target.x - port.x, target.y - port.y) < 24) addMinorRoad(port, target);
  }
  for (const localCenter of [...satelliteCities, ...newTowns]) {
    const target = pickEntities([...stations, ...markets, ...modernCities].map((p) => ({ ...p, score: 1 / (1 + Math.hypot(p.x - localCenter.x, p.y - localCenter.y)) })), { max: 1, minDistance: 1, threshold: 0 })[0];
    if (target && Math.hypot(target.x - localCenter.x, target.y - localCenter.y) < 30) addMinorRoad(localCenter, target);
  }
  for (const station of stations.slice(0, 12)) {
    const locals = pickEntities([...villages, ...markets, ...ports].map((p) => ({ ...p, score: 1 / (1 + Math.hypot(p.x - station.x, p.y - station.y)) })), { max: 1, minDistance: 1, threshold: 0 });
    for (const local of locals) if (Math.hypot(local.x - station.x, local.y - station.y) < 22) addMinorRoad(station, local);
  }
  for (const village of villages.slice(0, 16)) {
    const neighbor = pickEntities(villages.filter((v) => v !== village).map((v) => ({ ...v, score: 1 / (1 + Math.hypot(v.x - village.x, v.y - village.y)) })), { max: 1, minDistance: 1, threshold: 0 })[0];
    if (neighbor && Math.hypot(neighbor.x - village.x, neighbor.y - village.y) < 14) addMinorRoad(village, neighbor);
  }

  const combinedModernTransport = [...nationalRoads, ...externalRoads, ...externalExpressways, ...railways, ...branchRailways, ...externalRailways, ...expressways];
  const requiredTransportNodes = [capital, ...externalGateways, ...modernCities.filter((city) => (city.population || 0) >= 120000 || city.isPrefecturalCapital)];
  const connectedRequiredNodeCount = requiredTransportNodes.filter((node) => nearestPathCellDistance(node, combinedModernTransport) <= 7).length;
  const allExpresswayPaths = [...expressways, ...externalExpressways];
  const expresswayCells = allExpresswayPaths.flat();
  const expresswayAverageDensity = expresswayCells.length ? expresswayCells.reduce((sum, [x, y]) => sum + densityValue(x, y), 0) / expresswayCells.length : 0;
  const nationalRoadPopulationCoverageDebug = nationalRoadPopulationCoverage(nationalRoads, 8.5);
  const transportDebug = {
    urbanHierarchy: urbanHierarchyDebug,
    requiredNodeCount: requiredTransportNodes.length,
    connectedRequiredNodeCount,
    throughExpresswayAdded,
    nationalRoadDeadEndRepairs,
    railDeadEndRepairs,
    nationalRoadBranchDemotions,
    regionalNationalRoadsAdded,
    regionalNationalRoadFallbacks,
    metroRadialNationalRoadsAdded,
    metroRingRoadSegmentsAdded,
    internalNationalRoadFallbacks,
    postDemotionInternalNationalRoadFallbacks,
    nationalRoadPopulationCoverage: Number(nationalRoadPopulationCoverageDebug.ratio.toFixed(3)),
    nationalRoadUncoveredPopulation: Math.round(nationalRoadPopulationCoverageDebug.uncoveredPopulation || 0),
    internalNationalRoadCellCount: internalNationalRoadCellCount(),
    externalGatewayCount: externalGateways.length,
    externalNationalRoadCount: externalRoads.length,
    expresswayAverageDensity: Number(expresswayAverageDensity.toFixed(3)),
    expresswayPathCount: allExpresswayPaths.length,
    minorRoadCount: minorRoads.length,
    minorRoadTotalLength: Math.round(minorRoads.reduce((sum, path) => sum + pathLength(path), 0)),
    settlementModel: {
      villageCount: villages.length,
      marketCount: markets.length,
      averageUrbanFootprintCells: Number((modernCities.reduce((sum, city) => sum + (city.urbanFootprintCells || 0), 0) / Math.max(1, modernCities.length)).toFixed(1)),
      maxUrbanFootprintCells: Math.max(0, ...modernCities.map((city) => city.urbanFootprintCells || 0)),
    },
  };

  const newTownInfluence = influenceFromPoints(newTowns, 8, () => 1);
  const landuse = new Uint8Array(SIZE);
  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const mountain = elevation[i] > 0.62 || slope[i] > 0.46 || ridgeField[i] > 0.64;
      const farm = agriculture[i] > 0.26 && (plain[i] > 0.2 || valleyField[i] > 0.32 || basinField[i] > 0.25);
      let nearestCity = null;
      let nearestCityDistance = INF;
      for (const city of modernCities) {
        const d = Math.hypot(city.x - x, city.y - y);
        if (d < nearestCityDistance) { nearestCityDistance = d; nearestCity = city; }
      }
      const dCity = nearestCityDistance;
      const populationScale = nearestCity ? clamp(Math.log10(Math.max(10000, nearestCity.population)) - 4, 0.25, 2.2) : 0.5;
      const normalizedUrbanDistance = nearestCity ? dCity / Math.max(6, nearestCity.urbanRadius) : 99;
      const cityClusterBoost = nearestCity ? clamp(1 - normalizedUrbanDistance) * (0.18 + populationScale * 0.16) : 0;
      const density = populationDensity[i];
      const oldTownScore = oldCoreInfluence[i] * 0.64 + premodernInfluence[i] * 0.32 + plain[i] * 0.12 + density * 0.08;
      const terrainUrbanPenalty = slope[i] * 1.02 + ridgeField[i] * 0.55 + Math.max(0, elevation[i] - 0.56) * 0.56;
      const nodeCausalPull = Math.max(stationInfluence[i] * 0.18, premodernInfluence[i] * 0.13, coastalLowland[i] * river[i] * 0.12, valleyField[i] * 0.08);
      const satelliteEnvelope = satelliteInfluence[i] * 0.54;
      const urbanEnvelope = cityInfluence[i] * 0.58 + cityCoreInfluence[i] * 0.3 + satelliteEnvelope + density * 0.47 + stationInfluence[i] * 0.18 + oldCoreInfluence[i] * 0.14 + newTownInfluence[i] * 0.12 + cityClusterBoost + nodeCausalPull - terrainUrbanPenalty;
      const coreScore = cityCoreInfluence[i] * 0.74 + urbanEnvelope * 0.3 + density * 0.36 + satelliteInfluence[i] * 0.16 + stationInfluence[i] * 0.06 + railInfluence2[i] * 0.04 - slope[i] * 0.82 - ridgeField[i] * 0.28;
      const suburbScore = urbanEnvelope * 0.54 + density * 0.14 + satelliteInfluence[i] * 0.22 + stationInfluence[i] * 0.09 + roadInfluence[i] * 0.05 + railInfluence2[i] * 0.05 + plain[i] * 0.16 + valleyField[i] * 0.04 + populationScale * 0.05 + (coreScore < 0.58 ? 0.05 : 0) - slope[i] * 0.76 - ridgeField[i] * 0.22;
      const roadsideScore = interchangeInfluence[i] * 0.54 + logisticsInfluence[i] * 0.18 + roadInfluence[i] * 0.1 + plain[i] * 0.1 - cityInfluence[i] * 0.02;
      const isolatedCorridor = roadInfluence[i] > 0.22 && cityInfluence[i] < 0.08 && stationInfluence[i] < 0.08 && interchangeInfluence[i] < 0.18;
      const ruralScore = villageInfluence[i] * 0.3 + agriculture[i] * 0.38 + plain[i] * 0.18 - slope[i] * 0.08;

      if (mountain && !urbanFootprint[i] && !oldUrbanFootprint[i]) landuse[i] = LANDUSE.FOREST;
      else if (urbanCoreFootprint[i] || (coreScore > 0.72 && density > 0.52 && slope[i] < 0.26 && ridgeField[i] < 0.38)) landuse[i] = LANDUSE.CBD;
      else if (oldUrbanFootprint[i] || oldTownScore > 0.49) landuse[i] = LANDUSE.OLD_URBAN;
      else if (industrialInfluence[i] > 0.48 && !urbanCoreFootprint[i]) landuse[i] = LANDUSE.INDUSTRIAL;
      else if (logisticsInfluence[i] > 0.44 && !urbanCoreFootprint[i]) landuse[i] = LANDUSE.LOGISTICS;
      else if (newTownInfluence[i] > 0.42 && urbanEnvelope > 0.16) landuse[i] = LANDUSE.NEW_TOWN;
      else if (urbanFootprint[i]) landuse[i] = LANDUSE.SUBURB;
      else if (suburbScore > 0.235 && !isolatedCorridor && slope[i] < 0.32 && ridgeField[i] < 0.48 && (normalizedUrbanDistance < 1.42 || satelliteInfluence[i] > 0.24)) landuse[i] = LANDUSE.SUBURB;
      else if (roadsideScore > 0.5 && plain[i] > 0.18 && slope[i] < 0.34 && ridgeField[i] < 0.5 && !isolatedCorridor && (interchangeInfluence[i] > 0.24 || logisticsInfluence[i] > 0.16 || cityInfluence[i] > 0.09)) landuse[i] = LANDUSE.ROADSIDE;
      else if (farm) landuse[i] = LANDUSE.FARMLAND;
      else if (ruralSettlementFootprint[i] || ruralScore > 0.3) landuse[i] = LANDUSE.RURAL;
      else landuse[i] = LANDUSE.RURAL;
    }
  }

  function hasUrbanNeighborCluster(x, y, radius = 2, minUrban = 7) {
    let urban = 0;
    for (let dy = -radius; dy <= radius; dy++) {
      for (let dx = -radius; dx <= radius; dx++) {
        const nx = x + dx;
        const ny = y + dy;
        if (!inside(nx, ny)) continue;
        const lu = landuse[indexOf(nx, ny)];
        if (isUrbanResidentialLanduse(lu)) urban++;
      }
    }
    return urban >= minUrban;
  }

  function removeIsolatedUrbanPatches(maxCells = 22) {
    const seen = new Uint8Array(SIZE);
    const namedCenters = [...modernCities, ...(satelliteCities || []), ...markets, ...ports, ...newTowns, ...stations];
    const queue = [];
    for (let i = 0; i < SIZE; i++) {
      if (seen[i] || !isHumanRegionCell(i)) continue;
      const lu0 = landuse[i];
      if (!isBuiltLanduse(lu0)) continue;
      const component = [];
      let maxDensity = 0;
      queue.length = 0;
      queue.push(i);
      seen[i] = 1;
      for (let q = 0; q < queue.length; q++) {
        const cur = queue[q];
        component.push(cur);
        maxDensity = Math.max(maxDensity, populationDensity[cur]);
        const [x, y] = xyOf(cur);
        for (const [nx, ny] of neighbors8(x, y)) {
          const ni = indexOf(nx, ny);
          if (seen[ni] || !isHumanRegionCell(ni)) continue;
          if (!isBuiltLanduse(landuse[ni])) continue;
          seen[ni] = 1;
          queue.push(ni);
        }
      }
      if (component.length > maxCells) continue;
      let hasAnchor = false;
      for (const ci of component) {
        const [x, y] = xyOf(ci);
        if (distanceToNearest(namedCenters, x, y) <= 5.8) {
          hasAnchor = true;
          break;
        }
      }
      if (!hasAnchor) {
        for (const ci of component) landuse[ci] = agriculture[ci] > 0.34 ? LANDUSE.FARMLAND : LANDUSE.RURAL;
      }
    }
  }

  for (let pass = 0; pass < 2; pass++) removeIsolatedUrbanPatches(36);

  // CBD is no longer a marker. It is a DID-like contiguous high-density core:
  // first remove isolated core cells, then grow connected high-density cells
  // from each urban center according to population scale.
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (landuse[i] === LANDUSE.CBD && !hasUrbanNeighborCluster(x, y, 2, 8)) landuse[i] = LANDUSE.SUBURB;
    }
  }

  function growDidCore(center, city, salt) {
    if (!center || !city) return 0;
    const start = indexOf(center.x, center.y);
    if (!isHumanRegionCell(start)) return 0;
    if ((city.population || 0) < 220000) return 0;
    const targetCells = Math.round(clamp(2 + Math.sqrt(city.population || 80000) / 74, 4, 22));
    const maxRadius = clamp((city.coreRadius || 3) * 2.4 + Math.sqrt(city.population || 80000) / 260, 6, 16);
    const selected = new Set();
    const queued = new Set([start]);
    const heap = new MinHeap();
    heap.push({ i: start, f: -10 });
    let made = 0;

    while (heap.length > 0 && made < targetCells) {
      const cur = heap.pop();
      if (!cur || selected.has(cur.i)) continue;
      const [x, y] = xyOf(cur.i);
      const i = cur.i;
      const d = Math.hypot(x - center.x, y - center.y);
      const support = populationDensity[i] * 1.18 + cityInfluence[i] * 0.22 + stationInfluence[i] * 0.18 + plain[i] * 0.12 - slope[i] * 1.24 - ridgeField[i] * 0.54 - Math.max(0, elevation[i] - 0.58) * 0.50 - floodplain[i] * 0.08 - d / maxRadius * 0.22;
      if (d > maxRadius || support < 0.44 || !isHumanRegionCell(i)) continue;
      if (!(landuse[i] === LANDUSE.OLD_URBAN || landuse[i] === LANDUSE.CBD || landuse[i] === LANDUSE.SUBURB || landuse[i] === LANDUSE.NEW_TOWN || populationDensity[i] > 0.22 || stationInfluence[i] > 0.14)) continue;

      selected.add(i);
      landuse[i] = LANDUSE.CBD;
      made++;

      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        if (queued.has(ni) || selected.has(ni) || !isHumanRegionCell(ni)) continue;
        const nd = Math.hypot(nx - center.x, ny - center.y);
        if (nd > maxRadius + 1) continue;
        const score = populationDensity[ni] * 1.24 + cityInfluence[ni] * 0.22 + stationInfluence[ni] * 0.18 + plain[ni] * 0.12 - slope[ni] * 1.25 - ridgeField[ni] * 0.54 - nd / maxRadius * 0.22 + hash2(nx, ny, seed + salt) * 0.03;
        queued.add(ni);
        heap.push({ i: ni, f: -score });
      }
    }
    return made;
  }

  urbanCenters.forEach((center, n) => growDidCore(center, center.parent || modernCities[n], 9400 + n * 17));
  for (let pass = 0; pass < 3; pass++) removeIsolatedUrbanPatches(42);
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (landuse[i] === LANDUSE.CBD && !hasUrbanNeighborCluster(x, y, 2, 8)) landuse[i] = LANDUSE.SUBURB;
    }
  }


  return {
    ports,
    crossings,
    passes,
    settlementCluster,
    settlementScore,
    villages,
    markets,
    castles,
    premodernRoads,
    minorRoads,
    castleTowns,
    modernCities,
    populationDensity,
    railways,
    branchRailways,
    ringRailways,
    externalRailways,
    stations,
    industrialZones,
    nationalRoads,
    ringRoads,
    expressways,
    ringExpressways,
    icAccessRoads,
    externalRoads,
    externalExpressways,
    interchanges,
    logisticsParks,
    satelliteCities,
    newTowns,
    landuse,
    stationInfluence,
    roadInfluence,
    railInfluence2,
    villageInfluence,
    externalGateways,
    transportDebug,
    cityPopulationCap,
  };
}