clawhdf5/crates/clawhdf5/tests/vds_interop.rs

//! Virtual Dataset (VDS) reads checked against libhdf5 (through h5py).
//!
//! Each test has h5py build virtual datasets and their source files in a temp
//! directory, record what libhdf5 reads back (shape and values) next to them,
//! and then compares that with what clawhdf5 reads from the same files.
//!
//! Skipped when python3 or h5py is unavailable, unless
//! `CLAWHDF5_REQUIRE_INTEROP=1`.

use std::path::Path;
use std::process::Command;

use clawhdf5::File;

/// The Python interpreter to drive interop checks with (`CLAWHDF5_PYTHON`
/// lets these run against a virtualenv holding h5py).
fn python() -> String {
    std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string())
}

/// When `CLAWHDF5_REQUIRE_INTEROP=1` (set in CI), a missing Python dependency
/// is a test failure instead of a silent skip.
fn interop_required() -> bool {
    std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1")
}

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

macro_rules! skip_if_no_python {
    () => {
        if !python_available() {
            assert!(
                !interop_required(),
                "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with h5py is not available"
            );
            eprintln!("SKIP: python3 with h5py not available");
            return;
        }
    };
}

/// Prelude for every generator script: `expect(file, dset, tag)` records what
/// libhdf5 reads for `file:dset` as `<tag>.expect` (shape line, values line).
const PRELUDE: &str = r#"
import h5py, numpy as np
def expect(fn, dset, tag):
    with h5py.File(fn, "r") as f:
        d = f[dset]
        a = d[...]
        with open(tag + ".expect", "w") as out:
            out.write(" ".join(str(n) for n in d.shape) + "\n")
            out.write(" ".join(repr(float(v)) for v in a.ravel()) + "\n")
"#;

/// Run `body` (after [`PRELUDE`]) with `dir` as the working directory, so
/// relative source file names land next to the virtual file.
fn generate(dir: &Path, body: &str) {
    let script = format!("{PRELUDE}\n{body}");
    let out = Command::new(python())
        .args(["-c", &script])
        .current_dir(dir)
        .output()
        .expect("failed to run python");
    assert!(
        out.status.success(),
        "generator failed:\nSTDOUT: {}\nSTDERR: {}",
        String::from_utf8_lossy(&out.stdout),
        String::from_utf8_lossy(&out.stderr)
    );
}

/// What libhdf5 read for `tag`: (shape, values as f64).
fn expected(dir: &Path, tag: &str) -> (Vec<u64>, Vec<f64>) {
    let text = std::fs::read_to_string(dir.join(format!("{tag}.expect"))).unwrap();
    let mut lines = text.lines();
    let parse_line = |l: Option<&str>| -> Vec<String> {
        l.unwrap_or("")
            .split_whitespace()
            .map(str::to_string)
            .collect()
    };
    let shape = parse_line(lines.next())
        .iter()
        .map(|s| s.parse().unwrap())
        .collect();
    let values = parse_line(lines.next())
        .iter()
        .map(|s| s.parse().unwrap())
        .collect();
    (shape, values)
}

/// Assert clawhdf5 reads `file:dset` exactly as libhdf5 did for `tag`.
fn assert_matches_libhdf5(dir: &Path, file: &str, dset: &str, tag: &str) {
    let (shape, values) = expected(dir, tag);
    let f = File::open(dir.join(file)).unwrap();
    let ds = f.dataset(dset).unwrap();
    assert_eq!(
        ds.shape().unwrap(),
        shape,
        "{tag}: shape differs from libhdf5"
    );
    let got = ds
        .read_f64()
        .unwrap_or_else(|e| panic!("{tag}: read failed: {e}"));
    assert_eq!(got.len(), values.len(), "{tag}: element count differs");
    for (i, (g, e)) in got.iter().zip(&values).enumerate() {
        assert!(
            g == e || (g.is_nan() && e.is_nan()),
            "{tag}: element {i} is {g}, libhdf5 reads {e}\n ours: {got:?}\n libhdf5: {values:?}"
        );
    }
}

// ---------------------------------------------------------------------------
// Selection encodings
// ---------------------------------------------------------------------------

