crates/clawhdf5/Cargo.toml

@@ -48,8 +48,10 @@ bitshuffle = ["clawhdf5-format/bitshuffle"]
 bzip2 = ["clawhdf5-format/bzip2"]
 blosc = ["clawhdf5-format/blosc"]
 blosc2 = ["clawhdf5-format/blosc2"]
+# ZFP (32013), read-only.
+zfp = ["clawhdf5-format/zfp"]
 # Every plugin filter.
-plugin-filters = ["lzf", "bitshuffle", "bzip2", "blosc", "blosc2"]
+plugin-filters = ["lzf", "bitshuffle", "bzip2", "blosc", "blosc2", "zfp"]
 # Dataset::verify_provenance() — recompute a dataset's SHA-256 and compare
 # against its stored _provenance_sha256 attribute. On by default, matching
 # clawhdf5-format's own default-on `provenance` feature.

crates/clawhdf5/tests/plugin_filters_interop.rs

@@ -495,16 +495,15 @@ fn lzf_written_by_clawhdf5_reads_in_h5py() {
     });
 }
 
-/// ZFP is not implemented, and Blosc2 is not in a build without the
-/// `blosc2` feature: reading them must be a clear error naming the filter,
-/// never data.
+/// Blosc2 and ZFP are not in a build without their features: reading them
+/// must be a clear error naming the filter, never data.
 #[test]
