gwseq-io 0.2.0

Rust library for processing bigWig, bigBed, BAM and HiC files
Documentation
//! Time the read paths on a real file.
//!
//! `cargo run --release --example timeit -- <path> [parallel]`
//!
//! Not a criterion benchmark — one number per path on one file, for comparing
//! against another implementation or another machine.

use std::time::Instant;

use gwseq_io::bbi::{BbiReader, ProfileRequest, QuantifyRequest, ValuesRequest, Zoom};
use gwseq_io::genomic::Locs;

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let args: Vec<String> = std::env::args().skip(1).collect();
    let path = args.first().ok_or("usage: timeit <path> [parallel]")?;
    let parallel: i64 = args.get(1).map_or(Ok(-1), |a| a.parse())?;

    let reader = BbiReader::open(path, parallel, 1.0 / 3.0, None, None)?;
    let first = reader.chr_sizes().iter().next().ok_or("no chromosomes")?;
    let (name, size) = (first.id.clone(), first.size);

    // A thousand scattered 10 kb windows, which is the shape of a real request.
    let n = 1000usize;
    let step = (size - 20_000) / n as i64;
    let chr_ids = vec![name; n];
    let starts: Vec<i64> = (0..n as i64).map(|i| 10_000 + i * step).collect();
    let ends: Vec<i64> = starts.iter().map(|s| s + 10_000).collect();

    let run = |label: &str, f: &dyn Fn() -> Result<usize, gwseq_io::Error>| {
        // One warm pass so the block cache is not what is being timed.
        let _ = f();
        let started = Instant::now();
        let reps = 5;
        let mut count = 0;
        for _ in 0..reps {
            count = f().unwrap();
        }
        let each = started.elapsed().as_secs_f64() / reps as f64;
        println!("{label:28} {:8.1} ms   {count} values", each * 1000.0);
    };

    for (label, bin_size, zoom) in [
        ("read_values bin=1", 1.0, Zoom::Full),
        ("read_values bin=100", 100.0, Zoom::Full),
        ("read_values bin=100 zoom=auto", 100.0, Zoom::Auto),
    ] {
        run(label, &|| {
            let locs = Locs::spans(&chr_ids, &starts, &ends)?;
            Ok(reader
                .read_values(&ValuesRequest::new(locs).bin_size(bin_size).zoom(zoom))?
                .len())
        });
    }
    run("quantify", &|| {
        let locs = Locs::spans(&chr_ids, &starts, &ends)?;
        Ok(reader.quantify(&QuantifyRequest::new(locs))?.len())
    });
    run("profile bin=100", &|| {
        let locs = Locs::spans(&chr_ids, &starts, &ends)?;
        Ok(reader
            .profile(&ProfileRequest::new(locs).bin_size(100.0))?
            .len())
    });
    Ok(())
}