feat(08.1-01): implement CH contraction core — ordering, contraction, overlay witness, CSR (GREEN)

- contractGraph: greedy edgeDifference + λ·deletedNeighbours priority, lazy PQ with id tie-breaks, rank = contraction position (total order) - overlay-only witness search: contracted node + candidate absent, bounded by candidate distance and max(60, 4·localDegree) settles - directed shortcut emission (ascending from, then to), via midpoints, f32 CSR; dominated-insert guard; explosion tripwire fails loud at 4×N - deterministic: no PRNG/map-order; identical input -> identical output - 25-node grid: 41 shortcuts (1.64/node), 0 budget hits
2026-10-26 06:39:09 +08:00
parent a3e2c3253f
commit f75dc9e8c6
2 changed files with 545 additions and 64 deletions
@@ -0,0 +1,476 @@
+/*
+ * contraction.ts — the CH construction core (Phase 08.1, plan 01).
+ * ----------------------------------------------------------------------------
+ * One pure, typed-array module for everything build-time about the
+ * contraction hierarchy: deterministic node ordering (greedy
+ * edge-difference priority with the deleted-neighbour uniformity term,
+ * lazy binary PQ, ties to the lower node id), node contraction with an
+ * OVERLAY-ONLY witness search, and directed shortcut emission in CSR
+ * form. The fixture builder (test/fixtures/build-mini-graph.mjs), the
+ * real build script (scripts/build-road-graph.mjs, plan 02) and the
+ * graph.test.mjs ARG3 re-pin ALL compute their `ch:` blocks through
+ * this export — there is no second CH implementation anywhere
+ * (single-authority law).
+ *
+ * THE WEIGHT DOMAIN: CH runs on the REAL f32 seconds — no quantization.
+ * The hierarchy must preserve the exact edge weights the query reads
+ * back, because the query's reported seconds are re-summed from the
+ * unpacked original edges (ch-query.ts); a shortcut's stored weight is
+ * the f64 sum of its two constituents rounded once to f32 at CSR
+ * assembly — it steers the search only, it never enters reported
+ * seconds.
+ *
+ * TRAVERSAL LAW (mirrors astar.ts / landmark.ts verbatim — the A1
+ * precedent): an edge with bit0 set ("oneway-reverse", graph.ts) is
+ * NEVER relaxed — it never enters the overlay, so no witness search and
+ * no neighbour enumeration can traverse it, in either direction.
+ *
+ * THE OVERLAY LAW (research Pitfall 1 — the #1 correctness trap): the
+ * witness search sees ONLY the overlay graph — the node being
+ * contracted and every previously contracted node are absent (their
+ * incident edges are dead). A witness found through a contracted node
+ * does not exist in the reduced graph; honoring it would silently drop
+ * a required shortcut and break distance preservation. Bounding the
+ * search (candidate distance + settle budget) may add superfluous
+ * shortcuts but never removes a required one — correctness survives
+ * budget misses by construction (arXiv:1402.0402 §3).
+ *
+ * DETERMINISM: no randomness of any kind — no PRNG, no map-iteration
+ * order, no reliance on object key order. The priority PQ breaks ties
+ * by the lower node id, neighbour lists are sorted by node id, and
+ * identical input arrays give byte-identical output arrays (C7).
+ *
+ * Pure module: no I/O, no logging, no storage, zero value imports
+ * (T-06-05).
+ */
+
+/** The priority function's uniformity weight — the reference default:
+ *  p(u) = edgeDifference(u) + LAMBDA * deletedNeighbours(u). */
+export const CH_LAMBDA_DELETED_NEIGHBOURS = 1;
+/** The witness search's settle-budget floor (density-scaled above it). */
+export const CH_WPS_MAX_SETTLED = 60;
+/** The witness search's settle-budget scaling per unit of local degree. */
+export const CH_WPS_SETTLED_PER_DEGREE = 4;
+/** The shortcut-explosion tripwire: shortcuts per node — a build above
+ *  4×N fails loud (T-08.1-06: a blown hierarchy inflates the binary and
+ *  every query; the ratio names the diagnosis). */
+export const CH_SHORTCUT_EXPLOSION_RATIO = 4;
+
+/** The contraction compute input — the shape build scripts assemble. */
+export interface ContractionInput {
+  nodeCount: number;
+  /**
+   * Directed edges in any order. flags bit0 = oneway-reverse: relaxed by
+   * NO direction (mirrors the CSR's edgeFlags; see the traversal law).
+   */
+  edges: Array<{ from: number; to: number; seconds: number; flags: number }>;
+}
+
+/** The CH output: the contraction rank + the shortcut set in CSR form. */
+export interface ContractionResult {
+  /** rank[u] = u's position in the contraction order (0 = contracted
+   *  first). A STRICT total order — the query's strict uphill
+   *  comparisons depend on it (Pitfall 9). */
+  rank: Uint32Array;
+  /** CSR row pointers over the shortcut set — shortcutStart[u]..
+   *  shortcutStart[u+1] are u's out-shortcuts, length nodeCount + 1. */
+  shortcutStart: Uint32Array;
+  shortcutTarget: Uint32Array;
+  /** f32 round of w(from,via) + w(via,to) — steering weight only. */
+  shortcutSeconds: Float32Array;
+  /** The contracted midpoint of each shortcut — REQUIRED for the
+   *  query's recursive unpacking. */
+  shortcutVia: Uint32Array;
+  /** Build statistics (printed by the build script). */
+  stats: {
+    shortcutCount: number;
+    /** shortcuts per node — the explosion tripwire's input. */
+    ratioPerNode: number;
+    /** Witness searches run (priority evaluations + contractions). */
+    witnessSearches: number;
+    /** Witness searches bounded out by the settle budget. */
+    witnessBudgetHits: number;
+  };
+}
+
+function fail(context: string, message: string): never {
+  throw new Error(`${context}: ${message}`);
+}
+
+/* ————— overlay graph: one edge object shared by both adjacency lists ————— */
+
+interface OverlayEdge {
+  from: number;
+  to: number;
+  seconds: number;
+  alive: boolean;
+}
+
+/** Deterministic binary min-heap of (key, node) — ties break to the
+ *  LOWER node id at every level, so pop order is a pure function of the
+ *  key/node multiset (the determinism law). */
+class IdHeap {
+  private keys: Float64Array;
+  private nodes: Uint32Array;
+  size = 0;
+
+  constructor(capacity = 1024) {
+    this.keys = new Float64Array(capacity);
+    this.nodes = new Uint32Array(capacity);
+  }
+
+  private less(i: number, j: number): boolean {
+    const ki = this.keys[i]!;
+    const kj = this.keys[j]!;
+    return ki < kj || (ki === kj && this.nodes[i]! < this.nodes[j]!);
+  }
+
+  private swap(i: number, j: number): void {
+    const k = this.keys[i]!;
+    this.keys[i] = this.keys[j]!;
+    this.keys[j] = k;
+    const n = this.nodes[i]!;
+    this.nodes[i] = this.nodes[j]!;
+    this.nodes[j] = n;
+  }
+
+  push(key: number, node: number): void {
+    if (this.size === this.keys.length) this.grow();
+    let i = this.size;
+    this.size += 1;
+    this.keys[i] = key;
+    this.nodes[i] = node;
+    while (i > 0) {
+      const parent = (i - 1) >> 1;
+      if (!this.less(i, parent)) break;
+      this.swap(i, parent);
+      i = parent;
+    }
+  }
+
+  /** The min (key, node), or null when empty. */
+  pop(): { key: number; node: number } | null {
+    if (this.size === 0) return null;
+    const key = this.keys[0]!;
+    const node = this.nodes[0]!;
+    this.size -= 1;
+    if (this.size > 0) {
+      this.keys[0] = this.keys[this.size]!;
+      this.nodes[0] = this.nodes[this.size]!;
+      let i = 0;
+      for (;;) {
+        const left = 2 * i + 1;
+        if (left >= this.size) break;
+        const right = left + 1;
+        const child = right < this.size && this.less(right, left) ? right : left;
+        if (!this.less(child, i)) break;
+        this.swap(i, child);
+        i = child;
+      }
+    }
+    return { key, node };
+  }
+
+  private grow(): void {
+    const keys = new Float64Array(this.keys.length * 2);
+    keys.set(this.keys);
+    this.keys = keys;
+    const nodes = new Uint32Array(this.nodes.length * 2);
+    nodes.set(this.nodes);
+    this.nodes = nodes;
+  }
+}
+
+/** A neighbour (distinct node) with its minimum alive edge weight. */
+interface Neighbour {
+  node: number;
+  weight: number;
+}
+
+/**
+ * Distinct alive neighbours of the list's owner, ascending by node id,
+ * each carrying the MINIMUM alive edge weight (parallel edges collapse
+ * to the cheapest). `reverse=false` reads an out-adjacency (the
+ * neighbour is e.to); `reverse=true` reads an in-adjacency (the
+ * neighbour is e.from — the two lists share edge objects). Deterministic:
+ * collect, sort by (node, weight), keep the first per node.
+ */
+function collectNeighbours(adj: OverlayEdge[], reverse: boolean): Neighbour[] {
+  const pairs: Neighbour[] = [];
+  for (const e of adj) {
+    if (!e.alive) continue;
+    pairs.push({ node: reverse ? e.from : e.to, weight: e.seconds });
+  }
+  pairs.sort((a, b) => a.node - b.node || a.weight - b.weight);
+  const out: Neighbour[] = [];
+  for (const p of pairs) {
+    if (out.length > 0 && out[out.length - 1]!.node === p.node) continue;
+    out.push(p);
+  }
+  return out;
+}
+
+/** Count alive edges in an adjacency list (the incident-edge term). */
+function aliveCount(adj: OverlayEdge[]): number {
+  let n = 0;
+  for (const e of adj) if (e.alive) n += 1;
+  return n;
+}
+
+/**
+ * Build the contraction hierarchy for the whole graph: ordering,
+ * contraction, witness-bounded shortcut emission.
+ *
+ * Loud on malformed input (node ids out of range, non-finite/negative
+ * seconds, non-positive nodeCount) and on the shortcut-explosion
+ * tripwire — build-time failure is correct there; the fail-honest-null
+ * law belongs to decodeGraph.
+ */
+export function contractGraph(input: ContractionInput): ContractionResult {
+  const context = "contraction";
+  const n = input.nodeCount;
+  if (!Number.isInteger(n) || n < 1) fail(context, "nodeCount must be a positive integer");
+  // 1. Validate BEFORE any CSR work (the plan's step order): node ids in
+  //    range, finite non-negative seconds. Self-loops are skipped — no
+  //    cheapest path ever uses one and no shortcut pair can come from it.
+  const edgeCount = input.edges.length;
+  const overlay: OverlayEdge[] = [];
+  const outAdj: OverlayEdge[][] = new Array(n);
+  const inAdj: OverlayEdge[][] = new Array(n);
+  for (let i = 0; i < n; i++) {
+    outAdj[i] = [];
+    inAdj[i] = [];
+  }
+  for (let e = 0; e < edgeCount; e++) {
+    const edge = input.edges[e]!;
+    if (!Number.isInteger(edge.from) || edge.from < 0 || edge.from >= n) {
+      fail(context, `edge ${e} from-node ${edge.from} out of range`);
+    }
+    if (!Number.isInteger(edge.to) || edge.to < 0 || edge.to >= n) {
+      fail(context, `edge ${e} to-node ${edge.to} out of range`);
+    }
+    if (!Number.isFinite(edge.seconds) || edge.seconds < 0) {
+      fail(context, `edge ${edge.from}->${edge.to} has a non-finite/negative weight`);
+    }
+    if ((edge.flags & 1) !== 0) continue; // oneway-reverse: NEVER relaxed — not in the overlay
+    if (edge.from === edge.to) continue; // self-loop: no path can improve through it
+    const oe: OverlayEdge = { from: edge.from, to: edge.to, seconds: edge.seconds, alive: true };
+    overlay.push(oe);
+    outAdj[edge.from]!.push(oe);
+    inAdj[edge.to]!.push(oe);
+  }
+
+  // Original (input-graph) neighbour sets — the deleted-neighbour term
+  // counts neighbours ALREADY CONTRACTED regardless of overlay liveness
+  // (a contracted neighbour's edges are dead, so the overlay alone
+  // cannot answer it).
+  const origNeighbours: Array<Set<number>> = new Array(n);
+  for (let u = 0; u < n; u++) origNeighbours[u] = new Set();
+  for (const oe of overlay) {
+    origNeighbours[oe.from]!.add(oe.to);
+    origNeighbours[oe.to]!.add(oe.from);
+  }
+
+  // Witness search state: generation-stamped distances (no refills).
+  const dist = new Float64Array(n);
+  const visitGen = new Int32Array(n);
+  let gen = 0;
+
+  /**
+   * Bounded Dijkstra a -> b over ALIVE out-edges (the overlay: contracted
+   * nodes are unreachable because every incident edge is dead). Returns
+   * "witness" when b settles within `maxDist`; "budget" when the settle
+   * budget bound the search first; "none" otherwise (heap exhausted or
+   * the distance bound exceeded — no path within maxDist exists).
+   * `forbidden` is the node being contracted/evaluated: its edges are
+   * treated as absent even while they are still alive in the overlay —
+   * the search runs "with u removed" in BOTH call scopes (priority and
+   * contraction), so the trivial two-hop path THROUGH u can never
+   * witness u's own pairs.
+   */
+  function witness(
+    a: number,
+    b: number,
+    budgetSettles: number,
+    maxDist: number,
+    forbidden: number,
+  ): "witness" | "budget" | "none" {
+    gen += 1;
+    const heap = new IdHeap();
+    dist[a] = 0;
+    visitGen[a] = gen;
+    heap.push(0, a);
+    let settled = 0;
+    for (;;) {
+      const top = heap.pop();
+      if (top === null) return "none"; // exhausted: no path within maxDist
+      const u = top.node;
+      if (visitGen[u] === -gen) continue; // already settled (lazy deletion)
+      if (top.key > maxDist) return "none"; // distance bound: everything left is worse
+      visitGen[u] = -gen; // settle
+      settled += 1;
+      if (u === b) return "witness"; // settled within the candidate distance
+      if (settled >= budgetSettles) return "budget"; // bounded out — superfluous shortcut is SAFE
+      for (const e of outAdj[u]!) {
+        if (!e.alive) continue;
+        const v = e.to;
+        if (v === forbidden) continue; // the contracting node is REMOVED
+        if (visitGen[v] === -gen) continue;
+        const nd = top.key + e.seconds;
+        if (visitGen[v] !== gen || nd < dist[v]!) {
+          visitGen[v] = gen;
+          dist[v] = nd;
+          heap.push(nd, v);
+        }
+      }
+    }
+  }
+
+  /** Settle budget for contracting u: the density-scaled plan formula. */
+  function budgetFor(localDegree: number): number {
+    return Math.max(CH_WPS_MAX_SETTLED, CH_WPS_SETTLED_PER_DEGREE * localDegree);
+  }
+
+  const contracted = new Uint8Array(n);
+  const rank = new Uint32Array(n).fill(0xffffffff);
+  let witnessSearches = 0;
+  let witnessBudgetHits = 0;
+
+  /**
+   * priority(u) = edgeDifference(u) + LAMBDA * deletedNeighbours(u).
+   * edgeDifference runs the REAL witness searches (the honest cost of
+   * contracting u now) — recomputed fresh at every pop against the
+   * current overlay (the lazy-update law: a popped priority that
+   * IMPROVED is requeued; otherwise the node contracts at its popped
+   * key). Deterministic: neighbour lists are id-sorted, the witness
+   * heap ties break to the lower id.
+   */
+  function priorityOf(u: number): number {
+    const inN = collectNeighbours(inAdj[u]!, true);
+    const outN = collectNeighbours(outAdj[u]!, false);
+    const incident = aliveCount(inAdj[u]!) + aliveCount(outAdj[u]!);
+    let added = 0;
+    for (const a of inN) {
+      for (const b of outN) {
+        if (a.node === b.node) continue;
+        const candidate = a.weight + b.weight;
+        witnessSearches += 1;
+        const r = witness(a.node, b.node, budgetFor(inN.length + outN.length), candidate, u);
+        if (r === "budget") witnessBudgetHits += 1;
+        if (r !== "witness") added += 1;
+      }
+    }
+    let deleted = 0;
+    for (const v of origNeighbours[u]!) if (contracted[v] === 1) deleted += 1;
+    return added - incident + CH_LAMBDA_DELETED_NEIGHBOURS * deleted;
+  }
+
+  // The lazy priority queue: every node enters once; stale pops recompute
+  // and requeue on improvement (ties re-pop in id order — no loops: the
+  // recomputed value is not lower than itself).
+  const pq = new IdHeap(n > 0 ? n : 1);
+  for (let u = 0; u < n; u++) pq.push(priorityOf(u), u);
+
+  const shortcuts: Array<{ from: number; to: number; seconds: number; via: number }> = [];
+  let nextRank = 0;
+
+  for (;;) {
+    const top = pq.pop();
+    if (top === null) break;
+    const u = top.node;
+    if (contracted[u] === 1) continue; // defensive: no re-push ever happens
+    const fresh = priorityOf(u);
+    if (fresh < top.key) {
+      pq.push(fresh, u); // the priority IMPROVED — requeue (lazy update)
+      continue;
+    }
+    // 3. rank = the strictly increasing contraction position (never the
+    //    priority value) — a total order over all nodes.
+    contracted[u] = 1;
+    rank[u] = nextRank;
+    nextRank += 1;
+
+    // Neighbourhood BEFORE the removal (weights are the pre-kill minima).
+    const inN = collectNeighbours(inAdj[u]!, true);
+    const outN = collectNeighbours(outAdj[u]!, false);
+    // 6. Remove u and its incident overlay edges FIRST — every witness
+    //    search below runs with u absent (the overlay law), and every
+    //    previously contracted node is absent because its edges died at
+    //    its own contraction.
+    for (const e of outAdj[u]!) e.alive = false;
+    for (const e of inAdj[u]!) e.alive = false;
+
+    // 4./5. Candidate shortcuts a -> b, deterministic order: ascending a
+    //    (in-neighbours), then ascending b (out-neighbours).
+    const budget = budgetFor(inN.length + outN.length);
+    for (const a of inN) {
+      for (const b of outN) {
+        if (a.node === b.node) continue;
+        const candidate = a.weight + b.weight;
+        witnessSearches += 1;
+        const r = witness(a.node, b.node, budget, candidate, u);
+        if (r === "budget") witnessBudgetHits += 1;
+        if (r === "witness") continue;
+        // Dominated-insert guard: an alive a->b edge at or under the
+        // candidate already preserves the distance — emitting a parallel
+        // shortcut would only bloat the CSR. (The witness search itself
+        // usually catches this; the guard covers the bounded-out case.)
+        let dominated = false;
+        for (const e of outAdj[a.node]!) {
+          if (e.alive && e.to === b.node && e.seconds <= candidate) {
+            dominated = true;
+            break;
+          }
+        }
+        if (dominated) continue;
+        const se: OverlayEdge = { from: a.node, to: b.node, seconds: candidate, alive: true };
+        outAdj[a.node]!.push(se);
+        inAdj[b.node]!.push(se);
+        shortcuts.push({ from: a.node, to: b.node, seconds: candidate, via: u });
+      }
+    }
+  }
+
+  // 8. The explosion tripwire: fail loud naming the ratio and the count.
+  if (shortcuts.length > CH_SHORTCUT_EXPLOSION_RATIO * n) {
+    fail(
+      context,
+      `shortcut explosion — ${shortcuts.length} shortcuts exceeds ` +
+        `${CH_SHORTCUT_EXPLOSION_RATIO}x ${n} nodes ` +
+        `(check contraction order / witness budget)`,
+    );
+  }
+
+  // 7. CSR assembly: shortcuts were emitted ascending by (from, to); the
+  //    counting sort preserves that within each row.
+  const shortcutCount = shortcuts.length;
+  const shortcutStart = new Uint32Array(n + 1);
+  const shortcutTarget = new Uint32Array(shortcutCount);
+  const shortcutSeconds = new Float32Array(shortcutCount);
+  const shortcutVia = new Uint32Array(shortcutCount);
+  for (const s of shortcuts) shortcutStart[s.from + 1] = shortcutStart[s.from + 1]! + 1;
+  for (let i = 0; i < n; i++) shortcutStart[i + 1] = shortcutStart[i + 1]! + shortcutStart[i]!;
+  const cursor = shortcutStart.slice(0, n);
+  for (const s of shortcuts) {
+    const pos = cursor[s.from]!;
+    cursor[s.from] = pos + 1;
+    shortcutTarget[pos] = s.to;
+    shortcutSeconds[pos] = s.seconds; // the single f64 -> f32 rounding
+    shortcutVia[pos] = s.via;
+  }
+
+  // 9. Stats.
+  return {
+    rank,
+    shortcutStart,
+    shortcutTarget,
+    shortcutSeconds,
+    shortcutVia,
+    stats: {
+      shortcutCount,
+      ratioPerNode: shortcutCount / n,
+      witnessSearches,
+      witnessBudgetHits,
+    },
+  };
+}