diff --git a/crates/clawhdf5-format/src/chunk_grid.rs b/crates/clawhdf5-format/src/chunk_grid.rs new file mode 100644 index 0000000..9e03b9d --- /dev/null +++ b/crates/clawhdf5-format/src/chunk_grid.rs @@ -0,0 +1,200 @@ +//! Chunk-index linearisation shared by the Fixed Array and Extensible Array +//! chunk indexes (reader and writer). +//! +//! Both indexes store one element per chunk at a *linear* index, and the +//! library derives that index from the chunk's scaled coordinates +//! (`offset / chunk_dim`) using the dataset's **maximum** dimensions, not its +//! current ones (`H5D__farray_idx_get_addr` / `H5D__earray_idx_get_addr`, +//! via `layout->max_down_chunks`). A dataset whose current shape is smaller +//! than its maxshape therefore has gaps in the index, and laying it out by the +//! current shape puts every chunk after the first row in the wrong place. +//! +//! The Extensible Array adds one more step: its one unlimited dimension has no +//! finite chunk count, so the library *swizzles* the coordinates to make that +//! dimension the slowest-varying one (`H5VM_swizzle_coords`, which moves +//! `coords[unlim_dim]` to the front and shifts the dimensions before it right +//! by one) before linearising with `swizzled_max_down_chunks`. When the +//! unlimited dimension is already dimension 0 no swizzle happens. + +#[cfg(not(feature = "std"))] +extern crate alloc; + +#[cfg(not(feature = "std"))] +use alloc::{vec, vec::Vec}; + +use crate::error::FormatError; + +/// How a chunk index maps linear element indexes to chunk coordinates. +#[derive(Debug, Clone)] +pub(crate) struct ChunkGrid { + /// Spatial chunk dimensions, in dataset order. + chunk_dims: Vec, + /// Chunks per dimension covering the *current* extent, in dataset order. + cur_chunks: Vec, + /// Dataset dimension stored at each linearisation position (slowest + /// first). The identity except for a swizzled Extensible Array. + order: Vec, + /// Linear stride of each linearisation position. + down: Vec, +} + +impl ChunkGrid { + /// Grid for a Fixed Array index: row-major over the chunk counts of the + /// maximum dimensions (`max_dims`, falling back to the current dimensions + /// when the dataspace records none). + pub(crate) fn fixed_array( + cur_dims: &[u64], + max_dims: Option<&[u64]>, + chunk_dims: &[u64], + ) -> Result { + Self::build(cur_dims, max_dims, chunk_dims, None) + } + + /// Grid for an Extensible Array index: like the Fixed Array, but the + /// unlimited dimension (the one whose maximum is `H5S_UNLIMITED`) is moved + /// to the slowest-varying position first. + pub(crate) fn extensible_array( + cur_dims: &[u64], + max_dims: Option<&[u64]>, + chunk_dims: &[u64], + ) -> Result { + let unlim = max_dims.and_then(|m| m.iter().position(|&d| d == u64::MAX)); + Self::build(cur_dims, max_dims, chunk_dims, unlim) + } + + fn build( + cur_dims: &[u64], + max_dims: Option<&[u64]>, + chunk_dims: &[u64], + unlim: Option, + ) -> Result { + let rank = chunk_dims.len(); + if cur_dims.len() != rank || max_dims.is_some_and(|m| m.len() != rank) { + return Err(FormatError::ChunkedReadError( + "chunk index rank does not match the dataspace".into(), + )); + } + if chunk_dims.contains(&0) { + return Err(FormatError::ChunkedReadError( + "chunk dimension is zero".into(), + )); + } + let cur_chunks: Vec = cur_dims + .iter() + .zip(chunk_dims) + .map(|(&d, &c)| d.div_ceil(c)) + .collect(); + // Chunk counts of the maximum extent. An unlimited dimension has no + // finite count; it only ever sits in the slowest position, where its + // count never enters a stride. A (corrupt) maximum smaller than the + // current extent is widened so no allocated chunk becomes unreachable. + let max_chunks: Vec = (0..rank) + .map(|d| { + let max = max_dims.map_or(cur_dims[d], |m| m[d]); + if max == u64::MAX { + u64::MAX + } else { + max.div_ceil(chunk_dims[d]).max(cur_chunks[d]) + } + }) + .collect(); + + let mut order: Vec = (0..rank).collect(); + if let Some(u) = unlim { + order.remove(u); + order.insert(0, u); + } + let mut down = vec![1u64; rank]; + for p in (0..rank.saturating_sub(1)).rev() { + let next = max_chunks[order[p + 1]]; + if next == u64::MAX { + // Only reachable with more than one unlimited dimension, which + // neither index type can describe. + return Err(FormatError::ChunkedReadError( + "array chunk index with more than one unlimited dimension".into(), + )); + } + down[p] = down[p + 1].checked_mul(next).ok_or_else(|| { + FormatError::Overflow("chunk index linear stride overflows u64".into()) + })?; + } + Ok(Self { + chunk_dims: chunk_dims.to_vec(), + cur_chunks, + order, + down, + }) + } + + /// Dataset-space offsets of the chunk stored at linear `index`, or `None` + /// when that chunk lies outside the current extent (the index still has a + /// slot for it; the library ignores such chunks on read). + pub(crate) fn offsets(&self, index: u64) -> Option> { + let rank = self.chunk_dims.len(); + let mut offsets = vec![0u64; rank]; + let mut rem = index; + for p in 0..rank { + let d = self.order[p]; + let scaled = rem / self.down[p]; + rem %= self.down[p]; + if scaled >= self.cur_chunks[d] { + return None; + } + offsets[d] = scaled * self.chunk_dims[d]; + } + Some(offsets) + } + + /// Linear index of the chunk with scaled coordinates `scaled` + /// (`offset / chunk_dim` per dimension, in dataset order). + pub(crate) fn linear_index(&self, scaled: &[u64]) -> u64 { + self.order + .iter() + .zip(&self.down) + .map(|(&d, &stride)| scaled[d] * stride) + .sum() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn fixed_array_uses_max_dims() { + // shape (4, 6), chunks (2, 3), maxshape (20, 10): 10 x 4 chunk grid. + let g = ChunkGrid::fixed_array(&[4, 6], Some(&[20, 10]), &[2, 3]).unwrap(); + assert_eq!(g.offsets(0), Some(vec![0, 0])); + assert_eq!(g.offsets(1), Some(vec![0, 3])); + assert_eq!(g.offsets(2), None); // column chunk 2 is beyond the extent + assert_eq!(g.offsets(4), Some(vec![2, 0])); + assert_eq!(g.offsets(5), Some(vec![2, 3])); + assert_eq!(g.offsets(8), None); // row chunk 2 is beyond the extent + assert_eq!(g.linear_index(&[1, 1]), 5); + } + + #[test] + fn extensible_array_swizzles_unlimited_dim() { + // maxshape (10, None): dim 1 is unlimited and becomes slowest. + let g = ChunkGrid::extensible_array(&[4, 6], Some(&[10, u64::MAX]), &[2, 3]).unwrap(); + // max chunks of dim 0 = 5, so index = c1 * 5 + c0. + assert_eq!(g.linear_index(&[1, 0]), 1); + assert_eq!(g.linear_index(&[0, 1]), 5); + assert_eq!(g.offsets(5), Some(vec![0, 3])); + assert_eq!(g.offsets(6), Some(vec![2, 3])); + assert_eq!(g.offsets(2), None); + } + + #[test] + fn extensible_array_unlimited_first_is_row_major() { + let g = ChunkGrid::extensible_array(&[4, 6], Some(&[u64::MAX, 30]), &[2, 3]).unwrap(); + // max chunks of dim 1 = 10. + assert_eq!(g.linear_index(&[1, 1]), 11); + assert_eq!(g.offsets(11), Some(vec![2, 3])); + } + + #[test] + fn rejects_two_unlimited_dims_after_the_first() { + assert!(ChunkGrid::fixed_array(&[4, 6], Some(&[u64::MAX, u64::MAX]), &[2, 3]).is_err()); + } +} diff --git a/crates/clawhdf5-format/src/chunked_read.rs b/crates/clawhdf5-format/src/chunked_read.rs index 8dade54..f07ae27 100644 --- a/crates/clawhdf5-format/src/chunked_read.rs +++ b/crates/clawhdf5-format/src/chunked_read.rs @@ -593,6 +593,7 @@ pub fn list_chunks( file_data, &header, &dataspace.dimensions, + dataspace.max_dimensions.as_deref(), spatial_chunk_dims, elem_size as u32, offset_size, @@ -608,6 +609,7 @@ pub fn list_chunks( file_data, &header, &dataspace.dimensions, + dataspace.max_dimensions.as_deref(), spatial_chunk_dims, elem_size as u32, offset_size, diff --git a/crates/clawhdf5-format/src/chunked_write.rs b/crates/clawhdf5-format/src/chunked_write.rs index 6816546..682a1c6 100644 --- a/crates/clawhdf5-format/src/chunked_write.rs +++ b/crates/clawhdf5-format/src/chunked_write.rs @@ -4,10 +4,11 @@ extern crate alloc; #[cfg(not(feature = "std"))] -use alloc::{vec, vec::Vec}; +use alloc::{format, vec, vec::Vec}; use crate::checksum::jenkins_lookup3; use crate::chunk_cache::{CACHE_LINE_SIZE, align_to_cache_line}; +use crate::chunk_grid::ChunkGrid; use crate::ea_writer; use crate::error::FormatError; use crate::filter_pipeline::{ @@ -443,6 +444,27 @@ fn serialize_v4_fixed_array( element_size: u32, max_bits: u8, ) -> Vec { + let mut buf = layout_v4_chunked_prefix(chunk_dims, element_size); + + // chunk index type = 3 (Fixed Array) + buf.push(3); + + // max_dblk_page_nelmts_bits — must match FAHD max_nelmts_bits + buf.push(max_bits); + + // Fixed Array header address + match offset_size { + 4 => buf.extend_from_slice(&(fixed_array_address as u32).to_le_bytes()), + 8 => buf.extend_from_slice(&fixed_array_address.to_le_bytes()), + _ => {} + } + + buf +} + +/// The part of a v4 chunked layout message before the chunk index type: +/// version, class, flags and the chunk dimensions (plus the element size). +fn layout_v4_chunked_prefix(chunk_dims: &[u32], element_size: u32) -> Vec { let mut buf = Vec::new(); buf.push(4); // version buf.push(2); // class = chunked @@ -482,125 +504,143 @@ fn serialize_v4_fixed_array( 4 => buf.extend_from_slice(&element_size.to_le_bytes()), _ => {} } - - // chunk index type = 3 (Fixed Array) - buf.push(3); - - // max_dblk_page_nelmts_bits — must match FAHD max_nelmts_bits - buf.push(max_bits); - - // Fixed Array header address - match offset_size { - 4 => buf.extend_from_slice(&(fixed_array_address as u32).to_le_bytes()), - 8 => buf.extend_from_slice(&fixed_array_address.to_le_bytes()), - _ => {} - } - buf } +/// log2 of the elements per Fixed Array data block page (the library's +/// default, `H5D_FARRAY_MAX_DBLK_PAGE_NELMTS_BITS`). +const FA_PAGE_BITS: u8 = 10; + +pub(crate) fn push_addr(buf: &mut Vec, addr: u64, offset_size: u8) { + match offset_size { + 4 => buf.extend_from_slice(&(addr as u32).to_le_bytes()), + _ => buf.extend_from_slice(&addr.to_le_bytes()), + } +} + +/// Width of the chunk-size field of a filtered chunk index element. Must +/// match the library's `H5D_FARRAY_FILT_COMPUTE_CHUNK_SIZE_LEN` (the EA and +/// B-tree v2 indexes use the same formula): +/// `1 + ((log2(unfiltered chunk bytes) + 8) / 8)`, capped at 8. +pub(crate) fn filtered_chunk_size_len(slots: &[Option]) -> usize { + let max_raw = slots + .iter() + .flatten() + .map(|c| c.raw_size) + .max() + .unwrap_or(1); + let log2_val = if max_raw <= 1 { + 0 + } else { + 63 - max_raw.leading_zeros() + }; + (1 + ((log2_val + 8) / 8) as usize).min(8) +} + +/// Append one chunk index element: the chunk's address, plus its stored size +/// and filter mask when the dataset is filtered. `None` is an unallocated +/// chunk (undefined address, zero size and mask). +pub(crate) fn push_index_element( + buf: &mut Vec, + slot: Option<&WrittenChunk>, + offset_size: u8, + chunk_size_bytes: Option, +) { + match slot { + Some(c) => { + push_addr(buf, c.address, offset_size); + if let Some(n) = chunk_size_bytes { + buf.extend_from_slice(&c.compressed_size.to_le_bytes()[..n]); + buf.extend_from_slice(&c.filter_mask.to_le_bytes()); + } + } + None => { + buf.extend(core::iter::repeat_n(0xFF, offset_size as usize)); + if let Some(n) = chunk_size_bytes { + buf.extend(core::iter::repeat_n(0x00, n + 4)); + } + } + } +} + /// Build a complete Fixed Array at a known absolute address. +/// +/// `slots` holds one entry per element of the array, i.e. per chunk of the +/// dataset's *maximum* extent in the order [`crate::chunk_grid`] defines; +/// `None` marks a chunk that is not allocated. An array with more elements +/// than fit in one page (`2^FA_PAGE_BITS`) gets a paged data block: a +/// page-init bitmap after the prefix, then one checksummed page per +/// `2^FA_PAGE_BITS` elements, the last one short (`H5FA__dblock_create`). pub fn build_fixed_array_at( - chunks: &[WrittenChunk], + slots: &[Option], offset_size: u8, length_size: u8, has_filters: bool, fa_base_address: u64, ) -> Vec { let os = offset_size as usize; - let num_elements = chunks.len(); - - // For filtered chunks, compute chunk_size encoding width. - // Must match the HDF5 C library's H5D_FARRAY_FILT_COMPUTE_CHUNK_SIZE_LEN macro: - // chunk_size_len = 1 + ((H5VM_log2_gen(chunk.size) + 8) / 8) - // where chunk.size is the unfiltered chunk size in bytes (product of all chunk dims). - let chunk_size_bytes: usize = if has_filters { - let max_raw = chunks.iter().map(|c| c.raw_size).max().unwrap_or(1); - let log2_val = if max_raw <= 1 { - 0 - } else { - 63 - max_raw.leading_zeros() - }; - let len = 1 + ((log2_val + 8) / 8) as usize; - len.min(8) - } else { - 0 - }; - - let elem_size = if has_filters { - os + chunk_size_bytes + 4 - } else { - os - }; + let num_elements = slots.len(); + let chunk_size_bytes = has_filters.then(|| filtered_chunk_size_len(slots)); + let elem_size = os + chunk_size_bytes.map_or(0, |n| n + 4); let client_id: u8 = if has_filters { 1 } else { 0 }; // FAHD total size - let nelmts_field_size = length_size as usize; - let fahd_total_size = 4 + 1 + 1 + 1 + 1 + nelmts_field_size + os + 4; + let fahd_total_size = 4 + 1 + 1 + 1 + 1 + length_size as usize + os + 4; let fadb_address = fa_base_address + fahd_total_size as u64; - // Build FAHD let mut fahd = Vec::with_capacity(fahd_total_size); fahd.extend_from_slice(b"FAHD"); fahd.push(0); // version fahd.push(client_id); fahd.push(elem_size as u8); - - // max_nelmts_bits: use 10 as default (page_size = 1024), matching h5py convention - let max_bits: u8 = 10; - fahd.push(max_bits); - + fahd.push(FA_PAGE_BITS); match length_size { 4 => fahd.extend_from_slice(&(num_elements as u32).to_le_bytes()), - 8 => fahd.extend_from_slice(&(num_elements as u64).to_le_bytes()), _ => fahd.extend_from_slice(&(num_elements as u64).to_le_bytes()), } - - match offset_size { - 4 => fahd.extend_from_slice(&(fadb_address as u32).to_le_bytes()), - 8 => fahd.extend_from_slice(&fadb_address.to_le_bytes()), - _ => fahd.extend_from_slice(&fadb_address.to_le_bytes()), - } - - // Checksum + push_addr(&mut fahd, fadb_address, offset_size); let checksum = jenkins_lookup3(&fahd); fahd.extend_from_slice(&checksum.to_le_bytes()); - assert_eq!(fahd.len(), fahd_total_size); - // Build FADB + // FADB prefix let mut fadb = Vec::new(); fadb.extend_from_slice(b"FADB"); fadb.push(0); // version fadb.push(client_id); + push_addr(&mut fadb, fa_base_address, offset_size); - // header address - match offset_size { - 4 => fadb.extend_from_slice(&(fa_base_address as u32).to_le_bytes()), - 8 => fadb.extend_from_slice(&fa_base_address.to_le_bytes()), - _ => fadb.extend_from_slice(&fa_base_address.to_le_bytes()), - } - - // Element data - for chunk in chunks { - match offset_size { - 4 => fadb.extend_from_slice(&(chunk.address as u32).to_le_bytes()), - 8 => fadb.extend_from_slice(&chunk.address.to_le_bytes()), - _ => fadb.extend_from_slice(&chunk.address.to_le_bytes()), + let page_nelmts = 1usize << FA_PAGE_BITS; + if num_elements <= page_nelmts { + // Unpaged: the elements follow the prefix, one checksum over both. + for slot in slots { + push_index_element(&mut fadb, slot.as_ref(), offset_size, chunk_size_bytes); } - if has_filters { - // Write compressed size using chunk_size_bytes (variable width) - let cs_bytes = chunk.compressed_size.to_le_bytes(); - fadb.extend_from_slice(&cs_bytes[..chunk_size_bytes]); - fadb.extend_from_slice(&chunk.filter_mask.to_le_bytes()); + let fadb_checksum = jenkins_lookup3(&fadb); + fadb.extend_from_slice(&fadb_checksum.to_le_bytes()); + } else { + // Paged: every page is written, so every page-init bit is set + // (MSB-first, as `H5VM_bit_set` packs them). The prefix and bitmap + // share a checksum; each page carries its own. + let npages = num_elements.div_ceil(page_nelmts); + let mut bitmap = vec![0u8; npages.div_ceil(8)]; + for p in 0..npages { + bitmap[p / 8] |= 0x80 >> (p % 8); + } + fadb.extend_from_slice(&bitmap); + let prefix_checksum = jenkins_lookup3(&fadb); + fadb.extend_from_slice(&prefix_checksum.to_le_bytes()); + for page in slots.chunks(page_nelmts) { + let start = fadb.len(); + for slot in page { + push_index_element(&mut fadb, slot.as_ref(), offset_size, chunk_size_bytes); + } + let page_checksum = jenkins_lookup3(&fadb[start..]); + fadb.extend_from_slice(&page_checksum.to_le_bytes()); } } - // FADB checksum - let fadb_checksum = jenkins_lookup3(&fadb); - fadb.extend_from_slice(&fadb_checksum.to_le_bytes()); - let mut combined = fahd; combined.extend_from_slice(&fadb); combined @@ -667,7 +707,8 @@ pub fn build_chunked_data_from_precompressed( pre: &PrecompressedChunks, base_address: u64, maxshape: Option<&[u64]>, -) -> ChunkedDataResult { +) -> Result { + let index = ChunkIndexPlan::new(&pre.shape, maxshape, &pre.chunk_dims)?; let offset_size: u8 = 8; let length_size: u8 = 8; let num_chunks = pre.chunks.len(); @@ -693,71 +734,329 @@ pub fn build_chunked_data_from_precompressed( } let chunk_dims_u32: Vec = pre.chunk_dims.iter().map(|&d| d as u32).collect(); - let use_extensible = maxshape.is_some_and(|ms| ms.contains(&u64::MAX)); let aligned_idx = align_to_cache_line(data_buf.len()); if aligned_idx > data_buf.len() { data_buf.resize(aligned_idx, 0u8); } - let layout_message = if use_extensible { - let ea_address = base_address + data_buf.len() as u64; - let ea_bytes = ea_writer::build_extensible_array_at( - &written_chunks, - offset_size, - length_size, - pre.has_filters, - ea_address, - ); - data_buf.extend_from_slice(&ea_bytes); - ea_writer::serialize_v4_extensible_array( - &chunk_dims_u32, - ea_address, - offset_size, - element_size as u32, - ) - } else if num_chunks == 1 { - let chunk_addr = written_chunks[0].address; - let filtered_size = if pre.has_filters { - Some(written_chunks[0].compressed_size) - } else { - None - }; - let filter_mask = if pre.has_filters { Some(0u32) } else { None }; - serialize_v4_single_chunk( - &chunk_dims_u32, - chunk_addr, - filtered_size, - filter_mask, - offset_size, - element_size as u32, - ) - } else { - let fa_address = base_address + data_buf.len() as u64; - let fa_bytes = build_fixed_array_at( - &written_chunks, - offset_size, - length_size, - pre.has_filters, - fa_address, - ); - data_buf.extend_from_slice(&fa_bytes); - serialize_v4_fixed_array( - &chunk_dims_u32, - fa_address, - offset_size, - element_size as u32, - 10, // max_nelmts_bits — matches h5py convention - ) + let layout_message = match &index { + ChunkIndexPlan::ExtensibleArray(grid) => { + let ea_address = base_address + data_buf.len() as u64; + let slots = index_slots(grid, &pre.shape, &pre.chunk_dims, &written_chunks, None)?; + let ea_bytes = ea_writer::build_extensible_array_at( + &slots, + offset_size, + length_size, + pre.has_filters, + ea_address, + ); + data_buf.extend_from_slice(&ea_bytes); + ea_writer::serialize_v4_extensible_array( + &chunk_dims_u32, + ea_address, + offset_size, + element_size as u32, + ) + } + ChunkIndexPlan::SingleChunk => { + let chunk_addr = written_chunks[0].address; + let filtered_size = if pre.has_filters { + Some(written_chunks[0].compressed_size) + } else { + None + }; + let filter_mask = if pre.has_filters { Some(0u32) } else { None }; + serialize_v4_single_chunk( + &chunk_dims_u32, + chunk_addr, + filtered_size, + filter_mask, + offset_size, + element_size as u32, + ) + } + ChunkIndexPlan::FixedArray(grid, nslots) => { + let fa_address = base_address + data_buf.len() as u64; + let slots = index_slots( + grid, + &pre.shape, + &pre.chunk_dims, + &written_chunks, + Some(*nslots), + )?; + let fa_bytes = build_fixed_array_at( + &slots, + offset_size, + length_size, + pre.has_filters, + fa_address, + ); + data_buf.extend_from_slice(&fa_bytes); + serialize_v4_fixed_array( + &chunk_dims_u32, + fa_address, + offset_size, + element_size as u32, + FA_PAGE_BITS, + ) + } + ChunkIndexPlan::BTreeV2 => { + let bt_address = base_address + data_buf.len() as u64; + let records: Vec<(Vec, &WrittenChunk)> = written_chunks + .iter() + .enumerate() + .map(|(i, c)| (scaled_coords(&pre.shape, &pre.chunk_dims, i), c)) + .collect(); + let (bt_bytes, node_size) = build_btree_v2_chunk_index_at( + pre.shape.len(), + &records, + offset_size, + length_size, + pre.has_filters, + bt_address, + )?; + data_buf.extend_from_slice(&bt_bytes); + serialize_v4_btree_v2( + &chunk_dims_u32, + bt_address, + offset_size, + element_size as u32, + node_size, + ) + } }; - ChunkedDataResult { + Ok(ChunkedDataResult { data_bytes: data_buf, layout_message, pipeline_message: pre.pipeline_message.clone(), + }) +} + +/// Most slots a Fixed Array index may have before we refuse to build it: its +/// data block holds one element per chunk of the *maximum* extent, so a huge +/// finite maxshape with small chunks would otherwise exhaust memory. +const MAX_FIXED_ARRAY_SLOTS: u64 = 1 << 26; + +/// Which chunk index a dataset gets, following the library's choice in +/// `H5D__layout_set_latest_indexing`: version-2 B-tree for more than one +/// unlimited dimension, Extensible Array for exactly one, Fixed Array for a +/// finite maxshape, Single Chunk when the whole maximum extent is one chunk. +enum ChunkIndexPlan { + SingleChunk, + /// The grid and the number of array elements (chunks of the max extent). + FixedArray(ChunkGrid, usize), + ExtensibleArray(ChunkGrid), + BTreeV2, +} + +impl ChunkIndexPlan { + fn new( + shape: &[u64], + maxshape: Option<&[u64]>, + chunk_dims: &[u64], + ) -> Result { + let bad = |what: &str| FormatError::ChunkedReadError(format!("maxshape: {what}")); + if let Some(ms) = maxshape { + if ms.len() != shape.len() { + return Err(bad("rank differs from the shape")); + } + if ms.iter().zip(shape).any(|(&m, &s)| m < s) { + return Err(bad("smaller than the shape")); + } + } + let max = maxshape.unwrap_or(shape); + let nunlim = max.iter().filter(|&&d| d == u64::MAX).count(); + match nunlim { + 0 => { + let nslots = max + .iter() + .zip(chunk_dims) + .try_fold(1u64, |acc, (&m, &c)| acc.checked_mul(m.div_ceil(c.max(1)))) + .filter(|&n| n <= MAX_FIXED_ARRAY_SLOTS) + .ok_or_else(|| { + bad("too many chunks for a Fixed Array index; \ + use larger chunks or an unlimited dimension") + })?; + // A Single Chunk index needs that one chunk to exist; an + // empty dataset gets an all-unallocated Fixed Array instead. + let empty = shape.contains(&0); + if nslots == 1 && !empty { + Ok(Self::SingleChunk) + } else { + let grid = ChunkGrid::fixed_array(shape, Some(max), chunk_dims)?; + Ok(Self::FixedArray(grid, nslots as usize)) + } + } + 1 => Ok(Self::ExtensibleArray(ChunkGrid::extensible_array( + shape, + Some(max), + chunk_dims, + )?)), + _ => Ok(Self::BTreeV2), + } } } +/// Place each written chunk at its linear index in `grid`. `chunks` are in +/// row-major order over the chunks of the current extent (`split_into_chunks`). +/// `len` fixes the slot count (Fixed Array); otherwise it is one past the +/// highest index used. +fn index_slots( + grid: &ChunkGrid, + shape: &[u64], + chunk_dims: &[u64], + chunks: &[WrittenChunk], + len: Option, +) -> Result>, FormatError> { + let mut placed: Vec<(usize, &WrittenChunk)> = Vec::with_capacity(chunks.len()); + for (i, chunk) in chunks.iter().enumerate() { + let scaled = scaled_coords(shape, chunk_dims, i); + let idx = usize::try_from(grid.linear_index(&scaled)) + .map_err(|_| FormatError::Overflow("chunk index slot".into()))?; + placed.push((idx, chunk)); + } + let n = len.unwrap_or_else(|| placed.iter().map(|&(i, _)| i + 1).max().unwrap_or(0)); + let mut slots = vec![None; n]; + for (idx, chunk) in placed { + *slots + .get_mut(idx) + .ok_or_else(|| FormatError::Overflow("chunk index slot".into()))? = Some(chunk.clone()); + } + Ok(slots) +} + +/// Scaled coordinates (`offset / chunk_dim`) of the `i`-th chunk in the +/// row-major order `split_into_chunks` produces over the current extent. +fn scaled_coords(shape: &[u64], chunk_dims: &[u64], i: usize) -> Vec { + let rank = shape.len(); + let mut scaled = vec![0u64; rank]; + let mut rem = i as u64; + for d in (0..rank).rev() { + let n = shape[d].div_ceil(chunk_dims[d]); + scaled[d] = rem % n; + rem /= n; + } + scaled +} + +/// Node size the library gives a chunk index B-tree (`H5D_BT2_NODE_SIZE`), +/// with its split and merge percentages. +const BT2_NODE_SIZE: u32 = 2048; +const BT2_SPLIT_PERCENT: u8 = 100; +const BT2_MERGE_PERCENT: u8 = 40; +/// B-tree v2 record types for chunk indexes (`H5B2_CDSET_ID`, +/// `H5B2_CDSET_FILT_ID`). +const BT2_CHUNK_UNFILTERED: u8 = 10; +const BT2_CHUNK_FILTERED: u8 = 11; + +/// Build a version-2 B-tree chunk index (the library's index for datasets +/// with more than one unlimited dimension) at a known absolute address. +/// +/// `records` are `(scaled coordinates, chunk)` in lexicographic order of the +/// coordinates, which is the order the library's comparator +/// (`H5VM_vector_cmp_u`) keeps them in. The tree is a single leaf: the +/// library's 2048-byte node when the records fit, otherwise a leaf node +/// sized to hold them all (the root's record count is 16-bit, so at most +/// 65535 chunks). Returns the bytes and the node size the layout message +/// must record. +fn build_btree_v2_chunk_index_at( + rank: usize, + records: &[(Vec, &WrittenChunk)], + offset_size: u8, + length_size: u8, + has_filters: bool, + base_address: u64, +) -> Result<(Vec, u32), FormatError> { + let os = offset_size as usize; + let nrec = u16::try_from(records.len()).map_err(|_| { + FormatError::ChunkedReadError( + "more than 65535 chunks with more than one unlimited dimension: \ + use larger chunks" + .into(), + ) + })?; + let chunk_size_bytes = has_filters.then(|| { + let slots: Vec> = + records.iter().map(|(_, c)| Some((*c).clone())).collect(); + filtered_chunk_size_len(&slots) + }); + let record_size = os + chunk_size_bytes.map_or(0, |n| n + 4) + 8 * rank; + // Leaf: signature, version, type, records, checksum. + let leaf_len = 4 + 1 + 1 + records.len() * record_size + 4; + let node_size = u32::try_from(leaf_len) + .map_err(|_| FormatError::Overflow("B-tree v2 leaf size".into()))? + .max(BT2_NODE_SIZE); + let tree_type = if has_filters { + BT2_CHUNK_FILTERED + } else { + BT2_CHUNK_UNFILTERED + }; + + let hdr_len = 4 + 1 + 1 + 4 + 2 + 2 + 1 + 1 + os + 2 + length_size as usize + 4; + let leaf_address = base_address + hdr_len as u64; + + let mut out = Vec::with_capacity(hdr_len + node_size as usize); + out.extend_from_slice(b"BTHD"); + out.push(0); // version + out.push(tree_type); + out.extend_from_slice(&node_size.to_le_bytes()); + out.extend_from_slice(&(record_size as u16).to_le_bytes()); + out.extend_from_slice(&0u16.to_le_bytes()); // depth + out.push(BT2_SPLIT_PERCENT); + out.push(BT2_MERGE_PERCENT); + if records.is_empty() { + out.extend(core::iter::repeat_n(0xFF, os)); + } else { + push_addr(&mut out, leaf_address, offset_size); + } + out.extend_from_slice(&nrec.to_le_bytes()); + match length_size { + 4 => out.extend_from_slice(&(records.len() as u32).to_le_bytes()), + _ => out.extend_from_slice(&(records.len() as u64).to_le_bytes()), + } + let sum = jenkins_lookup3(&out); + out.extend_from_slice(&sum.to_le_bytes()); + debug_assert_eq!(out.len(), hdr_len); + if records.is_empty() { + return Ok((out, node_size)); + } + + let leaf_start = out.len(); + out.extend_from_slice(b"BTLF"); + out.push(0); // version + out.push(tree_type); + for (scaled, chunk) in records { + push_index_element(&mut out, Some(chunk), offset_size, chunk_size_bytes); + for &c in scaled { + out.extend_from_slice(&c.to_le_bytes()); + } + } + let sum = jenkins_lookup3(&out[leaf_start..]); + out.extend_from_slice(&sum.to_le_bytes()); + // The library reads whole nodes; pad the leaf out to the node size. + out.resize(leaf_start + node_size as usize, 0); + Ok((out, node_size)) +} + +/// Serialize a v4 layout message for a version-2 B-tree chunk index. +fn serialize_v4_btree_v2( + chunk_dims: &[u32], + btree_address: u64, + offset_size: u8, + element_size: u32, + node_size: u32, +) -> Vec { + let mut buf = layout_v4_chunked_prefix(chunk_dims, element_size); + buf.push(5); // chunk index type = 5 (version-2 B-tree) + buf.extend_from_slice(&node_size.to_le_bytes()); + buf.push(BT2_SPLIT_PERCENT); + buf.push(BT2_MERGE_PERCENT); + push_addr(&mut buf, btree_address, offset_size); + buf +} + /// Build chunked data with absolute addresses. /// If `maxshape` has unlimited dims, uses Extensible Array index. pub fn build_chunked_data_at( @@ -790,11 +1089,7 @@ pub fn build_chunked_data_at_ext( maxshape: Option<&[u64]>, ) -> Result { let pre = precompress_chunks(raw_data, shape, chunk_dims, element_size, options)?; - Ok(build_chunked_data_from_precompressed( - &pre, - base_address, - maxshape, - )) + build_chunked_data_from_precompressed(&pre, base_address, maxshape) } /// Write selected elements into an existing in-memory dataset buffer. @@ -1382,7 +1677,8 @@ mod tests { filter_mask: 0, }, ]; - let fa = build_fixed_array_at(&chunks, 8, 8, false, 0x2000); + let slots: Vec<_> = chunks.into_iter().map(Some).collect(); + let fa = build_fixed_array_at(&slots, 8, 8, false, 0x2000); // Should start with FAHD assert_eq!(&fa[0..4], b"FAHD"); // FAHD size = 4+1+1+1+1+8+8+4 = 28 @@ -1429,7 +1725,8 @@ mod tests { filter_mask: 0, }, ]; - let ea = ea_writer::build_extensible_array_at(&chunks, 8, 8, false, 0x2000); + let slots: Vec<_> = chunks.into_iter().map(Some).collect(); + let ea = ea_writer::build_extensible_array_at(&slots, 8, 8, false, 0x2000); assert_eq!(&ea[0..4], b"EAHD"); // Find EAIB after EAHD: 12 fixed + 6*8 stats + 8 addr + 4 checksum = 72 let aehd_size = 4 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 6 * 8 + 8 + 4; diff --git a/crates/clawhdf5-format/src/ea_writer.rs b/crates/clawhdf5-format/src/ea_writer.rs index d92fca2..e0f0375 100644 --- a/crates/clawhdf5-format/src/ea_writer.rs +++ b/crates/clawhdf5-format/src/ea_writer.rs @@ -7,7 +7,7 @@ extern crate alloc; use alloc::{vec, vec::Vec}; use crate::checksum::jenkins_lookup3; -use crate::chunked_write::WrittenChunk; +use crate::chunked_write::{WrittenChunk, filtered_chunk_size_len, push_addr, push_index_element}; /// Serialize a v4 Extensible Array layout message. pub(crate) fn serialize_v4_extensible_array( @@ -58,11 +58,11 @@ pub(crate) fn serialize_v4_extensible_array( buf.push(4); // EA creation parameters (must match AEHD and HDF5 C library defaults) - buf.push(32); // max_nelmts_bits - buf.push(4); // idx_blk_elmts - buf.push(4); // super_blk_min_data_ptrs - buf.push(16); // data_blk_min_elmts - buf.push(10); // max_dblk_page_nelmts_bits + buf.push(MAX_NELMTS_BITS); + buf.push(IDX_BLK_ELMTS); + buf.push(SUP_BLK_MIN_DATA_PTRS); + buf.push(DATA_BLK_MIN_ELMTS); + buf.push(MAX_DBLK_PAGE_NELMTS_BITS); // EA header address match offset_size { @@ -74,304 +74,281 @@ pub(crate) fn serialize_v4_extensible_array( buf } +// EA creation parameters — the HDF5 library's defaults for chunk indexes +// (`H5D_EARRAY_*`); the layout message above and the header must agree. +const MAX_NELMTS_BITS: u8 = 32; +const IDX_BLK_ELMTS: u8 = 4; +const SUP_BLK_MIN_DATA_PTRS: u8 = 4; +const DATA_BLK_MIN_ELMTS: u8 = 16; +const MAX_DBLK_PAGE_NELMTS_BITS: u8 = 10; + +/// One data block of the array: its first element (relative to the end of +/// the index block's own elements), element count, and address when it is +/// allocated. +struct DataBlock { + start: usize, + nelmts: usize, + addr: Option, +} + /// Build a complete Extensible Array at a known absolute address. /// -/// For simplicity, we put all elements inline in the index block when the -/// number of chunks is small (up to idx_blk_elmts), otherwise use inline + -/// direct data blocks. +/// `slots[i]` is the element at linear index `i` (see `chunk_grid`); `None` +/// marks an unallocated chunk. The first `IDX_BLK_ELMTS` elements live in +/// the index block, the rest in data blocks grouped by super block level +/// exactly as `H5EA__hdr_init` sizes them: level `u` has `2^(u/2)` data +/// blocks of `DATA_BLK_MIN_ELMTS * 2^ceil(u/2)` elements. The data blocks of +/// the first levels are addressed straight from the index block; later +/// levels go through a super block (EASB). Data blocks larger than a page +/// (`2^MAX_DBLK_PAGE_NELMTS_BITS` elements) are paged, with their page-init +/// bits kept in the owning super block. Only blocks holding a defined element +/// are allocated; the rest keep the undefined address, as in a file the +/// library wrote. pub fn build_extensible_array_at( - chunks: &[WrittenChunk], + slots: &[Option], offset_size: u8, length_size: u8, has_filters: bool, ea_base_address: u64, ) -> Vec { let os = offset_size as usize; - let num_elements = chunks.len(); - - // Compute element encoding size (same logic as Fixed Array) - let chunk_size_bytes: usize = if has_filters { - let max_raw = chunks.iter().map(|c| c.raw_size).max().unwrap_or(1); - let log2_val = if max_raw <= 1 { - 0 - } else { - 63 - max_raw.leading_zeros() - }; - let len = 1 + ((log2_val + 8) / 8) as usize; - len.min(8) - } else { - 0 - }; - - let elem_size = if has_filters { - os + chunk_size_bytes + 4 - } else { - os - }; - + let chunk_size_bytes = has_filters.then(|| filtered_chunk_size_len(slots)); + let elem_size = os + chunk_size_bytes.map_or(0, |n| n + 4); let client_id: u8 = if has_filters { 1 } else { 0 }; + let arr_off_size = (MAX_NELMTS_BITS as usize).div_ceil(8); + let page_nelmts = 1usize << MAX_DBLK_PAGE_NELMTS_BITS; + let idx_blk = IDX_BLK_ELMTS as usize; - // EA creation parameters — must match HDF5 C library defaults exactly - let max_nelmts_bits: u8 = 32; - let idx_blk_elmts: u8 = 4; - let min_dblk_nelmts: u8 = 16; - let super_blk_min_nelmts: u8 = 4; - let max_dblk_nelmts_bits: u8 = 10; + // Elements past the last defined one are never realised + // (`max_idx_set` is one past the highest index ever set). + let max_idx_set = slots.iter().rposition(Option::is_some).map_or(0, |i| i + 1); + let slots = &slots[..max_idx_set]; + let defined_in = |start: usize, n: usize| -> bool { + let lo = idx_blk.saturating_add(start).min(slots.len()); + let hi = idx_blk + .saturating_add(start) + .saturating_add(n) + .min(slots.len()); + slots[lo..hi].iter().any(Option::is_some) + }; - // EAHD size: fixed(12) + 6 stats(6*length_size) + addr(offset_size) + checksum(4) + // Super block levels: (ndblks, dblk_nelmts, first element). + let log2_dmin = (DATA_BLK_MIN_ELMTS as u32).trailing_zeros() as usize; + let nsblks = 1 + MAX_NELMTS_BITS as usize - log2_dmin; + let ndblk_addrs = 2 * (SUP_BLK_MIN_DATA_PTRS as usize - 1); + let mut levels: Vec<(usize, usize, usize)> = Vec::with_capacity(nsblks); + let mut start = 0usize; + for u in 0..nsblks { + let ndblks = 1usize << (u / 2); + let nelmts = (DATA_BLK_MIN_ELMTS as usize) << u.div_ceil(2); + levels.push((ndblks, nelmts, start)); + // Saturate: on 32-bit targets the last levels only need to compare + // as "beyond the end". + start = start.saturating_add(ndblks.saturating_mul(nelmts)); + } + // Levels whose data blocks the index block addresses directly. + let mut direct_levels = 0; + let mut n = 0; + while n < ndblk_addrs { + n += levels[direct_levels].0; + direct_levels += 1; + } + let nsblk_addrs = nsblks - direct_levels; + + let dblk_size = |nelmts: usize| -> usize { + let prefix = 4 + 1 + 1 + os + arr_off_size + 4; + if nelmts > page_nelmts { + prefix + (nelmts / page_nelmts) * (page_nelmts * elem_size + 4) + } else { + prefix + nelmts * elem_size + } + }; + let sblk_bitmap_len = |ndblks: usize, nelmts: usize| -> usize { + if nelmts > page_nelmts { + ndblks * (nelmts / page_nelmts).div_ceil(8) + } else { + 0 + } + }; + + // Plan addresses: header, index block, the direct data blocks, then each + // allocated super block followed by its allocated data blocks. let aehd_size = 4 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 6 * length_size as usize + os + 4; let aeib_address = ea_base_address + aehd_size as u64; + let aeib_size = 4 + 1 + 1 + os + idx_blk * elem_size + ndblk_addrs * os + nsblk_addrs * os + 4; + let mut cursor = aeib_address + aeib_size as u64; - // Determine how many elements go inline vs data blocks - let n_inline = (idx_blk_elmts as usize).min(num_elements); - let remaining_after_inline = num_elements.saturating_sub(n_inline); + let mut ndata_blks = 0u64; + let mut data_blk_size = 0u64; + let mut nsuper_blks = 0u64; + let mut super_blk_size = 0u64; + let mut realized = idx_blk as u64; - // Compute super block layout per HDF5 spec - let sblk_min = super_blk_min_nelmts as usize; - let log2_dblk_min = if min_dblk_nelmts <= 1 { - 0 - } else { - (min_dblk_nelmts as u32).trailing_zeros() as usize + let mut plan_dblk = |cursor: &mut u64, start: usize, nelmts: usize| -> DataBlock { + let addr = defined_in(start, nelmts).then(|| { + let a = *cursor; + let size = dblk_size(nelmts) as u64; + *cursor += size; + ndata_blks += 1; + data_blk_size += size; + realized += nelmts as u64; + a + }); + DataBlock { + start, + nelmts, + addr, + } }; - let nsblks = (max_nelmts_bits as usize).saturating_sub(log2_dblk_min) + 1; - // Direct data block addresses (from super blocks 0..sblk_min-1) - let mut dblk_sizes: Vec = Vec::new(); - for sblk_idx in 0..sblk_min.min(nsblks) { - let ndblks = 1usize << (sblk_idx / 2); - let dblk_nelmts = (min_dblk_nelmts as usize) * (1 << sblk_idx.div_ceil(2)); - for _ in 0..ndblks { - dblk_sizes.push(dblk_nelmts); + let mut direct: Vec = Vec::with_capacity(ndblk_addrs); + for &(ndblks, nelmts, first) in &levels[..direct_levels] { + for k in 0..ndblks { + direct.push(plan_dblk(&mut cursor, first + k * nelmts, nelmts)); } } - let n_direct_dblks = dblk_sizes.len(); - - // Super block addresses (for super blocks sblk_min..nsblks-1) - let n_sblk_addrs = nsblks.saturating_sub(sblk_min); - - // EAIB size - let aeib_size = 4 - + 1 - + 1 - + os - + idx_blk_elmts as usize * elem_size - + n_direct_dblks * os - + n_sblk_addrs * os - + 4; - - // Build AEHD - let mut aehd = Vec::with_capacity(aehd_size); - aehd.extend_from_slice(b"EAHD"); - aehd.push(0); // version - aehd.push(client_id); - aehd.push(elem_size as u8); - aehd.push(max_nelmts_bits); - aehd.push(idx_blk_elmts); - aehd.push(min_dblk_nelmts); - aehd.push(super_blk_min_nelmts); - aehd.push(max_dblk_nelmts_bits); - - // Count data blocks that will have chunks - let n_active_dblks: u64 = if remaining_after_inline > 0 { - let mut count = 0u64; - let mut ci = n_inline; - for &sz in &dblk_sizes { - if ci < num_elements { - count += 1; - ci += sz; - } + // (super block address, level, its data blocks) + let mut supers: Vec<(Option, usize, Vec)> = Vec::with_capacity(nsblk_addrs); + for (u, &(ndblks, nelmts, first)) in levels.iter().enumerate().skip(direct_levels) { + if !defined_in(first, ndblks.saturating_mul(nelmts)) { + supers.push((None, u, Vec::new())); + continue; } - count - } else { - 0 - }; - let blk_off_size = (max_nelmts_bits as usize).div_ceil(8); - let aedb_header_overhead = 4 + 1 + 1 + os + blk_off_size + 4; - let data_blk_total_size: u64 = if remaining_after_inline > 0 { - let mut total = 0u64; - let mut ci = n_inline; - for &sz in &dblk_sizes { - if ci < num_elements { - total += (aedb_header_overhead + sz * elem_size) as u64; - ci += sz; - } - } - total - } else { - 0 - }; - let max_idx_set: u64 = if remaining_after_inline > 0 { - let mut max_set = idx_blk_elmts as u64; - let mut ci = n_inline; - for &sz in &dblk_sizes { - if ci < num_elements { - max_set += sz as u64; - ci += sz; - } - } - max_set - } else { - idx_blk_elmts as u64 - }; + let sb_size = + 4 + 1 + 1 + os + arr_off_size + sblk_bitmap_len(ndblks, nelmts) + ndblks * os + 4; + let sb_addr = cursor; + cursor += sb_size as u64; + nsuper_blks += 1; + super_blk_size += sb_size as u64; + let dblks = (0..ndblks) + .map(|k| plan_dblk(&mut cursor, first + k * nelmts, nelmts)) + .collect(); + supers.push((Some(sb_addr), u, dblks)); + } + let slot = |i: usize| slots.get(i).and_then(Option::as_ref); let write_length = |buf: &mut Vec, val: u64| match length_size { 4 => buf.extend_from_slice(&(val as u32).to_le_bytes()), _ => buf.extend_from_slice(&val.to_le_bytes()), }; - let write_addr = |buf: &mut Vec, val: u64| match offset_size { - 4 => buf.extend_from_slice(&(val as u32).to_le_bytes()), - _ => buf.extend_from_slice(&val.to_le_bytes()), + let write_addr_opt = |buf: &mut Vec, addr: Option| match addr { + Some(a) => push_addr(buf, a, offset_size), + None => buf.extend(core::iter::repeat_n(0xFF, os)), + }; + let block_prefix = |buf: &mut Vec, sig: &[u8; 4], block_off: usize| { + buf.extend_from_slice(sig); + buf.push(0); // version + buf.push(client_id); + push_addr(buf, ea_base_address, offset_size); + buf.extend_from_slice(&(block_off as u64).to_le_bytes()[..arr_off_size]); + }; + // Serialise one data block (paged or not) onto `out`. + let write_dblk = |out: &mut Vec, db: &DataBlock| { + let at = out.len(); + block_prefix(out, b"EADB", db.start); + let first = idx_blk + db.start; + if db.nelmts > page_nelmts { + // Paged: the prefix carries only its own checksum; each page + // follows with one of its own. + let sum = jenkins_lookup3(&out[at..]); + out.extend_from_slice(&sum.to_le_bytes()); + for p in 0..db.nelmts / page_nelmts { + let page_at = out.len(); + for e in 0..page_nelmts { + let i = first + p * page_nelmts + e; + push_index_element(out, slot(i), offset_size, chunk_size_bytes); + } + let sum = jenkins_lookup3(&out[page_at..]); + out.extend_from_slice(&sum.to_le_bytes()); + } + } else { + for i in first..first + db.nelmts { + push_index_element(out, slot(i), offset_size, chunk_size_bytes); + } + let sum = jenkins_lookup3(&out[at..]); + out.extend_from_slice(&sum.to_le_bytes()); + } + debug_assert_eq!(out.len() - at, dblk_size(db.nelmts)); }; - write_length(&mut aehd, 0); - write_length(&mut aehd, 0); - write_length(&mut aehd, n_active_dblks); - write_length(&mut aehd, data_blk_total_size); - write_length(&mut aehd, num_elements as u64); - write_length(&mut aehd, max_idx_set); + // Header (EAHD). The six statistics are, in order: super blocks, their + // bytes, data blocks, their bytes, max index set, elements realised. + let mut out = Vec::with_capacity((cursor - ea_base_address) as usize); + out.extend_from_slice(b"EAHD"); + out.push(0); // version + out.push(client_id); + out.push(elem_size as u8); + out.push(MAX_NELMTS_BITS); + out.push(IDX_BLK_ELMTS); + out.push(DATA_BLK_MIN_ELMTS); + out.push(SUP_BLK_MIN_DATA_PTRS); + out.push(MAX_DBLK_PAGE_NELMTS_BITS); + write_length(&mut out, nsuper_blks); + write_length(&mut out, super_blk_size); + write_length(&mut out, ndata_blks); + write_length(&mut out, data_blk_size); + write_length(&mut out, max_idx_set as u64); + write_length(&mut out, realized); + push_addr(&mut out, aeib_address, offset_size); + let sum = jenkins_lookup3(&out); + out.extend_from_slice(&sum.to_le_bytes()); + debug_assert_eq!(out.len(), aehd_size); - write_addr(&mut aehd, aeib_address); - - let aehd_checksum = jenkins_lookup3(&aehd); - aehd.extend_from_slice(&aehd_checksum.to_le_bytes()); - debug_assert_eq!(aehd.len(), aehd_size); - - // Build AEIB - let mut aeib = Vec::with_capacity(aeib_size); - aeib.extend_from_slice(b"EAIB"); - aeib.push(0); - aeib.push(client_id); - - match offset_size { - 4 => aeib.extend_from_slice(&(ea_base_address as u32).to_le_bytes()), - 8 => aeib.extend_from_slice(&ea_base_address.to_le_bytes()), - _ => aeib.extend_from_slice(&ea_base_address.to_le_bytes()), + // Index block (EAIB): inline elements, data block and super block + // addresses. + let ib_start = out.len(); + out.extend_from_slice(b"EAIB"); + out.push(0); + out.push(client_id); + push_addr(&mut out, ea_base_address, offset_size); + for i in 0..idx_blk { + push_index_element(&mut out, slot(i), offset_size, chunk_size_bytes); } - - // Inline elements - #[allow(clippy::needless_range_loop)] - for i in 0..idx_blk_elmts as usize { - if i < n_inline { - write_chunk_element( - &mut aeib, - &chunks[i], - offset_size, - has_filters, - chunk_size_bytes, - ); - } else { - write_undefined_element(&mut aeib, offset_size, has_filters, chunk_size_bytes); - } + for db in &direct { + write_addr_opt(&mut out, db.addr); } + for (sb_addr, _, _) in &supers { + write_addr_opt(&mut out, *sb_addr); + } + let sum = jenkins_lookup3(&out[ib_start..]); + out.extend_from_slice(&sum.to_le_bytes()); + debug_assert_eq!(out.len() - ib_start, aeib_size); - // Data block addresses + build data blocks - let mut data_blocks_buf = Vec::new(); - let dblks_base = aeib_address + aeib_size as u64; - let mut dblk_cursor = dblks_base; - let mut chunk_idx = n_inline; - - for &nelmts in &dblk_sizes { - if chunk_idx >= num_elements { - match offset_size { - 4 => aeib.extend_from_slice(&u32::MAX.to_le_bytes()), - 8 => aeib.extend_from_slice(&u64::MAX.to_le_bytes()), - _ => aeib.extend_from_slice(&u64::MAX.to_le_bytes()), - } + for db in direct.iter().filter(|d| d.addr.is_some()) { + write_dblk(&mut out, db); + } + for (sb_addr, u, dblks) in &supers { + if sb_addr.is_none() { continue; } - - match offset_size { - 4 => aeib.extend_from_slice(&(dblk_cursor as u32).to_le_bytes()), - 8 => aeib.extend_from_slice(&dblk_cursor.to_le_bytes()), - _ => aeib.extend_from_slice(&dblk_cursor.to_le_bytes()), - } - - // Build EADB - let mut aedb = Vec::new(); - aedb.extend_from_slice(b"EADB"); - aedb.push(0); - aedb.push(client_id); - match offset_size { - 4 => aedb.extend_from_slice(&(ea_base_address as u32).to_le_bytes()), - 8 => aedb.extend_from_slice(&ea_base_address.to_le_bytes()), - _ => aedb.extend_from_slice(&ea_base_address.to_le_bytes()), - } - - let blk_off_size = (max_nelmts_bits as usize).div_ceil(8); - let blk_off_val = (chunk_idx - n_inline) as u64; - aedb.extend_from_slice(&blk_off_val.to_le_bytes()[..blk_off_size]); - - for slot in 0..nelmts { - if chunk_idx + slot < num_elements { - write_chunk_element( - &mut aedb, - &chunks[chunk_idx + slot], - offset_size, - has_filters, - chunk_size_bytes, - ); - } else { - write_undefined_element(&mut aedb, offset_size, has_filters, chunk_size_bytes); + let (ndblks, nelmts, first) = levels[*u]; + let sb_start = out.len(); + block_prefix(&mut out, b"EASB", first); + if nelmts > page_nelmts { + // Page-init bits, `npages` per data block, packed MSB-first + // (`H5VM_bit_set`): every page of an allocated data block is + // written. + let npages = nelmts / page_nelmts; + let mut bitmap = vec![0u8; sblk_bitmap_len(ndblks, nelmts)]; + for (k, db) in dblks.iter().enumerate() { + if db.addr.is_some() { + for p in 0..npages { + let bit = k * npages + p; + bitmap[bit / 8] |= 0x80 >> (bit % 8); + } + } } + out.extend_from_slice(&bitmap); } - - let aedb_checksum = jenkins_lookup3(&aedb); - aedb.extend_from_slice(&aedb_checksum.to_le_bytes()); - - dblk_cursor += aedb.len() as u64; - data_blocks_buf.extend_from_slice(&aedb); - chunk_idx += nelmts; - } - - // Super block addresses (all undefined) - for _ in 0..n_sblk_addrs { - match offset_size { - 4 => aeib.extend_from_slice(&u32::MAX.to_le_bytes()), - 8 => aeib.extend_from_slice(&u64::MAX.to_le_bytes()), - _ => aeib.extend_from_slice(&u64::MAX.to_le_bytes()), + for db in dblks { + write_addr_opt(&mut out, db.addr); + } + let sum = jenkins_lookup3(&out[sb_start..]); + out.extend_from_slice(&sum.to_le_bytes()); + for db in dblks.iter().filter(|d| d.addr.is_some()) { + write_dblk(&mut out, db); } } - - let aeib_checksum = jenkins_lookup3(&aeib); - aeib.extend_from_slice(&aeib_checksum.to_le_bytes()); - debug_assert_eq!(aeib.len(), aeib_size); - - let mut combined = aehd; - combined.extend_from_slice(&aeib); - combined.extend_from_slice(&data_blocks_buf); - combined -} - -fn write_chunk_element( - buf: &mut Vec, - chunk: &WrittenChunk, - offset_size: u8, - has_filters: bool, - chunk_size_bytes: usize, -) { - match offset_size { - 4 => buf.extend_from_slice(&(chunk.address as u32).to_le_bytes()), - 8 => buf.extend_from_slice(&chunk.address.to_le_bytes()), - _ => buf.extend_from_slice(&chunk.address.to_le_bytes()), - } - if has_filters { - let cs_bytes = chunk.compressed_size.to_le_bytes(); - buf.extend_from_slice(&cs_bytes[..chunk_size_bytes]); - buf.extend_from_slice(&chunk.filter_mask.to_le_bytes()); - } -} - -fn write_undefined_element( - buf: &mut Vec, - offset_size: u8, - has_filters: bool, - chunk_size_bytes: usize, -) { - let os = offset_size as usize; - // Use extend with repeat to avoid heap-allocating a temporary Vec on each call. - buf.extend(core::iter::repeat_n(0xFF, os)); - if has_filters { - buf.extend(core::iter::repeat_n(0x00, chunk_size_bytes)); - buf.extend_from_slice(&0u32.to_le_bytes()); - } + debug_assert_eq!(out.len() as u64, cursor - ea_base_address); + out } diff --git a/crates/clawhdf5-format/src/extensible_array.rs b/crates/clawhdf5-format/src/extensible_array.rs index 4416a76..ba1c822 100644 --- a/crates/clawhdf5-format/src/extensible_array.rs +++ b/crates/clawhdf5-format/src/extensible_array.rs @@ -9,6 +9,7 @@ extern crate alloc; #[cfg(not(feature = "std"))] use alloc::{format, vec, vec::Vec}; +use crate::chunk_grid::ChunkGrid; use crate::chunked_read::ChunkInfo; use crate::error::FormatError; @@ -203,8 +204,7 @@ fn read_element( offset_size: u8, chunk_byte_size: u64, linear_index: usize, - num_chunks_per_dim: &[u64], - chunk_dimensions: &[u32], + grid: &ChunkGrid, ) -> Result<(Option, usize), FormatError> { let os = offset_size as usize; @@ -220,7 +220,10 @@ fn read_element( return Ok((None, os)); } let address = read_offset(data, pos, offset_size)?; - let offsets = index_to_chunk_offsets(linear_index, num_chunks_per_dim, chunk_dimensions); + // A slot beyond the current extent is ignored, as the library does. + let Some(offsets) = grid.offsets(linear_index as u64) else { + return Ok((None, os)); + }; Ok(( Some(ChunkInfo { chunk_size: chunk_byte_size as u32, @@ -261,7 +264,9 @@ fn read_element( data[fm_off + 2], data[fm_off + 3], ]); - let offsets = index_to_chunk_offsets(linear_index, num_chunks_per_dim, chunk_dimensions); + let Some(offsets) = grid.offsets(linear_index as u64) else { + return Ok((None, elem_total)); + }; Ok(( Some(ChunkInfo { chunk_size: chunk_size as u32, @@ -274,27 +279,6 @@ fn read_element( } } -/// Convert a linear chunk index to N-dimensional chunk offsets in dataset space. -fn index_to_chunk_offsets( - index: usize, - num_chunks_per_dim: &[u64], - chunk_dimensions: &[u32], -) -> Vec { - let rank = num_chunks_per_dim.len(); - let mut offsets = vec![0u64; rank]; - let mut remaining = index as u64; - for d in (0..rank).rev() { - let nchunks = num_chunks_per_dim[d]; - if nchunks == 0 { - continue; - } - let chunk_idx = remaining % nchunks; - remaining /= nchunks; - offsets[d] = chunk_idx * chunk_dimensions[d] as u64; - } - offsets -} - /// Collect elements from a data block at the given offset. #[allow(clippy::too_many_arguments)] /// Layout of super block `u`, per the HDF5 spec: the number of data blocks it @@ -339,8 +323,7 @@ fn read_data_block_elements( offset_size: u8, chunk_byte_size: u64, start_index: usize, - num_chunks_per_dim: &[u64], - chunk_dimensions: &[u32], + grid: &ChunkGrid, page_init: &[u8], first_page: usize, ) -> Result, FormatError> { @@ -376,8 +359,7 @@ fn read_data_block_elements( offset_size, chunk_byte_size, first_index + i, - num_chunks_per_dim, - chunk_dimensions, + grid, )?; if let Some(ci) = info { chunks.push(ci); @@ -449,25 +431,19 @@ pub fn read_extensible_array_chunks( file_data: &[u8], header: &ExtensibleArrayHeader, dataset_dims: &[u64], + max_dims: Option<&[u64]>, chunk_dimensions: &[u32], element_size: u32, offset_size: u8, _length_size: u8, ) -> Result, FormatError> { - let rank = chunk_dimensions.len(); let os = offset_size as usize; - let mut num_chunks_per_dim = Vec::with_capacity(rank); - for d in 0..rank { - let ch_dim = chunk_dimensions[d] as u64; - if ch_dim == 0 { - return Err(FormatError::ChunkedReadError( - "chunk dimension is zero".into(), - )); - } - let ds_dim = dataset_dims[d]; - num_chunks_per_dim.push(ds_dim.div_ceil(ch_dim)); - } + // Linear indexes follow the maximum dimensions, with the unlimited + // dimension swizzled to the slowest position (see `chunk_grid`). + let dims_u64: Vec = chunk_dimensions.iter().map(|&d| d as u64).collect(); + let grid = ChunkGrid::extensible_array(dataset_dims, max_dims, &dims_u64)?; + let grid = &grid; let chunk_byte_size: u64 = chunk_dimensions.iter().map(|&d| d as u64).product::() * element_size as u64; @@ -557,8 +533,7 @@ pub fn read_extensible_array_chunks( offset_size, chunk_byte_size, i, - &num_chunks_per_dim, - chunk_dimensions, + grid, )?; if let Some(ci) = info { chunks.push(ci); @@ -594,8 +569,7 @@ pub fn read_extensible_array_chunks( offset_size, chunk_byte_size, global_index, - &num_chunks_per_dim, - chunk_dimensions, + grid, &[], 0, )?); @@ -625,8 +599,7 @@ pub fn read_extensible_array_chunks( offset_size, chunk_byte_size, global_index, - &num_chunks_per_dim, - chunk_dimensions, + grid, )?); } global_index = @@ -653,8 +626,7 @@ fn read_super_block( offset_size: u8, chunk_byte_size: u64, start_index: usize, - num_chunks_per_dim: &[u64], - chunk_dimensions: &[u32], + grid: &ChunkGrid, ) -> Result, FormatError> { let os = offset_size as usize; let sb_header_size = 4 + 1 + 1 + os + arr_off_size(header); @@ -710,8 +682,7 @@ fn read_super_block( offset_size, chunk_byte_size, global_idx, - num_chunks_per_dim, - chunk_dimensions, + grid, bitmap, i * npages, )?); @@ -735,35 +706,18 @@ mod tests { } #[test] fn index_to_offsets_1d() { - let num_chunks = vec![5u64]; - let chunk_dims = vec![20u32]; - assert_eq!(index_to_chunk_offsets(0, &num_chunks, &chunk_dims), vec![0]); - assert_eq!( - index_to_chunk_offsets(1, &num_chunks, &chunk_dims), - vec![20] - ); - assert_eq!( - index_to_chunk_offsets(4, &num_chunks, &chunk_dims), - vec![80] - ); + let g = ChunkGrid::fixed_array(&[100], None, &[20]).unwrap(); + assert_eq!(g.offsets(0).unwrap(), vec![0]); + assert_eq!(g.offsets(1).unwrap(), vec![20]); + assert_eq!(g.offsets(4).unwrap(), vec![80]); } #[test] fn index_to_offsets_2d() { - let num_chunks = vec![3u64, 2]; - let chunk_dims = vec![4u32, 3]; - assert_eq!( - index_to_chunk_offsets(0, &num_chunks, &chunk_dims), - vec![0, 0] - ); - assert_eq!( - index_to_chunk_offsets(1, &num_chunks, &chunk_dims), - vec![0, 3] - ); - assert_eq!( - index_to_chunk_offsets(2, &num_chunks, &chunk_dims), - vec![4, 0] - ); + let g = ChunkGrid::fixed_array(&[10, 6], None, &[4, 3]).unwrap(); + assert_eq!(g.offsets(0).unwrap(), vec![0, 0]); + assert_eq!(g.offsets(1).unwrap(), vec![0, 3]); + assert_eq!(g.offsets(2).unwrap(), vec![4, 0]); } #[test] @@ -830,7 +784,7 @@ mod tests { index_block_address: (usize::MAX - 4) as u64, }; let buf = vec![0u8; 64]; - let r = read_extensible_array_chunks(&buf, &header, &[100], &[20], 8, 8, 8); + let r = read_extensible_array_chunks(&buf, &header, &[100], None, &[20], 8, 8, 8); assert!(r.is_err()); } @@ -913,9 +867,17 @@ mod tests { let header = ExtensibleArrayHeader::parse(&file_data, aehd_offset, os, ls).unwrap(); let ds_dims = vec![40u64]; // 2 chunks × 20 elements let chunk_dims = vec![20u32]; - let chunks = - read_extensible_array_chunks(&file_data, &header, &ds_dims, &chunk_dims, 8, os, ls) - .unwrap(); + let chunks = read_extensible_array_chunks( + &file_data, + &header, + &ds_dims, + None, + &chunk_dims, + 8, + os, + ls, + ) + .unwrap(); assert_eq!(chunks.len(), 2); assert_eq!(chunks[0].address, base_addr); @@ -1023,9 +985,17 @@ mod tests { let header = ExtensibleArrayHeader::parse(&file_data, aehd_offset, os, ls).unwrap(); let ds_dims = vec![40u64]; let chunk_dims = vec![10u32]; - let chunks = - read_extensible_array_chunks(&file_data, &header, &ds_dims, &chunk_dims, 8, os, ls) - .unwrap(); + let chunks = read_extensible_array_chunks( + &file_data, + &header, + &ds_dims, + None, + &chunk_dims, + 8, + os, + ls, + ) + .unwrap(); assert_eq!(chunks.len(), 4); for (i, c) in chunks.iter().enumerate() { @@ -1047,10 +1017,8 @@ mod tests { #[test] fn read_element_unallocated() { let data = vec![0xFFu8; 16]; - let num_chunks = vec![5u64]; - let chunk_dims = vec![10u32]; - let (info, consumed) = - read_element(&data, 0, 0, 8, 8, 80, 0, &num_chunks, &chunk_dims).unwrap(); + let grid = ChunkGrid::fixed_array(&[50], None, &[10]).unwrap(); + let (info, consumed) = read_element(&data, 0, 0, 8, 8, 80, 0, &grid).unwrap(); assert!(info.is_none()); assert_eq!(consumed, 8); } @@ -1069,20 +1037,9 @@ mod tests { // Filter mask data[12..16].copy_from_slice(&0u32.to_le_bytes()); - let num_chunks = vec![5u64]; - let chunk_dims = vec![10u32]; - let (info, consumed) = read_element( - &data, - 0, - 1, - elem_size as u8, - os, - 80, - 2, - &num_chunks, - &chunk_dims, - ) - .unwrap(); + let grid = ChunkGrid::fixed_array(&[50], None, &[10]).unwrap(); + let (info, consumed) = + read_element(&data, 0, 1, elem_size as u8, os, 80, 2, &grid).unwrap(); let ci = info.unwrap(); assert_eq!(ci.address, 0x2000); assert_eq!(ci.chunk_size, 120); diff --git a/crates/clawhdf5-format/src/file_writer.rs b/crates/clawhdf5-format/src/file_writer.rs index c262cd4..644b64e 100644 --- a/crates/clawhdf5-format/src/file_writer.rs +++ b/crates/clawhdf5-format/src/file_writer.rs @@ -1124,7 +1124,12 @@ impl FileWriter { let is_chunked: Vec = all_ds .iter() .enumerate() - .map(|(i, d)| !is_vds[i] && (d.chunk_options.is_chunked() || d.maxshape.is_some())) + .map(|(i, d)| { + // Only a dataset that can grow needs chunks; a maxshape equal + // to the shape is as fixed as no maxshape at all. + let resizable = d.maxshape.as_ref().is_some_and(|m| *m != d.ds.dimensions); + !is_vds[i] && (d.chunk_options.is_chunked() || resizable) + }) .collect(); // Determine which datasets use compact storage let is_compact: Vec = all_ds @@ -1244,7 +1249,7 @@ impl FileWriter { &pre, dummy_cursor, d.maxshape.as_deref(), - ); + )?; dummy_cursor += result.data_bytes.len() as u64; let dense_blob = if ds_dense[i] { Some(build_dense_attrs(&d.attrs, 0)) @@ -1424,7 +1429,7 @@ impl FileWriter { .expect("chunked dataset missing precompressed cache"), base_address, d.maxshape.as_deref(), - ); + )?; cursor2 += result.data_bytes.len(); let oh = build_chunked_dataset_oh( &d.dt, diff --git a/crates/clawhdf5-format/src/fixed_array.rs b/crates/clawhdf5-format/src/fixed_array.rs index 4c79f03..8546224 100644 --- a/crates/clawhdf5-format/src/fixed_array.rs +++ b/crates/clawhdf5-format/src/fixed_array.rs @@ -6,6 +6,7 @@ extern crate alloc; #[cfg(not(feature = "std"))] use alloc::{format, vec, vec::Vec}; +use crate::chunk_grid::ChunkGrid; use crate::chunked_read::ChunkInfo; use crate::error::FormatError; @@ -151,13 +152,13 @@ pub fn read_fixed_array_chunks( file_data: &[u8], header: &FixedArrayHeader, dataset_dims: &[u64], + max_dims: Option<&[u64]>, chunk_dimensions: &[u32], element_size: u32, offset_size: u8, _length_size: u8, ) -> Result, FormatError> { let db_offset = header.data_block_address as usize; - let rank = chunk_dimensions.len(); // Parse data block header: FADB(4) + version(1) + client_id(1) + header_address(offset_size) let db_header_size = 4 + 1 + 1 + offset_size as usize; @@ -198,19 +199,10 @@ pub fn read_fixed_array_chunks( )) }; - // Compute chunk offsets based on index. - // Chunks are stored in row-major order within the dataset space. - let mut num_chunks_per_dim = Vec::with_capacity(rank); - for d_idx in 0..rank { - let ch_dim = chunk_dimensions[d_idx] as u64; - if ch_dim == 0 { - return Err(FormatError::ChunkedReadError( - "chunk dimension is zero".into(), - )); - } - let ds_dim = dataset_dims[d_idx]; - num_chunks_per_dim.push(ds_dim.div_ceil(ch_dim)); - } + // The index is laid out over the chunk grid of the *maximum* dimensions + // (row-major), so a dataset smaller than its maxshape has gaps. + let dims_u64: Vec = chunk_dimensions.iter().map(|&d| d as u64).collect(); + let grid = ChunkGrid::fixed_array(dataset_dims, max_dims, &dims_u64)?; let chunk_byte_size: u64 = chunk_dimensions.iter().map(|&d| d as u64).product::() * element_size as u64; @@ -226,7 +218,11 @@ pub fn read_fixed_array_chunks( header.element_size, chunk_byte_size, )? { - let offsets = index_to_chunk_offsets(i, &num_chunks_per_dim, chunk_dimensions); + // A slot beyond the current extent is ignored, as the + // library does. + let Some(offsets) = grid.offsets(i as u64) else { + return Ok(()); + }; chunks.push(ChunkInfo { chunk_size, filter_mask, @@ -367,27 +363,6 @@ fn parse_fa_element( } } -/// Convert a linear chunk index to N-dimensional chunk offsets in dataset space. -fn index_to_chunk_offsets( - index: usize, - num_chunks_per_dim: &[u64], - chunk_dimensions: &[u32], -) -> Vec { - let rank = num_chunks_per_dim.len(); - let mut offsets = vec![0u64; rank]; - let mut remaining = index as u64; - for d in (0..rank).rev() { - let nchunks = num_chunks_per_dim[d]; - if nchunks == 0 { - continue; - } - let chunk_idx = remaining % nchunks; - remaining /= nchunks; - offsets[d] = chunk_idx * chunk_dimensions[d] as u64; - } - offsets -} - /// Read a variable-length little-endian unsigned integer. fn read_variable_length(data: &[u8], size: usize) -> Result { if size > 8 || data.len() < size { @@ -416,44 +391,21 @@ mod tests { #[test] fn index_to_offsets_1d() { - let num_chunks = vec![5u64]; - let chunk_dims = vec![20u32]; - assert_eq!(index_to_chunk_offsets(0, &num_chunks, &chunk_dims), vec![0]); - assert_eq!( - index_to_chunk_offsets(1, &num_chunks, &chunk_dims), - vec![20] - ); - assert_eq!( - index_to_chunk_offsets(4, &num_chunks, &chunk_dims), - vec![80] - ); + let g = ChunkGrid::fixed_array(&[100], None, &[20]).unwrap(); + assert_eq!(g.offsets(0).unwrap(), vec![0]); + assert_eq!(g.offsets(1).unwrap(), vec![20]); + assert_eq!(g.offsets(4).unwrap(), vec![80]); } #[test] fn index_to_offsets_2d() { // 10x6 dataset with 4x3 chunks => ceil(10/4)=3, ceil(6/3)=2 => 6 chunks - let num_chunks = vec![3u64, 2]; - let chunk_dims = vec![4u32, 3]; - assert_eq!( - index_to_chunk_offsets(0, &num_chunks, &chunk_dims), - vec![0, 0] - ); - assert_eq!( - index_to_chunk_offsets(1, &num_chunks, &chunk_dims), - vec![0, 3] - ); - assert_eq!( - index_to_chunk_offsets(2, &num_chunks, &chunk_dims), - vec![4, 0] - ); - assert_eq!( - index_to_chunk_offsets(3, &num_chunks, &chunk_dims), - vec![4, 3] - ); - assert_eq!( - index_to_chunk_offsets(5, &num_chunks, &chunk_dims), - vec![8, 3] - ); + let g = ChunkGrid::fixed_array(&[10, 6], None, &[4, 3]).unwrap(); + assert_eq!(g.offsets(0).unwrap(), vec![0, 0]); + assert_eq!(g.offsets(1).unwrap(), vec![0, 3]); + assert_eq!(g.offsets(2).unwrap(), vec![4, 0]); + assert_eq!(g.offsets(3).unwrap(), vec![4, 3]); + assert_eq!(g.offsets(5).unwrap(), vec![8, 3]); } #[test] @@ -517,7 +469,7 @@ mod tests { let read = |f: &[u8], fahd: usize| -> Result, FormatError> { let h = FixedArrayHeader::parse(f, fahd, 8, 8)?; - read_fixed_array_chunks(f, &h, &[60], &[20], 8, 8, 8) + read_fixed_array_chunks(f, &h, &[60], None, &[20], 8, 8, 8) }; let (clean, fahd) = build(); @@ -562,7 +514,7 @@ mod tests { let db = 0x100usize; buf[db..db + 4].copy_from_slice(b"FADB"); let header = FixedArrayHeader::parse(&buf, fahd, 8, 8).unwrap(); - let r = read_fixed_array_chunks(&buf, &header, &[100], &[20], 8, 8, 8); + let r = read_fixed_array_chunks(&buf, &header, &[100], None, &[20], 8, 8, 8); assert!(r.is_err()); } @@ -579,7 +531,7 @@ mod tests { stamp_checksum(&mut buf, fahd, fahd + 24); buf[0x80..0x84].copy_from_slice(b"FADB"); let header = FixedArrayHeader::parse(&buf, fahd, 8, 8).unwrap(); - let r = read_fixed_array_chunks(&buf, &header, &[100], &[20], 8, 8, 8); + let r = read_fixed_array_chunks(&buf, &header, &[100], None, &[20], 8, 8, 8); assert!(r.is_err()); } @@ -602,7 +554,7 @@ mod tests { data_block_address: (usize::MAX - 4) as u64, }; let buf = vec![0u8; 64]; - let r = read_fixed_array_chunks(&buf, &header, &[100], &[20], 8, 8, 8); + let r = read_fixed_array_chunks(&buf, &header, &[100], None, &[20], 8, 8, 8); assert!(r.is_err()); } @@ -664,6 +616,7 @@ mod tests { &file_data, &header, &ds_dims, + None, &chunk_dims, 8, offset_size, @@ -740,6 +693,7 @@ mod tests { &file_data, &header, &ds_dims, + None, &chunk_dims, 8, offset_size, @@ -840,6 +794,7 @@ mod tests { &file_data, &header, &ds_dims, + None, &chunk_dims, 8, offset_size, diff --git a/crates/clawhdf5-format/src/lib.rs b/crates/clawhdf5-format/src/lib.rs index 54bd326..4a6c2a1 100644 --- a/crates/clawhdf5-format/src/lib.rs +++ b/crates/clawhdf5-format/src/lib.rs @@ -54,6 +54,7 @@ pub mod btree_v1; pub mod btree_v2; pub mod checksum; pub mod chunk_cache; +mod chunk_grid; pub mod chunk_index; pub mod chunked_read; pub mod chunked_write; diff --git a/crates/clawhdf5/tests/chunk_index_interop.rs b/crates/clawhdf5/tests/chunk_index_interop.rs new file mode 100644 index 0000000..b72b7e2 --- /dev/null +++ b/crates/clawhdf5/tests/chunk_index_interop.rs @@ -0,0 +1,593 @@ +//! Fixed Array / Extensible Array chunk-index interop with libhdf5 (via h5py). +//! +//! Both indexes place each chunk at a linear index computed from the +//! dataset's *maximum* dimensions, and the Extensible Array additionally +//! moves its unlimited dimension to the slowest-varying position. Getting +//! either wrong reads (or writes) every chunk after the first row in the +//! wrong place, silently, so these tests compare every value. +//! +//! Skipped when python3 with h5py is unavailable, unless +//! `CLAWHDF5_REQUIRE_INTEROP=1`. + +use std::process::Command; + +use clawhdf5::{File, FileBuilder}; + +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"]) + .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) -> String { + let output = Command::new(python()) + .args(["-c", script]) + .output() + .expect("failed to run python"); + if !output.status.success() { + panic!( + "Python script failed:\nSTDOUT: {}\nSTDERR: {}", + String::from_utf8_lossy(&output.stdout), + String::from_utf8_lossy(&output.stderr) + ); + } + String::from_utf8_lossy(&output.stdout).trim().to_string() +} + +/// Row-major `arange` of `shape`, cropped to `crop` (the current extent). +fn arange_cropped(full: &[usize], crop: &[usize]) -> Vec { + let n: usize = crop.iter().product(); + let mut out = Vec::with_capacity(n); + for flat in 0..n { + let mut rem = flat; + let mut src = 0usize; + let mut stride = 1usize; + let mut coords = vec![0usize; crop.len()]; + for d in (0..crop.len()).rev() { + coords[d] = rem % crop[d]; + rem /= crop[d]; + } + for d in (0..full.len()).rev() { + src += coords[d] * stride; + stride *= full[d]; + } + out.push(src as i32); + } + out +} + +/// One `i4` dataset, filled with `arange` over `full` and then resized to +/// `shape` (equal to `full` unless the case shrinks it). +struct Case { + name: &'static str, + full: Vec, + shape: Vec, + chunks: Vec, + maxshape: &'static str, + extra: &'static str, + index: &'static str, +} + +fn py_tuple(v: &[usize]) -> String { + let parts: Vec = v.iter().map(|x| x.to_string()).collect(); + format!("({},)", parts.join(",")) +} + +/// Have h5py (`libver="latest"`, so Fixed/Extensible Array indexes) write +/// every case to one file, then read each back and compare every value. +fn check_h5py_written(cases: &[Case]) { + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("h5py_chunk_index.h5"); + let path_str = path.display().to_string(); + + let mut script = + format!("import h5py, numpy as np\nf = h5py.File(r'{path_str}', 'w', libver='latest')\n"); + for c in cases { + script += &format!( + "d = f.create_dataset('{name}', data=np.arange({n}, dtype='i4').reshape({full}), \ + chunks={chunks}, maxshape={maxshape}{extra})\n\ + d.resize({shape})\n", + name = c.name, + n = c.full.iter().product::(), + full = py_tuple(&c.full), + chunks = py_tuple(&c.chunks), + maxshape = c.maxshape, + extra = c.extra, + shape = py_tuple(&c.shape), + ); + } + script += "f.close()\n"; + run_python(&script); + + let file = File::open(&path).unwrap(); + for c in cases { + let ds = file.dataset(c.name).unwrap(); + let shape: Vec = ds.shape().unwrap().iter().map(|&d| d as usize).collect(); + assert_eq!(shape, c.shape, "{}: shape", c.name); + let got = ds.read_i32().unwrap(); + let want = arange_cropped(&c.full, &c.shape); + let bad = got.iter().zip(&want).filter(|(a, b)| a != b).count(); + assert_eq!( + got, + want, + "{}: {bad} of {} values differ (index {})", + c.name, + want.len(), + c.index + ); + } +} + +/// h5py-written Extensible Array whose unlimited dimension is not the first, +/// with the current shape smaller than the finite maximum: the library +/// swizzles the unlimited dimension to the slowest position and strides the +/// rest by their maximum chunk counts. +#[test] +fn h5py_extensible_array_partial_extent_reads_correctly() { + skip_if_no_python!(); + check_h5py_written(&[ + // The `ea_fa_partial.h5` repro from the conformance sweep. + Case { + name: "ea_10_none", + full: vec![4, 6], + shape: vec![4, 6], + chunks: vec![2, 3], + maxshape: "(10, None)", + extra: "", + index: "EA, unlimited dim 1", + }, + Case { + name: "ea_none_10", + full: vec![4, 6], + shape: vec![4, 6], + chunks: vec![2, 3], + maxshape: "(None, 10)", + extra: "", + index: "EA, unlimited dim 0", + }, + Case { + name: "ea_3d_mid", + full: vec![3, 4, 5], + shape: vec![3, 4, 5], + chunks: vec![2, 3, 2], + maxshape: "(5, None, 7)", + extra: "", + index: "EA, unlimited dim 1 of 3", + }, + Case { + name: "ea_3d_last_gzip", + full: vec![3, 4, 5], + shape: vec![3, 4, 5], + chunks: vec![2, 3, 2], + maxshape: "(5, 9, None)", + extra: ", compression='gzip'", + index: "EA, unlimited dim 2 of 3, filtered", + }, + // Many chunks: crosses data blocks, super blocks and paging. + Case { + name: "ea_many", + full: vec![3, 1500], + shape: vec![3, 1500], + chunks: vec![1, 1], + maxshape: "(4, None)", + extra: "", + index: "EA, 4500 slots", + }, + // Shrunk after writing: chunks beyond the extent must be ignored. + Case { + name: "ea_shrunk", + full: vec![8, 9], + shape: vec![3, 4], + chunks: vec![2, 3], + maxshape: "(10, None)", + extra: "", + index: "EA, shrunk", + }, + ]); +} + +/// h5py-written Fixed Array with the current shape smaller than a finite +/// maxshape: the index has one slot per chunk of the *maximum* extent. +#[test] +fn h5py_fixed_array_partial_extent_reads_correctly() { + skip_if_no_python!(); + check_h5py_written(&[ + Case { + name: "fa_20_10", + full: vec![4, 6], + shape: vec![4, 6], + chunks: vec![2, 3], + maxshape: "(20, 10)", + extra: "", + index: "FA", + }, + Case { + name: "fa_3d_gzip", + full: vec![3, 4, 5], + shape: vec![3, 4, 5], + chunks: vec![2, 3, 2], + maxshape: "(6, 8, 10)", + extra: ", compression='gzip'", + index: "FA, filtered", + }, + // Paged (> 1024 slots) with most of them beyond the extent. + Case { + name: "fa_paged", + full: vec![30, 50], + shape: vec![30, 50], + chunks: vec![1, 1], + maxshape: "(40, 60)", + extra: "", + index: "FA, 2400 slots, paged", + }, + Case { + name: "fa_shrunk", + full: vec![8, 9], + shape: vec![5, 2], + chunks: vec![2, 3], + maxshape: "(20, 10)", + extra: "", + index: "FA, shrunk", + }, + ]); +} + +// =========================================================================== +// Files we write, read back by libhdf5 (h5py and h5dump) and by us +// =========================================================================== + +/// One `i4` dataset we write, filled with `arange` over `shape`. +struct WriteCase { + name: String, + shape: Vec, + chunks: Vec, + maxshape: Option>, + deflate: bool, +} + +fn wcase(name: &str, shape: &[u64], chunks: &[u64], maxshape: Option<&[u64]>) -> WriteCase { + WriteCase { + name: name.to_string(), + shape: shape.to_vec(), + chunks: chunks.to_vec(), + maxshape: maxshape.map(<[u64]>::to_vec), + deflate: false, + } +} + +fn h5dump_available() -> bool { + Command::new("h5dump") + .arg("--version") + .output() + .map(|o| o.status.success()) + .unwrap_or(false) +} + +/// Write every case into one file with our writer, then check that our own +/// reader, h5py and h5dump (when installed) all return every value. Only the +/// libhdf5 half is skipped without h5py. +fn check_we_write(cases: &[WriteCase]) { + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("ours_chunk_index.h5"); + let path_str = path.display().to_string(); + + let mut b = FileBuilder::new(); + for c in cases { + let n: u64 = c.shape.iter().product(); + let data: Vec = (0..n as i32).collect(); + let ds = b.create_dataset(&c.name); + ds.with_i32_data(&data) + .with_shape(&c.shape) + .with_chunks(&c.chunks); + if let Some(ms) = &c.maxshape { + ds.with_maxshape(ms); + } + if c.deflate { + ds.with_deflate(4); + } + } + b.write(&path).unwrap(); + + // Our reader. + let file = File::open(&path).unwrap(); + for c in cases { + let got = file.dataset(&c.name).unwrap().read_i32().unwrap(); + let n: u64 = c.shape.iter().product(); + let bad = got + .iter() + .enumerate() + .filter(|&(i, &v)| v != i as i32) + .count(); + assert!( + got.len() == n as usize && bad == 0, + "{}: our reader: {bad} of {n} values wrong", + c.name + ); + } + + // libhdf5 via h5py. + skip_if_no_python!(); + let mut script = + format!("import h5py, numpy as np\nbad = []\nf = h5py.File(r'{path_str}', 'r')\n"); + for c in cases { + let shape: Vec = c.shape.iter().map(u64::to_string).collect(); + let maxshape: Vec = c + .maxshape + .as_ref() + .unwrap_or(&c.shape) + .iter() + .map(|&d| { + if d == u64::MAX { + "None".to_string() + } else { + d.to_string() + } + }) + .collect(); + script += &format!( + "d = f['{name}']\n\ + want = np.arange({n}, dtype='i4').reshape(({shape},))\n\ + got = d[()]\n\ + if d.maxshape != ({maxshape},): bad.append(('{name}', 'maxshape', d.maxshape))\n\ + elif not np.array_equal(got, want): \ + bad.append(('{name}', int((got != want).sum()), 'of', got.size))\n", + name = c.name, + n = c.shape.iter().product::(), + shape = shape.join(","), + maxshape = maxshape.join(","), + ); + } + script += "print(bad if bad else 'OK')\n"; + let out = run_python(&script); + assert_eq!(out, "OK", "h5py disagrees"); + + // libhdf5's own tool, when installed. + if h5dump_available() { + let o = Command::new("h5dump").arg(&path).output().unwrap(); + let stderr = String::from_utf8_lossy(&o.stderr); + assert!( + o.status.success() && !stderr.to_lowercase().contains("error"), + "h5dump failed: {stderr}" + ); + } + + // Let libhdf5 grow every resizable dataset by two chunks per dimension + // (capped at the maxshape) and rewrite it, which updates our index in + // place and inserts new chunks into it. Then both readers must agree. + let script = format!( + r#" +import h5py, numpy as np +grown = {{}} +with h5py.File(r'{path_str}', 'r+') as f: + for name in f: + d = f[name] + if d.chunks is None: + continue + new = tuple(s + 2 * c if m is None else min(m, s + 2 * c) + for s, m, c in zip(d.shape, d.maxshape, d.chunks)) + if new == d.shape: + continue + old = d[()] + full = np.full(new, -7, 'i4') + full[tuple(slice(0, s) for s in old.shape)] = old + d.resize(new) + d[...] = full + grown[name] = (list(old.shape), list(new)) +with h5py.File(r'{path_str}', 'r') as f: + for name, (old, new) in grown.items(): + want = np.full(new, -7, 'i4') + want[tuple(slice(0, s) for s in old)] = np.arange(int(np.prod(old)), dtype='i4').reshape(old) + assert np.array_equal(f[name][()], want), name +for name, (old, new) in grown.items(): + print(name, ','.join(map(str, old)), ','.join(map(str, new))) +"# + ); + let out = run_python(&script); + let growable = cases + .iter() + .filter(|c| c.maxshape.as_ref().is_some_and(|m| *m != c.shape)) + .count(); + assert_eq!(out.lines().count(), growable, "libhdf5 grew: {out}"); + let dims = |s: &str| -> Vec { s.split(',').map(|x| x.parse().unwrap()).collect() }; + let file = File::open(&path).unwrap(); + for line in out.lines() { + let mut parts = line.split(' '); + let (name, old, new) = ( + parts.next().unwrap(), + dims(parts.next().unwrap()), + dims(parts.next().unwrap()), + ); + let got = file.dataset(name).unwrap().read_i32().unwrap(); + let n: usize = new.iter().product(); + let mut want = vec![-7i32; n]; + for (flat, w) in want.iter_mut().enumerate() { + let mut rem = flat; + let mut coords = vec![0usize; new.len()]; + for d in (0..new.len()).rev() { + coords[d] = rem % new[d]; + rem /= new[d]; + } + if coords.iter().zip(&old).all(|(c, o)| c < o) { + *w = coords.iter().zip(&old).fold(0, |acc, (c, o)| acc * o + c) as i32; + } + } + let bad = got.iter().zip(&want).filter(|(a, b)| a != b).count(); + assert!( + got.len() == n && bad == 0, + "{name}: after libhdf5 grew it, our reader got {bad} of {n} values wrong" + ); + } +} + +/// A Fixed Array with more than 1024 elements must be paged, or libhdf5 +/// rejects the data block's checksum. +#[test] +fn we_write_paged_fixed_array() { + let mut cases: Vec = [1023u64, 1024, 1025, 2048, 5000] + .iter() + .map(|&n| wcase(&format!("fa_{n}"), &[n * 4], &[4], None)) + .collect(); + // Filtered elements are wider; a 2-D grid pages the same way. + let mut filtered = wcase("fa_1500_deflate", &[1500 * 4], &[4], None); + filtered.deflate = true; + cases.push(filtered); + cases.push(wcase("fa_2d_1100", &[110, 40], &[1, 4], None)); + check_we_write(&cases); +} + +/// An Extensible Array holds 4 elements in its index block and 240 in the +/// data blocks the index block addresses; everything after that lives under +/// super blocks, and from ~131K elements on in paged data blocks. Chunks past +/// index 243 used to be written but never indexed (read back as fill by us +/// and by libhdf5). +#[test] +fn we_write_extensible_array_past_index_block() { + let unl: &[u64] = &[u64::MAX]; + let mut cases: Vec = [1u64, 4, 5, 243, 244, 245, 300, 1000, 5000] + .iter() + .map(|&n| wcase(&format!("ea_{n}"), &[n * 4], &[4], Some(unl))) + .collect(); + let mut filtered = wcase("ea_300_deflate", &[300 * 4], &[4], Some(unl)); + filtered.deflate = true; + cases.push(filtered); + // Several super blocks and paged data blocks (level 13, the first with + // data blocks over 1024 elements, starts at element 4 + 131056). + cases.push(wcase("ea_140000", &[140_000], &[1], Some(unl))); + check_we_write(&cases); +} + +/// A maxshape larger than the shape: the index must be laid out over the +/// chunks of the maximum extent (libhdf5 read our Fixed Array past its end: +/// "addr overflow"), and an Extensible Array whose unlimited dimension is not +/// the first must swizzle it to the slowest position (libhdf5 read our +/// `(20, None)` dataset scrambled). +#[test] +fn we_write_maxshape_larger_than_shape() { + const U: u64 = u64::MAX; + let mut cases = vec![ + // Fixed Array over the maximum extent. + wcase("fa2d_finite_max", &[20, 30], &[5, 5], Some(&[40, 60])), + wcase("fa1d_finite_max", &[40], &[4], Some(&[100])), + wcase("fa3d_edges", &[6, 7, 8], &[4, 3, 5], Some(&[10, 9, 20])), + wcase("fa_paged_max", &[30, 50], &[1, 1], Some(&[40, 60])), + wcase("fa_one_chunk_now", &[5], &[5], Some(&[50])), + // Extensible Array, unlimited dimension first (no swizzle) ... + wcase("ea2d_unl_fin", &[20, 30], &[5, 5], Some(&[U, 30])), + wcase("ea2d_unl_fin_max", &[20, 30], &[5, 5], Some(&[U, 60])), + // ... and not first (swizzled). + wcase("ea2d_fin_unl", &[20, 30], &[5, 5], Some(&[20, U])), + wcase("ea2d_fin_max_unl", &[20, 30], &[5, 5], Some(&[40, U])), + wcase("ea3d_mid", &[6, 7, 8], &[4, 3, 5], Some(&[10, U, 20])), + // Past the index block and into super blocks, swizzled. + wcase("ea2d_many", &[3, 2000], &[1, 1], Some(&[4, U])), + ]; + let mut filtered = wcase( + "ea3d_last_deflate", + &[6, 7, 8], + &[4, 3, 5], + Some(&[6, 8, U]), + ); + filtered.deflate = true; + cases.push(filtered); + check_we_write(&cases); +} + +/// More than one unlimited dimension needs a version-2 B-tree chunk index, +/// as the library uses; an Extensible Array for `(None, None)` made libhdf5 +/// refuse the whole file ("already found unlimited dimension"). +#[test] +fn we_write_btree_v2_for_several_unlimited_dims() { + const U: u64 = u64::MAX; + let mut cases = vec![ + wcase("unl_unl", &[20, 30], &[5, 5], Some(&[U, U])), + wcase("unl_fin_unl", &[6, 7, 8], &[4, 3, 5], Some(&[U, 9, U])), + // More records than the library's 2048-byte node holds (84 here). + wcase("unl_unl_2400", &[40, 60], &[1, 1], Some(&[U, U])), + wcase("unl_unl_empty", &[0, 0], &[4, 4], Some(&[U, U])), + ]; + let mut filtered = wcase("unl_unl_deflate", &[6, 7, 8], &[4, 3, 5], Some(&[U, U, U])); + filtered.deflate = true; + cases.push(filtered); + check_we_write(&cases); +} + +/// A single-leaf B-tree has a 16-bit record count; beyond it the writer +/// refuses rather than writing a tree libhdf5 would misread. +#[test] +fn btree_v2_index_past_one_leaf_is_refused() { + let mut b = FileBuilder::new(); + b.create_dataset("d") + .with_i32_data(&vec![0i32; 70_000]) + .with_shape(&[70_000, 1]) + .with_chunks(&[1, 1]) + .with_maxshape(&[u64::MAX, u64::MAX]); + let dir = tempfile::tempdir().unwrap(); + assert!(b.write(dir.path().join("too_many.h5")).is_err()); +} + +/// A maxshape equal to the shape cannot grow, so it needs no chunks: the +/// dataset stays contiguous (as h5py makes it) unless chunks are requested. +#[test] +fn maxshape_equal_to_shape_stays_contiguous() { + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("ms_eq.h5"); + let data: Vec = (0..40).collect(); + let mut b = FileBuilder::new(); + b.create_dataset("plain") + .with_i32_data(&data) + .with_shape(&[40]) + .with_maxshape(&[40]); + b.create_dataset("chunked") + .with_i32_data(&data) + .with_shape(&[40]) + .with_maxshape(&[40]) + .with_chunks(&[8]); + b.write(&path).unwrap(); + + let file = File::open(&path).unwrap(); + let plain = file.dataset("plain").unwrap(); + assert_eq!(plain.read_i32().unwrap(), data); + assert_eq!(plain.max_dimensions().unwrap(), Some(vec![40])); + assert!( + plain.read_raw_ref().unwrap().is_some(), + "maxshape == shape should be contiguous" + ); + let chunked = file.dataset("chunked").unwrap(); + assert_eq!(chunked.read_i32().unwrap(), data); + assert!(chunked.read_raw_ref().unwrap().is_none()); + + skip_if_no_python!(); + let out = run_python(&format!( + "import h5py, numpy as np\n\ + f = h5py.File(r'{}', 'r')\n\ + for n in ('plain', 'chunked'):\n\ + \x20 d = f[n]\n\ + \x20 assert np.array_equal(d[()], np.arange(40, dtype='i4')), n\n\ + \x20 print(n, d.chunks, d.maxshape)\n", + path.display() + )); + assert_eq!(out, "plain None (40,)\nchunked (8,) (40,)"); +}