miniprot (Rust port)

Current H800 end-to-end result: Rust does not yet beat the C oracle on the 2000-query RefSeq map. This figure tracks the remaining gap, not a claimed win.

This project is a Rust port of lh3/miniprot.

However, it is not a trivial rewrite of the original C code, and certainly not an AI-generated line-by-line translation of the upstream repository. In this port, I intentionally removed many legacy C patterns and favored implementations that are more direct, maintainable, and idiomatic in Rust.

• Preserve the core mathematical model and algorithmic ideas of the original miniprot.
• Reorganize the implementation in Rust instead of reproducing the original C control flow, memory layout, macro style, and historical baggage line by line.
• Prioritize clarity, maintainability, and performance over source-level similarity to the upstream code.

• This project does not aim to be a bit-for-bit clone of the C implementation.
• It does not guarantee that internal intermediate states, floating-point details, iteration order, tie-breaking behavior, parallel scheduling, or final outputs will exactly match the upstream implementation on every input.
• Preserving legacy C coding patterns is not itself a compatibility goal.

This project only aims to stay consistent with the original implementation at the mathematical and algorithmic level:

• the problem definition should remain aligned
• the core algorithmic intent should remain aligned
• the main scoring, chaining, and alignment logic should remain aligned

It does not guarantee computational identity or byte-for-byte identity. In other words, even when two implementations follow the same algorithmic ideas, the exact numerical path, edge-case behavior, output details, and sometimes even parts of the final result may differ. This is an intentional design boundary, not by itself a bug definition.

cargo build --release

Benchmarked against upstream C oracle (lh3/miniprot v0.18-r281). Both binaries are compiled with arch-native flags (-C target-cpu=native / -march=native -mtune=native) unless noted otherwise.

This is the primary end-to-end GPU integration benchmark. It measures full mapping, including index load, candidate mapping, extension DP, traceback, ranking, and PAF formatting.

Dataset:

• Reference index: /root/grch38_gpu_test.mpi, built from /autodl-fs/data/refseq/GCF_000001405.40_GRCh38.p14_genomic.fna.gz
• Query: first 2000 records from /autodl-fs/data/refseq/uniprot_le100aa/GCF_000001405.40_Homo_sapiens_tax9606_uniprotkb_len001-100.fasta.gz
• Threads: -t4
• CUDA build: RUSTFLAGS="-C target-cpu=native" CUDA_HOME=/usr/local/cuda-12.8 cargo build --release --features cuda
• The H800 was shared during the final run; nvidia-smi showed a Python process using about 6.4 GB VRAM and 46-94% GPU util, so CUDA timings are conservative.

--gpu is part of the production mapping path, but it is deliberately gated. It batches eligible extension-DP jobs across queries and sends them to CUDA only when the shape is profitable (aa <= 128, nt <= 8192, batch >=4096). On the 2000-query RefSeq run, the mapper found only 53 left-extension CUDA candidates, so the GPU path now exits early and CPU SIMD handles the alignment. A looser experimental gate (nt <= 16384) did run CUDA extension kernels (5182 left and 5233 right jobs) but slowed the same workload to 25.17s, despite byte-identical Rust PAF output. The current baseline therefore keeps the stricter gate.

The chart is intentionally not a victory chart: on this larger real mapping case the Rust port still trails the C oracle. The gap is in alignment/traceback; candidate mapping alone remains faster in Rust (-A -t4: Rust 14.19s vs C 15.23s on the same 2000-query dataset).

Dataset:

• Reference: /autodl-fs/data/refseq/GCF_000001405.40_GRCh38.p14_genomic.fna.gz (GRCh38.p14, 928 MB gzipped FASTA)
• Query: first 200 records from /autodl-fs/data/refseq/uniprot_le100aa/GCF_000001405.40_Homo_sapiens_tax9606_uniprotkb_len001-100.fasta.gz
• Median of 3 measured samples; index outputs and the shared map index were written to /root scratch so shared-filesystem write latency is not counted as implementation time.
• Mapping uses the same C-built .mpi for both binaries.

On the same AutoDL H800 machine, the CUDA prepared-batch DP benchmark shows the GPU path's compute headroom for the short-window, large-batch shape it is meant to accept. These are isolated DP microbenchmarks, not end-to-end mapping claims:

The final rerun happened while another process was actively using the H800. Even under contention, isolated CUDA DP remained much faster than CPU SIMD; the production mapper still avoids CUDA for real RefSeq extension windows unless the full batch shape is profitable end to end.

On this broader real-data query set, PAF output is not byte-identical to the C oracle. Sorted -A output differs by 1 line out of 274; sorted default output has 266 common lines, 10 C-only lines, and 3 Rust-only lines. This project does not treat byte-for-byte identity with upstream C as a correctness guarantee, but the divergence should be considered when using the C oracle as a strict regression checker.