/// Files written with the default (earliest) format bounds serialize every
/// VDS hyperslab as a version-1 *irregular* selection (4-byte block corners),
/// and a strided selection as many blocks. These were refused outright.
#[test]
fn vds_version1_irregular_hyperslab_selections() {
    skip_if_no_python!();
    let dir = tempfile::tempdir().unwrap();
    generate(
        dir.path(),
        r#"
with h5py.File("src.h5", "w") as s:
    s.create_dataset("a", data=np.arange(12.0))
    s.create_dataset("m", data=np.arange(20.0).reshape(4, 5))
with h5py.File("v1.h5", "w") as f:   # default libver: hyperslab version 1
    f.create_dataset("local", data=np.arange(10.0) * -1)
    lay = h5py.VirtualLayout(shape=(12,), dtype="f8")
    lay[0:4] = h5py.VirtualSource(".", "local", shape=(10,))[2:6]
    lay[4:10] = h5py.VirtualSource("src.h5", "a", shape=(12,))[::2]
    lay[10:12] = h5py.VirtualSource("src.h5", "a", shape=(12,))[10:12]
    f.create_virtual_dataset("strided", lay)
    lay = h5py.VirtualLayout(shape=(4, 6), dtype="f8")
    lay[:, 0:2] = h5py.VirtualSource("src.h5", "m", shape=(4, 5))[:, 3:5]
    lay[:, 2:6:2] = h5py.VirtualSource("src.h5", "m", shape=(4, 5))[:, 0:2]
    lay[:, 3:6:2] = h5py.VirtualSource("src.h5", "m", shape=(4, 5))[:, 2:4]
    f.create_virtual_dataset("grid", lay)
expect("v1.h5", "strided", "strided")
expect("v1.h5", "grid", "grid")
"#,
    );
    assert_matches_libhdf5(dir.path(), "v1.h5", "strided", "strided");
    assert_matches_libhdf5(dir.path(), "v1.h5", "grid", "grid");
}

// ---------------------------------------------------------------------------
// Mapping list encoding
// ---------------------------------------------------------------------------

/// With a 2.0 low version bound libhdf5 writes the mapping list as block
/// version 1: a flags byte per entry, and repeated names stored as the index
/// of the entry that first spelled them out. The flags byte was mistaken for
/// an empty (same-file) name.
#[test]
fn vds_mapping_block_version1_shared_names() {
    skip_if_no_python!();
    let dir = tempfile::tempdir().unwrap();
    generate(
        dir.path(),
        r#"
name = "a_rather_long_dataset_name"
with h5py.File("a_rather_long_source_file.h5", "w") as s:
    s.create_dataset(name, data=np.arange(12.0) + 100)
with h5py.File("shared.h5", "w", libver=("v200", "v200")) as f:
    f.create_dataset(name, data=np.arange(12.0) * -1)
    lay = h5py.VirtualLayout(shape=(4, 4), dtype="f8")
    for i in range(3):
        src = h5py.VirtualSource("a_rather_long_source_file.h5", name, shape=(12,))
        lay[i] = src[4 * i:4 * i + 4]
    lay[3] = h5py.VirtualSource(".", name, shape=(12,))[0:4]
    f.create_virtual_dataset("v", lay)
# the heap block must really be version 1 for this test to mean anything
raw = open("shared.h5", "rb").read()
gcol = raw.index(b"GCOL")
assert raw[gcol + 32] == 1, "expected a version-1 VDS mapping block"
expect("shared.h5", "v", "shared")
"#,
    );
    assert_matches_libhdf5(dir.path(), "shared.h5", "v", "shared");
}

// ---------------------------------------------------------------------------
// Fill value
// ---------------------------------------------------------------------------

