#!/usr/bin/env node
/**
 * scripts/coarsen-regions.mjs — coarsens the committed `public/regions.geo.json`
 * in place, so the map is genuinely coarser, not just finer-grained.
 *
 * Rule (per user request: "focus on joining cities to the region around it.
 * Regions shouldn't be smaller than around 100x100km, excluding islands"):
 *
 *   FLOOR   — no region should be smaller than ~100×100 km = 10,000 km²,
 *             EXCEPT islands.
 *   ISLAND  — a region is an island if it has NO same-country neighbor: no
 *             other region with the same `properties.parent` whose bounding box
 *             is within CLOSE_DEG (≈≤22 km) of its own. Small islands stay
 *             as-is regardless of size (Corsica, the Hawaiian counties,
 *             Okinawa, the Faroes, island micro-states, …).
 *   MERGE   — a region that is small (area < FLOOR) AND not an island is
 *             merged into its NEAREST same-country neighbor (nearest by
 *             centroid distance; tie-break → larger area, then lexicographic
 *             id). Standalone cities (FR-75 Paris, the Latvian/Slovenian
 *             municipalities, the UK's 231 micro-districts) collapse into the
 *             region around them.
 *   DISSOLVE— the merge is a real boolean union (`@turf/union`, polygon-
 *             clipping): the shared boundary between the merged region and its
 *             host is removed, producing a single clean polygon — not two
 *             overlapping shapes. (Fallback if the union fails on a pair:
 *             combine the rings into a MultiPolygon without dissolving; the
 *             count of such fallbacks is reported.)
 *   ITERATE — passes run smallest-area-first (small regions merge into their
 *             neighbors, which grow and then absorb more). After each pass the
 *             area/neighbor/centroid metrics are recomputed; it stops at the
 *             fixed point (no small non-island region left) or at MAX_ITER.
 *   IDENTITY— the host (the region merged INTO) keeps its own `id` and
 *             `name`; the absorbed region's id is dropped. `properties.parent`
 *             is unchanged, so the map code (WorldMap/App/protocol/map.ts)
 *             keeps working untouched.
 *
 * Deterministic: fixed processing order (area, then id), deterministic
 * neighbor choice (distance, then -area, then id), deterministic boolean
 * union. Re-running on the coarsened file reaches the fixed point
 * immediately and rewrites a byte-identical file.
 *
 * Area model: the SAME planar approximation the validation test uses
 * (`scripts/regions.test.ts`) — `x = lon·111.32·cos(midLat)`,
 * `y = lat·110.574` km, shoelace per OUTER ring (holes ignored), each ring at
 * its own mean latitude — so the script's merge decisions and the test's
 * assertions measure the same quantity.
 *
 * Longitude handling: every ring is unwrapped to a continuous [0,360)-derived
 * space before its bbox/centroid is computed (the shortest-difference
 * "np.unwrap" scheme), which makes the one full-360° antimeridian-wrapping
 * feature (Antarctica's mainland ring) and any seam-straddling ring behave
 * correctly; cross-feature longitude gaps use the 360°-cyclic minimum.
 *
 * Output: features sorted by (parent, id) (the same comparator as
 * `scripts/generate-regions.mjs`), coordinates rounded to 4 decimal places,
 * single line of JSON + trailing newline. If the file would exceed the
 * ~5 MB guard, a light Douglas–Peucker pass (tolerances escalated through a
 * small table) brings it under ~4 MB before writing.
 *
 * Run: `npm run coarsen:regions`
 */

import { readFileSync, writeFileSync } from 'node:fs';
import path from 'node:path';
import { fileURLToPath } from 'node:url';
import { union } from '@turf/union';
import { feature, featureCollection } from '@turf/helpers';

const here = path.dirname(fileURLToPath(import.meta.url));
const root = path.resolve(here, '..');
const FILE = path.join(root, 'public', 'regions.geo.json');

