GPU Multi-Source RAPTOR — Design

Status: design only, nothing executed. Target host: geocoder (RTX 5080, 16 GB VRAM, 13.2 GB free after H3 dataset). Authored 2026-06-10 on geo; sync this directory to geocoder before implementation (rsync -a /home/krish/gpu-pilot/raptor-design/ geocoder:/home/krish/gpu-pilot/raptor-design/).

Algorithm reference: Delling, Pajor, Werneck, "Round-Based Public Transit Routing" (Transportation Science, 2015). We implement standard RAPTOR semantics, batched over many origins.

External review: an independent design review (2026-06-10) is filed as REVIEW.md. Its verified numeric corrections are folded into this doc inline; its design deltas are summarized in §11. Read it before implementing — especially before any costing, multimodal, or tri-state work.

1. Problem statement

Compute a travel-time matrix T[origin][stop] = earliest arrival (seconds since departure) from every H3 res-9 cell centroid in the NYC tri-state grid to every transit stop, for a fixed departure time, with at most max_transfers transfers. Unreachable = sentinel.

Dimensions (NYC MTA Subway pilot):

Each origin is an independent single-source RAPTOR query. "Multi-source" means batching independent queries, not merging labels across sources. The replaced CPU pipeline (OTP polar sampling, 15 Java threads, 10 days) is the baseline.

2. Data layout (struct-of-arrays)

GPU memory holds (a) a small read-only timetable graph and (b) large per-origin state arrays. All times are int32 seconds since service-day midnight (GTFS allows >24:00:00, so values may exceed 86,400 — do not modulo). INF = 0x7F7F7F7F so a single-byte cudaMemset(ptr, 0x7F, bytes) fills arrays with a valid huge sentinel that is safe under atomicMin.

2.1 Read-only timetable (built host-side in cudf, uploaded once)

A RAPTOR "route" is a pattern: a maximal set of trips sharing an identical stop sequence. GTFS route_id is display-level only; the host preprocessor splits each GTFS route into patterns and sorts each pattern's trips by departure (FIFO/no-overtaking required for binary search — see risk R2).

The schedule is stored twice, in two layouts (8 MB total — trivial): dep_by_stop is transposed so the earliest-trip binary search at a fixed stop position reads a contiguous ascending array; arr_by_trip is trip-major so riding a boarded trip forward reads contiguous arrivals. Classic space-for-access-pattern trade; both slabs are int64-base-indexed so tri-state scale (100M+ stop_times) does not overflow.

Calendar handling: trips are filtered to one service date on the host (calendar + calendar_dates), so the kernel never sees service logic. service_calendar_t exists only as a host-side staging row.

2.2 Per-origin state (the big memory)

Layout rule: stop-major, origin-minor — tau[stop * origins_padded + origin]. A warp processes 32 consecutive origins, so every access at a fixed stop is one coalesced 128-byte transaction.

Marked bitmap layout: one u32 word per (origin-batch, stop); bit = lane. A warp tests "does any of my 32 origins mark this stop" with a single scalar load.

