use std::fs;
use std::io::{BufWriter, Write};
use std::path::PathBuf;
use std::process::Command;
use std::time::{Duration, Instant};
use slow5lib::aux::AuxMeta;
use slow5lib::aux::AuxValue;
use slow5lib::header::{Header, ReadGroup, RecordCompression, SignalCompression};
use slow5lib::record::Record;
use slow5lib::{Slow5IndexedReader, Slow5Reader, Slow5Writer};
struct Args {
file: PathBuf,
slow5tools: String,
n_correct: usize,
n_indexed: usize,
n_write: usize,
skip_correct: bool,
skip_seq: bool,
skip_par: bool,
skip_idx: bool,
skip_write: bool,
skip_write_par: bool,
skip_write_pipe: bool,
skip_c: bool,
}
impl Default for Args {
fn default() -> Self {
Self {
file: PathBuf::new(),
slow5tools: "slow5tools".into(),
n_correct: 1_000,
n_indexed: 10_000,
n_write: 10_000,
skip_correct: false,
skip_seq: false,
skip_par: false,
skip_idx: false,
skip_write: false,
skip_write_par: false,
skip_write_pipe: false,
skip_c: false,
}
}
}
fn parse_args() -> Args {
let mut args = Args::default();
let raw: Vec<String> = std::env::args().skip(1).collect();
let mut i = 0;
while i < raw.len() {
match raw[i].as_str() {
"--slow5tools" => {
i += 1;
args.slow5tools = raw[i].clone();
}
"--n-correct" => {
i += 1;
args.n_correct = raw[i].parse().expect("--n-correct integer");
}
"--n-indexed" => {
i += 1;
args.n_indexed = raw[i].parse().expect("--n-indexed integer");
}
"--n-write" => {
i += 1;
args.n_write = raw[i].parse().expect("--n-write integer");
}
"--skip-correct" => args.skip_correct = true,
"--skip-seq" => args.skip_seq = true,
"--skip-par" => args.skip_par = true,
"--skip-idx" => args.skip_idx = true,
"--skip-write" => args.skip_write = true,
"--skip-write-par" => args.skip_write_par = true,
"--skip-write-pipe" => args.skip_write_pipe = true,
"--skip-c" => args.skip_c = true,
other if !other.starts_with('-') => {
if args.file.as_os_str().is_empty() {
args.file = PathBuf::from(other);
} else {
eprintln!("unexpected argument: {other}");
std::process::exit(1);
}
}
other => {
eprintln!("unknown flag: {other}");
std::process::exit(1);
}
}
i += 1;
}
if args.file.as_os_str().is_empty() {
eprintln!("usage: harness <FILE> [OPTIONS]\n");
eprintln!("OPTIONS:");
eprintln!(" --slow5tools PATH path to slow5tools [slow5tools]");
eprintln!(" --n-correct N reads for correctness check [1000]");
eprintln!(" --n-indexed N reads for indexed benchmark [10000]");
eprintln!(" --n-write N reads for write benchmark [10000]");
eprintln!(
" --skip-correct / --skip-seq / --skip-par / --skip-idx / --skip-write / --skip-write-par / --skip-write-pipe / --skip-c"
);
std::process::exit(1);
}
args
}
fn signal_hash(signal: &[i16]) -> u64 {
signal.iter().fold(0u64, |h, &s| {
h.wrapping_mul(0x517cc1b727220a95)
.wrapping_add(s as u16 as u64)
})
}
fn strided_ids(all_ids: &[String], n: usize) -> Vec<String> {
if all_ids.len() <= n {
return all_ids.to_vec();
}
let stride = all_ids.len() / n;
(0..n).map(|i| all_ids[i * stride].clone()).collect()
}
fn write_id_list(ids: &[String], path: &str) {
let mut f = BufWriter::new(fs::File::create(path).expect("create id list"));
for id in ids {
writeln!(f, "{id}").expect("write id");
}
}
fn slow5tools_ok(st: &str) -> bool {
Command::new(st)
.arg("--version")
.output()
.map(|o| o.status.success())
.unwrap_or(false)
}
fn fmt_dur(d: Duration) -> String {
let s = d.as_secs_f64();
if s >= 3600.0 {
format!("{:.1}h", s / 3600.0)
} else if s >= 60.0 {
format!("{:.1}m", s / 60.0)
} else {
format!("{:.2}s", s)
}
}
fn row(label: &str, records: u64, bytes: u64, d: Duration, note: &str) {
let s = d.as_secs_f64();
let rps = records as f64 / s;
let gbps = bytes as f64 / s / 1e9;
let rps_str = if records > 0 {
format!("{rps:>10.0} reads/s")
} else {
" ".repeat(17)
};
let gbps_str = if bytes > 0 {
format!("{gbps:>6.2} GB/s")
} else {
" ".repeat(12)
};
println!(
" {:<22} {:>7} {rps_str} {gbps_str} {note}",
label,
fmt_dur(d)
);
}
fn separator() {
println!("{}", "-".repeat(72));
}
fn aux_approx_eq(a: &AuxValue, b: &AuxValue) -> bool {
match (a, b) {
(AuxValue::Float(x), AuxValue::Float(y)) => rel_close_f32(*x, *y),
(AuxValue::Double(x), AuxValue::Double(y)) => rel_close_f64(*x, *y),
_ => a == b,
}
}
fn rel_close_f32(a: f32, b: f32) -> bool {
if a == b {
return true;
}
let m = a.abs().max(b.abs());
if m == 0.0 {
return true;
}
((a - b) / m).abs() < 1e-4
}
fn rel_close_f64(a: f64, b: f64) -> bool {
if a == b {
return true;
}
let m = a.abs().max(b.abs());
if m == 0.0 {
return true;
}
((a - b) / m).abs() < 1e-6
}
fn run_correctness(args: &Args, ids: &[String]) {
let st = &args.slow5tools;
let stride = args.n_correct.max(1);
println!(
"\nCorrectness check ({} reads, stride ~{})",
ids.len(),
stride
);
separator();
println!(" Extracting via slow5tools get...");
let id_list = "/tmp/harness_correct_ids.txt";
write_id_list(ids, id_list);
let tmp = "/tmp/harness_correct_st.slow5";
let status = Command::new(st)
.args(["get", "--to", "slow5", "-l", id_list, "-o", tmp])
.arg(&args.file)
.stderr(std::process::Stdio::null())
.status()
.expect("spawn slow5tools get");
if !status.success() {
println!(" slow5tools get failed -- skipping correctness check");
return;
}
let mut st_reader = Slow5Reader::open(tmp).expect("open st correctness output");
let st_recs: std::collections::HashMap<String, Record> = st_reader
.records()
.map(|r| {
let rec = r.expect("st rec");
(rec.read_id.clone(), rec)
})
.collect();
let rust_reader = Slow5IndexedReader::open(&args.file).expect("open indexed");
let aux_names = rust_reader.header().aux_meta.names.clone();
let mut signal_mismatches = 0usize;
let mut primary_mismatches: Vec<String> = Vec::new();
let mut aux_mismatches: Vec<String> = Vec::new();
let mut running_hash = 0u64;
for id in ids {
let rust = rust_reader
.get(id)
.unwrap_or_else(|e| panic!("get {id}: {e}"));
let Some(st) = st_recs.get(id) else {
eprintln!(" warning: {id} missing from slow5tools output");
continue;
};
let rh = signal_hash(&rust.raw_signal);
let sh = signal_hash(&st.raw_signal);
running_hash ^= rh;
if rh != sh {
signal_mismatches += 1;
}
for (field, rust_val, st_val) in [
("digitisation", rust.digitisation, st.digitisation),
("offset", rust.offset, st.offset),
("range", rust.range, st.range),
("sampling_rate", rust.sampling_rate, st.sampling_rate),
] {
if !rel_close_f64(rust_val, st_val) {
primary_mismatches.push(format!("{id}/{field}: {rust_val} vs {st_val}"));
}
}
for name in &aux_names {
let rv = rust
.aux
.get(name.as_str())
.cloned()
.unwrap_or(AuxValue::Missing);
let sv = st
.aux
.get(name.as_str())
.cloned()
.unwrap_or(AuxValue::Missing);
if !aux_approx_eq(&rv, &sv) {
aux_mismatches.push(format!("{id}/{name}: rust={rv:?} st={sv:?}"));
}
}
}
let pass = |n: usize| if n == 0 { "PASS" } else { "FAIL" };
println!(
" Signal: {} (hash 0x{running_hash:016x}, {} mismatches)",
pass(signal_mismatches),
signal_mismatches
);
println!(
" Primary: {} ({} mismatches)",
pass(primary_mismatches.len()),
primary_mismatches.len()
);
println!(
" Aux fields: {} ({} fields, {} mismatches)",
pass(aux_mismatches.len()),
aux_names.len(),
aux_mismatches.len()
);
for m in primary_mismatches.iter().take(3) {
println!(" {m}");
}
for m in aux_mismatches.iter().take(5) {
println!(" {m}");
}
}
fn run_seq(args: &Args, file_bytes: u64, have_c: bool) {
println!("\nSequential read -- whole file, decompress all records");
println!(" Rust measures: I/O + block decompress + SVB-ZD decode + signal hash");
println!(" slow5tools measures: same + text encoding (extra overhead)");
separator();
let mut reader = Slow5Reader::open(&args.file).expect("open seq");
let mut records = 0u64;
let mut hash = 0u64;
let t0 = Instant::now();
for rec in reader.records() {
let rec = rec.expect("seq record");
hash ^= signal_hash(&rec.raw_signal);
records += 1;
}
let rust_dur = t0.elapsed();
row(
"Rust",
records,
file_bytes,
rust_dur,
&format!("hash 0x{hash:016x}"),
);
if have_c {
let t0 = Instant::now();
let status = Command::new(&args.slow5tools)
.args(["view", "--to", "slow5", "-t", "1"])
.arg(&args.file)
.stdout(std::process::Stdio::null())
.stderr(std::process::Stdio::null())
.status()
.expect("spawn slow5tools view");
let st_dur = t0.elapsed();
if status.success() {
row(
"slow5tools view -t 1",
records,
file_bytes,
st_dur,
"[+ text encoding]",
);
let ratio = st_dur.as_secs_f64() / rust_dur.as_secs_f64();
println!(" Ratio (st/rust): {ratio:.2}x");
}
}
}
#[cfg(feature = "rayon")]
fn run_par(args: &Args, file_bytes: u64) {
use rayon::prelude::*;
println!("\nParallel read -- whole file, rayon par_records");
println!(" Sequential I/O pass collects all raw records, then rayon decompresses in parallel");
separator();
let mut reader = Slow5Reader::open(&args.file).expect("open par");
let t0 = Instant::now();
let (records, hash): (u64, u64) = reader
.par_records()
.map(|r| {
let rec = r.expect("par record");
(1u64, signal_hash(&rec.raw_signal))
})
.reduce(|| (0, 0), |(ra, ha), (rb, hb)| (ra + rb, ha ^ hb));
let dur = t0.elapsed();
row(
"Rust par_records",
records,
file_bytes,
dur,
&format!("hash 0x{hash:016x}"),
);
}
fn run_idx(args: &Args, ids: &[String], have_c: bool) {
println!(
"\nIndexed random access -- {} reads, evenly strided across file",
ids.len()
);
println!(" Measures: index lookup + pread + block decompress + SVB-ZD decode");
separator();
let reader = Slow5IndexedReader::open(&args.file).expect("open indexed");
let mut hash = 0u64;
let t0 = Instant::now();
for id in ids {
let rec = reader
.get(id)
.unwrap_or_else(|e| panic!("indexed get {id}: {e}"));
hash ^= signal_hash(&rec.raw_signal);
}
let rust_dur = t0.elapsed();
row(
"Rust indexed",
ids.len() as u64,
0,
rust_dur,
&format!("hash 0x{hash:016x}"),
);
if have_c {
let id_list = "/tmp/harness_idx_ids.txt";
write_id_list(ids, id_list);
let t0 = Instant::now();
let status = Command::new(&args.slow5tools)
.args(["get", "--to", "slow5", "-l", id_list])
.arg(&args.file)
.stdout(std::process::Stdio::null())
.stderr(std::process::Stdio::null())
.status()
.expect("spawn slow5tools get");
let st_dur = t0.elapsed();
if status.success() {
row(
"slow5tools get",
ids.len() as u64,
0,
st_dur,
"[+ text encoding]",
);
let ratio = st_dur.as_secs_f64() / rust_dur.as_secs_f64();
println!(" Ratio (st/rust): {ratio:.2}x");
}
}
}
fn run_write(args: &Args, ids: &[String], have_c: bool) {
println!("\nWrite benchmark -- {} reads, zstd + SVB-ZD", ids.len());
println!(" Both tools: indexed read from source BLOW5, write new BLOW5 (zstd+svbzd)");
println!(" Measures: index lookup + pread + decompress + re-compress + write");
separator();
let id_list = "/tmp/harness_write_ids.txt";
write_id_list(ids, id_list);
let tmp_rust = "/tmp/harness_write_rust.blow5";
let header = blow5_header();
let src_reader = Slow5IndexedReader::open(&args.file).expect("open indexed for write");
let mut writer = Slow5Writer::create(tmp_rust, header).expect("create rust write output");
let t0 = Instant::now();
for id in ids {
let rec = src_reader
.get(id)
.unwrap_or_else(|e| panic!("write get {id}: {e}"));
writer.write(&rec).expect("write record");
}
writer.finish().expect("finish write");
let rust_dur = t0.elapsed();
let rust_bytes = fs::metadata(tmp_rust).map(|m| m.len()).unwrap_or(0);
row("Rust", ids.len() as u64, rust_bytes, rust_dur, "");
if have_c {
let tmp_st = "/tmp/harness_write_st.blow5";
let t0 = Instant::now();
let status = Command::new(&args.slow5tools)
.args([
"get", "--to", "blow5", "-c", "zstd", "-s", "svb-zd", "-t", "1", "-l", id_list,
"-o", tmp_st,
])
.arg(&args.file)
.stderr(std::process::Stdio::null())
.status()
.expect("spawn slow5tools get for write bench");
let st_dur = t0.elapsed();
if status.success() {
let st_bytes = fs::metadata(tmp_st).map(|m| m.len()).unwrap_or(0);
row(
"slow5tools get -t 1",
ids.len() as u64,
st_bytes,
st_dur,
"",
);
let ratio = st_dur.as_secs_f64() / rust_dur.as_secs_f64();
println!(" Ratio (st/rust): {ratio:.2}x");
}
}
}
#[cfg(feature = "rayon")]
fn run_write_par(args: &Args, ids: &[String]) {
println!(
"\nParallel write benchmark -- {} reads, zstd + SVB-ZD",
ids.len()
);
println!(" Pre-load records via indexed reads, then time compression + I/O only");
println!(" Sequential write uses the same pre-loaded records (isolates compression cost)");
separator();
let src_reader = Slow5IndexedReader::open(&args.file).expect("open indexed for write-par");
let records: Vec<Record> = ids
.iter()
.map(|id| {
src_reader
.get(id)
.unwrap_or_else(|e| panic!("write_par get {id}: {e}"))
})
.collect();
let tmp_seq = "/tmp/harness_write_par_seq.blow5";
let mut writer = Slow5Writer::create(tmp_seq, blow5_header()).expect("create seq output");
let t0 = Instant::now();
for rec in &records {
writer.write(rec).expect("write seq");
}
writer.finish().expect("finish seq");
let seq_dur = t0.elapsed();
let seq_bytes = fs::metadata(tmp_seq).map(|m| m.len()).unwrap_or(0);
row(
"Rust sequential ",
records.len() as u64,
seq_bytes,
seq_dur,
"",
);
let tmp_par = "/tmp/harness_write_par.blow5";
let mut writer = Slow5Writer::create(tmp_par, blow5_header()).expect("create par output");
let t0 = Instant::now();
writer.write_all_par(&records).expect("write_all_par");
writer.finish().expect("finish par");
let par_dur = t0.elapsed();
let par_bytes = fs::metadata(tmp_par).map(|m| m.len()).unwrap_or(0);
row(
"Rust parallel ",
records.len() as u64,
par_bytes,
par_dur,
"",
);
let ratio = seq_dur.as_secs_f64() / par_dur.as_secs_f64();
println!(" Speedup (seq/par): {ratio:.2}x");
}
#[cfg(feature = "rayon")]
fn run_write_pipeline(args: &Args, ids: &[String], have_c: bool) {
use rayon::prelude::*;
use slow5lib::writer::ParallelSlow5Writer;
let n_threads = std::thread::available_parallelism()
.map(|n| n.get())
.unwrap_or(4);
println!("\nPipeline write -- {} reads, zstd + SVB-ZD", ids.len());
println!(" Rayon workers: parallel indexed pread + compress across {n_threads} threads");
println!(" Background thread: sequential I/O, overlaps with compression");
separator();
let id_list = "/tmp/harness_write_ids.txt";
write_id_list(ids, id_list);
let tmp_rust = "/tmp/harness_write_pipeline_rust.blow5";
let src_reader = Slow5IndexedReader::open(&args.file).expect("open indexed for pipeline write");
let par_writer = ParallelSlow5Writer::create(tmp_rust, blow5_header(), n_threads * 4)
.expect("create pipeline writer");
let write_handle = par_writer.write_handle();
let t0 = Instant::now();
ids.par_iter().for_each(|id| {
let rec = src_reader
.get(id)
.unwrap_or_else(|e| panic!("pipeline get {id}: {e}"));
write_handle.write(&rec).expect("pipeline write");
});
drop(write_handle);
par_writer.finish().expect("pipeline finish");
let rust_dur = t0.elapsed();
let rust_bytes = fs::metadata(tmp_rust).map(|m| m.len()).unwrap_or(0);
row("Rust pipeline", ids.len() as u64, rust_bytes, rust_dur, "");
if have_c {
let n_str = n_threads.to_string();
let tmp_st = "/tmp/harness_write_pipeline_st.blow5";
let t0 = Instant::now();
let status = Command::new(&args.slow5tools)
.args([
"get", "--to", "blow5", "-c", "zstd", "-s", "svb-zd", "-t", &n_str, "-l", id_list,
"-o", tmp_st,
])
.arg(&args.file)
.stderr(std::process::Stdio::null())
.status()
.expect("spawn slow5tools get for pipeline bench");
let st_dur = t0.elapsed();
if status.success() {
let st_bytes = fs::metadata(tmp_st).map(|m| m.len()).unwrap_or(0);
row(
&format!("slow5tools get -t {n_threads}"),
ids.len() as u64,
st_bytes,
st_dur,
"",
);
let ratio = st_dur.as_secs_f64() / rust_dur.as_secs_f64();
println!(" Ratio (st/rust): {ratio:.2}x");
}
}
}
fn blow5_header() -> Header {
Header {
version: (0, 2, 0),
num_read_groups: 1,
record_compression: RecordCompression::Zstd,
signal_compression: SignalCompression::SvbZd,
read_groups: vec![ReadGroup::default()],
aux_meta: AuxMeta::default(),
}
}
fn main() {
let args = parse_args();
let file_bytes = fs::metadata(&args.file).expect("stat file").len();
let gb = file_bytes as f64 / 1e9;
let have_c = !args.skip_c && slow5tools_ok(&args.slow5tools);
println!("slow5lib harness");
println!("================");
println!("File: {} ({:.1} GB)", args.file.display(), gb);
println!(
"slow5tools: {}",
if have_c {
args.slow5tools.as_str()
} else {
"not available (C comparisons skipped)"
}
);
print!("\nLoading index... ");
std::io::stdout().flush().ok();
let t0 = Instant::now();
let reader = Slow5IndexedReader::open(&args.file).expect("open indexed reader");
let all_ids: Vec<String> = reader.read_ids().map(str::to_string).collect();
println!(
"{} reads ({:.2}s)",
all_ids.len(),
t0.elapsed().as_secs_f64()
);
let correct_ids = strided_ids(&all_ids, args.n_correct);
let indexed_ids = strided_ids(&all_ids, args.n_indexed);
let write_ids = strided_ids(&all_ids, args.n_write);
if !args.skip_correct && have_c {
run_correctness(&args, &correct_ids);
} else if !args.skip_correct {
println!("\nCorrectness check: skipped (slow5tools not available)");
}
if !args.skip_seq {
run_seq(&args, file_bytes, have_c);
}
#[cfg(feature = "rayon")]
if !args.skip_par {
run_par(&args, file_bytes);
}
#[cfg(not(feature = "rayon"))]
if !args.skip_par {
println!("\nParallel read: compile with --features rayon to enable");
}
if !args.skip_idx {
run_idx(&args, &indexed_ids, have_c);
}
if !args.skip_write {
run_write(&args, &write_ids, have_c);
}
#[cfg(feature = "rayon")]
if !args.skip_write_par {
run_write_par(&args, &write_ids);
}
#[cfg(not(feature = "rayon"))]
if !args.skip_write_par {
println!("\nParallel write: compile with --features rayon to enable");
}
#[cfg(feature = "rayon")]
if !args.skip_write_pipe {
run_write_pipeline(&args, &write_ids, have_c);
}
#[cfg(not(feature = "rayon"))]
if !args.skip_write_pipe {
println!("\nPipeline write: compile with --features rayon to enable");
}
println!();
}