// --- parameters ---------------------------------------------------------------
const AREA_FLOOR_KM2 = 10_000; // ~100×100 km floor for non-island regions
const CLOSE_DEG = 0.2; // "spatially close" bbox gap (both lat and lon ≤ this; ≈≤22 km)
const MAX_ITER = 20; // safety cap on fixed-point passes (spec: ~15)
const SIZE_GUARD_KB = 5_120; // keep the committed file under ~5 MB
const SIZE_TARGET_KB = 4_096; // …and land under ~4 MB once simplified
const DP_STEPS = [0.05, 0.1, 0.2, 0.3, 0.5, 0.75, 1.0]; // Douglas–Peucker tolerances (degrees)
const R4 = 4; // coordinate precision (decimal places, ≈11 m)
const SLIVER_EPS_KM2 = 1e-6; // drop only truly-degenerate (≈zero-area) clipper slivers

// --- area model (identical to scripts/regions.test.ts) -------------------------
const KM_LON = 111.32;
const KM_LAT = 110.574;
const rad = (deg) => (deg * Math.PI) / 180;

/** Shoelace area (km²) of one ring, measured at the ring's own mean latitude. */
function ringAreaKm2(ring) {
  let lat = 0;
  for (const [, la] of ring) lat += la;
  const kx = KM_LON * Math.cos(rad(ring.length ? lat / ring.length : 0));
  let a = 0;
  const n = ring.length;
  for (let i = 0; i < n; i++) {
    const [lon1, lat1] = ring[i];
    const [lon2, lat2] = ring[(i + 1) % n];
    a += lon1 * kx * (lat2 * KM_LAT) - lon2 * kx * (lat1 * KM_LAT);
  }
  return Math.abs(a / 2);
}

/** Sum of a geometry's outer-ring areas (holes are lakes — ignored by design). */
function geomAreaKm2(geom) {
  if (!geom) return 0;
  const polys = geom.type === 'Polygon' ? [geom.coordinates] : geom.coordinates;
  let total = 0;
  for (const poly of polys) total += ringAreaKm2(poly[0]);
  return total;
}

// --- longitude unwrapping (antimeridian-safe) ----------------------------------

/**
 * Unwrap a ring's longitudes to a continuous real line (the "np.unwrap"
 * shortest-difference scheme, anchored at the first vertex's [0,360) value).
 * Returns { lons, full } — `full` when the unwrapped span is ≥ 300°, i.e. the
 * ring is a full-360° antimeridian band (Antarctica's mainland ring).
 */
function unwrapRing(ring) {
  const n = ring.length;
  const norm = (x) => ((x % 360) + 360) % 360; // → [0,360)
  const u = new Array(n);
  u[0] = norm(ring[0][0]);
  for (let i = 1; i < n; i++) {
    // The PREVIOUS unwrapped value may lie outside [0,360) (it accumulates),
    // so normalize it back before measuring the difference; otherwise a bogus
    // ~180° jump is inserted and the ring's span balloons to ~180°.
    const a = norm(u[i - 1]);
    const b = norm(ring[i][0]);
    const d = (((b - a) % 360) + 540) % 360 - 180; // shortest signed diff in (-180, 180]
    u[i] = u[i - 1] + d;
  }
  let min = Infinity;
  let max = -Infinity;
  for (const v of u) {
    if (v < min) min = v;
    if (v > max) max = v;
  }
  return { lons: u, full: max - min >= 300 };
}

/**
 * Cyclic 1-D gap (degrees) between two tight longitude intervals in their own
 * unwrapped [0,360)-derived spaces. Tries the three 360° shifts of `b` so
 * different unwrap anchors reconcile; a `full` interval is close to everything.
 */
function lonGap(intA, intB) {
  if (intA.full || intB.full) return 0;
  let best = Infinity;
  for (const k of [-1, 0, 1]) {
    const a2 = intB.a + 360 * k;
    const b2 = intB.b + 360 * k;
    const gap = intA.b < a2 ? a2 - intA.b : b2 < intA.a ? intA.a - b2 : 0;
    if (gap < best) best = gap;
  }
  return best;
}

// --- geometry helpers ----------------------------------------------------------
const rN = (v, N) => Math.round(v * 10 ** N) / 10 ** N;