• Arch-native SIMD DP — NEON on aarch64, SSE4.1 on x86_64. Both match the C implementation's SIMD kernel (nasw-sse.c) with native intrinsics, not an SSE→NEON translation shim.
• Splice fast path — SIMD inner loop skips 12 splice-state operations per column when no splice signals are present, matching the C scalar path optimization.
• Fast approximate log2 — bit-manipulation log2 in chain scoring, identical to C's mp_log2, avoiding math library calls in the hot loop.
• zlib-rs — pure-Rust gzip backend for FASTA decompression during index build, with a linear unpack buffer replacing per-base closures.
• Radix sort — LSB radix sort for k-mer anchors, ~3–5× faster than comparison sort for the typical input sizes.
• Rayon parallelism — work-stealing thread pool for both index building and per-query mapping.
• Thread-local buffer reuse — profile matrix and SIMD scratch arrays reused across DP calls via thread-local storage, cutting ~1 GB of temporary allocations per benchmark run.
• GPU-accelerated DP — Metal and CUDA compute kernels with int16 arrays, pointer rotation, score row cache, optimised launch sizing, and prepared CUDA batches. Metal surpasses NEON SIMD at batch ≥4096 on Apple M2; CUDA on H800 surpasses CPU SIMD in isolated splice-free DP microbenchmarks. The end-to-end mapper uses CUDA only behind a profitability gate. Cross-platform via wgpu (Vulkan/Metal/DX12). See GPU optimization log for full experiment history.

GPU backends for batched DP computation. At batch ≥4096, Metal GPU surpasses NEON SIMD on Apple M2. On NVIDIA H800, the CUDA backend surpasses CPU SIMD in the splice-free DP sweep and reaches ~15× at batch 8192. Production --gpu mapping is stricter because real extension DP includes varied window lengths, splice penalties, host/device transfers, and CPU traceback work.

At batch 4096: Metal 80ms total vs NEON 110ms — GPU 1.4x faster. At batch 8192: Metal 169ms vs NEON 219ms — GPU 1.3x faster. Kernel-only per-call cost down to 18us (vs NEON 27us). GPU beats NEON because marginal per-call compute (11-24us) is below NEON's 27us, and fixed dispatch overhead (~40ms) is amortized at large batch.

CUDA is opt-in and requires the CUDA toolkit at build time:

CUDA_HOME=/usr/local/cuda-12.8 cargo test --release --features cuda bench_batch_size_sweep -- --nocapture

Defaults are tuned for H800 (CUDA_ARCH=sm_90, CUDA_THREADS=32) and can be overridden via environment variables.

At batch 8192: CUDA 21.74ms total vs CPU SIMD 341.96ms. The first batch pays CUDA context initialization, so small-batch latency is not representative. For repeated 8192-call batches, the CUDASW++-style prepared-batch path uploads nas/aas/params/matrix once and then reruns the resident workload with a best-observed steady-state average of 10.91ms (best iteration: 10.42ms), 32.0x faster than CPU SIMD on the same uncontended workload.

• int16 DP arrays — half stack footprint (1.8KB/thread), reduces VRAM spill
• Pointer rotation — swap 8 row pointers instead of copying 8×(al+1) shorts per row
• Score row cache — precompute score_row[j] once per nucleotide (every 3 rows), eliminating 2D VRAM lookup per column in the inner loop
• Thread-local aa buffer — copy aas[] to register array, avoid VRAM random read
• Conditional row_max — skip max tracking in non-extension mode
• TG=32 — per-threadgroup stack 58KB, fits in M2 GPU register file
• CUDA block size 32 — fastest H800 setting for large per-thread local arrays
• CUDA non-extension kernel — launches a specialized kernel for ext=false, removing h_best, row_max, and extension branches from the hot path
• Reusable CUDA buffers — keeps device buffers between batches and grows them only when needed, reducing repeated-batch allocation overhead
• Prepared CUDA batches — keeps a repeated batch's input and scoring matrix resident on the H800, reducing steady-state 8192-call batches from 13.16ms to 10.91ms average
• Production CUDA extension hook — batches eligible left/right extension DP across queries, uses donor/acceptor splice penalties, and falls back to CPU SIMD for small batches, long windows, CIGAR/traceback, or failed CUDA runs
• No-copy input buffers — new_buffer_with_bytes_no_copy wraps host slices directly on Apple Silicon unified memory
• Dead store elimination — remove per-row zero-fill that gets overwritten

GPU benchmarks:

cargo test --release --lib gpu_bench -- --nocapture
cargo test --release --lib bench_batch_size_sweep -- --nocapture
CUDA_HOME=/usr/local/cuda-12.8 cargo test --release --features cuda bench_batch_size_sweep -- --nocapture
CUDA_HOME=/usr/local/cuda-12.8 cargo test --release --features cuda bench_cuda_repeated_batch -- --nocapture

Performance baseline is the current Rust implementation. GPU regressions measured against current Metal shader, not against C oracle. CPU NEON SIMD is the single-query baseline; GPU DP is the high-batch baseline. The Metal and wgpu benchmark correctness claim covers the splice-free non-extension DP path (ext=false). CUDA additionally has the gated production extension-DP path used by --gpu; real end-to-end speed still decides whether that path is enabled for a batch.

CPU benchmarks (vs C oracle):

# index build
miniprot -t 4 -d human.mpi GCF_000001405.40_GRCh38.p14_genomic.fna.gz
# mapping
/usr/bin/time -p miniprot -t 4 human.mpi queries.fa > /dev/null

The oracle and parity tests depend on the upstream C implementation and its bundled fixtures. They are not required to build the Rust binary.

Use either of these setups:

git submodule update --init --recursive
cargo test

or point MINIPROT_C_ORACLE at an existing upstream miniprot binary.

If you need the original miniprot, use the upstream project directly: lh3/miniprot.

If you want a Rust reimplementation of miniprot that removes much of the legacy C style and is written as Rust rather than as translated C, that is what this repository is for.