diff --git a/BENCHMARKS.md b/BENCHMARKS.md index e8ba435..44cb0f7 100644 --- a/BENCHMARKS.md +++ b/BENCHMARKS.md @@ -28,6 +28,90 @@ --- +## 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 | + +### 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 0c75b95..705e7ac 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,45 @@ ## 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. +- **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 + 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 + 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 + `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/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 8ba9332..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,116 @@ 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; + +/// 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 +614,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 +768,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 +810,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); } @@ -777,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); + } + } } } @@ -985,29 +1086,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 +1128,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 +1240,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 +1282,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() { @@ -2284,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/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/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index e6c4271..0e6773f 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( @@ -858,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 @@ -885,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); @@ -975,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); @@ -1044,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); @@ -1126,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); @@ -1407,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), @@ -1414,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/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/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-format/src/filters.rs b/crates/clawhdf5-format/src/filters.rs index 2e99d70..9979d3f 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.as_chunks_mut::().0.iter_mut().enumerate() { + for (byte, plane) in element.iter_mut().zip(&planes) { + *byte = plane[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" + ); + } + } + } } 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-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=" = 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/crates/clawhdf5/src/reader.rs b/crates/clawhdf5/src/reader.rs index fb17557..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)?) } @@ -458,6 +483,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/h5py_interop_tests.rs b/crates/clawhdf5/tests/h5py_interop_tests.rs index 2402a07..e177b72 100644 --- a/crates/clawhdf5/tests/h5py_interop_tests.rs +++ b/crates/clawhdf5/tests/h5py_interop_tests.rs @@ -913,3 +913,128 @@ 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"] + ); +} + +// --------------------------------------------------------------------------- +// 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] + ); +} 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] + ); +} diff --git a/docs/known-issues.md b/docs/known-issues.md index 5eecaac..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 @@ -116,16 +124,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