/// Elements no mapping supplies read as the virtual dataset's fill value, not
/// as 0: unmapped regions, a missing source file, a missing source dataset.
/// A source's own unallocated chunks read as *its* fill value.
#[test]
fn vds_unmapped_regions_read_as_fill_value() {
    skip_if_no_python!();
    let dir = tempfile::tempdir().unwrap();
    generate(
        dir.path(),
        r#"
for i in range(3):
    with h5py.File(f"src_{i}.h5", "w") as s:
        s.create_dataset("data", data=np.arange(10.0) + i * 100)
with h5py.File("sparse_src.h5", "w") as s:
    d = s.create_dataset("data", shape=(10,), chunks=(5,), dtype="f8", fillvalue=42.0)
    d[0:5] = np.arange(5.0) + 1000   # the second chunk is never written
for libver in ["earliest", "latest"]:
    with h5py.File(f"fill_{libver}.h5", "w", libver=libver) as f:
        f.create_dataset("local", data=np.arange(10.0) * -1)
        lay = h5py.VirtualLayout(shape=(6, 10), dtype="f8")
        for i in range(3):
            lay[i] = h5py.VirtualSource(f"src_{i}.h5", "data", shape=(10,))
        lay[3] = h5py.VirtualSource("no_such_file.h5", "data", shape=(10,))
        lay[4] = h5py.VirtualSource("src_0.h5", "no_such_dataset", shape=(10,))
        # row 5 is not mapped at all
        f.create_virtual_dataset("files", lay, fillvalue=-1.0)
        lay = h5py.VirtualLayout(shape=(20,), dtype="f8")
        lay[0:10] = h5py.VirtualSource(".", "local", shape=(10,))
        f.create_virtual_dataset("same_file", lay, fillvalue=7.0)
        lay = h5py.VirtualLayout(shape=(12,), dtype="f8")
        lay[1:11] = h5py.VirtualSource("sparse_src.h5", "data", shape=(10,))
        f.create_virtual_dataset("sparse_source", lay, fillvalue=-3.5)
        lay = h5py.VirtualLayout(shape=(3, 4), dtype="i4")
        lay[1, :] = h5py.VirtualSource(".", "ints", shape=(4,))
        f.create_dataset("ints", data=np.arange(4, dtype="i4") + 1)
        f.create_virtual_dataset("int_fill", lay, fillvalue=-99)
    for name in ["files", "same_file", "sparse_source", "int_fill"]:
        expect(f"fill_{libver}.h5", name, f"{name}_{libver}")
"#,
    );
    for libver in ["earliest", "latest"] {
        let file = format!("fill_{libver}.h5");
        for name in ["files", "same_file", "sparse_source", "int_fill"] {
            assert_matches_libhdf5(dir.path(), &file, name, &format!("{name}_{libver}"));
        }
    }

    // A selection read goes through the same fill-aware assembly.
    let f = File::open(dir.path().join("fill_latest.h5")).unwrap();
    let sel = clawhdf5::Selection::slice(std::slice::from_ref(&(8..14)));
    let got = f
        .dataset("same_file")
        .unwrap()
        .read_f64_selection(&sel)
        .unwrap();
    assert_eq!(got, vec![-8.0, -9.0, 7.0, 7.0, 7.0, 7.0]);
}

/// A source name that would leave the virtual file's directory is refused
/// with an error; it used to be skipped and read silently as fill.
#[test]
fn vds_source_outside_directory_is_an_error_not_fill() {
    skip_if_no_python!();
    let dir = tempfile::tempdir().unwrap();
    std::fs::create_dir(dir.path().join("sub")).unwrap();
    generate(
        dir.path(),
        r#"
with h5py.File("src.h5", "w") as s:
    s.create_dataset("data", data=np.arange(4.0))
with h5py.File("sub/up.h5", "w", libver="latest") as f:
    lay = h5py.VirtualLayout(shape=(4,), dtype="f8")
    lay[:] = h5py.VirtualSource("../src.h5", "data", shape=(4,))
    f.create_virtual_dataset("v", lay, fillvalue=-1.0)
with h5py.File("nested.h5", "w", libver="latest") as f:
    lay = h5py.VirtualLayout(shape=(4,), dtype="f8")
    lay[:] = h5py.VirtualSource("sub/inner.h5", "data", shape=(4,))
    f.create_virtual_dataset("v", lay, fillvalue=-1.0)
with h5py.File("sub/inner.h5", "w") as s:
    s.create_dataset("data", data=np.arange(4.0) + 10)
expect("nested.h5", "v", "nested")
"#,
    );
    // libhdf5 resolves "../src.h5" (and would read [0, 1, 2, 3]); we refuse
    // to leave the directory, and say so.
    let f = File::open(dir.path().join("sub/up.h5")).unwrap();
    let err = f.dataset("v").unwrap().read_f64().unwrap_err();
    assert!(
        err.to_string().contains("not followed"),
        "unexpected error: {err}"
    );
    // A relative name below the virtual file's directory resolves there.
    assert_matches_libhdf5(dir.path(), "nested.h5", "v", "nested");
}

