test(benchmark): transform costs + GPU benchmarks for a reformulation cost model (stacked on #11)

[dougnukem]

Stacked on #11. Asks whether reformulating fastp for a GPU (or other batch accelerator) pays off, counting the cost of transforming data into and out of the accelerator's format. Write-up with results: benchmark/reformulations.md. Data: benchmark/results/reformulation-2026-09/.

Tools

• tools/soa_pack.cpp: CPU cost of index, pack (SoA / 2-bit) and emit, plus a representative QC kernel, on real reads. --dump writes a batch for the GPU check.
• tools/gpu_bench.py: PCIe, the same kernel on the GPU (checked read for read against the CPU), nvCOMP batched BGZF inflate, single-stream gzip (nvlzcat lookahead, nvCOMP ≥ 5.2), and GPU gzip compression.
• scripts/reformulation_model.py: projected CPU, GPU and wall time per reformulation from measured inputs. CPU: D1 BGZF, D2 parallel gzip reader, D3 -z 1, D4 uncompressed streaming, D5 flat-batch parser. GPU: A offload QC, B GPU-resident, C fused into a GPU aligner. QC-speedup sensitivity rows included.
• scripts/run_reformulation_suite.sh: every CPU-side measurement in one command.

Measured (n2d-highmem-48 + one L4, nvCOMP 5.3):

• Transform costs are small: index + pack + emit is 80–120 ns/read, against fastp's 670–1100 ns parse. 2-bit packing isn't worth it (4–5× the cost, ~30% fewer bytes).
• L4 QC kernel: 5.6 ns/read with 0 mismatches vs CPU over 1M reads, but moving the batch over PCIe costs 24.5 ns/read, 4× the kernel.
• L4 single-stream gzip inflate (nvlzcat): 1.26 GB/s. That's below 16 CPU threads of rapidgzip (3.4 GB/s). The nvCOMP Python default is 31 MB/s.
• L4 BGZF inflate 7.3 GB/s vs bgzip at 16 CPU threads 9.0 GB/s.
• L4 gzip at ratio 4.61 runs at 2.1 GB/s. At fastp's -z 4 ratio (4.88) it runs at 0.12 GB/s, about one CPU core.

Conclusion (details in reformulations.md):

• Largest wins are CPU changes: fix and speed up trimBySequence's gap search (33% of worker CPU, plus a correctness bug), a cheaper parser (1.2–2.1× on reader-bound data), right-sized -w, and uncompressed streaming to the aligner.
• A pays off for CPU-hours (−70–75% CPU) only if fastp's real per-read work runs ≥ ~10–15× a CPU thread on the GPU. That's unmeasured: the analogue kernel covers ~1/8 of the work.
• B is limited by single-stream inflate and GPU gzip's ratio.
• C is the one real GPU opportunity, but it needs a pre-alignment hook that Parabricks doesn't have.

Self-review fixes already in:

• nvCOMP Python API usage: enums, device sizes, freeing GPU memory, and the non-lookahead gzip path measured on a 100 MB member.
• soa_pack truncated-record over-read.
• B uses only the single-stream rate for ordinary gzip.
• SE BGZF pool formula.
• Python-API encode rows excluded from the model.

Remaining

• Deciding experiment for A: port the fixed trimBySequence plus statRead to the GPU and measure the real speedup.
• Measurement builds used branch measure/reformulation-build (this stack + fix: scale reader backpressure with worker thread count (#721) OpenGene/fastp#723) because upstream 1.3.7 deadlocks at -w 48.

🤖 Generated with Claude Code