gwseq-io 0.2.0

//! Properties, over inputs `proptest` chooses.
//!
//! The unit tests check the cases someone thought of; these check the ones
//! nobody did. Each is a statement that has to hold for *every* input of a
//! shape, and proptest's job is to find the input where it does not and then
//! shrink it to the smallest one that still fails — which is the part worth
//! having, because a 4 000-value counterexample teaches nothing and the
//! three-value one it shrinks to teaches everything.
//!
//! The properties here are the ones with no easier way to check them:
//!
//! - **Round trip.** Whatever is written comes back. Given the writer picks
//!   between three section encodings by a cost model, and the reader has a
//!   different code path for each, this is where a mis-encoded section shows.
//! - **Composition.** Binning at `n` and then at `m` is binning at `n × m`;
//!   walking a chromosome in windows is reading it whole. These are the
//!   library agreeing with itself, and they fail when two code paths that
//!   should be one have drifted.
//! - **Refusal.** A request that cannot be served is refused, not answered
//!   wrongly and not panicked on.
//!
//! Run more than the default 256 cases with `PROPTEST_CASES=5000`.

use proptest::prelude::*;

use gwseq_io::bbi::{
    BbiKind, BbiReader, BbiWriter, BbiWriterOptions, QuantifyRequest, SectionPolicy, ValuesRequest,
};
use gwseq_io::genomic::{BinMode, ChrMap, Locs, Reduce};

const CHR: &str = "chr1";

/// A temp file named after the property and the case, removed on drop.
struct Scratch(std::path::PathBuf);

impl Scratch {
    fn new(tag: &str) -> Self {
        use std::sync::atomic::{AtomicU64, Ordering};
        static N: AtomicU64 = AtomicU64::new(0);
        let mut path = std::env::temp_dir();
        path.push(format!(
            "gwseq-prop-{tag}-{}-{}.bigwig",
            std::process::id(),
            N.fetch_add(1, Ordering::Relaxed)
        ));
        Self(path)
    }
    fn as_str(&self) -> &str {
        self.0.to_str().expect("utf-8 temp path")
    }
}

impl Drop for Scratch {
    fn drop(&mut self) {
        let _ = std::fs::remove_file(&self.0);
    }
}

/// Values a bigWig can hold: finite, and spread over enough orders of magnitude
/// that a rounding difference in the wrong place is visible.
fn values(len: impl Into<proptest::collection::SizeRange>) -> impl Strategy<Value = Vec<f32>> {
    proptest::collection::vec(
        prop_oneof![
            3 => -1000.0f32..1000.0,
            1 => -1.0e-6f32..1.0e-6,
            1 => Just(0.0f32),
        ],
        len,
    )
}

fn write(path: &str, start: i64, span: i64, vals: &[f32], size: i64, parallel: i64) {
    write_with(
        path,
        start,
        span,
        vals,
        size,
        parallel,
        SectionPolicy::default(),
    );
}

fn write_with(
    path: &str,
    start: i64,
    span: i64,
    vals: &[f32],
    size: i64,
    parallel: i64,
    section_policy: SectionPolicy,
) {
    let mut w = BbiWriter::create(
        path,
        BbiWriterOptions {
            kind: BbiKind::BigWig,
            chr_sizes: Some(ChrMap::from_entries([(CHR.to_string(), size)])),
            parallel,
            section_policy,
            ..Default::default()
        },
    )
    .expect("a writable file");
    w.write_values(CHR, start, span, vals).expect("valid run");
    w.close().expect("a finished file");
}

fn read(path: &str, start: i64, end: i64, bin_size: f64) -> Vec<f32> {
    let reader = BbiReader::open(path, 1, 1.0 / 3.0, None, None).expect("a readable file");
    let locs = Locs::spans(&[CHR.to_string()], &[start], &[end]).expect("a well-formed request");
    reader
        .read_values(&ValuesRequest::new(locs).bin_size(bin_size))
        .expect("a valid request")
        .iter()
        .copied()
        .collect()
}