/** Round a ring to 4 dp, drop consecutive duplicates, and (re-)close it. */
function roundRing(ring, N = R4) {
  const out = [];
  for (const [x, y] of ring) {
    const p = [rN(x, N), rN(y, N)];
    const last = out[out.length - 1];
    if (!last || last[0] !== p[0] || last[1] !== p[1]) out.push(p);
  }
  if (
    out.length > 1 &&
    out[0][0] === out[out.length - 1][0] &&
    out[0][1] === out[out.length - 1][1]
  )
    out.pop();
  if (out.length < 3) return out; // degenerate; caller drops it
  return [...out, [...out[0]]];
}

/**
 * Round every ring of a geometry to 4 dp, drop degenerate rings (<3 pts) and
 * empty polygons, keep the Polygon/MultiPolygon type (a single-part
 * MultiPolygon collapses to Polygon). Returns null if nothing usable survived.
 */
function roundAll(geom) {
  if (!geom || (geom.type !== 'Polygon' && geom.type !== 'MultiPolygon')) return null;
  const polys = (geom.type === 'Polygon' ? [geom.coordinates] : geom.coordinates)
    .map((poly) => poly.map((ring) => roundRing(ring)))
    .filter((poly) => poly.length >= 1 && poly.every((r) => r.length >= 3));
  if (!polys.length) return null;
  if (polys.length === 1 && geom.type === 'Polygon') return { type: 'Polygon', coordinates: polys[0] };
  return { type: 'MultiPolygon', coordinates: polys };
}

/**
 * Normalize a union output (or fallback MultiPolygon): roundAll, then drop
 * truly-degenerate (≈zero-area, < SLIVER_EPS_KM2) polygons the clipper can
 * leave behind at shared boundaries. It never drops *every* part — a merge of
 * two sub-resolution (zero-area) micro-regions must still collapse into one
 * (degenerate) region rather than vanish. Returns null only when the input
 * itself had no usable ring (e.g. `null`).
 */
function normalizeGeom(geom) {
  const rounded = roundAll(geom);
  if (!rounded) return null;
  const all = rounded.type === 'Polygon' ? [rounded.coordinates] : rounded.coordinates;
  const kept = all.filter((poly) => ringAreaKm2(poly[0]) >= SLIVER_EPS_KM2);
  const polys = kept.length ? kept : all;
  if (polys.length === 1 && rounded.type === 'Polygon') return { type: 'Polygon', coordinates: polys[0] };
  return { type: 'MultiPolygon', coordinates: polys };
}

/** Douglas–Peucker of an OPEN point sequence (keeps endpoints). */
function dpOpen(pts, tol) {
  if (pts.length < 3) return pts;
  let maxD = 0;
  let idx = 0;
  const a = pts[0];
  const b = pts[pts.length - 1];
  for (let i = 1; i < pts.length - 1; i++) {
    const dx = b[0] - a[0];
    const dy = b[1] - a[1];
    const l2 = dx * dx + dy * dy;
    let t = l2 === 0 ? 0 : ((pts[i][0] - a[0]) * dx + (pts[i][1] - a[1]) * dy) / l2;
    t = Math.max(0, Math.min(1, t));
    const d = Math.hypot(pts[i][0] - (a[0] + t * dx), pts[i][1] - (a[1] + t * dy));
    if (d > maxD) {
      maxD = d;
      idx = i;
    }
  }
  if (maxD > tol) {
    const left = dpOpen(pts.slice(0, idx + 1), tol);
    const right = dpOpen(pts.slice(idx), tol);
    return left.concat(right.slice(1));
  }
  return [a, b];
}

/** Douglas–Peucker a CLOSED ring at tolerance `tol` (0 = skip). */
function dpRing(ring, tol) {
  if (tol <= 0 || ring.length < 4) return ring;
  let pts = ring;
  if (pts[0][0] === pts[pts.length - 1][0] && pts[0][1] === pts[pts.length - 1][1]) pts = pts.slice(0, -1);
  const s = dpOpen(pts, tol);
  if (s.length < 3) return ring;
  return [...s, [...s[0]]];
}