/// Variable-length and reference elements are addresses into their own file.
/// Copied raw from an external source they would be decoded against the
/// virtual dataset's file and name another object, so they are refused.
#[test]
fn vds_external_variable_length_source_is_an_error_not_foreign_addresses() {
    skip_if_no_python!();
    let dir = tempfile::tempdir().unwrap();
    generate(
        dir.path(),
        r#"
st = h5py.string_dtype()
with h5py.File("src.h5", "w") as s:
    s.create_dataset("names", data=np.array(["alpha", "beta", "gamma"], dtype=object), dtype=st)
    s.create_dataset("refs", data=[s.ref, s.ref], dtype=h5py.ref_dtype)
with h5py.File("v.h5", "w", libver="latest") as f:
    f.create_dataset("pad", data=np.arange(64.0))
    lay = h5py.VirtualLayout(shape=(3,), dtype=st)
    lay[:] = h5py.VirtualSource("src.h5", "names", shape=(3,))
    f.create_virtual_dataset("names", lay)
    lay = h5py.VirtualLayout(shape=(2,), dtype=h5py.ref_dtype)
    lay[:] = h5py.VirtualSource("src.h5", "refs", shape=(2,))
    f.create_virtual_dataset("refs", lay)
"#,
    );
    let f = File::open(dir.path().join("v.h5")).unwrap();
    for name in ["names", "refs"] {
        let err = f
            .dataset(name)
            .unwrap()
            .read_selection(&clawhdf5_format::selection::Selection::All)
            .expect_err("raw addresses from another file must not be returned");
        assert!(
            err.to_string().contains("from another file"),
            "{name}: unexpected error: {err}"
        );
    }
}

// ---------------------------------------------------------------------------
// Unlimited and printf-style mappings
// ---------------------------------------------------------------------------

/// Helpers for building unlimited VDS mappings through h5py's low-level API.
const UNLIMITED_HELPERS: &str = r#"
U = h5py.h5s.UNLIMITED
def space(dims, maxdims, start=None, count=None, stride=None, block=None):
    s = h5py.h5s.create_simple(dims, maxdims)
    if start is not None:
        s.select_hyperslab(start, count, stride, block)
    return s
def make_vds(fn, name, dims, maxdims, maps, fill, libver="latest", mode="w"):
    # maps: [(vsel_kwargs, source_file, source_dataset, source_space)]
    with h5py.File(fn, mode, libver=libver) as f:
        dcpl = h5py.h5p.create(h5py.h5p.DATASET_CREATE)
        for vsel, sfile, sdset, sspace in maps:
            dcpl.set_virtual(space(dims, maxdims, **vsel), sfile.encode(), sdset.encode(), sspace)
        dcpl.set_fill_value(np.array(fill, dtype="f8"))
        h5py.h5d.create(f.id, name.encode(), h5py.h5t.IEEE_F64LE,
                        h5py.h5s.create_simple(dims, maxdims), dcpl=dcpl)
"#;

