diff --git a/CHANGELOG.md b/CHANGELOG.md index 99f86f2..bb70d01 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,27 @@ `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 diff --git a/crates/clawhdf5-format/src/bulk_alloc.rs b/crates/clawhdf5-format/src/bulk_alloc.rs index 22a8bf2..085f125 100644 --- a/crates/clawhdf5-format/src/bulk_alloc.rs +++ b/crates/clawhdf5-format/src/bulk_alloc.rs @@ -19,6 +19,7 @@ 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, diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index 1c53418..633ecd7 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -531,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()) @@ -540,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) } @@ -757,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()); + 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); @@ -781,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 @@ -799,9 +840,8 @@ 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); @@ -941,19 +981,8 @@ 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); @@ -986,6 +1015,12 @@ 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 = crate::bulk_alloc::vec_for_bulk(count); for i in 0..count { @@ -1013,9 +1048,8 @@ 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) { @@ -1103,19 +1137,8 @@ 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); 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 05ebd46..23840b8 100644 --- a/crates/clawhdf5-format/src/lib.rs +++ b/crates/clawhdf5-format/src/lib.rs @@ -94,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/docs/known-issues.md b/docs/known-issues.md index bf8f77a..60fc9d8 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. Measured on tank with `concurrent_read` against h5py +**Status:** open for chunked full reads; 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"): - Full reads of chunked datasets from several threads through one `File` stop scaling at about 4 threads (880 MB/s on deflate data vs 4424 MB/s @@ -17,6 +18,12 @@ deleting it. scale to 1244 MB/s, so the `File`'s shared chunk cache is the suspect. - 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