/** Round (and optionally DP-simplify) every ring of a geometry. */
function transformGeom(geom, tol) {
  const polys = (geom.type === 'Polygon' ? [geom.coordinates] : geom.coordinates)
    .map((poly) => poly.map((ring) => roundRing(ring, R4)).map((ring) => (tol > 0 ? dpRing(ring, tol) : ring)))
    .filter((poly) => poly.length >= 1 && poly.every((r) => r.length >= 3));
  if (!polys.length) return null;
  if (polys.length === 1 && geom.type === 'Polygon') return { type: 'Polygon', coordinates: polys[0] };
  return { type: 'MultiPolygon', coordinates: polys };
}

// --- region record -------------------------------------------------------------

/**
 * A region under management. `lonInt` is the tight longitude interval in
 * [0,360)-derived unwrapped space (or `full` for a 360° band); `(clat, clon)`
 * is the area-weighted centroid of the LARGEST outer ring (the region's
 * "location"), unwrapped for its ring then wrapped to [-180,180) for
 * cyclic-distance use.
 */
class Region {
  constructor(f) {
    this.id = f.id;
    this.name = f.properties.name;
    this.parent = f.properties.parent;
    this.geom = f.geometry;
    this.rawGeom = null; // last pre-normalization union output (safety net for build())
    this.dead = false;
    this.metrics();
  }

  metrics() {
    this.area = geomAreaKm2(this.geom);
    const polys = this.geom.type === 'Polygon' ? [this.geom.coordinates] : this.geom.coordinates;
    let latMin = Infinity;
    let latMax = -Infinity;
    let lonA = Infinity;
    let lonB = -Infinity;
    let lonFull = false;
    let bestRing = null;
    let bestArea = -1;
    for (const poly of polys) {
      const outer = poly[0];
      for (const [, la] of outer) {
        if (la < latMin) latMin = la;
        if (la > latMax) latMax = la;
      }
      const { lons, full } = unwrapRing(outer);
      if (full) lonFull = true;
      let a = Infinity;
      let b = -Infinity;
      for (const v of lons) {
        if (v < a) a = v;
        if (v > b) b = v;
      }
      lonA = Math.min(lonA, a);
      lonB = Math.max(lonB, b);
      const ar = ringAreaKm2(outer);
      if (ar > bestArea) {
        bestArea = ar;
        bestRing = outer;
      }
    }
    this.lonInt = { a: lonA, b: lonB, full: lonFull };
    this.latMin = latMin;
    this.latMax = latMax;
    // area-weighted centroid of the largest outer ring, in unwrapped space
    const { lons } = unwrapRing(bestRing);
    const lats = bestRing.map((p) => p[1]);
    const n = bestRing.length;
    let A = 0;
    let cx = 0;
    let cy = 0;
    for (let i = 0; i < n; i++) {
      const x1 = lons[i];
      const y1 = lats[i];
      const x2 = lons[(i + 1) % n];
      const y2 = lats[(i + 1) % n];
      const f = x1 * y2 - x2 * y1;
      A += f;
      cx += (x1 + x2) * f;
      cy += (y1 + y2) * f;
    }
    A /= 2;
    if (Math.abs(A) > 1e-9) {
      cx /= 6 * A;
      cy /= 6 * A;
    } else {
      cx = (lonA + lonB) / 2;
      cy = (latMin + latMax) / 2;
    }
    this.clon = (((cx + 180) % 360) + 360) % 360 - 180;
    this.clat = cy;
  }

  /** Cyclic planar distance (km) to another region's centroid. */
  distKm(o) {
    let dlon = Math.abs(this.clon - o.clon) % 360;
    if (dlon > 180) dlon = 360 - dlon;
    const midLat = rad((this.clat + o.clat) / 2);
    return Math.hypot(dlon * KM_LON * Math.cos(midLat), (this.clat - o.clat) * KM_LAT);
  }

  /** Bounding-box closeness to another same-space region (see CLOSE_DEG). */
  closeTo(o) {
    const latGap = this.latMax < o.latMin ? o.latMin - this.latMax : o.latMax < this.latMin ? this.latMin - o.latMax : 0;
    if (latGap > CLOSE_DEG) return false;
    return lonGap(this.lonInt, o.lonInt) <= CLOSE_DEG;
  }
}