/// printf-style names: block `j` of the virtual selection comes from the
/// source named with `j` in place of `%b` (`%%` is a literal `%`), probing
/// j = 0, 1, ... until the first missing source. libhdf5 also recomputes the
/// extent from what it finds, so the stored dataspace is not the shape.
#[test]
fn vds_printf_source_names() {
    skip_if_no_python!();
    let dir = tempfile::tempdir().unwrap();
    let body = format!(
        "{UNLIMITED_HELPERS}{}",
        r#"
for i in [0, 1, 2, 4]:   # 3 is missing: 4 is past the first gap and unused
    with h5py.File(f"vds_src_{i}.h5", "w") as s:
        s.create_dataset("data", data=np.arange(10.0) + i * 100)
    with h5py.File(f"p%c_{i}.h5", "w") as s:
        s.create_dataset("data", data=np.arange(10.0) - i * 100)
for libver in ["earliest", "latest"]:
    fn = f"printf_{libver}.h5"
    make_vds(fn, "files", (10,), (U,),
             [(dict(start=(0,), count=(U,), stride=(10,), block=(10,)), "vds_src_%b.h5", "data",
               space((10,), (10,), (0,), (1,), (1,), (10,)))], -1.0, libver)
    # interleaved blocks with gaps between them, and an escaped percent sign
    make_vds(fn, "escaped", (4,), (U,),
             [(dict(start=(1,), count=(U,), stride=(6,), block=(4,)), "p%%c_%b.h5", "data",
               space((10,), (10,), (2,), (1,), (1,), (4,)))], -5.0, libver, "a")
    # printf in the dataset name, same file, 2-D frames
    with h5py.File(fn, "a") as f:
        for i in range(3):
            f.create_dataset(f"frame_{i}", data=np.arange(6.0).reshape(2, 3) + 10 * i)
    with h5py.File(fn, "a", libver=libver) as f:
        dcpl = h5py.h5p.create(h5py.h5p.DATASET_CREATE)
        dcpl.set_virtual(space((1, 2, 3), (U, 2, 3), (0, 0, 0), (U, 1, 1), (1, 1, 1), (1, 2, 3)),
                         b".", b"frame_%b", space((2, 3), (2, 3)))
        h5py.h5d.create(f.id, b"frames", h5py.h5t.IEEE_F64LE,
                        h5py.h5s.create_simple((1, 2, 3), (U, 2, 3)), dcpl=dcpl)
    for name in ["files", "escaped", "frames"]:
        expect(fn, name, f"{name}_{libver}")
"#
    );
    generate(dir.path(), &body);
    for libver in ["earliest", "latest"] {
        for name in ["files", "escaped", "frames"] {
            let file = format!("printf_{libver}.h5");
            assert_matches_libhdf5(dir.path(), &file, name, &format!("{name}_{libver}"));
        }
    }
    // The shape libhdf5 reports: three 10-element blocks.
    let f = File::open(dir.path().join("printf_latest.h5")).unwrap();
    assert_eq!(f.dataset("files").unwrap().shape().unwrap(), vec![30]);
}

/// Unlimited source and virtual selections: each mapping covers as much as
/// its source's current extent fills (a partial last block included), the
/// extent is the largest of them but never smaller than the limited
/// mappings need, and a missing source contributes nothing.
#[test]
fn vds_unlimited_mappings_follow_source_extents() {
    skip_if_no_python!();
    let dir = tempfile::tempdir().unwrap();
    let body = format!(
        "{UNLIMITED_HELPERS}{}",
        r#"
with h5py.File("grow.h5", "w") as s:
    s.create_dataset("a", data=np.arange(7.0) + 1, maxshape=(None,))
    s.create_dataset("b", data=np.arange(5.0) + 100, maxshape=(None,))
    s.create_dataset("rows", data=np.arange(12.0).reshape(4, 3) + 50, maxshape=(None, 3))
unlim_src = lambda: space((1,), (U,), (0,), (U,), (1,), (1,))
for libver in ["earliest", "latest"]:
    fn = f"unlim_{libver}.h5"
    make_vds(fn, "interleaved", (1,), (U,), [
        # blocks of 3 every 4: 7 source elements end mid-block
        (dict(start=(0,), count=(U,), stride=(4,), block=(3,)), "grow.h5", "a", unlim_src()),
        (dict(start=(3,), count=(U,), stride=(4,), block=(1,)), "grow.h5", "b", unlim_src()),
        (dict(start=(0,), count=(U,), stride=(1,), block=(1,)), "missing.h5", "a", unlim_src()),
    ], -2.0, libver)
    make_vds(fn, "rows", (6, 3), (U, 3), [
        # an unlimited *block*, plus a limited mapping reaching row 5
        (dict(start=(0, 0), count=(1, 1), stride=(1, 1), block=(U, 3)), "grow.h5", "rows",
         space((1, 3), (U, 3), (0, 0), (1, 1), (1, 1), (U, 3))),
        (dict(start=(5, 0), count=(1, 1), stride=(1, 1), block=(1, 3)), "grow.h5", "rows",
         space((4, 3), (U, 3), (1, 0), (1, 1), (1, 1), (1, 3))),
    ], -4.0, libver, "a")
    for name in ["interleaved", "rows"]:
        expect(fn, name, f"{name}_{libver}")
"#
    );
    generate(dir.path(), &body);
    for libver in ["earliest", "latest"] {
        for name in ["interleaved", "rows"] {
            let file = format!("unlim_{libver}.h5");
            assert_matches_libhdf5(dir.path(), &file, name, &format!("{name}_{libver}"));
        }
    }
    // "a" (7 elements, blocks of 3 every 4) ends at 9; "b" (5 elements from
    // 3, every 4) at 20. "rows" fills 4 rows but a limited mapping needs 6.
    let f = File::open(dir.path().join("unlim_latest.h5")).unwrap();
    assert_eq!(f.dataset("interleaved").unwrap().shape().unwrap(), vec![20]);
    assert_eq!(f.dataset("rows").unwrap().shape().unwrap(), vec![6, 3]);
}

