fix(format): read unmapped VDS elements as the virtual dataset's fill value

Elements of a virtual dataset that no mapping supplies (unmapped regions,
a missing source file, a missing source dataset) read as 0 instead of the
fill value libhdf5 returns — silent wrong data for any VDS created with a
non-zero fillvalue (read-matrix cases 0471/0472: -1 and 7 read as 0). A
missing source dataset was an error; libhdf5 reads it as fill.

Move VDS assembly into a new vds module following H5Dvirtual.c:
vds::read_virtual_dataset takes the dataset's fill value and a
VdsFileResolver that can refuse a name, and reports how many elements were
unmapped. Sources are read with their own fill value, and a source whose
datatype differs from the virtual dataset's is an error (libhdf5 converts).
File passes the dataset's fill value, resolves source names against the
virtual file's directory, and refuses names that leave it with an error
instead of reading them as fill. read_selection on a VDS goes through the
same fill-aware path.

The raw-read API (read_raw_data_full*) has no fill value, so it now errors
for a VDS with unmapped elements instead of guessing zeros.

Tests: vds_interop::vds_unmapped_regions_read_as_fill_value (external,
same-file, missing file/dataset, sparse source with its own fill, int
fill; earliest and latest format) and
vds_source_outside_directory_is_an_error_not_fill, both against h5py;
integration_test::v4_virtual_dataset_raw_api_refuses_to_guess_the_fill_value.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>

CHANGELOG.md

@@ -275,6 +275,21 @@
 ### Correctness
 - `clawhdf5-format` virtual datasets (VDS), checked against HDF5 2.0 through
   h5py (`crates/clawhdf5/tests/vds_interop.rs`):
