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<u8>, 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<WrittenChunk>]) -> 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<u8>,
+    slot: Option<&WrittenChunk>,
+    offset_size: u8,
+    chunk_size_bytes: Option<usize>,
+) {
+    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<WrittenChunk>],
     offset_size: u8,
     length_size: u8,
     has_filters: bool,
     fa_base_address: u64,
 ) -> Vec<u8> {
     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<Option<WrittenChunk>> = 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

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<u64>,
+    chunks: Vec<u64>,
+    maxshape: Option<Vec<u64>>,
+    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<i32> = (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<String> = c.shape.iter().map(u64::to_string).collect();
+        let maxshape: Vec<String> = 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::<u64>(),
+            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<WriteCase> = [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);
+}