hspZ: GPU-accelerated HSPs backend
hspZ v0.0.4

C. Hit discovery

Count seed occurrences, preserve MAX_HITS chunk boundaries, and expand reference/query anchors in oracle order.

The cumulative hit array is both a total and an address map. It is therefore part of exactness, not temporary bookkeeping.

seed count       0   3   1   0   5
inclusive prefix 0   3   4   4   9

seed 1 writes    [2] [1] [0]     ← backwards from prefix 3
seed 2 writes                [3]
seed 4 writes    [8] [7] [6] [5] [4]
C.1 — Count each seed's reference occurrencesStep 1 of 5
React Flow

WHAT GOES IN: Query seed offsets and the reference index table. WHAT HAPPENS: `find_num_hits` subtracts adjacent cumulative bucket ends for each seed. WHAT COMES OUT: One reference-hit count per query seed. INVARIANT: Empty buckets contribute zero without changing seed order.

Readable journey

Step 1

C.1 — Count each seed's reference occurrences

what goes in
Query seed offsets and the reference index table.
what happens
`find_num_hits` subtracts adjacent cumulative bucket ends for each seed.
what comes out
One reference-hit count per query seed.
invariant
Empty buckets contribute zero without changing seed order.
  1. Query seeds → find_num_hits: packed k-mer per seed
  2. find_num_hits → Index table: bucket end - previous end
  3. find_num_hits → Hit counts: count in seed order
Step 2

C.2 — Turn counts into global prefixes

what goes in
Per-seed hit counts on the device.
what happens
Device blocks scan locally, the host exclusively scans small block sums, then offsets are added back. A batch whose total exceeds the 32-bit hit index is refused, never wrapped.
what comes out
A global inclusive prefix and the total number of anchors.
invariant
Only block sums cross the bus; when MAX_HITS cuts a batch, only the 256-count block holding each boundary is copied as well.
  1. Device counts → Block scan: 256 counts per block
  2. Block scan → Host block prefix: one sum per block
  3. Host block prefix → Block scan: exclusive block offsets
  4. Block scan → Global prefix: inclusive seed prefixes
Step 3

C.3 — Cut at MAX_HITS boundaries

what goes in
Global hit prefixes, total hits, and the resolved MAX_HITS cap.
what happens
The common case uses one chunk; oversized streams reproduce KegAlign's lower-bound walk over seed prefixes, located through the block sums (`HSPZ_CHUNK_WALK=full` copies the whole prefix array instead).
what comes out
Contiguous seed ranges and their exact hit ranges.
invariant
A seed is never split across chunks, so a chunk may exceed MAX_HITS by one seed bucket.
  1. Total hits → Lower-bound walk: MAX_HITS cap
  2. Lower-bound walk → Seed prefixes: find preceding seed boundary
  3. Seed prefixes → Lower-bound walk: cumulative hit position
  4. Lower-bound walk → Hit chunks: seed range + hit range
Step 4

C.4 — Choose sparse or dense expansion

what goes in
One chunk's seed count and hit count.
what happens
Sparse launches assign one thread per seed; at 16 or more hits per seed, one warp walks a seed bucket coalescently.
what comes out
The same packed anchors under either execution mapping.
invariant
Mapping changes memory traffic, not prefix-derived destinations.
  1. Hit chunk → Density choice: hits / seeds
  2. Density choice → Thread-per-seed: sparse path
  3. Density choice → Warp-per-seed: density at least 16
  4. Thread-per-seed → Anchors: prefix destinations
  5. Warp-per-seed → Anchors: same destinations
Step 5

C.5 — Expand buckets backwards

what goes in
Seed prefixes, reference bucket positions, query positions, and chunk offsets.
what happens
Every reference occurrence becomes a packed anchor written backwards from that seed's exclusive prefix.
what comes out
An anchor array in the raw order expected by the oracle.
invariant
Reversing the per-seed bucket walk is load-bearing because this order feeds later MAX_HITS and dedup scopes.
  1. One query seed → Reference bucket: bucket start..end
  2. Reference bucket → Exclusive prefix: occurrence j
  3. Exclusive prefix → Anchor array: prefix - 1 - chunk start - j