feat: extend same-file VDS assembly to N dimensions

Generalize Selection iteration from 1-D to arbitrary rank: iter_linear(dims)
enumerates a selection's row-major linear indices over a dataspace of the
given shape (ALL, NONE, regular hyperslabs, points), which is the order HDF5
uses to pair virtual and source selections.

read_virtual_data now passes the full virtual/source dimensions instead of a
single extent, so multi-dimensional block mappings scatter to the correct
non-contiguous linear positions. read_named_dataset_raw returns the source
dataset's dimensions. The rank-1 restriction is removed; only external-file
sources remain unsupported.

Tests: 2-D integration fixture (vds_2d_same_file.h5: two 2x2 sources placed
as non-contiguous blocks in a 4x4 virtual) plus N-D iter_linear unit tests
(block, strided, ALL, rank-mismatch). The 1-D path is unchanged.

Co-Authored-By: Claude Opus 4.8 <[email protected]>

CHANGELOG.md

@@ -15,17 +15,18 @@
   float members all read end-to-end, validated against HDF5 2.0.
 
 ### New Features
-- `clawhdf5-format`: assemble **1-D, same-file Virtual Datasets (VDS)**.
+- `clawhdf5-format`: assemble **same-file Virtual Datasets (VDS)** of any rank.
   Previously a virtual layout returned `UnsupportedVersion`. The reader now
   decodes the global-heap mapping block (reverse-engineered against HDF5 2.0:
   `version · nused · [source-file · source-dataset · source-selection ·
   virtual-selection]* · checksum`, including the block-version-1 same-file
   marker), decodes the `H5S` source/virtual dataspace **selections** (ALL,
   NONE, and version-3 regular hyperslabs), reads each same-file source dataset,
-  and scatters its selected elements into the virtual buffer; unmapped regions
+  and scatters its selected elements into the virtual buffer in row-major order
-  are left at the zero fill value. External-file sources and N-dimensional
+  (so multi-dimensional block mappings land correctly); unmapped regions are
-  selections return a clean unsupported error. The previous `parse_vds_mappings`
+  left at the zero fill value. External-file sources return a clean unsupported
-  used a guessed layout that did not match real files and is replaced.
+  error. The previous `parse_vds_mappings` used a guessed layout that did not
+  match real files and is replaced.
 
 ### Bug Fixes
 - `clawhdf5-format`: read **paged Fixed Array** chunk indexes. A filtered,

crates/clawhdf5-format/src/data_read.rs