+  - **Wrong data:** elements no mapping supplies — unmapped regions, and
+    mappings whose source file or dataset is missing — read as 0 instead of
+    the virtual dataset's fill value (e.g. h5py `fillvalue=-1`). Assembly moved
+    to the new `vds` module: `vds::read_virtual_dataset` takes the fill value
+    and a resolver that can refuse a name (`VdsFileResolver`), and `File`
+    passes the dataset's fill value. A missing source *dataset* read as an
+    error; it is fill now, as in libhdf5. Source datasets are read with their
+    own fill value for unallocated chunks, and a source whose datatype differs
+    from the virtual dataset's is an error (libhdf5 converts; we do not).
+    `File` now refuses a source name that leaves the virtual file's directory
+    (`../x.h5`, absolute paths), or any external source of a `File::from_bytes`
+    file, with an error — these used to read as fill.
+    **Behaviour change:** the raw-read API (`read_raw_data_full*`), which has
+    no fill value, now returns an error for a virtual dataset with unmapped
+    elements instead of zeros.
   - Hyperslab selection versions 1 and 2 were refused ("only version-3
     hyperslab selections are supported"). Version 1 is what libhdf5 writes for
     every VDS created with the default format bounds (h5py's default), so

crates/clawhdf5-format/src/data_read.rs

@@ -191,14 +191,9 @@ fn read_raw_data_full_impl(
             offset_size,
             length_size,
         ),
-        DataLayout::Virtual {
+        DataLayout::Virtual { .. } => read_virtual_data(
-            global_heap_address,
-            global_heap_index,
-            ..
-        } => read_virtual_data(
             file_data,
-            *global_heap_address,
+            layout,
-            *global_heap_index,
             dataspace,
             datatype,
             offset_size,
@@ -465,158 +460,54 @@ pub fn read_raw_data_selection(
     }
 }
 
-/// Assemble a **Virtual Dataset (VDS)** from its source mappings.
+/// Assemble a **Virtual Dataset (VDS)** through the raw-read API, which has no
+/// access to the dataset's fill value message.
 ///
-/// Supports virtual datasets of any rank. Same-file sources are read directly;
+/// Delegates to [`crate::vds::read_virtual_dataset`]. Because the fill value
-/// **external-file** sources are read through the caller-supplied `resolver`,
+/// is unknown here, a virtual dataset with any element no mapping supplies
-/// which maps a stored source file name to that file's bytes. Each mapping's
+/// (an unmapped region, or a missing source file or dataset) is an error
-/// selected source elements are scattered into the virtual buffer at the
+/// rather than a guess at the fill value; so is one whose extent libhdf5
-/// positions given by the virtual selection (both enumerated in row-major
+/// would report differently from the stored dataspace (unlimited mappings).
-/// order, as HDF5 pairs them). Unmapped regions are left at the zero fill value.
+/// Use [`crate::vds::read_virtual_dataset`] to read those.
-///
-/// A mapping whose external source file the resolver cannot supply (`None`) is
-/// skipped, leaving its region at fill — matching HDF5's tolerance of missing
-/// sources. An external source with no resolver at all is a hard error.
 #[allow(clippy::too_many_arguments)]
 fn read_virtual_data(
     file_data: &[u8],
-    global_heap_address: Option<u64>,
+    layout: &DataLayout,
-    global_heap_index: u32,
     dataspace: &Dataspace,
     datatype: &Datatype,
     offset_size: u8,
     length_size: u8,
     resolver: Option<&VdsSourceResolver>,
 ) -> Result<Vec<u8>, FormatError> {
-    use crate::data_layout::parse_vds_mappings;
+    let wrapped =
-    use crate::global_heap::GlobalHeapCollection;
+        resolver.map(|r| move |name: &str| -> Result<Option<Vec<u8>>, FormatError> { Ok(r(name)) });
-    use crate::selection::Selection;
+    let wrapped_ref = wrapped.as_ref().map(|w| w as &crate::vds::VdsFileResolver);
-
+    let v = crate::vds::read_virtual_dataset(
-    let elem_size = datatype.type_size() as usize;
-    let mut out = crate::chunked_read::alloc_output(crate::chunked_read::checked_byte_len(
-        dataspace.checked_num_elements()?,
-        elem_size,
-    )?)?;
-
-    let virtual_dims = &dataspace.dimensions;
-
-    let addr = global_heap_address.ok_or_else(|| {
-        FormatError::ChunkedReadError("virtual dataset has no mapping global heap".into())
-    })?;
-    let coll = GlobalHeapCollection::parse(file_data, addr as usize, length_size)?;
-    let obj =
-        coll.get_object(global_heap_index as u16)
-            .ok_or(FormatError::GlobalHeapObjectNotFound {
-                collection_address: addr,
-                index: global_heap_index as u16,
-            })?;
-    let mappings = parse_vds_mappings(&obj.data, length_size)?;
-
-    for m in &mappings {
-        let same_file = m.source_file.is_empty() || m.source_file == ".";
-
-        // Resolve the bytes of the file holding this source dataset.
-        let external;
-        let src_file_data: &[u8] = if same_file {
-            file_data
-        } else {
-            let r = resolver.ok_or_else(|| {
-                FormatError::ChunkedReadError(
-                    "external-file virtual dataset sources require a file resolver".into(),
-                )
-            })?;
-            match r(&m.source_file) {
-                Some(bytes) => {
-                    external = bytes;
-                    &external
-                }
-                // Source file unavailable: leave this region at fill value.
-                None => continue,
-            }
-        };
-
-        let (vsel, _) = Selection::decode_serialized(&m.virtual_selection)?;
-        let (ssel, _) = Selection::decode_serialized(&m.source_selection)?;
-
-        let (src_raw, src_dims) =
-            read_named_dataset_raw(src_file_data, &m.source_dataset, offset_size, length_size)?;
-
-        let vidx = vsel.iter_linear(virtual_dims)?;
-        let sidx = ssel.iter_linear(&src_dims)?;
-        if vidx.len() != sidx.len() {
-            return Err(FormatError::ChunkedReadError(
-                "virtual/source selection element counts differ".into(),
-            ));
-        }
-
-        for (&v, &s) in vidx.iter().zip(sidx.iter()) {
-            let (vo, so) = (v as usize * elem_size, s as usize * elem_size);
-            if vo + elem_size > out.len() || so + elem_size > src_raw.len() {
-                return Err(FormatError::ChunkedReadError(
-                    "virtual dataset selection out of bounds".into(),
-                ));
-            }
-            out[vo..vo + elem_size].copy_from_slice(&src_raw[so..so + elem_size]);
-        }
-    }
-
-    Ok(out)
-}
-
-/// 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>, Vec<u64>), FormatError> {
-    use crate::filter_pipeline::FilterPipeline;
-    use crate::group_v2::resolve_path_any;
-    use crate::message_type::MessageType;
-    use crate::object_header::ObjectHeader;
-    use crate::signature::find_signature;
-    use crate::superblock::Superblock;
-
-    let sig = find_signature(file_data)?;
-    let sb = Superblock::parse(file_data, sig)?;
-    let addr = resolve_path_any(file_data, &sb, path)?;
-    let hdr = ObjectHeader::parse(file_data, addr as usize, sb.offset_size, sb.length_size)?;
-
-    let find = |t: MessageType| hdr.messages.iter().find(|m| m.msg_type == t);
-    let ds_msg = find(MessageType::Dataspace)
-        .ok_or_else(|| FormatError::ChunkedReadError("VDS source has no dataspace".into()))?;
-    let dataspace = Dataspace::parse(&ds_msg.data, sb.length_size)?;
-    let dt_msg = find(MessageType::Datatype)
-        .ok_or_else(|| FormatError::ChunkedReadError("VDS source has no datatype".into()))?;
-    let (datatype, _) = Datatype::parse(&dt_msg.data)?;
-    let dl_msg = find(MessageType::DataLayout)
-        .ok_or_else(|| FormatError::ChunkedReadError("VDS source has no data layout".into()))?;
-    let layout = DataLayout::parse(&dl_msg.data, sb.offset_size, sb.length_size)?;
-    // A virtual dataset whose source is itself another virtual dataset could
-    // form a cycle (A -> B -> A) and recurse into a stack overflow. Nested
-    // virtual sources are exotic and unsupported, so stop here cleanly.
-    if matches!(layout, DataLayout::Virtual { .. }) {
-        return Err(FormatError::ChunkedReadError(
-            "virtual dataset source is itself virtual (unsupported)".into(),
-        ));
-    }
-    let pipeline = find(MessageType::FilterPipeline)
-        .map(|m| FilterPipeline::parse(&m.data))
-        .transpose()?;
-
-    let raw = read_raw_data_full(
         file_data,
-        &layout,
+        layout,
-        &dataspace,
+        dataspace,
-        &datatype,
+        datatype,
-        pipeline.as_ref(),
+        None,
-        sb.offset_size,
+        offset_size,
-        sb.length_size,
+        length_size,
+        wrapped_ref,
     )?;
-    Ok((raw, dataspace.dimensions.clone()))
+    if v.dims != dataspace.dimensions {
+        return Err(FormatError::ChunkedReadError(
+            "virtual dataset extent differs from its stored dataspace; \
+             read it with vds::read_virtual_dataset"
+                .into(),
+        ));
+    }
+    if v.unmapped > 0 {
+        return Err(FormatError::ChunkedReadError(
+            "virtual dataset has elements no source supplies, which read as its \
+             fill value; read it with vds::read_virtual_dataset and the fill value"
+                .into(),
+        ));
+    }
+    Ok(v.data)
 }
-
 /// Extract selected elements from a full dataset buffer.
 pub fn extract_selection_from_buffer(
     full_data: &[u8],

crates/clawhdf5-format/src/lib.rs

@@ -100,6 +100,7 @@ pub mod signature;
 pub mod superblock;
 pub mod symbol_table;
 pub mod type_builders;
+pub mod vds;
 pub mod vl_data;
 
 #[cfg(feature = "provenance")]

crates/clawhdf5-format/src/vds.rs

@@ -0,0 +1,497 @@
+//! Virtual Dataset (VDS) assembly, following libhdf5's `H5Dvirtual.c`.
+//!
+//! A virtual dataset stores no data of its own: a list of mappings (kept in
+//! the global heap) pairs a selection of the virtual dataspace with a
+//! selection of a *source* dataset, in the same file (`"."`) or another one.
+//! Reading it means reading each source and scattering the selected source
+//! elements into the virtual buffer, pairing the two selections element by
+//! element in row-major order. Elements no mapping supplies — unmapped
+//! regions, and mappings whose source file or dataset does not exist — read
+//! as the virtual dataset's **fill value**, as in libhdf5.
+//!
+//! Source files other than the virtual file itself are obtained through a
+//! caller-supplied [`VdsFileResolver`], since this crate has no filesystem.
+
+#[cfg(not(feature = "std"))]
+use alloc::{format, string::String, vec, vec::Vec};
+
+use crate::data_layout::{DataLayout, VdsMapping, parse_vds_mappings};
+use crate::dataspace::Dataspace;
+use crate::datatype::Datatype;
+use crate::error::FormatError;
+use crate::selection::{SerializedSelection, UNLIMITED};
+
+/// Resolves the name of an external VDS source file, as stored in the
+/// mapping, to that file's bytes.
+///
+/// `Ok(None)` means the file does not exist; its mappings then read as the
+/// fill value, as libhdf5 does for a missing source. `Err` refuses the name
+/// (e.g. a path the caller will not follow) and fails the read, so that a
+/// refused source is never passed off as fill.
+pub type VdsFileResolver<'a> = dyn Fn(&str) -> Result<Option<Vec<u8>>, FormatError> + 'a;
+
+/// A fully assembled virtual dataset.
+#[derive(Debug, Clone, PartialEq)]
+pub struct VirtualData {
+    /// The virtual dataset's extent.
+    pub dims: Vec<u64>,
+    /// Raw element bytes, row-major, in the virtual dataset's datatype.
+    pub data: Vec<u8>,
+    /// Number of elements no mapping supplied; they hold the fill value.
+    pub unmapped: u64,
+}
+
+fn vds_err(msg: impl Into<String>) -> FormatError {
+    FormatError::ChunkedReadError(msg.into())
+}
+
+/// One mapping with its selections decoded.
+struct Mapping {
+    file: String,
+    dataset: String,
+    vsel: SerializedSelection,
+    ssel: SerializedSelection,
+}
+
+/// Load and decode the mapping list of a virtual layout.
+fn load_mappings(
+    file_data: &[u8],
+    layout: &DataLayout,
+    length_size: u8,
+) -> Result<Vec<Mapping>, FormatError> {
+    let DataLayout::Virtual {
+        global_heap_address,
+        global_heap_index,
+        ..
+    } = layout
+    else {
+        return Err(vds_err("not a virtual dataset layout"));
+    };
+    let Some(addr) = *global_heap_address else {
+        return Ok(Vec::new());
+    };
+    let coll =
+        crate::global_heap::GlobalHeapCollection::parse(file_data, addr as usize, length_size)?;
+    let index = u16::try_from(*global_heap_index)
+        .map_err(|_| vds_err("VDS mapping heap index out of range"))?;
+    let obj = coll
+        .get_object(index)
+        .ok_or(FormatError::GlobalHeapObjectNotFound {
+            collection_address: addr,
+            index,
+        })?;
+    parse_vds_mappings(&obj.data, length_size)?
+        .into_iter()
+        .map(|m: VdsMapping| {
+            let (vsel, _) = SerializedSelection::decode(&m.virtual_selection)?;
+            let (ssel, _) = SerializedSelection::decode(&m.source_selection)?;
+            Ok(Mapping {
+                file: m.source_file,
+                dataset: m.source_dataset,
+                vsel,
+                ssel,
+            })
+        })
+        .collect()
+}
+
+/// The virtual dataset's extent as libhdf5 reports it (`H5Dget_space`).
+///
+/// For a virtual dataset whose mappings are all of fixed size this is the
+/// stored dataspace.
+pub fn virtual_dataset_extent(
+    file_data: &[u8],
+    layout: &DataLayout,
+    dataspace: &Dataspace,
+    _offset_size: u8,
+    length_size: u8,
+    _resolver: Option<&VdsFileResolver>,
+) -> Result<Vec<u64>, FormatError> {
+    let mappings = load_mappings(file_data, layout, length_size)?;
+    check_fixed(&mappings)?;
+    Ok(dataspace.dimensions.clone())
+}
+
+fn check_fixed(mappings: &[Mapping]) -> Result<(), FormatError> {
+    if mappings
+        .iter()
+        .any(|m| m.vsel.unlimited_dim().is_some() || m.ssel.unlimited_dim().is_some())
+    {
+        return Err(vds_err(
+            "unlimited virtual dataset mappings are not supported",
+        ));
+    }
+    Ok(())
+}
+
+/// Read a whole virtual dataset.
+///
+/// `fill` is the virtual dataset's fill value (from its fill value message;
+/// `None` for the default of zeros); every element no mapping supplies holds
+/// it. External source files are read through `resolver`; without one, a
+/// mapping to another file is an error.
+#[allow(clippy::too_many_arguments)]
+pub fn read_virtual_dataset(
+    file_data: &[u8],
+    layout: &DataLayout,
+    dataspace: &Dataspace,
+    datatype: &Datatype,
+    fill: Option<&[u8]>,
+    offset_size: u8,
+    length_size: u8,
+    resolver: Option<&VdsFileResolver>,
+) -> Result<VirtualData, FormatError> {
+    let mappings = load_mappings(file_data, layout, length_size)?;
+    check_fixed(&mappings)?;
+    let dims = virtual_dataset_extent(
+        file_data,
+        layout,
+        dataspace,
+        offset_size,
+        length_size,
+        resolver,
+    )?;
+
+    let elem_size = datatype.type_size() as usize;
+    let total = dims
+        .iter()
+        .try_fold(1u64, |acc, &d| acc.checked_mul(d))
+        .ok_or_else(|| FormatError::Overflow("virtual dataset extent".into()))?;
+    let mut data = crate::chunked_read::alloc_output(crate::chunked_read::checked_byte_len(
+        total, elem_size,
+    )?)?;
+    if let Some(fill) = fill.filter(|f| f.len() == elem_size && f.iter().any(|&b| b != 0)) {
+        for element in data.chunks_exact_mut(elem_size) {
+            element.copy_from_slice(fill);
+        }
+    }
+    let mut mapped = vec![false; usize::try_from(total).map_err(|_| vds_err("VDS too large"))?];
+
+    let mut sources = Sources::new(file_data, resolver);
+    for m in &mappings {
+        let Some(src) = sources.dataset(&m.file, &m.dataset, datatype)? else {
+            continue; // missing source file or dataset: fill
+        };
+        let vidx = selection_indices(&m.vsel, &dims, None)?;
+        let sidx = selection_indices(&m.ssel, &src.dims, None)?;
+        scatter(&mut data, &mut mapped, &src.raw, &vidx, &sidx, elem_size)?;
+    }
+
+    let unmapped = mapped.iter().filter(|&&m| !m).count() as u64;
+    Ok(VirtualData {
+        dims,
+        data,
+        unmapped,
+    })
+}
+
+/// Copy source element `sidx[i]` to virtual element `vidx[i]` for every `i`.
+fn scatter(
+    out: &mut [u8],
+    mapped: &mut [bool],
+    src: &[u8],
+    vidx: &[u64],
+    sidx: &[u64],
+    elem_size: usize,
+) -> Result<(), FormatError> {
+    if vidx.len() != sidx.len() {
+        return Err(vds_err("virtual/source selection element counts differ"));
+    }
+    for (&v, &s) in vidx.iter().zip(sidx) {
+        let (vo, so) = (v as usize * elem_size, s as usize * elem_size);
+        if vo + elem_size > out.len() || so + elem_size > src.len() {
+            return Err(vds_err("virtual dataset selection out of bounds"));
+        }
+        out[vo..vo + elem_size].copy_from_slice(&src[so..so + elem_size]);
+        mapped[v as usize] = true;
+    }
+    Ok(())
+}
+
+/// Row-major linear indices of the elements `sel` selects in a dataspace of
+/// shape `dims`, in the order libhdf5 iterates them (row-major).
+///
+/// `clip` = `(dim, limit)` drops every coordinate `>= limit` in `dim` — the
+/// clipping libhdf5 applies to an unlimited selection
+/// (`H5S_hyper_clip_unlim`); an unlimited selection must be clipped.
+fn selection_indices(
+    sel: &SerializedSelection,
+    dims: &[u64],
+    clip: Option<(usize, u64)>,
+) -> Result<Vec<u64>, FormatError> {
+    let overflow = || FormatError::Overflow("VDS selection index overflow".into());
+    let rank = dims.len();
+    let total = dims
+        .iter()
+        .try_fold(1u64, |acc, &d| acc.checked_mul(d))
+        .ok_or_else(overflow)?;
+    let mut row_stride = vec![1u64; rank];
+    for d in (0..rank.saturating_sub(1)).rev() {
+        row_stride[d] = row_stride[d + 1]
+            .checked_mul(dims[d + 1])
+            .ok_or_else(overflow)?;
+    }
+    if sel.rank().is_some_and(|r| r != rank) {
+        return Err(vds_err("VDS selection rank does not match dataspace rank"));
+    }
+    let limit = |d: usize| match clip {
+        Some((cd, l)) if cd == d => l,
+        _ => u64::MAX,
+    };
+
+    match sel {
+        SerializedSelection::All => Ok((0..total).collect()),
+        SerializedSelection::None => Ok(Vec::new()),
+        SerializedSelection::Regular {
+            start,
+            stride,
+            count,
+            block,
+        } => {
+            // Selected coordinates along each dimension, in order.
+            let mut per_dim: Vec<Vec<u64>> = Vec::with_capacity(rank);
+            for d in 0..rank {
+                let lim = limit(d);
+                if (count[d] == UNLIMITED || block[d] == UNLIMITED) && lim == u64::MAX {
+                    return Err(vds_err("unlimited VDS selection was not clipped"));
+                }
+                let mut coords = Vec::new();
+                let mut ci = 0u64;
+                'blocks: while ci < count[d] {
+                    let base = ci
+                        .checked_mul(stride[d])
+                        .and_then(|o| start[d].checked_add(o))
+                        .ok_or_else(overflow)?;
+                    if base >= lim {
+                        break;
+                    }
+                    let mut bi = 0u64;
+                    while bi < block[d] {
+                        let coord = base.checked_add(bi).ok_or_else(overflow)?;
+                        if coord >= lim {
+                            break 'blocks;
+                        }
+                        // Past the extent is malformed; bail before the list
+                        // can grow without bound.
+                        if coord >= dims[d] {
+                            return Err(vds_err("VDS selection exceeds the dataspace extent"));
+                        }
+                        coords.push(coord);
+                        bi += 1;
+                    }
+                    ci += 1;
+                }
+                per_dim.push(coords);
+            }
+            if per_dim.iter().any(|c| c.is_empty()) {
+                return Ok(Vec::new());
+            }
+            let n = per_dim
+                .iter()
+                .try_fold(1usize, |acc, c| acc.checked_mul(c.len()))
+                .ok_or_else(overflow)?;
+            let mut out = Vec::with_capacity(n);
+            let mut idx = vec![0usize; rank];
+            loop {
+                let lin: u64 = (0..rank).map(|d| per_dim[d][idx[d]] * row_stride[d]).sum();
+                out.push(lin);
+                // Mixed-radix increment, last dimension fastest.
+                let mut d = rank;
+                loop {
+                    if d == 0 {
+                        return Ok(out);
+                    }
+                    d -= 1;
+                    idx[d] += 1;
+                    if idx[d] < per_dim[d].len() {
+                        break;
+                    }
+                    idx[d] = 0;
+                }
+            }
+        }
+        SerializedSelection::Blocks {
+            rank: _,
+            starts,
+            ends,
+        } => {
+            // libhdf5 serializes the union as disjoint blocks, so their volumes
+            // never add up to more than the dataspace.
+            let mut volume = 0u64;
+            for (s, e) in starts.chunks_exact(rank).zip(ends.chunks_exact(rank)) {
+                for d in 0..rank {
+                    if e[d] >= dims[d] {
+                        return Err(vds_err("VDS selection exceeds the dataspace extent"));
+                    }
+                }
+                let v = s
+                    .iter()
+                    .zip(e)
+                    .try_fold(1u64, |acc, (&s, &e)| acc.checked_mul(e - s + 1))
+                    .ok_or_else(overflow)?;
+                volume = volume.checked_add(v).ok_or_else(overflow)?;
+                if volume > total {
+                    return Err(vds_err("VDS selection blocks overlap"));
+                }
+            }
+            let mut out = Vec::with_capacity(volume as usize);
+            for (s, e) in starts.chunks_exact(rank).zip(ends.chunks_exact(rank)) {
+                let mut cur = s.to_vec();
+                'block: loop {
+                    if (0..rank).all(|d| cur[d] < limit(d)) {
+                        out.push((0..rank).map(|d| cur[d] * row_stride[d]).sum());
+                    }
+                    for d in (0..rank).rev() {
+                        if cur[d] < e[d] {
+                            cur[d] += 1;
+                            continue 'block;
+                        }
+                        cur[d] = s[d];
+                    }
+                    break;
+                }
+            }
+            out.sort_unstable();
+            out.dedup();
+            Ok(out)
+        }
+    }
+}
+
+/// A source dataset's decoded contents.
+struct SourceData {
+    dims: Vec<u64>,
+    raw: Vec<u8>,
+}
+
+/// Source files and datasets, fetched on demand. The most recently used
+/// external file is kept, since consecutive mappings usually share one.
+struct Sources<'a, 'r> {
+    file_data: &'a [u8],
+    resolver: Option<&'r VdsFileResolver<'r>>,
+    cached_file: Option<(String, Option<Vec<u8>>)>,
+}
+
+impl<'a, 'r> Sources<'a, 'r> {
+    fn new(file_data: &'a [u8], resolver: Option<&'r VdsFileResolver<'r>>) -> Self {
+        Sources {
+            file_data,
+            resolver,
+            cached_file: None,
+        }
+    }
+
+    /// The bytes of source file `name`, or `None` if it does not exist.
+    fn file(&mut self, name: &str) -> Result<Option<&[u8]>, FormatError> {
+        if name == "." {
+            return Ok(Some(self.file_data));
+        }
+        if self.cached_file.as_ref().is_none_or(|(n, _)| n != name) {
+            let resolver = self.resolver.ok_or_else(|| {
+                vds_err("external-file virtual dataset sources require a file resolver")
+            })?;
+            self.cached_file = Some((String::from(name), resolver(name)?));
+        }
+        Ok(self.cached_file.as_ref().and_then(|(_, b)| b.as_deref()))
+    }
+
+    /// Read source dataset `path` from file `file`, or `None` when either
+    /// does not exist. Its datatype must be the virtual dataset's: libhdf5
+    /// converts between types here, which is not supported.
+    fn dataset(
+        &mut self,
+        file: &str,
+        path: &str,
+        datatype: &Datatype,
+    ) -> Result<Option<SourceData>, FormatError> {
+        let Some(bytes) = self.file(file)? else {
+            return Ok(None);
+        };
+        read_source(bytes, path, datatype)
+    }
+}
+
+/// Read source dataset `path` of the file in `file_data` in full (its own
+/// fill value applied to unallocated chunks), or `None` if it does not exist.
+fn read_source(
+    file_data: &[u8],
+    path: &str,
+    datatype: &Datatype,
+) -> Result<Option<SourceData>, FormatError> {
+    use crate::filter_pipeline::FilterPipeline;
+    use crate::message_type::MessageType;
+    use crate::object_header::ObjectHeader;
+    use crate::shared_message::message_data_with_sohm;
+
+    let sig = crate::signature::find_signature(file_data)?;
+    let sb = crate::superblock::Superblock::parse(file_data, sig)?;
+    let (os, ls) = (sb.offset_size, sb.length_size);
+    let addr = match crate::group_v2::resolve_path_any(file_data, &sb, path) {
+        Ok(a) => a,
+        Err(FormatError::PathNotFound(_)) => return Ok(None),
+        Err(e) => return Err(e),
+    };
+    let hdr = ObjectHeader::parse(file_data, addr as usize, os, ls)?;
+    let msg = |t: MessageType| {
+        hdr.messages
+            .iter()
+            .find(|m| m.msg_type == t)
+            .ok_or_else(|| vds_err(format!("VDS source {path} has no {t:?} message")))
+    };
+    let dataspace = Dataspace::parse(
+        &message_data_with_sohm(file_data, msg(MessageType::Dataspace)?, os, ls)?,
+        ls,
+    )?;
+    let (src_type, _) = Datatype::parse(&message_data_with_sohm(
+        file_data,
+        msg(MessageType::Datatype)?,
+        os,
+        ls,
+    )?)?;
+    if &src_type != datatype {
+        return Err(vds_err(format!(
+            "VDS source {path} has a different datatype from the virtual dataset \
+             (type conversion is not supported)"
+        )));
+    }
+    let layout = DataLayout::parse(&msg(MessageType::DataLayout)?.data, os, ls)?;
+    // A source that is itself virtual could form a cycle (A -> B -> A) and
+    // recurse without bound. Nested virtual sources are not supported.
+    if matches!(layout, DataLayout::Virtual { .. }) {
+        return Err(vds_err(
+            "virtual dataset source is itself virtual (unsupported)",
+        ));
+    }
+    let pipeline = hdr
+        .messages
+        .iter()
+        .find(|m| m.msg_type == MessageType::FilterPipeline)
+        .map(|m| {
+            message_data_with_sohm(file_data, m, os, ls).and_then(|d| FilterPipeline::parse(&d))
+        })
+        .transpose()?;
+    let raw = crate::fill_value::read_full_with_fill(
+        &hdr.messages,
+        file_data,
+        &layout,
+        &dataspace,
+        src_type.type_size() as usize,
+        os,
+        ls,
+        || {
+            crate::data_read::read_raw_data_full(
+                file_data,
+                &layout,
+                &dataspace,
+                &src_type,
+                pipeline.as_ref(),
+                os,
+                ls,
+            )
+        },
+    )?;
+    Ok(Some(SourceData {
+        dims: dataspace.dimensions,
+        raw,
+    }))
+}

crates/clawhdf5-format/tests/integration_test.rs

@@ -83,6 +83,45 @@ fn read_chunked_dataset(file_data: &[u8], dataset_path: &str) -> (Vec<u8>, Datat
     (raw, datatype, dataspace)
 }
 
+/// Helper: read a virtual dataset with `vds::read_virtual_dataset`, giving it
+/// the dataset's own fill value (same-file sources only).
+fn read_virtual_fixture(file_data: &[u8], path: &str) -> (Vec<u8>, Datatype) {
+    let sig = find_signature(file_data).unwrap();
+    let sb = Superblock::parse(file_data, sig).unwrap();
+    let addr = resolve_path_any(file_data, &sb, path).unwrap();
+    let hdr =
+        ObjectHeader::parse(file_data, addr as usize, sb.offset_size, sb.length_size).unwrap();
+    let msg = |t: MessageType| hdr.messages.iter().find(|m| m.msg_type == t).unwrap();
+    let ds = Dataspace::parse(&msg(MessageType::Dataspace).data, sb.length_size).unwrap();
+    let (dt, _) = Datatype::parse(&msg(MessageType::Datatype).data).unwrap();
+    let layout = DataLayout::parse(
+        &msg(MessageType::DataLayout).data,
+        sb.offset_size,
+        sb.length_size,
+    )
+    .unwrap();
+    let fill = clawhdf5_format::fill_value::dataset_fill_value_in(
+        file_data,
+        &hdr.messages,
+        sb.offset_size,
+        sb.length_size,
+    )
+    .unwrap();
+    let v = clawhdf5_format::vds::read_virtual_dataset(
+        file_data,
+        &layout,
+        &ds,
+        &dt,
+        fill.as_deref(),
+        sb.offset_size,
+        sb.length_size,
+        None,
+    )
+    .unwrap();
+    assert_eq!(v.dims, ds.dimensions);
+    (v.data, dt)
+}
+
 /// Helper: read any dataset (contiguous or chunked) as f64.
 fn read_dataset_f64_any(bytes: &[u8], path: &str) -> Vec<f64> {
     let sig = find_signature(bytes).unwrap();
@@ -672,7 +711,7 @@ fn v4_virtual_dataset_same_file_read() {
     //   virt[4:8]  <- (unmapped)                              => fill 0
     //   virt[8:12] <- src_b[0:4]  (ALL)                       => 20,21,22,23
     let file_data = include_bytes!("fixtures/vds_same_file.h5");
-    let (raw, datatype, _) = read_chunked_dataset(file_data, "virt");
+    let (raw, datatype) = read_virtual_fixture(file_data, "virt");
     let values = read_as_i32(&raw, &datatype).unwrap();
     assert_eq!(
         values,
@@ -681,6 +720,38 @@ fn v4_virtual_dataset_same_file_read() {
     );
 }
 
+#[test]
+fn v4_virtual_dataset_raw_api_refuses_to_guess_the_fill_value() {
+    // The raw read API has no fill value message, so a virtual dataset with an
+    // unmapped region is an error there instead of zeros that may be wrong.
+    let file_data = include_bytes!("fixtures/vds_same_file.h5");
+    let sig = find_signature(file_data).unwrap();
+    let sb = Superblock::parse(file_data, sig).unwrap();
+    let addr = resolve_path_any(file_data, &sb, "virt").unwrap();
+    let hdr =
+        ObjectHeader::parse(file_data, addr as usize, sb.offset_size, sb.length_size).unwrap();
+    let msg = |t: MessageType| hdr.messages.iter().find(|m| m.msg_type == t).unwrap();
+    let ds = Dataspace::parse(&msg(MessageType::Dataspace).data, sb.length_size).unwrap();
+    let (dt, _) = Datatype::parse(&msg(MessageType::Datatype).data).unwrap();
+    let layout = DataLayout::parse(
+        &msg(MessageType::DataLayout).data,
+        sb.offset_size,
+        sb.length_size,
+    )
+    .unwrap();
+    let err = read_raw_data_full(
+        file_data,
+        &layout,
+        &ds,
+        &dt,
+        None,
+        sb.offset_size,
+        sb.length_size,
+    )
+    .unwrap_err();
+    assert!(err.to_string().contains("fill value"), "{err}");
+}
+
 #[test]
 fn v4_virtual_dataset_2d_same_file_read() {
     // A 4x4 virtual dataset assembled from two 2x2 same-file sources placed as
@@ -689,7 +760,7 @@ fn v4_virtual_dataset_2d_same_file_read() {
     //   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 (raw, datatype) = read_virtual_fixture(file_data, "virt");
     let values = read_as_i32(&raw, &datatype).unwrap();
     assert_eq!(
         values,

crates/clawhdf5/src/reader.rs

@@ -485,6 +485,7 @@ impl<'f> Dataset<'f> {
             self.file.length_size(),
         )?;
         let fill_matters = !clawhdf5_format::fill_value::has_storage(&dl)
+            || matches!(dl, DataLayout::Virtual { .. })
             || (matches!(dl, DataLayout::Chunked { .. })
                 && !clawhdf5_format::fill_value::is_default(fill.as_deref()));
         if fill_matters {
@@ -837,24 +838,9 @@ impl<'f> Dataset<'f> {
         let pipeline = self.filter_pipeline()?;
 
         // Virtual datasets are assembled from source datasets; the per-file
-        // chunk cache does not apply. Route them through the resolver path so
+        // chunk cache does not apply.
-        // external sibling files resolve relative to this file's directory.
         if matches!(dl, DataLayout::Virtual { .. }) {
-            let base_dir = self.file.base_dir.clone();
+            return self.read_virtual(&dl, &ds, &dt);
-            let resolver = move |name: &str| -> Option<Vec<u8>> {
-                let dir = base_dir.as_ref()?;
-                std::fs::read(dir.join(sibling_file_name(name)?)).ok()
-            };
-            return Ok(data_read::read_raw_data_full_with_resolver(
-                self.file.data.as_bytes(),
-                &dl,
-                &ds,
-                &dt,
-                pipeline.as_ref(),
-                self.file.offset_size(),
-                self.file.length_size(),
-                Some(&resolver),
-            )?);
         }
 
         // Unallocated storage reads as the dataset's fill value.
@@ -880,6 +866,62 @@ impl<'f> Dataset<'f> {
             },
         )
     }
+
+    /// Resolver for external Virtual Dataset source files: names are
+    /// resolved against the directory of the file that holds the virtual
+    /// dataset, as libhdf5 does. A missing file is `Ok(None)` (its mappings
+    /// read as the fill value); a name that would leave that directory is
+    /// refused with an error rather than read as fill.
+    fn vds_resolver(&self) -> impl Fn(&str) -> Result<Option<Vec<u8>>, FormatError> + use<> {
+        let base_dir = self.file.base_dir.clone();
+        move |name: &str| {
+            let Some(dir) = base_dir.as_ref() else {
+                return Err(FormatError::ChunkedReadError(format!(
+                    "virtual dataset source file {name:?} cannot be resolved for an in-memory file"
+                )));
+            };
+            let rel = sibling_file_name(name).ok_or_else(|| {
+                FormatError::ChunkedReadError(format!(
+                    "virtual dataset source file {name:?} is outside the virtual file's \
+                     directory and is not followed"
+                ))
+            })?;
+            match std::fs::read(dir.join(rel)) {
+                Ok(bytes) => Ok(Some(bytes)),
+                Err(e) if e.kind() == std::io::ErrorKind::NotFound => Ok(None),
+                Err(e) => Err(FormatError::ChunkedReadError(format!(
+                    "cannot read virtual dataset source file {name:?}: {e}"
+                ))),
+            }
+        }
+    }
+
+    /// Read a whole virtual dataset; unmapped elements hold its fill value.
+    fn read_virtual(
+        &self,
+        dl: &DataLayout,
+        ds: &Dataspace,
+        dt: &Datatype,
+    ) -> Result<Vec<u8>, Error> {
+        let fill = clawhdf5_format::fill_value::dataset_fill_value_in(
+            self.file.data.as_bytes(),
+            &self.header.messages,
+            self.file.offset_size(),
+            self.file.length_size(),
+        )?;
+        let resolver = self.vds_resolver();
+        let v = clawhdf5_format::vds::read_virtual_dataset(
+            self.file.data.as_bytes(),
+            dl,
+            ds,
+            dt,
+            fill.as_deref(),
+            self.file.offset_size(),
+            self.file.length_size(),
+            Some(&resolver),
+        )?;
+        Ok(v.data)
+    }
 }
 
 // ---------------------------------------------------------------------------

crates/clawhdf5/tests/vds_interop.rs

@@ -190,3 +190,102 @@ expect("shared.h5", "v", "shared")
     );
     assert_matches_libhdf5(dir.path(), "shared.h5", "v", "shared");
 }
+
+// ---------------------------------------------------------------------------
+// Fill value
+// ---------------------------------------------------------------------------
+
+/// Elements no mapping supplies read as the virtual dataset's fill value, not
+/// as 0: unmapped regions, a missing source file, a missing source dataset.
+/// A source's own unallocated chunks read as *its* fill value.
+#[test]
+fn vds_unmapped_regions_read_as_fill_value() {
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    generate(
+        dir.path(),
+        r#"
+for i in range(3):
+    with h5py.File(f"src_{i}.h5", "w") as s:
+        s.create_dataset("data", data=np.arange(10.0) + i * 100)
+with h5py.File("sparse_src.h5", "w") as s:
+    d = s.create_dataset("data", shape=(10,), chunks=(5,), dtype="f8", fillvalue=42.0)
+    d[0:5] = np.arange(5.0) + 1000   # the second chunk is never written
+for libver in ["earliest", "latest"]:
+    with h5py.File(f"fill_{libver}.h5", "w", libver=libver) as f:
+        f.create_dataset("local", data=np.arange(10.0) * -1)
+        lay = h5py.VirtualLayout(shape=(6, 10), dtype="f8")
+        for i in range(3):
+            lay[i] = h5py.VirtualSource(f"src_{i}.h5", "data", shape=(10,))
+        lay[3] = h5py.VirtualSource("no_such_file.h5", "data", shape=(10,))
+        lay[4] = h5py.VirtualSource("src_0.h5", "no_such_dataset", shape=(10,))
+        # row 5 is not mapped at all
+        f.create_virtual_dataset("files", lay, fillvalue=-1.0)
+        lay = h5py.VirtualLayout(shape=(20,), dtype="f8")
+        lay[0:10] = h5py.VirtualSource(".", "local", shape=(10,))
+        f.create_virtual_dataset("same_file", lay, fillvalue=7.0)
+        lay = h5py.VirtualLayout(shape=(12,), dtype="f8")
+        lay[1:11] = h5py.VirtualSource("sparse_src.h5", "data", shape=(10,))
+        f.create_virtual_dataset("sparse_source", lay, fillvalue=-3.5)
+        lay = h5py.VirtualLayout(shape=(3, 4), dtype="i4")
+        lay[1, :] = h5py.VirtualSource(".", "ints", shape=(4,))
+        f.create_dataset("ints", data=np.arange(4, dtype="i4") + 1)
+        f.create_virtual_dataset("int_fill", lay, fillvalue=-99)
+    for name in ["files", "same_file", "sparse_source", "int_fill"]:
+        expect(f"fill_{libver}.h5", name, f"{name}_{libver}")
+"#,
+    );
+    for libver in ["earliest", "latest"] {
+        let file = format!("fill_{libver}.h5");
+        for name in ["files", "same_file", "sparse_source", "int_fill"] {
+            assert_matches_libhdf5(dir.path(), &file, name, &format!("{name}_{libver}"));
+        }
+    }
+
+    // A selection read goes through the same fill-aware assembly.
+    let f = File::open(dir.path().join("fill_latest.h5")).unwrap();
+    let sel = clawhdf5::Selection::slice(std::slice::from_ref(&(8..14)));
+    let got = f
+        .dataset("same_file")
+        .unwrap()
+        .read_f64_selection(&sel)
+        .unwrap();
+    assert_eq!(got, vec![-8.0, -9.0, 7.0, 7.0, 7.0, 7.0]);
+}
+
+/// A source name that would leave the virtual file's directory is refused
+/// with an error; it used to be skipped and read silently as fill.
+#[test]
+fn vds_source_outside_directory_is_an_error_not_fill() {
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    std::fs::create_dir(dir.path().join("sub")).unwrap();
+    generate(
+        dir.path(),
+        r#"
+with h5py.File("src.h5", "w") as s:
+    s.create_dataset("data", data=np.arange(4.0))
+with h5py.File("sub/up.h5", "w", libver="latest") as f:
+    lay = h5py.VirtualLayout(shape=(4,), dtype="f8")
+    lay[:] = h5py.VirtualSource("../src.h5", "data", shape=(4,))
+    f.create_virtual_dataset("v", lay, fillvalue=-1.0)
+with h5py.File("nested.h5", "w", libver="latest") as f:
+    lay = h5py.VirtualLayout(shape=(4,), dtype="f8")
+    lay[:] = h5py.VirtualSource("sub/inner.h5", "data", shape=(4,))
+    f.create_virtual_dataset("v", lay, fillvalue=-1.0)
+with h5py.File("sub/inner.h5", "w") as s:
+    s.create_dataset("data", data=np.arange(4.0) + 10)
+expect("nested.h5", "v", "nested")
+"#,
+    );
+    // libhdf5 resolves "../src.h5" (and would read [0, 1, 2, 3]); we refuse
+    // to leave the directory, and say so.
+    let f = File::open(dir.path().join("sub/up.h5")).unwrap();
+    let err = f.dataset("v").unwrap().read_f64().unwrap_err();
+    assert!(
+        err.to_string().contains("not followed"),
+        "unexpected error: {err}"
+    );
+    // A relative name below the virtual file's directory resolves there.
+    assert_matches_libhdf5(dir.path(), "nested.h5", "v", "nested");
+}

docs/known-issues.md

@@ -63,13 +63,15 @@ segfault or abort.
 
 ## Gaps found by the 2026-09-25 HDF5 audit (open)
 
-**Status:** open. These fail with an error; none returns wrong data, except
+**Status:** open. These fail with an error; none returns wrong data (the VDS
-the VDS item, which is marked.
+fill-value item that did is fixed).
 
 - **Layout message versions 1 and 2** (HDF5 1.6-era files): 84 of the 686
   sweep files, `InvalidLayoutVersion`. This is the largest single gap.
 - **Virtual datasets:**
-  - **Wrong data:** unmapped regions read as 0 instead of the fill value.
+  - ~~**Wrong data:** unmapped regions read as 0 instead of the fill value.~~
+    Fixed 2026-09-25: unmapped elements and missing sources read as the
+    virtual dataset's fill value.
   - `%b` printf-style source names are not expanded.
   - ~~Hyperslab selection versions 1 and 2 are refused.~~ Fixed 2026-09-25:
     versions 1-3 and irregular hyperslabs are decoded.