Footprint summary (the prompt's question): read-only graph ~30 MB; state ~7.6 GB at platform-level stops (S=1,000), ~3.9 GB at parent-station level (S=500). Per-origin cost is S x 14.25 B (3 i32 tau = 12 B, u16 output = 2 B, 2 mark bits = 0.25 B) ≈ 14.25 KB at S=1,000 — an earlier 16.3 KB figure double-counted travel_u16 as i32 (REVIEW.md §4.1). Origins per 13.2 GB chunk (85% usable, state only — keep headroom for graph + compaction temporaries): ~780K at S=1,000; ~1.5M at S=500. Res-9 fits in one chunk either way; res-10 (3.75M) needs ~5-6 chunks at S=1,000. At S=50K (tri-state) the per-origin cost is 712.5 KB -> ~15K origins/chunk (not 21K as previously stated). Output per departure sample: 535,488 x 1,000 x 2 B = 1.07 GB raw u16, far less in Parquet (RLE on INF-heavy rows); long-format reachable-only output is smaller still. Output unit: internal i32 seconds (GTFS granularity — ms adds nothing); exported column travel_sec u16 capped 65,534 = 18.2 h — acceptable for the subway pilot, a silent-clamp hazard for tri-state/overnight; width policy decision required before scaling (§11). 0xFFFF = unreachable. Consumers multiply by 1,000 if they want ms.

3. Multi-source RAPTOR on SIMD

Why no priority queue: RAPTOR is label-correcting by trip count, not by cost. Round k computes tau_k(p) = earliest arrival using ≤ k trips. Within a round, every update is min(old, candidate) — idempotent and commutative — so any execution order inside a round is correct, which is exactly the property SIMD wants. There is no node-duplication detection because state is a dense array indexed by (stop, origin), not an open list. This is the structural opposite of GPU A* (Zhou & Zeng 2015), which fights the queue and dedup hash; we have neither.

Per round k (1..K), three phases, grid-synchronized by separate kernel launches:

1. Copy-forward: tau_curr <- tau_prev (device-to-device, ~2 GB, ~3 ms). Preserves exact paper semantics (tau_k initialized to tau_{k-1}). A documented relaxation ("best-hopping": board from tau_best, skip the copy) breaks per-round Pareto snapshots and is not proven to converge to the same tau_best: a label improved during the current round can become boardable by another pattern in the same round — two boardings in one round, order-dependent (REVIEW.md §4.7). Off in v1; do not enable on empirical matching alone.
2. Route scan (raptor_round_scan): for each (origin-batch, pattern) work item, scan the pattern's stops in order. Per lane (= origin): maintain trip (current boarded trip, -1 = off). At each stop position i: (a) if on board, arr = arr_by_trip[trip][i]; if arr < tau_best[stop][o] then atomicMin(tau_curr), atomicMin(tau_best) (local pruning per the paper), atomicOr the mark bit; (b) if tau_prev[stop][o] ≤ departure of the current trip at i (or not on board), binary-search the earliest catchable trip and upgrade. This is the paper's inner loop verbatim, vectorized 32-wide.
3. Transfer relax (raptor_round_relax): for each (origin-batch, stop) with any mark from phase 2, push tau_curr[stop] + walk over footpath CSR with atomicMin, marking targets. Footpaths must be transitively closed in preprocessing (paper assumption). Revised per REVIEW.md §4.8: the "closure makes the read+extend-one-mask race benign" argument holds only if closure times are shortest-valid and every chained transfer is unconditionally interchangeable with its shortcut — v1 instead uses separate input/output mark buffers (trip_marked snapshot -> relax -> next_marked); one extra bitmap buys determinism.

Where atomics are required: two patterns may write the same (stop, origin) cell concurrently (hub stops served by many routes), and 32 origins share one mark word. atomicMin(int32) and atomicOr(u32) are native and cheap. No atomics are needed across origins on tau (each origin owns its column).

Marked-stop maintenance: scan reads marked_in (round k-1 result incl. transfers), writes marked_out; transfers read+extend marked_out; buffers swap each round. A device any_marked flag (4-byte D2H read per round) gives early exit when the frontier empties.

Work-skipping: a warp pre-passes the pattern's stop list, ORing mark words; if the union mask is 0, all 32 origins skip the pattern (return — whole warp exits, no divergence). v1 scans all (batch, pattern) pairs with this skip; a compacted frontier queue (Gunrock-style) is a v2 optimization, likely unnecessary because rounds are memory-bound, not launch-bound.

Back-of-envelope: full scan touches ~500 patterns x ~25 stops x 535K origins x ~12 B ≈ 80 GB/round; at ~850 GB/s effective on the 5080 that is ~100 ms/round, so ~0.5-1 s per departure sample for all 535K origins (vs 10 days CPU for the whole product). Estimates only; week-2 task is to measure.

4. Block/warp/thread granularity

Decision: one warp per (origin-batch-of-32, pattern); lane = origin. Block = 256 threads = 8 work items; grid = ceil(batches x patterns / 8).

Why, against alternatives: - One thread per origin, private full RAPTOR — no atomics, but every memory access diverges (different marked sets, different trips), so nothing coalesces; SIMD efficiency collapses. - One block per pattern, threads = stops — the scan carries a sequential dependency (trip state flows stop to stop); parallel-over-stops needs segmented-scan gymnastics. Too clever for v1. - Warp per (batch, pattern) — pattern topology and schedule reads are warp-uniform (one broadcast transaction); tau/mark accesses are perfectly coalesced (32 consecutive origins); divergence is confined to the binary search and gather reads of per-lane trips, which are short and bounded. Register pressure stays low (~30 registers/lane).

5. First-mile (walk) integration

Options: (a) precompute walk matrix [H3 cell x stops ≤ 1 km] with cuGraph SSSP over the OSM walk network, persist Parquet, load into Round 0; (b) inline crow-fly lookup inside Round 0 via spatial index.

Recommendation: (a), precomputed network walk matrix. Reasons: (1) NYC crow-fly lies — rivers, rail yards, and highway trenches make euclidean x 1.3 wrong by 10-30 minutes exactly where it matters; (2) the matrix is reusable across every departure sample, every GTFS refresh, and the last-mile direction — it amortizes to near-zero; (3) Round 0 becomes a trivial CSR read (origin -> (stop, walk_sec) pairs), keeping the kernel free of geometry; (4) the build is offline: ~1,000 SSSP runs seeded at station entrances (reverse direction), truncated at 1 km / 15 min. Note: cuGraph is not in the Phase 1 venv (cudf/cuspatial/cupy only) — installing cugraph-cu12==25.04 is an implementer task. Until then the pipeline ships a haversine x 1.35 / 1.33 m/s fallback (flagged walk_source='crowfly_fallback') so kernel bring-up is not blocked. Round 0 init: tau_0[stop][o] = depart_sec + walk_sec, mark stop; boarding first vehicles happens in round 1 — paper-exact.

6. Last-mile / cell-to-cell

Same matrix, transposed: stop -> cells within 1 km. The kernel's committed output stays [origin x stop]; cell-to-cell closes downstream as T_cell(o, d) = min over s (T[o][s] + walk(s, d)) — a cudf join/groupby-min over the long-format output, not kernel work. This keeps the kernel's output contract stable while letting hexmaps products choose destination granularity (stop, station, or cell). Same (a)-style precompute recommendation applies.

7. Multi-criteria handling

Recommendation: single-criterion arrival arrays, with per-round snapshots as the optional Pareto product. RAPTOR's round structure is the (arrival, transfers) Pareto frontier: tau_k(p) for k = 0..K is the full staircase. So we need no Pareto sets per stop — if a product wants "reachable within ≤ 1 transfer" isochrones, snapshot tau_curr after each round. Budget correction (REVIEW.md §4.2): resident i32 snapshots are 6 x 2.14 GB = 12.85 GB for k=0..5 — they do NOT fit beside the 7.6 GB state in 13.2 GB. Default is streaming each round's slab to host (zero extra resident cost); the alternative is finalizing each snapshot to u16 first (6 x 1.07 = 6.4 GB resident). v1 product output is tau_best only. If "transfer count of the best journey" is wanted later without snapshots: pack (arrival << 8) | round into u64 and atomicMin(u64) — lexicographic min keeps the pair consistent; doubles state memory; deferred to v1.1.

8. Validation plan

Three tiers, strictly ordered: - T1 — synthetic micro-feed, exact. 4-6 stops, 2 patterns, hand-traceable schedules (incl. an overtaking case and a >24:00 trip). GPU result must equal hand-computed labels bit-for-bit, per round. - T2 — CPU reference RAPTOR, exact. Plain-Python RAPTOR over the same preprocessed arrays (same pattern split, same service day). 1,000 random origins x all stops: GPU == CPU exactly — same algorithm, same data, integer arithmetic, no tolerance. - T3 — OTP cross-check, tolerance. Query OTP-northeast (http://gtfs3:8091/otp/gtfs/v1 GraphQL) for ~20 OD pairs, expect agreement within ±3 min (OTP differs on footpath models/boarding slack). Known-truth pairs, scoped to feed coverage: - Times Sq-42 St -> 34 St-Penn (1/2/3): 2-6 min, 0 transfers. - Times Sq -> Grand Central (S shuttle): 2-5 min. - 96 St -> South Ferry (1): 30-38 min. - Coney Island-Stillwell -> Times Sq (D/N/Q): 45-60 min. - Times Sq -> Howard Beach-JFK (A): 45-60 min. The prompt's "Times Sq -> JFK ≈ 50 min via A + AirTrain" includes the PANYNJ AirTrain feed, which is not in the MTA subway feed — validate the subway leg now, full pair when multi-feed lands. - 08854 -> 10001 (67-78 min via NJT NEC P&R): out of scope for subway-only; becomes the tri-state acceptance test when NJT rail joins (phase D).

9. Risk register

10. Sprint plan

Week 1 — correctness on small data (implementer agent): 1. Pattern builder in cudf (split + overtaking detect + trip sort) with asserts; service-date filter; transfer transitive closure. 2. Schedule slab builder (both layouts, int64 bases); dense reindex of stops/trips. 3. CPU reference RAPTOR (plain Python, ~100 lines) over the same arrays. 4. Compile kernels via cupy.RawModule (NVRTC); run T1 micro-feed; exact-match per round. 5. MTA feed, 1,000-origin subset; T2 exact match; fix divergences.

Week 2 — scale, walk legs, performance: 1. Full res-9 (535K origins) single-chunk run; finalize kernel -> u16 -> long-format Parquet writer. 2. Walk matrix: entrances from stops/pathways, cuGraph install + SSSP job, persist Parquet; wire into Round 0; keep crow-fly flag. 3. T3 OTP cross-check (~20 pairs); investigate any >3 min deltas. 4. Perf pass: per-round cudaEvent timings, Nsight memory-throughput check, warp-skip effectiveness; only then optimize. 5. Res-10 chunked driver (6 chunks); departure-window loop (e.g., 7-9 AM every 10 min); handoff notes.

11. Review deltas accepted into the plan (2026-06-10)

Full review: REVIEW.md. Numeric errata are already folded in above (§2.2, §7, R4). Design changes accepted:

1. Deterministic transfer relax — separate input/output mark buffers replace the read+extend-one-mask scheme (§3 phase 3; REVIEW §4.8).
2. Best-hopping demoted — from "converges, off in v1" to "unproven, hazardous; requires formal proof + adversarial same-round-chaining tests before any enablement" (§3 phase 1; REVIEW §4.7).
3. Boarding slack — RAPTOR_BOARD_SLACK_SEC = 0 is test-only. v1 ships a nonzero default (30-60 s) plus transfers.txt min_transfer_time; pathways-aware station topology is a later phase (REVIEW §4.5, §4.6).
4. Output width policy — u16 seconds saturates at 18.2 h; fine for the subway pilot, must be revisited (u16 minutes / u32 s / saturation flag) before tri-state or overnight service (REVIEW §4.4).
5. Tiled export — no full-matrix cp.nonzero; compaction is tiled/streamed (REVIEW §4.3).
6. Benchmark honesty — the acceptance comparison is GPU vs in-process optimized CPU RAPTOR over the same arrays, not vs the OTP API pipeline; report phase-by-phase timings (REVIEW §4.10, §4.11).
7. Scope framing — this kernel computes transit time only. The product ("how long AND how much, A to B, multimodal") additionally needs journey reconstruction, a separate cost engine (Fares v2, EPA MPG / EIA prices, tolls, parking, ride-hail), mode/asset state, backward arrive-by search, and time-based (not hex-ring) access limits. Those live in TODO.md phase E and REVIEW §5-6 — the kernel's output contract stays [origin x stop] earliest-arrival and is upstream input to that layer, not the product itself.
8. Architecture direction at scale — prefer sparse H3->stop access + reusable stop-to-stop profiles (rRAPTOR/ULTRA family) over re-running per H3 origin; never materialize hex-to-hex (535K² ≈ 573 GB/sample) (REVIEW §3, §5.5).