@@ -385,13 +385,14 @@ pub fn read_raw_data_selection(
 
 /// Assemble a **Virtual Dataset (VDS)** from its source mappings.
 ///
-/// Supports **1-D, same-file** virtual datasets: each mapping's source dataset
+/// Supports **same-file** virtual datasets of any rank: each mapping's source
-/// is read from the same file and its selected elements are scattered into the
+/// dataset is read from the same file and its selected elements are scattered
-/// virtual buffer at the positions given by the virtual selection. Unmapped
+/// into the virtual buffer at the positions given by the virtual selection
-/// regions are left at the zero fill value.
+/// (both enumerated in row-major order, as HDF5 pairs them). Unmapped regions
+/// are left at the zero fill value.
 ///
-/// External-file sources and N-dimensional selections are reported as
+/// External-file sources are reported as unsupported rather than silently
-/// unsupported rather than silently producing wrong data.
+/// producing wrong data.
 fn read_virtual_data(
     file_data: &[u8],
     global_heap_address: Option<u64>,
@@ -409,16 +410,7 @@ fn read_virtual_data(
     let total_elems = dataspace.num_elements() as usize;
     let mut out = vec![0u8; total_elems.saturating_mul(elem_size)];
 
-    if dataspace.dimensions.len() > 1 {
+    let virtual_dims = &dataspace.dimensions;
-        return Err(FormatError::ChunkedReadError(
-            "N-dimensional virtual datasets are not supported".into(),
-        ));
-    }
-    let virtual_extent = dataspace
-        .dimensions
-        .first()
-        .copied()
-        .unwrap_or(total_elems as u64);
 
     let addr = global_heap_address.ok_or_else(|| {
         FormatError::ChunkedReadError("virtual dataset has no mapping global heap".into())
@@ -443,11 +435,11 @@ fn read_virtual_data(
         let (vsel, _) = Selection::decode_serialized(&m.virtual_selection)?;
         let (ssel, _) = Selection::decode_serialized(&m.source_selection)?;
 
-        let (src_raw, src_elems) =
+        let (src_raw, src_dims) =
             read_named_dataset_raw(file_data, &m.source_dataset, offset_size, length_size)?;
 
-        let vidx = vsel.iter_linear_1d(virtual_extent)?;
+        let vidx = vsel.iter_linear(virtual_dims)?;
-        let sidx = ssel.iter_linear_1d(src_elems as u64)?;
+        let sidx = ssel.iter_linear(&src_dims)?;
         if vidx.len() != sidx.len() {
             return Err(FormatError::ChunkedReadError(
                 "virtual/source selection element counts differ".into(),
@@ -468,14 +460,14 @@ fn read_virtual_data(
     Ok(out)
 }
 
-/// Read a named dataset's raw (decoded) bytes and element count, navigating
+/// Read a named dataset's raw (decoded) bytes and its dimensions, navigating
 /// from the superblock. Used to pull VDS source datasets out of the same file.
 fn read_named_dataset_raw(
     file_data: &[u8],
     path: &str,
     _offset_size: u8,
     _length_size: u8,
-) -> Result<(Vec<u8>, usize), FormatError> {
+) -> Result<(Vec<u8>, Vec<u64>), FormatError> {
     use crate::filter_pipeline::FilterPipeline;
     use crate::group_v2::resolve_path_any;
     use crate::message_type::MessageType;
@@ -511,7 +503,7 @@ fn read_named_dataset_raw(
         sb.offset_size,
         sb.length_size,
     )?;
-    Ok((raw, dataspace.num_elements() as usize))
+    Ok((raw, dataspace.dimensions.clone()))
 }
 
 /// Extract selected elements from a full dataset buffer.

crates/clawhdf5-format/src/selection.rs

@@ -258,13 +258,32 @@ impl Selection {
     }
 
     /// Enumerate the selected element indices of a **1-D** dataspace of the
-    /// given `extent`, in selection order.
+    /// given `extent`, in row-major selection order.
     ///
-    /// Returns an error for selections of rank != 1 (N-dimensional VDS
+    /// Convenience wrapper over [`Selection::iter_linear`] for rank-1 spaces.
-    /// assembly is not supported in this build).
     pub fn iter_linear_1d(&self, extent: u64) -> Result<Vec<u64>, FormatError> {
+        self.iter_linear(&[extent])
+    }
+
+    /// Enumerate the **row-major linear indices** of the selected elements of a
+    /// dataspace with shape `dims`, in row-major (C) iteration order.
+    ///
+    /// This is the order HDF5 uses to pair a virtual selection with a source
+    /// selection in a Virtual Dataset, so the i-th index returned here for the
+    /// virtual selection corresponds to the i-th index for the source
+    /// selection. Hyperslab/point selections whose rank differs from
+    /// `dims.len()` are rejected.
+    pub fn iter_linear(&self, dims: &[u64]) -> Result<Vec<u64>, FormatError> {
+        let total: u64 = dims.iter().product();
+        // Row-major strides: row_stride[d] = product(dims[d+1..]).
+        let rank = dims.len();
+        let mut row_stride = vec![1u64; rank];
+        for d in (0..rank.saturating_sub(1)).rev() {
+            row_stride[d] = row_stride[d + 1] * dims[d + 1];
+        }
+
         match self {
-            Selection::All => Ok((0..extent).collect()),
+            Selection::All => Ok((0..total).collect()),
             Selection::None => Ok(Vec::new()),
             Selection::Hyperslab {
                 start,
@@ -272,23 +291,54 @@ impl Selection {
                 count,
                 block,
             } => {
-                if start.len() != 1 {
+                if start.len() != rank {
                     return Err(FormatError::ChunkedReadError(
-                        "only 1-D VDS hyperslab selections are supported".into(),
+                        "VDS selection rank does not match dataspace rank".into(),
                     ));
                 }
-                let (s, st, c, b) = (start[0], stride[0], count[0], block[0]);
+                // Selected coordinates along each dimension, in order.
-                let mut out = Vec::new();
+                let mut per_dim: Vec<Vec<u64>> = Vec::with_capacity(rank);
-                for ci in 0..c {
+                for d in 0..rank {
-                    let base = s + ci * st;
+                    let mut coords = Vec::new();
-                    for bi in 0..b {
+                    for ci in 0..count[d] {
-                        let idx = base + bi;
+                        let base = start[d] + ci * stride[d];
-                        if idx >= extent {
+                        for bi in 0..block[d] {
-                            return Err(FormatError::ChunkedReadError(
+                            let coord = base + bi;
-                                "VDS hyperslab selection exceeds dataspace extent".into(),
+                            if coord >= dims[d] {
-                            ));
+                                return Err(FormatError::ChunkedReadError(
+                                    "VDS hyperslab selection exceeds dataspace extent".into(),
+                                ));
+                            }
+                            coords.push(coord);
                         }
-                        out.push(idx);
+                    }
+                    per_dim.push(coords);
+                }
+                if per_dim.iter().any(|c| c.is_empty()) {
+                    return Ok(Vec::new());
+                }
+                // Cartesian product in row-major order (dim 0 slowest-varying).
+                let out_len: usize = per_dim.iter().map(|c| c.len()).product();
+                let mut out = Vec::with_capacity(out_len);
+                let mut idx = vec![0usize; rank];
+                loop {
+                    let mut lin = 0u64;
+                    for d in 0..rank {
+                        lin += per_dim[d][idx[d]] * row_stride[d];
+                    }
+                    out.push(lin);
+                    // Increment the mixed-radix counter, last dimension fastest.
+                    let mut carry = true;
+                    for d in (0..rank).rev() {
+                        idx[d] += 1;
+                        if idx[d] < per_dim[d].len() {
+                            carry = false;
+                            break;
+                        }
+                        idx[d] = 0;
+                    }
+                    if carry {
+                        break;
                     }
                 }
                 Ok(out)
@@ -296,17 +346,21 @@ impl Selection {
             Selection::Points(pts) => {
                 let mut out = Vec::with_capacity(pts.len());
                 for p in pts {
-                    if p.len() != 1 {
+                    if p.len() != rank {
                         return Err(FormatError::ChunkedReadError(
-                            "only 1-D VDS point selections are supported".into(),
+                            "VDS point selection rank does not match dataspace rank".into(),
                         ));
                     }
-                    if p[0] >= extent {
+                    let mut lin = 0u64;
-                        return Err(FormatError::ChunkedReadError(
+                    for d in 0..rank {
-                            "VDS point selection exceeds dataspace extent".into(),
+                        if p[d] >= dims[d] {
-                        ));
+                            return Err(FormatError::ChunkedReadError(
+                                "VDS point selection exceeds dataspace extent".into(),
+                            ));
+                        }
+                        lin += p[d] * row_stride[d];
                     }
-                    out.push(p[0]);
+                    out.push(lin);
                 }
                 Ok(out)
             }
@@ -553,4 +607,57 @@ mod tests {
         let bytes = [0x02u8, 0, 0, 0, 0x03, 0, 0, 0, 0x00, 0x02, 0x01, 0, 0, 0];
         assert!(Selection::decode_serialized(&bytes).is_err());
     }
+
+    #[test]
+    fn iter_linear_2d_block_row_major() {
+        // A 2x2 block at the top-left of a 4x4 space => linear 0,1,4,5.
+        let sel = Selection::Hyperslab {
+            start: vec![0, 0],
+            stride: vec![1, 1],
+            count: vec![1, 1],
+            block: vec![2, 2],
+        };
+        assert_eq!(sel.iter_linear(&[4, 4]).unwrap(), vec![0, 1, 4, 5]);
+
+        // The same block shifted to the bottom-right => 10,11,14,15.
+        let sel2 = Selection::Hyperslab {
+            start: vec![2, 2],
+            stride: vec![1, 1],
+            count: vec![1, 1],
+            block: vec![2, 2],
+        };
+        assert_eq!(sel2.iter_linear(&[4, 4]).unwrap(), vec![10, 11, 14, 15]);
+    }
+
+    #[test]
+    fn iter_linear_2d_strided() {
+        // start=(0,0) stride=(2,2) count=(2,2) block=(1,1) over 4x4 =>
+        // coords (0,0)(0,2)(2,0)(2,2) => linear 0,2,8,10.
+        let sel = Selection::Hyperslab {
+            start: vec![0, 0],
+            stride: vec![2, 2],
+            count: vec![2, 2],
+            block: vec![1, 1],
+        };
+        assert_eq!(sel.iter_linear(&[4, 4]).unwrap(), vec![0, 2, 8, 10]);
+    }
+
+    #[test]
+    fn iter_linear_all_2d() {
+        assert_eq!(
+            Selection::All.iter_linear(&[2, 3]).unwrap(),
+            (0..6).collect::<Vec<_>>()
+        );
+    }
+
+    #[test]
+    fn iter_linear_rank_mismatch_rejected() {
+        let sel = Selection::Hyperslab {
+            start: vec![0],
+            stride: vec![1],
+            count: vec![1],
+            block: vec![2],
+        };
+        assert!(sel.iter_linear(&[4, 4]).is_err());
+    }
 }

crates/clawhdf5-format/tests/integration_test.rs

@@ -678,6 +678,23 @@ fn v4_virtual_dataset_same_file_read() {
     );
 }
 
+#[test]
+fn v4_virtual_dataset_2d_same_file_read() {
+    // A 4x4 virtual dataset assembled from two 2x2 same-file sources placed as
+    // non-contiguous blocks (exercises N-dimensional row-major scatter):
+    //   virt[0:2,0:2] <- src_a = [[1,2],[3,4]]
+    //   virt[2:4,2:4] <- src_b = [[5,6],[7,8]]
+    //   everything else -> fill 0
+    let file_data = include_bytes!("fixtures/vds_2d_same_file.h5");
+    let (raw, datatype, _) = read_chunked_dataset(file_data, "virt");
+    let values = read_as_i32(&raw, &datatype).unwrap();
+    assert_eq!(
+        values,
+        vec![1, 2, 0, 0, 3, 4, 0, 0, 0, 0, 5, 6, 0, 0, 7, 8],
+        "2-D VDS block scatter mismatch"
+    );
+}
+
 #[test]
 fn v4_paged_fixed_array_read() {
     // 1025 chunks of 16 int32s, gzip-filtered => Fixed Array index whose data