-fn unimplemented_filters_are_a_clear_error() {
+fn filters_left_out_of_the_build_are_a_clear_error() {
     if !have_python("h5py, hdf5plugin") {
         return;
     }
     let dir = tempfile::tempdir().unwrap();
-    let path = dir.path().join("unimplemented.h5");
+    let path = dir.path().join("left_out.h5");
     run_python(
         r#"
 import sys
@@ -517,7 +516,10 @@ with h5py.File(sys.argv[1], 'w') as f:
         &[path.to_str().unwrap()],
     );
     let file = File::open(&path).unwrap();
-    let mut missing = vec![("zfp", 32013u16, "ZFP")];
+    let mut missing = Vec::new();
+    if !cfg!(feature = "zfp") {
+        missing.push(("zfp", 32013u16, "ZFP"));
+    }
     if !cfg!(feature = "blosc2") {
         missing.push(("blosc2", 32026, "Blosc2"));
     }
@@ -526,7 +528,7 @@ with h5py.File(sys.argv[1], 'w') as f:
             .dataset(name)
             .unwrap()
             .read_selection(&Selection::All)
-            .expect_err("an unimplemented filter must not read");
+            .expect_err("a filter left out of the build must not read");
         let msg = err.to_string();
         assert!(
             msg.contains(&id.to_string()) && msg.contains(label),

crates/clawhdf5/tests/zfp_interop.rs

@@ -0,0 +1,277 @@
+//! ZFP (H5Z-ZFP, filter 32013) against libhdf5 + libzfp.
+//!
+//! h5py with hdf5plugin (H5Z-ZFP 1.1.1, zfp 1.0.1) writes datasets in every
+//! ZFP mode (fixed rate, precision and accuracy, reversible, expert) for
+//! each type ZFP supports (int32, int64, float, double), in 1 to 4
+//! dimensions, with chunks that are partial at the dataset's edges, blocks
+//! that are partial at the chunks' edges, and chunks with unit dimensions
+//! (a lower-dimensional ZFP field). Next to each it stores what h5py reads
+//! back, unfiltered, and clawhdf5 must read the ZFP dataset bit for bit
+//! equal to that: the decoder is deterministic, so lossy modes have one
+//! right answer.
+//!
+//! Skipped when python3 with h5py and hdf5plugin is unavailable, unless
+//! `CLAWHDF5_REQUIRE_INTEROP=1`.
+#![cfg(feature = "zfp")]
+
+use std::collections::BTreeMap;
+use std::process::Command;
+
+use clawhdf5::File;
+use clawhdf5_format::selection::Selection;
+
+fn python() -> String {
+    std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string())
+}
+
+fn have_python() -> bool {
+    let ok = Command::new(python())
+        .args(["-c", "import h5py, hdf5plugin"])
+        .output()
+        .map(|o| o.status.success())
+        .unwrap_or(false);
+    if ok {
+        return true;
+    }
+    assert!(
+        !std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1"),
+        "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with h5py and hdf5plugin is not available"
+    );
+    eprintln!("SKIP: python3 with h5py and hdf5plugin not available");
+    false
+}
+
+fn run_python(script: &str, args: &[&str]) -> String {
+    let output = Command::new(python())
+        .arg("-c")
+        .arg(script)
+        .args(args)
+        .output()
+        .expect("failed to run python");
+    assert!(
+        output.status.success(),
+        "Python script failed:\nSTDOUT: {}\nSTDERR: {}",
+        String::from_utf8_lossy(&output.stdout),
+        String::from_utf8_lossy(&output.stderr)
+    );
+    String::from_utf8_lossy(&output.stdout).trim().to_string()
+}
+
+/// Writes `f{i}` (ZFP) and `r{i}` (h5py's reading of `f{i}`, unfiltered)
+/// for every mode x dtype x shape x data kind H5Z-ZFP accepts; each `f{i}`
+/// has a `case` attribute. Prints the number of pairs.
+const GENERATE: &str = r#"
+import sys
+import numpy as np, h5py, hdf5plugin
+path = sys.argv[1]
+MODES = [
+    ('rate 2.5', dict(rate=2.5)),
+    ('rate 8', dict(rate=8)),
+    ('rate 16', dict(rate=16)),
+    ('rate 31', dict(rate=31)),
+    ('rate 64', dict(rate=64)),
+    ('precision 6', dict(precision=6)),
+    ('precision 20', dict(precision=20)),
+    ('precision 64', dict(precision=64)),
+    ('accuracy 0.5', dict(accuracy=0.5)),
+    ('accuracy 1e-3', dict(accuracy=1e-3)),
+    ('accuracy 1e-12', dict(accuracy=1e-12)),
+    ('reversible', dict(reversible=True)),
+    ('expert', dict(minbits=24, maxbits=600, maxprec=30, minexp=-20)),
+    # maxbits past ZFP_MAX_BITS: H5Z-ZFP stores it as fixed precision.
+    ('expert maxbits 20000', dict(minbits=1, maxbits=20000, maxprec=40, minexp=-1074)),
+    # A budget of 10 bits a block: one bit left after a float's exponent.
+    # (Doubles are left out: 12 bits of exponent and flag overrun it, and
+    # libzfp's encoder then writes past its buffer.)
+    ('expert maxbits 10', dict(minbits=5, maxbits=10, maxprec=64, minexp=-100)),
+    ('expert minbits 900', dict(minbits=900, maxbits=4000, maxprec=64, minexp=-60)),
+]
+DTYPES = ['<f4', '<f8', '<i4', '<i8']
+SHAPES = [
+    ((37,), (16,)),
+    ((13, 22), (5, 8)),
+    ((9, 7, 6), (4, 5, 6)),
+    ((5, 6, 7, 3), (3, 5, 6, 3)),
+    ((4, 1, 30), (2, 1, 30)),
+    ((3, 10, 1, 9), (2, 10, 1, 9)),
+    ((64, 64), (64, 64)),
+]
+KINDS = ['smooth', 'noise', 'wide', 'zeros', 'special']
+rng = np.random.default_rng(13)
+
+def data(dt, shape, kind):
+    n = int(np.prod(shape))
+    d = np.dtype(dt)
+    idx = np.arange(n, dtype=np.float64)
+    if d.kind == 'f':
+        if kind == 'smooth':
+            v = 100 * np.sin(idx / 5.0) + idx / 3.0
+        elif kind == 'noise':
+            v = rng.normal(size=n) * 10.0 ** rng.uniform(-3, 3, n)
+        elif kind == 'wide':
+            lo, hi = (-150, 120) if d.itemsize == 4 else (-1075, 1000)
+            v = rng.choice([-1.0, 1.0], n) * np.exp2(rng.integers(lo, hi, n).astype(np.float64))
+            v[rng.random(n) < 0.1] = 0.0
+        elif kind == 'zeros':
+            v = np.zeros(n)
+            v[: n // 3] = 0.0
+        else:
+            v = rng.normal(size=n) * 1e3
+            v[rng.random(n) < 0.05] = np.inf
+            v[rng.random(n) < 0.05] = -np.inf
+            v[rng.random(n) < 0.05] = np.nan
+            v[rng.random(n) < 0.05] = -0.0
+        with np.errstate(over='ignore'):
+            return v.astype(dt).reshape(shape)
+    info = np.iinfo(d)
+    if kind == 'smooth':
+        v = ((idx * 7) % 1000 - 500).astype(np.int64)
+    elif kind == 'noise':
+        v = rng.integers(info.min, info.max, n, dtype=np.int64, endpoint=True)
+    elif kind == 'wide':
+        bits = rng.integers(0, d.itemsize * 8 - 2, n)
+        v = rng.integers(-(2 ** 20), 2 ** 20, n) << np.minimum(bits, d.itemsize * 8 - 22)
+    elif kind == 'zeros':
+        v = np.zeros(n, dtype=np.int64)
+    else:
+        v = rng.choice([info.min, info.max, 0, -1, 1], n)
+    return v.astype(dt).reshape(shape)
+
+# Written first and read back after the file is closed: h5py reads a chunk
+# still in libhdf5's chunk cache without decoding it.
+i = 0
+with h5py.File(path, 'w') as f:
+    for label, kw in MODES:
+        for dt in DTYPES:
+            if label == 'expert maxbits 10' and dt == '<f8':
+                continue
+            for shape, chunks in SHAPES:
+                for kind in KINDS:
+                    v = data(dt, shape, kind)
+                    name = f'f{i}'
+                    try:
+                        ds = f.create_dataset(name, data=v, chunks=chunks, **hdf5plugin.Zfp(**kw))
+                    except Exception:
+                        continue
+                    ds.attrs['case'] = f'{label}|{dt}|{shape}|{chunks}|{kind}'
+                    i += 1
+with h5py.File(path, 'a') as f:
+    for k in range(i):
+        f.create_dataset(f'r{k}', data=f[f'f{k}'][()])
+print(i)
+"#;
+
+#[test]
+fn every_mode_reads_bit_exact() {
+    if !have_python() {
+        return;
+    }
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("zfp.h5");
+    let n: usize = run_python(GENERATE, &[path.to_str().unwrap()])
+        .parse()
+        .unwrap();
+    let file = File::open(&path).unwrap();
+    // mode -> dtypes read
+    let mut seen: BTreeMap<String, Vec<String>> = BTreeMap::new();
+    let mut failures = Vec::new();
+    for i in 0..n {
+        let ds = file.dataset(&format!("f{i}")).unwrap();
+        let case = match ds.attrs().unwrap().get("case") {
+            Some(clawhdf5::AttrValue::String(s)) => s.clone(),
+            other => panic!("f{i}: case attribute {other:?}"),
+        };
+        let want = file
+            .dataset(&format!("r{i}"))
+            .unwrap()
+            .read_selection(&Selection::All)
+            .unwrap();
+        match ds.read_selection(&Selection::All) {
+            Ok(got) if got == want => {
+                let mut parts = case.split('|');
+                let mode = parts.next().unwrap().to_string();
+                let dt = parts.next().unwrap().to_string();
+                let e = seen.entry(mode).or_default();
+                if !e.contains(&dt) {
+                    e.push(dt);
+                }
+            }
+            Ok(got) => {
+                let first = got.iter().zip(&want).position(|(a, b)| a != b);
+                failures.push(format!("f{i} {case}: differs from byte {first:?}"));
+            }
+            Err(e) => failures.push(format!("f{i} {case}: {e}")),
+        }
+    }
+    assert!(
+        failures.is_empty(),
+        "{} of {n} datasets:\n{}",
+        failures.len(),
+        failures.join("\n")
+    );
+    // Every mode was exercised on every type H5Z-ZFP accepts it for.
+    for (mode, dts) in &seen {
+        assert!(dts.len() >= 2, "{mode}: only {dts:?}");
+    }
+    assert_eq!(seen.len(), 16, "{:?}", seen.keys());
+    eprintln!("{n} ZFP datasets bit-exact; modes and types: {seen:?}");
+}
+
+/// A header written on a big-endian machine: H5Z-ZFP finds the magic only
+/// after byte-swapping the `cd_values`, and then byte-swaps the decoded
+/// values, since the dataset's datatype is big-endian there. The file is
+/// made by swapping the header words of a little-endian one in place, so
+/// the datatype stays little-endian and libhdf5 reads the values swapped;
+/// clawhdf5 must read the same bytes.
+#[test]
+fn big_endian_header_swaps_the_values() {
+    if !have_python() {
+        return;
+    }
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("zfp_be.h5");
+    let n: usize = run_python(
+        r#"
+import sys, struct
+import numpy as np, h5py, hdf5plugin
+path = sys.argv[1]
+cases = [('<f4', dict(rate=12)), ('<f8', dict(precision=30)), ('<i4', dict(reversible=True)),
+         ('<i8', dict(accuracy=4))]
+cds = []
+with h5py.File(path, 'w') as f:
+    for k, (dt, kw) in enumerate(cases):
+        v = (np.arange(150) * 13 % 97 - 40).reshape(10, 15).astype(dt)
+        ds = f.create_dataset(f'f{k}', data=v, chunks=(4, 8), **hdf5plugin.Zfp(**kw))
+        cds.append(ds.id.get_create_plist().get_filter(0)[2])
+raw = bytearray(open(path, 'rb').read())
+for cd in cds:
+    le = struct.pack(f'<{len(cd)}I', *cd)
+    swapped = le[:4] + struct.pack(f'>{len(cd) - 1}I', *cd[1:])
+    at = raw.find(le)
+    assert at > 0 and raw.find(le, at + 1) < 0, 'cd_values not found once'
+    raw[at:at + len(le)] = swapped
+open(path, 'wb').write(raw)
+with h5py.File(path, 'a') as f:
+    for k in range(len(cases)):
+        f.create_dataset(f'r{k}', data=f[f'f{k}'][()])
+        f.create_dataset(f'o{k}', data=f[f'f{k}'][()].byteswap())
+print(len(cases))
+"#,
+        &[path.to_str().unwrap()],
+    )
+    .parse()
+    .unwrap();
+    let file = File::open(&path).unwrap();
+    for k in 0..n {
+        let read = |name: String| {
+            file.dataset(&name)
+                .unwrap()
+                .read_selection(&Selection::All)
+                .unwrap_or_else(|e| panic!("{name}: {e}"))
+        };
+        let got = read(format!("f{k}"));
+        assert!(got == read(format!("r{k}")), "f{k}: differs from h5py");
+        // Byte-swapped back, the values are the right ones.
+        assert!(got != read(format!("o{k}")), "f{k}: not swapped");
+    }
+}

scripts/ci-test.sh

@@ -67,7 +67,7 @@ run_step "cargo clippy --all-targets" cargo clippy \
 
 # 3. Clippy over clawhdf5-format's optional features, which the default
 #    workspace build never compiles (szip is left out: it needs libaec).
-#    plugin-filters = bitshuffle, bzip2, blosc, blosc2 (and the default-on lzf).
+#    plugin-filters = bitshuffle, bzip2, blosc, blosc2, zfp (and the default-on lzf).
 run_step "cargo clippy (format feature matrix)" cargo clippy \
     -p clawhdf5-format \
     --all-targets \
@@ -78,7 +78,7 @@ run_step "cargo clippy (format feature matrix)" cargo clippy \
 # dependencies (bitshuffle and blosc share code).
 plugin_filters_alone() {
     local f
-    for f in bitshuffle bzip2 blosc blosc2; do
+    for f in bitshuffle bzip2 blosc blosc2 zfp; do
         echo "--- $f"
         cargo clippy -p clawhdf5-format --all-targets --features "$f" -- -D warnings || return 1
     done
@@ -208,9 +208,10 @@ if "$PYTHON" -c "import h5py" >/dev/null 2>&1 || [ "${CLAWHDF5_REQUIRE_INTEROP:-
     # libhdf5's registered plugins) compile and run too.
     run_step "h5py interop (format, ignored tests)" cargo test \
         -p clawhdf5-format --features lz4,zstd --test writer_h5py_tests -- --include-ignored
-    # LZF, bitshuffle, bzip2 and Blosc both ways against h5py + hdf5plugin.
+    # LZF, bitshuffle, bzip2 and Blosc both ways against h5py + hdf5plugin;
+    # Blosc2 and ZFP (read-only) against what h5py reads.
     run_step "h5py interop (plugin filters)" cargo test \
-        -p clawhdf5 --features plugin-filters --test plugin_filters_interop
+        -p clawhdf5 --features plugin-filters --test plugin_filters_interop --test zfp_interop
 else
     echo ""
     echo "==> [h5py interop] SKIPPED: no h5py in $PYTHON"