diff --git a/CHANGELOG.md b/CHANGELOG.md index 9030c3c..ad55d54 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,44 @@ its lookups cost a few percent at 16 threads. Throughput with the default pool is unchanged, and still short of an h5py process pool. +### Contiguous read speed (2026-09-26) +- **Large read buffers are backed by transparent huge pages.** A full read + of a contiguous dataset was one `memcpy` from the mapped file, yet ran at + a quarter of h5py's speed on one thread: the fresh output `Vec` took a + page fault (and a kernel page clear) for every 4 KiB page it was written + to, 16384 of them for 64 MiB, and those cost several times the copy. + numpy, and so h5py, asks for transparent huge pages on every allocation of + 4 MiB or more; clawhdf5-format's read buffers now do too + (`madvise(MADV_HUGEPAGE)` on Linux, `libc` added as a Linux-only + dependency; a no-op elsewhere or when THP is disabled). It applies to the + typed readers' output (`read_f32`, `read_f64`, `read_i32`, `read_i64`, + `read_u64`, both byte orders), the raw contiguous read and the chunk + assembly buffer. Values are unchanged; new h5py comparison + `crates/clawhdf5/tests/contiguous_read_interop.rs` covers every 1-8-byte + integer and float type in both byte orders, ranks 1-4, and datasets past + the 4 MiB threshold. +- **Hyperslab and point reads of contiguous data copy runs, not elements.** + A 256 x 256 hyperslab of a contiguous `f32` dataset read at an eighth of + h5py's speed: the selection's bounding box was copied out of the file, + then walked element by element (a recursive call and two bounds checks per + element) into a second buffer, which `read_f32_selection` converted into + a third. Selections of contiguous data are now copied straight from the + file, one `memcpy` per run of elements that is contiguous in the file + (a block along the last dimension, blocks that touch, and whole rows when + the inner dimensions are selected in full, merged), with no zero-filled + intermediate; a selection covering most of the dataset no longer makes a + full copy first. The typed selection readers (`read_f32_selection`, + `read_f64_selection`, `read_i32_selection`, `read_i64_selection`) copy + directly into their output when the dataset stores that type natively, + and convert as before otherwise (big-endian, other widths). The chunked + paths use the same run-based extraction. New public + `clawhdf5_format::data_read::read_selection_native` and the sealed + `NativeElement` trait (also used by the `read_as_*` fast paths, which + gained one for native `u64`). Values are unchanged: checked against h5py + by `contiguous_read_interop.rs` (strided, blocked, adjacent-block and + whole-row hyperslabs, points, empty selections; every type, both byte + orders, ranks 1-4). + ### Plugin filters (2026-09-26) - **LZF, bitshuffle, bzip2 and Blosc read and write, in pure Rust.** Files written by h5py with `compression="lzf"`, or with hdf5plugin's diff --git a/conformance/probe/Cargo.lock b/conformance/probe/Cargo.lock index 964e0ab..36cfe85 100644 --- a/conformance/probe/Cargo.lock +++ b/conformance/probe/Cargo.lock @@ -64,6 +64,7 @@ dependencies = [ "bzip2", "flate2", "libaec-sys", + "libc", "lz4_flex", "pco", "portable-atomic", diff --git a/crates/clawhdf5-format/Cargo.toml b/crates/clawhdf5-format/Cargo.toml index 75fd233..1d1998c 100644 --- a/crates/clawhdf5-format/Cargo.toml +++ b/crates/clawhdf5-format/Cargo.toml @@ -30,6 +30,10 @@ ruzstd = { version = "0.9", optional = true } bzip2 = { version = "0.6", optional = true } snap = { version = "1", optional = true } +[target.'cfg(target_os = "linux")'.dependencies] +# madvise(MADV_HUGEPAGE) for large read buffers (see src/bulk_alloc.rs). +libc = { version = "0.2", default-features = false } + [dev-dependencies] half = { workspace = true } serde_json = "1" diff --git a/crates/clawhdf5-format/src/bulk_alloc.rs b/crates/clawhdf5-format/src/bulk_alloc.rs new file mode 100644 index 0000000..085f125 --- /dev/null +++ b/crates/clawhdf5-format/src/bulk_alloc.rs @@ -0,0 +1,79 @@ +//! Large output buffers backed by transparent huge pages where the OS offers +//! them. +//! +//! A fresh multi-megabyte `Vec` is mapped lazily by the kernel: the first +//! write to each 4 KiB page takes a page fault, and the kernel zeroes the page +//! before handing it over. For a 64 MiB read that is 16384 faults, and they +//! cost far more than the copy that fills the buffer — single-threaded +//! contiguous reads ran at about a quarter of h5py's speed because of them. +//! numpy (so h5py) avoids this by asking for transparent huge pages +//! (`madvise(MADV_HUGEPAGE)`) on every allocation of 4 MiB or more, which +//! turns 512 faults into one; this module does the same. +//! +//! The advice only changes how the pages are backed, never their contents, so +//! it is harmless when it cannot be honoured (THP disabled, not Linux, a +//! region that is part of the heap): the buffer is then exactly what it would +//! have been without it. + +#[cfg(not(feature = "std"))] +use alloc::vec::Vec; + +/// Buffers smaller than this are left alone (numpy uses the same threshold). +#[cfg(any(target_os = "linux", test))] +pub(crate) const HUGE_PAGE_THRESHOLD: usize = 4 << 20; + +/// Advise the kernel to back `[ptr, ptr + len)` with transparent huge pages, +/// when `len` is large enough to benefit. Call it before the first write so +/// the faults happen at huge-page granularity. +#[inline] +pub(crate) fn advise_huge_pages(ptr: *const u8, len: usize) { + #[cfg(target_os = "linux")] + if len >= HUGE_PAGE_THRESHOLD { + const PAGE: usize = 4096; + let start = (ptr as usize).next_multiple_of(PAGE); + let end = (ptr as usize + len) & !(PAGE - 1); + if end > start { + // SAFETY: `[start, end)` lies inside an allocation of `len` bytes + // at `ptr` that the caller owns, and is page aligned as madvise + // requires. MADV_HUGEPAGE does not change the memory's contents or + // validity; on failure (EINVAL when THP is compiled out, etc.) the + // region is simply left as it was, so the result is ignored. + unsafe { + libc::madvise(start as *mut libc::c_void, end - start, libc::MADV_HUGEPAGE); + } + } + } + #[cfg(not(target_os = "linux"))] + let _ = (ptr, len); +} + +/// `Vec::with_capacity(count)` for a buffer about to be filled in bulk, with +/// huge-page advice when it is large (see the module docs). +#[inline] +pub(crate) fn vec_for_bulk(count: usize) -> Vec { + let v: Vec = Vec::with_capacity(count); + advise_huge_pages( + v.as_ptr().cast::(), + v.capacity().saturating_mul(core::mem::size_of::()), + ); + v +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn bulk_vec_is_an_ordinary_vec() { + for count in [0usize, 1, 1000, HUGE_PAGE_THRESHOLD / 4 + 3] { + let mut v: Vec = vec_for_bulk(count); + assert!(v.capacity() >= count); + v.extend((0..count as u32).map(|i| i.wrapping_mul(2654435761))); + assert!( + v.iter() + .enumerate() + .all(|(i, &x)| x == (i as u32).wrapping_mul(2654435761)) + ); + } + } +} diff --git a/crates/clawhdf5-format/src/chunked_read.rs b/crates/clawhdf5-format/src/chunked_read.rs index 9169ab7..ffc1379 100644 --- a/crates/clawhdf5-format/src/chunked_read.rs +++ b/crates/clawhdf5-format/src/chunked_read.rs @@ -278,6 +278,8 @@ pub(crate) fn alloc_output(len: usize) -> Result, FormatError> { if ptr.is_null() { return Err(failed()); } + // Before anything writes to it, so a large buffer faults in huge pages. + crate::bulk_alloc::advise_huge_pages(ptr, len); // 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). diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index b8f2ffe..633ecd7 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -180,7 +180,9 @@ fn read_raw_data_full_impl( }); } ensure_len(file_data, addr, sz)?; - Ok(file_data[addr..addr + sz].to_vec()) + let mut out = crate::bulk_alloc::vec_for_bulk(sz); + out.extend_from_slice(&file_data[addr..addr + sz]); + Ok(out) } DataLayout::Chunked { .. } => read_chunked_data( file_data, @@ -529,6 +531,11 @@ pub fn extract_selection_from_buffer( block, } => { let rank = dims.len(); + if [start.len(), stride.len(), count.len(), block.len()] != [rank; 4] { + return Err(FormatError::SelectionOutOfBounds(format!( + "hyperslab rank does not match dataset rank {rank}" + ))); + } let output_elements = count .iter() .zip(block.iter()) @@ -538,96 +545,40 @@ pub fn extract_selection_from_buffer( crate::chunked_read::checked_byte_len(output_elements, elem_size)?, )?; - // Compute dataset strides (row-major) - let mut ds_strides = vec![1usize; rank]; - for i in (0..rank.saturating_sub(1)).rev() { - ds_strides[i] = ds_strides[i + 1] * dims[i + 1] as usize; - } - - // Compute output shape and strides - let output_dims: Vec = count - .iter() - .zip(block.iter()) - .map(|(&c, &b)| (c * b) as usize) - .collect(); - let mut out_strides = vec![1usize; rank]; - for i in (0..rank.saturating_sub(1)).rev() { - out_strides[i] = out_strides[i + 1] * output_dims[i + 1]; - } - - // Iterate over all selected elements - // For each block in the hyperslab, copy the elements - let mut out_linear = 0usize; - let _block_coords = vec![0u64; rank]; - - #[allow(clippy::too_many_arguments)] - fn iterate_hyperslab( - d: usize, - rank: usize, - start: &[u64], - stride: &[u64], - count: &[u64], - block: &[u64], - dims: &[u64], - ds_strides: &[usize], - elem_size: usize, - full_data: &[u8], - output: &mut [u8], - out_linear: &mut usize, - current_ds_offset: usize, - ) { - if d == rank { - // Copy one element - let src = current_ds_offset * elem_size; - let dst = *out_linear * elem_size; - if src + elem_size <= full_data.len() && dst + elem_size <= output.len() { - output[dst..dst + elem_size] - .copy_from_slice(&full_data[src..src + elem_size]); - } - *out_linear += 1; - return; - } - - for bi in 0..count[d] { - let block_start = start[d] + bi * stride[d]; - for bj in 0..block[d] { - let coord = block_start + bj; - if coord < dims[d] { - iterate_hyperslab( - d + 1, - rank, - start, - stride, - count, - block, - dims, - ds_strides, - elem_size, - full_data, - output, - out_linear, - current_ds_offset + coord as usize * ds_strides[d], - ); + // One copy per run of elements contiguous in `full_data` + // (`gather`'s runs). Coordinates past the extent are skipped and + // runs past the end of `full_data` left as zeros, element by + // element, as this extractor always did; validated selections + // never hit either. + let mut out_at = 0usize; + crate::gather::hyperslab_runs(dims, start, stride, count, block, |first, n| { + let big = |v: u64| usize::try_from(v).unwrap_or(usize::MAX); + let (first, n) = (big(first), big(n)); + let len = n.saturating_mul(elem_size); + let src = first.saturating_mul(elem_size); + let out_end = out_at.saturating_add(len); + if let (Some(from), Some(to)) = ( + full_data.get(src..src.saturating_add(len)), + output.get_mut(out_at..out_end), + ) { + to.copy_from_slice(from); + } else { + for k in 0..n { + let s = first.saturating_add(k).saturating_mul(elem_size); + let o = out_at.saturating_add(k.saturating_mul(elem_size)); + if o >= output.len() { + break; + } + if let (Some(from), Some(to)) = ( + full_data.get(s..s.saturating_add(elem_size)), + output.get_mut(o..o.saturating_add(elem_size)), + ) { + to.copy_from_slice(from); } } } - } - - iterate_hyperslab( - 0, - rank, - start, - stride, - count, - block, - dims, - &ds_strides, - elem_size, - full_data, - &mut output, - &mut out_linear, - 0, - ); + out_at = out_end; + }); Ok(output) } @@ -755,22 +706,76 @@ 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. +mod sealed { + pub trait Sealed {} +} + +/// A numeric type whose values can be copied straight out of a dataset's +/// bytes when the dataset stores exactly that type in the target's byte +/// order: `u8`, `i32`, `i64`, `u64`, `f32` and `f64`. +/// +/// # Safety +/// +/// Implementors have no padding and no invalid bit patterns, so a buffer of +/// them may be filled by copying bytes. The trait is sealed. +pub unsafe trait NativeElement: sealed::Sealed + Copy + 'static { + /// Whether `datatype`'s stored bytes are this type's native in-memory + /// representation (same size, byte order, signedness, full precision, + /// IEEE layout), so reading needs a copy and no conversion. + fn is_native(datatype: &Datatype) -> bool; +} + +/// A full-width fixed-point type of `size` bytes and the given signedness in +/// the target's byte order. +fn is_native_int(datatype: &Datatype, size: u32, want_signed: bool) -> bool { + let order = if cfg!(target_endian = "little") { + DatatypeByteOrder::LittleEndian + } else { + DatatypeByteOrder::BigEndian + }; + matches!( + datatype, + Datatype::FixedPoint { size: s, signed, byte_order, .. } + if *s == size && *signed == want_signed && (size == 1 || *byte_order == order) + ) && is_full_width(datatype) +} + +macro_rules! native_element { + ($($t:ty => |$dt:ident| $check:expr;)*) => {$( + impl sealed::Sealed for $t {} + // SAFETY: a primitive integer or float: no padding, and every bit + // pattern is a valid value. + unsafe impl NativeElement for $t { + fn is_native($dt: &Datatype) -> bool { + $check + } + } + )*}; +} + +native_element! { + u8 => |dt| is_native_int(dt, 1, false); + i32 => |dt| is_native_int(dt, 4, true); + i64 => |dt| is_native_int(dt, 8, true); + u64 => |dt| is_native_int(dt, 8, false); + f32 => |dt| cfg!(target_endian = "little") && is_native_le_float(dt, FloatFormat::Single); + f64 => |dt| cfg!(target_endian = "little") && is_native_le_float(dt, FloatFormat::Double); +} + +/// Copy `count` values of `T` out of `raw`, which holds them in `T`'s native +/// representation (see [`NativeElement::is_native`]), 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 { +fn native_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); + assert!(bytes <= raw.len(), "native_to_vec: source too short"); + let mut result: Vec = crate::bulk_alloc::vec_for_bulk(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 + // bytes; `raw` holds at least `bytes` bytes (asserted); the regions + // cannot overlap because `result` was just allocated. `T: NativeElement` + // 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); @@ -779,6 +784,44 @@ fn native_le_to_vec(raw: &[u8], count: usize) -> Vec { result } +/// Read `selection` of a dataset whose raw bytes (all of them, row-major, of +/// shape `dims`) are `raw` — typically a contiguous dataset's bytes borrowed +/// from the file — straight into a `Vec`, copying each contiguous run of +/// selected elements once. +/// +/// Returns `Ok(None)` when `datatype` is not `T`'s native representation +/// ([`NativeElement::is_native`]); the caller then converts through +/// [`read_raw_data_selection`] and the `read_as_*` functions. The selection is +/// validated like every selection read: out-of-range coordinates are +/// [`FormatError::SelectionOutOfBounds`]. +pub fn read_selection_native( + raw: &[u8], + dims: &[u64], + datatype: &Datatype, + selection: &crate::selection::Selection, +) -> Result>, FormatError> { + if !T::is_native(datatype) { + return Ok(None); + } + let elem_size = core::mem::size_of::(); + let total = dims + .iter() + .try_fold(1u64, |acc, &d| acc.checked_mul(d)) + .ok_or_else(|| FormatError::Overflow("dataset shape overflows".into()))?; + let expected = crate::chunked_read::checked_byte_len(total, elem_size)?; + if raw.len() != expected { + return Err(FormatError::DataSizeMismatch { + expected, + actual: raw.len(), + }); + } + if let crate::selection::Selection::All = selection { + return Ok(Some(native_to_vec(raw, expected / elem_size))); + } + crate::partial_read::validate(selection, dims)?; + crate::gather::gather::(raw, dims, elem_size, selection).map(Some) +} + /// Convert raw bytes to `f64` values. pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result, FormatError> { // Array datatypes read as a flat sequence of their base elements, and @@ -797,13 +840,12 @@ pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result, FormatEr let count = raw.len() / elem_size; // Fast path: native-endian f64 — single bulk memcpy - #[cfg(target_endian = "little")] - if is_native_le_float(datatype, FloatFormat::Double) { - return Ok(native_le_to_vec::(raw, count)); + if f64::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); } let order = get_byte_order(datatype); - let mut result = Vec::with_capacity(count); + let mut result = crate::bulk_alloc::vec_for_bulk(count); if let Datatype::FloatingPoint { .. } = datatype { let format = FloatFormat::of(datatype)?; for chunk in raw.chunks_exact(elem_size) { @@ -939,23 +981,12 @@ pub fn read_as_i64(raw: &[u8], datatype: &Datatype) -> Result, FormatEr let count = raw.len() / elem_size; // Fast path: native LE i64 — single bulk memcpy - #[cfg(target_endian = "little")] - if elem_size == 8 - && is_full_width(datatype) - && matches!( - datatype, - Datatype::FixedPoint { - byte_order: DatatypeByteOrder::LittleEndian, - signed: true, - .. - } - ) - { - return Ok(native_le_to_vec::(raw, count)); + if i64::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); } let order = get_byte_order(datatype); - let mut result = Vec::with_capacity(count); + let mut result = crate::bulk_alloc::vec_for_bulk(count); for i in 0..count { let chunk = &raw[i * elem_size..(i + 1) * elem_size]; result.push(decode_scalar(chunk, datatype, &order)?.to_i64()); @@ -984,8 +1015,14 @@ pub fn read_as_u64(raw: &[u8], datatype: &Datatype) -> Result, FormatEr }); } let count = raw.len() / elem_size; + + // Fast path: native u64 — single bulk memcpy + if u64::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); + } + let order = get_byte_order(datatype); - let mut result = Vec::with_capacity(count); + let mut result = crate::bulk_alloc::vec_for_bulk(count); for i in 0..count { let chunk = &raw[i * elem_size..(i + 1) * elem_size]; result.push(decode_scalar(chunk, datatype, &order)?.to_u64()); @@ -1011,21 +1048,23 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result, FormatEr let count = raw.len() / elem_size; // Fast path: native-endian f32 — single bulk memcpy - #[cfg(target_endian = "little")] - if is_native_le_float(datatype, FloatFormat::Single) { - return Ok(native_le_to_vec::(raw, count)); + if f32::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); } // Little-endian IEEE half precision (numpy float16): widen directly. if is_native_le_float(datatype, FloatFormat::Half) { let (halves, _) = raw[..count * 2].as_chunks::<2>(); - return Ok(halves - .iter() - .map(|&b| f16_bits_to_f32(u16::from_le_bytes(b))) - .collect()); + let mut result = crate::bulk_alloc::vec_for_bulk(count); + result.extend( + halves + .iter() + .map(|&b| f16_bits_to_f32(u16::from_le_bytes(b))), + ); + return Ok(result); } let order = get_byte_order(datatype); - let mut result = Vec::with_capacity(count); + let mut result = crate::bulk_alloc::vec_for_bulk(count); if let Datatype::FloatingPoint { .. } = datatype { let format = FloatFormat::of(datatype)?; for chunk in raw.chunks_exact(elem_size) { @@ -1098,23 +1137,12 @@ pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result, FormatEr let count = raw.len() / elem_size; // Fast path: native LE i32 — single bulk memcpy - #[cfg(target_endian = "little")] - if elem_size == 4 - && is_full_width(datatype) - && matches!( - datatype, - Datatype::FixedPoint { - byte_order: DatatypeByteOrder::LittleEndian, - signed: true, - .. - } - ) - { - return Ok(native_le_to_vec::(raw, count)); + if i32::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); } let order = get_byte_order(datatype); - let mut result = Vec::with_capacity(count); + let mut result = crate::bulk_alloc::vec_for_bulk(count); for i in 0..count { let chunk = &raw[i * elem_size..(i + 1) * elem_size]; result.push(decode_scalar(chunk, datatype, &order)?.to_i32()); diff --git a/crates/clawhdf5-format/src/gather.rs b/crates/clawhdf5-format/src/gather.rs new file mode 100644 index 0000000..22d2ae7 --- /dev/null +++ b/crates/clawhdf5-format/src/gather.rs @@ -0,0 +1,343 @@ +//! Copying a selection out of a row-major buffer one contiguous run at a time. +//! +//! A selection's elements, in output order, fall into runs that are adjacent +//! in the source: a whole block along the last dimension, blocks that touch +//! (`stride == block`), and whole rows when the inner dimensions are selected +//! in full. Copying run by run turns a 256 x 256 hyperslab of a 1024-wide +//! dataset into 256 `memcpy`s of 1 KiB, where the old extractor recursed and +//! bounds-checked once per element. + +#[cfg(not(feature = "std"))] +use alloc::{vec, vec::Vec}; + +use crate::data_read::NativeElement; +use crate::error::FormatError; +use crate::selection::Selection; + +/// Row-major element strides of `dims` (the last dimension has stride 1). +fn strides(dims: &[u64]) -> Vec { + let mut s = vec![1u64; dims.len()]; + for d in (0..dims.len().saturating_sub(1)).rev() { + s[d] = s[d + 1].wrapping_mul(dims[d + 1]); + } + s +} + +/// Merges adjacent runs before handing them on. +struct Coalesce { + start: u64, + len: u64, + emit: F, +} + +impl Coalesce { + #[inline] + fn push(&mut self, start: u64, len: u64) { + if len == 0 { + return; + } + if self.len > 0 && self.start.wrapping_add(self.len) == start { + self.len += len; + return; + } + self.flush(); + self.start = start; + self.len = len; + } + + fn flush(&mut self) { + if self.len > 0 { + (self.emit)(self.start, self.len); + self.len = 0; + } + } +} + +/// Call `emit(first_element, element_count)` for each run of a hyperslab's +/// elements that is contiguous in a row-major dataset of shape `dims`, in +/// the order the selection returns them. Adjacent runs are merged. +/// +/// Coordinates at or past a dimension's extent are skipped, as the +/// element-wise extractor always did; callers that want them to be an error +/// validate the selection first. The four vectors must have `dims.len()` +/// entries. +pub(crate) fn hyperslab_runs( + dims: &[u64], + start: &[u64], + stride: &[u64], + count: &[u64], + block: &[u64], + emit: impl FnMut(u64, u64), +) { + let rank = dims.len(); + let mut out = Coalesce { + start: 0, + len: 0, + emit, + }; + if rank == 0 { + out.push(0, 1); + out.flush(); + return; + } + if (0..rank).any(|d| count[d] == 0 || block[d] == 0) { + return; + } + let strides = strides(dims); + let last = rank - 1; + // Odometer over the outer dimensions: (block index, offset in block). + let mut ci = vec![0u64; last]; + let mut bi = vec![0u64; last]; + 'outer: loop { + // Base offset of this row, or skip it if a coordinate is out of range. + let mut base = 0u64; + let mut in_range = true; + for d in 0..last { + let coord = start[d] + .saturating_add(ci[d].saturating_mul(stride[d])) + .saturating_add(bi[d]); + if coord >= dims[d] { + in_range = false; + break; + } + base = base.wrapping_add(coord.wrapping_mul(strides[d])); + } + if in_range && (stride[last] == block[last] || count[last] == 1) { + // Blocks that touch (the common unit-stride case: block 1, + // stride 1) are one range; don't split it into per-element runs. + let s = start[last]; + let e = s + .saturating_add(count[last].saturating_mul(block[last])) + .min(dims[last]); + if s < e { + out.push(base.wrapping_add(s), e - s); + } + } else if in_range { + for c in 0..count[last] { + let s = start[last].saturating_add(c.saturating_mul(stride[last])); + if s >= dims[last] { + continue; + } + let e = s.saturating_add(block[last]).min(dims[last]); + out.push(base.wrapping_add(s), e - s); + } + } + // Advance the odometer, last outer dimension fastest. + let mut d = last; + loop { + if d == 0 { + break 'outer; + } + d -= 1; + bi[d] += 1; + if bi[d] < block[d] { + break; + } + bi[d] = 0; + ci[d] += 1; + if ci[d] < count[d] { + break; + } + ci[d] = 0; + } + } + out.flush(); +} + +/// The selected elements of `src` — a row-major dataset of shape `dims` and +/// `elem_size`-byte elements — copied into a fresh `Vec`, one `memcpy` per +/// contiguous run, with no zero-filling of the output first. +/// +/// For `T` other than `u8`, `elem_size` must equal `size_of::()`. The +/// selection must be a validated hyperslab, point list or `None` (`All` is the +/// caller's to handle); `src` must hold exactly the dataset. Anything that +/// would read outside `src` is an error, never a partial result. +pub(crate) fn gather( + src: &[u8], + dims: &[u64], + elem_size: usize, + selection: &Selection, +) -> Result, FormatError> { + let t_size = core::mem::size_of::(); + if elem_size == 0 || (t_size != 1 && t_size != elem_size) { + return Err(FormatError::DataSizeMismatch { + expected: t_size, + actual: elem_size, + }); + } + let n_elements = match selection { + Selection::None => 0, + Selection::Hyperslab { count, block, .. } => count + .iter() + .zip(block) + .try_fold(1u64, |acc, (&c, &b)| acc.checked_mul(c.checked_mul(b)?)) + .ok_or_else(|| FormatError::Overflow("hyperslab count x block overflows".into()))?, + Selection::Points(points) => points.len() as u64, + Selection::All => { + return Err(FormatError::SelectionOutOfBounds( + "gather does not take Selection::All".into(), + )); + } + }; + let out_bytes = crate::chunked_read::checked_byte_len(n_elements, elem_size)?; + let out_len = out_bytes / t_size; + let mut out: Vec = crate::bulk_alloc::vec_for_bulk(out_len); + let dst = out.as_mut_ptr().cast::(); + let mut written = 0usize; + let mut failed = false; + let mut copy_run = |first: u64, n: u64| { + if failed { + return; + } + let range = usize::try_from(first) + .ok() + .and_then(|f| f.checked_mul(elem_size)) + .zip( + usize::try_from(n) + .ok() + .and_then(|n| n.checked_mul(elem_size)), + ) + .and_then(|(at, len)| Some((at, len, at.checked_add(len)?))); + match range { + Some((at, len, end)) if end <= src.len() && written + len <= out_bytes => { + // SAFETY: `src[at..end]` is in bounds (checked above), and + // `dst + written .. + len` lies within `out`'s capacity of + // `out_bytes` bytes (checked above); `out` is a fresh + // allocation, so the regions do not overlap. + unsafe { + core::ptr::copy_nonoverlapping(src.as_ptr().add(at), dst.add(written), len) + }; + written += len; + } + _ => failed = true, + } + }; + let mut bad_point = false; + match selection { + Selection::Hyperslab { + start, + stride, + count, + block, + } => { + let rank = dims.len(); + if [start.len(), stride.len(), count.len(), block.len()] != [rank; 4] { + return Err(FormatError::SelectionOutOfBounds( + "hyperslab rank does not match dataset rank".into(), + )); + } + hyperslab_runs(dims, start, stride, count, block, &mut copy_run); + } + Selection::Points(points) => { + let strides = strides(dims); + let mut runs = Coalesce { + start: 0, + len: 0, + emit: &mut copy_run, + }; + for p in points { + if p.len() != dims.len() || p.iter().zip(dims).any(|(c, n)| c >= n) { + bad_point = true; + break; + } + let at = p + .iter() + .zip(&strides) + .fold(0u64, |acc, (c, s)| acc.wrapping_add(c.wrapping_mul(*s))); + runs.push(at, 1); + } + runs.flush(); + } + Selection::None | Selection::All => {} + } + if failed || bad_point || written != out_bytes { + return Err(FormatError::SelectionOutOfBounds( + "selection addresses elements outside the dataset".into(), + )); + } + // SAFETY: all `out_bytes` bytes, i.e. `out_len` values of `T`, were + // written above, and every bit pattern is a valid `T` (`NativeElement`). + unsafe { out.set_len(out_len) }; + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn runs(dims: &[u64], sel: [&[u64]; 4]) -> Vec<(u64, u64)> { + let mut v = Vec::new(); + hyperslab_runs(dims, sel[0], sel[1], sel[2], sel[3], |s, n| v.push((s, n))); + v + } + + #[test] + fn runs_merge_blocks_and_whole_rows() { + // A box: one run per row. + assert_eq!( + runs(&[4, 10], [&[1, 2], &[1, 1], &[2, 3], &[1, 1]]), + vec![(12, 3), (22, 3)] + ); + // Whole rows: one run. + assert_eq!( + runs(&[4, 10], [&[1, 0], &[1, 1], &[3, 10], &[1, 1]]), + vec![(10, 30)] + ); + // stride == block: blocks merge. + assert_eq!( + runs(&[1, 10], [&[0, 1], &[1, 2], &[1, 4], &[1, 2]]), + vec![(1, 8)] + ); + // Strided with blocks along both dimensions. + assert_eq!( + runs(&[6, 10], [&[0, 1], &[3, 4], &[2, 2], &[2, 2]]), + vec![ + (1, 2), + (5, 2), + (11, 2), + (15, 2), + (31, 2), + (35, 2), + (41, 2), + (45, 2) + ] + ); + // Empty. + assert!(runs(&[4, 10], [&[0, 0], &[1, 1], &[0, 3], &[1, 1]]).is_empty()); + // Scalar. + assert_eq!(runs(&[], [&[], &[], &[], &[]]), vec![(0, 1)]); + } + + #[test] + fn gather_matches_element_order_and_rejects_out_of_range() { + let dims = [3u64, 4]; + let src: Vec = (0..12u16).flat_map(|v| v.to_le_bytes()).collect(); + let sel = Selection::Hyperslab { + start: vec![0, 1], + stride: vec![2, 2], + count: vec![2, 2], + block: vec![1, 1], + }; + let got: Vec = gather(&src, &dims, 2, &sel).unwrap(); + let want: Vec = [1u16, 3, 9, 11] + .iter() + .flat_map(|v| v.to_le_bytes()) + .collect(); + assert_eq!(got, want); + let pts = Selection::Points(vec![vec![2, 3], vec![0, 0], vec![0, 1]]); + let got: Vec = gather(&src, &dims, 2, &pts).unwrap(); + let want: Vec = [11u16, 0, 1].iter().flat_map(|v| v.to_le_bytes()).collect(); + assert_eq!(got, want); + // Past the extent, or a source shorter than the dataset: an error. + let bad = Selection::Points(vec![vec![3, 0]]); + assert!(gather::(&src, &dims, 2, &bad).is_err()); + let past = Selection::Hyperslab { + start: vec![2, 0], + stride: vec![1, 1], + count: vec![2, 4], + block: vec![1, 1], + }; + assert!(gather::(&src, &dims, 2, &past).is_err()); + assert!(gather::(&src[..20], &dims, 2, &pts).is_err()); + } +} diff --git a/crates/clawhdf5-format/src/lib.rs b/crates/clawhdf5-format/src/lib.rs index 905ce84..23840b8 100644 --- a/crates/clawhdf5-format/src/lib.rs +++ b/crates/clawhdf5-format/src/lib.rs @@ -61,6 +61,7 @@ pub mod attribute; pub mod attribute_info; pub mod btree_v1; pub mod btree_v2; +mod bulk_alloc; pub mod checksum; pub mod chunk_cache; mod chunk_grid; @@ -93,6 +94,7 @@ mod filters_szip; pub mod fixed_array; pub mod float16; pub mod fractal_heap; +mod gather; pub mod global_heap; pub mod group_info; pub mod group_v1; diff --git a/crates/clawhdf5-format/src/partial_read.rs b/crates/clawhdf5-format/src/partial_read.rs index 7d73599..9865c46 100644 --- a/crates/clawhdf5-format/src/partial_read.rs +++ b/crates/clawhdf5-format/src/partial_read.rs @@ -3,11 +3,13 @@ //! //! [`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. +//! large dataset took about as long as reading all of it. A contiguous +//! dataset's selection is now copied straight out of the file, one `memcpy` +//! per contiguous run of selected elements (`crate::gather`). For chunked +//! data the selection's bounding box is materialised — only the chunks that +//! overlap it — and the 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; @@ -250,10 +252,33 @@ pub fn read_selection( if dims.is_empty() || elem_size == 0 { return Ok(None); } + let total = dataspace.checked_num_elements()?; + // Contiguous data is addressable in place: copy the selection's runs + // straight out of it, whatever fraction of the dataset it covers, with no + // intermediate box (and no full copy for a large selection). + if let ( + DataLayout::Contiguous { + address: Some(address), + .. + }, + Selection::Hyperslab { .. } | Selection::Points(_), + ) = (layout, selection) + { + validate(selection, dims)?; + 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(), + })?; + return crate::gather::gather::(data, dims, elem_size, selection).map(Some); + } 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)) @@ -265,30 +290,6 @@ pub fn read_selection( 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(_), .. diff --git a/crates/clawhdf5/src/reader.rs b/crates/clawhdf5/src/reader.rs index 83cbfa8..74c268d 100644 --- a/crates/clawhdf5/src/reader.rs +++ b/crates/clawhdf5/src/reader.rs @@ -569,9 +569,7 @@ impl<'f> Dataset<'f> { &self, selection: &clawhdf5_format::selection::Selection, ) -> Result, Error> { - let raw = self.read_selection(selection)?; - let dt = self.datatype()?; - Ok(data_read::read_as_f64(&raw, &dt)?) + self.read_typed_selection(selection, data_read::read_as_f64, || self.read_f64()) } /// Read selected elements as `f32` values. @@ -579,9 +577,7 @@ impl<'f> Dataset<'f> { &self, selection: &clawhdf5_format::selection::Selection, ) -> Result, Error> { - let raw = self.read_selection(selection)?; - let dt = self.datatype()?; - Ok(data_read::read_as_f32(&raw, &dt)?) + self.read_typed_selection(selection, data_read::read_as_f32, || self.read_f32()) } /// Read selected elements as `i32` values. @@ -589,9 +585,7 @@ impl<'f> Dataset<'f> { &self, selection: &clawhdf5_format::selection::Selection, ) -> Result, Error> { - let raw = self.read_selection(selection)?; - let dt = self.datatype()?; - Ok(data_read::read_as_i32(&raw, &dt)?) + self.read_typed_selection(selection, data_read::read_as_i32, || self.read_i32()) } /// Read selected elements as `i64` values. @@ -599,9 +593,34 @@ impl<'f> Dataset<'f> { &self, selection: &clawhdf5_format::selection::Selection, ) -> Result, Error> { - let raw = self.read_selection(selection)?; + self.read_typed_selection(selection, data_read::read_as_i64, || self.read_i64()) + } + + /// The typed selection readers. `All` is a full read. A contiguous dataset + /// that stores `T` natively is copied from the file straight into the + /// `Vec`, one copy per contiguous run of selected elements; anything + /// else reads the selection's bytes and converts them with `convert`. + fn read_typed_selection( + &self, + selection: &clawhdf5_format::selection::Selection, + convert: fn(&[u8], &Datatype) -> Result, FormatError>, + full: impl FnOnce() -> Result, Error>, + ) -> Result, Error> { + if matches!(selection, clawhdf5_format::selection::Selection::All) { + return full(); + } let dt = self.datatype()?; - Ok(data_read::read_as_i64(&raw, &dt)?) + if T::is_native(&dt) + && let Ok(Some(raw)) = self.read_raw_ref() + { + let dims = self.dataspace()?.dimensions; + if let Some(values) = data_read::read_selection_native::(raw, &dims, &dt, selection)? + { + return Ok(values); + } + } + let raw = self.read_selection(selection)?; + Ok(convert(&raw, &dt)?) } /// Zero-copy read of contiguous raw data. diff --git a/crates/clawhdf5/tests/contiguous_read_interop.rs b/crates/clawhdf5/tests/contiguous_read_interop.rs new file mode 100644 index 0000000..188856e --- /dev/null +++ b/crates/clawhdf5/tests/contiguous_read_interop.rs @@ -0,0 +1,403 @@ +//! Reads of contiguous datasets — full reads and hyperslab/point selections, +//! through every typed reader — checked against h5py/libhdf5 for every +//! integer and float width, both byte orders, ranks 1 to 4, and datasets +//! larger than the huge-page threshold (4 MiB) the read buffers use. +//! +//! h5py writes the file and, for each selection, reads it with libhdf5's own +//! hyperslab/point selection (`select_hyperslab` with stride and block, +//! `select_elements`) and saves the raw bytes it gets back; the byte-level +//! [`Dataset::read_selection`] must return exactly those bytes, and the typed +//! readers the same values. Skipped when python3 with h5py is unavailable, +//! unless `CLAWHDF5_REQUIRE_INTEROP=1`. + +use std::path::Path; +use std::process::Command; + +use clawhdf5::File; +use clawhdf5_format::selection::Selection; + +fn python() -> String { + std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string()) +} + +fn interop_required() -> bool { + std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1") +} + +fn python_available() -> bool { + Command::new(python()) + .args(["-c", "import h5py, numpy"]) + .output() + .map(|o| o.status.success()) + .unwrap_or(false) +} + +macro_rules! skip_if_no_python { + () => { + if !python_available() { + assert!( + !interop_required(), + "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with h5py is not available" + ); + eprintln!("SKIP: python3 with h5py not available"); + return; + } + }; +} + +fn run_python(script: &str) { + let output = Command::new(python()) + .args(["-c", script]) + .output() + .expect("failed to run python"); + assert!( + output.status.success(), + "python failed:\n{}", + String::from_utf8_lossy(&output.stderr) + ); +} + +/// numpy type codes, with the modulus of the value pattern: every value is an +/// integer exactly representable in the type and in every typed reader's +/// output (f16 is exact below 2048, f32 below 2^24). +const DTYPES: [(&str, i64, bool); 11] = [ + ("i1", 201, true), + ("u1", 251, false), + ("i2", 2039, true), + ("u2", 2039, false), + ("i4", 1_000_003, true), + ("u4", 1_000_003, false), + ("i8", 1_000_003, true), + ("u8", 1_000_003, false), + ("f2", 2039, true), + ("f4", 1_000_003, true), + ("f8", 1_000_003, true), +]; + +/// Value of element `i` of a dataset of type `code` (the same formula as the +/// Python side): a permutation-revealing pattern, centred on 0 when signed. +fn value(code: &str, i: u64) -> i64 { + let (_, m, signed) = DTYPES.iter().find(|d| d.0 == code).unwrap(); + let v = ((i as i128 * 7919) % *m as i128) as i64; + if *signed { v - m / 2 } else { v } +} + +const SHAPES: [&[u64]; 4] = [&[1000], &[37, 53], &[7, 11, 13], &[3, 5, 7, 9]]; + +/// Datasets past the 4 MiB huge-page threshold, as (type, shape). +const BIG: [(&str, [u64; 2]); 4] = [ + ("f4", [1100, 1024]), + ("i4", [1100, 1024]), + ("f8", [600, 1024]), + ("i8", [600, 1024]), +]; + +fn datasets() -> Vec<(String, String, Vec)> { + let mut out = Vec::new(); + for (code, _, _) in DTYPES { + for (tag, _) in [("le", '<'), ("be", '>')] { + for shape in SHAPES { + out.push(( + format!("{code}{tag}_r{}", shape.len()), + code.to_string(), + shape.to_vec(), + )); + } + } + } + for (code, shape) in BIG { + for tag in ["le", "be"] { + out.push((format!("{code}{tag}_big"), code.to_string(), shape.to_vec())); + } + } + out +} + +fn write_file(path: &Path) { + let script = format!( + r#" +import h5py, numpy as np +M = {{'i1': 201, 'u1': 251, 'i2': 2039, 'u2': 2039, 'i4': 1000003, 'u4': 1000003, + 'i8': 1000003, 'u8': 1000003, 'f2': 2039, 'f4': 1000003, 'f8': 1000003}} +def values(code, n): + v = (np.arange(n, dtype=np.int64) * 7919) % M[code] + if code[0] != 'u': + v -= M[code] // 2 + return v +shapes = [(1000,), (37, 53), (7, 11, 13), (3, 5, 7, 9)] +big = [('f4', (1100, 1024)), ('i4', (1100, 1024)), ('f8', (600, 1024)), ('i8', (600, 1024))] +with h5py.File("{path}", "w") as f: + for code in M: + for tag, e in (('le', '<'), ('be', '>')): + for shape in shapes: + n = int(np.prod(shape)) + f.create_dataset(f"{{code}}{{tag}}_r{{len(shape)}}", + data=values(code, n).astype(e + code).reshape(shape)) + for code, shape in big: + for tag, e in (('le', '<'), ('be', '>')): + n = int(np.prod(shape)) + f.create_dataset(f"{{code}}{{tag}}_big", + data=values(code, n).astype(e + code).reshape(shape)) +"#, + path = path.display() + ); + run_python(&script); +} + +#[test] +fn full_reads_match_h5py_for_every_type_order_and_size() { + skip_if_no_python!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("contig.h5"); + write_file(&path); + let file = File::open(&path).unwrap(); + for (name, code, shape) in datasets() { + let ds = file.dataset(&name).unwrap(); + assert!(ds.read_raw_ref().unwrap().is_some(), "{name} is contiguous"); + let n: u64 = shape.iter().product(); + let want: Vec = (0..n).map(|i| value(&code, i)).collect(); + assert_eq!( + ds.read_f64().unwrap(), + want.iter().map(|&v| v as f64).collect::>(), + "{name} read_f64" + ); + assert_eq!( + ds.read_f32().unwrap(), + want.iter().map(|&v| v as f32).collect::>(), + "{name} read_f32" + ); + assert_eq!(ds.read_i64().unwrap(), want, "{name} read_i64"); + assert_eq!( + ds.read_i32().unwrap(), + want.iter().map(|&v| v as i32).collect::>(), + "{name} read_i32" + ); + // libhdf5 saturates negative values to 0 when reading as unsigned. + assert_eq!( + ds.read_u64().unwrap(), + want.iter().map(|&v| v.max(0) as u64).collect::>(), + "{name} read_u64" + ); + } +} + +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) + } +} + +/// A hyperslab from per-dimension `(start, stride, count, block)`. +fn slab(dims: &[(u64, u64, u64, u64)]) -> Selection { + Selection::Hyperslab { + start: dims.iter().map(|d| d.0).collect(), + stride: dims.iter().map(|d| d.1).collect(), + count: dims.iter().map(|d| d.2).collect(), + block: dims.iter().map(|d| d.3).collect(), + } +} + +/// Selections of every shape the read paths distinguish, all valid for `dims`. +fn selections(rng: &mut Rng, dims: &[u64]) -> Vec { + let mut out = Vec::new(); + // Unit-stride box. + out.push(slab( + &dims + .iter() + .map(|&n| { + let c = 1 + rng.below(n); + (rng.below(n - c + 1), 1, c, 1) + }) + .collect::>(), + )); + // Strided (block 1), blocked (stride > block), and adjacent blocks + // (stride == block, which reads like a box): (block, stride - block). + for (block, gap) in [(1, 1), (2, 1), (2, 0)] { + out.push(slab( + &dims + .iter() + .map(|&n| { + let b = (block + rng.below(2)).min(n); + let st = b + gap + rng.below(2) * gap; + let s = rng.below(n - b + 1); + let c = 1 + rng.below((n - s - b) / st + 1); + (s, st, c, b) + }) + .collect::>(), + )); + } + // Whole inner rows (one run across rows), and the whole dataset. + let r0 = rng.below(dims[0]); + let mut rows = vec![(r0, 1, 1 + rng.below(dims[0] - r0), 1)]; + rows.extend(dims[1..].iter().map(|&n| (0, 1, n, 1))); + out.push(slab(&rows)); + out.push(slab( + &dims.iter().map(|&n| (0, 1, n, 1)).collect::>(), + )); + // One element. + out.push(slab( + &dims + .iter() + .map(|&n| (rng.below(n), 1, 1, 1)) + .collect::>(), + )); + // Distinct points in no particular order. + let mut points: Vec> = Vec::new(); + for _ in 0..1 + rng.below(15) { + let p: Vec = dims.iter().map(|&n| rng.below(n)).collect(); + if !points.contains(&p) { + points.push(p); + } + } + out.push(Selection::Points(points)); + out +} + +fn join(v: &[u64]) -> String { + v.iter().map(u64::to_string).collect::>().join(",") +} + +/// One element of type `code`, given as its raw file-order bytes, as an +/// integer (every value in these files is one). +fn decode(code: &str, big_endian: bool, bytes: &[u8]) -> i64 { + let mut b = bytes.to_vec(); + if big_endian { + b.reverse(); + } + let mut w = [0u8; 8]; + w[..b.len()].copy_from_slice(&b); + let u = u64::from_le_bytes(w); + match code { + "i1" => u as u8 as i8 as i64, + "i2" => u as u16 as i16 as i64, + "i4" => u as u32 as i32 as i64, + "i8" => u as i64, + "u1" | "u2" | "u4" | "u8" => u as i64, + "f2" => clawhdf5_format::float16::f16_bits_to_f32(u as u16) as i64, + "f4" => f32::from_bits(u as u32) as i64, + "f8" => f64::from_bits(u) as i64, + _ => unreachable!(), + } +} + +#[test] +fn selection_reads_match_h5py_for_every_type_order_and_rank() { + skip_if_no_python!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("contig.h5"); + write_file(&path); + + let mut rng = Rng(2026); + let mut cases: Vec<(String, String, Selection)> = Vec::new(); + for (name, code, shape) in datasets() { + for sel in selections(&mut rng, &shape) { + cases.push((name.clone(), code.clone(), sel)); + } + // Empty: a zero count, and Selection::None. + let mut empty = shape.iter().map(|&n| (0, 1, n, 1)).collect::>(); + empty[shape.len() - 1].2 = 0; + cases.push((name.clone(), code.clone(), slab(&empty))); + cases.push((name, code, Selection::None)); + } + let mut spec = String::new(); + for (k, (name, _, sel)) in cases.iter().enumerate() { + let line = match sel { + Selection::Hyperslab { count, .. } if count.contains(&0) => "N".to_string(), + Selection::Hyperslab { + start, + stride, + count, + block, + } => format!( + "H {};{};{};{}", + join(start), + join(stride), + join(count), + join(block) + ), + Selection::Points(points) => format!( + "P {}", + points.iter().map(|p| join(p)).collect::>().join(";") + ), + Selection::None => "N".to_string(), + Selection::All => unreachable!(), + }; + spec.push_str(&format!("{k} {name} {line}\n")); + } + let spec_path = dir.path().join("cases.txt"); + std::fs::write(&spec_path, spec).unwrap(); + run_python(&format!( + r#" +import h5py, numpy as np +with h5py.File("{path}", "r") as f: + for line in open("{spec}"): + k, name, kind, *rest = line.split() + d = f[name] + space = d.id.get_space() + if kind == 'H': + start, stride, count, block = (tuple(int(x) for x in part.split(',')) + for part in rest[0].split(';')) + space.select_hyperslab(start, count, stride, block) + elif kind == 'P': + pts = np.array([[int(x) for x in p.split(',')] for p in rest[0].split(';')], + dtype=np.uint64) + space.select_elements(pts) + else: + space.select_none() + n = space.get_select_npoints() + out = np.empty(n, dtype=d.dtype) + if n: + d.id.read(h5py.h5s.create_simple((n,)), space, out) + open("{dir}/sel_" + k + ".bin", "wb").write(out.tobytes()) +"#, + path = path.display(), + spec = spec_path.display(), + dir = dir.path().display(), + )); + + let file = File::open(&path).unwrap(); + for (k, (name, code, sel)) in cases.iter().enumerate() { + let ds = file.dataset(name).unwrap(); + let want_bytes = std::fs::read(dir.path().join(format!("sel_{k}.bin"))).unwrap(); + let got_bytes = ds.read_selection(sel).unwrap(); + assert!( + got_bytes == want_bytes, + "{name} {sel:?}: raw bytes differ from libhdf5's ({} vs {} bytes)", + got_bytes.len(), + want_bytes.len() + ); + let size = ds.raw_datatype().unwrap().type_size() as usize; + let want: Vec = want_bytes + .chunks_exact(size) + .map(|e| decode(code, name.contains("be_"), e)) + .collect(); + assert_eq!( + ds.read_f64_selection(sel).unwrap(), + want.iter().map(|&v| v as f64).collect::>(), + "{name} {sel:?} as f64" + ); + assert_eq!( + ds.read_f32_selection(sel).unwrap(), + want.iter().map(|&v| v as f32).collect::>(), + "{name} {sel:?} as f32" + ); + assert_eq!( + ds.read_i64_selection(sel).unwrap(), + want, + "{name} {sel:?} as i64" + ); + assert_eq!( + ds.read_i32_selection(sel).unwrap(), + want.iter().map(|&v| v as i32).collect::>(), + "{name} {sel:?} as i32" + ); + } +} diff --git a/docs/known-issues.md b/docs/known-issues.md index f8d7798..b5e3155 100644 --- a/docs/known-issues.md +++ b/docs/known-issues.md @@ -9,7 +9,8 @@ deleting it. ## Concurrent and contiguous read performance (measured 2026-09-26) -**Status:** open; one cause of the first bullet fixed (2026-09-26). Measured on +**Status:** open for chunked full reads (one cause fixed 2026-09-26); the +contiguous item is fixed (2026-09-26). Measured on tank with `concurrent_read` against h5py 3.16 / HDF5 2.0 (`BENCHMARKS.md`, "Concurrent reads"): - **Partly fixed 2026-09-26.** Full reads of chunked datasets from several threads @@ -34,6 +35,12 @@ tank with `concurrent_read` against h5py 3.16 / HDF5 2.0 (`BENCHMARKS.md`, the `f32` copy of it, and a new buffer per decoded chunk). - Contiguous datasets read 4x slower than h5py on one thread (2.5 vs 9.8 GB/s full, 0.12x for 256 x 256 hyperslabs). + **Fixed 2026-09-26** (not yet re-measured for `BENCHMARKS.md`): full + reads were dominated by 4 KiB page faults on the fresh output buffer, + which is now backed by transparent huge pages as numpy's is; hyperslab + reads copied the selection three times, element by element, and now copy + each contiguous run once, straight from the file into the output (see + `CHANGELOG.md`). The chunked-read scaling item above is still open. Values are correct; this is speed only. ## Silent wrong data found by the 2026-09-25 HDF5 audit