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] + ); +}