From c6a7bbfc67930a8a05c8f14b1f8400f3f86627e6 Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 19 Sep 2026 13:57:14 -0700 Subject: [PATCH 1/7] perf(format): partial selection reads; out-of-range selections are errors MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit read_raw_data_selection computed which chunks a selection intersects, threw the answer away, decoded the entire dataset and picked elements out of it — for contiguous layouts too. A 64x64 window of a 64 MB deflate dataset cost 105 ms, about half a full read; every selection cost the same whatever its size. New partial_read module: materialise only the selection's bounding box — the overlapping rows of a contiguous dataset (straight from the file bytes) or the overlapping chunks (only those are decompressed) — then run the existing extractor over that buffer with the selection translated to the box origin, so extraction semantics are exactly the full-read ones. It declines (falling back to the old path) for All/None, compact/virtual/storage-less layouts, and boxes covering more than half the dataset. That window now takes 0.39 ms, one row 2.7 ms, one column 5.2 ms. Selections are validated against the dataset shape first. They were not: a hyperslab past an edge came back padded with zeros and a point with an out-of-range column wrapped into the next row, returning the wrong element with no error. Now FormatError::SelectionOutOfBounds (also rank mismatch and overlapping blocks); the facade's fill-aware path validates too. Tests: equivalence against a reference extraction from a full read over 60 random hyperslabs/point lists per layout (contiguous, chunked, deflate) for ranks 1-3. New read_harness bench binary with before/after in BENCHMARKS.md. Co-Authored-By: Claude Fable 5.1 --- BENCHMARKS.md | 63 +++ CHANGELOG.md | 15 + crates/clawhdf5-bench/Cargo.toml | 6 + crates/clawhdf5-bench/src/bin/read_harness.rs | 176 +++++++++ crates/clawhdf5-format/src/data_read.rs | 18 + crates/clawhdf5-format/src/error.rs | 6 + crates/clawhdf5-format/src/lib.rs | 1 + crates/clawhdf5-format/src/partial_read.rs | 359 ++++++++++++++++++ crates/clawhdf5/src/reader.rs | 1 + .../tests/partial_read_equivalence.rs | 197 ++++++++++ 10 files changed, 842 insertions(+) create mode 100644 crates/clawhdf5-bench/src/bin/read_harness.rs create mode 100644 crates/clawhdf5-format/src/partial_read.rs create mode 100644 crates/clawhdf5/tests/partial_read_equivalence.rs diff --git a/BENCHMARKS.md b/BENCHMARKS.md index e8ba435..cc118cb 100644 --- a/BENCHMARKS.md +++ b/BENCHMARKS.md @@ -28,6 +28,69 @@ --- +## Read harness + +Produced by `cargo run --release -p clawhdf5-bench --bin read_harness`: a 4096 x +2048 `f64` dataset (64 MB) written three ways, read in full and through four +hyperslab selections, each from a fresh file handle. The last column is the +point: does a selection cost what the *selection* costs? + +### Baseline (v2.4.0): every selection decodes the whole dataset + +4096 x 2048 f64 (64 MB per dataset), chunks 256 x 256, file 129 MB + +| layout | read | selected | time ms | MB/s of selection | vs full read | +|---|---|---:|---:|---:|---:| +| chunked + deflate | full (first) | 64 MB | 181.8 | 352 | | +| chunked + deflate | full (repeat) | 64 MB | 162.1 | 395 | 1.00x | +| chunked + deflate | 64 x 64 window (1 chunk) | 0.03 MB | 104.89 | 0 | 0.577x | +| chunked + deflate | 512 x 512 window (4-9 chunks) | 2.00 MB | 110.26 | 18 | 0.606x | +| chunked + deflate | one row | 0.02 MB | 105.81 | 0 | 0.582x | +| chunked + deflate | one column | 0.03 MB | 108.37 | 0 | 0.596x | +| chunked | full (first) | 64 MB | 97.4 | 657 | | +| chunked | full (repeat) | 64 MB | 86.7 | 738 | 1.00x | +| chunked | 64 x 64 window (1 chunk) | 0.03 MB | 40.97 | 1 | 0.420x | +| chunked | 512 x 512 window (4-9 chunks) | 2.00 MB | 44.66 | 45 | 0.458x | +| chunked | one row | 0.02 MB | 30.88 | 1 | 0.317x | +| chunked | one column | 0.03 MB | 30.27 | 1 | 0.311x | +| contiguous | full (first) | 64 MB | 57.6 | 1112 | | +| contiguous | full (repeat) | 64 MB | 53.5 | 1195 | 1.00x | +| contiguous | 64 x 64 window (1 chunk) | 0.03 MB | 30.97 | 1 | 0.538x | +| contiguous | 512 x 512 window (4-9 chunks) | 2.00 MB | 31.64 | 63 | 0.550x | +| contiguous | one row | 0.02 MB | 31.90 | 0 | 0.554x | +| contiguous | one column | 0.03 MB | 29.36 | 1 | 0.510x | + +### After: partial reads + +Only the rows of a contiguous dataset, or the chunks, that overlap the +selection's bounding box are read/decoded. A 64 x 64 window of the compressed +dataset: **105 -> 0.39 ms**; one row: **106 -> 2.7 ms**; one column: +**108 -> 5.2 ms**. (Absolute full-read times differ between the two runs +because the machine's speed drifted; compare the *vs full read* column.) + +4096 x 2048 f64 (64 MB per dataset), chunks 256 x 256, file 129 MB + +| layout | read | selected | time ms | MB/s of selection | vs full read | +|---|---|---:|---:|---:|---:| +| chunked + deflate | full (first) | 64 MB | 112.5 | 569 | | +| chunked + deflate | full (repeat) | 64 MB | 104.5 | 612 | 1.00x | +| chunked + deflate | 64 x 64 window (1 chunk) | 0.03 MB | 0.39 | 81 | 0.003x | +| chunked + deflate | 512 x 512 window (4-9 chunks) | 2.00 MB | 4.85 | 412 | 0.043x | +| chunked + deflate | one row | 0.02 MB | 2.69 | 6 | 0.024x | +| chunked + deflate | one column | 0.03 MB | 5.23 | 6 | 0.046x | +| chunked | full (first) | 64 MB | 70.0 | 915 | | +| chunked | full (repeat) | 64 MB | 61.7 | 1037 | 1.00x | +| chunked | 64 x 64 window (1 chunk) | 0.03 MB | 0.06 | 541 | 0.001x | +| chunked | 512 x 512 window (4-9 chunks) | 2.00 MB | 1.99 | 1005 | 0.028x | +| chunked | one row | 0.02 MB | 0.05 | 285 | 0.001x | +| chunked | one column | 0.03 MB | 0.45 | 69 | 0.006x | +| contiguous | full (first) | 64 MB | 60.3 | 1062 | | +| contiguous | full (repeat) | 64 MB | 56.4 | 1134 | 1.00x | +| contiguous | 64 x 64 window (1 chunk) | 0.03 MB | 0.08 | 396 | 0.001x | +| contiguous | 512 x 512 window (4-9 chunks) | 2.00 MB | 2.12 | 944 | 0.035x | +| contiguous | one row | 0.02 MB | 0.03 | 576 | 0.000x | +| contiguous | one column | 0.03 MB | 2.55 | 12 | 0.042x | + ## Search harness baseline (v2.3.0) Produced by `cargo run --release -p clawhdf5-bench --bin search_harness -- --full` diff --git a/CHANGELOG.md b/CHANGELOG.md index 0c75b95..41950b1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,21 @@ ## Unreleased +### HDF5 Read Path +- **Selection reads cost what the selection costs.** `read_*_selection` decoded + the *entire* dataset and then picked elements out, so a 64 x 64 window of a + 64 MB compressed dataset took 105 ms - about as long as reading all of it. + Now only the rows (contiguous) or chunks that overlap the selection's + bounding box are read and decompressed: that window takes 0.39 ms, one row + 2.7 ms, one column 5.2 ms. Results are identical to the full-read path + (equivalence-tested over random hyperslabs and point lists, ranks 1-3, + contiguous / chunked / deflate). New `read_harness` bench binary. +- **Out-of-range selections are errors.** They used to return data: a hyperslab + past an edge came back padded with zeros, and a point whose column was out of + range wrapped into the next row and returned that element. Now + `FormatError::SelectionOutOfBounds` (also for a rank mismatch or overlapping + blocks). + ### Search - `clawhdf5-ann`: **faster index builds.** Back-link pruning is 90% of a build's distance evaluations; the bulk build now inserts in batches and diff --git a/crates/clawhdf5-bench/Cargo.toml b/crates/clawhdf5-bench/Cargo.toml index f6f737d..d93530d 100644 --- a/crates/clawhdf5-bench/Cargo.toml +++ b/crates/clawhdf5-bench/Cargo.toml @@ -13,6 +13,10 @@ path = "src/bin/longmemeval_bench.rs" name = "memory_arena" path = "src/bin/memory_arena.rs" +[[bin]] +name = "read_harness" +path = "src/bin/read_harness.rs" + [[bin]] name = "search_harness" path = "src/bin/search_harness.rs" @@ -53,6 +57,8 @@ harness = false [dependencies] clawhdf5-agent = { path = "../clawhdf5-agent" } clawhdf5-ann = { path = "../clawhdf5-ann" } +clawhdf5 = { path = "../clawhdf5" } +clawhdf5-format = { path = "../clawhdf5-format" } clawhdf5-io = { path = "../clawhdf5-io" } mpi = { version = "0.8", optional = true } serde = { workspace = true } diff --git a/crates/clawhdf5-bench/src/bin/read_harness.rs b/crates/clawhdf5-bench/src/bin/read_harness.rs new file mode 100644 index 0000000..b09879a --- /dev/null +++ b/crates/clawhdf5-bench/src/bin/read_harness.rs @@ -0,0 +1,176 @@ +//! HDF5 read-path measurement harness: full reads vs. hyperslab selections on +//! a chunked 2-D dataset, compressed and uncompressed, plus a contiguous one. +//! +//! The question it answers for every read-path change: does the cost of a +//! selection scale with the *selection*, or with the whole dataset? +//! +//! ```text +//! cargo run --release -p clawhdf5-bench --bin read_harness +//! cargo run --release -p clawhdf5-bench --bin read_harness -- --large # 512 MB +//! ``` + +use std::time::{Duration, Instant}; + +use clawhdf5::{File, FileBuilder}; +use clawhdf5_format::selection::Selection; + +const CHUNK: u64 = 256; + +struct Layout { + name: &'static str, + chunked: bool, + deflate: bool, +} + +const LAYOUTS: [Layout; 3] = [ + Layout { + name: "chunked + deflate", + chunked: true, + deflate: true, + }, + Layout { + name: "chunked", + chunked: true, + deflate: false, + }, + Layout { + name: "contiguous", + chunked: false, + deflate: false, + }, +]; + +/// Smooth-ish, compressible data whose value encodes its position, so a read +/// can be verified exactly. +fn value(row: u64, col: u64) -> f64 { + (row * 100_003 + col) as f64 * 0.5 +} + +fn write_file(path: &std::path::Path, rows: u64, cols: u64) { + let data: Vec = (0..rows) + .flat_map(|r| (0..cols).map(move |c| value(r, c))) + .collect(); + let mut builder = FileBuilder::new(); + for (i, layout) in LAYOUTS.iter().enumerate() { + let ds = builder.create_dataset(&format!("d{i}")); + ds.with_f64_data(&data).with_shape(&[rows, cols]); + if layout.chunked { + ds.with_chunks(&[CHUNK, CHUNK]); + } + if layout.deflate { + ds.with_deflate(4); + } + } + builder.write(path).unwrap(); +} + +fn median(mut samples: Vec) -> Duration { + samples.sort(); + samples[samples.len() / 2] +} + +fn time(reps: usize, mut f: impl FnMut() -> T) -> Duration { + median( + (0..reps) + .map(|_| { + let t = Instant::now(); + std::hint::black_box(f()); + t.elapsed() + }) + .collect(), + ) +} + +fn slab(start: [u64; 2], count: [u64; 2]) -> Selection { + Selection::Hyperslab { + start: start.to_vec(), + stride: vec![1, 1], + count: count.to_vec(), + block: vec![1, 1], + } +} + +fn main() { + let large = std::env::args().any(|a| a == "--large"); + let (rows, cols) = if large { (8192, 8192) } else { (4096, 2048) }; + let total_mb = (rows * cols * 8) as f64 / (1 << 20) as f64; + if cfg!(debug_assertions) { + eprintln!("warning: debug build — numbers are meaningless. Use --release."); + } + + let dir = tempfile::TempDir::new().unwrap(); + let path = dir.path().join("read_harness.h5"); + write_file(&path, rows, cols); + let file_mb = std::fs::metadata(&path).unwrap().len() as f64 / (1 << 20) as f64; + + println!("## Read harness"); + println!( + "\n{rows} x {cols} f64 ({total_mb:.0} MB per dataset), chunks {CHUNK} x {CHUNK}, file {file_mb:.0} MB\n" + ); + + // (label, selection, elements selected) + let selections: Vec<(&str, Selection, u64)> = vec![ + ( + "64 x 64 window (1 chunk)", + slab([300, 300], [64, 64]), + 64 * 64, + ), + ( + "512 x 512 window (4-9 chunks)", + slab([1000, 700], [512, 512]), + 512 * 512, + ), + ("one row", slab([rows / 2, 0], [1, cols]), cols), + ("one column", slab([0, cols / 2], [rows, 1]), rows), + ]; + + println!("| layout | read | selected | time ms | MB/s of selection | vs full read |"); + println!("|---|---|---:|---:|---:|---:|"); + for (i, layout) in LAYOUTS.iter().enumerate() { + // Fresh handle per layout so one dataset's cached chunks don't help + // (or evict) another's. + let file = File::open(&path).unwrap(); + let ds = file.dataset(&format!("d{i}")).unwrap(); + + let full_cold = time(1, || ds.read_f64().unwrap()); + let full = time(3, || ds.read_f64().unwrap()); + println!( + "| {} | full (first) | {total_mb:.0} MB | {:.1} | {:.0} | |", + layout.name, + full_cold.as_secs_f64() * 1e3, + total_mb / full_cold.as_secs_f64() + ); + println!( + "| {} | full (repeat) | {total_mb:.0} MB | {:.1} | {:.0} | 1.00x |", + layout.name, + full.as_secs_f64() * 1e3, + total_mb / full.as_secs_f64() + ); + + for (label, selection, elements) in &selections { + // A fresh handle again: measure the selection on its own, not + // served from chunks the full read just cached. + let file = File::open(&path).unwrap(); + let ds = file.dataset(&format!("d{i}")).unwrap(); + let got = ds.read_f64_selection(selection).unwrap(); + assert_eq!(got.len() as u64, *elements, "{label}"); + if let Selection::Hyperslab { start, .. } = selection { + assert_eq!(got[0], value(start[0], start[1]), "{label}: wrong data"); + } + let took = time(5, || { + let file = File::open(&path).unwrap(); + let ds = file.dataset(&format!("d{i}")).unwrap(); + ds.read_f64_selection(selection).unwrap() + }); + let mb = (*elements * 8) as f64 / (1 << 20) as f64; + println!( + "| {} | {label} | {:.2} MB | {:.2} | {:.0} | {:.3}x |", + layout.name, + mb, + took.as_secs_f64() * 1e3, + mb / took.as_secs_f64(), + took.as_secs_f64() / full_cold.as_secs_f64() + ); + } + } +} diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index e6c4271..13b4953 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -307,6 +307,24 @@ pub fn read_raw_data_selection( ) -> Result, FormatError> { use crate::selection::Selection; + crate::partial_read::validate(selection, &dataspace.dimensions)?; + + // Read only what the selection's bounding box touches when that is + // possible; everything below is the decode-everything-then-pick path, + // kept for the cases `partial_read` declines. + if let Some(selected) = crate::partial_read::read_selection( + file_data, + layout, + dataspace, + datatype.type_size() as usize, + pipeline, + offset_size, + length_size, + selection, + )? { + return Ok(selected); + } + match selection { Selection::All => { return read_raw_data_full( diff --git a/crates/clawhdf5-format/src/error.rs b/crates/clawhdf5-format/src/error.rs index a8fa0ec..6d10c3d 100644 --- a/crates/clawhdf5-format/src/error.rs +++ b/crates/clawhdf5-format/src/error.rs @@ -117,6 +117,9 @@ pub enum FormatError { /// A message is marked shared but was parsed without access to the file, /// so the reference to the real message could not be followed. UnresolvedSharedMessage, + /// A selection does not fit the dataset it was applied to (wrong rank, or + /// it reaches past a dimension's extent). + SelectionOutOfBounds(String), /// The dataset's raw data is stored in external files (External Data /// Files message), which this reader does not follow. ExternalDataFilesUnsupported, @@ -333,6 +336,9 @@ impl fmt::Display for FormatError { f, "dataset raw data is stored in external file(s), which is not supported" ), + FormatError::SelectionOutOfBounds(msg) => { + write!(f, "selection out of bounds: {msg}") + } FormatError::UnresolvedSharedMessage => write!( f, "message is shared but no file data was available to resolve it" diff --git a/crates/clawhdf5-format/src/lib.rs b/crates/clawhdf5-format/src/lib.rs index dae44b6..83e9740 100644 --- a/crates/clawhdf5-format/src/lib.rs +++ b/crates/clawhdf5-format/src/lib.rs @@ -89,6 +89,7 @@ pub mod object_header; pub mod object_header_writer; #[cfg(feature = "parallel")] pub mod parallel_read; +pub mod partial_read; pub mod profiling; pub mod property_list; pub mod selection; diff --git a/crates/clawhdf5-format/src/partial_read.rs b/crates/clawhdf5-format/src/partial_read.rs new file mode 100644 index 0000000..7cdd2de --- /dev/null +++ b/crates/clawhdf5-format/src/partial_read.rs @@ -0,0 +1,359 @@ +//! Selection reads that cost what the selection costs, not what the dataset +//! costs. +//! +//! [`crate::data_read::read_raw_data_selection`] used to decode the *entire* +//! dataset and then pick elements out of it, so reading a 64x64 window of a +//! large dataset took about as long as reading all of it. Here the selection's +//! bounding box is materialised instead — only the rows of a contiguous +//! dataset, or only the chunks, that overlap it — and the existing extractor +//! runs over that small buffer with the selection translated to the box's +//! origin. Extraction semantics are therefore exactly the full-read ones. + +#[cfg(not(feature = "std"))] +use alloc::string as alloc_or_std; +#[cfg(not(feature = "std"))] +use alloc::{format, vec, vec::Vec}; +#[cfg(feature = "std")] +use std::string as alloc_or_std; + +use crate::chunked_read::{alloc_output, checked_byte_len, list_chunks}; +use crate::data_layout::DataLayout; +use crate::data_read::extract_selection_from_buffer; +use crate::dataspace::Dataspace; +use crate::error::FormatError; +use crate::filter_pipeline::FilterPipeline; +use crate::filters::decompress_chunk; +use crate::selection::Selection; + +/// The smallest axis-aligned box containing every selected element, as +/// `(start, extent)` per dimension. `None` when there is nothing to gain or +/// the selection is not valid for `dims` (the caller's full path then reports +/// the error exactly as before). +fn bounding_box(selection: &Selection, dims: &[u64]) -> Option<(Vec, Vec)> { + match selection { + Selection::Hyperslab { + start, + stride, + count, + block, + } => { + let rank = dims.len(); + if [start.len(), stride.len(), count.len(), block.len()] != [rank; 4] { + return None; + } + let mut extent = Vec::with_capacity(rank); + for d in 0..rank { + if count[d] == 0 || block[d] == 0 { + return None; + } + // Last selected index + 1, relative to start. + let span = (count[d] - 1) + .checked_mul(stride[d])? + .checked_add(block[d])?; + if start[d].checked_add(span)? > dims[d] { + return None; + } + extent.push(span); + } + Some((start.clone(), extent)) + } + Selection::Points(points) => { + let rank = dims.len(); + let first = points.first()?; + if first.len() != rank { + return None; + } + let (mut lo, mut hi) = (first.clone(), first.clone()); + for p in points { + if p.len() != rank { + return None; + } + for d in 0..rank { + if p[d] >= dims[d] { + return None; + } + lo[d] = lo[d].min(p[d]); + hi[d] = hi[d].max(p[d]); + } + } + let extent = lo.iter().zip(&hi).map(|(l, h)| h - l + 1).collect(); + Some((lo, extent)) + } + Selection::All | Selection::None => None, + } +} + +/// Check that `selection` addresses only elements that exist in a dataset of +/// shape `dims`. Without this an out-of-range selection read *something*: a +/// hyperslab past the edge came back padded with zeros, and a point whose +/// column was out of range wrapped into the next row. +pub fn validate(selection: &Selection, dims: &[u64]) -> Result<(), FormatError> { + let rank = dims.len(); + let bad = |msg: alloc_or_std::String| Err(FormatError::SelectionOutOfBounds(msg)); + match selection { + Selection::All | Selection::None => Ok(()), + Selection::Hyperslab { + start, + stride, + count, + block, + } => { + if [start.len(), stride.len(), count.len(), block.len()] != [rank; 4] { + return bad(format!("hyperslab rank does not match dataset rank {rank}")); + } + for d in 0..rank { + if count[d] == 0 || block[d] == 0 { + continue; // selects nothing along this dimension + } + let end = (count[d] - 1) + .checked_mul(stride[d]) + .and_then(|v| v.checked_add(block[d])) + .and_then(|v| v.checked_add(start[d])); + if !end.is_some_and(|end| end <= dims[d]) { + return bad(format!( + "dimension {d}: start {} stride {} count {} block {} exceeds extent {}", + start[d], stride[d], count[d], block[d], dims[d] + )); + } + if block[d] > stride[d] && count[d] > 1 { + return bad(format!( + "dimension {d}: block {} larger than stride {} (overlapping blocks)", + block[d], stride[d] + )); + } + } + Ok(()) + } + Selection::Points(points) => { + for p in points { + if p.len() != rank { + return bad(format!("point {p:?} does not match dataset rank {rank}")); + } + if let Some(d) = (0..rank).find(|&d| p[d] >= dims[d]) { + return bad(format!( + "point {p:?}: coordinate {} exceeds extent {} of dimension {d}", + p[d], dims[d] + )); + } + } + Ok(()) + } + } +} + +/// The same selection expressed relative to `origin`. +fn translate(selection: &Selection, origin: &[u64]) -> Selection { + match selection { + Selection::Hyperslab { + start, + stride, + count, + block, + } => Selection::Hyperslab { + start: start.iter().zip(origin).map(|(s, o)| s - o).collect(), + stride: stride.clone(), + count: count.clone(), + block: block.clone(), + }, + Selection::Points(points) => Selection::Points( + points + .iter() + .map(|p| p.iter().zip(origin).map(|(c, o)| c - o).collect()) + .collect(), + ), + other => other.clone(), + } +} + +/// Copy the part of a source region that overlaps the box into `out` (which +/// is the box, row-major). +/// +/// The source region starts at `src_origin` in dataset coordinates, has shape +/// `src_shape`, and its elements are in `src` row-major. One `memcpy` per +/// overlapping row of the last dimension. +#[allow(clippy::too_many_arguments)] +fn copy_overlap( + src: &[u8], + src_origin: &[u64], + src_shape: &[u64], + out: &mut [u8], + box_start: &[u64], + box_extent: &[u64], + elem_size: usize, +) { + let rank = box_start.len(); + // Overlap in dataset coordinates. + let mut lo = vec![0u64; rank]; + let mut hi = vec![0u64; rank]; + for d in 0..rank { + lo[d] = src_origin[d].max(box_start[d]); + hi[d] = (src_origin[d] + src_shape[d]).min(box_start[d] + box_extent[d]); + if lo[d] >= hi[d] { + return; + } + } + let strides = |shape: &[u64]| { + let mut s = vec![1u64; rank]; + for d in (0..rank.saturating_sub(1)).rev() { + s[d] = s[d + 1] * shape[d + 1]; + } + s + }; + let (src_strides, out_strides) = (strides(src_shape), strides(box_extent)); + let last = rank - 1; + let run = ((hi[last] - lo[last]) as usize) * elem_size; + + let mut idx = lo.clone(); + loop { + let src_at: u64 = (0..rank) + .map(|d| (idx[d] - src_origin[d]) * src_strides[d]) + .sum(); + let out_at: u64 = (0..rank) + .map(|d| (idx[d] - box_start[d]) * out_strides[d]) + .sum(); + let (s, o) = (src_at as usize * elem_size, out_at as usize * elem_size); + if let (Some(from), Some(to)) = (src.get(s..s + run), out.get_mut(o..o + run)) { + to.copy_from_slice(from); + } + // Advance over every dimension but the last. + let mut d = last; + loop { + if d == 0 { + return; + } + d -= 1; + idx[d] += 1; + if idx[d] < hi[d] { + break; + } + idx[d] = lo[d]; + } + } +} + +/// Read `selection` without materialising the whole dataset, when that is +/// possible and worthwhile. `Ok(None)` means "use the full-read path": an +/// `All`/`None`/invalid selection, a layout this doesn't handle (compact, +/// virtual, storage-less), or a bounding box covering most of the dataset. +#[allow(clippy::too_many_arguments)] +pub fn read_selection( + file_data: &[u8], + layout: &DataLayout, + dataspace: &Dataspace, + elem_size: usize, + pipeline: Option<&FilterPipeline>, + offset_size: u8, + length_size: u8, + selection: &Selection, +) -> Result>, FormatError> { + let dims = &dataspace.dimensions; + if dims.is_empty() || elem_size == 0 { + return Ok(None); + } + let Some((box_start, box_extent)) = bounding_box(selection, dims) else { + return Ok(None); + }; + let total = dataspace.checked_num_elements()?; + let box_elements = box_extent + .iter() + .try_fold(1u64, |acc, &e| acc.checked_mul(e)) + .ok_or_else(|| FormatError::Overflow("selection bounding box overflows".into()))?; + // A box covering most of the dataset gains nothing over the full path. + if box_elements.saturating_mul(2) > total { + return Ok(None); + } + let mut boxed = alloc_output(checked_byte_len(box_elements, elem_size)?)?; + + match layout { + DataLayout::Contiguous { + address: Some(address), + .. + } => { + let base = usize::try_from(*address) + .map_err(|_| FormatError::Overflow("data address exceeds usize".into()))?; + let data = file_data + .get(base..) + .and_then(|d| d.get(..checked_byte_len(total, elem_size).ok()?)) + .ok_or(FormatError::UnexpectedEof { + expected: base, + available: file_data.len(), + })?; + let origin = vec![0u64; dims.len()]; + copy_overlap( + data, + &origin, + dims, + &mut boxed, + &box_start, + &box_extent, + elem_size, + ); + } + DataLayout::Chunked { + btree_address: Some(_), + .. + } => { + let (chunks, chunk_dims) = list_chunks( + file_data, + layout, + dataspace, + elem_size, + offset_size, + length_size, + )?; + let rank = dims.len(); + let chunk_shape: Vec = chunk_dims.iter().map(|&d| d as u64).collect(); + let chunk_bytes = crate::chunked_read::checked_chunk_byte_len(&chunk_dims, elem_size)?; + for chunk in &chunks { + if chunk.offsets.len() < rank || chunk.address == u64::MAX { + continue; + } + let origin = &chunk.offsets[..rank]; + let overlaps = (0..rank).all(|d| { + origin[d] < box_start[d] + box_extent[d] + && origin[d].saturating_add(chunk_shape[d]) > box_start[d] + }); + if !overlaps { + continue; + } + let at = usize::try_from(chunk.address) + .map_err(|_| FormatError::Overflow("chunk address exceeds usize".into()))?; + let raw = at + .checked_add(chunk.chunk_size as usize) + .and_then(|end| file_data.get(at..end)) + .ok_or(FormatError::UnexpectedEof { + expected: at.saturating_add(chunk.chunk_size as usize), + available: file_data.len(), + })?; + // Mirrors the full-read path: a non-zero filter mask means the + // chunk was stored unfiltered. + let decoded; + let data: &[u8] = match pipeline { + Some(pl) if chunk.filter_mask == 0 => { + decoded = decompress_chunk(raw, pl, chunk_bytes, elem_size as u32)?; + &decoded + } + _ => raw, + }; + copy_overlap( + data, + origin, + &chunk_shape, + &mut boxed, + &box_start, + &box_extent, + elem_size, + ); + } + } + _ => return Ok(None), + } + + extract_selection_from_buffer( + &boxed, + &box_extent, + elem_size, + &translate(selection, &box_start), + ) + .map(Some) +} diff --git a/crates/clawhdf5/src/reader.rs b/crates/clawhdf5/src/reader.rs index fb17557..d6d6032 100644 --- a/crates/clawhdf5/src/reader.rs +++ b/crates/clawhdf5/src/reader.rs @@ -458,6 +458,7 @@ impl<'f> Dataset<'f> { || (matches!(dl, DataLayout::Chunked { .. }) && !clawhdf5_format::fill_value::is_default(fill.as_deref())); if fill_matters { + clawhdf5_format::partial_read::validate(selection, &ds.dimensions)?; let full = self.read_raw()?; return Ok(data_read::extract_selection_from_buffer( &full, diff --git a/crates/clawhdf5/tests/partial_read_equivalence.rs b/crates/clawhdf5/tests/partial_read_equivalence.rs new file mode 100644 index 0000000..a84b220 --- /dev/null +++ b/crates/clawhdf5/tests/partial_read_equivalence.rs @@ -0,0 +1,197 @@ +//! Selection reads must return exactly what a full read followed by element +//! extraction returns — for every layout, rank and selection shape — while +//! touching only what the selection needs. + +use clawhdf5::{File, FileBuilder}; +use clawhdf5_format::selection::Selection; + +struct Rng(u64); +impl Rng { + fn next(&mut self) -> u64 { + self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15); + let mut z = self.0; + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) + } + fn below(&mut self, n: u64) -> u64 { + self.next() % n.max(1) + } +} + +/// Row-major reference extraction from a full read. +fn reference(full: &[i32], dims: &[u64], selection: &Selection) -> Vec { + let strides: Vec = (0..dims.len()) + .map(|d| dims[d + 1..].iter().product()) + .collect(); + let at = + |coord: &[u64]| full[coord.iter().zip(&strides).map(|(c, s)| c * s).sum::() as usize]; + match selection { + Selection::Points(points) => points.iter().map(|p| at(p)).collect(), + Selection::Hyperslab { + start, + stride, + count, + block, + } => { + // Selected indices per dimension, then their cartesian product. + let per_dim: Vec> = (0..dims.len()) + .map(|d| { + (0..count[d]) + .flat_map(|c| (0..block[d]).map(move |b| (c, b))) + .map(|(c, b)| start[d] + c * stride[d] + b) + .collect() + }) + .collect(); + let mut out = Vec::new(); + let mut idx = vec![0usize; dims.len()]; + loop { + let coord: Vec = idx + .iter() + .enumerate() + .map(|(d, &i)| per_dim[d][i]) + .collect(); + out.push(at(&coord)); + let mut d = dims.len(); + loop { + if d == 0 { + return out; + } + d -= 1; + idx[d] += 1; + if idx[d] < per_dim[d].len() { + break; + } + idx[d] = 0; + } + } + } + _ => unreachable!(), + } +} + +fn random_hyperslab(rng: &mut Rng, dims: &[u64]) -> Selection { + let mut start = Vec::new(); + let mut stride = Vec::new(); + let mut count = Vec::new(); + let mut block = Vec::new(); + for &dim in dims { + let b = 1 + rng.below(3); + let st = b + rng.below(4); // stride >= block: no overlap + let s = rng.below(dim - b + 1); + let max_count = (dim - s - b) / st + 1; + let c = 1 + rng.below(max_count.min(6)); + start.push(s); + stride.push(st); + count.push(c); + block.push(b); + } + Selection::Hyperslab { + start, + stride, + count, + block, + } +} + +#[test] +fn selection_reads_match_full_reads_for_every_layout() { + let dir = tempfile::tempdir().unwrap(); + let mut rng = Rng(7); + // (dims, chunk dims) + let shapes: [(&[u64], &[u64]); 3] = [ + (&[97], &[10]), + (&[41, 53], &[8, 9]), + (&[11, 13, 17], &[4, 5, 6]), + ]; + for (dims, chunks) in shapes { + let n: u64 = dims.iter().product(); + let data: Vec = (0..n as i32).map(|v| v * 3 - 7).collect(); + + let path = dir.path().join(format!("r{}.h5", dims.len())); + let mut builder = FileBuilder::new(); + builder + .create_dataset("contiguous") + .with_i32_data(&data) + .with_shape(dims); + builder + .create_dataset("chunked") + .with_i32_data(&data) + .with_shape(dims) + .with_chunks(chunks); + builder + .create_dataset("deflated") + .with_i32_data(&data) + .with_shape(dims) + .with_chunks(chunks) + .with_deflate(3); + builder.write(&path).unwrap(); + + let file = File::open(&path).unwrap(); + for name in ["contiguous", "chunked", "deflated"] { + let ds = file.dataset(name).unwrap(); + let full = ds.read_i32().unwrap(); + assert_eq!(full, data, "{name} full read"); + + for case in 0..60 { + let selection = if case % 5 == 4 { + let points = (0..1 + rng.below(12)) + .map(|_| dims.iter().map(|&d| rng.below(d)).collect()) + .collect(); + Selection::Points(points) + } else { + random_hyperslab(&mut rng, dims) + }; + assert_eq!( + ds.read_i32_selection(&selection).unwrap(), + reference(&full, dims, &selection), + "{name} rank {} case {case}: {selection:?}", + dims.len() + ); + } + } + } +} + +#[test] +fn out_of_bounds_selections_are_errors() { + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("oob.h5"); + let mut builder = FileBuilder::new(); + builder + .create_dataset("d") + .with_i32_data(&(0..100).collect::>()) + .with_shape(&[10, 10]) + .with_chunks(&[4, 4]); + builder.write(&path).unwrap(); + let file = File::open(&path).unwrap(); + let ds = file.dataset("d").unwrap(); + let beyond = Selection::Hyperslab { + start: vec![8, 8], + stride: vec![1, 1], + count: vec![5, 5], + block: vec![1, 1], + }; + use clawhdf5::Error; + use clawhdf5_format::error::FormatError; + let is_oob = |s: &Selection| { + matches!( + ds.read_i32_selection(s), + Err(Error::Format(FormatError::SelectionOutOfBounds(_))) + ) + }; + // Used to come back padded with zeros. + assert!(is_oob(&beyond)); + // Row out of range. + assert!(is_oob(&Selection::Points(vec![vec![10, 0]]))); + // Column out of range: used to wrap into the next row and return its value. + assert!(is_oob(&Selection::Points(vec![vec![0, 12]]))); + // Wrong rank. + assert!(is_oob(&Selection::Points(vec![vec![3]]))); + // In range is fine. + assert_eq!( + ds.read_i32_selection(&Selection::Points(vec![vec![9, 9]])) + .unwrap(), + [99] + ); +} From d668e45ab54f1db82a8f538119761371ae87f72c Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 19 Sep 2026 13:59:18 -0700 Subject: [PATCH 2/7] feat(format): read chunked datasets indexed by a version-2 B-tree MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit With libver='latest', a chunked dataset with two or more unlimited dimensions indexes its chunks with a v2 B-tree (layout v4, index type 5). Reading one failed with "unsupported chunked layout version=4, index_type=Some(5)". read_btree_v2_chunks decodes record types 10 (address + scaled offsets) and 11 (address, stored size, filter mask, scaled offsets). The width of the stored-size field is taken from the record size the tree header declares rather than re-deriving the library's formula. Scaled offsets are multiplied back by the chunk dimensions with overflow checks. The chunk-index dispatch existed four times (uncached, cached, sweep and indexed readers). The three copies outside list_chunks now call it, so every read path — and fill-value handling and partial reads — supports every index type from one place. h5py interop test: plain, gzip+shuffle, a 2500-chunk tree with internal nodes, a sparse dataset with a fill value, and a strided hyperslab. Co-Authored-By: Claude Fable 5.1 --- CHANGELOG.md | 4 + crates/clawhdf5-format/src/chunked_read.rs | 400 ++++++++------------ crates/clawhdf5/tests/h5py_interop_tests.rs | 66 ++++ docs/known-issues.md | 15 +- 4 files changed, 229 insertions(+), 256 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 41950b1..c1eed95 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,10 @@ 2.7 ms, one column 5.2 ms. Results are identical to the full-read path (equivalence-tested over random hyperslabs and point lists, ranks 1-3, contiguous / chunked / deflate). New `read_harness` bench binary. +- **Datasets indexed by a version-2 B-tree now read** (layout v4, chunk index + type 5 — what `libver='latest'` uses for two or more unlimited dimensions; + previously "unsupported chunked layout"). The four copies of the chunk-index + dispatch are now one shared function, so every read path gets it. - **Out-of-range selections are errors.** They used to return data: a hyperslab past an edge came back padded with zeros, and a point whose column was out of range wrapped into the next row and returned that element. Now diff --git a/crates/clawhdf5-format/src/chunked_read.rs b/crates/clawhdf5-format/src/chunked_read.rs index 8ba9332..0dbf6d4 100644 --- a/crates/clawhdf5-format/src/chunked_read.rs +++ b/crates/clawhdf5-format/src/chunked_read.rs @@ -362,6 +362,112 @@ pub fn generate_implicit_chunks( } /// Read a chunked dataset, decompressing chunks as needed. +/// B-tree v2 record types used for chunk indexing. +const BT2_CHUNK_UNFILTERED: u8 = 10; +const BT2_CHUNK_FILTERED: u8 = 11; + +/// Chunks indexed by a version-2 B-tree (layout v4, index type 5). +/// +/// Record layouts (all little endian): +/// * type 10, unfiltered: address, then one 8-byte *scaled* offset per +/// dimension (offset / chunk dimension); +/// * type 11, filtered: address, stored chunk size (a variable number of +/// bytes), 4-byte filter mask, then the scaled offsets. +/// +/// The width of the stored-size field depends on the largest possible chunk; +/// rather than re-derive the library's formula it is taken from the record +/// size the tree header declares, which is what actually governs the bytes. +fn read_btree_v2_chunks( + file_data: &[u8], + addr: u64, + chunk_dims: &[usize], + elem_size: usize, + offset_size: u8, + length_size: u8, +) -> Result, FormatError> { + use crate::btree_v2::{BTreeV2Header, collect_btree_v2_records}; + + let bad = |what: &str| FormatError::ChunkedReadError(format!("B-tree v2 chunk index: {what}")); + let header = BTreeV2Header::parse(file_data, addr as usize, offset_size, length_size)?; + let rank = chunk_dims.len(); + let os = offset_size as usize; + let record_size = header.record_size as usize; + let size_len = match header.tree_type { + BT2_CHUNK_UNFILTERED => { + if record_size != os + 8 * rank { + return Err(bad("unexpected record size for unfiltered chunks")); + } + 0 + } + BT2_CHUNK_FILTERED => { + let fixed = os + 4 + 8 * rank; + let size_len = record_size + .checked_sub(fixed) + .ok_or_else(|| bad("record too small"))?; + if !(1..=8).contains(&size_len) { + return Err(bad("implausible chunk-size field width")); + } + size_len + } + _ => return Err(bad("tree is not a chunk index")), + }; + let unfiltered_bytes = checked_chunk_byte_len(chunk_dims, elem_size)?; + let unfiltered_bytes = + u32::try_from(unfiltered_bytes).map_err(|_| bad("chunk larger than 4 GiB"))?; + + let records = collect_btree_v2_records(file_data, &header, offset_size, length_size)?; + let mut chunks = Vec::with_capacity(records.len()); + for record in &records { + let data = record.data.as_slice(); + if data.len() < record_size { + return Err(bad("truncated record")); + } + let address = read_offset(data, 0, offset_size)?; + let mut pos = os; + let (chunk_size, filter_mask) = if size_len == 0 { + (unfiltered_bytes, 0) + } else { + let mut size = 0u64; + for (i, &b) in data[pos..pos + size_len].iter().enumerate() { + size |= u64::from(b) << (8 * i); + } + pos += size_len; + let mask = u32::from_le_bytes([data[pos], data[pos + 1], data[pos + 2], data[pos + 3]]); + pos += 4; + ( + u32::try_from(size).map_err(|_| bad("stored chunk larger than 4 GiB"))?, + mask, + ) + }; + let mut offsets = Vec::with_capacity(rank); + for &dim in chunk_dims { + let scaled = u64::from_le_bytes([ + data[pos], + data[pos + 1], + data[pos + 2], + data[pos + 3], + data[pos + 4], + data[pos + 5], + data[pos + 6], + data[pos + 7], + ]); + pos += 8; + offsets.push( + scaled + .checked_mul(dim as u64) + .ok_or_else(|| bad("chunk offset overflows"))?, + ); + } + chunks.push(ChunkInfo { + chunk_size, + filter_mask, + offsets, + address, + }); + } + Ok(chunks) +} + /// Every allocated chunk of a chunked dataset, for any supported chunk index, /// plus the spatial chunk dimensions. Chunks the file never allocated (sparse /// datasets) are simply absent from the list. @@ -487,6 +593,18 @@ pub fn list_chunks( length_size, )? } + (4, Some(5)) => { + // Version-2 B-tree: what the library uses for a dataset with two + // or more unlimited dimensions. + read_btree_v2_chunks( + file_data, + addr, + &chunk_dims, + elem_size, + offset_size, + length_size, + )? + } (v, idx) => { return Err(FormatError::ChunkedReadError(format!( "unsupported chunked layout version={v}, index_type={idx:?}" @@ -629,29 +747,12 @@ pub fn read_chunked_data_cached( length_size: u8, cache: &ChunkCache, ) -> Result, FormatError> { - let ( - chunk_dimensions, - version, - chunk_index_type, - addr_opt, - single_filtered_size, - single_filter_mask, - ) = match layout { + let (chunk_dimensions, addr_opt) = match layout { DataLayout::Chunked { chunk_dimensions, btree_address, - version, - chunk_index_type, - single_chunk_filtered_size, - single_chunk_filter_mask, - } => ( - chunk_dimensions, - *version, - *chunk_index_type, - *btree_address, - *single_chunk_filtered_size, - *single_chunk_filter_mask, - ), + .. + } => (chunk_dimensions, *btree_address), _ => { return Err(FormatError::ChunkedReadError( "expected chunked layout".into(), @@ -688,69 +789,14 @@ pub fn read_chunked_data_cached( // Populate chunk index on first access if !cache.has_index() { - let chunks = match (version, chunk_index_type) { - (3, _) => collect_chunk_info(file_data, addr, ndims, offset_size, length_size)?, - (4, Some(1)) => { - let chunk_byte_size = checked_chunk_byte_len(&chunk_dims, elem_size)?; - let (csize, fmask) = if let Some(fs) = single_filtered_size { - (fs as u32, single_filter_mask.unwrap_or(0)) - } else { - (chunk_byte_size as u32, 0) - }; - vec![ChunkInfo { - chunk_size: csize, - filter_mask: fmask, - offsets: vec![0u64; rank], - address: addr, - }] - } - (4, Some(2)) => { - let spatial_chunk_dims: &[u32] = &chunk_dimensions[..rank]; - generate_implicit_chunks( - addr, - &dataspace.dimensions, - spatial_chunk_dims, - elem_size as u32, - ) - } - (4, Some(3)) => { - let spatial_chunk_dims: &[u32] = &chunk_dimensions[..rank]; - let header = - FixedArrayHeader::parse(file_data, addr as usize, offset_size, length_size)?; - read_fixed_array_chunks( - file_data, - &header, - &dataspace.dimensions, - spatial_chunk_dims, - elem_size as u32, - offset_size, - length_size, - )? - } - (4, Some(4)) => { - let spatial_chunk_dims: &[u32] = &chunk_dimensions[..rank]; - let header = ExtensibleArrayHeader::parse( - file_data, - addr as usize, - offset_size, - length_size, - )?; - read_extensible_array_chunks( - file_data, - &header, - &dataspace.dimensions, - spatial_chunk_dims, - elem_size as u32, - offset_size, - length_size, - )? - } - (v, idx) => { - return Err(FormatError::ChunkedReadError(format!( - "unsupported chunked layout version={v}, index_type={idx:?}" - ))); - } - }; + let (chunks, _) = list_chunks( + file_data, + layout, + dataspace, + elem_size, + offset_size, + length_size, + )?; cache.populate_index(&chunks, rank); } @@ -985,29 +1031,12 @@ pub fn read_chunked_data_sweep( cache: &ChunkCache, sweep: &mut SweepContext, ) -> Result, FormatError> { - let ( - chunk_dimensions, - version, - chunk_index_type, - addr_opt, - single_filtered_size, - single_filter_mask, - ) = match layout { + let (chunk_dimensions, addr_opt) = match layout { DataLayout::Chunked { chunk_dimensions, btree_address, - version, - chunk_index_type, - single_chunk_filtered_size, - single_chunk_filter_mask, - } => ( - chunk_dimensions, - *version, - *chunk_index_type, - *btree_address, - *single_chunk_filtered_size, - *single_chunk_filter_mask, - ), + .. + } => (chunk_dimensions, *btree_address), _ => { return Err(FormatError::ChunkedReadError( "expected chunked layout".into(), @@ -1044,69 +1073,14 @@ pub fn read_chunked_data_sweep( // Populate chunk index on first access if !cache.has_index() { - let chunks = match (version, chunk_index_type) { - (3, _) => collect_chunk_info(file_data, addr, ndims, offset_size, length_size)?, - (4, Some(1)) => { - let chunk_byte_size = checked_chunk_byte_len(&chunk_dims, elem_size)?; - let (csize, fmask) = if let Some(fs) = single_filtered_size { - (fs as u32, single_filter_mask.unwrap_or(0)) - } else { - (chunk_byte_size as u32, 0) - }; - vec![ChunkInfo { - chunk_size: csize, - filter_mask: fmask, - offsets: vec![0u64; rank], - address: addr, - }] - } - (4, Some(2)) => { - let spatial_chunk_dims: &[u32] = &chunk_dimensions[..rank]; - generate_implicit_chunks( - addr, - &dataspace.dimensions, - spatial_chunk_dims, - elem_size as u32, - ) - } - (4, Some(3)) => { - let spatial_chunk_dims: &[u32] = &chunk_dimensions[..rank]; - let header = - FixedArrayHeader::parse(file_data, addr as usize, offset_size, length_size)?; - read_fixed_array_chunks( - file_data, - &header, - &dataspace.dimensions, - spatial_chunk_dims, - elem_size as u32, - offset_size, - length_size, - )? - } - (4, Some(4)) => { - let spatial_chunk_dims: &[u32] = &chunk_dimensions[..rank]; - let header = ExtensibleArrayHeader::parse( - file_data, - addr as usize, - offset_size, - length_size, - )?; - read_extensible_array_chunks( - file_data, - &header, - &dataspace.dimensions, - spatial_chunk_dims, - elem_size as u32, - offset_size, - length_size, - )? - } - (v, idx) => { - return Err(FormatError::ChunkedReadError(format!( - "unsupported chunked layout version={v}, index_type={idx:?}" - ))); - } - }; + let (chunks, _) = list_chunks( + file_data, + layout, + dataspace, + elem_size, + offset_size, + length_size, + )?; cache.populate_index(&chunks, rank); } @@ -1211,29 +1185,12 @@ pub fn read_chunked_data_indexed( length_size: u8, cache: &ChunkCache, ) -> Result, FormatError> { - let ( - chunk_dimensions, - version, - chunk_index_type, - addr_opt, - single_filtered_size, - single_filter_mask, - ) = match layout { + let (chunk_dimensions, addr_opt) = match layout { DataLayout::Chunked { chunk_dimensions, btree_address, - version, - chunk_index_type, - single_chunk_filtered_size, - single_chunk_filter_mask, - } => ( - chunk_dimensions, - *version, - *chunk_index_type, - *btree_address, - *single_chunk_filtered_size, - *single_chunk_filter_mask, - ), + .. + } => (chunk_dimensions, *btree_address), _ => { return Err(FormatError::ChunkedReadError( "expected chunked layout".into(), @@ -1270,69 +1227,14 @@ pub fn read_chunked_data_indexed( // Build chunk index on first access if !cache.has_chunk_index() { - let chunks = match (version, chunk_index_type) { - (3, _) => collect_chunk_info(file_data, addr, ndims, offset_size, length_size)?, - (4, Some(1)) => { - let chunk_byte_size = checked_chunk_byte_len(&chunk_dims, elem_size)?; - let (csize, fmask) = if let Some(fs) = single_filtered_size { - (fs as u32, single_filter_mask.unwrap_or(0)) - } else { - (chunk_byte_size as u32, 0) - }; - vec![ChunkInfo { - chunk_size: csize, - filter_mask: fmask, - offsets: vec![0u64; rank], - address: addr, - }] - } - (4, Some(2)) => { - let spatial_chunk_dims: &[u32] = &chunk_dimensions[..rank]; - generate_implicit_chunks( - addr, - &dataspace.dimensions, - spatial_chunk_dims, - elem_size as u32, - ) - } - (4, Some(3)) => { - let spatial_chunk_dims: &[u32] = &chunk_dimensions[..rank]; - let header = - FixedArrayHeader::parse(file_data, addr as usize, offset_size, length_size)?; - read_fixed_array_chunks( - file_data, - &header, - &dataspace.dimensions, - spatial_chunk_dims, - elem_size as u32, - offset_size, - length_size, - )? - } - (4, Some(4)) => { - let spatial_chunk_dims: &[u32] = &chunk_dimensions[..rank]; - let header = ExtensibleArrayHeader::parse( - file_data, - addr as usize, - offset_size, - length_size, - )?; - read_extensible_array_chunks( - file_data, - &header, - &dataspace.dimensions, - spatial_chunk_dims, - elem_size as u32, - offset_size, - length_size, - )? - } - (v, idx) => { - return Err(FormatError::ChunkedReadError(format!( - "unsupported chunked layout version={v}, index_type={idx:?}" - ))); - } - }; + let (chunks, _) = list_chunks( + file_data, + layout, + dataspace, + elem_size, + offset_size, + length_size, + )?; cache.populate_chunk_index(&chunks, rank); // Also populate the legacy index for compatibility if !cache.has_index() { diff --git a/crates/clawhdf5/tests/h5py_interop_tests.rs b/crates/clawhdf5/tests/h5py_interop_tests.rs index 2402a07..b5692d3 100644 --- a/crates/clawhdf5/tests/h5py_interop_tests.rs +++ b/crates/clawhdf5/tests/h5py_interop_tests.rs @@ -913,3 +913,69 @@ with h5py.File("{dst_str}", "r") as f: "[(1, 2.5), (3, 4.5)] ('a', 'b') [18446744073709551615, 0, 9223372036854775808] uint64" ); } + +// --------------------------------------------------------------------------- +// h5py writes datasets indexed by a version-2 B-tree -> clawhdf5 reads +// --------------------------------------------------------------------------- + +/// With `libver='latest'`, a chunked dataset with two or more unlimited +/// dimensions indexes its chunks with a version-2 B-tree (layout v4, index +/// type 5). These used to fail with "unsupported chunked layout". +#[test] +fn h5py_btree_v2_chunk_index_clawhdf5_reads() { + skip_if_no_python!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("bt2.h5"); + let path_str = path.display().to_string(); + let script = format!( + r#" +import h5py, numpy as np +with h5py.File("{path_str}", "w", libver="latest") as f: + a = np.arange(60 * 45, dtype="> = out + .lines() + .map(|l| { + let (name, list) = l.split_once(' ').unwrap(); + (name, parse_int_list(list)) + }) + .collect(); + + let file = File::open(&path).unwrap(); + let small: Vec = (0..60 * 45).collect(); + assert_eq!(file.dataset("plain").unwrap().read_i32().unwrap(), small); + assert_eq!(file.dataset("gz").unwrap().read_i32().unwrap(), small); + let deep: Vec = (0..200 * 200).collect(); + assert_eq!(file.dataset("deep").unwrap().read_i32().unwrap(), deep); + assert_eq!( + file.dataset("sparse").unwrap().read_i32().unwrap(), + expected["sparse"] + ); + // Partial read through the same index: rows 37,50,..,128 x cols 5,36,..,160. + let slab = clawhdf5_format::selection::Selection::Hyperslab { + start: vec![37, 5], + stride: vec![13, 31], + count: vec![8, 6], + block: vec![1, 1], + }; + assert_eq!( + file.dataset("deep") + .unwrap() + .read_i32_selection(&slab) + .unwrap(), + expected["slab"] + ); +} diff --git a/docs/known-issues.md b/docs/known-issues.md index 5eecaac..d3c8993 100644 --- a/docs/known-issues.md +++ b/docs/known-issues.md @@ -116,16 +116,17 @@ not reported). ## B-tree v2 chunk index (layout v4, index type 5) is not supported -**Status:** open. +**Status:** fixed 2026-09-19. **Summary:** a chunked dataset with **two or more unlimited dimensions** written -with `libver='latest'` indexes its chunks with a version-2 B-tree. Reading it -fails with `ChunkedReadError("unsupported chunked layout version=4, -index_type=Some(5)")`. Single-chunk, implicit, fixed-array and -extensible-array indexes (and the v3 B-tree v1) are supported. +with `libver='latest'` indexes its chunks with a version-2 B-tree, and reading it +failed with `unsupported chunked layout version=4, index_type=Some(5)`. -**Repro:** `f.create_dataset("d", shape=(5, 7), chunks=(2, 3), maxshape=(None, None))` -with `h5py.File(..., libver='latest')`. +**Fix:** record types 10 (unfiltered) and 11 (filtered) are decoded — address, +stored size, filter mask, scaled offsets — through the shared chunk-listing +function, so full reads, cached reads, partial reads and fill-value handling +all work. Covered by an h5py interop test (plain, gzip+shuffle, a 2500-chunk +tree with internal nodes, a sparse dataset with a fill value, a hyperslab). ## External links and external raw data are not followed From 0addf328bc43130752c799b1e498453f07d67bdf Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 19 Sep 2026 14:05:15 -0700 Subject: [PATCH 3/7] perf(format): parallel cached decode and fewer copies on full reads MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Same-moment A/B on a 64 MB f64 dataset: chunked+deflate 110 -> 69 ms, chunked 72 -> 60 ms, contiguous 56 -> 30 ms. - read_chunked_data_cached — the path the facade uses — decompressed chunks one at a time; only the uncached reader was parallel. Cache misses are now decoded in bounded batches (128), in parallel with the `parallel` feature. - Every chunk was pushed into the 16 MiB chunk cache, which a larger dataset just churns (insert, evict moments later). Chunks are cached only when the whole dataset fits (new ChunkCache::max_bytes). - Unfiltered chunks went file -> Vec -> aligned cache buffer -> output. They are copied straight from the file bytes. - The facade's typed reads convert a contiguous dataset straight from the borrowed file bytes instead of copying it into a Vec first. - The native little-endian fast paths allocated vec![0; n] and then overwrote it; they now fill an uninitialised buffer in one copy (native_le_to_vec). alloc_output requests zeroed memory from the allocator instead of reserving and filling. The unit test that expected unfiltered chunks to land in the decompressed cache now asserts the new design (index reused, cache not involved). Co-Authored-By: Claude Fable 5.1 --- BENCHMARKS.md | 21 +++ CHANGELOG.md | 7 + crates/clawhdf5-format/src/chunk_cache.rs | 5 + crates/clawhdf5-format/src/chunked_read.rs | 155 ++++++++++++++------- crates/clawhdf5-format/src/data_read.rs | 54 +++---- crates/clawhdf5/src/reader.rs | 35 ++++- 6 files changed, 197 insertions(+), 80 deletions(-) diff --git a/BENCHMARKS.md b/BENCHMARKS.md index cc118cb..44cb0f7 100644 --- a/BENCHMARKS.md +++ b/BENCHMARKS.md @@ -91,6 +91,27 @@ because the machine's speed drifted; compare the *vs full read* column.) | contiguous | one row | 0.02 MB | 0.03 | 576 | 0.000x | | contiguous | one column | 0.03 MB | 2.55 | 12 | 0.042x | +### After: parallel cached decode, fewer copies (full reads) + +Full-read times, old and new binaries run alternately at the same moment (this +machine's absolute speed drifts over a long session, so only same-moment +comparisons mean anything): + +| layout (64 MB `f64`) | before | after | +|---|---:|---:| +| chunked + deflate | 110 ms | 69 ms | +| chunked | 72 ms | 60 ms | +| contiguous | 56 ms | 30 ms | + +What changed: the facade's cached read path decompressed chunks one at a time +(only the uncached reader was parallel) and pushed every chunk through a 16 MiB +cache that a 64 MB read simply churns; it now decodes cache misses in parallel +batches and caches only datasets that fit. Unfiltered chunks are copied +straight from the file bytes instead of via two intermediate buffers. A +contiguous dataset is converted straight from the file bytes (one copy instead +of two), and the native-endian conversions no longer zero a buffer they are +about to overwrite. + ## Search harness baseline (v2.3.0) Produced by `cargo run --release -p clawhdf5-bench --bin search_harness -- --full` diff --git a/CHANGELOG.md b/CHANGELOG.md index c1eed95..b600b8b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,13 @@ 2.7 ms, one column 5.2 ms. Results are identical to the full-read path (equivalence-tested over random hyperslabs and point lists, ranks 1-3, contiguous / chunked / deflate). New `read_harness` bench binary. +- **Faster full reads** (same-moment A/B, 64 MB `f64`): chunked + deflate + 110 -> 69 ms, chunked 72 -> 60 ms, contiguous 56 -> 30 ms. The facade's + cached read path now decompresses cache misses in parallel batches (it was + sequential; only the uncached reader was parallel) and caches only datasets + that fit the chunk cache; unfiltered chunks are copied straight from the file + bytes; a contiguous dataset is converted straight from the file bytes; and + the native-endian conversions no longer zero a buffer before overwriting it. - **Datasets indexed by a version-2 B-tree now read** (layout v4, chunk index type 5 — what `libver='latest'` uses for two or more unlimited dimensions; previously "unsupported chunked layout"). The four copies of the chunk-index diff --git a/crates/clawhdf5-format/src/chunk_cache.rs b/crates/clawhdf5-format/src/chunk_cache.rs index 369c8ac..89703f4 100644 --- a/crates/clawhdf5-format/src/chunk_cache.rs +++ b/crates/clawhdf5-format/src/chunk_cache.rs @@ -374,6 +374,11 @@ impl ChunkCache { // ----- Index operations ----- + /// The most decompressed bytes this cache will hold. + pub fn max_bytes(&self) -> usize { + self.inner.lock().map(|g| g.max_bytes).unwrap_or(0) + } + /// Bind the cache to the dataset at chunk-index address `addr`. /// /// The cache is shared per file across all of its datasets. If the cache diff --git a/crates/clawhdf5-format/src/chunked_read.rs b/crates/clawhdf5-format/src/chunked_read.rs index 0dbf6d4..8dade54 100644 --- a/crates/clawhdf5-format/src/chunked_read.rs +++ b/crates/clawhdf5-format/src/chunked_read.rs @@ -165,12 +165,29 @@ pub(crate) fn checked_chunk_byte_len( /// process when the allocation fails; a size taken from the file must surface /// as an error instead. pub(crate) fn alloc_output(len: usize) -> Result, FormatError> { - let mut out = Vec::new(); - out.try_reserve_exact(len).map_err(|_| { - FormatError::Overflow(format!("cannot allocate {len} bytes for dataset output")) - })?; - out.resize(len, 0); - Ok(out) + if len == 0 { + return Ok(Vec::new()); + } + let failed = + || FormatError::Overflow(format!("cannot allocate {len} bytes for dataset output")); + let layout = core::alloc::Layout::array::(len).map_err(|_| failed())?; + // Ask the allocator for zeroed memory instead of reserving and then + // writing zeros: for a large buffer the OS hands out already-zero pages + // lazily, where an explicit fill touches every page up front — and most of + // the buffer is about to be overwritten with chunk data anyway. + // + // SAFETY (both arms): `layout` has non-zero size (len > 0) and alignment 1. + #[cfg(feature = "std")] + let ptr = unsafe { std::alloc::alloc_zeroed(layout) }; + #[cfg(not(feature = "std"))] + let ptr = unsafe { alloc::alloc::alloc_zeroed(layout) }; + if ptr.is_null() { + return Err(failed()); + } + // SAFETY: `ptr` came from the global allocator with the layout of + // `[u8; len]`, which is exactly what `Vec` with capacity `len` frees; + // all `len` bytes are initialised (zero). + Ok(unsafe { Vec::from_raw_parts(ptr, len, len) }) } fn read_offset(data: &[u8], pos: usize, size: u8) -> Result { @@ -362,6 +379,10 @@ pub fn generate_implicit_chunks( } /// Read a chunked dataset, decompressing chunks as needed. +/// Chunks decompressed together before being copied out, bounding the extra +/// memory a parallel full read holds at once. +const DECODE_BATCH: usize = 128; + /// B-tree v2 record types used for chunk indexing. const BT2_CHUNK_UNFILTERED: u8 = 10; const BT2_CHUNK_FILTERED: u8 = 11; @@ -823,52 +844,86 @@ pub fn read_chunked_data_cached( let chunk_total_bytes = checked_chunk_byte_len(&chunk_dims, elem_size)?; - for chunk_info in &chunks { - let coord: Vec = chunk_info.offsets.iter().take(rank).copied().collect(); - - // Try decompressed cache first - let decompressed = if let Some(cached) = cache.get_decompressed_aligned(&coord) { - cached - } else { - // Decompress from file - let c_addr = chunk_info.address as usize; - let size = chunk_info.chunk_size as usize; - ensure_len(file_data, c_addr, size)?; - let raw_chunk = &file_data[c_addr..c_addr + size]; - let dec = if let Some(pl) = pipeline { - if chunk_info.filter_mask == 0 { - decompress_chunk(raw_chunk, pl, chunk_total_bytes, elem_size as u32)? - } else { - raw_chunk.to_vec() - } - } else { - raw_chunk.to_vec() - }; - cache.put_decompressed(coord, dec) - }; - + let mut place = |data: &[u8], chunk_info: &ChunkInfo| { + if rank == 0 { + let copy_len = data.len().min(output.len()); + output[..copy_len].copy_from_slice(&data[..copy_len]); + return; + } let chunk_offsets: Vec = chunk_info .offsets .iter() .take(rank) .map(|&o| o as usize) .collect(); + copy_chunk_to_output( + data, + &mut output, + &chunk_offsets, + &chunk_dims, + &ds_dims, + &ds_strides, + &chunk_strides, + elem_size, + rank, + ); + }; + let raw_bytes = |chunk_info: &ChunkInfo| -> Result<&[u8], FormatError> { + let c_addr = chunk_info.address as usize; + let size = chunk_info.chunk_size as usize; + ensure_len(file_data, c_addr, size)?; + Ok(&file_data[c_addr..c_addr + size]) + }; - if rank == 0 { - let copy_len = decompressed.len().min(output.len()); - output[..copy_len].copy_from_slice(&decompressed[..copy_len]); - } else { - copy_chunk_to_output( - &decompressed, - &mut output, - &chunk_offsets, - &chunk_dims, - &ds_dims, - &ds_strides, - &chunk_strides, - elem_size, - rank, - ); + // Chunks stored as-is (no pipeline, or the filter mask says this chunk + // skipped it) are copied straight from the file bytes: they are already in + // memory, so routing them through a Vec and then an aligned cache buffer + // was two extra copies of the whole dataset for nothing. + let stored_raw = |c: &ChunkInfo| pipeline.is_none() || c.filter_mask != 0; + let mut misses: Vec<&ChunkInfo> = Vec::new(); + for chunk_info in &chunks { + if stored_raw(chunk_info) { + place(raw_bytes(chunk_info)?, chunk_info); + continue; + } + let coord: Vec = chunk_info.offsets.iter().take(rank).copied().collect(); + match cache.get_decompressed_aligned(&coord) { + Some(cached) => place(&cached, chunk_info), + None => misses.push(chunk_info), + } + } + + // Decompress what the cache didn't have, a bounded batch at a time — in + // parallel with the `parallel` feature (this path, the one the facade + // uses, was sequential; only the uncached reader was parallel). Chunks are + // cached only when the whole dataset fits: pushing a larger dataset + // through the cache just evicts each chunk moments after inserting it. + let cache_them = total_bytes <= cache.max_bytes(); + if let Some(pl) = pipeline { + let decode = |c: &&ChunkInfo| -> Result, FormatError> { + decompress_chunk(raw_bytes(c)?, pl, chunk_total_bytes, elem_size as u32) + }; + for batch in misses.chunks(DECODE_BATCH) { + #[cfg(feature = "parallel")] + let decoded: Vec, FormatError>> = if batch.len() >= 4 { + use rayon::prelude::*; + batch.par_iter().map(decode).collect() + } else { + batch.iter().map(decode).collect() + }; + #[cfg(not(feature = "parallel"))] + let decoded: Vec, FormatError>> = batch.iter().map(decode).collect(); + + for (chunk_info, data) in batch.iter().zip(decoded) { + let data = data?; + if cache_them { + let coord: Vec = chunk_info.offsets.iter().take(rank).copied().collect(); + let cached = cache.put_decompressed(coord, data); + place(&cached, chunk_info); + } else { + place(&data, chunk_info); + } + } } } @@ -2186,21 +2241,23 @@ mod tests { } #[test] - fn cached_read_second_call_uses_cache() { + fn cached_read_second_call_reuses_the_index() { let values: Vec = (0..20).map(|i| i as f64).collect(); let (file_data, layout, dataspace) = build_1d_chunked_file(&values, 10); let datatype = make_f64_type(); let cache = ChunkCache::new(); - // First read — populates index + decompressed cache + // First read — populates the chunk index. These chunks are stored + // unfiltered, so they are copied straight from the file bytes and the + // decompressed-chunk cache is (deliberately) not involved. let raw1 = read_chunked_data_cached( &file_data, &layout, &dataspace, &datatype, None, 8, 8, &cache, ) .unwrap(); assert!(cache.has_index()); - assert!(cache.cached_chunk_count() > 0); + assert_eq!(cache.cached_chunk_count(), 0); - // Second read — should hit the decompressed cache + // Second read — reuses the cached index let raw2 = read_chunked_data_cached( &file_data, &layout, &dataspace, &datatype, None, 8, 8, &cache, ) diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index 13b4953..4aa83a4 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -876,6 +876,30 @@ fn get_size(dt: &Datatype) -> usize { dt.type_size() as usize } +/// Reinterpret little-endian bytes as `count` native values of `T` on a +/// little-endian target, in one copy. +/// +/// The buffer is allocated uninitialised and filled by the copy. It used to be +/// `vec![0; count]` first, which for a large dataset meant writing every page +/// twice (zero it, then overwrite it) — about as expensive as the copy itself. +#[cfg(target_endian = "little")] +fn native_le_to_vec(raw: &[u8], count: usize) -> Vec { + let bytes = count * core::mem::size_of::(); + debug_assert!(bytes <= raw.len()); + let mut result: Vec = Vec::with_capacity(count); + // SAFETY: `result` has capacity for `count` values of `T`, i.e. `bytes` + // bytes; `raw` holds at least `bytes` bytes (callers derive `count` from + // `raw.len() / size_of::()`); the regions cannot overlap because + // `result` was just allocated. Every `T` used here (f32/f64/i32/i64) is + // valid for any bit pattern, so after the copy all `count` values are + // initialised and `set_len` is sound. + unsafe { + core::ptr::copy_nonoverlapping(raw.as_ptr(), result.as_mut_ptr().cast::(), bytes); + result.set_len(count); + } + result +} + /// Convert raw bytes to `f64` values. pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result, FormatError> { // Array datatypes (e.g. an array-typed compound member) are read as a flat @@ -903,14 +927,7 @@ pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result, FormatEr .. } ) { - let mut result = vec![0.0f64; count]; - // SAFETY: On LE platforms, f64 in-memory representation matches LE bytes. - // We copy raw bytes directly into the f64 buffer. - // SAFETY: The byte slice is properly aligned for this type and the length is divisible by size_of::(). - unsafe { - core::ptr::copy_nonoverlapping(raw.as_ptr(), result.as_mut_ptr() as *mut u8, raw.len()); - } - return Ok(result); + return Ok(native_le_to_vec::(raw, count)); } let order = get_byte_order(datatype); @@ -993,12 +1010,7 @@ pub fn read_as_i64(raw: &[u8], datatype: &Datatype) -> Result, FormatEr } ) { - let mut result = vec![0i64; count]; - // SAFETY: The byte slice is properly aligned for this type and the length is divisible by size_of::(). - unsafe { - core::ptr::copy_nonoverlapping(raw.as_ptr(), result.as_mut_ptr() as *mut u8, raw.len()); - } - return Ok(result); + return Ok(native_le_to_vec::(raw, count)); } let order = get_byte_order(datatype); @@ -1062,12 +1074,7 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result, FormatEr .. } ) { - let mut result = vec![0.0f32; count]; - // SAFETY: The byte slice is properly aligned for this type and the length is divisible by size_of::(). - unsafe { - core::ptr::copy_nonoverlapping(raw.as_ptr(), result.as_mut_ptr() as *mut u8, raw.len()); - } - return Ok(result); + return Ok(native_le_to_vec::(raw, count)); } let order = get_byte_order(datatype); @@ -1144,12 +1151,7 @@ pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result, FormatEr } ) { - let mut result = vec![0i32; count]; - // SAFETY: The byte slice is properly aligned for this type and the length is divisible by size_of::(). - unsafe { - core::ptr::copy_nonoverlapping(raw.as_ptr(), result.as_mut_ptr() as *mut u8, raw.len()); - } - return Ok(result); + return Ok(native_le_to_vec::(raw, count)); } let order = get_byte_order(datatype); diff --git a/crates/clawhdf5/src/reader.rs b/crates/clawhdf5/src/reader.rs index d6d6032..64f96ef 100644 --- a/crates/clawhdf5/src/reader.rs +++ b/crates/clawhdf5/src/reader.rs @@ -381,8 +381,13 @@ impl<'f> Dataset<'f> { /// Read all data as `f64` values. pub fn read_f64(&self) -> Result, Error> { - let raw = self.read_raw()?; let dt = self.datatype()?; + // A contiguous dataset is converted straight from the file bytes; going + // through `read_raw` first copied the whole dataset an extra time. + if let Ok(Some(bytes)) = self.read_raw_ref() { + return Ok(data_read::read_as_f64(bytes, &dt)?); + } + let raw = self.read_raw()?; Ok(data_read::read_as_f64(&raw, &dt)?) } @@ -393,29 +398,49 @@ impl<'f> Dataset<'f> { /// /// Read all data as `f32` values. pub fn read_f32(&self) -> Result, Error> { - let raw = self.read_raw()?; let dt = self.datatype()?; + // A contiguous dataset is converted straight from the file bytes; going + // through `read_raw` first copied the whole dataset an extra time. + if let Ok(Some(bytes)) = self.read_raw_ref() { + return Ok(data_read::read_as_f32(bytes, &dt)?); + } + let raw = self.read_raw()?; Ok(data_read::read_as_f32(&raw, &dt)?) } /// Read all data as `i32` values. pub fn read_i32(&self) -> Result, Error> { - let raw = self.read_raw()?; let dt = self.datatype()?; + // A contiguous dataset is converted straight from the file bytes; going + // through `read_raw` first copied the whole dataset an extra time. + if let Ok(Some(bytes)) = self.read_raw_ref() { + return Ok(data_read::read_as_i32(bytes, &dt)?); + } + let raw = self.read_raw()?; Ok(data_read::read_as_i32(&raw, &dt)?) } /// Read all data as `i64` values. pub fn read_i64(&self) -> Result, Error> { - let raw = self.read_raw()?; let dt = self.datatype()?; + // A contiguous dataset is converted straight from the file bytes; going + // through `read_raw` first copied the whole dataset an extra time. + if let Ok(Some(bytes)) = self.read_raw_ref() { + return Ok(data_read::read_as_i64(bytes, &dt)?); + } + let raw = self.read_raw()?; Ok(data_read::read_as_i64(&raw, &dt)?) } /// Read all data as `u64` values. pub fn read_u64(&self) -> Result, Error> { - let raw = self.read_raw()?; let dt = self.datatype()?; + // A contiguous dataset is converted straight from the file bytes; going + // through `read_raw` first copied the whole dataset an extra time. + if let Ok(Some(bytes)) = self.read_raw_ref() { + return Ok(data_read::read_as_u64(bytes, &dt)?); + } + let raw = self.read_raw()?; Ok(data_read::read_as_u64(&raw, &dt)?) } From f3d63dbdcd268b8ceb9afd3c7344b3f755472326 Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 19 Sep 2026 14:06:10 -0700 Subject: [PATCH 4/7] perf(format): un-shuffle by interleaving fixed-width byte planes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit shuffle_decompress — on the read path of every compressed dataset, since shuffle is applied automatically before compression — was the naive `result[i * es + j] = data[j * n + i]`: a multiply and two bounds checks per byte. It now interleaves fixed-width arrays of byte planes for element sizes 2/4/8/16 (bounds checks hoisted, vectorisable), with a chunked generic fallback. The write-side shuffle was already optimised; this was the asymmetry the survey flagged. Modest wall-clock effect now that decode is parallel (chunked+deflate full read ~70 -> ~66 ms). Round-trip test over element sizes 1-24 and several lengths. Co-Authored-By: Claude Fable 5.1 --- crates/clawhdf5-format/src/filters.rs | 46 +++++++++++++++++++++++++-- 1 file changed, 43 insertions(+), 3 deletions(-) diff --git a/crates/clawhdf5-format/src/filters.rs b/crates/clawhdf5-format/src/filters.rs index 2e99d70..0a0225f 100644 --- a/crates/clawhdf5-format/src/filters.rs +++ b/crates/clawhdf5-format/src/filters.rs @@ -845,9 +845,33 @@ fn shuffle_decompress(data: &[u8], element_size: usize) -> Result, Forma let num_elements = data.len() / element_size; let mut result = vec![0u8; data.len()]; - for i in 0..num_elements { - for j in 0..element_size { - result[i * element_size + j] = data[j * num_elements + i]; + // The shuffled stream is `element_size` byte planes of `num_elements` + // bytes each; un-shuffling interleaves them. This is on the read path of + // every compressed dataset (shuffle is applied automatically before + // compression). The naive `result[i * es + j] = data[j * n + i]` form does + // a multiply and two bounds checks per byte and defeats vectorisation; + // fixed-width plane arrays sliced to a common length let the compiler + // hoist the checks and emit interleaves for the common 4- and 8-byte + // element sizes. + fn interleave(data: &[u8], n: usize, out: &mut [u8]) { + let planes: [&[u8]; W] = core::array::from_fn(|j| &data[j * n..(j + 1) * n]); + for (i, element) in out.chunks_exact_mut(W).enumerate() { + for j in 0..W { + element[j] = planes[j][i]; + } + } + } + match element_size { + 2 => interleave::<2>(data, num_elements, &mut result), + 4 => interleave::<4>(data, num_elements, &mut result), + 8 => interleave::<8>(data, num_elements, &mut result), + 16 => interleave::<16>(data, num_elements, &mut result), + _ => { + for (i, element) in result.chunks_exact_mut(element_size).enumerate() { + for (j, byte) in element.iter_mut().enumerate() { + *byte = data[j * num_elements + i]; + } + } } } @@ -1848,4 +1872,20 @@ mod tests { }; assert!(decompress_chunk(&data, &pipeline, 16, 1).is_err()); } + #[test] + fn unshuffle_inverts_shuffle_for_every_element_size() { + for element_size in [1usize, 2, 3, 4, 5, 8, 12, 16, 24] { + for elements in [0usize, 1, 2, 7, 64, 1000] { + let original: Vec = (0..element_size * elements) + .map(|i| (i * 31 + 7) as u8) + .collect(); + let shuffled = shuffle_compress(&original, element_size).unwrap(); + assert_eq!( + shuffle_decompress(&shuffled, element_size).unwrap(), + original, + "element_size {element_size}, {elements} elements" + ); + } + } + } } From b36c6ec2af76ad172145a43d26d00ce1bb4e3e26 Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 19 Sep 2026 14:06:36 -0700 Subject: [PATCH 5/7] style(format): as_chunks_mut in the un-shuffle interleave (clippy) Co-Authored-By: Claude Fable 5.1 --- crates/clawhdf5-format/src/filters.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/crates/clawhdf5-format/src/filters.rs b/crates/clawhdf5-format/src/filters.rs index 0a0225f..9979d3f 100644 --- a/crates/clawhdf5-format/src/filters.rs +++ b/crates/clawhdf5-format/src/filters.rs @@ -855,9 +855,9 @@ fn shuffle_decompress(data: &[u8], element_size: usize) -> Result, Forma // element sizes. fn interleave(data: &[u8], n: usize, out: &mut [u8]) { let planes: [&[u8]; W] = core::array::from_fn(|j| &data[j * n..(j + 1) * n]); - for (i, element) in out.chunks_exact_mut(W).enumerate() { - for j in 0..W { - element[j] = planes[j][i]; + for (i, element) in out.as_chunks_mut::().0.iter_mut().enumerate() { + for (byte, plane) in element.iter_mut().zip(&planes) { + *byte = plane[i]; } } } From 05c665a898c084481cd88274ee9c6c84bc72a45a Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 19 Sep 2026 14:21:04 -0700 Subject: [PATCH 6/7] feat(format): choose chunk dimensions automatically for large datasets MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Requesting a filter without chunk dimensions made the whole dataset a single chunk. Any read, even one row, then decompresses everything, and a large dataset cannot be decoded in parallel — which also made the new partial reads pointless for such files. auto_chunk_dims keeps datasets up to 1 MiB as one chunk (unchanged behaviour) and splits larger ones by halving the dimensions in turn, so chunks keep roughly the dataset's proportions, until a chunk is at most 1 MiB — h5py's approach. An empty (unlimited, unwritten) dimension is treated as 1024. The writer passes the element size through resolve_chunk_dims_for; the old resolve_chunk_dims assumes 8-byte elements. Explicit with_chunks always wins. Interop test: h5py reads an auto-chunked 13 MB deflate dataset, sees chunks between 128 KiB and 1 MiB, and a small dataset still has one chunk. Co-Authored-By: Claude Fable 5.1 --- CHANGELOG.md | 8 ++ crates/clawhdf5-format/src/chunked_write.rs | 87 +++++++++++++++++++-- crates/clawhdf5-format/src/file_writer.rs | 4 +- crates/clawhdf5/tests/h5py_interop_tests.rs | 59 ++++++++++++++ 4 files changed, 152 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b600b8b..8e8ee83 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,14 @@ type 5 — what `libver='latest'` uses for two or more unlimited dimensions; previously "unsupported chunked layout"). The four copies of the chunk-index dispatch are now one shared function, so every read path gets it. +- **Automatic chunk sizes.** Asking for compression (or any filter) without + `with_chunks` used to store the whole dataset as one chunk, so any read had + to decompress everything and nothing could be decoded in parallel. Datasets up + to 1 MiB stay a single chunk, as before; larger ones are split by halving the + dimensions in turn until a chunk is at most 1 MiB (the approach h5py takes). + **Behaviour change:** large compressed datasets written without explicit + chunk dimensions get a different (standard, h5py-readable) layout. Explicit + `with_chunks` is unaffected. - **Out-of-range selections are errors.** They used to return data: a hyperslab past an edge came back padded with zeros, and a point whose column was out of range wrapped into the next row and returned that element. Now diff --git a/crates/clawhdf5-format/src/chunked_write.rs b/crates/clawhdf5-format/src/chunked_write.rs index 405312a..e6a9e44 100644 --- a/crates/clawhdf5-format/src/chunked_write.rs +++ b/crates/clawhdf5-format/src/chunked_write.rs @@ -49,6 +49,38 @@ pub struct ChunkOptions { pub pcodec: bool, } +/// Largest chunk the automatic choice produces, in bytes. +const AUTO_CHUNK_TARGET_BYTES: u64 = 1 << 20; + +/// Extent assumed for a dimension that is currently empty (an unlimited +/// dimension not yet written to) — the same stand-in h5py uses. +const AUTO_CHUNK_EMPTY_DIM: u64 = 1024; + +/// Choose chunk dimensions for a dataset nobody specified them for. +/// +/// Asking for compression (or any filter) without chunk dimensions used to +/// make the whole dataset one chunk. That defeats the point of chunking: any +/// read — even a single row — must decompress everything, and a large dataset +/// cannot be decompressed in parallel. Datasets up to the target size stay a +/// single chunk, exactly as before; larger ones are split by halving the +/// dimensions in turn (so chunks keep roughly the dataset's proportions, the +/// approach h5py takes) until a chunk fits the target. +pub fn auto_chunk_dims(shape: &[u64], elem_size: usize) -> Vec { + let mut dims: Vec = shape + .iter() + .map(|&d| if d == 0 { AUTO_CHUNK_EMPTY_DIM } else { d }) + .collect(); + let elem = elem_size.max(1) as u64; + let bytes = |dims: &[u64]| dims.iter().fold(elem, |acc, &d| acc.saturating_mul(d)); + let mut axis = 0; + while bytes(&dims) > AUTO_CHUNK_TARGET_BYTES && dims.iter().any(|&d| d > 1) { + let i = axis % dims.len(); + dims[i] = dims[i].div_ceil(2); + axis += 1; + } + dims +} + impl ChunkOptions { /// Whether any chunking option is enabled. pub fn is_chunked(&self) -> bool { @@ -135,11 +167,17 @@ impl ChunkOptions { /// Determine chunk dimensions, using user-specified or auto-computing. pub fn resolve_chunk_dims(&self, shape: &[u64]) -> Vec { - if let Some(ref dims) = self.chunk_dims { - dims.clone() - } else { - // Auto chunk: use the full dataset shape (single chunk) - shape.to_vec() + // Without the element size, assume 8 bytes (the widest common scalar); + // the writer uses `resolve_chunk_dims_for`. + self.resolve_chunk_dims_for(shape, 8) + } + + /// Chunk dimensions for a dataset of `shape` whose elements are `elem_size` + /// bytes: the caller's if given, otherwise chosen automatically. + pub fn resolve_chunk_dims_for(&self, shape: &[u64], elem_size: usize) -> Vec { + match self.chunk_dims { + Some(ref dims) => dims.clone(), + None => auto_chunk_dims(shape, elem_size), } } } @@ -1143,6 +1181,45 @@ mod tests { assert_eq!(dims, vec![100, 50]); } + #[test] + fn auto_chunking_splits_only_large_datasets() { + let bytes = |dims: &[u64], elem: u64| dims.iter().product::() * elem; + // Up to the target: one chunk, as before. + assert_eq!(auto_chunk_dims(&[100, 50], 8), [100, 50]); + assert_eq!(auto_chunk_dims(&[131_072], 8), [131_072]); // exactly 1 MiB + // Larger: split, keeping proportions, never above the target. + let big = auto_chunk_dims(&[4096, 2048], 8); + assert!(bytes(&big, 8) <= AUTO_CHUNK_TARGET_BYTES, "{big:?}"); + assert!(bytes(&big, 8) > AUTO_CHUNK_TARGET_BYTES / 4, "{big:?}"); + assert_eq!(big[0] / big[1], 2, "proportions kept: {big:?}"); + // Every dimension stays within the dataset and at least 1. + for shape in [ + vec![10_000_000u64], + vec![3, 5_000_000], + vec![1, 1, 9_000_000], + vec![7; 9], + ] { + let dims = auto_chunk_dims(&shape, 4); + assert!( + dims.iter().zip(&shape).all(|(c, s)| *c >= 1 && c <= s), + "{shape:?} -> {dims:?}" + ); + assert!( + bytes(&dims, 4) <= AUTO_CHUNK_TARGET_BYTES, + "{shape:?} -> {dims:?}" + ); + } + // An empty (unlimited, unwritten) dimension still gets a usable chunk. + let growable = auto_chunk_dims(&[0, 128], 8); + assert!(growable[0] >= 1 && bytes(&growable, 8) <= AUTO_CHUNK_TARGET_BYTES); + // Explicit dimensions always win. + let explicit = ChunkOptions { + chunk_dims: Some(vec![10, 10]), + ..Default::default() + }; + assert_eq!(explicit.resolve_chunk_dims_for(&[4096, 2048], 8), [10, 10]); + } + #[test] fn chunk_options_pipeline_deflate() { // Auto-shuffle is applied before compression by default (matches h5py). diff --git a/crates/clawhdf5-format/src/file_writer.rs b/crates/clawhdf5-format/src/file_writer.rs index 644c483..c57691b 100644 --- a/crates/clawhdf5-format/src/file_writer.rs +++ b/crates/clawhdf5-format/src/file_writer.rs @@ -1221,8 +1221,10 @@ impl FileWriter { precompressed: None, }); } else if is_chunked[i] { - let chunk_dims = d.chunk_options.resolve_chunk_dims(&d.ds.dimensions); let elem_size = d.dt.type_size() as usize; + let chunk_dims = d + .chunk_options + .resolve_chunk_dims_for(&d.ds.dimensions, elem_size); // Compress once in Pass 1; cache the result so Pass 2 can skip // re-compression and just rebuild the index with real addresses. let pre = precompress_chunks( diff --git a/crates/clawhdf5/tests/h5py_interop_tests.rs b/crates/clawhdf5/tests/h5py_interop_tests.rs index b5692d3..e177b72 100644 --- a/crates/clawhdf5/tests/h5py_interop_tests.rs +++ b/crates/clawhdf5/tests/h5py_interop_tests.rs @@ -979,3 +979,62 @@ with h5py.File("{path_str}", "r") as f: expected["slab"] ); } + +// --------------------------------------------------------------------------- +// clawhdf5 auto-chunks a large compressed dataset -> h5py reads +// --------------------------------------------------------------------------- + +/// Compression without explicit chunk dimensions used to store the whole +/// dataset as a single chunk. Large datasets are now split automatically; +/// h5py must read the result and see sensibly sized chunks. +#[test] +fn clawhdf5_auto_chunked_dataset_h5py_reads() { + skip_if_no_python!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("auto_chunk.h5"); + let path_str = path.display().to_string(); + + let (rows, cols) = (1500u64, 1100u64); // 13.2 MB of f64 + let data: Vec = (0..rows * cols).map(|i| (i % 9973) as f64 * 0.25).collect(); + let mut builder = FileBuilder::new(); + builder + .create_dataset("big") + .with_f64_data(&data) + .with_shape(&[rows, cols]) + .with_deflate(4); + builder + .create_dataset("small") + .with_f64_data(&data[..600]) + .with_shape(&[20, 30]) + .with_deflate(4); + builder.write(&path).unwrap(); + + let out = run_python_output(&format!( + r#" +import h5py, numpy as np +with h5py.File("{path_str}", "r") as f: + big, small = f["big"], f["small"] + expect = (np.arange(1500 * 1100) % 9973) * 0.25 + ok = bool(np.array_equal(big[...].ravel(), expect)) and bool(np.array_equal(small[...].ravel(), expect[:600])) + chunk_bytes = int(np.prod(big.chunks)) * 8 + print(ok, chunk_bytes <= 1 << 20, chunk_bytes >= 1 << 17, small.chunks == (20, 30), big.compression) +"# + )); + assert_eq!(out.trim(), "True True True True gzip"); + + // And it reads back here, in full and partially. + let file = File::open(&path).unwrap(); + let ds = file.dataset("big").unwrap(); + assert_eq!(ds.read_f64().unwrap(), data); + let row = clawhdf5_format::selection::Selection::Hyperslab { + start: vec![777, 0], + stride: vec![1, 1], + count: vec![1, cols], + block: vec![1, 1], + }; + let start = (777 * cols) as usize; + assert_eq!( + ds.read_f64_selection(&row).unwrap(), + data[start..start + cols as usize] + ); +} From 52cfcf20b2113a1f733e12076c827312f544e0aa Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 19 Sep 2026 14:24:04 -0700 Subject: [PATCH 7/7] feat(format): parse H5T_STD_REF references and decode object references HDF5 1.12 revised the reference datatype (class 7) in datatype message version 4: reference types 2-4 are the new H5T_STD_REF object / dataset-region / attribute references. Datatype::parse rejected them with InvalidReferenceType, so any dataset of that type was unreadable. h5py cannot write this type, which is why it had never been tested. A real file was produced by calling the libhdf5 bundled in the h5py wheel through ctypes (H5T_STD_REF_g, H5Rcreate_object, H5Dwrite); the 2 KB result is committed as tests/fixtures/std_ref_hdf5_2_0.h5 with its generator, gen_std_ref.py. - ReferenceType gains Object2, DatasetRegion2 and Attribute, accepted only from datatype version 4. - read_object_references decodes Object2 elements: type(1) flags(1) token_size(1) token, zero-padded to the element size; the token is the target's object header address. A null reference decodes to the undefined address; an external reference, a wrong type byte or a token that doesn't fit is an error. The fixture test follows both references and checks they resolve to the objects they were created from. Co-Authored-By: Claude Fable 5.1 --- CHANGELOG.md | 5 + crates/clawhdf5-format/src/data_read.rs | 60 +++++++++++ crates/clawhdf5-format/src/datatype.rs | 44 +++++++- .../tests/fixtures/gen_std_ref.py | 44 ++++++++ .../tests/fixtures/std_ref_hdf5_2_0.h5 | Bin 0 -> 2104 bytes .../clawhdf5-format/tests/reference_tests.rs | 94 ++++++++++++++++++ docs/known-issues.md | 20 ++-- 7 files changed, 258 insertions(+), 9 deletions(-) create mode 100644 crates/clawhdf5-format/tests/fixtures/gen_std_ref.py create mode 100644 crates/clawhdf5-format/tests/fixtures/std_ref_hdf5_2_0.h5 diff --git a/CHANGELOG.md b/CHANGELOG.md index 8e8ee83..705e7ac 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,11 @@ type 5 — what `libver='latest'` uses for two or more unlimited dimensions; previously "unsupported chunked layout"). The four copies of the chunk-index dispatch are now one shared function, so every read path gets it. +- **`H5T_STD_REF` references** (HDF5 1.12+, datatype message version 4) parse: + `ReferenceType` gains `Object2`, `DatasetRegion2` and `Attribute`, and + `read_object_references` decodes the new object references. Previously any + dataset of this type failed with `InvalidReferenceType(2)`. Tested against a + file written by HDF5 2.0 itself (fixture + generator script committed). - **Automatic chunk sizes.** Asking for compression (or any filter) without `with_chunks` used to store the whole dataset as one chunk, so any read had to decompress everything and nothing could be decoded in parallel. Datasets up diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index 4aa83a4..0e6773f 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -1427,6 +1427,26 @@ pub fn read_object_references( } Ok(result) } + Datatype::Reference { + ref_type: crate::datatype::ReferenceType::Object2, + size, + } => { + let elem_size = *size as usize; + if elem_size == 0 { + return Ok(Vec::new()); + } + if !raw.len().is_multiple_of(elem_size) { + return Err(FormatError::DataSizeMismatch { + expected: 0, + actual: raw.len(), + }); + } + raw.chunks_exact(elem_size) + .map(|element| { + decode_std_object_ref(element).map(|address| ObjectReference { address }) + }) + .collect() + } _ => Err(FormatError::TypeMismatch { expected: "Reference(Object)", actual: datatype_name(datatype), @@ -1434,6 +1454,46 @@ pub fn read_object_references( } } +/// Decode one `H5T_STD_REF` object reference as stored in a dataset: +/// `type(1) flags(1) token_size(1) token(token_size)`, zero-padded to the +/// element size. For a reference within the same file the token is the target +/// object's header address. An all-zero element is a null reference and +/// decodes to the undefined address (`u64::MAX`). +fn decode_std_object_ref(element: &[u8]) -> Result { + const STD_REF_OBJECT: u8 = 2; + const FLAG_EXTERNAL: u8 = 0x01; + if element.iter().all(|&b| b == 0) { + return Ok(u64::MAX); + } + let [ref_type, flags, token_size, token @ ..] = element else { + return Err(FormatError::UnexpectedEof { + expected: 3, + available: element.len(), + }); + }; + if *ref_type != STD_REF_OBJECT { + return Err(FormatError::InvalidReferenceType(*ref_type)); + } + if flags & FLAG_EXTERNAL != 0 { + // Carries a file name as well; nothing here follows those. + return Err(FormatError::TypeMismatch { + expected: "object reference within this file", + actual: "external object reference", + }); + } + let n = *token_size as usize; + if n == 0 || n > 8 || n > token.len() { + return Err(FormatError::UnexpectedEof { + expected: 3 + n, + available: element.len(), + }); + } + Ok(token[..n] + .iter() + .rev() + .fold(0u64, |addr, &byte| (addr << 8) | u64::from(byte))) +} + /// Read region references from raw bytes. /// /// Region references encode a dataset selection (hyperslab, point list, etc.) diff --git a/crates/clawhdf5-format/src/datatype.rs b/crates/clawhdf5-format/src/datatype.rs index 19242c0..2a83e0b 100644 --- a/crates/clawhdf5-format/src/datatype.rs +++ b/crates/clawhdf5-format/src/datatype.rs @@ -36,8 +36,18 @@ pub enum CharacterSet { /// Reference type. #[derive(Debug, Clone, PartialEq)] pub enum ReferenceType { + /// Legacy object reference: the target's object header address. Object, + /// Legacy dataset region reference. DatasetRegion, + /// `H5T_STD_REF` object reference (HDF5 1.12+, datatype message version + /// 4): a small header followed by an object token. Decoded by + /// `data_read::read_object_references`. + Object2, + /// `H5T_STD_REF` dataset region reference. + DatasetRegion2, + /// `H5T_STD_REF` attribute reference. + Attribute, } /// A member of a compound datatype. @@ -424,9 +434,15 @@ impl Datatype { 7 => { // Reference let ref_type_val = bf0 & 0x0F; - let ref_type = match ref_type_val { - 0 => ReferenceType::Object, - 1 => ReferenceType::DatasetRegion, + // Datatype message version 4 (HDF5 1.12) revised this class: + // types 2-4 are the new `H5T_STD_REF` references, and the high + // nibble of the first flag byte carries their encoding version. + let ref_type = match (ref_type_val, version) { + (0, _) => ReferenceType::Object, + (1, _) => ReferenceType::DatasetRegion, + (2, 4..) => ReferenceType::Object2, + (3, 4..) => ReferenceType::DatasetRegion2, + (4, 4..) => ReferenceType::Attribute, _ => return Err(FormatError::InvalidReferenceType(ref_type_val)), }; Ok((Datatype::Reference { size, ref_type }, pos)) @@ -1563,6 +1579,28 @@ mod tests { assert_eq!(err, FormatError::InvalidCharacterSet(2)); } + #[test] + fn test_reference_v4_std_ref_from_hdf5_2_0() { + // Datatype message of an H5T_STD_REF dataset written by HDF5 2.0: + // class 7, version 4, type 2 (object), encoding version 1, 18 bytes. + let bytes = [0x47, 0x12, 0x00, 0x00, 0x12, 0x00, 0x00, 0x00]; + let (dt, consumed) = Datatype::parse(&bytes).unwrap(); + assert_eq!(consumed, 8); + assert_eq!( + dt, + Datatype::Reference { + size: 18, + ref_type: ReferenceType::Object2 + } + ); + // The new types are only valid from datatype version 4. + let old_version = [0x37, 0x12, 0x00, 0x00, 0x12, 0x00, 0x00, 0x00]; + assert_eq!( + Datatype::parse(&old_version).unwrap_err(), + FormatError::InvalidReferenceType(2) + ); + } + #[test] fn test_error_invalid_reference_type() { let buf = build_dt_header(7, 1, [5, 0, 0], 8); diff --git a/crates/clawhdf5-format/tests/fixtures/gen_std_ref.py b/crates/clawhdf5-format/tests/fixtures/gen_std_ref.py new file mode 100644 index 0000000..f73a8c5 --- /dev/null +++ b/crates/clawhdf5-format/tests/fixtures/gen_std_ref.py @@ -0,0 +1,44 @@ +"""Generate std_ref_hdf5_2_0.h5: a dataset of H5T_STD_REF (the reference +datatype introduced in HDF5 1.12, datatype message version 4) holding two +object references — to /target (a dataset) and /grp (a group). + +h5py has no API for this type, so the file is written by calling the libhdf5 +bundled in the h5py wheel directly through ctypes. Written with h5py 3.16.0 / +HDF5 2.0.0. Re-run only if the fixture ever needs regenerating: + + python gen_std_ref.py std_ref_hdf5_2_0.h5 +""" +import ctypes +import glob +import os +import sys + +import h5py +import numpy as np + +libdir = os.path.join(os.path.dirname(os.path.dirname(h5py.__file__)), "h5py.libs") +libs = [p for p in glob.glob(os.path.join(libdir, "libhdf5*.so*")) if "_hl" not in os.path.basename(p)] +lib = ctypes.CDLL(libs[0]) +lib.H5open() +hid = ctypes.c_int64 +std_ref = hid.in_dll(lib, "H5T_STD_REF_g").value + +lib.H5Screate_simple.restype = hid +lib.H5Screate_simple.argtypes = [ctypes.c_int, ctypes.POINTER(ctypes.c_uint64), ctypes.POINTER(ctypes.c_uint64)] +lib.H5Dcreate2.restype = hid +lib.H5Dcreate2.argtypes = [hid, ctypes.c_char_p, hid, hid, hid, hid, hid] +lib.H5Rcreate_object.argtypes = [hid, ctypes.c_char_p, hid, ctypes.c_void_p] +lib.H5Dwrite.argtypes = [hid, hid, hid, hid, hid, ctypes.c_void_p] +lib.H5Dclose.argtypes = [hid] + +with h5py.File(sys.argv[1], "w", libver="latest") as f: + f.create_dataset("target", data=np.arange(5, dtype="uwDa20ohR0SouGLA&$=vRjEEjs+`?#P4hBYdU}S(&5OUOrK%+(l=G;9n#?Yjj wglAtGnl2x8^=Jr;hQP=P0Z@euDxH~u7*uI8F>t`DZ6FT>fFiJJ8zR5}01qQlYXATM literal 0 HcmV?d00001 diff --git a/crates/clawhdf5-format/tests/reference_tests.rs b/crates/clawhdf5-format/tests/reference_tests.rs index af44b15..b1dcf40 100644 --- a/crates/clawhdf5-format/tests/reference_tests.rs +++ b/crates/clawhdf5-format/tests/reference_tests.rs @@ -316,3 +316,97 @@ print('ok') // Clean up let _ = std::fs::remove_file(&path); } + +// --------------------------------------------------------------------------- +// H5T_STD_REF (HDF5 1.12+ references, datatype message version 4) +// --------------------------------------------------------------------------- + +/// `fixtures/std_ref_hdf5_2_0.h5` (see `gen_std_ref.py`) holds a dataset of +/// `H5T_STD_REF` with two object references, written by HDF5 2.0 itself. The +/// datatype used to be rejected with `InvalidReferenceType(2)`. +#[test] +fn std_ref_object_references_from_hdf5_2_0() { + use clawhdf5_format::data_layout::DataLayout; + use clawhdf5_format::dataspace::Dataspace; + use clawhdf5_format::group_v2::resolve_path_any; + use clawhdf5_format::message_type::MessageType; + use clawhdf5_format::object_header::ObjectHeader; + use clawhdf5_format::signature::find_signature; + use clawhdf5_format::superblock::Superblock; + + let bytes: &[u8] = include_bytes!("fixtures/std_ref_hdf5_2_0.h5"); + let sb = Superblock::parse(bytes, find_signature(bytes).unwrap()).unwrap(); + let (os, ls) = (sb.offset_size, sb.length_size); + + let refs_addr = resolve_path_any(bytes, &sb, "refs").unwrap(); + let header = ObjectHeader::parse(bytes, refs_addr as usize, os, ls).unwrap(); + let message = |t: MessageType| { + &header + .messages + .iter() + .find(|m| m.msg_type == t) + .unwrap() + .data + }; + + let (datatype, _) = Datatype::parse(message(MessageType::Datatype)).unwrap(); + assert_eq!( + datatype, + Datatype::Reference { + size: 18, + ref_type: ReferenceType::Object2 + } + ); + let dataspace = Dataspace::parse(message(MessageType::Dataspace), ls).unwrap(); + let layout = DataLayout::parse(message(MessageType::DataLayout), os, ls).unwrap(); + let raw = + clawhdf5_format::data_read::read_raw_data(bytes, &layout, &dataspace, &datatype).unwrap(); + assert_eq!(raw.len(), 2 * 18); + + // The references point at the objects they were created from. + let refs = read_object_references(&raw, &datatype, os).unwrap(); + let addresses: Vec = refs.iter().map(|r| r.address).collect(); + assert_eq!( + addresses, + [ + resolve_path_any(bytes, &sb, "target").unwrap(), + resolve_path_any(bytes, &sb, "grp").unwrap(), + ] + ); + // And what they point at is a real object header. + for address in addresses { + ObjectHeader::parse(bytes, address as usize, os, ls).unwrap(); + } +} + +#[test] +fn std_ref_decoding_rejects_malformed_elements() { + let dt = Datatype::Reference { + size: 18, + ref_type: ReferenceType::Object2, + }; + let mut good = vec![0u8; 18]; + good[..4].copy_from_slice(&[2, 0, 8, 0xb3]); + assert_eq!( + read_object_references(&good, &dt, 8).unwrap()[0].address, + 0xb3 + ); + + // Null reference. + assert_eq!( + read_object_references(&[0u8; 18], &dt, 8).unwrap()[0].address, + u64::MAX + ); + for (what, patch) in [ + ("wrong reference type", (0usize, 3u8)), + ("external flag", (1, 1)), + ("token longer than the element", (2, 200)), + ("zero-length token", (2, 0)), + ] { + let mut bad = good.clone(); + bad[patch.0] = patch.1; + assert!(read_object_references(&bad, &dt, 8).is_err(), "{what}"); + } + // Not a whole number of elements. + assert!(read_object_references(&good[..17], &dt, 8).is_err()); +} diff --git a/docs/known-issues.md b/docs/known-issues.md index d3c8993..c762b1a 100644 --- a/docs/known-issues.md +++ b/docs/known-issues.md @@ -62,13 +62,21 @@ only files using the native type through the C API / h5py low-level API hit this ## Revised reference datatype (class 7, version 4) is not parsed -**Status:** open, unconfirmed against a real file. +**Status:** fixed 2026-09-19 for object references; region and attribute +references are recognised but not decoded. -**Summary:** HDF5 1.12+ `H5T_STD_REF` references use datatype version 4 with -reference types 2–4 (object2 / region2 / attribute), which `Datatype::parse` -rejects with `InvalidReferenceType`. h5py still writes the legacy v1 -object/region references, which read correctly, so no reproducing file has been -generated yet; one written with the C API (`H5T_STD_REF`) is needed. +**Summary:** HDF5 1.12+ `H5T_STD_REF` references use datatype message version 4 +with reference types 2-4 (object / region / attribute), which `Datatype::parse` +rejected with `InvalidReferenceType`. h5py still writes the legacy references, +so no file had been available to test against. + +**Fix:** a real file was produced by driving the libhdf5 bundled in the h5py +wheel through ctypes (`tests/fixtures/gen_std_ref.py` -> +`std_ref_hdf5_2_0.h5`). The three new types parse as +`ReferenceType::{Object2, DatasetRegion2, Attribute}`, and +`read_object_references` decodes `Object2` elements (type, flags, token size, +token = target object header address). External references (flag bit 0) and +the region/attribute payloads are errors rather than misreads. ## `clawhdf5-gpu` `gpu_tests` can hang under the default parallel test runner