proptest! {
    /// A run written at span 1 reads back value for value, whatever the values
    /// were and however the writer chose to encode the sections holding them.
    #[test]
    fn a_written_run_reads_back_bit_for_bit(vals in values(1..400usize)) {
        let path = Scratch::new("roundtrip");
        let size = vals.len() as i64;
        write(path.as_str(), 0, 1, &vals, size, 1);
        let back = read(path.as_str(), 0, size, 1.0);
        prop_assert_eq!(back.len(), vals.len());
        for (i, (got, want)) in back.iter().zip(&vals).enumerate() {
            prop_assert_eq!(got.to_bits(), want.to_bits(), "value {}", i);
        }
    }

    /// The same, at a span the writer cannot express as fixedStep for free —
    /// every value covers `span` bases, so the section carries starts as well.
    #[test]
    fn a_run_of_any_span_reads_back_over_the_bases_it_covers(
        vals in values(1..120usize),
        span in 1i64..40,
    ) {
        let path = Scratch::new("span");
        let size = vals.len() as i64 * span;
        write(path.as_str(), 0, span, &vals, size, 1);
        let back = read(path.as_str(), 0, size, span as f64);
        prop_assert_eq!(back.len(), vals.len());
        for (i, (got, want)) in back.iter().zip(&vals).enumerate() {
            prop_assert_eq!(got.to_bits(), want.to_bits(), "value {} at span {}", i, span);
        }
    }

    /// The thread count is a performance knob, not a semantic one: the same
    /// input written by any number of workers is the same file.
    #[test]
    fn the_written_bytes_do_not_depend_on_the_thread_count(vals in values(200..1200usize)) {
        let a = Scratch::new("par1");
        let b = Scratch::new("par4");
        let size = vals.len() as i64;
        write(a.as_str(), 0, 1, &vals, size, 1);
        write(b.as_str(), 0, 1, &vals, size, 4);
        prop_assert_eq!(
            std::fs::read(&a.0).unwrap(),
            std::fs::read(&b.0).unwrap(),
            "the deflate pipeline changed the bytes"
        );
    }

    /// Summing a bin is summing its bases, so binning twice is binning once at
    /// the product. Checked with `sum` rather than `mean` because it is exact
    /// under regrouping in a way an average of averages is not.
    #[test]
    fn binning_twice_is_binning_once_at_the_product(
        vals in values(240..600usize),
        n in 2i64..8,
        m in 2i64..8,
    ) {
        let path = Scratch::new("compose");
        // A whole number of `n*m` bins, so no partial bin is involved and the
        // two groupings cover exactly the same bases.
        let bins = (vals.len() as i64) / (n * m);
        prop_assume!(bins >= 2);
        let end = bins * n * m;
        write(path.as_str(), 0, 1, &vals, vals.len() as i64, 1);

        let reader = BbiReader::open(path.as_str(), 1, 1.0 / 3.0, None, None).unwrap();
        let request = |bin: f64| {
            let locs = Locs::spans(&[CHR.to_string()], &[0], &[end]).unwrap();
            ValuesRequest::new(locs).bin_size(bin).bin_mode(BinMode::Sum)
        };
        let coarse = reader.read_values(&request((n * m) as f64)).unwrap();
        let fine = reader.read_values(&request(n as f64)).unwrap();

        prop_assert_eq!(coarse.len() as i64, bins);
        prop_assert_eq!(fine.len() as i64, bins * m);
        let fine: Vec<f32> = fine.iter().copied().collect();
        for (b, want) in coarse.iter().enumerate() {
            let got: f64 = fine[b * m as usize..(b + 1) * m as usize]
                .iter()
                .map(|v| *v as f64)
                .sum();

            // Regrouping an f32 sum is not exact, and the bound is not on the
            // result: summing k values in f32 carries an error of about
            // k * eps * sum|x|, so a bin whose values cancel has a small
            // answer and a large tolerance. Comparing against |want| instead
            // would fail on exactly those bins for no reason.
            let bin: &[f32] = &fine[b * m as usize..(b + 1) * m as usize];
            let magnitude: f64 = bin.iter().map(|v| v.abs() as f64).sum();
            let bound = (n * m) as f64 * f32::EPSILON as f64 * magnitude.max(1.0);
            prop_assert!(
                (got - *want as f64).abs() <= bound,
                "bin {}: {} regrouped vs {} (bound {})", b, got, want, bound
            );
        }
    }

    /// A walk of a chromosome is a read of it, in windows. The windows tile the
    /// same grid, so concatenating them is the whole read exactly — not
    /// approximately, since no bin straddles a window boundary.
    #[test]
    fn a_walk_concatenates_to_the_whole_read(
        vals in values(500..2000usize),
        bin_size in 1i64..64,
        span in 200i64..900,
    ) {
        let path = Scratch::new("walk");
        let size = vals.len() as i64;
        prop_assume!(size / bin_size >= 2);
        write(path.as_str(), 0, 1, &vals, size, 1);
        let reader = BbiReader::open(path.as_str(), 1, 1.0 / 3.0, None, None).unwrap();

        let locs = Locs::spans(&[CHR.to_string()], &[0], &[size]).unwrap();
        let whole = reader
            .read_values(&ValuesRequest::new(locs).bin_size(bin_size as f64))
            .unwrap();
        let plan = ValuesRequest::new(Locs::chromosomes(vec![CHR.to_string()]))
            .bin_size(bin_size as f64);
        let walked: Vec<f32> = reader
            .iter_all_values(&plan, span)
            .unwrap()
            .flat_map(|w| w.unwrap().to_vec())
            .collect();

        prop_assert_eq!(
            walked.iter().map(|v| v.to_bits()).collect::<Vec<_>>(),
            whole.iter().map(|v| v.to_bits()).collect::<Vec<_>>()
        );
    }

    /// `quantify` reduces the same numbers `read_values` returns, so `sum` over
    /// one locus is the sum of its bins and `max` is their maximum. If these
    /// two ever disagree, one of the reductions is reading different data.
    #[test]
    fn quantify_reduces_what_read_values_returns(vals in values(100..800usize)) {
        let path = Scratch::new("quantify");
        let size = vals.len() as i64;
        write(path.as_str(), 0, 1, &vals, size, 1);
        let reader = BbiReader::open(path.as_str(), 1, 1.0 / 3.0, None, None).unwrap();
        let locs = || Locs::spans(&[CHR.to_string()], &[0], &[size]).unwrap();

        let per_base = reader
            .read_values(&ValuesRequest::new(locs()).bin_size(1.0))
            .unwrap();
        let got_max = reader
            .quantify(&QuantifyRequest::new(locs()).reduce(Reduce::Max))
            .unwrap();
        let want_max = per_base.iter().copied().fold(f32::NEG_INFINITY, f32::max);
        prop_assert_eq!(got_max[0].to_bits(), want_max.to_bits());

        let got_count = reader
            .quantify(&QuantifyRequest::new(locs()).reduce(Reduce::Count))
            .unwrap();
        prop_assert_eq!(got_count[0], size as f32);
    }

    /// A window past the end of a chromosome is answered with the default, not
    /// refused and not read out of bounds.
    #[test]
    fn a_window_past_the_end_is_all_default(
        vals in values(50..200usize),
        beyond in 1i64..5000,
    ) {
        let path = Scratch::new("past");
        let size = vals.len() as i64;
        write(path.as_str(), 0, 1, &vals, size, 1);
        let reader = BbiReader::open(path.as_str(), 1, 1.0 / 3.0, None, None).unwrap();
        let locs = Locs::spans(&[CHR.to_string()], &[size + beyond], &[size + beyond + 100])
            .unwrap();
        let out = reader
            .read_values(&ValuesRequest::new(locs).bin_size(1.0).def_value(-7.0))
            .unwrap();
        let out: Vec<f32> = out.iter().copied().collect();
        prop_assert!(out.iter().all(|v| *v == -7.0), "{:?}", &out[..4]);
    }

    /// The L1 norm is `Σ|x|`, which is the sum only for data that never goes
    /// negative — and the reduction is there for data that does.
    #[test]
    fn l1norm_is_the_sum_of_absolute_values(vals in values(50..500usize)) {
        let path = Scratch::new("l1");
        let size = vals.len() as i64;
        write(path.as_str(), 0, 1, &vals, size, 1);
        let reader = BbiReader::open(path.as_str(), 1, 1.0 / 3.0, None, None).unwrap();
        let locs = || Locs::spans(&[CHR.to_string()], &[0], &[size]).unwrap();
        let got = reader
            .quantify(&QuantifyRequest::new(locs()).reduce(Reduce::L1Norm))
            .unwrap()[0] as f64;
        let want: f64 = vals.iter().map(|v| v.abs() as f64).sum();
        let bound = vals.len() as f64 * f32::EPSILON as f64 * want.max(1.0);
        prop_assert!((got - want).abs() <= bound, "{} against {}", got, want);
    }

    /// A window is snapped to the genome-wide grid of `bin_size`, at both
    /// ends and by the same amount — which is what keeps bin `n` meaning the
    /// same offset for every locus in a request, and so what makes a profile
    /// mean anything.
    #[test]
    fn a_window_keeps_its_width_however_it_is_offset_from_the_grid(
        vals in values(2000..4000usize),
        bin_size in 2i64..97,
        offset in 0i64..500,
    ) {
        let path = Scratch::new("grid");
        let size = vals.len() as i64;
        write(path.as_str(), 0, 1, &vals, size, 1);
        let reader = BbiReader::open(path.as_str(), 1, 1.0 / 3.0, None, None).unwrap();

        let width = bin_size * 10;
        prop_assume!(offset + width < size);
        let locs = Locs::spans(&[CHR.to_string()], &[offset], &[offset + width]).unwrap();
        let out = reader
            .read_values(&ValuesRequest::new(locs).bin_size(bin_size as f64))
            .unwrap();
        // Ten bins whatever the offset: both ends move by `offset % bin_size`,
        // so the span between them is untouched.
        prop_assert_eq!(out.len(), 10);

        // And the values are the ones on the grid, not the ones from `offset`.
        let snapped = offset - offset.rem_euclid(bin_size);
        let want: Vec<f32> = (0..10)
            .map(|b| {
                let lo = (snapped + b * bin_size) as usize;
                let hi = lo + bin_size as usize;
                (vals[lo..hi].iter().map(|v| *v as f64).sum::<f64>() / bin_size as f64) as f32
            })
            .collect();
        for (b, (got, expect)) in out.iter().zip(&want).enumerate() {
            let scale = (expect.abs() as f64).max(1.0);
            prop_assert!(
                (*got as f64 - *expect as f64).abs() <= scale * 1e-5,
                "bin {}: {} against {}", b, got, expect
            );
        }
    }

    /// A fractional bin size is refused by every read path, not just by the
    /// whole-file walks.
    #[test]
    fn a_fractional_bin_size_is_refused(
        vals in values(50..200usize),
        whole in 1i64..50,
        frac in 1u32..1000,
    ) {
        let bin_size = whole as f64 + frac as f64 / 1000.0;
        prop_assume!(bin_size.fract() != 0.0);
        let path = Scratch::new("frac");
        let size = vals.len() as i64;
        write(path.as_str(), 0, 1, &vals, size, 1);
        let reader = BbiReader::open(path.as_str(), 1, 1.0 / 3.0, None, None).unwrap();
        let locs = Locs::spans(&[CHR.to_string()], &[0], &[size]).unwrap();
        let err = reader
            .read_values(&ValuesRequest::new(locs).bin_size(bin_size))
            .unwrap_err()
            .to_string();
        prop_assert!(err.contains("whole number of base pairs"), "{}", err);
    }

    /// The section-encoding policy decides how a file is *laid out* and
    /// nothing else. Every rule in `SectionPolicy` has to produce a file that
    /// reads back the same values — otherwise the size comparison in
    /// `examples/section_policy.rs` is between files that are not
    /// interchangeable, and the smallest one is worthless.
    ///
    /// This is what no size comparison can check and no ordinary correctness
    /// test would think to: every policy writes a *valid* bigWig, so a reader
    /// cannot tell them apart, and that is precisely why the values have to be
    /// asserted equal rather than assumed to be.
    #[test]
    fn every_section_policy_writes_the_same_values(
        vals in values(200..1200usize),
        span in 1i64..12,
        break_every in 3usize..40,
        break_span in 1i64..30,
    ) {
        // A run whose span changes every so often, which is what makes the
        // policies disagree about where a section ends.
        let mut items = Vec::new();
        let mut at = 0i64;
        for (i, v) in vals.iter().enumerate() {
            let s = if i % break_every == 0 { break_span } else { span };
            items.push((at, s, *v));
            at += s;
        }
        let size = at + 10;

        let mut baseline: Option<Vec<u32>> = None;
        for policy in [
            SectionPolicy::Cost,
            SectionPolicy::Room,
            SectionPolicy::Runt,
            SectionPolicy::Adaptive,
            SectionPolicy::Split,
            SectionPolicy::Widen,
        ] {
            let path = Scratch::new("policy");
            {
                let mut w = BbiWriter::create(
                    path.as_str(),
                    BbiWriterOptions {
                        kind: BbiKind::BigWig,
                        chr_sizes: Some(ChrMap::from_entries([(CHR.to_string(), size)])),
                        parallel: 1,
                        section_policy: policy,
                        ..Default::default()
                    },
                )
                .expect("a writable file");
                for (s, sp, v) in &items {
                    w.write_value(CHR, *s, *s + *sp, *v).expect("valid value");
                }
                w.close().expect("a finished file");
            }
            let reader = BbiReader::open(path.as_str(), 1, 1.0 / 3.0, None, None).unwrap();
            let locs = Locs::spans(&[CHR.to_string()], &[0], &[size]).unwrap();
            let bits: Vec<u32> = reader
                .read_values(&ValuesRequest::new(locs).bin_size(1.0))
                .unwrap()
                .iter()
                .map(|v| v.to_bits())
                .collect();
            match &baseline {
                None => baseline = Some(bits),
                Some(first) => prop_assert_eq!(&bits, first, "{:?} read back differently", policy),
            }
        }
    }

    /// The same for a run handed over whole, which takes the writer's fast
    /// path and the flush-before-a-run decision with it.
    #[test]
    fn a_run_reads_back_the_same_under_every_policy(
        vals in values(100..800usize),
        span in 1i64..15,
    ) {
        let size = vals.len() as i64 * span;
        let mut baseline: Option<Vec<u32>> = None;
        for policy in [
            SectionPolicy::Cost,
            SectionPolicy::Adaptive,
            SectionPolicy::Split,
            SectionPolicy::Widen,
        ] {
            let path = Scratch::new("runpolicy");
            write_with(path.as_str(), 0, span, &vals, size, 1, policy);
            let bits: Vec<u32> = read(path.as_str(), 0, size, span as f64)
                .iter()
                .map(|v| v.to_bits())
                .collect();
            match &baseline {
                None => baseline = Some(bits),
                Some(first) => prop_assert_eq!(&bits, first, "{:?} read back differently", policy),
            }
        }
    }

    /// An end before its start is refused when the request meets a file, at
    /// any coordinates, with an error rather than a panic or a wrong answer.
    ///
    /// Not at `Locs::spans`, which is a constructor and knows nothing about
    /// the file: a locus is only well-formed or not against the chromosome it
    /// names, so that is where the check belongs and where a caller meets it.
    #[test]
    fn an_inverted_window_is_refused(
        vals in values(20..60usize),
        start in 0i64..1000,
        back in 1i64..1000,
    ) {
        let path = Scratch::new("inverted");
        write(path.as_str(), 0, 1, &vals, 2000, 1);
        let reader = BbiReader::open(path.as_str(), 1, 1.0 / 3.0, None, None).unwrap();
        let locs = Locs::spans(&[CHR.to_string()], &[start], &[start - back]).unwrap();
        let result = reader.read_values(&ValuesRequest::new(locs).bin_size(1.0));
        prop_assert!(
            result.is_err(),
            "{}:{}-{} was answered, not refused", CHR, start, start - back
        );
    }
}