/**
 * A "degenerate" region is a zero-area artifact: after 4-dp rounding its only
 * ring has <3 distinct points (a line/point, e.g. a micro-municipality that
 * v0.13's Douglas–Peucker collapsed to `[A, B, A]`). It contributes no area,
 * so it can never be dissolved by the clipper (which returns null) nor kept
 * through ring-dropping (which loses the 2-point ring). Merging one is purely
 * an absorption: the host keeps its own geometry.
 */
const isDegenerate = (r) => r.area < SLIVER_EPS_KM2;

// --- main ----------------------------------------------------------------------
function main() {
  const fc0 = JSON.parse(readFileSync(FILE, 'utf8'));
  const regions = fc0.features.map((f) => new Region(f));
  const byParent = new Map();
  for (const r of regions) {
    if (!byParent.has(r.parent)) byParent.set(r.parent, []);
    byParent.get(r.parent).push(r);
  }

  const t0 = Date.now();
  console.log(`input: ${fc0.features.length} regions across ${byParent.size} countries (${(Buffer.byteLength(readFileSync(FILE)) / 1024).toFixed(0)} KB)`);
  console.log(`parameters: floor ${AREA_FLOOR_KM2.toLocaleString('en-US')} km2, close ${CLOSE_DEG} deg, max ${MAX_ITER} passes`);

  const merges = []; // { from, fromName, fromArea, to, toName, toAreaAfter }
  let fallbackUnions = 0;

  // The fixed-point loop. Each pass: every region that is small AND has a
  // close same-country neighbor merges into the nearest such neighbor
  // (nearest by centroid distance; tie → larger area → lexicographic id).
  for (let pass = 0; pass < MAX_ITER; pass++) {
    let skippedIsland = 0;
    let skippedGrew = 0;
    const items = [];
    for (const r of regions) {
      if (r.dead || r.area >= AREA_FLOOR_KM2) continue;
      const siblings = byParent.get(r.parent);
      for (const s of siblings) {
        if (s !== r && !s.dead && r.closeTo(s)) {
          items.push(r);
          break;
        }
      }
    }
    if (!items.length) {
      console.log(`after pass ${pass}: fixed point reached (no small non-island regions left)`);
      break;
    }
    items.sort((a, b) => (a.area !== b.area ? a.area - b.area : a.id < b.id ? -1 : 1));
    let merged = 0;
    for (const r of items) {
      if (r.dead) continue;
      if (r.area >= AREA_FLOOR_KM2) {
        skippedGrew++;
        continue; // absorbed a smaller region earlier this pass and crossed the floor
      }
      const siblings = byParent.get(r.parent);
      // Close same-country neighbors as of NOW (mid-pass): some may have died
      // earlier in this pass.
      const close = [];
      for (const s of siblings) if (s !== r && !s.dead && r.closeTo(s)) close.push(s);
      if (!close.length) {
        skippedIsland++;
        continue; // lost its last close neighbor this pass → island now → exempt
      }
      // Prefer a REAL (non-degenerate) host so real regions dissolve cleanly;
      // fall back to any close neighbor only when no real one is nearby.
      const realCands = close.filter((s) => !isDegenerate(s));
      const pool = realCands.length ? realCands : close;
      let host = null;
      for (const s of pool) {
        if (
          !host ||
          r.distKm(s) < r.distKm(host) ||
          (r.distKm(s) === r.distKm(host) &&
            (s.area !== host.area ? s.area > host.area : s.id < host.id))
        )
          host = s;
      }
      if (isDegenerate(r)) {
        // r is a zero-area artifact (line/point): it contributes no area, so the
        // merge is a pure absorption — the host keeps its own geometry (a real
        // host is untouched; if the host is also degenerate a cluster of dots
        // collapses to one). Never route degenerate input through the clipper
        // (returns null) or ring-dropping (would lose the 2-point ring the map
        // already renders).
        r.dead = true;
        merged++;
        merges.push({ from: r.id, fromName: r.name, fromArea: r.area, to: host.id, toName: host.name, toAreaAfter: host.area });
        continue;
      }
      if (isDegenerate(host)) {
        // r is real but the chosen host is a zero-area artifact: the host takes
        // on r's geometry (keeping its own id/name). No clipper needed.
        host.geom = r.geom;
        host.rawGeom = r.geom;
        host.metrics();
        r.dead = true;
        merged++;
        merges.push({ from: r.id, fromName: r.name, fromArea: r.area, to: host.id, toName: host.name, toAreaAfter: host.area });
        continue;
      }
      // Both real: dissolve via a boolean union; the host keeps its id/name
      // and the item vanishes. @turf/union@7 dispatches through geomEach, which
      // only understands a FeatureCollection / single Feature / bare geometry —
      // a raw array of features has no `.type` and throws "Unknown Geometry
      // Type". Wrap it.
      let unionGeom = null;
      let rawUnion = null;
      try {
        const u = union(featureCollection([feature(r.geom), feature(host.geom)]));
        rawUnion = u ? u.geometry : null;
        unionGeom = normalizeGeom(rawUnion);
      } catch (err) {
        console.error(`  [warn] union failed for ${r.id} into ${host.id} (${err.message}); combining without dissolve`);
        fallbackUnions++;
      }
      if (!unionGeom) {
        // fallback: no dissolve — keep both shapes in one MultiPolygon
        const polys = [...(r.geom.type === 'Polygon' ? [r.geom.coordinates] : r.geom.coordinates), ...(host.geom.type === 'Polygon' ? [host.geom.coordinates] : host.geom.coordinates)];
        rawUnion = { type: 'MultiPolygon', coordinates: polys };
        unionGeom = normalizeGeom(rawUnion);
        if (!unionGeom) fallbackUnions++;
      }
      if (!unionGeom) {
        console.error(`  [error] could not merge ${r.id} into ${host.id}; leaving ${r.id} as-is`);
        continue;
      }
      host.geom = unionGeom;
      host.rawGeom = rawUnion;
      host.metrics();
      r.dead = true;
      merged++;
      merges.push({ from: r.id, fromName: r.name, fromArea: r.area, to: host.id, toName: host.name, toAreaAfter: host.area });
    }
    const alive = regions.filter((r) => !r.dead).length;
    console.log(
      `pass ${pass + 1}: ${items.length} small non-island region(s) → ${merged} merged` +
        ` (${skippedGrew} crossed the floor, ${skippedIsland} became island) → ${alive} regions remain`,
    );
    if (merged === 0) {
      // items exist but nothing could merge (all became islands this pass) → stable
      console.log(`pass ${pass + 1}: no merges possible — stable state`);
      break;
    }
  }

  // --- assemble output ----------------------------------------------------------
  const alive = regions.filter((r) => !r.dead);
  const sortById = (a, b) =>
    a.parent !== b.parent ? (a.parent < b.parent ? -1 : 1) : a.id.localeCompare(b.id, 'en', { numeric: true });
  alive.sort(sortById);

  const build = (tol) => {
    const feats = [];
    for (const r of alive) {
      // Never drop a region: if the simplified/normalized geometry degenerated,
      // fall back to the raw (pre-normalization) union output, then to the
      // untouched input geometry.
      const g =
        transformGeom(r.geom, tol) ??
        (r.rawGeom ? transformGeom(r.rawGeom, 0) ?? r.rawGeom : r.geom);
      feats.push({ type: 'Feature', id: r.id, properties: { name: r.name, parent: r.parent }, geometry: g });
    }
    return { type: 'FeatureCollection', features: feats };
  };

  let out = build(0);
  let json = JSON.stringify(out);
  let sizeKB = Buffer.byteLength(json) / 1024;
  console.log(`\ncoarsened: ${fc0.features.length} → ${out.features.length} regions (${fc0.features.length - out.features.length} merged away)`);

  let usedTol = 0;
  if (sizeKB > SIZE_GUARD_KB) {
    for (const t of DP_STEPS) {
      out = build(t);
      json = JSON.stringify(out);
      sizeKB = Buffer.byteLength(json) / 1024;
      usedTol = t;
      console.log(`DP ${t}: ${sizeKB.toFixed(1)} KB`);
      if (sizeKB < SIZE_TARGET_KB) break;
    }
  }
  writeFileSync(FILE, json + '\n');

  // --- summary -------------------------------------------------------------------
  const areaOf = new Map();
  const counts = new Map();
  let maxA = 0;
  let maxR = null;
  for (const f of out.features) {
    const a = geomAreaKm2(f.geometry);
    areaOf.set(f.id, a);
    counts.set(f.properties.parent, (counts.get(f.properties.parent) ?? 0) + 1);
    if (a > maxA) {
      maxA = a;
      maxR = f;
    }
  }
  const ccounts = [...counts.values()].sort((a, b) => a - b);
  const median =
    ccounts.length % 2 ? ccounts[(ccounts.length - 1) / 2] : (ccounts[ccounts.length / 2 - 1] + ccounts[ccounts.length / 2]) / 2;
  const maxCountry = [...counts.entries()].sort((a, b) => b[1] - a[1])[0];

  // islands = small regions with no close same-country neighbor in the FINAL set
  const finalByParent = new Map();
  for (const f of out.features) {
    if (!finalByParent.has(f.properties.parent)) finalByParent.set(f.properties.parent, []);
    finalByParent.get(f.properties.parent).push(f);
  }
  // rebuild final metrics (light: only need island flag on the < FLOOR set)
  const finalRegions = out.features.map((f) => new Region(f));
  const finalByParent2 = new Map();
  for (const r of finalRegions) {
    if (!finalByParent2.has(r.parent)) finalByParent2.set(r.parent, []);
    finalByParent2.get(r.parent).push(r);
  }
  let smallIslands = 0;
  let smallNonIslands = 0;
  const islandSamples = [];
  for (const r of finalRegions) {
    if (r.area >= AREA_FLOOR_KM2) continue;
    let hasNeighbor = false;
    for (const s of finalByParent2.get(r.parent)) {
      if (s !== r && r.closeTo(s)) {
        hasNeighbor = true;
        break;
      }
    }
    if (hasNeighbor) smallNonIslands++;
    else {
      smallIslands++;
      if (islandSamples.length < 12) islandSamples.push(`${r.id} ${r.name} (${Math.round(r.area)} km2)`);
    }
  }

  console.log(`wrote ${FILE} — ${sizeKB.toFixed(1)} KB (${(sizeKB / 1024).toFixed(2)} MB)` + (usedTol ? ` [DP @ ${usedTol}°]` : ''));
  console.log(`per-country count: min ${ccounts[0]} / median ${median} / max ${ccounts[ccounts.length - 1]} (max: ${maxCountry[0]})`);
  console.log(`max region: ${maxR.id} (${maxR.properties.name}) ≈ ${Math.round(maxA).toLocaleString('en-US')} km2`);
  console.log(`small (<${AREA_FLOOR_KM2.toLocaleString('en-US')} km2) in final set: ${smallIslands + smallNonIslands} — ${smallIslands} islands (exempt) / ${smallNonIslands} non-island residuals`);
  if (islandSamples.length) console.log(`island samples: ${islandSamples.join(', ')}`);
  console.log(`fallback (non-dissolved) unions: ${fallbackUnions}`);
  console.log(`merge log (${merges.length} merges): first 15:`);
  for (const m of merges.slice(0, 15)) {
    console.log(`  ${m.from} (${m.fromName}, ${Math.round(m.fromArea)} km2) → ${m.to} (${m.toName}, now ${Math.round(m.toAreaAfter)} km2)`);
  }
  const logPath = '/tmp/opencode/coarsen-log.txt';
  writeFileSync(logPath, merges.map((m) => `${m.from}\t${m.fromName}\t${Math.round(m.fromArea)}\t${m.to}\t${m.toName}\t${Math.round(m.toAreaAfter)}`).join('\n') + '\n');
  console.log(`full merge log → ${logPath}; elapsed ${((Date.now() - t0) / 1000).toFixed(1)}s`);
}

main();