clawhdf5/crates/clawhdf5/tests/scaleoffset_interop.rs

//! Scale-offset datasets written by h5py (libhdf5) read exactly as h5py
//! reads them.
//!
//! h5py writes the whole matrix the filter has: every integer type (`i1` ..
//! `u8`) and `f4`/`f8`, both byte orders, with and without a fill value
//! (the type's minimum, maximum, or another value), with random, constant,
//! all-fill and extreme data, and every interesting `scaleoffset` setting
//! (0 = let libhdf5 choose the bits, a few bits, full width less one, full
//! width; decimal scale factors 0..7 for floats) — about 1480 datasets. For
//! each one h5py's decoded values are stored uncompressed next to it, and
//! the raw bytes clawhdf5 decodes must equal them.
//!
//! Until 2026-09-26 clawhdf5 silently returned wrong values for 332 of
//! these, in every release that decoded scale-offset: ordinary `u8`, `u4`,
//! `i8` and `f4` data with a wide range (where libhdf5 stores the elements
//! as they are, at full width), chunks whose `minval` field was recorded at
//! another size than 8 bytes, and chunks with `minbits` 0 and a fill value.
//!
//! Skipped when python3 with h5py is unavailable, unless
//! `CLAWHDF5_REQUIRE_INTEROP=1`.

use std::process::Command;

use clawhdf5::File;

fn python() -> String {
    std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string())
}

fn python_available() -> bool {
    Command::new(python())
        .args(["-c", "import h5py, numpy"])
        .output()
        .map(|o| o.status.success())
        .unwrap_or(false)
}

const GENERATE: &str = r#"
import sys, h5py, numpy as np
rng = np.random.default_rng(7)
out, expect = sys.argv[1], sys.argv[2]
made = []
with h5py.File(out, "w") as f:
    def make(name, a, so, fv, desc):
        try:
            d = f.create_dataset(name, data=a, chunks=(16,), scaleoffset=so, fillvalue=fv)
        except Exception:
            return  # a combination libhdf5 refuses to write
        d.attrs["desc"] = desc
        made.append(name)
    i = 0
    for dt in ["i1", "u1", "i2", "u2", "i4", "u4", "i8", "u8"]:
        for bo in "<>":
            t = np.dtype(bo + dt)
            info = np.iinfo(t)
            for so in [0, 1, 3, 8 * t.itemsize - 1, 8 * t.itemsize]:
                for fill in [None, "min", "max", "mid"]:
                    for pat in ["rand", "const", "allfill", "extreme"]:
                        n = 64
                        fv = {None: None, "min": info.min, "max": info.max, "mid": t.type(7)}[fill]
                        if pat == "rand":
                            lo = max(info.min, -50) if so else info.min
                            hi = min(info.max, 50) if so else info.max
                            a = rng.integers(lo, hi, size=n, endpoint=True,
                                             dtype=np.int64 if t.kind == "i" else np.uint64).astype(t)
                        elif pat == "const":
                            a = np.full(n, 3, t)
                        elif pat == "extreme":
                            a = np.array([info.min, info.max] * (n // 2), t)
                        else:
                            if fv is None:
                                continue
                            a = np.full(n, fv, t)
                        make(f"d{i}", a, so, fv, f"{bo}{dt} so={so} fill={fill} pat={pat}")
                        i += 1
    for dt in ["f4", "f8"]:
        for bo in "<>":
            t = np.dtype(bo + dt)
            for so in [0, 1, 2, 4, 7]:
                for fill in [None, -1.5, 0.0]:
                    for pat in ["rand", "const", "allfill", "neg", "big"]:
                        n = 50
                        if pat == "rand":
                            a = rng.normal(size=n).astype(t) * 10
                        elif pat == "const":
                            a = np.full(n, 2.25, t)
                        elif pat == "neg":
                            a = -np.abs(rng.normal(size=n)).astype(t) * 1000
                        elif pat == "big":
                            a = rng.normal(size=n).astype(t) * 1e6
                        else:
                            if fill is None:
                                continue
                            a = np.full(n, fill, t)
                        make(f"d{i}", a, so, fill, f"{bo}{dt} so={so} fill={fill} pat={pat}")
                        i += 1
# What libhdf5 decodes, stored uncompressed in the same datatype.
with h5py.File(out, "r") as f, h5py.File(expect, "w") as e:
    for name in made:
        e.create_dataset(name, data=f[name][()])
print(len(made))
"#;

#[test]
fn every_scale_offset_dataset_reads_as_h5py_reads_it() {
    if !python_available() {
        assert!(
            !std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1"),
            "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with h5py is not available"
        );
        eprintln!("SKIP: python3 with h5py not available");
        return;
    }
    let dir = tempfile::tempdir().unwrap();
    let written = dir.path().join("scaleoffset.h5");
    let expected = dir.path().join("expected.h5");
    let out = Command::new(python())
        .args(["-c", GENERATE])
        .arg(&written)
        .arg(&expected)
        .output()
        .expect("failed to run python");
    assert!(
        out.status.success(),
        "python failed:\n{}",
        String::from_utf8_lossy(&out.stderr)
    );
    let count: usize = String::from_utf8_lossy(&out.stdout).trim().parse().unwrap();
    assert!(count > 1400, "only {count} datasets written");

    let file = File::open(&written).unwrap();
    let reference = File::open(&expected).unwrap();
    let mut names = reference.root().datasets().unwrap();
    names.sort();
    assert_eq!(names.len(), count);
    let mut wrong = Vec::new();
    for name in &names {
        let want = reference.read_multi(&[name]).unwrap().remove(0);
        let got = file.read_multi(&[name]).map(|mut v| v.remove(0));
        if got.as_ref().ok() != Some(&want) {
            let desc = file.dataset(name).unwrap().attrs().unwrap().remove("desc");
            wrong.push(format!("{name} {desc:?}: {:?}", got.map(|g| g.len())));
        }
    }
    assert!(
        wrong.is_empty(),
        "{} of {count} scale-offset datasets differ from h5py:\n{}",
        wrong.len(),
        wrong.join("\n")
    );
}