/// libhdf5's own VDS test files (HDF5 `tools/test/testfiles/vds`):
/// printf-style Eiger frames (with a source past the first gap that must be
/// ignored), a printf mapping in the 1.10 format, and Percival's four
/// interleaved unlimited sources of different lengths.
#[test]
fn vds_libhdf5_test_files() {
    skip_if_no_python!();
    let dir = tempfile::tempdir().unwrap();
    let fixtures = Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures/vds");
    for entry in std::fs::read_dir(&fixtures).unwrap() {
        let path = entry.unwrap().path();
        if path.extension().is_some_and(|e| e == "h5") {
            std::fs::copy(&path, dir.path().join(path.file_name().unwrap())).unwrap();
        }
    }
    let cases = [
        ("vds-eiger.h5", "/VDS-Eiger"),
        ("4_vds.h5", "/vds_dset"),
        ("vds-percival-unlim-maxmin.h5", "/VDS-Percival-unlim-maxmin"),
    ];
    let mut body = String::new();
    for (i, (file, dset)) in cases.iter().enumerate() {
        body.push_str(&format!("expect({file:?}, {dset:?}, \"case{i}\")\n"));
    }
    generate(dir.path(), &body);
    for (i, (file, dset)) in cases.iter().enumerate() {
        assert_matches_libhdf5(dir.path(), file, dset, &format!("case{i}"));
    }
    // Stored as 20 frames; only f-0.h5 is found before the first gap.
    let f = File::open(dir.path().join("vds-eiger.h5")).unwrap();
    assert_eq!(
        f.dataset("VDS-Eiger").unwrap().shape().unwrap(),
        vec![5, 10, 10]
    );
}

/// A source file with a metadata cache image (libhdf5's
/// `h5clear_mdc_image.h5`, whose root group's header exists only in the
/// image) is read through the image, as libhdf5 opens it. Source files were
/// read from their own bytes, and this one failed on the root group.
#[test]
fn vds_source_file_with_a_metadata_cache_image() {
    skip_if_no_python!();
    let dir = tempfile::tempdir().unwrap();
    let fixture = Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures/h5clear_mdc_image.h5");
    std::fs::copy(&fixture, dir.path().join("src.h5")).unwrap();
    generate(
        dir.path(),
        r#"
with h5py.File("vds.h5", "w", libver="latest") as f:
    lay = h5py.VirtualLayout(shape=(3, 100), dtype="i4")
    lay[0:3, :] = h5py.VirtualSource("src.h5", "DSET", shape=(50, 100))[5:8, :]
    f.create_virtual_dataset("v", lay)
expect("vds.h5", "v", "image_source")
"#,
    );
    assert_matches_libhdf5(dir.path(), "vds.h5", "v", "image_source");
}