diff --git a/crates/clawhdf5-format/src/chunked_write.rs b/crates/clawhdf5-format/src/chunked_write.rs index 6816546..c0d02e1 100644 --- a/crates/clawhdf5-format/src/chunked_write.rs +++ b/crates/clawhdf5-format/src/chunked_write.rs @@ -499,108 +499,140 @@ fn serialize_v4_fixed_array( 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; + +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 @@ -734,8 +766,9 @@ pub fn build_chunked_data_from_precompressed( ) } else { let fa_address = base_address + data_buf.len() as u64; + let slots: Vec> = written_chunks.iter().cloned().map(Some).collect(); let fa_bytes = build_fixed_array_at( - &written_chunks, + &slots, offset_size, length_size, pre.has_filters, @@ -747,7 +780,7 @@ pub fn build_chunked_data_from_precompressed( fa_address, offset_size, element_size as u32, - 10, // max_nelmts_bits — matches h5py convention + FA_PAGE_BITS, ) }; @@ -1382,7 +1415,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 diff --git a/crates/clawhdf5/tests/chunk_index_interop.rs b/crates/clawhdf5/tests/chunk_index_interop.rs index 9c030a3..cce2df6 100644 --- a/crates/clawhdf5/tests/chunk_index_interop.rs +++ b/crates/clawhdf5/tests/chunk_index_interop.rs @@ -11,7 +11,7 @@ use std::process::Command; -use clawhdf5::File; +use clawhdf5::{File, FileBuilder}; fn python() -> String { std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string()) @@ -254,3 +254,139 @@ fn h5py_fixed_array_partial_extent_reads_correctly() { }, ]); } + +// =========================================================================== +// 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}" + ); + } +} + +/// 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); +}