From 608de9ac9b2b9d35beb83e77c920c3e52483d6f0 Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Fri, 18 Sep 2026 04:51:53 +0000
Subject: [PATCH 01/11] feat(dataform): add chunk-shape policy for chunked
 feature storage

HDF5 reads whole chunks even when only a few of their elements are
selected, so the chunk shape is what decides the cost of a partial read.
Add choose_chunk_shape(), which derives a chunk shape from a byte budget
(1 MiB by default, per h5py guidance) instead of hard-coding an extent
along any one axis.

Both axes callers actually slice are chunked: axis 0 (samples) and axis 1
(channels/units). Trailing spatial axes are kept whole -- they are small
and never sliced, and splitting them would only shrink the chunks. The
budget is split evenly between the two chunked axes, with each absorbing
the headroom the other leaves, so a short sample axis does not force a
narrow feature chunk.
---
 bdpy/dataform/_feature_chunking.py      | 108 ++++++++++++++++++++++
 tests/dataform/test_feature_chunking.py | 117 ++++++++++++++++++++++++
 2 files changed, 225 insertions(+)
 create mode 100644 bdpy/dataform/_feature_chunking.py
 create mode 100644 tests/dataform/test_feature_chunking.py

diff --git a/bdpy/dataform/_feature_chunking.py b/bdpy/dataform/_feature_chunking.py
new file mode 100644
index 0000000..44b1b45
--- /dev/null
+++ b/bdpy/dataform/_feature_chunking.py
@@ -0,0 +1,108 @@
+"""Chunk-shape policy for chunked feature storage.
+
+HDF5 reads whole chunks even when only a few elements of a chunk are selected,
+so the chunk shape - not the file format alone - is what decides the cost of a
+partial read. This module turns an array shape and dtype into a chunk shape
+under a byte budget, rather than hard-coding an extent along any one axis.
+
+This file is a part of BdPy.
+"""
+
+import math
+from typing import Sequence, Tuple
+
+import numpy as np
+
+__all__ = ["DEFAULT_TARGET_CHUNK_BYTES", "choose_chunk_shape"]
+
+#: Default per-chunk byte budget. h5py's own guidance puts a useful chunk
+#: somewhere between 10 KiB and 1 MiB; features are large and read in slabs, so
+#: we sit at the top of that range.
+DEFAULT_TARGET_CHUNK_BYTES = 1 << 20
+
+
+def _prod(values: Sequence[int]) -> int:
+    out = 1
+    for v in values:
+        out *= int(v)
+    return out
+
+
+def choose_chunk_shape(
+    shape: Sequence[int],
+    dtype: np.dtype,
+    target_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
+    n_samples_known: bool = True,
+) -> Tuple[int, ...]:
+    """Choose an HDF5 chunk shape for a feature array under a byte budget.
+
+    The two axes callers actually slice are axis 0 (samples) and axis 1 (the
+    outermost feature axis, typically channels or units), so both are chunked.
+    Trailing axes - the spatial axes of a convolutional feature map, say
+    ``(13, 13)`` - are kept whole: they are small and are never sliced in
+    practice, and splitting them would only shrink the chunks.
+
+    Within the budget the two chunked axes are split roughly evenly in element
+    count, so that neither a sample-wise read nor a channel-wise read degenerates
+    into reading the entire dataset.
+
+    Parameters
+    ----------
+    shape : sequence of int
+        Full array shape, ``(n_samples, *feature_shape)``.
+    dtype : numpy.dtype
+        Array dtype; only its itemsize matters.
+    target_bytes : int, optional
+        Upper bound on the size of one chunk in bytes
+        (default: ``DEFAULT_TARGET_CHUNK_BYTES``).
+    n_samples_known : bool, optional
+        Whether ``shape[0]`` is the final number of samples. Pass ``False`` when
+        writing to a resizable dataset whose sample count is not yet known, so
+        that the sample-axis extent is not capped by the current size.
+
+    Returns
+    -------
+    tuple of int
+        Chunk shape, same length as ``shape``. Every element is at least 1 and,
+        when ``n_samples_known`` is True, at most the corresponding entry of
+        ``shape``.
+
+    Examples
+    --------
+    >>> choose_chunk_shape((1200, 1000), np.dtype(np.float32))
+    (512, 512)
+    >>> choose_chunk_shape((1200, 256, 13, 13), np.dtype(np.float32))
+    (39, 39, 13, 13)
+    >>> choose_chunk_shape((50, 1000), np.dtype(np.float32))
+    (50, 1000)
+    """
+    shape = tuple(int(s) for s in shape)
+    if len(shape) < 2:
+        raise ValueError(
+            "Feature arrays need at least a sample axis and a feature axis; "
+            "got shape {}".format(shape)
+        )
+    if target_bytes < 1:
+        raise ValueError("target_bytes must be positive, got {}".format(target_bytes))
+
+    itemsize = np.dtype(dtype).itemsize
+
+    # The whole array fits in one chunk: nothing to gain from splitting it.
+    if n_samples_known and _prod(shape) * itemsize <= target_bytes:
+        return shape
+
+    # Bytes taken by one (sample, feature-0) cell, i.e. one element of the two
+    # chunked axes with the trailing axes kept whole.
+    inner = max(1, _prod(shape[2:])) * itemsize
+    budget_cells = max(1, target_bytes // inner)
+
+    # Split the budget evenly between the two chunked axes, then let each axis
+    # absorb the headroom the other did not need: a short sample axis should not
+    # force a narrow feature chunk, and vice versa.
+    c1 = min(shape[1], max(1, math.isqrt(budget_cells)))
+    c0 = max(1, budget_cells // c1)
+    if n_samples_known:
+        c0 = min(shape[0], c0)
+    c1 = min(shape[1], max(c1, budget_cells // c0))
+
+    return (c0, c1, *shape[2:])
diff --git a/tests/dataform/test_feature_chunking.py b/tests/dataform/test_feature_chunking.py
new file mode 100644
index 0000000..579f70c
--- /dev/null
+++ b/tests/dataform/test_feature_chunking.py
@@ -0,0 +1,117 @@
+import unittest
+
+import numpy as np
+
+from bdpy.dataform._feature_chunking import (
+    DEFAULT_TARGET_CHUNK_BYTES,
+    choose_chunk_shape,
+)
+
+
+def _nbytes(chunk, dtype):
+    return int(np.prod(chunk)) * np.dtype(dtype).itemsize
+
+
+class TestChooseChunkShape(unittest.TestCase):
+    """Chunk shape decides partial-read cost, so pin the policy down."""
+
+    shapes = [
+        (1200, 1000),
+        (1200, 256, 13, 13),
+        (50, 1000),
+        (1, 256, 13, 13),
+        (3, 4096, 7, 7),
+        (1000000, 1),
+        (2, 1000000),
+        (17, 63, 5),
+    ]
+
+    def test_invariants(self):
+        for dtype in (np.float32, np.float64, np.int16):
+            for shape in self.shapes:
+                chunk = choose_chunk_shape(shape, np.dtype(dtype))
+                with self.subTest(shape=shape, dtype=dtype):
+                    self.assertEqual(len(chunk), len(shape))
+                    # Never zero, never larger than the dataset.
+                    for c, s in zip(chunk, shape):
+                        self.assertGreaterEqual(c, 1)
+                        self.assertLessEqual(c, s)
+                    # Trailing axes are kept whole.
+                    self.assertEqual(chunk[2:], tuple(shape[2:]))
+
+    def test_chunk_stays_within_budget(self):
+        dtype = np.dtype(np.float32)
+        for shape in self.shapes:
+            chunk = choose_chunk_shape(shape, dtype)
+            with self.subTest(shape=shape):
+                if _nbytes(shape, dtype) <= DEFAULT_TARGET_CHUNK_BYTES:
+                    # Small enough to store as a single chunk.
+                    self.assertEqual(chunk, tuple(shape))
+                else:
+                    self.assertLessEqual(
+                        _nbytes(chunk, dtype), DEFAULT_TARGET_CHUNK_BYTES
+                    )
+
+    def test_budget_is_actually_used(self):
+        # A chunk far below budget means needless read amplification: every
+        # oversized array should fill at least a quarter of the budget.
+        dtype = np.dtype(np.float32)
+        for shape in self.shapes:
+            if _nbytes(shape, dtype) <= DEFAULT_TARGET_CHUNK_BYTES:
+                continue
+            chunk = choose_chunk_shape(shape, dtype)
+            with self.subTest(shape=shape):
+                self.assertGreater(
+                    _nbytes(chunk, dtype), DEFAULT_TARGET_CHUNK_BYTES // 4
+                )
+
+    def test_both_sliceable_axes_are_chunked(self):
+        # The point of #144: a channel slice must not have to read every channel.
+        shape = (1200, 256, 13, 13)
+        chunk = choose_chunk_shape(shape, np.dtype(np.float32))
+        self.assertLess(chunk[0], shape[0])
+        self.assertLess(chunk[1], shape[1])
+
+    def test_whole_array_below_budget_is_one_chunk(self):
+        shape = (50, 1000)
+        self.assertEqual(
+            choose_chunk_shape(shape, np.dtype(np.float32)), shape
+        )
+
+    def test_single_sample(self):
+        chunk = choose_chunk_shape((1, 4096, 7, 7), np.dtype(np.float32))
+        self.assertEqual(chunk[0], 1)
+
+    def test_target_bytes_is_respected(self):
+        dtype = np.dtype(np.float32)
+        shape = (1200, 256, 13, 13)
+        small = choose_chunk_shape(shape, dtype, target_bytes=64 * 1024)
+        large = choose_chunk_shape(shape, dtype, target_bytes=4 * 1024 * 1024)
+        self.assertLessEqual(_nbytes(small, dtype), 64 * 1024)
+        self.assertLessEqual(_nbytes(large, dtype), 4 * 1024 * 1024)
+        self.assertLess(_nbytes(small, dtype), _nbytes(large, dtype))
+
+    def test_unknown_sample_count_is_not_capped(self):
+        # Resizable datasets start empty; the sample extent must come from the
+        # budget, not from the current (zero) size.
+        chunk = choose_chunk_shape(
+            (0, 256, 13, 13), np.dtype(np.float32), n_samples_known=False
+        )
+        self.assertGreater(chunk[0], 1)
+
+    def test_deterministic(self):
+        shape = (1200, 256, 13, 13)
+        dtype = np.dtype(np.float32)
+        self.assertEqual(
+            choose_chunk_shape(shape, dtype), choose_chunk_shape(shape, dtype)
+        )
+
+    def test_rejects_bad_input(self):
+        with self.assertRaises(ValueError):
+            choose_chunk_shape((100,), np.dtype(np.float32))
+        with self.assertRaises(ValueError):
+            choose_chunk_shape((100, 10), np.dtype(np.float32), target_bytes=0)
+
+
+if __name__ == "__main__":
+    unittest.main()

From 74230ad2b0ecf1a0db79a9ed65ada8652987ca9a Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Fri, 18 Sep 2026 04:55:54 +0000
Subject: [PATCH 02/11] feat(dataform): read chunked HDF5 features through a
 FeatureStore backend

Features was welded to one on-disk layout: a directory per layer holding
one .mat file per stimulus. That layout has a file boundary only on the
sample axis, so reading a few channels still costs a full read of every
stimulus (issue #144).

Introduce a FeatureStore interface so the layout is a backend detail:

- MatFeatureStore carries today's behaviour unchanged (per-label paths,
  key 'feat', _mat_v73.loadmat_key, the multiprocessing fan-out).
- HDF5FeatureStore reads <layer>.h5 holding /features and /labels, with
  /features explicitly chunked, and pushes feature-axis slices down to
  h5py so only the covered chunks are read. Root attributes bdpy_format
  and bdpy_format_version identify the file and are validated on open.

Features holds one store per dpath and dispatches reads to whichever
store owns each label, so a mix of layouts across dpath entries works and
a future backend needs no change here. The format is auto-detected per
directory and can be forced with format='mat'/'hdf5'.

New read API, identical on both backends:

    features.get(layer, label=..., feature_slice=np.s_[128:256])
    features.iter_chunks(layer, axis=1)

feature_slice is a logical slice rather than the get_chunk() of the
issue, so the HDF5 chunk concept stays out of the public API; on the
legacy backend it loads in full and slices, on HDF5 it is a real partial
read. iter_chunks yields (slice, block) so callers can place results back
without tracking offsets, and defaults to the on-disk chunk extent so
every element is read exactly once.

Existing behaviour is untouched: labels, index, layers, feature_index,
get, get_features, statistic, the single-slot layer cache and
statistics.pkl all behave exactly as before, including caller-order row
selection. Sliced reads deliberately bypass the layer cache, which exists
to hold a whole layer.
---
 bdpy/dataform/_feature_store.py | 580 ++++++++++++++++++++++++++++++++
 bdpy/dataform/features.py       | 295 +++++++++++-----
 2 files changed, 788 insertions(+), 87 deletions(-)
 create mode 100644 bdpy/dataform/_feature_store.py

diff --git a/bdpy/dataform/_feature_store.py b/bdpy/dataform/_feature_store.py
new file mode 100644
index 0000000..38e2112
--- /dev/null
+++ b/bdpy/dataform/_feature_store.py
@@ -0,0 +1,580 @@
+"""Storage backends for DNN feature collections.
+
+``Features`` used to be welded to one on-disk layout: a directory per layer
+holding one ``.mat`` file per stimulus. That layout has a file boundary only on
+the sample axis, so reading a few channels still costs a full read of every
+stimulus. This module puts a backend interface between ``Features`` and the
+bytes, so the legacy layout and the chunked HDF5 layout are interchangeable and
+a future backend needs no change to ``Features`` itself.
+
+Every backend answers the same question -- "give me these labels, sliced this
+way along the feature axes" -- and differs only in how much it has to read to do
+it. The legacy backend loads whole stimulus files and slices in memory; the HDF5
+backend pushes the slice down to h5py and reads only the chunks it covers.
+
+This file is a part of BdPy.
+"""
+
+import glob
+import os
+from abc import ABC, abstractmethod
+from functools import partial
+from multiprocessing import Pool
+from typing import Dict, Iterator, List, Optional, Sequence, Tuple, Union
+
+import h5py
+import numpy as np
+
+from . import _mat_v73
+from ._feature_chunking import DEFAULT_TARGET_CHUNK_BYTES
+
+__all__ = [
+    "FeatureStore",
+    "HDF5FeatureStore",
+    "MatFeatureStore",
+    "detect_format",
+]
+
+#: Root attribute marking a file as bdpy chunked feature storage.
+FORMAT_ATTR = "bdpy_format"
+#: Root attribute holding the schema version.
+FORMAT_VERSION_ATTR = "bdpy_format_version"
+#: Value of ``FORMAT_ATTR`` for feature files.
+FORMAT_NAME = "features"
+#: Highest schema version this module can read.
+SUPPORTED_FORMAT_VERSION = 1
+#: Dataset holding the feature array, shape ``(n_samples, *feature_shape)``.
+FEATURES_DATASET = "features"
+#: Dataset holding the stimulus labels, shape ``(n_samples,)``.
+LABELS_DATASET = "labels"
+#: Extension of a chunked feature file.
+HDF5_EXT = "h5"
+
+# A slice spec for the feature axes: what ``numpy.s_[...]`` produces.
+FeatureSlice = Union[slice, int, Sequence[int], np.ndarray, Tuple, None]
+
+
+def _load_array_with_key(key: str, path: str) -> np.ndarray:
+    # v5 .mat via scipy, v7.3 (HDF5) via h5py; avoids hdf5storage on the load
+    # path, which breaks under NumPy 2.0 (see bdpy/dataform/_mat_v73.py).
+    # NOTE: key comes first so that partial(..., 'feat') is a Pool.map callable.
+    return _mat_v73.loadmat_key(path, key)
+
+
+def _determine_num_parallel(num_files: int) -> int:
+    # NOTE: optimal number of parallel processes is not clear. It could depend
+    # on several factors such as the number of files, the size of files, the
+    # number of cores, etc. For now, we use a simple heuristic based on the
+    # number of files.
+    num_parallel: int
+    if num_files < 16:
+        num_parallel = 1
+    elif num_files < 64:
+        num_parallel = 16
+    else:
+        num_parallel = 64
+    return num_parallel
+
+
+def _normalize_feature_slice(feature_slice: FeatureSlice) -> Tuple:
+    """Normalize a feature-axis slice spec to a tuple of per-axis indexers."""
+    if feature_slice is None:
+        return ()
+    if isinstance(feature_slice, tuple):
+        return feature_slice
+    return (feature_slice,)
+
+
+def _is_fancy(indexer: object) -> bool:
+    """Whether an indexer is a fancy (list/array) index rather than a slice."""
+    if isinstance(indexer, (slice, int, np.integer)):
+        return False
+    return isinstance(indexer, (list, tuple, np.ndarray, range))
+
+
+class FeatureStore(ABC):
+    """Interface to a collection of DNN features on disk.
+
+    A store exposes layers and stimulus labels, and reads features for a set of
+    labels with an optional slice along the feature axes. Implementations differ
+    only in how much I/O a given request costs.
+    """
+
+    @property
+    @abstractmethod
+    def layers(self) -> List[str]:
+        """List of DNN layers held by this store."""
+
+    @property
+    @abstractmethod
+    def labels(self) -> List[str]:
+        """List of stimulus labels, in the store's own order."""
+
+    @abstractmethod
+    def shape(self, layer: str) -> Tuple[int, ...]:
+        """Full shape ``(n_samples, *feature_shape)`` of `layer`."""
+
+    @abstractmethod
+    def dtype(self, layer: str) -> np.dtype:
+        """Dtype of the features in `layer`."""
+
+    @abstractmethod
+    def read(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]] = None,
+        feature_slice: FeatureSlice = None,
+    ) -> np.ndarray:
+        """Read features from `layer`.
+
+        Parameters
+        ----------
+        layer : str
+            DNN layer.
+        labels : sequence of str, optional
+            Stimulus labels to read. ``None`` reads every label in store order.
+            Otherwise rows come back **in the order given**, and repeated labels
+            yield repeated rows.
+        feature_slice : slice, int, array-like or tuple, optional
+            Index applied to the feature axes (axes 1 and up), as produced by
+            ``numpy.s_[...]``. ``None`` reads the whole feature tensor.
+
+        Returns
+        -------
+        numpy.ndarray
+            Array of shape ``(n_labels, *sliced_feature_shape)``.
+        """
+
+    def chunk_extent(self, layer: str, axis: int) -> Optional[int]:
+        """On-disk chunk extent of `layer` along `axis`, if the store has one.
+
+        Returns ``None`` when the backend has no chunking to align with, in
+        which case :meth:`iter_chunks` falls back to a byte budget.
+        """
+        return None
+
+    def iter_chunks(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]] = None,
+        feature_slice: FeatureSlice = None,
+        axis: int = 1,
+        size: Optional[int] = None,
+    ) -> Iterator[Tuple[slice, np.ndarray]]:
+        """Iterate over `layer` in slabs along `axis`.
+
+        Yields ``(sl, block)`` pairs, where `sl` is the slice applied to `axis`
+        of the selected array and `block` is that slab. The slice is yielded so
+        that callers can place results back without tracking offsets themselves.
+
+        Parameters
+        ----------
+        layer : str
+            DNN layer.
+        labels : sequence of str, optional
+            Stimulus labels, as in :meth:`read`.
+        feature_slice : slice, int, array-like or tuple, optional
+            Index applied to the feature axes before iterating.
+        axis : int, optional
+            Axis of the *selected* array to iterate over. Axis 0 is the sample
+            axis; the default, axis 1, is the outermost feature axis, which is
+            the one partial feature reads are about.
+        size : int, optional
+            Number of elements per slab. Defaults to the on-disk chunk extent
+            along `axis` when the backend has one, so that iteration is
+            chunk-aligned and every element is read exactly once; otherwise a
+            size derived from the same byte budget is used.
+
+        Yields
+        ------
+        tuple of (slice, numpy.ndarray)
+            The slice applied to `axis`, and the corresponding slab.
+        """
+        selected_shape = self._selected_shape(layer, labels, feature_slice)
+        if axis < 0:
+            axis += len(selected_shape)
+        if not 0 <= axis < len(selected_shape):
+            raise ValueError(
+                "axis {} is out of range for selected shape {}".format(
+                    axis, selected_shape
+                )
+            )
+
+        length = selected_shape[axis]
+        if size is None:
+            size = self._default_iter_size(layer, axis, selected_shape)
+        if size < 1:
+            raise ValueError("size must be positive, got {}".format(size))
+
+        base = _normalize_feature_slice(feature_slice)
+        for start in range(0, length, size):
+            sl = slice(start, min(start + size, length))
+            if axis == 0:
+                block_labels = (
+                    list(self.labels) if labels is None else list(labels)
+                )[sl]
+                yield sl, self.read(layer, block_labels, feature_slice)
+            else:
+                yield sl, self.read(
+                    layer, labels, _compose_slice(base, axis - 1, sl)
+                )
+
+    def _selected_shape(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]],
+        feature_slice: FeatureSlice,
+    ) -> Tuple[int, ...]:
+        """Shape the selection would have, without reading the data."""
+        full = self.shape(layer)
+        n_samples = full[0] if labels is None else len(labels)
+        feature_shape = _sliced_shape(full[1:], feature_slice)
+        return (n_samples, *feature_shape)
+
+    def _default_iter_size(
+        self, layer: str, axis: int, selected_shape: Tuple[int, ...]
+    ) -> int:
+        extent = self.chunk_extent(layer, axis)
+        if extent is not None:
+            return max(1, min(extent, selected_shape[axis]))
+        # No on-disk chunking to align with: fall back to the byte budget.
+        itemsize = np.dtype(self.dtype(layer)).itemsize
+        other = 1
+        for i, n in enumerate(selected_shape):
+            if i != axis:
+                other *= int(n)
+        per_element = max(1, other * itemsize)
+        return max(1, min(selected_shape[axis], DEFAULT_TARGET_CHUNK_BYTES // per_element))
+
+
+def _compose_slice(base: Tuple, axis: int, sl: slice) -> Tuple:
+    """Return `base` with `sl` intersected into its entry for `axis`.
+
+    `axis` indexes the *feature* axes (axis 0 here is array axis 1), and `sl`
+    is relative to the result of applying `base`.
+    """
+    out = list(base) + [slice(None)] * max(0, axis + 1 - len(base))
+    current = out[axis]
+    if isinstance(current, slice) and current == slice(None):
+        out[axis] = sl
+    elif isinstance(current, slice):
+        # Compose the two slices by indexing a range with both.
+        out[axis] = _compose_two_slices(current, sl)
+    else:
+        out[axis] = np.asarray(current)[sl]
+    return tuple(out)
+
+
+def _compose_two_slices(outer: slice, inner: slice) -> slice:
+    """Compose ``x[outer][inner]`` into a single slice where possible."""
+    # Only the simple (positive-step) case arises from iter_chunks, which
+    # always produces contiguous forward slices.
+    ostep = 1 if outer.step is None else outer.step
+    if ostep < 0:
+        # Fall back to an explicit index array for reversed slices.
+        raise ValueError("cannot compose a reversed slice; pass an explicit size")
+    ostart = 0 if outer.start is None else outer.start
+    start = ostart + (0 if inner.start is None else inner.start) * ostep
+    stop = ostart + (inner.stop if inner.stop is not None else 0) * ostep
+    if outer.stop is not None:
+        stop = min(stop, outer.stop)
+    return slice(start, stop, ostep)
+
+
+def _sliced_shape(
+    feature_shape: Tuple[int, ...], feature_slice: FeatureSlice
+) -> Tuple[int, ...]:
+    """Shape of ``numpy.empty(feature_shape)[feature_slice]``, cheaply."""
+    indexers = _normalize_feature_slice(feature_slice)
+    if not indexers:
+        return feature_shape
+    # A broadcast dummy allocates nothing but lets NumPy do the index
+    # arithmetic, including negative steps, ellipses and newaxis.
+    dummy = np.broadcast_to(np.int8(0), feature_shape)
+    return tuple(int(n) for n in dummy[indexers].shape)
+
+
+class MatFeatureStore(FeatureStore):
+    """Legacy per-stimulus ``.mat`` feature directory.
+
+    Layout is ``<dpath>/<layer>/<label>.<ext>``, each file holding one sample
+    under `key` with a leading sample axis. There is no file boundary on the
+    feature axes, so `feature_slice` is applied after the full stimulus files
+    have been loaded and concatenated -- the read cost is the same as before,
+    and the slice only saves the caller from doing it themselves.
+
+    Parameters
+    ----------
+    dpath : str
+        Feature directory.
+    ext : str, optional
+        Feature file extension (default: ``'mat'``).
+    key : str, optional
+        Variable name inside each file (default: ``'feat'``).
+    """
+
+    def __init__(self, dpath: str, ext: str = "mat", key: str = "feat"):
+        self._dpath = dpath
+        self._ext = ext
+        self._key = key
+        self._layers = self._collect_layers()
+        self._labels = self._collect_labels()
+        self._file_table: Dict[str, Dict[str, str]] = {
+            layer: {
+                label: os.path.join(dpath, layer, label + "." + ext)
+                for label in self._labels
+            }
+            for layer in self._layers
+        }
+
+    @property
+    def layers(self) -> List[str]:
+        return self._layers
+
+    @property
+    def labels(self) -> List[str]:
+        return self._labels
+
+    def path(self, layer: str, label: str) -> str:
+        """Path of the file holding `label` in `layer`."""
+        return self._file_table[layer][label]
+
+    def shape(self, layer: str) -> Tuple[int, ...]:
+        if not self._labels:
+            raise RuntimeError("No features found in {}".format(self._dpath))
+        sample = _load_array_with_key(self._key, self.path(layer, self._labels[0]))
+        return (len(self._labels), *sample.shape[1:])
+
+    def dtype(self, layer: str) -> np.dtype:
+        if not self._labels:
+            raise RuntimeError("No features found in {}".format(self._dpath))
+        sample = _load_array_with_key(self._key, self.path(layer, self._labels[0]))
+        return sample.dtype
+
+    def read(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]] = None,
+        feature_slice: FeatureSlice = None,
+    ) -> np.ndarray:
+        if labels is None:
+            labels = self._labels
+        paths = [self.path(layer, label) for label in labels]
+
+        num_parallel = _determine_num_parallel(len(paths))
+        load = partial(_load_array_with_key, self._key)
+        if num_parallel == 1:
+            arrays = list(map(load, paths))
+        else:
+            with Pool(processes=num_parallel) as pool:
+                arrays = pool.map(load, paths)
+        features = np.concatenate(arrays, axis=0)
+
+        indexers = _normalize_feature_slice(feature_slice)
+        if indexers:
+            features = features[(slice(None), *indexers)]
+        return features
+
+    def _collect_layers(self) -> List[str]:
+        return sorted(
+            d
+            for d in os.listdir(self._dpath)
+            if os.path.isdir(os.path.join(self._dpath, d))
+        )
+
+    def _collect_labels(self) -> List[str]:
+        labels: List[str] = []
+        for layer in self._layers:
+            layer_dir = os.path.join(self._dpath, layer)
+            layer_dir = layer_dir.replace("[", "[[]")  # Use glob.escape for Python 3.4 or later
+            files = glob.glob(os.path.join(layer_dir, "*." + self._ext))
+            labels_t = sorted(
+                os.path.splitext(os.path.basename(f))[0] for f in files
+            )
+            if not labels:
+                labels = labels_t
+            elif labels != labels_t:
+                raise RuntimeError("Invalid feature file in %s " % self._dpath)
+        return labels
+
+
+class HDF5FeatureStore(FeatureStore):
+    """Chunked HDF5 feature directory, one ``<layer>.h5`` file per layer.
+
+    Each file holds ``/features`` with shape ``(n_samples, *feature_shape)`` and
+    ``/labels`` with the matching stimulus labels. ``/features`` is explicitly
+    chunked, so a slice along the feature axes reads only the chunks it covers
+    instead of the whole layer.
+
+    Parameters
+    ----------
+    dpath : str
+        Directory holding ``<layer>.h5`` files.
+    """
+
+    def __init__(self, dpath: str):
+        self._dpath = dpath
+        self._layers = sorted(
+            os.path.splitext(os.path.basename(p))[0]
+            for p in glob.glob(os.path.join(glob.escape(dpath), "*." + HDF5_EXT))
+        )
+        if not self._layers:
+            raise RuntimeError("No .{} feature file found in {}".format(HDF5_EXT, dpath))
+        self._labels: Optional[List[str]] = None
+        self._label_index: Optional[Dict[str, int]] = None
+        self._validate_all()
+
+    @property
+    def layers(self) -> List[str]:
+        return self._layers
+
+    @property
+    def labels(self) -> List[str]:
+        if self._labels is None:
+            with self._open(self._layers[0]) as f:
+                self._labels = _decode_labels(f[LABELS_DATASET][()])
+            self._label_index = {
+                label: i for i, label in enumerate(self._labels)
+            }
+        return self._labels
+
+    def path(self, layer: str) -> str:
+        """Path of the file holding `layer`."""
+        return os.path.join(self._dpath, layer + "." + HDF5_EXT)
+
+    def shape(self, layer: str) -> Tuple[int, ...]:
+        with self._open(layer) as f:
+            return tuple(int(n) for n in f[FEATURES_DATASET].shape)
+
+    def dtype(self, layer: str) -> np.dtype:
+        with self._open(layer) as f:
+            dtype: np.dtype = np.dtype(f[FEATURES_DATASET].dtype)
+        return dtype
+
+    def chunk_extent(self, layer: str, axis: int) -> Optional[int]:
+        with self._open(layer) as f:
+            chunks = f[FEATURES_DATASET].chunks
+        if chunks is None or axis >= len(chunks):
+            return None
+        return int(chunks[axis])
+
+    def read(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]] = None,
+        feature_slice: FeatureSlice = None,
+    ) -> np.ndarray:
+        indexers = _normalize_feature_slice(feature_slice)
+
+        with self._open(layer) as f:
+            dset = f[FEATURES_DATASET]
+
+            if labels is None:
+                return self._read_rows(dset, slice(None), indexers)
+
+            rows = self._row_indices(labels)
+            # h5py wants a strictly increasing index list and allows at most one
+            # fancy index per selection. Read each distinct row once in order,
+            # then restore the caller's order (and any repeats) with `inverse`.
+            uniq, inverse = np.unique(rows, return_inverse=True)
+            if uniq.size == len(rows) and np.array_equal(uniq, rows):
+                return self._read_rows(dset, list(uniq), indexers)
+            block = self._read_rows(dset, list(uniq), indexers)
+            return block[inverse]
+
+    @staticmethod
+    def _read_rows(
+        dset: h5py.Dataset, rows: Union[slice, List[int]], indexers: Tuple
+    ) -> np.ndarray:
+        """Read `rows` from `dset`, pushing as much of `indexers` into h5py as it allows."""
+        if not indexers:
+            return np.asarray(dset[rows])
+        if any(_is_fancy(ix) for ix in indexers):
+            # A fancy row list plus a fancy feature index would be two fancy
+            # indices in one selection, which h5py rejects. Read the feature
+            # axes whole and apply the fancy part in NumPy.
+            data = np.asarray(dset[rows])
+            return data[(slice(None), *indexers)]
+        # Plain slices are not fancy, so they ride along in the same selection
+        # and h5py reads only the chunks they cover.
+        return np.asarray(dset[(rows, *indexers)])
+
+    def _row_indices(self, labels: Sequence[str]) -> np.ndarray:
+        _ = self.labels  # populate the label -> row index map
+        assert self._label_index is not None
+        try:
+            return np.array([self._label_index[label] for label in labels], dtype=int)
+        except KeyError as exc:
+            raise KeyError(
+                "Label {} not found in {}".format(exc.args[0], self._dpath)
+            ) from None
+
+    def _open(self, layer: str) -> h5py.File:
+        return h5py.File(self.path(layer), "r")
+
+    def _validate_all(self) -> None:
+        for layer in self._layers:
+            with self._open(layer) as f:
+                _validate_format(f, self.path(layer))
+
+
+def _validate_format(f: h5py.File, path: str) -> None:
+    """Check that `f` is bdpy feature storage of a version we can read."""
+    fmt = f.attrs.get(FORMAT_ATTR)
+    if isinstance(fmt, bytes):
+        fmt = fmt.decode("utf-8")
+    if fmt != FORMAT_NAME:
+        raise RuntimeError(
+            "{} is not bdpy feature storage ({} = {!r}); expected {!r}".format(
+                path, FORMAT_ATTR, fmt, FORMAT_NAME
+            )
+        )
+    version = int(f.attrs.get(FORMAT_VERSION_ATTR, 0))
+    if version > SUPPORTED_FORMAT_VERSION:
+        raise RuntimeError(
+            "{} uses feature storage format version {}, but this version of "
+            "bdpy supports up to version {}. Please upgrade bdpy.".format(
+                path, version, SUPPORTED_FORMAT_VERSION
+            )
+        )
+    for name in (FEATURES_DATASET, LABELS_DATASET):
+        if name not in f:
+            raise RuntimeError("{} has no /{} dataset".format(path, name))
+
+
+def _decode_labels(raw: np.ndarray) -> List[str]:
+    """Decode an HDF5 string dataset into a list of str (h5py 3.x gives bytes)."""
+    return [
+        item.decode("utf-8") if isinstance(item, bytes) else str(item)
+        for item in np.asarray(raw).ravel().tolist()
+    ]
+
+
+def detect_format(dpath: str, ext: str = "mat") -> str:
+    """Guess which storage format `dpath` holds.
+
+    Returns ``'hdf5'`` when the directory holds ``<layer>.h5`` files and no
+    layer subdirectories, and ``'mat'`` otherwise. A directory holding both is
+    read as ``'mat'``, the historical layout; pass an explicit format to
+    override.
+
+    Parameters
+    ----------
+    dpath : str
+        Feature directory.
+    ext : str, optional
+        Extension of the legacy per-stimulus files (default: ``'mat'``).
+
+    Returns
+    -------
+    str
+        ``'hdf5'`` or ``'mat'``.
+    """
+    escaped = glob.escape(dpath)
+    has_subdirs = any(
+        os.path.isdir(os.path.join(dpath, d)) for d in os.listdir(dpath)
+    )
+    has_h5 = bool(glob.glob(os.path.join(escaped, "*." + HDF5_EXT)))
+    if has_h5 and not has_subdirs:
+        return "hdf5"
+    return "mat"
diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py
index 1a181f5..44fc63a 100644
--- a/bdpy/dataform/features.py
+++ b/bdpy/dataform/features.py
@@ -15,12 +15,19 @@
 import warnings
 from functools import partial
 from multiprocessing import Pool
-from typing import Any, Dict, List, Optional, Union
+from typing import Any, Dict, Iterator, List, Optional, Sequence, Tuple, Union
 
 import hdf5storage
 import numpy as np
 
 from . import _mat_v73
+from ._feature_store import (
+    FeatureSlice,
+    FeatureStore,
+    HDF5FeatureStore,
+    MatFeatureStore,
+    detect_format,
+)
 
 # Deprecation notice emitted by the MATLAB-compatible write path. The actual
 # switch to bdpy-native plain HDF5 (and the drop of the hdf5storage write
@@ -58,12 +65,29 @@ def _determine_num_parallel(num_files: int) -> int:
 class Features(object):
     """DNN features class.
 
+    Reads a collection of DNN features from disk. Two on-disk layouts are
+    supported and behave identically through this API:
+
+    - the legacy per-stimulus layout, ``<dpath>/<layer>/<label>.mat``;
+    - chunked HDF5 storage, ``<dpath>/<layer>.h5`` (see
+      :mod:`bdpy.dataform.feature_hdf5`), which additionally supports partial
+      reads along the feature axes without loading a whole layer.
+
     Parameters
     ----------
     dpath: str or list
         (List of) DNN feature directory(ies)
     ext: str
-        DNN feature file extension (default: mat)
+        DNN feature file extension of the legacy per-stimulus layout
+        (default: mat). Ignored by chunked HDF5 storage.
+    feature_index: str, optional
+        Path to a ``.mat`` file holding a per-layer unit index under the key
+        ``index``. When given, features are flattened and restricted to those
+        units.
+    format: str
+        Storage format: ``'auto'`` (default), ``'mat'`` or ``'hdf5'``. With
+        ``'auto'`` each directory is inspected independently, so a mix of
+        layouts across `dpath` entries is fine.
 
     Attributes
     ----------
@@ -77,17 +101,20 @@ class Features(object):
 
     def __init__(
             self, dpath: Union[str, List[str]] = [],
-            ext: str = 'mat', feature_index: Optional[str] = None
+            ext: str = 'mat', feature_index: Optional[str] = None,
+            format: str = 'auto'
         ):
         if not isinstance(dpath, list):
             dpath = [dpath]
         self.__dpath = dpath
+        self.__format = format
 
-        self.__feature_file_table: Dict[str, Dict[str, str]] = {}  # Stimulus feature file tables
+        self.__stores: List[FeatureStore] = []
+        self.__label_store: Dict[str, FeatureStore] = {}  # label -> owning store
         self.__labels: List[str] = []  # Stimulus labels
         self.__index: List[int] = []  # Stimulus index (one-based)
         self.__layers: List[str] = []  # DNN layers
-        self.__collect_feature_files(ext=ext)
+        self.__collect_features(ext=ext)
 
         self.__c_feature_name: Optional[str] = None  # Loaded layer
         self.__features: Optional[np.ndarray] = None  # Loaded features
@@ -126,7 +153,26 @@ def layers(self):
     def feature_index(self):
         return self.__feature_index
 
-    def get(self, layer: str, label: Union[str, List[str], None] = None) -> np.ndarray:
+    def shape(self, layer: str) -> Tuple[int, ...]:
+        """Return the full shape of `layer` without reading the features.
+
+        Parameters
+        ----------
+        layer: str
+            DNN layer
+
+        Returns
+        -------
+        tuple of int
+            ``(n_samples, *feature_shape)``
+        """
+        feature_shape = self.__stores[0].shape(layer)[1:]
+        return (len(self.__labels), *feature_shape)
+
+    def get(
+            self, layer: str, label: Union[str, List[str], None] = None,
+            feature_slice: FeatureSlice = None
+        ) -> np.ndarray:
         """Return features in `layer`.
 
         Parameters
@@ -134,40 +180,110 @@ def get(self, layer: str, label: Union[str, List[str], None] = None) -> np.ndarr
         layer: str
             DNN layer
         label: str or list
-            Sample label(s)
+            Sample label(s). Rows come back in the order given.
+        feature_slice: slice, int, array-like or tuple, optional
+            Index applied to the feature axes (axes 1 and up), as produced by
+            ``numpy.s_[...]``. With chunked HDF5 storage this is a genuine
+            partial read; with the legacy layout the files are loaded in full
+            and then sliced.
 
         Returns
         -------
         numpy.ndarray, shape=(n_samples, shape_layers)
             DNN features
+
+        Examples
+        --------
+        >>> features.get('conv5')                                  # doctest: +SKIP
+        >>> features.get('conv5', label=['img0001', 'img0002'])     # doctest: +SKIP
+        >>> features.get('conv5', feature_slice=np.s_[128:256])     # doctest: +SKIP
         """
-        if label is None:
+        if label is None and feature_slice is None:
             return self.get_features(layer)
 
-        if isinstance(label, str):
+        if label is None:
+            labels = self.__labels
+        elif isinstance(label, str):
             labels = [label]
         else:
-            labels = label
+            labels = list(label)
 
-        features: np.ndarray
-        num_labels = len(labels)
-        num_parallel = _determine_num_parallel(num_labels)
-        path_iterator = map(lambda label: self.__feature_file_table[layer][label], labels)
-        if num_parallel == 1:
-            features = np.concatenate(list(map(partial(_load_array_with_key, 'feat'), path_iterator)), axis=0)
-        else:
-            with Pool(processes=num_parallel) as pool:
-                features = np.concatenate(pool.map(partial(_load_array_with_key, 'feat'), path_iterator), axis=0)
+        features = self.__read(layer, labels, feature_slice)
+        return self.__apply_feature_index(features, layer)
 
+    def iter_chunks(
+            self, layer: str, label: Union[str, List[str], None] = None,
+            feature_slice: FeatureSlice = None, axis: int = 1,
+            size: Optional[int] = None
+        ) -> Iterator[Tuple[slice, np.ndarray]]:
+        """Iterate over `layer` in slabs along `axis`.
+
+        This is the streaming counterpart of :meth:`get`: it never holds more
+        than one slab in memory, so a layer far larger than RAM can be processed
+        end to end.
+
+        Parameters
+        ----------
+        layer: str
+            DNN layer
+        label: str or list, optional
+            Sample label(s). ``None`` iterates over every label.
+        feature_slice: slice, int, array-like or tuple, optional
+            Index applied to the feature axes before iterating.
+        axis: int
+            Axis of the selected array to iterate over. Axis 0 is the sample
+            axis; the default, axis 1, is the outermost feature axis.
+        size: int, optional
+            Elements per slab. Defaults to the on-disk chunk extent along
+            `axis`, so that each element is read exactly once.
+
+        Yields
+        ------
+        tuple of (slice, numpy.ndarray)
+            The slice applied to `axis`, and the corresponding slab.
+
+        Raises
+        ------
+        RuntimeError
+            If a unit index (`feature_index`) is in use, which flattens the
+            feature axes and so has no meaningful per-axis iteration.
+
+        Examples
+        --------
+        >>> for sl, block in features.iter_chunks('conv5'):  # doctest: +SKIP
+        ...     out[:, sl] = transform(block)
+        """
         if self.__feat_index_table is not None:
-            # Select features by index
-            self.__feature_index = self.__feat_index_table[layer]
-            n_sample = features.shape[0]
-            n_feat = np.array(features.shape[1:]).prod()
+            raise RuntimeError(
+                'iter_chunks is not supported together with feature_index, '
+                'which flattens the feature axes'
+            )
 
-            features = features.reshape([n_sample, n_feat], order='C')[:, self.__feature_index]
+        if label is None:
+            labels = self.__labels
+        elif isinstance(label, str):
+            labels = [label]
+        else:
+            labels = list(label)
 
-        return features
+        store = self.__store_for(labels)
+        if store is not None:
+            yield from store.iter_chunks(
+                layer, labels, feature_slice, axis=axis, size=size
+            )
+            return
+
+        # Labels span several directories, so no single store can stream them.
+        # Fall back to slicing a full read, which still yields the same blocks.
+        features = self.__read(layer, labels, feature_slice)
+        length = features.shape[axis]
+        if size is None:
+            size = length
+        for start in range(0, length, size):
+            sl = slice(start, min(start + size, length))
+            index: List[Any] = [slice(None)] * features.ndim
+            index[axis] = sl
+            yield sl, features[tuple(index)]
 
     def statistic(self, statistic: str = 'mean', layer: Optional[str] = None):
         # NOTE: return type is ambiguous. currently, it is inferred as Unkown | Any
@@ -220,79 +336,84 @@ def get_features(self, layer: str) -> np.ndarray:
             assert isinstance(self.__features, np.ndarray)
             return self.__features  # self.__features could be None
 
-        features: np.ndarray
-        num_labels = len(self.__labels)
-        num_parallel = _determine_num_parallel(num_labels)
-        path_iterator = map(lambda label: self.__feature_file_table[layer][label], self.__labels)
-        if num_parallel == 1:
-            features = np.concatenate(list(map(partial(_load_array_with_key, 'feat'), path_iterator)), axis=0)
-        else:
-            with Pool(processes=num_parallel) as pool:
-                features = np.concatenate(pool.map(partial(_load_array_with_key, 'feat'), path_iterator), axis=0)
-        self.__features = features
-
+        self.__features = self.__read(layer, self.__labels)
         self.__c_feature_name = layer
-
-        if self.__feat_index_table is not None:
-            # Select features by index
-            self.__feature_index = self.__feat_index_table[layer]
-            n_sample = self.__features.shape[0]
-            n_feat = np.array(self.__features.shape[1:]).prod()
-
-            self.__features = self.__features.reshape([n_sample, n_feat], order='C')[:, self.__feature_index]
+        self.__features = self.__apply_feature_index(self.__features, layer)
 
         return self.__features
 
-    def __collect_feature_files(self, ext='mat'):
-        dpath_lst = self.__dpath
-
-        # List-up layers and stimulus labels
-        label_dir = {}
-        for dpath in dpath_lst:
-            # List-up layers
-            self.__layers = self.__get_layers(dpath)
+    def __read(
+            self, layer: str, labels: Sequence[str],
+            feature_slice: FeatureSlice = None
+        ) -> np.ndarray:
+        """Read `labels` from `layer`, dispatching to the owning store(s)."""
+        store = self.__store_for(labels)
+        if store is not None:
+            return store.read(layer, labels, feature_slice)
+        # Labels are spread over several directories: read each run of
+        # consecutive labels from its own store, preserving the caller's order.
+        blocks: List[np.ndarray] = []
+        run: List[str] = []
+        run_store: Optional[FeatureStore] = None
+        for label in labels:
+            owner = self.__label_store[label]
+            if owner is not run_store and run:
+                assert run_store is not None
+                blocks.append(run_store.read(layer, run, feature_slice))
+                run = []
+            run_store = owner
+            run.append(label)
+        if run:
+            assert run_store is not None
+            blocks.append(run_store.read(layer, run, feature_slice))
+        return np.concatenate(blocks, axis=0)
+
+    def __store_for(self, labels: Sequence[str]) -> Optional[FeatureStore]:
+        """The single store owning every label, or None if they are spread out."""
+        if len(self.__stores) == 1:
+            return self.__stores[0]
+        owners = {id(self.__label_store[label]) for label in labels}
+        if len(owners) == 1:
+            return self.__label_store[labels[0]]
+        return None
 
-            # List-up stimulus labels
-            labels_in_dir = self.__get_labels(dpath, self.__layers, ext=ext)
-            label_dir.update({label: dpath for label in labels_in_dir})
-            self.__labels += labels_in_dir
+    def __apply_feature_index(self, features: np.ndarray, layer: str) -> np.ndarray:
+        if self.__feat_index_table is None:
+            return features
+        # Select features by index
+        self.__feature_index = self.__feat_index_table[layer]
+        n_sample = features.shape[0]
+        n_feat = np.array(features.shape[1:]).prod()
+        return features.reshape([n_sample, n_feat], order='C')[:, self.__feature_index]
+
+    def __collect_features(self, ext='mat'):
+        for dpath in self.__dpath:
+            store = self.__make_store(dpath, ext)
+            if self.__layers and store.layers != self.__layers:
+                raise RuntimeError('Invalid layers in %s' % dpath)
+            self.__layers = store.layers
+            self.__stores.append(store)
+            for label in store.labels:
+                self.__label_store[label] = store
+            self.__labels += store.labels
 
         # NOTE: type incompatibility here. Is it OK to cast to list?
         self.__index = np.arange(len(self.__labels)) + 1
 
-        # List-up feature files
-        for lay in self.__layers:
-            self.__feature_file_table.update(
-                {
-                    lay:
-                    {
-                        label:
-                         os.path.join(label_dir[label], lay, label + '.' + ext)
-                         for label in self.__labels
-                    }
-                })
-
         return None
 
-    def __get_layers(self, dpath: str):
-        layers = sorted([d for d in os.listdir(dpath) if os.path.isdir(os.path.join(dpath, d))])
-        if self.__layers and (layers != self.__layers):
-            raise RuntimeError('Invalid layers in %s' % dpath)
-        return layers
-
-    def __get_labels(self, dpath: str, layers: List[str], ext: str = 'mat'):
-        labels: List[str] = []
-        for lay in layers:
-            lay_dir = os.path.join(dpath, lay)
-            lay_dir = lay_dir.replace('[', '[[]') # Use glob.escape for Python 3.4 or later
-            files = glob.glob(os.path.join(lay_dir, '*.' + ext))
-            labels_t = sorted([os.path.splitext(os.path.basename(f))[0] for f in files])
-            if not labels:
-                labels = labels_t
-            else:
-                if labels != labels_t:
-                    raise RuntimeError('Invalid feature file in %s ' % dpath)
-        return labels
+    def __make_store(self, dpath: str, ext: str) -> FeatureStore:
+        fmt = self.__format
+        if fmt == 'auto':
+            fmt = detect_format(dpath, ext=ext)
+        if fmt == 'hdf5':
+            return HDF5FeatureStore(dpath)
+        if fmt == 'mat':
+            return MatFeatureStore(dpath, ext=ext)
+        raise ValueError(
+            "Unknown feature storage format {!r}; expected "
+            "'auto', 'mat' or 'hdf5'".format(fmt)
+        )
 
 
 class DecodedFeatures(object):

From 6c43af73e2ca9ce2cf65c90f638e4fa66896339b Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Fri, 18 Sep 2026 05:02:04 +0000
Subject: [PATCH 03/11] feat(dataform): add HDF5 feature writer, converter and
 tests

Reading the chunked format is no use without a way to produce it, so add
the write side and the migration path:

- save_features() writes a whole layer at once.
- FeatureWriter appends incrementally onto resizable datasets, which is
  what feature extraction needs since it yields one stimulus at a time.
- convert_features_to_hdf5() migrates an existing .mat tree. It reads
  through MatFeatureStore, so the output is by construction what the
  legacy reader sees, and streams in batches so a layer is never held in
  memory in full.

All three are exported from bdpy.dataform. save_feature() is left alone;
it belongs to the separate MATLAB-write deprecation.

Tests cover the part that matters most: every combination of label
selection and feature slicing is run against both backends and compared,
and against ground truth, so the two layouts are indistinguishable
through the Features API -- including out-of-order labels, repeated
labels, fancy feature indices and mixed-layout dpaths. Partial reads are
asserted by spying on the selection handed to h5py, since byte counts are
not portable.

Note on feature_index: its happy path is deliberately left untested. Such
index files are struct .mat files, which _mat_v73.loadmat_key cannot
load (dense arrays only), so the feature has been broken since the
hdf5storage -> h5py read-path change of issue #106. That regression is
independent of storage backends and is left untouched here rather than
silently changed; the test file records it.
---
 bdpy/dataform/__init__.py            |   1 +
 bdpy/dataform/_feature_store.py      |  12 +-
 bdpy/dataform/feature_hdf5.py        | 409 +++++++++++++++++++++++++++
 docs/dataform_features.md            | 124 ++++++++
 tests/dataform/test_feature_hdf5.py  | 246 ++++++++++++++++
 tests/dataform/test_feature_store.py | 273 ++++++++++++++++++
 tests/dataform/test_features.py      |  96 ++++++-
 7 files changed, 1155 insertions(+), 6 deletions(-)
 create mode 100644 bdpy/dataform/feature_hdf5.py
 create mode 100644 tests/dataform/test_feature_hdf5.py
 create mode 100644 tests/dataform/test_feature_store.py

diff --git a/bdpy/dataform/__init__.py b/bdpy/dataform/__init__.py
index 30d66a6..fe9ea43 100644
--- a/bdpy/dataform/__init__.py
+++ b/bdpy/dataform/__init__.py
@@ -5,6 +5,7 @@
 """
 
 from .datastore import *
+from .feature_hdf5 import FeatureWriter, convert_features_to_hdf5, save_features
 from .features import *
 from .kvs import SQLite3KeyValueStore
 from .pd import *
diff --git a/bdpy/dataform/_feature_store.py b/bdpy/dataform/_feature_store.py
index 38e2112..1db4fa4 100644
--- a/bdpy/dataform/_feature_store.py
+++ b/bdpy/dataform/_feature_store.py
@@ -473,13 +473,19 @@ def read(
                 return self._read_rows(dset, slice(None), indexers)
 
             rows = self._row_indices(labels)
+            if rows.size == dset.shape[0] and np.array_equal(
+                rows, np.arange(dset.shape[0])
+            ):
+                # Asking for every row in order: a plain slice reads
+                # contiguously, where a full index list would not.
+                return self._read_rows(dset, slice(None), indexers)
             # h5py wants a strictly increasing index list and allows at most one
             # fancy index per selection. Read each distinct row once in order,
             # then restore the caller's order (and any repeats) with `inverse`.
             uniq, inverse = np.unique(rows, return_inverse=True)
-            if uniq.size == len(rows) and np.array_equal(uniq, rows):
-                return self._read_rows(dset, list(uniq), indexers)
-            block = self._read_rows(dset, list(uniq), indexers)
+            block = self._read_rows(dset, [int(i) for i in uniq], indexers)
+            if uniq.size == rows.size and np.array_equal(uniq, rows):
+                return block
             return block[inverse]
 
     @staticmethod
diff --git a/bdpy/dataform/feature_hdf5.py b/bdpy/dataform/feature_hdf5.py
new file mode 100644
index 0000000..8a5b3a4
--- /dev/null
+++ b/bdpy/dataform/feature_hdf5.py
@@ -0,0 +1,409 @@
+"""Chunked HDF5 storage for DNN features.
+
+The legacy feature layout writes one ``.mat`` file per stimulus, so the only
+file boundary is the sample axis and a partial read along the feature axes is
+impossible. This module writes the chunked layout instead: one ``<layer>.h5``
+per layer, holding the whole layer in a dataset that is explicitly chunked on
+both the sample axis and the outermost feature axis.
+
+Schema version 1::
+
+    <layer>.h5
+        /features   (n_samples, *feature_shape)
+        /labels     (n_samples,)   variable-length UTF-8
+        attrs:
+            bdpy_format         = "features"
+            bdpy_format_version = 1
+            layer               = "<layer name>"
+
+One file per layer rather than one file for everything: layers differ in shape
+and in optimal chunk shape, and separate files keep regeneration, copying and
+parallel writing per-layer.
+
+Compression is off by default. The point of this format is partial-read latency,
+and every compressed chunk costs a decompression on the way out; pass
+``compression=`` when size matters more than speed.
+
+This file is a part of BdPy.
+
+Examples
+--------
+Write a whole layer at once::
+
+    save_features('features/conv5.h5', array, labels)
+
+Write incrementally, as feature extraction produces one stimulus at a time::
+
+    with FeatureWriter('features/conv5.h5', feature_shape=(256, 13, 13),
+                       dtype=np.float32) as writer:
+        for label, feature in extract():
+            writer.append(feature, label)
+
+Migrate an existing ``.mat`` tree::
+
+    convert_features_to_hdf5('features_mat', 'features_h5')
+"""
+
+import os
+from types import TracebackType
+from typing import Iterable, Optional, Sequence, Tuple, Type
+
+import h5py
+import numpy as np
+
+from ._feature_chunking import DEFAULT_TARGET_CHUNK_BYTES, choose_chunk_shape
+from ._feature_store import (
+    FEATURES_DATASET,
+    FORMAT_ATTR,
+    FORMAT_NAME,
+    FORMAT_VERSION_ATTR,
+    HDF5_EXT,
+    LABELS_DATASET,
+    SUPPORTED_FORMAT_VERSION,
+    MatFeatureStore,
+)
+
+__all__ = [
+    "FeatureWriter",
+    "convert_features_to_hdf5",
+    "save_features",
+]
+
+#: Number of samples the resizable dataset grows by at a time.
+_GROW_BLOCK = 64
+
+
+def _string_dtype() -> np.dtype:
+    dtype: np.dtype = h5py.string_dtype(encoding="utf-8")
+    return dtype
+
+
+def _write_header(f: h5py.File, layer: Optional[str]) -> None:
+    f.attrs[FORMAT_ATTR] = FORMAT_NAME
+    f.attrs[FORMAT_VERSION_ATTR] = SUPPORTED_FORMAT_VERSION
+    if layer is not None:
+        f.attrs["layer"] = layer
+
+
+def save_features(
+    path: str,
+    features: np.ndarray,
+    labels: Sequence[str],
+    layer: Optional[str] = None,
+    chunks: Optional[Tuple[int, ...]] = None,
+    target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
+    compression: Optional[str] = None,
+    dtype: Optional[np.dtype] = None,
+) -> None:
+    """Write a whole layer to a chunked HDF5 feature file.
+
+    Parameters
+    ----------
+    path : str
+        Output file. Overwritten if it exists.
+    features : numpy.ndarray
+        Feature array of shape ``(n_samples, *feature_shape)``.
+    labels : sequence of str
+        One stimulus label per sample, in the same order as `features`.
+    layer : str, optional
+        Layer name recorded in the file. Defaults to the file's basename.
+    chunks : tuple of int, optional
+        Explicit chunk shape. Defaults to :func:`choose_chunk_shape`.
+    target_chunk_bytes : int, optional
+        Per-chunk byte budget used when `chunks` is not given.
+    compression : str, optional
+        h5py compression filter, e.g. ``'lzf'`` or ``'gzip'``. ``None``
+        (default) stores the features uncompressed, which keeps partial reads
+        fast.
+    dtype : numpy.dtype, optional
+        Dtype to store. Defaults to the dtype of `features`.
+
+    Raises
+    ------
+    ValueError
+        If `labels` does not have one entry per sample, or `features` has fewer
+        than two axes.
+    """
+    features = np.asarray(features)
+    if features.ndim < 2:
+        raise ValueError(
+            "features need a sample axis and at least one feature axis; "
+            "got shape {}".format(features.shape)
+        )
+    labels = list(labels)
+    if len(labels) != features.shape[0]:
+        raise ValueError(
+            "got {} labels for {} samples".format(len(labels), features.shape[0])
+        )
+
+    dtype = np.dtype(features.dtype if dtype is None else dtype)
+    if chunks is None:
+        chunks = choose_chunk_shape(
+            features.shape, dtype, target_bytes=target_chunk_bytes
+        )
+    if layer is None:
+        layer = os.path.splitext(os.path.basename(path))[0]
+
+    _makedirs_for(path)
+    with h5py.File(path, "w") as f:
+        _write_header(f, layer)
+        f.create_dataset(
+            FEATURES_DATASET,
+            data=features.astype(dtype, copy=False),
+            chunks=chunks,
+            compression=compression,
+        )
+        f.create_dataset(LABELS_DATASET, data=labels, dtype=_string_dtype())
+
+
+class FeatureWriter:
+    """Incremental writer for a chunked HDF5 feature file.
+
+    Feature extraction produces one stimulus at a time, so the datasets are
+    created resizable (``maxshape=(None, *feature_shape)``) and grown in blocks
+    as samples arrive. Use it as a context manager, or call :meth:`close`.
+
+    Parameters
+    ----------
+    path : str
+        Output file. Overwritten if it exists.
+    feature_shape : sequence of int
+        Shape of a single sample's features, without the sample axis.
+    dtype : numpy.dtype
+        Dtype to store.
+    layer : str, optional
+        Layer name recorded in the file. Defaults to the file's basename.
+    n_samples : int, optional
+        Expected number of samples, if known. Used only to pick a chunk shape.
+    chunks : tuple of int, optional
+        Explicit chunk shape, overriding the byte budget.
+    target_chunk_bytes : int, optional
+        Per-chunk byte budget used when `chunks` is not given.
+    compression : str, optional
+        h5py compression filter. ``None`` (default) stores uncompressed.
+
+    Examples
+    --------
+    >>> with FeatureWriter('conv5.h5', (256, 13, 13), np.float32) as w:  # doctest: +SKIP
+    ...     w.append(feature, 'img0001')
+    """
+
+    def __init__(
+        self,
+        path: str,
+        feature_shape: Sequence[int],
+        dtype: np.dtype,
+        layer: Optional[str] = None,
+        n_samples: Optional[int] = None,
+        chunks: Optional[Tuple[int, ...]] = None,
+        target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
+        compression: Optional[str] = None,
+    ):
+        self._feature_shape = tuple(int(s) for s in feature_shape)
+        if not self._feature_shape:
+            raise ValueError("feature_shape must have at least one axis")
+        self._dtype = np.dtype(dtype)
+        self._n = 0
+
+        if chunks is None:
+            chunks = choose_chunk_shape(
+                (n_samples if n_samples is not None else 0, *self._feature_shape),
+                self._dtype,
+                target_bytes=target_chunk_bytes,
+                n_samples_known=n_samples is not None,
+            )
+        if layer is None:
+            layer = os.path.splitext(os.path.basename(path))[0]
+
+        _makedirs_for(path)
+        self._file: Optional[h5py.File] = h5py.File(path, "w")
+        _write_header(self._file, layer)
+        self._features = self._file.create_dataset(
+            FEATURES_DATASET,
+            shape=(0, *self._feature_shape),
+            maxshape=(None, *self._feature_shape),
+            dtype=self._dtype,
+            chunks=chunks,
+            compression=compression,
+        )
+        self._labels = self._file.create_dataset(
+            LABELS_DATASET,
+            shape=(0,),
+            maxshape=(None,),
+            dtype=_string_dtype(),
+            chunks=(max(1, _GROW_BLOCK),),
+        )
+
+    @property
+    def n_samples(self) -> int:
+        """Number of samples written so far."""
+        return self._n
+
+    def append(self, feature: np.ndarray, label: str) -> None:
+        """Append a single sample.
+
+        Parameters
+        ----------
+        feature : numpy.ndarray
+            One sample, either ``feature_shape`` or ``(1, *feature_shape)``.
+        label : str
+            Its stimulus label.
+        """
+        feature = np.asarray(feature)
+        if feature.shape == self._feature_shape:
+            feature = feature[np.newaxis]
+        self.extend(feature, [label])
+
+    def extend(self, features: np.ndarray, labels: Sequence[str]) -> None:
+        """Append several samples at once.
+
+        Parameters
+        ----------
+        features : numpy.ndarray
+            Samples of shape ``(n, *feature_shape)``.
+        labels : sequence of str
+            One label per sample.
+
+        Raises
+        ------
+        ValueError
+            If the feature shape does not match the writer's, or the label count
+            does not match the sample count.
+        RuntimeError
+            If the writer is already closed.
+        """
+        if self._file is None:
+            raise RuntimeError("FeatureWriter is closed")
+        features = np.asarray(features)
+        labels = list(labels)
+        if features.shape[1:] != self._feature_shape:
+            raise ValueError(
+                "expected features of shape (n, {}), got {}".format(
+                    ", ".join(str(s) for s in self._feature_shape), features.shape
+                )
+            )
+        if len(labels) != features.shape[0]:
+            raise ValueError(
+                "got {} labels for {} samples".format(len(labels), features.shape[0])
+            )
+        if not labels:
+            return
+
+        new_n = self._n + len(labels)
+        self._features.resize(new_n, axis=0)
+        self._labels.resize(new_n, axis=0)
+        self._features[self._n:new_n] = features.astype(self._dtype, copy=False)
+        self._labels[self._n:new_n] = labels
+        self._n = new_n
+
+    def close(self) -> None:
+        """Close the underlying file. Idempotent."""
+        if self._file is not None:
+            self._file.close()
+            self._file = None
+
+    def __enter__(self) -> "FeatureWriter":
+        """Enter the context manager."""
+        return self
+
+    def __exit__(
+        self,
+        exc_type: Optional[Type[BaseException]],
+        exc: Optional[BaseException],
+        tb: Optional[TracebackType],
+    ) -> None:
+        """Close the file on leaving the context manager."""
+        self.close()
+
+
+def convert_features_to_hdf5(
+    src_dir: str,
+    dst_dir: str,
+    layers: Optional[Iterable[str]] = None,
+    ext: str = "mat",
+    key: str = "feat",
+    overwrite: bool = False,
+    batch_size: int = 64,
+    target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
+    compression: Optional[str] = None,
+    verbose: bool = False,
+) -> None:
+    """Convert a legacy per-stimulus feature directory to chunked HDF5.
+
+    Reads `src_dir` through the same backend :class:`~bdpy.dataform.Features`
+    uses for the legacy layout, so the output is by construction what the legacy
+    reader sees. Samples are streamed in batches, so a layer is never held in
+    memory in full.
+
+    Parameters
+    ----------
+    src_dir : str
+        Legacy directory, ``<src_dir>/<layer>/<label>.<ext>``.
+    dst_dir : str
+        Output directory; ``<dst_dir>/<layer>.h5`` is written per layer. Created
+        if missing.
+    layers : iterable of str, optional
+        Layers to convert. Defaults to every layer found in `src_dir`.
+    ext : str, optional
+        Extension of the legacy files (default: ``'mat'``).
+    key : str, optional
+        Variable name inside the legacy files (default: ``'feat'``).
+    overwrite : bool, optional
+        Overwrite an existing output file instead of skipping it
+        (default: False).
+    batch_size : int, optional
+        Number of stimulus files read per batch.
+    target_chunk_bytes : int, optional
+        Per-chunk byte budget.
+    compression : str, optional
+        h5py compression filter. ``None`` (default) stores uncompressed.
+    verbose : bool, optional
+        Print progress.
+
+    Raises
+    ------
+    KeyError
+        If a requested layer is not present in `src_dir`.
+    """
+    store = MatFeatureStore(src_dir, ext=ext, key=key)
+    selected = list(store.layers) if layers is None else list(layers)
+    missing = [layer for layer in selected if layer not in store.layers]
+    if missing:
+        raise KeyError(
+            "Layer(s) {} not found in {}".format(", ".join(missing), src_dir)
+        )
+
+    os.makedirs(dst_dir, exist_ok=True)
+    all_labels = store.labels
+
+    for layer in selected:
+        out_path = os.path.join(dst_dir, layer + "." + HDF5_EXT)
+        if os.path.exists(out_path) and not overwrite:
+            if verbose:
+                print("{} already exists. Skipped.".format(out_path))
+            continue
+
+        full_shape = store.shape(layer)
+        writer = FeatureWriter(
+            out_path,
+            feature_shape=full_shape[1:],
+            dtype=store.dtype(layer),
+            layer=layer,
+            n_samples=full_shape[0],
+            target_chunk_bytes=target_chunk_bytes,
+            compression=compression,
+        )
+        try:
+            for start in range(0, len(all_labels), batch_size):
+                batch = all_labels[start:start + batch_size]
+                writer.extend(store.read(layer, batch), batch)
+        finally:
+            writer.close()
+
+        if verbose:
+            print("Saved {} ({} samples).".format(out_path, full_shape[0]))
+
+
+def _makedirs_for(path: str) -> None:
+    parent = os.path.dirname(os.path.abspath(path))
+    os.makedirs(parent, exist_ok=True)
diff --git a/docs/dataform_features.md b/docs/dataform_features.md
index cd07c19..f5d0131 100644
--- a/docs/dataform_features.md
+++ b/docs/dataform_features.md
@@ -51,4 +51,128 @@ decoded_features.statistic('mean', layer='fc8', subject='sub-01', roi='VC')
 
 ```
 
+## Chunked HDF5 feature storage
 
+The layout above stores one file per stimulus, so the only file boundary is the
+sample axis: reading a few channels of a layer still costs a full read of every
+stimulus. For large DNN features, `Features` can instead read **chunked HDF5**
+storage, one file per layer:
+
+```
+features/
+  conv1_1.h5    # /features (n_stimuli, *feature_shape) + /labels (n_stimuli,)
+  conv2_1.h5
+  fc8.h5
+```
+
+`/features` is explicitly chunked along both the sample axis and the outermost
+feature axis, so a slice reads only the chunks it covers.
+
+### Reading
+
+Nothing changes for existing code -- `Features` detects the layout per
+directory, and `.mat` trees keep working exactly as before:
+
+``` python
+features = Features('/path/to/features')   # either layout
+feat = features.get(layer='conv5')
+```
+
+Pass `format='mat'` or `format='hdf5'` to skip detection. A directory holding
+both layouts is read as `.mat` unless you say otherwise.
+
+### Partial reads
+
+`feature_slice` selects along the feature axes (axis 1 and up). On chunked HDF5
+this is a real partial read; on the legacy layout the files are loaded in full
+and then sliced, so the same code works on both.
+
+``` python
+import numpy as np
+
+# Channels 128-255 only, without loading the rest of the layer
+feat = features.get(layer='conv5', feature_slice=np.s_[128:256])
+
+# Combine with label selection; rows come back in the order given
+feat = features.get(
+    layer='conv5',
+    label=['stimulus-0003', 'stimulus-0001'],
+    feature_slice=np.s_[128:256],
+)
+
+# Full shape without reading anything
+n_stimuli, *feature_shape = features.shape('conv5')
+```
+
+To stream a layer that does not fit in memory, iterate. `iter_chunks` yields
+`(slice, block)` so results can be placed back without tracking offsets, and by
+default uses the on-disk chunk extent, so every element is read exactly once:
+
+``` python
+out = np.empty(features.shape('conv5'))
+for sl, block in features.iter_chunks('conv5', axis=1):
+    out[:, sl] = transform(block)
+```
+
+`axis=0` iterates over stimuli instead of features.
+
+### Writing
+
+Write a whole layer at once:
+
+``` python
+from bdpy.dataform import save_features
+
+save_features('features/conv5.h5', array, labels)
+```
+
+Or incrementally, which is what feature extraction needs since it produces one
+stimulus at a time:
+
+``` python
+from bdpy.dataform import FeatureWriter
+
+with FeatureWriter('features/conv5.h5', feature_shape=(256, 13, 13),
+                   dtype=np.float32) as writer:
+    for label, feature in extract():
+        writer.append(feature, label)
+```
+
+### Migrating existing features
+
+``` python
+from bdpy.dataform import convert_features_to_hdf5
+
+convert_features_to_hdf5('/path/to/features_mat', '/path/to/features_h5')
+```
+
+The converter reads through the same backend `Features` uses for the legacy
+layout, and streams in batches, so a layer is never held in memory in full.
+
+### Chunk shape
+
+Chunk shape, not the file format alone, is what decides the cost of a partial
+read: HDF5 reads whole chunks even when only a few of their elements are
+selected. By default bdpy derives it from a 1 MiB budget, splitting it between
+the sample axis and the outermost feature axis and keeping trailing spatial axes
+whole -- e.g. `(1200, 256, 13, 13)` float32 becomes chunks of
+`(39, 39, 13, 13)`.
+
+That means reading a *single* stimulus costs one chunk row (~39 stimuli of I/O),
+which matters for `FeaturesDataset`-style per-sample access. Tune it if your
+access pattern is skewed:
+
+``` python
+save_features(path, array, labels, target_chunk_bytes=256 * 1024)  # smaller chunks
+save_features(path, array, labels, chunks=(1, 256, 13, 13))        # per-stimulus
+```
+
+Compression is off by default, since every compressed chunk costs a
+decompression on the way out. Pass `compression='gzip'` or `compression='lzf'`
+when size matters more than read speed.
+
+### Format
+
+Files are marked with root attributes `bdpy_format="features"` and
+`bdpy_format_version=1`, and are rejected on read if those are missing or
+newer than the running bdpy understands.
diff --git a/tests/dataform/test_feature_hdf5.py b/tests/dataform/test_feature_hdf5.py
new file mode 100644
index 0000000..2b661cf
--- /dev/null
+++ b/tests/dataform/test_feature_hdf5.py
@@ -0,0 +1,246 @@
+import os
+import tempfile
+import unittest
+import warnings
+
+import h5py
+import numpy as np
+from numpy.testing import assert_array_equal
+
+from bdpy.dataform import Features
+from bdpy.dataform._feature_store import (
+    FORMAT_ATTR,
+    FORMAT_NAME,
+    FORMAT_VERSION_ATTR,
+    SUPPORTED_FORMAT_VERSION,
+    HDF5FeatureStore,
+)
+from bdpy.dataform.feature_hdf5 import (
+    FeatureWriter,
+    convert_features_to_hdf5,
+    save_features,
+)
+
+from .test_features import prepare_mat_features
+
+
+class TestSaveFeatures(unittest.TestCase):
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.labels = ['img%04d' % i for i in range(12)]
+        self.data = np.random.rand(12, 32, 5, 5).astype(np.float32)
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_writes_schema_v1(self):
+        save_features(self.path, self.data, self.labels)
+
+        with h5py.File(self.path, 'r') as f:
+            self.assertEqual(f.attrs[FORMAT_ATTR], FORMAT_NAME)
+            self.assertEqual(
+                int(f.attrs[FORMAT_VERSION_ATTR]), SUPPORTED_FORMAT_VERSION
+            )
+            self.assertEqual(f.attrs['layer'], 'conv5')
+            assert_array_equal(f['features'][()], self.data)
+            self.assertEqual(
+                [s.decode('utf-8') for s in f['labels'][()]], self.labels
+            )
+
+    def test_features_dataset_is_chunked_and_uncompressed(self):
+        # Chunking is the whole point; compression would tax every partial read.
+        save_features(self.path, self.data, self.labels)
+        with h5py.File(self.path, 'r') as f:
+            self.assertIsNotNone(f['features'].chunks)
+            self.assertIsNone(f['features'].compression)
+
+    def test_explicit_chunks_and_compression(self):
+        save_features(
+            self.path, self.data, self.labels,
+            chunks=(4, 8, 5, 5), compression='gzip',
+        )
+        with h5py.File(self.path, 'r') as f:
+            self.assertEqual(f['features'].chunks, (4, 8, 5, 5))
+            self.assertEqual(f['features'].compression, 'gzip')
+        # Compression must not change what comes back out.
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_dtype_is_preserved(self):
+        for dtype in (np.float32, np.float64, np.int32):
+            with self.subTest(dtype=dtype):
+                path = os.path.join(self.tmpdir.name, 'l_%s.h5' % np.dtype(dtype).name)
+                save_features(path, self.data.astype(dtype), self.labels)
+                with h5py.File(path, 'r') as f:
+                    self.assertEqual(f['features'].dtype, np.dtype(dtype))
+
+    def test_rejects_mismatched_labels(self):
+        with self.assertRaises(ValueError):
+            save_features(self.path, self.data, self.labels[:-1])
+
+    def test_rejects_1d_features(self):
+        with self.assertRaises(ValueError):
+            save_features(self.path, np.arange(5), ['a'] * 5)
+
+
+class TestFeatureWriter(unittest.TestCase):
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+        self.labels = ['img%04d' % i for i in range(150)]
+        self.data = np.random.rand(150, 16, 3, 3).astype(np.float32)
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_append_round_trip(self):
+        # 150 samples crosses several resize blocks, which is where an
+        # incremental writer would lose data if the resize were wrong.
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            for label, feature in zip(self.labels, self.data):
+                writer.append(feature, label)
+            self.assertEqual(writer.n_samples, len(self.labels))
+
+        store = HDF5FeatureStore(self.tmpdir.name)
+        self.assertEqual(store.labels, self.labels)
+        assert_array_equal(store.read('conv5'), self.data)
+
+    def test_append_accepts_leading_sample_axis(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            writer.append(self.data[0][np.newaxis], self.labels[0])
+            writer.append(self.data[1], self.labels[1])
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data[:2]
+        )
+
+    def test_extend(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            writer.extend(self.data[:40], self.labels[:40])
+            writer.extend(self.data[40:], self.labels[40:])
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_writes_valid_header(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32, layer='conv5') as writer:
+            writer.append(self.data[0], self.labels[0])
+        with h5py.File(self.path, 'r') as f:
+            self.assertEqual(f.attrs[FORMAT_ATTR], FORMAT_NAME)
+            self.assertEqual(f.attrs['layer'], 'conv5')
+            self.assertIsNotNone(f['features'].chunks)
+
+    def test_rejects_wrong_shape(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            with self.assertRaises(ValueError):
+                writer.extend(np.zeros((2, 8, 3, 3)), ['a', 'b'])
+
+    def test_rejects_label_count_mismatch(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            with self.assertRaises(ValueError):
+                writer.extend(np.zeros((2, 16, 3, 3)), ['a'])
+
+    def test_write_after_close_raises(self):
+        writer = FeatureWriter(self.path, (16, 3, 3), np.float32)
+        writer.close()
+        writer.close()  # idempotent
+        with self.assertRaises(RuntimeError):
+            writer.append(self.data[0], self.labels[0])
+
+
+class TestConvertFeaturesToHDF5(unittest.TestCase):
+    """The converter must reproduce exactly what the legacy reader sees."""
+
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.matdir = os.path.join(self.tmpdir.name, 'mat')
+        self.h5dir = os.path.join(self.tmpdir.name, 'h5')
+        os.makedirs(self.matdir)
+        self.labels = ['img%04d' % i for i in range(10)]
+        self.stacked = prepare_mat_features(
+            self.matdir, ['conv5', 'fc8'], self.labels,
+            [(1, 16, 3, 3), (1, 50)],
+        )
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_converted_features_match_source(self):
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+
+        from_mat = Features(self.matdir)
+        from_h5 = Features(self.h5dir)
+        self.assertEqual(from_h5.layers, from_mat.layers)
+        self.assertEqual(from_h5.labels, from_mat.labels)
+        for layer in from_mat.layers:
+            with self.subTest(layer=layer):
+                assert_array_equal(from_h5.get(layer), self.stacked[layer])
+                assert_array_equal(from_h5.get(layer), from_mat.get(layer))
+
+    def test_selected_layers_only(self):
+        convert_features_to_hdf5(self.matdir, self.h5dir, layers=['fc8'])
+        self.assertEqual(Features(self.h5dir).layers, ['fc8'])
+
+    def test_unknown_layer_raises(self):
+        with self.assertRaises(KeyError):
+            convert_features_to_hdf5(self.matdir, self.h5dir, layers=['nope'])
+
+    def test_existing_file_is_skipped_unless_overwrite(self):
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+        marker = os.path.join(self.h5dir, 'fc8.h5')
+        mtime = os.path.getmtime(marker)
+        os.utime(marker, (mtime - 100, mtime - 100))
+
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+        self.assertEqual(os.path.getmtime(marker), mtime - 100)
+
+        convert_features_to_hdf5(self.matdir, self.h5dir, overwrite=True)
+        self.assertGreater(os.path.getmtime(marker), mtime - 100)
+
+    def test_small_batches_produce_the_same_file(self):
+        # Batching is an implementation detail; it must not affect the result.
+        convert_features_to_hdf5(self.matdir, self.h5dir, batch_size=3)
+        assert_array_equal(Features(self.h5dir).get('conv5'), self.stacked['conv5'])
+
+
+class TestFormatValidation(unittest.TestCase):
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_foreign_hdf5_is_rejected(self):
+        # A file that merely happens to have /features is not bdpy storage.
+        with h5py.File(self.path, 'w') as f:
+            f.create_dataset('features', data=np.zeros((2, 3)))
+            f.create_dataset('labels', data=['a', 'b'],
+                             dtype=h5py.string_dtype(encoding='utf-8'))
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+    def test_future_version_is_rejected_with_a_clear_message(self):
+        save_features(self.path, np.zeros((2, 3)), ['a', 'b'])
+        with h5py.File(self.path, 'a') as f:
+            f.attrs[FORMAT_VERSION_ATTR] = SUPPORTED_FORMAT_VERSION + 1
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('upgrade bdpy', str(ctx.exception))
+
+    def test_missing_dataset_is_rejected(self):
+        save_features(self.path, np.zeros((2, 3)), ['a', 'b'])
+        with h5py.File(self.path, 'a') as f:
+            del f['labels']
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+    def test_empty_directory_is_rejected(self):
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/dataform/test_feature_store.py b/tests/dataform/test_feature_store.py
new file mode 100644
index 0000000..ff7c8b1
--- /dev/null
+++ b/tests/dataform/test_feature_store.py
@@ -0,0 +1,273 @@
+"""Backend equivalence for feature storage.
+
+The legacy .mat layout and chunked HDF5 storage must be indistinguishable
+through the Features API: same labels, same order, same arrays, for every
+combination of label selection and feature slicing. These tests run each query
+against both backends and compare.
+"""
+
+import os
+import tempfile
+import unittest
+import warnings
+
+import h5py
+import numpy as np
+from numpy.testing import assert_array_equal
+
+from bdpy.dataform import Features, convert_features_to_hdf5, save_features
+from bdpy.dataform._feature_store import (
+    HDF5FeatureStore,
+    MatFeatureStore,
+    detect_format,
+)
+
+from .test_features import prepare_mat_features
+
+LAYERS = ['conv5', 'fc8']
+SHAPES = [(1, 24, 5, 5), (1, 60)]
+LABELS = ['img%04d' % i for i in range(11)]
+
+# Every way a caller might ask for features. Each is run against both backends.
+QUERIES = [
+    ('all', 'conv5', None, None),
+    ('single label', 'conv5', 'img0003', None),
+    ('label list', 'conv5', ['img0005', 'img0001', 'img0009'], None),
+    ('unsorted labels', 'conv5', ['img0009', 'img0000'], None),
+    ('repeated labels', 'conv5', ['img0009', 'img0001', 'img0009', 'img0000'], None),
+    ('feature slice', 'conv5', None, np.s_[8:16]),
+    ('slice with step', 'conv5', None, np.s_[::2]),
+    ('multi-axis slice', 'conv5', None, np.s_[8:16, 1:4, :]),
+    ('integer index', 'conv5', None, np.s_[5]),
+    ('fancy index', 'conv5', None, np.s_[[3, 1, 7]]),
+    ('ellipsis', 'conv5', None, np.s_[..., 1:3]),
+    ('labels and slice', 'conv5', ['img0009', 'img0001', 'img0009'], np.s_[8:16]),
+    ('2d layer', 'fc8', None, np.s_[10:40]),
+    ('2d layer with labels', 'fc8', ['img0002', 'img0000'], np.s_[10:40]),
+]
+
+
+class _BackendPair(unittest.TestCase):
+    """Builds the same features as a .mat tree and as chunked HDF5."""
+
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.matdir = os.path.join(self.tmpdir.name, 'mat')
+        self.h5dir = os.path.join(self.tmpdir.name, 'h5')
+        os.makedirs(self.matdir)
+        self.stacked = prepare_mat_features(self.matdir, LAYERS, LABELS, SHAPES)
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+        self.from_mat = Features(self.matdir)
+        self.from_h5 = Features(self.h5dir)
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+
+class TestBackendEquivalence(_BackendPair):
+    def test_metadata_matches(self):
+        self.assertEqual(self.from_h5.layers, self.from_mat.layers)
+        self.assertEqual(self.from_h5.labels, self.from_mat.labels)
+        assert_array_equal(self.from_h5.index, self.from_mat.index)
+        for layer in LAYERS:
+            self.assertEqual(self.from_h5.shape(layer), self.from_mat.shape(layer))
+
+    def test_queries_match_across_backends(self):
+        for name, layer, label, feature_slice in QUERIES:
+            with self.subTest(query=name):
+                from_mat = self.from_mat.get(
+                    layer, label=label, feature_slice=feature_slice
+                )
+                from_h5 = self.from_h5.get(
+                    layer, label=label, feature_slice=feature_slice
+                )
+                self.assertEqual(from_h5.shape, from_mat.shape)
+                assert_array_equal(from_h5, from_mat)
+
+    def test_queries_match_ground_truth(self):
+        # Both backends agreeing is not enough if both are wrong the same way.
+        for name, layer, label, feature_slice in QUERIES:
+            with self.subTest(query=name):
+                if label is None:
+                    rows = slice(None)
+                elif isinstance(label, str):
+                    rows = [LABELS.index(label)]
+                else:
+                    rows = [LABELS.index(s) for s in label]
+                expected = self.stacked[layer][rows]
+                if feature_slice is not None:
+                    indexers = (
+                        feature_slice
+                        if isinstance(feature_slice, tuple)
+                        else (feature_slice,)
+                    )
+                    expected = expected[(slice(None), *indexers)]
+                assert_array_equal(
+                    self.from_h5.get(layer, label=label, feature_slice=feature_slice),
+                    expected,
+                )
+
+    def test_statistic_matches(self):
+        for statistic in ('mean', 'std', 'std, ddof=0'):
+            for layer in LAYERS:
+                with self.subTest(statistic=statistic, layer=layer):
+                    assert_array_equal(
+                        self.from_h5.statistic(statistic, layer=layer),
+                        self.from_mat.statistic(statistic, layer=layer),
+                    )
+
+
+class TestIterChunks(_BackendPair):
+    def test_reassembles_to_a_full_read(self):
+        for axis in (0, 1):
+            for features in (self.from_mat, self.from_h5):
+                with self.subTest(axis=axis, backend=type(features).__name__):
+                    blocks = list(features.iter_chunks('conv5', axis=axis, size=4))
+                    joined = np.concatenate([b for _, b in blocks], axis=axis)
+                    assert_array_equal(joined, self.stacked['conv5'])
+
+    def test_yielded_slices_address_the_right_block(self):
+        full = self.from_h5.get('conv5')
+        for sl, block in self.from_h5.iter_chunks('conv5', axis=1, size=5):
+            assert_array_equal(block, full[:, sl])
+
+    def test_slices_cover_the_axis_without_overlap(self):
+        length = self.from_h5.shape('conv5')[1]
+        covered = []
+        for sl, _ in self.from_h5.iter_chunks('conv5', axis=1, size=5):
+            covered.extend(range(*sl.indices(length)))
+        self.assertEqual(covered, list(range(length)))
+
+    def test_composes_with_a_feature_slice(self):
+        expected = self.stacked['conv5'][:, 4:20]
+        blocks = list(
+            self.from_h5.iter_chunks('conv5', feature_slice=np.s_[4:20], axis=1, size=6)
+        )
+        assert_array_equal(np.concatenate([b for _, b in blocks], axis=1), expected)
+
+    def test_restricts_to_labels(self):
+        labels = ['img0009', 'img0001']
+        blocks = list(self.from_h5.iter_chunks('conv5', label=labels, axis=1, size=5))
+        joined = np.concatenate([b for _, b in blocks], axis=1)
+        assert_array_equal(joined, self.stacked['conv5'][[9, 1]])
+
+    def test_default_size_is_chunk_aligned(self):
+        store = HDF5FeatureStore(self.h5dir)
+        extent = store.chunk_extent('conv5', 1)
+        sizes = {
+            sl.stop - sl.start
+            for sl, _ in self.from_h5.iter_chunks('conv5', axis=1)
+        }
+        # Every block but the last is a whole chunk along the axis.
+        self.assertLessEqual(max(sizes), max(1, extent))
+
+    def test_rejects_bad_axis(self):
+        with self.assertRaises(ValueError):
+            list(self.from_h5.iter_chunks('conv5', axis=9))
+
+
+class TestPartialReads(_BackendPair):
+    def test_slice_is_pushed_down_to_h5py(self):
+        # The point of the format: h5py must receive the narrowed selection, so
+        # that only the covering chunks are read, rather than a full read that
+        # NumPy then slices.
+        seen = []
+        original = h5py.Dataset.__getitem__
+
+        def spy(dataset, args):
+            seen.append(args)
+            return original(dataset, args)
+
+        h5py.Dataset.__getitem__ = spy
+        try:
+            self.from_h5.get('conv5', feature_slice=np.s_[8:16])
+        finally:
+            h5py.Dataset.__getitem__ = original
+
+        self.assertEqual(len(seen), 1)
+        rows, feature_index = seen[0]
+        self.assertEqual(rows, slice(None))
+        self.assertEqual(feature_index, slice(8, 16, None))
+
+    def test_labels_are_read_once_and_in_order(self):
+        # Duplicated and out-of-order labels must reach h5py as a single
+        # strictly increasing index list, which is all h5py accepts.
+        seen = []
+        original = h5py.Dataset.__getitem__
+
+        def spy(dataset, args):
+            seen.append(args)
+            return original(dataset, args)
+
+        h5py.Dataset.__getitem__ = spy
+        try:
+            out = self.from_h5.get('conv5', label=['img0009', 'img0001', 'img0009'])
+        finally:
+            h5py.Dataset.__getitem__ = original
+
+        self.assertEqual(len(seen), 1)
+        rows = seen[0] if not isinstance(seen[0], tuple) else seen[0][0]
+        self.assertEqual(rows, [1, 9])
+        assert_array_equal(out, self.stacked['conv5'][[9, 1, 9]])
+
+    def test_dataset_is_chunked_on_both_sliceable_axes(self):
+        with h5py.File(os.path.join(self.h5dir, 'conv5.h5'), 'r') as f:
+            chunks = f['features'].chunks
+            shape = f['features'].shape
+        self.assertIsNotNone(chunks)
+        self.assertEqual(chunks[2:], shape[2:])  # spatial axes kept whole
+
+
+class TestFormatDetection(unittest.TestCase):
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.matdir = os.path.join(self.tmpdir.name, 'mat')
+        self.h5dir = os.path.join(self.tmpdir.name, 'h5')
+        os.makedirs(self.matdir)
+        self.stacked = prepare_mat_features(self.matdir, LAYERS, LABELS, SHAPES)
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_detects_each_layout(self):
+        self.assertEqual(detect_format(self.matdir), 'mat')
+        self.assertEqual(detect_format(self.h5dir), 'hdf5')
+
+    def test_explicit_format_overrides_detection(self):
+        self.assertIsInstance(
+            Features(self.h5dir, format='hdf5')._Features__stores[0],
+            HDF5FeatureStore,
+        )
+        self.assertIsInstance(
+            Features(self.matdir, format='mat')._Features__stores[0],
+            MatFeatureStore,
+        )
+
+    def test_unknown_format_raises(self):
+        with self.assertRaises(ValueError):
+            Features(self.h5dir, format='parquet')
+
+    def test_mixed_directories(self):
+        # One dpath per layout, read through a single Features.
+        h5_only = os.path.join(self.tmpdir.name, 'h5b')
+        os.makedirs(h5_only)
+        other_labels = ['other%04d' % i for i in range(4)]
+        for layer, shape in zip(LAYERS, SHAPES):
+            data = np.random.rand(len(other_labels), *shape[1:])
+            save_features(
+                os.path.join(h5_only, layer + '.h5'), data, other_labels, layer=layer
+            )
+
+        features = Features([self.matdir, h5_only])
+        self.assertEqual(features.labels, LABELS + other_labels)
+        # A query spanning both directories keeps the caller's order.
+        got = features.get('conv5', label=['other0002', 'img0003', 'other0000'])
+        assert_array_equal(got[1], self.stacked['conv5'][3])
+        self.assertEqual(got.shape[0], 3)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/dataform/test_features.py b/tests/dataform/test_features.py
index 8b8bbd4..421e6fe 100644
--- a/tests/dataform/test_features.py
+++ b/tests/dataform/test_features.py
@@ -13,13 +13,13 @@
 from bdpy.dataform.features import Features, save_feature
 
 
-def _prepare_mock_data(
+def prepare_mat_features(
         tmpdir: str,
         mock_layer_names: List[str],
         mock_image_names: List[str],
         mock_shapes: List[Tuple[int, ...]]
     ) -> dict:
-    """Prepare mock data for testing.
+    """Write a legacy per-stimulus feature directory for testing.
 
     Files are written as MATLAB v7.3 (HDF5) so the loader's h5py path (used for
     NumPy 2.0 compatibility) is exercised. The stacked arrays are returned so
@@ -57,7 +57,7 @@ def setUp(self):
         ]
         self.mock_shapes = [(1, 1000), (1, 256, 13, 13)]
         self.feature_dir = tempfile.TemporaryDirectory()
-        stacked = _prepare_mock_data(
+        stacked = prepare_mat_features(
             self.feature_dir.name,
             self.mock_layer_names,
             self.mock_image_names,
@@ -122,6 +122,96 @@ def test_features_get_label(self):
         )
 
 
+class TestFeaturesPartialRead(unittest.TestCase):
+    """feature_slice / iter_chunks on the legacy .mat backend.
+
+    The legacy layout cannot read partially, so these only check that the
+    result is the same as slicing a full read -- which is what makes the API
+    backend-independent.
+    """
+
+    def setUp(self):
+        self.layers = ['fc8', 'conv5']
+        self.labels = [
+            'n01443537_22563',
+            'n01443537_22564',
+            'n01677366_18182',
+            'n04572121_3262',
+        ]
+        self.shapes = [(1, 100), (1, 32, 5, 5)]
+        self.feature_dir = tempfile.TemporaryDirectory()
+        self.stacked = prepare_mat_features(
+            self.feature_dir.name, self.layers, self.labels, self.shapes
+        )
+
+    def tearDown(self):
+        self.feature_dir.cleanup()
+
+    def test_shape_without_reading(self):
+        feat = Features(self.feature_dir.name)
+        self.assertEqual(feat.shape('conv5'), self.stacked['conv5'].shape)
+        self.assertEqual(feat.shape('fc8'), self.stacked['fc8'].shape)
+
+    def test_feature_slice(self):
+        feat = Features(self.feature_dir.name)
+        assert_array_equal(
+            feat.get('conv5', feature_slice=np.s_[8:16]),
+            self.stacked['conv5'][:, 8:16],
+        )
+        assert_array_equal(
+            feat.get('fc8', feature_slice=np.s_[10:50]),
+            self.stacked['fc8'][:, 10:50],
+        )
+
+    def test_feature_slice_with_labels(self):
+        feat = Features(self.feature_dir.name)
+        labels = [self.labels[2], self.labels[0]]
+        assert_array_equal(
+            feat.get('conv5', label=labels, feature_slice=np.s_[8:16]),
+            self.stacked['conv5'][[2, 0]][:, 8:16],
+        )
+
+    def test_iter_chunks_reassembles(self):
+        feat = Features(self.feature_dir.name)
+        blocks = list(feat.iter_chunks('conv5', axis=1, size=7))
+        assert_array_equal(
+            np.concatenate([b for _, b in blocks], axis=1), self.stacked['conv5']
+        )
+
+    def test_unsliced_get_still_uses_the_layer_cache(self):
+        # The cache holds a whole layer; a sliced read must not replace it.
+        feat = Features(self.feature_dir.name)
+        cached = feat.get_features('conv5')
+        feat.get('conv5', feature_slice=np.s_[0:2])
+        self.assertIs(feat.get_features('conv5'), cached)
+
+
+class TestFeaturesFeatureIndex(unittest.TestCase):
+    """Unit-index selection.
+
+    NOTE: the happy path is not covered here. `feature_index` files are struct
+    .mat files, which the current reader (_mat_v73.loadmat_key) cannot load --
+    it handles dense arrays only, so a struct raises TypeError. That is a
+    pre-existing regression from the hdf5storage -> h5py read-path change
+    (issue #106), independent of the storage backends, and is left untouched
+    here rather than silently changed. Only the unambiguous case is asserted.
+    """
+
+    def setUp(self):
+        self.labels = ['img0001', 'img0002', 'img0003']
+        self.feature_dir = tempfile.TemporaryDirectory()
+        prepare_mat_features(
+            self.feature_dir.name, ['fc8'], self.labels, [(1, 20)]
+        )
+
+    def tearDown(self):
+        self.feature_dir.cleanup()
+
+    def test_missing_index_file_raises(self):
+        with self.assertRaises(RuntimeError):
+            Features(self.feature_dir.name, feature_index='/no/such/file.mat')
+
+
 class TestSaveFeature(unittest.TestCase):
     def test_save_feature_warns_future(self):
         # save_feature writes a MATLAB-compatible v7.3 .mat file, a path that is

From 186204893452491ca33c0851150a9c1dd2e8670d Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Fri, 18 Sep 2026 05:03:57 +0000
Subject: [PATCH 04/11] perf(dataform): trim HDF5 chunk edge padding

HDF5 allocates whole chunks, so an extent that divides its axis unevenly
pads the edge chunks. For (200, 256, 13, 13) float32 that cost 41.2 MiB
on disk for 33.0 MiB of data.

Shrink each extent to the smallest one needing the same number of chunks.
That can only lower the chunk size, so the byte budget still holds, and
it brings the same file to 34.1 MiB. Covered by a test capping the
overhead at 10%.
---
 bdpy/dataform/_feature_chunking.py      | 24 +++++++++++++++++++++---
 docs/dataform_features.md               |  9 ++++++---
 tests/dataform/test_feature_chunking.py | 12 ++++++++++++
 3 files changed, 39 insertions(+), 6 deletions(-)

diff --git a/bdpy/dataform/_feature_chunking.py b/bdpy/dataform/_feature_chunking.py
index 44b1b45..c3a057e 100644
--- a/bdpy/dataform/_feature_chunking.py
+++ b/bdpy/dataform/_feature_chunking.py
@@ -44,7 +44,9 @@ def choose_chunk_shape(
 
     Within the budget the two chunked axes are split roughly evenly in element
     count, so that neither a sample-wise read nor a channel-wise read degenerates
-    into reading the entire dataset.
+    into reading the entire dataset. Each extent is then shrunk to the smallest
+    one needing the same number of chunks, which keeps HDF5 from padding the
+    edge chunks and inflating the file.
 
     Parameters
     ----------
@@ -70,9 +72,9 @@ def choose_chunk_shape(
     Examples
     --------
     >>> choose_chunk_shape((1200, 1000), np.dtype(np.float32))
-    (512, 512)
+    (400, 500)
     >>> choose_chunk_shape((1200, 256, 13, 13), np.dtype(np.float32))
-    (39, 39, 13, 13)
+    (39, 37, 13, 13)
     >>> choose_chunk_shape((50, 1000), np.dtype(np.float32))
     (50, 1000)
     """
@@ -105,4 +107,20 @@ def choose_chunk_shape(
         c0 = min(shape[0], c0)
     c1 = min(shape[1], max(c1, budget_cells // c0))
 
+    # HDF5 allocates whole chunks, so a chunk extent that divides its axis
+    # unevenly pads the edge chunks and inflates the file. Shrinking the extent
+    # to the smallest one that needs the same number of chunks removes most of
+    # that padding, and can only lower the chunk size, so the budget still holds.
+    c1 = _snap(shape[1], c1)
+    if n_samples_known:
+        c0 = _snap(shape[0], c0)
+
     return (c0, c1, *shape[2:])
+
+
+def _snap(length: int, chunk: int) -> int:
+    """Smallest extent covering `length` in the same number of chunks as `chunk`."""
+    if length < 1 or chunk >= length:
+        return chunk
+    n_chunks = -(-length // chunk)  # ceil
+    return -(-length // n_chunks)   # ceil
diff --git a/docs/dataform_features.md b/docs/dataform_features.md
index f5d0131..9d299c9 100644
--- a/docs/dataform_features.md
+++ b/docs/dataform_features.md
@@ -156,7 +156,8 @@ read: HDF5 reads whole chunks even when only a few of their elements are
 selected. By default bdpy derives it from a 1 MiB budget, splitting it between
 the sample axis and the outermost feature axis and keeping trailing spatial axes
 whole -- e.g. `(1200, 256, 13, 13)` float32 becomes chunks of
-`(39, 39, 13, 13)`.
+`(39, 37, 13, 13)`. Extents are chosen to divide each axis as evenly as
+possible, so HDF5 does not pad the edge chunks.
 
 That means reading a *single* stimulus costs one chunk row (~39 stimuli of I/O),
 which matters for `FeaturesDataset`-style per-sample access. Tune it if your
@@ -168,8 +169,10 @@ save_features(path, array, labels, chunks=(1, 256, 13, 13))        # per-stimulu
 ```
 
 Compression is off by default, since every compressed chunk costs a
-decompression on the way out. Pass `compression='gzip'` or `compression='lzf'`
-when size matters more than read speed.
+decompression on the way out. Note that this makes the new format *larger* on
+disk than a legacy `.mat` tree, whose files `hdf5storage` compresses by default.
+Pass `compression='gzip'` or `compression='lzf'` when size matters more than
+read speed.
 
 ### Format
 
diff --git a/tests/dataform/test_feature_chunking.py b/tests/dataform/test_feature_chunking.py
index 579f70c..2ad6719 100644
--- a/tests/dataform/test_feature_chunking.py
+++ b/tests/dataform/test_feature_chunking.py
@@ -99,6 +99,18 @@ def test_unknown_sample_count_is_not_capped(self):
         )
         self.assertGreater(chunk[0], 1)
 
+    def test_edge_padding_is_small(self):
+        # HDF5 allocates whole chunks, so uneven extents pad the edges and
+        # inflate the file. Keep that overhead under 10%.
+        dtype = np.dtype(np.float32)
+        for shape in self.shapes:
+            chunk = choose_chunk_shape(shape, dtype)
+            allocated = 1
+            for s, c in zip(shape, chunk):
+                allocated *= -(-s // c) * c
+            with self.subTest(shape=shape):
+                self.assertLessEqual(allocated, int(np.prod(shape) * 1.1) + 1)
+
     def test_deterministic(self):
         shape = (1200, 256, 13, 13)
         dtype = np.dtype(np.float32)

From 2f19f81868ae6b22e8d59d93ff5c17ff0d5ddda9 Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Thu, 24 Sep 2026 02:22:59 +0000
Subject: [PATCH 05/11] fix(dataform): reject ambiguous slicing, make writes
 atomic, verify labels

Three review findings on the chunked feature storage, all of which returned
wrong data or lost it silently.

feature_slice + feature_index returned the wrong units. get() sliced the
feature axes first, then applied a unit index that addresses the flattened
FULL feature space to the already-narrowed array: silently the wrong units
when the index fit, IndexError when it did not. Refuse the combination with
ValueError instead of guessing. iter_chunks' existing feature_index guard
moves from RuntimeError to ValueError so the two read the same. The standalone
feature_index defects are untouched and remain out of scope.

A failed conversion published a truncated file. FeatureWriter opened the final
path with mode "w", truncating it up front, and the converter closed the writer
from a finally block, so an error mid-layer left a short file exactly where the
skip-if-exists check would take it for a finished one -- permanently, across
re-runs. Writes now stage in a hidden .partial sibling and publish with a single
os.replace, and close()/abort() split commit from discard, with __exit__ routing
to abort() when the body raises. The converter uses the context manager, so a
failure leaves nothing behind and re-running just works. An overwrite keeps the
old file readable until the swap. save_features and FeatureWriter also stop
clobbering: an existing target raises FileExistsError unless overwrite=True.
Staging opens with mode "x" rather than going through tempfile.mkstemp, so
published files keep umask permissions instead of being tightened to 0600, and
the scratch name is dotted and non-.h5 so a concurrent reader's *.h5 scan never
sees it. A weakref finalizer removes the scratch file if a writer is dropped.

Layers were not checked for agreeing labels. HDF5FeatureStore read labels from
the first layer only and reused that label -> row map against every layer's
/features, so a layer written in a different order returned another stimulus'
rows. Validate at construction that every layer holds the same labels in the
same order, reporting which layer differs and how, matching the invariant
MatFeatureStore already enforces.

Tests cover all three, and were checked to fail against the code before this
commit.
---
 bdpy/dataform/_feature_store.py      |  69 +++++++--
 bdpy/dataform/feature_hdf5.py        | 186 +++++++++++++++++++---
 bdpy/dataform/features.py            |  23 ++-
 docs/dataform_features.md            |  19 +++
 tests/dataform/test_feature_hdf5.py  | 220 +++++++++++++++++++++++++++
 tests/dataform/test_feature_store.py |  69 +++++++++
 tests/dataform/test_features.py      |  38 ++++-
 7 files changed, 587 insertions(+), 37 deletions(-)

diff --git a/bdpy/dataform/_feature_store.py b/bdpy/dataform/_feature_store.py
index 1db4fa4..a6810a4 100644
--- a/bdpy/dataform/_feature_store.py
+++ b/bdpy/dataform/_feature_store.py
@@ -420,8 +420,8 @@ def __init__(self, dpath: str):
         )
         if not self._layers:
             raise RuntimeError("No .{} feature file found in {}".format(HDF5_EXT, dpath))
-        self._labels: Optional[List[str]] = None
-        self._label_index: Optional[Dict[str, int]] = None
+        self._labels: List[str] = []
+        self._label_index: Dict[str, int] = {}
         self._validate_all()
 
     @property
@@ -430,12 +430,6 @@ def layers(self) -> List[str]:
 
     @property
     def labels(self) -> List[str]:
-        if self._labels is None:
-            with self._open(self._layers[0]) as f:
-                self._labels = _decode_labels(f[LABELS_DATASET][()])
-            self._label_index = {
-                label: i for i, label in enumerate(self._labels)
-            }
         return self._labels
 
     def path(self, layer: str) -> str:
@@ -506,8 +500,6 @@ def _read_rows(
         return np.asarray(dset[(rows, *indexers)])
 
     def _row_indices(self, labels: Sequence[str]) -> np.ndarray:
-        _ = self.labels  # populate the label -> row index map
-        assert self._label_index is not None
         try:
             return np.array([self._label_index[label] for label in labels], dtype=int)
         except KeyError as exc:
@@ -519,9 +511,32 @@ def _open(self, layer: str) -> h5py.File:
         return h5py.File(self.path(layer), "r")
 
     def _validate_all(self) -> None:
+        """Validate every layer file and pin down the shared label sequence.
+
+        A row index built from one layer is used against every layer's
+        ``/features``, so the layers must agree on the labels *and on their
+        order*; otherwise a label lookup would silently return another
+        stimulus' row. The legacy ``.mat`` backend enforces the same invariant
+        in :meth:`MatFeatureStore._collect_labels`.
+        """
+        reference: Optional[List[str]] = None
+        reference_layer = ""
         for layer in self._layers:
             with self._open(layer) as f:
                 _validate_format(f, self.path(layer))
+                labels = _decode_labels(f[LABELS_DATASET][()])
+            if reference is None:
+                reference, reference_layer = labels, layer
+            elif labels != reference:
+                raise RuntimeError(
+                    _label_mismatch_message(
+                        self._dpath, reference_layer, reference, layer, labels
+                    )
+                )
+
+        assert reference is not None  # self._layers is non-empty
+        self._labels = reference
+        self._label_index = {label: i for i, label in enumerate(reference)}
 
 
 def _validate_format(f: h5py.File, path: str) -> None:
@@ -548,6 +563,40 @@ def _validate_format(f: h5py.File, path: str) -> None:
             raise RuntimeError("{} has no /{} dataset".format(path, name))
 
 
+def _label_mismatch_message(
+    dpath: str,
+    reference_layer: str,
+    reference: List[str],
+    layer: str,
+    labels: List[str],
+) -> str:
+    """Describe *how* two layers' label sequences differ, not just that they do."""
+    if len(labels) != len(reference):
+        detail = "{} has {} labels, {} has {}".format(
+            reference_layer, len(reference), layer, len(labels)
+        )
+    elif sorted(labels) == sorted(reference):
+        first = next(
+            i for i, (a, b) in enumerate(zip(reference, labels)) if a != b
+        )
+        detail = (
+            "same labels in a different order; first difference at index {}: "
+            "{} has {!r}, {} has {!r}".format(
+                first, reference_layer, reference[first], layer, labels[first]
+            )
+        )
+    else:
+        only_ref = sorted(set(reference) - set(labels))[:3]
+        only_this = sorted(set(labels) - set(reference))[:3]
+        detail = "different labels; only in {}: {}, only in {}: {}".format(
+            reference_layer, only_ref or "-", layer, only_this or "-"
+        )
+    return (
+        "Inconsistent labels across layers in {}: {}. Every layer must hold the "
+        "same stimulus labels in the same order.".format(dpath, detail)
+    )
+
+
 def _decode_labels(raw: np.ndarray) -> List[str]:
     """Decode an HDF5 string dataset into a list of str (h5py 3.x gives bytes)."""
     return [
diff --git a/bdpy/dataform/feature_hdf5.py b/bdpy/dataform/feature_hdf5.py
index 8a5b3a4..ef39795 100644
--- a/bdpy/dataform/feature_hdf5.py
+++ b/bdpy/dataform/feature_hdf5.py
@@ -45,6 +45,8 @@
 """
 
 import os
+import uuid
+import weakref
 from types import TracebackType
 from typing import Iterable, Optional, Sequence, Tuple, Type
 
@@ -85,6 +87,62 @@ def _write_header(f: h5py.File, layer: Optional[str]) -> None:
         f.attrs["layer"] = layer
 
 
+#: Suffix of the scratch file a write is staged in.
+_PARTIAL_SUFFIX = ".partial"
+#: Cap on the copied basename so the staged name stays inside NAME_MAX.
+_MAX_STEM = 64
+
+
+def _staging_path(path: str) -> str:
+    """Name a scratch file next to `path`, invisible to feature-directory scans.
+
+    The name is dotted and does not end in ``.h5``, so neither
+    ``glob('*.h5')`` in :class:`~bdpy.dataform._feature_store.HDF5FeatureStore`
+    nor :func:`~bdpy.dataform._feature_store.detect_format` picks it up as a
+    layer while a write is in flight. The random component lets two processes
+    stage the same layer without colliding.
+    """
+    directory, name = os.path.split(os.path.abspath(path))
+    return os.path.join(
+        directory,
+        ".{}.{}{}".format(name[:_MAX_STEM], uuid.uuid4().hex, _PARTIAL_SUFFIX),
+    )
+
+
+def _prepare_target(path: str, overwrite: bool) -> str:
+    """Refuse to clobber, then name a scratch file next to `path`.
+
+    Writing goes to a scratch file in the *same directory*, so the final
+    :func:`os.replace` is an atomic rename on the same filesystem: either the
+    finished file appears at `path`, or nothing does. A half-written file must
+    never be left where a reader -- or the converter's skip-if-exists check --
+    would take it for a complete one.
+
+    The caller opens the returned path with mode ``"x"``, which both creates it
+    exclusively and lets HDF5 apply the process umask, so published files keep
+    the permissions they had before staging existed.
+    """
+    if os.path.exists(path) and not overwrite:
+        raise FileExistsError(
+            "{} already exists. Pass overwrite=True to replace it.".format(path)
+        )
+    _makedirs_for(path)
+    return _staging_path(path)
+
+
+def _promote(tmp_path: str, path: str) -> None:
+    """Move the finished scratch file into place (atomic on one filesystem)."""
+    os.replace(tmp_path, path)
+
+
+def _discard(tmp_path: str) -> None:
+    """Remove the scratch file, ignoring an already-vanished one."""
+    try:
+        os.remove(tmp_path)
+    except OSError:
+        pass
+
+
 def save_features(
     path: str,
     features: np.ndarray,
@@ -94,13 +152,15 @@ def save_features(
     target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
     compression: Optional[str] = None,
     dtype: Optional[np.dtype] = None,
+    overwrite: bool = False,
 ) -> None:
     """Write a whole layer to a chunked HDF5 feature file.
 
     Parameters
     ----------
     path : str
-        Output file. Overwritten if it exists.
+        Output file. An existing file is kept unless `overwrite` is True. The
+        write is atomic: on failure nothing is left at `path`.
     features : numpy.ndarray
         Feature array of shape ``(n_samples, *feature_shape)``.
     labels : sequence of str
@@ -117,12 +177,16 @@ def save_features(
         fast.
     dtype : numpy.dtype, optional
         Dtype to store. Defaults to the dtype of `features`.
+    overwrite : bool, optional
+        Replace an existing file at `path` (default: False).
 
     Raises
     ------
     ValueError
         If `labels` does not have one entry per sample, or `features` has fewer
         than two axes.
+    FileExistsError
+        If `path` exists and `overwrite` is False.
     """
     features = np.asarray(features)
     if features.ndim < 2:
@@ -144,16 +208,21 @@ def save_features(
     if layer is None:
         layer = os.path.splitext(os.path.basename(path))[0]
 
-    _makedirs_for(path)
-    with h5py.File(path, "w") as f:
-        _write_header(f, layer)
-        f.create_dataset(
-            FEATURES_DATASET,
-            data=features.astype(dtype, copy=False),
-            chunks=chunks,
-            compression=compression,
-        )
-        f.create_dataset(LABELS_DATASET, data=labels, dtype=_string_dtype())
+    tmp_path = _prepare_target(path, overwrite)
+    try:
+        with h5py.File(tmp_path, "x") as f:
+            _write_header(f, layer)
+            f.create_dataset(
+                FEATURES_DATASET,
+                data=features.astype(dtype, copy=False),
+                chunks=chunks,
+                compression=compression,
+            )
+            f.create_dataset(LABELS_DATASET, data=labels, dtype=_string_dtype())
+    except BaseException:
+        _discard(tmp_path)
+        raise
+    _promote(tmp_path, path)
 
 
 class FeatureWriter:
@@ -161,12 +230,18 @@ class FeatureWriter:
 
     Feature extraction produces one stimulus at a time, so the datasets are
     created resizable (``maxshape=(None, *feature_shape)``) and grown in blocks
-    as samples arrive. Use it as a context manager, or call :meth:`close`.
+    as samples arrive.
+
+    Writing goes to a temporary file next to `path`; :meth:`close` moves it into
+    place and :meth:`abort` throws it away, so a failed write leaves nothing at
+    `path`. **Prefer the context manager**, which aborts when the block raises --
+    a bare ``try/finally: writer.close()`` would publish a half-written file.
 
     Parameters
     ----------
     path : str
-        Output file. Overwritten if it exists.
+        Output file. An existing file is kept unless `overwrite` is True, and
+        nothing is written there until :meth:`close` succeeds.
     feature_shape : sequence of int
         Shape of a single sample's features, without the sample axis.
     dtype : numpy.dtype
@@ -181,6 +256,13 @@ class FeatureWriter:
         Per-chunk byte budget used when `chunks` is not given.
     compression : str, optional
         h5py compression filter. ``None`` (default) stores uncompressed.
+    overwrite : bool, optional
+        Replace an existing file at `path` (default: False).
+
+    Raises
+    ------
+    FileExistsError
+        If `path` exists and `overwrite` is False.
 
     Examples
     --------
@@ -198,6 +280,7 @@ def __init__(
         chunks: Optional[Tuple[int, ...]] = None,
         target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
         compression: Optional[str] = None,
+        overwrite: bool = False,
     ):
         self._feature_shape = tuple(int(s) for s in feature_shape)
         if not self._feature_shape:
@@ -215,8 +298,18 @@ def __init__(
         if layer is None:
             layer = os.path.splitext(os.path.basename(path))[0]
 
-        _makedirs_for(path)
-        self._file: Optional[h5py.File] = h5py.File(path, "w")
+        # Fail before doing any work if the target is occupied, then write to
+        # a scratch file so `path` stays untouched until close() succeeds.
+        self._path = path
+        self._aborted = False
+        self._tmp_path: Optional[str] = _prepare_target(path, overwrite)
+        self._file: Optional[h5py.File] = h5py.File(self._tmp_path, "x")
+        # If the writer is dropped without close() or abort(), the scratch file
+        # would otherwise linger next to the output. weakref.finalize (rather
+        # than __del__) runs exactly once, is detachable at publish time, and
+        # closes over the path instead of the writer, so it can never resurrect
+        # the object or outlive a name it no longer owns.
+        self._finalizer = weakref.finalize(self, _discard, self._tmp_path)
         _write_header(self._file, layer)
         self._features = self._file.create_dataset(
             FEATURES_DATASET,
@@ -297,10 +390,47 @@ def extend(self, features: np.ndarray, labels: Sequence[str]) -> None:
         self._n = new_n
 
     def close(self) -> None:
-        """Close the underlying file. Idempotent."""
+        """Finish the file and move it into place. Idempotent.
+
+        This is the commit: call it only once the data is complete, because it
+        publishes whatever has been written so far. To throw a partial file
+        away, call :meth:`abort` instead. ``try: ... finally: writer.close()``
+        is therefore wrong -- it would publish a truncated file -- so use the
+        writer as a context manager, which routes to :meth:`abort` when the
+        body raises.
+
+        Raises
+        ------
+        RuntimeError
+            If the writer was already aborted.
+        """
+        if self._aborted:
+            raise RuntimeError(
+                "FeatureWriter was aborted; nothing was written to {}".format(
+                    self._path
+                )
+            )
         if self._file is not None:
             self._file.close()
             self._file = None
+        if self._tmp_path is not None:
+            _promote(self._tmp_path, self._path)
+            self._tmp_path = None
+            self._finalizer.detach()  # the scratch file is now the output
+
+    def abort(self) -> None:
+        """Discard the partial file without touching the target path.
+
+        Idempotent, and a no-op once :meth:`close` has published the file --
+        it must never remove a published output.
+        """
+        if self._file is not None:
+            self._file.close()
+            self._file = None
+        if self._tmp_path is not None:
+            self._finalizer()  # runs _discard exactly once
+            self._tmp_path = None
+            self._aborted = True
 
     def __enter__(self) -> "FeatureWriter":
         """Enter the context manager."""
@@ -312,8 +442,13 @@ def __exit__(
         exc: Optional[BaseException],
         tb: Optional[TracebackType],
     ) -> None:
-        """Close the file on leaving the context manager."""
-        self.close()
+        """Publish the file on success, discard it if the block raised."""
+        if self._tmp_path is None:
+            return  # the body already called close() or abort() itself
+        if exc_type is None:
+            self.close()
+        else:
+            self.abort()
 
 
 def convert_features_to_hdf5(
@@ -341,7 +476,8 @@ def convert_features_to_hdf5(
         Legacy directory, ``<src_dir>/<layer>/<label>.<ext>``.
     dst_dir : str
         Output directory; ``<dst_dir>/<layer>.h5`` is written per layer. Created
-        if missing.
+        if missing. Each layer is written atomically, so a failed conversion
+        leaves no file behind and can simply be re-run.
     layers : iterable of str, optional
         Layers to convert. Defaults to every layer found in `src_dir`.
     ext : str, optional
@@ -384,7 +520,11 @@ def convert_features_to_hdf5(
             continue
 
         full_shape = store.shape(layer)
-        writer = FeatureWriter(
+        # The context manager is load-bearing: if a batch fails, it discards the
+        # partial file instead of publishing it. A published partial file would
+        # be taken for a finished one by the skip check above, and the layer
+        # would stay silently truncated across re-runs.
+        with FeatureWriter(
             out_path,
             feature_shape=full_shape[1:],
             dtype=store.dtype(layer),
@@ -392,13 +532,11 @@ def convert_features_to_hdf5(
             n_samples=full_shape[0],
             target_chunk_bytes=target_chunk_bytes,
             compression=compression,
-        )
-        try:
+            overwrite=overwrite,
+        ) as writer:
             for start in range(0, len(all_labels), batch_size):
                 batch = all_labels[start:start + batch_size]
                 writer.extend(store.read(layer, batch), batch)
-        finally:
-            writer.close()
 
         if verbose:
             print("Saved {} ({} samples).".format(out_path, full_shape[0]))
diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py
index 44fc63a..7b8d049 100644
--- a/bdpy/dataform/features.py
+++ b/bdpy/dataform/features.py
@@ -197,7 +197,26 @@ def get(
         >>> features.get('conv5')                                  # doctest: +SKIP
         >>> features.get('conv5', label=['img0001', 'img0002'])     # doctest: +SKIP
         >>> features.get('conv5', feature_slice=np.s_[128:256])     # doctest: +SKIP
+
+        Raises
+        ------
+        ValueError
+            If `feature_slice` is combined with a unit index (`feature_index`).
+            The index addresses the flattened *full* feature space, so it is
+            meaningless against an already-sliced array.
         """
+        if feature_slice is not None and self.__feat_index_table is not None:
+            # The unit index addresses the flattened full feature space, while
+            # feature_slice has already narrowed it. Applying one to the other
+            # would silently select the wrong units (or raise IndexError),
+            # so refuse the combination instead of guessing.
+            raise ValueError(
+                'feature_slice cannot be combined with feature_index: the unit '
+                'index addresses the full feature space, not a slice of it. '
+                'Call get() without feature_slice and slice the result, or '
+                'construct Features without feature_index.'
+            )
+
         if label is None and feature_slice is None:
             return self.get_features(layer)
 
@@ -244,7 +263,7 @@ def iter_chunks(
 
         Raises
         ------
-        RuntimeError
+        ValueError
             If a unit index (`feature_index`) is in use, which flattens the
             feature axes and so has no meaningful per-axis iteration.
 
@@ -254,7 +273,7 @@ def iter_chunks(
         ...     out[:, sl] = transform(block)
         """
         if self.__feat_index_table is not None:
-            raise RuntimeError(
+            raise ValueError(
                 'iter_chunks is not supported together with feature_index, '
                 'which flattens the feature axes'
             )
diff --git a/docs/dataform_features.md b/docs/dataform_features.md
index 9d299c9..9d8486e 100644
--- a/docs/dataform_features.md
+++ b/docs/dataform_features.md
@@ -100,6 +100,11 @@ feat = features.get(
     feature_slice=np.s_[128:256],
 )
 
+# NOTE: feature_slice cannot be combined with a unit index (`feature_index`).
+# The index addresses the flattened full feature space, so applying it to an
+# already-sliced array would select the wrong units; the combination raises
+# ValueError instead.
+
 # Full shape without reading anything
 n_stimuli, *feature_shape = features.shape('conv5')
 ```
@@ -126,6 +131,14 @@ from bdpy.dataform import save_features
 save_features('features/conv5.h5', array, labels)
 ```
 
+Writing never clobbers: an existing file raises `FileExistsError` unless you pass
+`overwrite=True`. Writes are also atomic -- the file is built in a hidden
+`.<name>.<id>.partial` sibling and moved into place only once it is complete --
+so a failed write leaves nothing behind and can simply be re-run, and an
+overwrite keeps the old file readable until the moment it is replaced. (A
+process killed outright cannot clean up after itself, so a stray `.partial` file
+may be left for you to delete.)
+
 Or incrementally, which is what feature extraction needs since it produces one
 stimulus at a time:
 
@@ -148,6 +161,12 @@ convert_features_to_hdf5('/path/to/features_mat', '/path/to/features_h5')
 
 The converter reads through the same backend `Features` uses for the legacy
 layout, and streams in batches, so a layer is never held in memory in full.
+Existing `<layer>.h5` files are skipped unless `overwrite=True`, and a layer that
+fails to convert leaves no file behind, so re-running picks up where it stopped.
+
+All layers in a directory must hold the same stimulus labels in the same order;
+a mismatch is rejected when the directory is opened, as it is for the legacy
+layout.
 
 ### Chunk shape
 
diff --git a/tests/dataform/test_feature_hdf5.py b/tests/dataform/test_feature_hdf5.py
index 2b661cf..b93845e 100644
--- a/tests/dataform/test_feature_hdf5.py
+++ b/tests/dataform/test_feature_hdf5.py
@@ -1,4 +1,7 @@
+import gc
+import glob
 import os
+import stat
 import tempfile
 import unittest
 import warnings
@@ -9,6 +12,7 @@
 
 from bdpy.dataform import Features
 from bdpy.dataform._feature_store import (
+    MatFeatureStore,
     FORMAT_ATTR,
     FORMAT_NAME,
     FORMAT_VERSION_ATTR,
@@ -24,6 +28,12 @@
 from .test_features import prepare_mat_features
 
 
+def _current_umask():
+    mask = os.umask(0)
+    os.umask(mask)
+    return mask
+
+
 class TestSaveFeatures(unittest.TestCase):
     def setUp(self):
         self.tmpdir = tempfile.TemporaryDirectory()
@@ -205,6 +215,216 @@ def test_small_batches_produce_the_same_file(self):
         assert_array_equal(Features(self.h5dir).get('conv5'), self.stacked['conv5'])
 
 
+class TestAtomicWriteAndOverwrite(unittest.TestCase):
+    """A failed write must leave nothing behind, and must never clobber.
+
+    The converter treats an existing <layer>.h5 as finished and skips it, so a
+    half-written file published at the target path would make a truncated layer
+    permanent across re-runs.
+    """
+
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.matdir = os.path.join(self.tmpdir.name, 'mat')
+        self.h5dir = os.path.join(self.tmpdir.name, 'h5')
+        os.makedirs(self.matdir)
+        self.labels = ['img%04d' % i for i in range(10)]
+        self.stacked = prepare_mat_features(
+            self.matdir, ['conv5'], self.labels, [(1, 16, 3, 3)]
+        )
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+        self.data = np.random.rand(10, 8).astype(np.float32)
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def _leftovers(self, directory):
+        return sorted(os.listdir(directory))
+
+    # --- overwrite guard -------------------------------------------------
+
+    def test_save_features_refuses_to_clobber(self):
+        save_features(self.path, self.data, self.labels)
+        with self.assertRaises(FileExistsError):
+            save_features(self.path, self.data * 2, self.labels)
+        # The original must be untouched.
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_save_features_overwrite_replaces(self):
+        save_features(self.path, self.data, self.labels)
+        save_features(self.path, self.data * 2, self.labels, overwrite=True)
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data * 2
+        )
+
+    def test_writer_refuses_to_clobber(self):
+        save_features(self.path, self.data, self.labels)
+        with self.assertRaises(FileExistsError):
+            FeatureWriter(self.path, (8,), np.float32)
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_writer_overwrite_replaces(self):
+        save_features(self.path, self.data, self.labels)
+        with FeatureWriter(self.path, (8,), np.float32, overwrite=True) as w:
+            w.extend(self.data * 3, self.labels)
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data * 3
+        )
+
+    # --- atomicity -------------------------------------------------------
+
+    def test_target_is_absent_until_close(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.append(self.data[0], self.labels[0])
+        self.assertFalse(os.path.exists(self.path))
+        writer.close()
+        self.assertTrue(os.path.exists(self.path))
+
+    def test_abort_leaves_nothing(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.append(self.data[0], self.labels[0])
+        writer.abort()
+        self.assertFalse(os.path.exists(self.path))
+        self.assertEqual(self._leftovers(self.tmpdir.name), ['mat'])
+
+    def test_context_manager_aborts_on_exception(self):
+        with self.assertRaises(ZeroDivisionError):
+            with FeatureWriter(self.path, (8,), np.float32) as writer:
+                writer.append(self.data[0], self.labels[0])
+                raise ZeroDivisionError
+        self.assertFalse(os.path.exists(self.path))
+        self.assertEqual(self._leftovers(self.tmpdir.name), ['mat'])
+
+    def test_close_and_abort_are_idempotent(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.extend(self.data, self.labels)
+        writer.close()
+        writer.close()
+        writer.abort()  # must not delete the published file
+        self.assertTrue(os.path.exists(self.path))
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_failed_save_features_leaves_nothing(self):
+        # Mismatched labels raise after the target has been prepared.
+        with self.assertRaises(ValueError):
+            save_features(self.path, self.data, self.labels[:-1])
+        self.assertFalse(os.path.exists(self.path))
+        self.assertEqual(self._leftovers(self.tmpdir.name), ['mat'])
+
+    def test_published_file_is_group_readable(self):
+        # Staging must not tighten permissions: these files live on shared lab
+        # storage, and 0600 would lock collaborators out.
+        save_features(self.path, self.data, self.labels)
+        mode = stat.S_IMODE(os.stat(self.path).st_mode)
+        expected = 0o666 & ~_current_umask()
+        self.assertEqual(mode, expected)
+
+    def test_staging_file_is_not_mistaken_for_a_layer(self):
+        # A directory being written into is also a directory someone may read.
+        # The scratch file must be invisible to the *.h5 scan, or a concurrent
+        # reader would either invent a bogus layer or fail validation.
+        writer = FeatureWriter(os.path.join(self.h5dir, 'conv5.h5'), (8,), np.float32)
+        try:
+            writer.append(self.data[0], self.labels[0])
+            self.assertEqual(glob.glob(os.path.join(self.h5dir, '*.h5')), [])
+            with self.assertRaises(RuntimeError):
+                HDF5FeatureStore(self.h5dir)  # "no .h5 feature file found"
+        finally:
+            writer.abort()
+
+    def test_failed_overwrite_preserves_the_previous_file(self):
+        # The old file stays readable for the whole write and is swapped only at
+        # the end, so a failure mid-overwrite loses neither old nor new.
+        save_features(self.path, self.data, self.labels)
+        with self.assertRaises(ZeroDivisionError):
+            with FeatureWriter(self.path, (8,), np.float32, overwrite=True) as w:
+                w.extend(self.data * 9, self.labels)
+                raise ZeroDivisionError
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_dropped_writer_does_not_leak_a_scratch_file(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.append(self.data[0], self.labels[0])
+        del writer
+        gc.collect()
+        self.assertEqual(self._leftovers(self.tmpdir.name), ['mat'])
+
+    def test_close_after_abort_is_refused(self):
+        # Silently doing nothing would let a caller believe it published.
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.abort()
+        with self.assertRaises(RuntimeError):
+            writer.close()
+        self.assertFalse(os.path.exists(self.path))
+
+    def test_zero_samples_round_trips(self):
+        # An empty layer is representable; refusing it would add a failure mode
+        # on the success path for something the caller can check themselves.
+        with FeatureWriter(self.path, (8,), np.float32):
+            pass
+        store = HDF5FeatureStore(self.tmpdir.name)
+        self.assertEqual(store.labels, [])
+        self.assertEqual(store.read('conv5').shape, (0, 8))
+
+    # --- the converter ---------------------------------------------------
+
+    def _fail_on_second_batch(self):
+        """Patch MatFeatureStore.read to blow up partway through a layer."""
+        original = MatFeatureStore.read
+        state = {'calls': 0}
+
+        def flaky(store, layer, labels=None, feature_slice=None):
+            state['calls'] += 1
+            if state['calls'] > 1:
+                raise OSError('simulated read failure')
+            return original(store, layer, labels, feature_slice)
+
+        return original, flaky
+
+    def test_failed_conversion_leaves_no_file(self):
+        original, flaky = self._fail_on_second_batch()
+        MatFeatureStore.read = flaky
+        try:
+            with self.assertRaises(OSError):
+                convert_features_to_hdf5(
+                    self.matdir, self.h5dir, batch_size=3
+                )
+        finally:
+            MatFeatureStore.read = original
+
+        out_path = os.path.join(self.h5dir, 'conv5.h5')
+        self.assertFalse(os.path.exists(out_path))
+        # No scratch file left behind either.
+        self.assertEqual(self._leftovers(self.h5dir), [])
+
+    def test_rerun_after_failure_succeeds(self):
+        # The regression this guards: a published partial file would be taken
+        # for a finished one by the skip check and never repaired.
+        original, flaky = self._fail_on_second_batch()
+        MatFeatureStore.read = flaky
+        try:
+            with self.assertRaises(OSError):
+                convert_features_to_hdf5(
+                    self.matdir, self.h5dir, batch_size=3
+                )
+        finally:
+            MatFeatureStore.read = original
+
+        convert_features_to_hdf5(self.matdir, self.h5dir, batch_size=3)
+        assert_array_equal(
+            Features(self.h5dir).get('conv5'), self.stacked['conv5']
+        )
+
+
 class TestFormatValidation(unittest.TestCase):
     def setUp(self):
         self.tmpdir = tempfile.TemporaryDirectory()
diff --git a/tests/dataform/test_feature_store.py b/tests/dataform/test_feature_store.py
index ff7c8b1..e82e44c 100644
--- a/tests/dataform/test_feature_store.py
+++ b/tests/dataform/test_feature_store.py
@@ -219,6 +219,75 @@ def test_dataset_is_chunked_on_both_sliceable_axes(self):
         self.assertEqual(chunks[2:], shape[2:])  # spatial axes kept whole
 
 
+class TestCrossLayerLabelConsistency(unittest.TestCase):
+    """Every layer must hold the same labels in the same order.
+
+    A row index built from one layer is used against every layer's /features,
+    so a layer whose labels are ordered differently would make a label lookup
+    silently return another stimulus' row. The legacy .mat backend enforces the
+    same invariant, and the HDF5 backend must match it.
+    """
+
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def _write(self, layer, labels, n_features=8):
+        data = np.random.rand(len(labels), n_features)
+        save_features(
+            os.path.join(self.tmpdir.name, layer + '.h5'), data, labels,
+            layer=layer,
+        )
+        return data
+
+    def test_consistent_layers_are_accepted(self):
+        self._write('conv5', LABELS)
+        self._write('fc8', LABELS)
+        self.assertEqual(Features(self.tmpdir.name).labels, LABELS)
+
+    def test_different_order_is_rejected(self):
+        self._write('conv5', LABELS)
+        self._write('fc8', list(reversed(LABELS)))
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        message = str(ctx.exception)
+        self.assertIn('different order', message)
+        self.assertIn('fc8', message)
+
+    def test_different_content_is_rejected(self):
+        self._write('conv5', LABELS)
+        self._write('fc8', ['other%04d' % i for i in range(len(LABELS))])
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('different labels', str(ctx.exception))
+
+    def test_different_length_is_rejected(self):
+        self._write('conv5', LABELS)
+        self._write('fc8', LABELS[:-2])
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('labels', str(ctx.exception))
+
+    def test_features_construction_also_rejects(self):
+        # The guard must fire through the public entry point too.
+        self._write('conv5', LABELS)
+        self._write('fc8', list(reversed(LABELS)))
+        with self.assertRaises(RuntimeError):
+            Features(self.tmpdir.name)
+
+    def test_legacy_backend_rejects_the_same_way(self):
+        # Parity with the .mat layout, which raises when layers disagree.
+        matdir = os.path.join(self.tmpdir.name, 'mat')
+        os.makedirs(matdir)
+        prepare_mat_features(matdir, ['conv5'], LABELS, [(1, 8)])
+        prepare_mat_features(matdir, ['fc8'], LABELS[:-2], [(1, 8)])
+        with self.assertRaises(RuntimeError):
+            Features(matdir)
+
+
 class TestFormatDetection(unittest.TestCase):
     def setUp(self):
         warnings.simplefilter('ignore', FutureWarning)
diff --git a/tests/dataform/test_features.py b/tests/dataform/test_features.py
index 421e6fe..2bd68d5 100644
--- a/tests/dataform/test_features.py
+++ b/tests/dataform/test_features.py
@@ -200,13 +200,49 @@ class TestFeaturesFeatureIndex(unittest.TestCase):
     def setUp(self):
         self.labels = ['img0001', 'img0002', 'img0003']
         self.feature_dir = tempfile.TemporaryDirectory()
-        prepare_mat_features(
+        self.stacked = prepare_mat_features(
             self.feature_dir.name, ['fc8'], self.labels, [(1, 20)]
         )
 
     def tearDown(self):
         self.feature_dir.cleanup()
 
+    def _features_with_index(self):
+        """Features with a unit index installed.
+
+        NOTE: the index table is injected directly rather than loaded from a
+        file. Loading is broken on `dev` (struct .mat files are unreadable by
+        _mat_v73.loadmat_key) and fixing that is out of scope here, but the
+        guards below must still be covered.
+        """
+        feat = Features(self.feature_dir.name)
+        feat._Features__feat_index_table = {'fc8': np.array([0, 5, 11, 19])}
+        return feat
+
+    def test_feature_slice_with_feature_index_is_refused(self):
+        # The unit index addresses the flattened FULL feature space, so applying
+        # it to an already-sliced array would silently pick the wrong units.
+        feat = self._features_with_index()
+        with self.assertRaises(ValueError) as ctx:
+            feat.get('fc8', feature_slice=np.s_[0:10])
+        self.assertIn('feature_index', str(ctx.exception))
+
+        # ... including when labels are also given.
+        with self.assertRaises(ValueError):
+            feat.get('fc8', label=self.labels[0], feature_slice=np.s_[0:10])
+
+    def test_unsliced_get_still_works_with_feature_index(self):
+        # The guard must not break the supported path.
+        feat = self._features_with_index()
+        assert_array_equal(
+            feat.get('fc8'), self.stacked['fc8'][:, [0, 5, 11, 19]]
+        )
+
+    def test_iter_chunks_with_feature_index_is_refused(self):
+        feat = self._features_with_index()
+        with self.assertRaises(ValueError):
+            list(feat.iter_chunks('fc8'))
+
     def test_missing_index_file_raises(self):
         with self.assertRaises(RuntimeError):
             Features(self.feature_dir.name, feature_index='/no/such/file.mat')

From 9a8c84f2fc3d913a9883f2f57c1d3fa062c66752 Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Tue, 29 Sep 2026 07:18:26 +0000
Subject: [PATCH 06/11] fix(dataform): narrow feature_slice, validate schema,
 dedupe labels

A second review found six defects. Three of them came from one decision --
generalising feature_slice towards full NumPy indexing -- so they are
removed by narrowing the API rather than by teaching it NumPy's index
arithmetic.

feature_slice is now basic forward indexing only: positive-step slices
(negative start/stop included), integers, one Ellipsis, and tuples of
those. Everything else raises ValueError on both backends. h5py handles
exactly the allowed set with NumPy's semantics but rejects a negative
step, so anything wider either makes the layouts disagree or needs NumPy
reimplemented to keep them together. np.s_[::-1] used to work on .mat and
raise inside h5py; now both refuse it the same way. bool is rejected
before the integer check, since isinstance(True, int) is True in Python.

iter_chunks loses its feature_slice parameter. Composing an arbitrary
index with the per-slab slice is the same arithmetic problem, and it was
wrong: iter_chunks(feature_slice=np.s_[:-1], size=5) returned 46 feature
columns where the slice itself has 19, and an integer feature_slice
raised IndexError. Slice the blocks as they come out instead. This
deletes _compose_slice, _compose_two_slices, _sliced_shape and _is_fancy
outright.

Duplicate labels are now invalid. Labels map to row indices once, so
["a", "a"] resolved both rows to the last one and silently dropped the
first. The legacy layout cannot express duplicates -- the file name is
the label -- so rejecting them keeps the layouts equivalent. Refused on
read, in save_features, and in FeatureWriter both within one extend()
call and across earlier calls, checked before writing so a refused batch
is a no-op.

Schema validation is consolidated and tightened: the version attribute
must be present, integral and in 1..SUPPORTED_FORMAT_VERSION (missing and
zero were accepted before, contradicting the docs), /features and /labels
must have the expected rank, and their row counts must agree.

Format detection keys off the files that are actually present instead of
the bare existence of a subdirectory, so a stray directory no longer
flips a chunked tree to .mat. A directory holding both layouts now raises
instead of silently choosing.

save_features becomes a thin wrapper over FeatureWriter. Writing a whole
layer is the degenerate case of writing it incrementally, and keeping two
write paths is what let a zero-sample layer succeed through one and fail
through the other. The unified path costs nothing measurable: the same
layer is 34.08 MiB with identical chunks; only the sample axis becomes
extensible.

Also: drop the layer attribute and the layer= arguments, since the file
name carries it and nothing validated the attribute; rename _GROW_BLOCK
to _LABELS_CHUNK to match what it does; de-duplicate _load_array_with_key
and _determine_num_parallel; and state the multi-directory limit of
iter_chunks in its docstring instead of claiming one slab unconditionally.

choose_chunk_shape now rejects empty feature axes and a negative sample
axis, and clamps a zero-sample layer to a chunk of at least 1, which HDF5
requires. Its budget is documented as an upper bound only where the
budget is reachable, which is not so when the trailing axes alone exceed
it. Its tests assert invariants and no longer pin chunk values, padding
ratios or budget usage -- those are a performance policy that should stay
free to change.

Every added test was run against the code before this commit to confirm
it fails there.
---
 bdpy/dataform/_feature_chunking.py      |  36 ++-
 bdpy/dataform/_feature_store.py         | 324 ++++++++++++++++--------
 bdpy/dataform/feature_hdf5.py           |  91 ++++---
 bdpy/dataform/features.py               |  52 ++--
 docs/dataform_features.md               |  33 ++-
 tests/dataform/test_feature_chunking.py |  77 +++---
 tests/dataform/test_feature_hdf5.py     | 153 ++++++++++-
 tests/dataform/test_feature_store.py    | 117 ++++++++-
 8 files changed, 631 insertions(+), 252 deletions(-)

diff --git a/bdpy/dataform/_feature_chunking.py b/bdpy/dataform/_feature_chunking.py
index c3a057e..24eaf74 100644
--- a/bdpy/dataform/_feature_chunking.py
+++ b/bdpy/dataform/_feature_chunking.py
@@ -55,8 +55,11 @@ def choose_chunk_shape(
     dtype : numpy.dtype
         Array dtype; only its itemsize matters.
     target_bytes : int, optional
-        Upper bound on the size of one chunk in bytes
-        (default: ``DEFAULT_TARGET_CHUNK_BYTES``).
+        Size budget for one chunk, in bytes
+        (default: ``DEFAULT_TARGET_CHUNK_BYTES``). It is an upper bound only
+        when the budget can be met at all: the trailing axes are kept whole, so
+        a shape whose trailing axes alone exceed `target_bytes` necessarily
+        produces a larger chunk.
     n_samples_known : bool, optional
         Whether ``shape[0]`` is the final number of samples. Pass ``False`` when
         writing to a resizable dataset whose sample count is not yet known, so
@@ -65,9 +68,16 @@ def choose_chunk_shape(
     Returns
     -------
     tuple of int
-        Chunk shape, same length as ``shape``. Every element is at least 1 and,
-        when ``n_samples_known`` is True, at most the corresponding entry of
-        ``shape``.
+        Chunk shape, same length as ``shape``. Every element is at least 1, and
+        no larger than the corresponding entry of ``shape`` for axes of
+        non-zero size (a layer with zero samples still needs a chunk extent of
+        at least 1 on the sample axis, since HDF5 forbids a zero chunk dim).
+
+    Raises
+    ------
+    ValueError
+        If `shape` has fewer than two axes, any feature axis has size 0, the
+        sample axis is negative, or `target_bytes` is not positive.
 
     Examples
     --------
@@ -86,12 +96,22 @@ def choose_chunk_shape(
         )
     if target_bytes < 1:
         raise ValueError("target_bytes must be positive, got {}".format(target_bytes))
+    if shape[0] < 0:
+        raise ValueError("sample axis cannot be negative, got {}".format(shape))
+    if any(s < 1 for s in shape[1:]):
+        # A zero-width feature axis carries nothing and cannot be chunked.
+        # A zero-sample layer, by contrast, is meaningful and is allowed below.
+        raise ValueError(
+            "feature axes must all be non-empty, got shape {}".format(shape)
+        )
 
     itemsize = np.dtype(dtype).itemsize
 
     # The whole array fits in one chunk: nothing to gain from splitting it.
+    # Clamp to at least 1 per axis -- a zero-sample layer is legal but HDF5
+    # rejects a chunk with a zero dimension.
     if n_samples_known and _prod(shape) * itemsize <= target_bytes:
-        return shape
+        return tuple(max(1, s) for s in shape)
 
     # Bytes taken by one (sample, feature-0) cell, i.e. one element of the two
     # chunked axes with the trailing axes kept whole.
@@ -103,7 +123,7 @@ def choose_chunk_shape(
     # force a narrow feature chunk, and vice versa.
     c1 = min(shape[1], max(1, math.isqrt(budget_cells)))
     c0 = max(1, budget_cells // c1)
-    if n_samples_known:
+    if n_samples_known and shape[0] > 0:
         c0 = min(shape[0], c0)
     c1 = min(shape[1], max(c1, budget_cells // c0))
 
@@ -112,7 +132,7 @@ def choose_chunk_shape(
     # to the smallest one that needs the same number of chunks removes most of
     # that padding, and can only lower the chunk size, so the budget still holds.
     c1 = _snap(shape[1], c1)
-    if n_samples_known:
+    if n_samples_known and shape[0] > 0:
         c0 = _snap(shape[0], c0)
 
     return (c0, c1, *shape[2:])
diff --git a/bdpy/dataform/_feature_store.py b/bdpy/dataform/_feature_store.py
index a6810a4..0bdd3a6 100644
--- a/bdpy/dataform/_feature_store.py
+++ b/bdpy/dataform/_feature_store.py
@@ -50,14 +50,20 @@
 #: Extension of a chunked feature file.
 HDF5_EXT = "h5"
 
-# A slice spec for the feature axes: what ``numpy.s_[...]`` produces.
-FeatureSlice = Union[slice, int, Sequence[int], np.ndarray, Tuple, None]
+# A slice spec for the feature axes. Deliberately narrow: basic indexing only,
+# i.e. what ``numpy.s_[128:256]`` or ``numpy.s_[8:16, 1:4]`` produces. See
+# _validate_feature_slice for what is rejected and why.
+BasicIndexer = Union[slice, int, "ellipsis"]  # noqa: F821
+FeatureSlice = Union[BasicIndexer, Tuple[BasicIndexer, ...], None]
 
 
 def _load_array_with_key(key: str, path: str) -> np.ndarray:
-    # v5 .mat via scipy, v7.3 (HDF5) via h5py; avoids hdf5storage on the load
-    # path, which breaks under NumPy 2.0 (see bdpy/dataform/_mat_v73.py).
-    # NOTE: key comes first so that partial(..., 'feat') is a Pool.map callable.
+    """Load one array from a ``.mat`` file. Shared with :mod:`bdpy.dataform.features`.
+
+    v5 ``.mat`` via scipy, v7.3 (HDF5) via h5py; avoids hdf5storage on the load
+    path, which breaks under NumPy 2.0 (see :mod:`bdpy.dataform._mat_v73`).
+    `key` comes first so ``partial(..., 'feat')`` is a ``Pool.map`` callable.
+    """
     return _mat_v73.loadmat_key(path, key)
 
 
@@ -77,19 +83,75 @@ def _determine_num_parallel(num_files: int) -> int:
 
 
 def _normalize_feature_slice(feature_slice: FeatureSlice) -> Tuple:
-    """Normalize a feature-axis slice spec to a tuple of per-axis indexers."""
+    """Normalize and validate a feature-axis slice spec.
+
+    Returns a tuple of per-axis indexers; ``None`` means "no slice" and gives
+    the empty tuple. Note this makes ``numpy.newaxis`` (which *is* ``None``)
+    inexpressible, so a bare newaxis reads as "no slice"; inside a tuple it is
+    rejected explicitly.
+    """
     if feature_slice is None:
         return ()
-    if isinstance(feature_slice, tuple):
-        return feature_slice
-    return (feature_slice,)
+    indexers = feature_slice if isinstance(feature_slice, tuple) else (feature_slice,)
+    _validate_feature_slice(indexers)
+    return indexers
+
 
+def _validate_feature_slice(indexers: Tuple) -> None:
+    """Reject anything outside basic, forward indexing.
+
+    `feature_slice` exists to read part of a layer, not to reimplement NumPy
+    indexing. Keeping it to basic forward indexing is what lets every backend
+    mean the same thing by it: h5py handles slices with a positive step
+    (negative start/stop included), integers and one ``Ellipsis`` exactly as
+    NumPy does, but rejects a negative step outright. Allowing more would make
+    the legacy and HDF5 backends diverge, or require reimplementing NumPy's
+    semantics to keep them together.
+
+    Parameters
+    ----------
+    indexers : tuple
+        Per-axis indexers, as returned by :func:`_normalize_feature_slice`.
 
-def _is_fancy(indexer: object) -> bool:
-    """Whether an indexer is a fancy (list/array) index rather than a slice."""
-    if isinstance(indexer, (slice, int, np.integer)):
-        return False
-    return isinstance(indexer, (list, tuple, np.ndarray, range))
+    Raises
+    ------
+    ValueError
+        If any indexer is outside the supported set.
+    """
+    advice = (
+        " feature_slice supports basic forward indexing only (slices with a "
+        "positive step, integers and a single Ellipsis). Read without "
+        "feature_slice and index the resulting array with NumPy instead."
+    )
+    if sum(ix is Ellipsis for ix in indexers) > 1:
+        raise ValueError("feature_slice may contain at most one Ellipsis." + advice)
+    for ix in indexers:
+        # bool is a subclass of int in Python, so it must be rejected *before*
+        # the integer check or True would silently be read as index 1.
+        if isinstance(ix, (bool, np.bool_)):
+            raise ValueError(
+                "feature_slice does not support boolean indexing." + advice
+            )
+        if ix is Ellipsis or isinstance(ix, (int, np.integer)):
+            continue
+        if isinstance(ix, slice):
+            step = ix.step
+            if step is not None and step < 1:
+                raise ValueError(
+                    "feature_slice does not support a step below 1 "
+                    "(got {}).".format(step) + advice
+                )
+            continue
+        if ix is None:
+            raise ValueError(
+                "feature_slice does not support numpy.newaxis." + advice
+            )
+        raise ValueError(
+            "feature_slice does not support {} indexing.".format(
+                type(ix).__name__
+            )
+            + advice
+        )
 
 
 class FeatureStore(ABC):
@@ -157,7 +219,6 @@ def iter_chunks(
         self,
         layer: str,
         labels: Optional[Sequence[str]] = None,
-        feature_slice: FeatureSlice = None,
         axis: int = 1,
         size: Optional[int] = None,
     ) -> Iterator[Tuple[slice, np.ndarray]]:
@@ -167,18 +228,22 @@ def iter_chunks(
         of the selected array and `block` is that slab. The slice is yielded so
         that callers can place results back without tracking offsets themselves.
 
+        There is deliberately no `feature_slice` here. Composing an arbitrary
+        slice with the per-slab slice means reimplementing NumPy's index
+        arithmetic, which is both the source of subtle wrong results and more
+        generality than reading a layer in slabs needs. Slice the blocks as they
+        come out instead.
+
         Parameters
         ----------
         layer : str
             DNN layer.
         labels : sequence of str, optional
             Stimulus labels, as in :meth:`read`.
-        feature_slice : slice, int, array-like or tuple, optional
-            Index applied to the feature axes before iterating.
         axis : int, optional
-            Axis of the *selected* array to iterate over. Axis 0 is the sample
-            axis; the default, axis 1, is the outermost feature axis, which is
-            the one partial feature reads are about.
+            Axis to iterate over. Axis 0 is the sample axis; the default, axis
+            1, is the outermost feature axis, which is the one partial feature
+            reads are about.
         size : int, optional
             Number of elements per slab. Defaults to the on-disk chunk extent
             along `axis` when the backend has one, so that iteration is
@@ -190,7 +255,7 @@ def iter_chunks(
         tuple of (slice, numpy.ndarray)
             The slice applied to `axis`, and the corresponding slab.
         """
-        selected_shape = self._selected_shape(layer, labels, feature_slice)
+        selected_shape = self._selected_shape(layer, labels)
         if axis < 0:
             axis += len(selected_shape)
         if not 0 <= axis < len(selected_shape):
@@ -206,30 +271,28 @@ def iter_chunks(
         if size < 1:
             raise ValueError("size must be positive, got {}".format(size))
 
-        base = _normalize_feature_slice(feature_slice)
         for start in range(0, length, size):
             sl = slice(start, min(start + size, length))
             if axis == 0:
                 block_labels = (
                     list(self.labels) if labels is None else list(labels)
                 )[sl]
-                yield sl, self.read(layer, block_labels, feature_slice)
+                yield sl, self.read(layer, block_labels)
             else:
-                yield sl, self.read(
-                    layer, labels, _compose_slice(base, axis - 1, sl)
-                )
+                # axis >= 1 indexes the feature axes, where axis 1 of the array
+                # is entry 0 of the feature-slice tuple.
+                indexers = (slice(None),) * (axis - 1) + (sl,)
+                yield sl, self.read(layer, labels, indexers)
 
     def _selected_shape(
         self,
         layer: str,
-        labels: Optional[Sequence[str]],
-        feature_slice: FeatureSlice,
+        labels: Optional[Sequence[str]] = None,
     ) -> Tuple[int, ...]:
         """Shape the selection would have, without reading the data."""
         full = self.shape(layer)
         n_samples = full[0] if labels is None else len(labels)
-        feature_shape = _sliced_shape(full[1:], feature_slice)
-        return (n_samples, *feature_shape)
+        return (n_samples, *full[1:])
 
     def _default_iter_size(
         self, layer: str, axis: int, selected_shape: Tuple[int, ...]
@@ -247,53 +310,6 @@ def _default_iter_size(
         return max(1, min(selected_shape[axis], DEFAULT_TARGET_CHUNK_BYTES // per_element))
 
 
-def _compose_slice(base: Tuple, axis: int, sl: slice) -> Tuple:
-    """Return `base` with `sl` intersected into its entry for `axis`.
-
-    `axis` indexes the *feature* axes (axis 0 here is array axis 1), and `sl`
-    is relative to the result of applying `base`.
-    """
-    out = list(base) + [slice(None)] * max(0, axis + 1 - len(base))
-    current = out[axis]
-    if isinstance(current, slice) and current == slice(None):
-        out[axis] = sl
-    elif isinstance(current, slice):
-        # Compose the two slices by indexing a range with both.
-        out[axis] = _compose_two_slices(current, sl)
-    else:
-        out[axis] = np.asarray(current)[sl]
-    return tuple(out)
-
-
-def _compose_two_slices(outer: slice, inner: slice) -> slice:
-    """Compose ``x[outer][inner]`` into a single slice where possible."""
-    # Only the simple (positive-step) case arises from iter_chunks, which
-    # always produces contiguous forward slices.
-    ostep = 1 if outer.step is None else outer.step
-    if ostep < 0:
-        # Fall back to an explicit index array for reversed slices.
-        raise ValueError("cannot compose a reversed slice; pass an explicit size")
-    ostart = 0 if outer.start is None else outer.start
-    start = ostart + (0 if inner.start is None else inner.start) * ostep
-    stop = ostart + (inner.stop if inner.stop is not None else 0) * ostep
-    if outer.stop is not None:
-        stop = min(stop, outer.stop)
-    return slice(start, stop, ostep)
-
-
-def _sliced_shape(
-    feature_shape: Tuple[int, ...], feature_slice: FeatureSlice
-) -> Tuple[int, ...]:
-    """Shape of ``numpy.empty(feature_shape)[feature_slice]``, cheaply."""
-    indexers = _normalize_feature_slice(feature_slice)
-    if not indexers:
-        return feature_shape
-    # A broadcast dummy allocates nothing but lets NumPy do the index
-    # arithmetic, including negative steps, ellipses and newaxis.
-    dummy = np.broadcast_to(np.int8(0), feature_shape)
-    return tuple(int(n) for n in dummy[indexers].shape)
-
-
 class MatFeatureStore(FeatureStore):
     """Legacy per-stimulus ``.mat`` feature directory.
 
@@ -486,17 +502,13 @@ def read(
     def _read_rows(
         dset: h5py.Dataset, rows: Union[slice, List[int]], indexers: Tuple
     ) -> np.ndarray:
-        """Read `rows` from `dset`, pushing as much of `indexers` into h5py as it allows."""
+        """Read `rows` from `dset`, pushing `indexers` down into the h5py selection."""
         if not indexers:
             return np.asarray(dset[rows])
-        if any(_is_fancy(ix) for ix in indexers):
-            # A fancy row list plus a fancy feature index would be two fancy
-            # indices in one selection, which h5py rejects. Read the feature
-            # axes whole and apply the fancy part in NumPy.
-            data = np.asarray(dset[rows])
-            return data[(slice(None), *indexers)]
-        # Plain slices are not fancy, so they ride along in the same selection
-        # and h5py reads only the chunks they cover.
+        # Basic forward indexers are not fancy, so they ride along in the same
+        # selection as the row list and h5py reads only the chunks they cover.
+        # _validate_feature_slice has already rejected anything else, which is
+        # what keeps this to a single line.
         return np.asarray(dset[(rows, *indexers)])
 
     def _row_indices(self, labels: Sequence[str]) -> np.ndarray:
@@ -523,8 +535,7 @@ def _validate_all(self) -> None:
         reference_layer = ""
         for layer in self._layers:
             with self._open(layer) as f:
-                _validate_format(f, self.path(layer))
-                labels = _decode_labels(f[LABELS_DATASET][()])
+                labels = _validate_format(f, self.path(layer))
             if reference is None:
                 reference, reference_layer = labels, layer
             elif labels != reference:
@@ -539,8 +550,36 @@ def _validate_all(self) -> None:
         self._label_index = {label: i for i, label in enumerate(reference)}
 
 
-def _validate_format(f: h5py.File, path: str) -> None:
-    """Check that `f` is bdpy feature storage of a version we can read."""
+def _validate_format(f: h5py.File, path: str) -> List[str]:
+    """Check that `f` is well-formed bdpy feature storage, and return its labels.
+
+    Everything a reader relies on is checked in one place: the format marker,
+    a version in the range this build understands, both datasets present with
+    the expected rank, one label per feature row, and no duplicate labels.
+
+    Duplicate labels matter more than they look. Labels are mapped to row
+    indices once, so a repeated label would silently resolve every occurrence
+    to the last row and drop the earlier one. The legacy ``.mat`` layout cannot
+    express duplicates at all -- the file name *is* the label -- so rejecting
+    them keeps the two layouts equivalent.
+
+    Parameters
+    ----------
+    f : h5py.File
+        Open feature file.
+    path : str
+        Its path, for error messages.
+
+    Returns
+    -------
+    list of str
+        The file's labels, already decoded.
+
+    Raises
+    ------
+    RuntimeError
+        If any of the above does not hold.
+    """
     fmt = f.attrs.get(FORMAT_ATTR)
     if isinstance(fmt, bytes):
         fmt = fmt.decode("utf-8")
@@ -550,7 +589,26 @@ def _validate_format(f: h5py.File, path: str) -> None:
                 path, FORMAT_ATTR, fmt, FORMAT_NAME
             )
         )
-    version = int(f.attrs.get(FORMAT_VERSION_ATTR, 0))
+
+    if FORMAT_VERSION_ATTR not in f.attrs:
+        raise RuntimeError(
+            "{} has no {} attribute, so it is not valid bdpy feature "
+            "storage.".format(path, FORMAT_VERSION_ATTR)
+        )
+    raw_version = f.attrs[FORMAT_VERSION_ATTR]
+    try:
+        version = int(raw_version)
+    except (TypeError, ValueError):
+        raise RuntimeError(
+            "{} has a malformed {} ({!r}); expected an integer.".format(
+                path, FORMAT_VERSION_ATTR, raw_version
+            )
+        ) from None
+    if version < 1:
+        raise RuntimeError(
+            "{} declares feature storage format version {}; versions start at "
+            "1.".format(path, version)
+        )
     if version > SUPPORTED_FORMAT_VERSION:
         raise RuntimeError(
             "{} uses feature storage format version {}, but this version of "
@@ -558,10 +616,50 @@ def _validate_format(f: h5py.File, path: str) -> None:
                 path, version, SUPPORTED_FORMAT_VERSION
             )
         )
+
     for name in (FEATURES_DATASET, LABELS_DATASET):
         if name not in f:
             raise RuntimeError("{} has no /{} dataset".format(path, name))
 
+    features, labels_dset = f[FEATURES_DATASET], f[LABELS_DATASET]
+    if features.ndim < 2:
+        raise RuntimeError(
+            "{}: /{} must have a sample axis and at least one feature axis, "
+            "got shape {}".format(path, FEATURES_DATASET, features.shape)
+        )
+    if labels_dset.ndim != 1:
+        raise RuntimeError(
+            "{}: /{} must be one-dimensional, got shape {}".format(
+                path, LABELS_DATASET, labels_dset.shape
+            )
+        )
+
+    labels = _decode_labels(labels_dset[()])
+    if len(labels) != features.shape[0]:
+        raise RuntimeError(
+            "{}: /{} has {} rows but /{} has {} entries".format(
+                path, FEATURES_DATASET, features.shape[0],
+                LABELS_DATASET, len(labels),
+            )
+        )
+    duplicates = _duplicates(labels)
+    if duplicates:
+        raise RuntimeError(
+            "{}: /{} contains duplicate labels ({}). Labels identify rows, so "
+            "they must be unique.".format(
+                path, LABELS_DATASET, ", ".join(repr(d) for d in duplicates[:3])
+            )
+        )
+    return labels
+
+
+def _duplicates(labels: Sequence[str]) -> List[str]:
+    """Labels appearing more than once, in first-seen order."""
+    seen: Dict[str, int] = {}
+    for label in labels:
+        seen[label] = seen.get(label, 0) + 1
+    return [label for label, n in seen.items() if n > 1]
+
 
 def _label_mismatch_message(
     dpath: str,
@@ -606,12 +704,17 @@ def _decode_labels(raw: np.ndarray) -> List[str]:
 
 
 def detect_format(dpath: str, ext: str = "mat") -> str:
-    """Guess which storage format `dpath` holds.
+    """Work out which storage layout `dpath` holds.
 
-    Returns ``'hdf5'`` when the directory holds ``<layer>.h5`` files and no
-    layer subdirectories, and ``'mat'`` otherwise. A directory holding both is
-    read as ``'mat'``, the historical layout; pass an explicit format to
-    override.
+    The two layouts are told apart by the files that are actually there, not by
+    the mere presence of a subdirectory: a legacy tree is one whose
+    subdirectories really contain ``*.<ext>`` feature files, so an unrelated
+    subdirectory next to ``<layer>.h5`` files does not turn a chunked directory
+    into a legacy one.
+
+    A directory holding both is ambiguous and raises rather than picking one
+    silently -- guessing wrong means reading different data than the caller
+    meant. Pass an explicit `format` to resolve it.
 
     Parameters
     ----------
@@ -624,12 +727,31 @@ def detect_format(dpath: str, ext: str = "mat") -> str:
     -------
     str
         ``'hdf5'`` or ``'mat'``.
+
+    Raises
+    ------
+    RuntimeError
+        If the directory holds both layouts, or neither.
     """
     escaped = glob.escape(dpath)
-    has_subdirs = any(
-        os.path.isdir(os.path.join(dpath, d)) for d in os.listdir(dpath)
-    )
     has_h5 = bool(glob.glob(os.path.join(escaped, "*." + HDF5_EXT)))
-    if has_h5 and not has_subdirs:
+    has_legacy = any(
+        glob.glob(os.path.join(glob.escape(os.path.join(dpath, d)), "*." + ext))
+        for d in os.listdir(dpath)
+        if os.path.isdir(os.path.join(dpath, d))
+    )
+
+    if has_h5 and has_legacy:
+        raise RuntimeError(
+            "{} holds both chunked .{} files and a legacy .{} tree, so the "
+            "storage format is ambiguous. Pass format='hdf5' or format='mat' "
+            "to choose one.".format(dpath, HDF5_EXT, ext)
+        )
+    if has_h5:
         return "hdf5"
-    return "mat"
+    if has_legacy:
+        return "mat"
+    raise RuntimeError(
+        "No features found in {}: expected either <layer>.{} files or "
+        "<layer>/<label>.{} subdirectories.".format(dpath, HDF5_EXT, ext)
+    )
diff --git a/bdpy/dataform/feature_hdf5.py b/bdpy/dataform/feature_hdf5.py
index ef39795..3a8ec58 100644
--- a/bdpy/dataform/feature_hdf5.py
+++ b/bdpy/dataform/feature_hdf5.py
@@ -14,7 +14,6 @@
         attrs:
             bdpy_format         = "features"
             bdpy_format_version = 1
-            layer               = "<layer name>"
 
 One file per layer rather than one file for everything: layers differ in shape
 and in optimal chunk shape, and separate files keep regeneration, copying and
@@ -48,7 +47,7 @@
 import uuid
 import weakref
 from types import TracebackType
-from typing import Iterable, Optional, Sequence, Tuple, Type
+from typing import Iterable, Optional, Sequence, Set, Tuple, Type
 
 import h5py
 import numpy as np
@@ -63,6 +62,7 @@
     LABELS_DATASET,
     SUPPORTED_FORMAT_VERSION,
     MatFeatureStore,
+    _duplicates,
 )
 
 __all__ = [
@@ -71,8 +71,8 @@
     "save_features",
 ]
 
-#: Number of samples the resizable dataset grows by at a time.
-_GROW_BLOCK = 64
+#: Chunk length of the resizable labels dataset.
+_LABELS_CHUNK = 64
 
 
 def _string_dtype() -> np.dtype:
@@ -80,11 +80,10 @@ def _string_dtype() -> np.dtype:
     return dtype
 
 
-def _write_header(f: h5py.File, layer: Optional[str]) -> None:
+def _write_header(f: h5py.File) -> None:
+    """Stamp the format marker and version. The layer name is the file name."""
     f.attrs[FORMAT_ATTR] = FORMAT_NAME
     f.attrs[FORMAT_VERSION_ATTR] = SUPPORTED_FORMAT_VERSION
-    if layer is not None:
-        f.attrs["layer"] = layer
 
 
 #: Suffix of the scratch file a write is staged in.
@@ -147,7 +146,6 @@ def save_features(
     path: str,
     features: np.ndarray,
     labels: Sequence[str],
-    layer: Optional[str] = None,
     chunks: Optional[Tuple[int, ...]] = None,
     target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
     compression: Optional[str] = None,
@@ -165,8 +163,6 @@ def save_features(
         Feature array of shape ``(n_samples, *feature_shape)``.
     labels : sequence of str
         One stimulus label per sample, in the same order as `features`.
-    layer : str, optional
-        Layer name recorded in the file. Defaults to the file's basename.
     chunks : tuple of int, optional
         Explicit chunk shape. Defaults to :func:`choose_chunk_shape`.
     target_chunk_bytes : int, optional
@@ -200,29 +196,20 @@ def save_features(
             "got {} labels for {} samples".format(len(labels), features.shape[0])
         )
 
-    dtype = np.dtype(features.dtype if dtype is None else dtype)
-    if chunks is None:
-        chunks = choose_chunk_shape(
-            features.shape, dtype, target_bytes=target_chunk_bytes
-        )
-    if layer is None:
-        layer = os.path.splitext(os.path.basename(path))[0]
-
-    tmp_path = _prepare_target(path, overwrite)
-    try:
-        with h5py.File(tmp_path, "x") as f:
-            _write_header(f, layer)
-            f.create_dataset(
-                FEATURES_DATASET,
-                data=features.astype(dtype, copy=False),
-                chunks=chunks,
-                compression=compression,
-            )
-            f.create_dataset(LABELS_DATASET, data=labels, dtype=_string_dtype())
-    except BaseException:
-        _discard(tmp_path)
-        raise
-    _promote(tmp_path, path)
+    # One write path, not two: writing a whole layer is the degenerate case of
+    # writing it incrementally, and keeping them separate is what let a bug
+    # (a zero-sample layer) land in only one of them.
+    with FeatureWriter(
+        path,
+        feature_shape=features.shape[1:],
+        dtype=features.dtype if dtype is None else dtype,
+        n_samples=features.shape[0],
+        chunks=chunks,
+        target_chunk_bytes=target_chunk_bytes,
+        compression=compression,
+        overwrite=overwrite,
+    ) as writer:
+        writer.extend(features, labels)
 
 
 class FeatureWriter:
@@ -246,8 +233,6 @@ class FeatureWriter:
         Shape of a single sample's features, without the sample axis.
     dtype : numpy.dtype
         Dtype to store.
-    layer : str, optional
-        Layer name recorded in the file. Defaults to the file's basename.
     n_samples : int, optional
         Expected number of samples, if known. Used only to pick a chunk shape.
     chunks : tuple of int, optional
@@ -275,7 +260,6 @@ def __init__(
         path: str,
         feature_shape: Sequence[int],
         dtype: np.dtype,
-        layer: Optional[str] = None,
         n_samples: Optional[int] = None,
         chunks: Optional[Tuple[int, ...]] = None,
         target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
@@ -287,6 +271,10 @@ def __init__(
             raise ValueError("feature_shape must have at least one axis")
         self._dtype = np.dtype(dtype)
         self._n = 0
+        # Labels identify rows, so they must be unique across the whole file --
+        # not just within one call. A repeated label would make every lookup
+        # resolve to the last row and silently drop the earlier one.
+        self._seen: Set[str] = set()
 
         if chunks is None:
             chunks = choose_chunk_shape(
@@ -295,9 +283,6 @@ def __init__(
                 target_bytes=target_chunk_bytes,
                 n_samples_known=n_samples is not None,
             )
-        if layer is None:
-            layer = os.path.splitext(os.path.basename(path))[0]
-
         # Fail before doing any work if the target is occupied, then write to
         # a scratch file so `path` stays untouched until close() succeeds.
         self._path = path
@@ -310,7 +295,7 @@ def __init__(
         # closes over the path instead of the writer, so it can never resurrect
         # the object or outlive a name it no longer owns.
         self._finalizer = weakref.finalize(self, _discard, self._tmp_path)
-        _write_header(self._file, layer)
+        _write_header(self._file)
         self._features = self._file.create_dataset(
             FEATURES_DATASET,
             shape=(0, *self._feature_shape),
@@ -324,7 +309,7 @@ def __init__(
             shape=(0,),
             maxshape=(None,),
             dtype=_string_dtype(),
-            chunks=(max(1, _GROW_BLOCK),),
+            chunks=(_LABELS_CHUNK,),
         )
 
     @property
@@ -379,6 +364,9 @@ def extend(self, features: np.ndarray, labels: Sequence[str]) -> None:
             raise ValueError(
                 "got {} labels for {} samples".format(len(labels), features.shape[0])
             )
+        # Check before writing anything, so a rejected batch leaves the file
+        # exactly as it was.
+        self._reject_duplicates(labels)
         if not labels:
             return
 
@@ -387,8 +375,28 @@ def extend(self, features: np.ndarray, labels: Sequence[str]) -> None:
         self._labels.resize(new_n, axis=0)
         self._features[self._n:new_n] = features.astype(self._dtype, copy=False)
         self._labels[self._n:new_n] = labels
+        self._seen.update(labels)
         self._n = new_n
 
+    def _reject_duplicates(self, labels: Sequence[str]) -> None:
+        """Refuse labels repeated within this batch or already written."""
+        within = _duplicates(labels)
+        if within:
+            raise ValueError(
+                "duplicate labels in this batch: {}. Labels identify rows, so "
+                "they must be unique.".format(
+                    ", ".join(repr(d) for d in within[:3])
+                )
+            )
+        already = [label for label in labels if label in self._seen]
+        if already:
+            raise ValueError(
+                "labels already written to {}: {}. Labels identify rows, so "
+                "they must be unique.".format(
+                    self._path, ", ".join(repr(d) for d in already[:3])
+                )
+            )
+
     def close(self) -> None:
         """Finish the file and move it into place. Idempotent.
 
@@ -528,7 +536,6 @@ def convert_features_to_hdf5(
             out_path,
             feature_shape=full_shape[1:],
             dtype=store.dtype(layer),
-            layer=layer,
             n_samples=full_shape[0],
             target_chunk_bytes=target_chunk_bytes,
             compression=compression,
diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py
index 7b8d049..c5ec3f0 100644
--- a/bdpy/dataform/features.py
+++ b/bdpy/dataform/features.py
@@ -26,6 +26,8 @@
     FeatureStore,
     HDF5FeatureStore,
     MatFeatureStore,
+    _determine_num_parallel,
+    _load_array_with_key,
     detect_format,
 )
 
@@ -41,27 +43,6 @@
 )
 
 
-def _load_array_with_key(key: str, path: str) -> np.ndarray:
-    # v5 .mat via scipy, v7.3 (HDF5) via h5py; avoids hdf5storage on the load
-    # path, which breaks under NumPy 2.0 (see bdpy/dataform/_mat_v73.py).
-    return _mat_v73.loadmat_key(path, key)
-
-
-def _determine_num_parallel(num_files: int) -> int:
-    # NOTE: optimal number of parallel processes is not clear. It could depend
-    # on several factors such as the number of files, the size of files, the
-    # number of cores, etc. For now, we use a simple heuristic based on the
-    # number of files.
-    num_parallel: int
-    if num_files < 16:
-        num_parallel = 1
-    elif num_files < 64:
-        num_parallel = 16
-    else:
-        num_parallel = 64
-    return num_parallel
-
-
 class Features(object):
     """DNN features class.
 
@@ -232,14 +213,18 @@ def get(
 
     def iter_chunks(
             self, layer: str, label: Union[str, List[str], None] = None,
-            feature_slice: FeatureSlice = None, axis: int = 1,
-            size: Optional[int] = None
+            axis: int = 1, size: Optional[int] = None
         ) -> Iterator[Tuple[slice, np.ndarray]]:
         """Iterate over `layer` in slabs along `axis`.
 
-        This is the streaming counterpart of :meth:`get`: it never holds more
-        than one slab in memory, so a layer far larger than RAM can be processed
-        end to end.
+        This is the streaming counterpart of :meth:`get`: with chunked HDF5
+        storage it holds only one slab at a time, so a layer far larger than RAM
+        can be processed end to end.
+
+        There is no `feature_slice` here on purpose -- slice the blocks as they
+        come out instead. Composing an arbitrary slice with the per-slab slice
+        means reimplementing NumPy's index arithmetic, which is more generality
+        than reading a layer in slabs needs.
 
         Parameters
         ----------
@@ -247,8 +232,6 @@ def iter_chunks(
             DNN layer
         label: str or list, optional
             Sample label(s). ``None`` iterates over every label.
-        feature_slice: slice, int, array-like or tuple, optional
-            Index applied to the feature axes before iterating.
         axis: int
             Axis of the selected array to iterate over. Axis 0 is the sample
             axis; the default, axis 1, is the outermost feature axis.
@@ -256,6 +239,13 @@ def iter_chunks(
             Elements per slab. Defaults to the on-disk chunk extent along
             `axis`, so that each element is read exactly once.
 
+        Notes
+        -----
+        The one-slab-at-a-time guarantee holds when the requested labels all
+        live in a single feature directory, which is the usual case. When they
+        are spread across several `dpath` entries, no single store can stream
+        them, so the selection is read in full and then split.
+
         Yields
         ------
         tuple of (slice, numpy.ndarray)
@@ -287,14 +277,12 @@ def iter_chunks(
 
         store = self.__store_for(labels)
         if store is not None:
-            yield from store.iter_chunks(
-                layer, labels, feature_slice, axis=axis, size=size
-            )
+            yield from store.iter_chunks(layer, labels, axis=axis, size=size)
             return
 
         # Labels span several directories, so no single store can stream them.
         # Fall back to slicing a full read, which still yields the same blocks.
-        features = self.__read(layer, labels, feature_slice)
+        features = self.__read(layer, labels)
         length = features.shape[axis]
         if size is None:
             size = length
diff --git a/docs/dataform_features.md b/docs/dataform_features.md
index 9d8486e..6e1b780 100644
--- a/docs/dataform_features.md
+++ b/docs/dataform_features.md
@@ -78,8 +78,10 @@ features = Features('/path/to/features')   # either layout
 feat = features.get(layer='conv5')
 ```
 
-Pass `format='mat'` or `format='hdf5'` to skip detection. A directory holding
-both layouts is read as `.mat` unless you say otherwise.
+Detection keys off the files that are actually present, so an unrelated
+subdirectory next to `<layer>.h5` files does not make a directory look legacy.
+A directory holding *both* layouts is ambiguous and raises rather than picking
+one silently; pass `format='mat'` or `format='hdf5'` to resolve it.
 
 ### Partial reads
 
@@ -104,6 +106,20 @@ feat = features.get(
 # The index addresses the flattened full feature space, so applying it to an
 # already-sliced array would select the wrong units; the combination raises
 # ValueError instead.
+```
+
+`feature_slice` is **basic forward indexing**: slices with a positive step
+(negative `start`/`stop` are fine), integers, a single `Ellipsis`, and tuples of
+those. Anything else -- a negative step, fancy indexing with a list or array,
+booleans, `np.newaxis` -- raises `ValueError` on *both* backends. The restriction
+is what lets the two layouts mean the same thing by the same index; read without
+`feature_slice` and index the result with NumPy when you need more.
+
+``` python
+features.get('conv5', feature_slice=np.s_[128:256])     # ok
+features.get('conv5', feature_slice=np.s_[8:16, 1:4])   # ok
+features.get('conv5', feature_slice=np.s_[::-1])        # ValueError
+features.get('conv5', feature_slice=np.s_[[3, 1, 7]])   # ValueError
 
 # Full shape without reading anything
 n_stimuli, *feature_shape = features.shape('conv5')
@@ -119,7 +135,11 @@ for sl, block in features.iter_chunks('conv5', axis=1):
     out[:, sl] = transform(block)
 ```
 
-`axis=0` iterates over stimuli instead of features.
+`axis=0` iterates over stimuli instead of features. `iter_chunks` takes no
+`feature_slice` -- slice each `block` as it comes out instead. The
+one-slab-at-a-time guarantee holds when the requested labels all live in one
+feature directory; when they are spread across several `dpath` entries no single
+store can stream them, so the selection is read in full and then split.
 
 ### Writing
 
@@ -164,9 +184,10 @@ layout, and streams in batches, so a layer is never held in memory in full.
 Existing `<layer>.h5` files are skipped unless `overwrite=True`, and a layer that
 fails to convert leaves no file behind, so re-running picks up where it stopped.
 
-All layers in a directory must hold the same stimulus labels in the same order;
-a mismatch is rejected when the directory is opened, as it is for the legacy
-layout.
+All layers in a directory must hold the same stimulus labels in the same order,
+and labels must be unique within a layer; both are rejected when the directory
+is opened, as they are for the legacy layout (where the file name *is* the
+label, so duplicates cannot arise). Writing a duplicate label raises too.
 
 ### Chunk shape
 
diff --git a/tests/dataform/test_feature_chunking.py b/tests/dataform/test_feature_chunking.py
index 2ad6719..d46fec1 100644
--- a/tests/dataform/test_feature_chunking.py
+++ b/tests/dataform/test_feature_chunking.py
@@ -32,52 +32,29 @@ def test_invariants(self):
                 chunk = choose_chunk_shape(shape, np.dtype(dtype))
                 with self.subTest(shape=shape, dtype=dtype):
                     self.assertEqual(len(chunk), len(shape))
-                    # Never zero, never larger than the dataset.
+                    # Never zero; never larger than the dataset on axes that
+                    # actually have a size. A zero-sample layer still needs a
+                    # chunk extent of at least 1 -- HDF5 forbids a zero dim.
                     for c, s in zip(chunk, shape):
                         self.assertGreaterEqual(c, 1)
-                        self.assertLessEqual(c, s)
-                    # Trailing axes are kept whole.
-                    self.assertEqual(chunk[2:], tuple(shape[2:]))
-
-    def test_chunk_stays_within_budget(self):
+                        if s > 0:
+                            self.assertLessEqual(c, s)
+
+    def test_chunk_stays_within_budget_when_the_budget_is_reachable(self):
+        # The budget is an upper bound only when some chunk can meet it. The
+        # current policy keeps the trailing axes whole, so a shape whose
+        # trailing axes alone exceed the budget cannot: assert only where the
+        # smallest chunk that policy can produce does fit.
         dtype = np.dtype(np.float32)
         for shape in self.shapes:
             chunk = choose_chunk_shape(shape, dtype)
+            smallest = _nbytes((1, 1, *shape[2:]), dtype)
             with self.subTest(shape=shape):
-                if _nbytes(shape, dtype) <= DEFAULT_TARGET_CHUNK_BYTES:
-                    # Small enough to store as a single chunk.
-                    self.assertEqual(chunk, tuple(shape))
-                else:
+                if smallest <= DEFAULT_TARGET_CHUNK_BYTES:
                     self.assertLessEqual(
                         _nbytes(chunk, dtype), DEFAULT_TARGET_CHUNK_BYTES
                     )
 
-    def test_budget_is_actually_used(self):
-        # A chunk far below budget means needless read amplification: every
-        # oversized array should fill at least a quarter of the budget.
-        dtype = np.dtype(np.float32)
-        for shape in self.shapes:
-            if _nbytes(shape, dtype) <= DEFAULT_TARGET_CHUNK_BYTES:
-                continue
-            chunk = choose_chunk_shape(shape, dtype)
-            with self.subTest(shape=shape):
-                self.assertGreater(
-                    _nbytes(chunk, dtype), DEFAULT_TARGET_CHUNK_BYTES // 4
-                )
-
-    def test_both_sliceable_axes_are_chunked(self):
-        # The point of #144: a channel slice must not have to read every channel.
-        shape = (1200, 256, 13, 13)
-        chunk = choose_chunk_shape(shape, np.dtype(np.float32))
-        self.assertLess(chunk[0], shape[0])
-        self.assertLess(chunk[1], shape[1])
-
-    def test_whole_array_below_budget_is_one_chunk(self):
-        shape = (50, 1000)
-        self.assertEqual(
-            choose_chunk_shape(shape, np.dtype(np.float32)), shape
-        )
-
     def test_single_sample(self):
         chunk = choose_chunk_shape((1, 4096, 7, 7), np.dtype(np.float32))
         self.assertEqual(chunk[0], 1)
@@ -99,18 +76,6 @@ def test_unknown_sample_count_is_not_capped(self):
         )
         self.assertGreater(chunk[0], 1)
 
-    def test_edge_padding_is_small(self):
-        # HDF5 allocates whole chunks, so uneven extents pad the edges and
-        # inflate the file. Keep that overhead under 10%.
-        dtype = np.dtype(np.float32)
-        for shape in self.shapes:
-            chunk = choose_chunk_shape(shape, dtype)
-            allocated = 1
-            for s, c in zip(shape, chunk):
-                allocated *= -(-s // c) * c
-            with self.subTest(shape=shape):
-                self.assertLessEqual(allocated, int(np.prod(shape) * 1.1) + 1)
-
     def test_deterministic(self):
         shape = (1200, 256, 13, 13)
         dtype = np.dtype(np.float32)
@@ -118,6 +83,22 @@ def test_deterministic(self):
             choose_chunk_shape(shape, dtype), choose_chunk_shape(shape, dtype)
         )
 
+    def test_zero_sample_layer_still_gets_a_usable_chunk(self):
+        # HDF5 forbids a zero chunk dimension, so an empty layer needs >= 1.
+        chunk = choose_chunk_shape((0, 8), np.dtype(np.float32))
+        self.assertEqual(chunk[0], 1)
+        self.assertGreaterEqual(chunk[1], 1)
+
+    def test_rejects_empty_feature_axis(self):
+        for bad in ((5, 0), (5, 3, 0)):
+            with self.subTest(shape=bad):
+                with self.assertRaises(ValueError):
+                    choose_chunk_shape(bad, np.dtype(np.float32))
+
+    def test_rejects_negative_sample_axis(self):
+        with self.assertRaises(ValueError):
+            choose_chunk_shape((-1, 4), np.dtype(np.float32))
+
     def test_rejects_bad_input(self):
         with self.assertRaises(ValueError):
             choose_chunk_shape((100,), np.dtype(np.float32))
diff --git a/tests/dataform/test_feature_hdf5.py b/tests/dataform/test_feature_hdf5.py
index b93845e..d6827d7 100644
--- a/tests/dataform/test_feature_hdf5.py
+++ b/tests/dataform/test_feature_hdf5.py
@@ -13,6 +13,7 @@
 from bdpy.dataform import Features
 from bdpy.dataform._feature_store import (
     MatFeatureStore,
+    SUPPORTED_FORMAT_VERSION,
     FORMAT_ATTR,
     FORMAT_NAME,
     FORMAT_VERSION_ATTR,
@@ -52,7 +53,6 @@ def test_writes_schema_v1(self):
             self.assertEqual(
                 int(f.attrs[FORMAT_VERSION_ATTR]), SUPPORTED_FORMAT_VERSION
             )
-            self.assertEqual(f.attrs['layer'], 'conv5')
             assert_array_equal(f['features'][()], self.data)
             self.assertEqual(
                 [s.decode('utf-8') for s in f['labels'][()]], self.labels
@@ -134,11 +134,10 @@ def test_extend(self):
         )
 
     def test_writes_valid_header(self):
-        with FeatureWriter(self.path, (16, 3, 3), np.float32, layer='conv5') as writer:
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
             writer.append(self.data[0], self.labels[0])
         with h5py.File(self.path, 'r') as f:
             self.assertEqual(f.attrs[FORMAT_ATTR], FORMAT_NAME)
-            self.assertEqual(f.attrs['layer'], 'conv5')
             self.assertIsNotNone(f['features'].chunks)
 
     def test_rejects_wrong_shape(self):
@@ -215,6 +214,154 @@ def test_small_batches_produce_the_same_file(self):
         assert_array_equal(Features(self.h5dir).get('conv5'), self.stacked['conv5'])
 
 
+class TestSchemaValidation(unittest.TestCase):
+    """Everything a reader depends on is checked when the file is opened."""
+
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+        save_features(self.path, np.zeros((3, 4)), ['a', 'b', 'c'])
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def _mutate(self, fn):
+        with h5py.File(self.path, 'a') as f:
+            fn(f)
+
+    def test_missing_version_is_rejected(self):
+        # Previously read as version 0 and accepted, contradicting the docs.
+        self._mutate(lambda f: f.attrs.__delitem__(FORMAT_VERSION_ATTR))
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+    def test_malformed_version_is_rejected(self):
+        self._mutate(lambda f: f.attrs.__setitem__(FORMAT_VERSION_ATTR, 'v1'))
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+    def test_version_below_one_is_rejected(self):
+        for bad in (0, -1):
+            with self.subTest(version=bad):
+                self._mutate(lambda f, b=bad: f.attrs.__setitem__(FORMAT_VERSION_ATTR, b))
+                with self.assertRaises(RuntimeError):
+                    HDF5FeatureStore(self.tmpdir.name)
+
+    def test_future_version_is_rejected(self):
+        self._mutate(
+            lambda f: f.attrs.__setitem__(
+                FORMAT_VERSION_ATTR, SUPPORTED_FORMAT_VERSION + 1
+            )
+        )
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('upgrade bdpy', str(ctx.exception))
+
+    def test_row_count_mismatch_is_rejected(self):
+        def shrink(f):
+            del f['labels']
+            f.create_dataset('labels', data=['a', 'b'],
+                             dtype=h5py.string_dtype(encoding='utf-8'))
+        self._mutate(shrink)
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('rows', str(ctx.exception))
+
+    def test_wrong_ndim_is_rejected(self):
+        def flatten_features(f):
+            del f['features']
+            f.create_dataset('features', data=np.zeros(3))
+        self._mutate(flatten_features)
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+        save_features(self.path, np.zeros((2, 4)), ['a', 'b'], overwrite=True)
+
+        def widen_labels(f):
+            del f['labels']
+            f.create_dataset('labels', data=[['a'], ['b']],
+                             dtype=h5py.string_dtype(encoding='utf-8'))
+        self._mutate(widen_labels)
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+    def test_duplicate_labels_are_rejected_on_read(self):
+        # Labels map to row indices, so a repeat would resolve both rows to the
+        # last one and silently drop the first. The .mat layout cannot express
+        # duplicates at all, so neither should this one.
+        def duplicate(f):
+            del f['labels']
+            f.create_dataset('labels', data=['a', 'a', 'c'],
+                             dtype=h5py.string_dtype(encoding='utf-8'))
+        self._mutate(duplicate)
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('duplicate', str(ctx.exception))
+
+
+class TestDuplicateLabelsOnWrite(unittest.TestCase):
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_save_features_rejects_duplicates(self):
+        with self.assertRaises(ValueError):
+            save_features(self.path, np.zeros((2, 4)), ['a', 'a'])
+        self.assertFalse(os.path.exists(self.path))
+
+    def test_writer_rejects_duplicates_within_one_call(self):
+        with FeatureWriter(self.path, (4,), np.float32) as writer:
+            with self.assertRaises(ValueError) as ctx:
+                writer.extend(np.zeros((2, 4)), ['a', 'a'])
+            self.assertIn('this batch', str(ctx.exception))
+            writer.abort()
+
+    def test_writer_rejects_duplicates_across_calls(self):
+        with FeatureWriter(self.path, (4,), np.float32) as writer:
+            writer.append(np.zeros(4), 'a')
+            with self.assertRaises(ValueError) as ctx:
+                writer.append(np.ones(4), 'a')
+            self.assertIn('already written', str(ctx.exception))
+            writer.abort()
+
+    def test_rejected_batch_leaves_earlier_rows_intact(self):
+        # The check runs before writing, so a refused batch is a no-op.
+        with FeatureWriter(self.path, (4,), np.float32) as writer:
+            writer.extend(np.ones((2, 4)), ['a', 'b'])
+            with self.assertRaises(ValueError):
+                writer.extend(np.zeros((2, 4)), ['c', 'a'])
+            self.assertEqual(writer.n_samples, 2)
+            writer.extend(np.full((1, 4), 7.0), ['c'])
+        store = HDF5FeatureStore(self.tmpdir.name)
+        self.assertEqual(store.labels, ['a', 'b', 'c'])
+
+
+class TestZeroSampleLayer(unittest.TestCase):
+    """Both write paths agree on an empty layer."""
+
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_save_features_accepts_zero_samples(self):
+        path = os.path.join(self.tmpdir.name, 'conv5.h5')
+        save_features(path, np.empty((0, 8), dtype=np.float32), [])
+        store = HDF5FeatureStore(self.tmpdir.name)
+        self.assertEqual(store.labels, [])
+        self.assertEqual(store.read('conv5').shape, (0, 8))
+
+    def test_writer_accepts_zero_samples(self):
+        path = os.path.join(self.tmpdir.name, 'fc8.h5')
+        with FeatureWriter(path, (8,), np.float32):
+            pass
+        self.assertEqual(HDF5FeatureStore(self.tmpdir.name).read('fc8').shape, (0, 8))
+
+
 class TestAtomicWriteAndOverwrite(unittest.TestCase):
     """A failed write must leave nothing behind, and must never clobber.
 
diff --git a/tests/dataform/test_feature_store.py b/tests/dataform/test_feature_store.py
index e82e44c..a70d9c5 100644
--- a/tests/dataform/test_feature_store.py
+++ b/tests/dataform/test_feature_store.py
@@ -39,7 +39,6 @@
     ('slice with step', 'conv5', None, np.s_[::2]),
     ('multi-axis slice', 'conv5', None, np.s_[8:16, 1:4, :]),
     ('integer index', 'conv5', None, np.s_[5]),
-    ('fancy index', 'conv5', None, np.s_[[3, 1, 7]]),
     ('ellipsis', 'conv5', None, np.s_[..., 1:3]),
     ('labels and slice', 'conv5', ['img0009', 'img0001', 'img0009'], np.s_[8:16]),
     ('2d layer', 'fc8', None, np.s_[10:40]),
@@ -139,12 +138,13 @@ def test_slices_cover_the_axis_without_overlap(self):
             covered.extend(range(*sl.indices(length)))
         self.assertEqual(covered, list(range(length)))
 
-    def test_composes_with_a_feature_slice(self):
-        expected = self.stacked['conv5'][:, 4:20]
-        blocks = list(
-            self.from_h5.iter_chunks('conv5', feature_slice=np.s_[4:20], axis=1, size=6)
-        )
-        assert_array_equal(np.concatenate([b for _, b in blocks], axis=1), expected)
+    def test_blocks_can_be_sliced_by_the_caller(self):
+        # iter_chunks deliberately takes no feature_slice; slicing the blocks as
+        # they come out is the supported way to combine the two.
+        wanted = np.zeros_like(self.stacked['conv5'][:, :, 1:2])
+        for sl, block in self.from_h5.iter_chunks('conv5', axis=1, size=6):
+            wanted[:, sl] = block[:, :, 1:2]
+        assert_array_equal(wanted, self.stacked['conv5'][:, :, 1:2])
 
     def test_restricts_to_labels(self):
         labels = ['img0009', 'img0001']
@@ -219,6 +219,76 @@ def test_dataset_is_chunked_on_both_sliceable_axes(self):
         self.assertEqual(chunks[2:], shape[2:])  # spatial axes kept whole
 
 
+class TestFeatureSliceValidation(_BackendPair):
+    """feature_slice is basic forward indexing, identically on both backends.
+
+    The point of narrowing it: whatever is accepted must mean the same thing on
+    both layouts, and whatever is rejected must be rejected by both. Previously
+    a negative step worked on .mat and raised inside h5py, which is exactly the
+    divergence this table guards against.
+    """
+
+    accepted = [
+        ('plain slice', np.s_[4:16]),
+        ('negative stop', np.s_[:-1]),
+        ('negative start', np.s_[-8:]),
+        ('positive step', np.s_[::2]),
+        ('integer', np.s_[5]),
+        ('ellipsis', np.s_[..., 1:3]),
+        ('multi-axis', np.s_[4:16, 1:3]),
+        ('open slice', np.s_[:]),
+    ]
+
+    rejected = [
+        ('negative step', np.s_[::-1]),
+        ('reversed with bounds', np.s_[16:4:-1]),
+        ('zero step', slice(None, None, 0)),
+        ('list', np.s_[[3, 1, 7]]),
+        ('ndarray', np.array([1, 2])),
+        ('range', range(3)),
+        ('bool', True),
+        ('numpy bool', np.True_),
+        ('newaxis in tuple', (np.newaxis, slice(None))),
+        ('two ellipses', (Ellipsis, Ellipsis)),
+    ]
+
+    def test_accepted_agree_across_backends_and_numpy(self):
+        for name, spec in self.accepted:
+            with self.subTest(index=name):
+                from_mat = self.from_mat.get('conv5', feature_slice=spec)
+                from_h5 = self.from_h5.get('conv5', feature_slice=spec)
+                indexers = spec if isinstance(spec, tuple) else (spec,)
+                expected = self.stacked['conv5'][(slice(None), *indexers)]
+                assert_array_equal(from_mat, expected)
+                assert_array_equal(from_h5, expected)
+
+    def test_rejected_by_both_backends(self):
+        for name, spec in self.rejected:
+            with self.subTest(index=name):
+                for backend, features in (
+                    ('mat', self.from_mat), ('hdf5', self.from_h5),
+                ):
+                    with self.subTest(backend=backend):
+                        with self.assertRaises(ValueError):
+                            features.get('conv5', feature_slice=spec)
+
+    def test_bool_is_not_read_as_an_integer(self):
+        # isinstance(True, int) is True in Python, so a naive integer check
+        # would silently read True as index 1.
+        with self.assertRaises(ValueError):
+            self.from_h5.get('conv5', feature_slice=True)
+
+    def test_rejection_message_points_at_the_alternative(self):
+        with self.assertRaises(ValueError) as ctx:
+            self.from_h5.get('conv5', feature_slice=np.s_[[1, 2]])
+        self.assertIn('NumPy', str(ctx.exception))
+
+    def test_iter_chunks_takes_no_feature_slice(self):
+        for features in (self.from_mat, self.from_h5):
+            with self.assertRaises(TypeError):
+                list(features.iter_chunks('conv5', feature_slice=np.s_[0:4]))
+
+
 class TestCrossLayerLabelConsistency(unittest.TestCase):
     """Every layer must hold the same labels in the same order.
 
@@ -237,10 +307,7 @@ def tearDown(self):
 
     def _write(self, layer, labels, n_features=8):
         data = np.random.rand(len(labels), n_features)
-        save_features(
-            os.path.join(self.tmpdir.name, layer + '.h5'), data, labels,
-            layer=layer,
-        )
+        save_features(os.path.join(self.tmpdir.name, layer + '.h5'), data, labels)
         return data
 
     def test_consistent_layers_are_accepted(self):
@@ -305,6 +372,32 @@ def test_detects_each_layout(self):
         self.assertEqual(detect_format(self.matdir), 'mat')
         self.assertEqual(detect_format(self.h5dir), 'hdf5')
 
+    def test_unrelated_subdirectory_does_not_change_detection(self):
+        # Detection keys off real feature files, not the bare existence of a
+        # subdirectory, so a stray directory cannot flip a chunked tree to .mat.
+        os.makedirs(os.path.join(self.h5dir, 'notes'))
+        self.assertEqual(detect_format(self.h5dir), 'hdf5')
+        self.assertEqual(Features(self.h5dir).layers, LAYERS)
+
+    def test_both_layouts_present_is_ambiguous(self):
+        # Guessing here means silently reading different data than intended.
+        for layer in LAYERS:
+            os.makedirs(os.path.join(self.h5dir, layer))
+            open(os.path.join(self.h5dir, layer, 'x.mat'), 'w').close()
+        with self.assertRaises(RuntimeError) as ctx:
+            detect_format(self.h5dir)
+        self.assertIn('ambiguous', str(ctx.exception))
+        with self.assertRaises(RuntimeError):
+            Features(self.h5dir)
+        # An explicit format resolves it.
+        self.assertEqual(Features(self.h5dir, format='hdf5').layers, LAYERS)
+
+    def test_empty_directory_is_rejected(self):
+        empty = os.path.join(self.tmpdir.name, 'empty')
+        os.makedirs(empty)
+        with self.assertRaises(RuntimeError):
+            detect_format(empty)
+
     def test_explicit_format_overrides_detection(self):
         self.assertIsInstance(
             Features(self.h5dir, format='hdf5')._Features__stores[0],
@@ -327,7 +420,7 @@ def test_mixed_directories(self):
         for layer, shape in zip(LAYERS, SHAPES):
             data = np.random.rand(len(other_labels), *shape[1:])
             save_features(
-                os.path.join(h5_only, layer + '.h5'), data, other_labels, layer=layer
+                os.path.join(h5_only, layer + '.h5'), data, other_labels
             )
 
         features = Features([self.matdir, h5_only])

From ce1e62345ed353433868737b7ea78b51044022f1 Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Tue, 29 Sep 2026 08:16:20 +0000
Subject: [PATCH 07/11] fix(dataform): close input-validation gaps in
 conversion and iteration

A third review found two ways to get a wrong result without an error, plus
validation that was looser than documented.

convert_features_to_hdf5(batch_size=-1) published an empty layer as a
finished one. range(0, n, -1) is empty, so the loop wrote no samples and
close() then committed that file; the converter's skip-if-exists check
would treat it as complete on every later run. batch_size is now checked
before anything is opened.

iter_chunks validated its arguments only on the path that goes through
FeatureStore. Labels spread across several dpath entries take a fallback
in Features instead, where size=-1 silently yielded nothing and an
out-of-range axis raised IndexError rather than ValueError. The checks
move to Features.iter_chunks, ahead of the dispatch, so both paths reject
the same inputs the same way; a negative axis is normalised there too.

Validation is tightened where it was weaker than claimed:

- bdpy_format_version is required to be an integer scalar rather than
  coerced with int(), which let 1.5 round down into a supported version.
  bool is rejected explicitly, being a subclass of int.
- feature_slice now checks slice components, not just the sign of the
  step: start, stop and step must each be None or an integer, so
  slice(1.5, 3) fails with ValueError instead of reaching NumPy or h5py.

The FeatureSlice alias contained a "ellipsis" forward reference, which is
not a Python name, so typing.get_type_hints() raised NameError on any
annotated member. It now uses builtins.ellipsis under TYPE_CHECKING and
type(Ellipsis) at runtime (types.EllipsisType is 3.10+, the floor here is
3.8), and includes np.integer to match what is accepted.

Two documentation mismatches: the feature_slice docstrings still
advertised array-like indexing, which is rejected, and the docs claimed
np.newaxis raises when a bare np.newaxis is None and so reads as "no
slice"; only a newaxis inside a tuple is refused.

Finally, a test-cleanup miss from the previous commit: an assertion in
test_feature_store.py still pinned "trailing axes are kept whole", a
heuristic detail that was supposed to stop being fixed, while its name
promised a check on the feature axis that it never made. The chunk-policy
assertion moves to test_feature_chunking.py against a representative
large shape, where both sliceable axes really are split -- that one is
issue #144's core promise, not tuning. The n_samples_known=False test now
compares two unknown-size calls that differ only in the current sample
count, which is the contract itself rather than a particular extent.

Every added test was run against the code before this commit; 18 of them
fail there.
---
 bdpy/dataform/_feature_store.py         | 48 ++++++++++++++-----
 bdpy/dataform/feature_hdf5.py           |  9 +++-
 bdpy/dataform/features.py               | 27 +++++++++--
 docs/dataform_features.md               | 12 +++--
 tests/dataform/test_feature_chunking.py | 26 ++++++++---
 tests/dataform/test_feature_hdf5.py     | 23 ++++++++--
 tests/dataform/test_feature_store.py    | 61 +++++++++++++++++++++++--
 7 files changed, 173 insertions(+), 33 deletions(-)

diff --git a/bdpy/dataform/_feature_store.py b/bdpy/dataform/_feature_store.py
index 0bdd3a6..664137a 100644
--- a/bdpy/dataform/_feature_store.py
+++ b/bdpy/dataform/_feature_store.py
@@ -20,7 +20,7 @@
 from abc import ABC, abstractmethod
 from functools import partial
 from multiprocessing import Pool
-from typing import Dict, Iterator, List, Optional, Sequence, Tuple, Union
+from typing import TYPE_CHECKING, Dict, Iterator, List, Optional, Sequence, Tuple, Union
 
 import h5py
 import numpy as np
@@ -53,7 +53,15 @@
 # A slice spec for the feature axes. Deliberately narrow: basic indexing only,
 # i.e. what ``numpy.s_[128:256]`` or ``numpy.s_[8:16, 1:4]`` produces. See
 # _validate_feature_slice for what is rejected and why.
-BasicIndexer = Union[slice, int, "ellipsis"]  # noqa: F821
+# `types.EllipsisType` is 3.10+, and a bare "ellipsis" forward reference does not
+# resolve at runtime (it is not a Python name), which makes
+# typing.get_type_hints() raise. typeshed does define builtins.ellipsis, so the
+# checker gets the real type and the interpreter gets an equivalent object.
+if TYPE_CHECKING:
+    from builtins import ellipsis
+else:
+    ellipsis = type(Ellipsis)
+BasicIndexer = Union[slice, int, np.integer, ellipsis]
 FeatureSlice = Union[BasicIndexer, Tuple[BasicIndexer, ...], None]
 
 
@@ -135,11 +143,22 @@ def _validate_feature_slice(indexers: Tuple) -> None:
         if ix is Ellipsis or isinstance(ix, (int, np.integer)):
             continue
         if isinstance(ix, slice):
-            step = ix.step
-            if step is not None and step < 1:
+            for part, value in (
+                ("start", ix.start), ("stop", ix.stop), ("step", ix.step),
+            ):
+                if value is None:
+                    continue
+                if isinstance(value, (bool, np.bool_)) or not isinstance(
+                    value, (int, np.integer)
+                ):
+                    raise ValueError(
+                        "feature_slice slice {} must be an integer or None, "
+                        "got {!r}.".format(part, value) + advice
+                    )
+            if ix.step is not None and ix.step < 1:
                 raise ValueError(
                     "feature_slice does not support a step below 1 "
-                    "(got {}).".format(step) + advice
+                    "(got {}).".format(ix.step) + advice
                 )
             continue
         if ix is None:
@@ -197,9 +216,12 @@ def read(
             Stimulus labels to read. ``None`` reads every label in store order.
             Otherwise rows come back **in the order given**, and repeated labels
             yield repeated rows.
-        feature_slice : slice, int, array-like or tuple, optional
+        feature_slice : slice, int, Ellipsis or tuple of those, optional
             Index applied to the feature axes (axes 1 and up), as produced by
-            ``numpy.s_[...]``. ``None`` reads the whole feature tensor.
+            ``numpy.s_[...]``. Basic forward indexing only: slices with a
+            positive step, integers, and at most one ``Ellipsis``. Fancy
+            indexing, a negative step, booleans and ``numpy.newaxis`` raise
+            ``ValueError``. ``None`` reads the whole feature tensor.
 
         Returns
         -------
@@ -596,14 +618,18 @@ def _validate_format(f: h5py.File, path: str) -> List[str]:
             "storage.".format(path, FORMAT_VERSION_ATTR)
         )
     raw_version = f.attrs[FORMAT_VERSION_ATTR]
-    try:
-        version = int(raw_version)
-    except (TypeError, ValueError):
+    # Require an actual integer scalar rather than coercing: int(1.5) would
+    # silently round a malformed version down to a supported one. bool is
+    # checked first because it is a subclass of int.
+    if isinstance(raw_version, (bool, np.bool_)) or not isinstance(
+        raw_version, (int, np.integer)
+    ):
         raise RuntimeError(
             "{} has a malformed {} ({!r}); expected an integer.".format(
                 path, FORMAT_VERSION_ATTR, raw_version
             )
-        ) from None
+        )
+    version = int(raw_version)
     if version < 1:
         raise RuntimeError(
             "{} declares feature storage format version {}; versions start at "
diff --git a/bdpy/dataform/feature_hdf5.py b/bdpy/dataform/feature_hdf5.py
index 3a8ec58..3d01560 100644
--- a/bdpy/dataform/feature_hdf5.py
+++ b/bdpy/dataform/feature_hdf5.py
@@ -496,7 +496,7 @@ def convert_features_to_hdf5(
         Overwrite an existing output file instead of skipping it
         (default: False).
     batch_size : int, optional
-        Number of stimulus files read per batch.
+        Number of stimulus files read per batch. Must be positive.
     target_chunk_bytes : int, optional
         Per-chunk byte budget.
     compression : str, optional
@@ -506,9 +506,16 @@ def convert_features_to_hdf5(
 
     Raises
     ------
+    ValueError
+        If `batch_size` is not positive.
     KeyError
         If a requested layer is not present in `src_dir`.
     """
+    if batch_size < 1:
+        # range(0, n, -1) is empty, so a non-positive batch size would write no
+        # samples at all and then publish that empty file as a finished layer.
+        raise ValueError("batch_size must be positive, got {}".format(batch_size))
+
     store = MatFeatureStore(src_dir, ext=ext, key=key)
     selected = list(store.layers) if layers is None else list(layers)
     missing = [layer for layer in selected if layer not in store.layers]
diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py
index c5ec3f0..3d1fe62 100644
--- a/bdpy/dataform/features.py
+++ b/bdpy/dataform/features.py
@@ -162,9 +162,12 @@ def get(
             DNN layer
         label: str or list
             Sample label(s). Rows come back in the order given.
-        feature_slice: slice, int, array-like or tuple, optional
+        feature_slice: slice, int, Ellipsis or tuple of those, optional
             Index applied to the feature axes (axes 1 and up), as produced by
-            ``numpy.s_[...]``. With chunked HDF5 storage this is a genuine
+            ``numpy.s_[...]``. Basic forward indexing only: slices with a
+            positive step, integers, and at most one ``Ellipsis``. Fancy
+            indexing, a negative step, booleans and ``numpy.newaxis`` raise
+            ``ValueError``. With chunked HDF5 storage this is a genuine
             partial read; with the legacy layout the files are loaded in full
             and then sliced.
 
@@ -255,7 +258,8 @@ def iter_chunks(
         ------
         ValueError
             If a unit index (`feature_index`) is in use, which flattens the
-            feature axes and so has no meaningful per-axis iteration.
+            feature axes and so has no meaningful per-axis iteration; if `axis`
+            is out of range; or if `size` is given and is not positive.
 
         Examples
         --------
@@ -275,6 +279,23 @@ def iter_chunks(
         else:
             labels = list(label)
 
+        # Validate here rather than in the store, so that both dispatch paths
+        # below reject the same inputs the same way. The multi-store fallback
+        # does not go through FeatureStore.iter_chunks and would otherwise
+        # silently yield nothing for size < 1, and raise IndexError rather than
+        # ValueError for an out-of-range axis.
+        ndim = len(self.shape(layer))
+        if axis < 0:
+            axis += ndim
+        if not 0 <= axis < ndim:
+            raise ValueError(
+                'axis {} is out of range for a {}-dimensional selection'.format(
+                    axis, ndim
+                )
+            )
+        if size is not None and size < 1:
+            raise ValueError('size must be positive, got {}'.format(size))
+
         store = self.__store_for(labels)
         if store is not None:
             yield from store.iter_chunks(layer, labels, axis=axis, size=size)
diff --git a/docs/dataform_features.md b/docs/dataform_features.md
index 6e1b780..ffc12ec 100644
--- a/docs/dataform_features.md
+++ b/docs/dataform_features.md
@@ -110,16 +110,20 @@ feat = features.get(
 
 `feature_slice` is **basic forward indexing**: slices with a positive step
 (negative `start`/`stop` are fine), integers, a single `Ellipsis`, and tuples of
-those. Anything else -- a negative step, fancy indexing with a list or array,
-booleans, `np.newaxis` -- raises `ValueError` on *both* backends. The restriction
-is what lets the two layouts mean the same thing by the same index; read without
-`feature_slice` and index the result with NumPy when you need more.
+those. Anything else -- a negative step, a non-integer slice bound, fancy
+indexing with a list or array, booleans, or `np.newaxis` *inside a tuple* --
+raises `ValueError` on *both* backends. (A bare `np.newaxis` is simply `None`,
+which is the "no slice" default, so it reads the whole feature tensor rather
+than raising.) The restriction is what lets the two layouts mean the same thing
+by the same index; read without `feature_slice` and index the result with NumPy
+when you need more.
 
 ``` python
 features.get('conv5', feature_slice=np.s_[128:256])     # ok
 features.get('conv5', feature_slice=np.s_[8:16, 1:4])   # ok
 features.get('conv5', feature_slice=np.s_[::-1])        # ValueError
 features.get('conv5', feature_slice=np.s_[[3, 1, 7]])   # ValueError
+features.get('conv5', feature_slice=slice(1.5, 3))      # ValueError
 
 # Full shape without reading anything
 n_stimuli, *feature_shape = features.shape('conv5')
diff --git a/tests/dataform/test_feature_chunking.py b/tests/dataform/test_feature_chunking.py
index d46fec1..3dfb1b2 100644
--- a/tests/dataform/test_feature_chunking.py
+++ b/tests/dataform/test_feature_chunking.py
@@ -68,13 +68,27 @@ def test_target_bytes_is_respected(self):
         self.assertLessEqual(_nbytes(large, dtype), 4 * 1024 * 1024)
         self.assertLess(_nbytes(small, dtype), _nbytes(large, dtype))
 
-    def test_unknown_sample_count_is_not_capped(self):
-        # Resizable datasets start empty; the sample extent must come from the
-        # budget, not from the current (zero) size.
-        chunk = choose_chunk_shape(
-            (0, 256, 13, 13), np.dtype(np.float32), n_samples_known=False
+    def test_unknown_sample_count_ignores_the_current_size(self):
+        # The contract of n_samples_known=False is exactly that the current
+        # sample-axis size is not consulted, so two calls that differ only in
+        # that size must agree. Stated this way the test does not depend on
+        # which extent the policy happens to pick.
+        dtype = np.dtype(np.float32)
+        empty = choose_chunk_shape((0, 256, 13, 13), dtype, n_samples_known=False)
+        huge = choose_chunk_shape(
+            (10 ** 6, 256, 13, 13), dtype, n_samples_known=False
         )
-        self.assertGreater(chunk[0], 1)
+        self.assertEqual(empty, huge)
+
+    def test_both_sliceable_axes_are_chunked(self):
+        # The promise of issue #144: reading a slice of channels must not have
+        # to read every channel, which requires the feature axis to be split as
+        # well as the sample axis. Asserted on a representative layer, since a
+        # small one legitimately fits in a single chunk.
+        shape = (1200, 256, 13, 13)
+        chunk = choose_chunk_shape(shape, np.dtype(np.float32))
+        self.assertLess(chunk[0], shape[0])
+        self.assertLess(chunk[1], shape[1])
 
     def test_deterministic(self):
         shape = (1200, 256, 13, 13)
diff --git a/tests/dataform/test_feature_hdf5.py b/tests/dataform/test_feature_hdf5.py
index d6827d7..dab9738 100644
--- a/tests/dataform/test_feature_hdf5.py
+++ b/tests/dataform/test_feature_hdf5.py
@@ -236,9 +236,15 @@ def test_missing_version_is_rejected(self):
             HDF5FeatureStore(self.tmpdir.name)
 
     def test_malformed_version_is_rejected(self):
-        self._mutate(lambda f: f.attrs.__setitem__(FORMAT_VERSION_ATTR, 'v1'))
-        with self.assertRaises(RuntimeError):
-            HDF5FeatureStore(self.tmpdir.name)
+        # Anything that is not an integer scalar, including a float that would
+        # round down to a supported version under int().
+        for bad in ('v1', 1.5, 1.0, True, np.True_):
+            with self.subTest(version=bad):
+                self._mutate(
+                    lambda f, b=bad: f.attrs.__setitem__(FORMAT_VERSION_ATTR, b)
+                )
+                with self.assertRaises(RuntimeError):
+                    HDF5FeatureStore(self.tmpdir.name)
 
     def test_version_below_one_is_rejected(self):
         for bad in (0, -1):
@@ -537,6 +543,17 @@ def flaky(store, layer, labels=None, feature_slice=None):
 
         return original, flaky
 
+    def test_non_positive_batch_size_is_refused(self):
+        # range(0, n, -1) is empty, so this used to write nothing and then
+        # publish the empty file as a finished layer.
+        for bad in (-1, 0):
+            with self.subTest(batch_size=bad):
+                with self.assertRaises(ValueError):
+                    convert_features_to_hdf5(
+                        self.matdir, self.h5dir, batch_size=bad
+                    )
+                self.assertFalse(os.path.exists(self.h5dir))
+
     def test_failed_conversion_leaves_no_file(self):
         original, flaky = self._fail_on_second_batch()
         MatFeatureStore.read = flaky
diff --git a/tests/dataform/test_feature_store.py b/tests/dataform/test_feature_store.py
index a70d9c5..c46624f 100644
--- a/tests/dataform/test_feature_store.py
+++ b/tests/dataform/test_feature_store.py
@@ -162,6 +162,39 @@ def test_default_size_is_chunk_aligned(self):
         # Every block but the last is a whole chunk along the axis.
         self.assertLessEqual(max(sizes), max(1, extent))
 
+    def test_multi_directory_validates_like_a_single_one(self):
+        # The multi-store path does not go through FeatureStore.iter_chunks, so
+        # it used to skip these checks entirely: size=-1 yielded nothing and a
+        # bad axis raised IndexError instead of ValueError.
+        other = os.path.join(self.tmpdir.name, 'other')
+        os.makedirs(other)
+        other_labels = ['other%04d' % i for i in range(3)]
+        for layer, shape in zip(LAYERS, SHAPES):
+            save_features(
+                os.path.join(other, layer + '.h5'),
+                np.random.rand(len(other_labels), *shape[1:]),
+                other_labels,
+            )
+        spread = Features([self.h5dir, other])
+        both = [LABELS[0], other_labels[0]]
+        # Spanning two directories is what takes the fallback path.
+        self.assertEqual(len(spread.get('conv5', label=both)), 2)
+
+        for kwargs in ({'size': -1}, {'size': 0}, {'axis': 9}, {'axis': -9}):
+            with self.subTest(**kwargs):
+                with self.assertRaises(ValueError):
+                    list(spread.iter_chunks('conv5', label=both, **kwargs))
+                with self.assertRaises(ValueError):
+                    list(self.from_h5.iter_chunks('conv5', **kwargs))
+
+    def test_negative_axis_is_accepted(self):
+        by_negative = list(self.from_h5.iter_chunks('conv5', axis=-3, size=4))
+        by_positive = list(self.from_h5.iter_chunks('conv5', axis=1, size=4))
+        self.assertEqual(len(by_negative), len(by_positive))
+        for (sl_a, a), (sl_b, b) in zip(by_negative, by_positive):
+            self.assertEqual(sl_a, sl_b)
+            assert_array_equal(a, b)
+
     def test_rejects_bad_axis(self):
         with self.assertRaises(ValueError):
             list(self.from_h5.iter_chunks('conv5', axis=9))
@@ -211,12 +244,14 @@ def spy(dataset, args):
         self.assertEqual(rows, [1, 9])
         assert_array_equal(out, self.stacked['conv5'][[9, 1, 9]])
 
-    def test_dataset_is_chunked_on_both_sliceable_axes(self):
+    def test_dataset_is_chunked(self):
+        # This fixture is small enough to fit in one chunk, so only the fact
+        # that chunking is enabled can be asserted here. That both sliceable
+        # axes are actually split is a property of the chunk-shape policy and
+        # is tested against a representative large shape in
+        # tests/dataform/test_feature_chunking.py.
         with h5py.File(os.path.join(self.h5dir, 'conv5.h5'), 'r') as f:
-            chunks = f['features'].chunks
-            shape = f['features'].shape
-        self.assertIsNotNone(chunks)
-        self.assertEqual(chunks[2:], shape[2:])  # spatial axes kept whole
+            self.assertIsNotNone(f['features'].chunks)
 
 
 class TestFeatureSliceValidation(_BackendPair):
@@ -250,6 +285,10 @@ class TestFeatureSliceValidation(_BackendPair):
         ('numpy bool', np.True_),
         ('newaxis in tuple', (np.newaxis, slice(None))),
         ('two ellipses', (Ellipsis, Ellipsis)),
+        ('float start', slice(1.5, 3)),
+        ('float stop', slice(1, 3.5)),
+        ('float step', slice(None, None, 2.0)),
+        ('bool start', slice(True, 3)),
     ]
 
     def test_accepted_agree_across_backends_and_numpy(self):
@@ -289,6 +328,18 @@ def test_iter_chunks_takes_no_feature_slice(self):
                 list(features.iter_chunks('conv5', feature_slice=np.s_[0:4]))
 
 
+class TestFeatureSliceAnnotation(unittest.TestCase):
+    def test_type_hints_resolve(self):
+        # FeatureSlice once contained an "ellipsis" forward reference, which is
+        # not a Python name, so any tool resolving annotations blew up.
+        import typing
+
+        from bdpy.dataform._feature_store import FeatureStore
+
+        typing.get_type_hints(FeatureStore.read)
+        typing.get_type_hints(MatFeatureStore.read)
+
+
 class TestCrossLayerLabelConsistency(unittest.TestCase):
     """Every layer must hold the same labels in the same order.
 

From b4acc01e8c33c7894436313148bf5d38a78cc96a Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Tue, 29 Sep 2026 09:10:32 +0000
Subject: [PATCH 08/11] docs(dataform): correct newaxis wording and drop a
 duplicate import

The feature_slice docstrings on FeatureStore.read and Features.get still
said numpy.newaxis raises ValueError. It does not: numpy.newaxis *is*
None, which is the "no slice" default, so a bare newaxis reads the whole
feature tensor and only a newaxis inside a tuple is refused. The Markdown
docs already said this; now all three agree, and both docstrings also
mention the non-integer slice bound check.

tests/dataform/test_feature_hdf5.py imported SUPPORTED_FORMAT_VERSION
twice. Ruff never caught it because pyproject excludes tests from linting;
an AST sweep over tests/dataform confirms this was the only duplicate.

No behaviour change: 349 passed either way.
---
 bdpy/dataform/_feature_store.py     |  8 +++++---
 bdpy/dataform/features.py           | 12 +++++++-----
 tests/dataform/test_feature_hdf5.py |  3 +--
 3 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/bdpy/dataform/_feature_store.py b/bdpy/dataform/_feature_store.py
index 664137a..9ebddc5 100644
--- a/bdpy/dataform/_feature_store.py
+++ b/bdpy/dataform/_feature_store.py
@@ -219,9 +219,11 @@ def read(
         feature_slice : slice, int, Ellipsis or tuple of those, optional
             Index applied to the feature axes (axes 1 and up), as produced by
             ``numpy.s_[...]``. Basic forward indexing only: slices with a
-            positive step, integers, and at most one ``Ellipsis``. Fancy
-            indexing, a negative step, booleans and ``numpy.newaxis`` raise
-            ``ValueError``. ``None`` reads the whole feature tensor.
+            positive step and integer bounds, integers, and at most one
+            ``Ellipsis``. Fancy indexing, a negative step, a non-integer slice
+            bound, booleans, and ``numpy.newaxis`` *inside a tuple* raise
+            ``ValueError``. A bare ``numpy.newaxis`` is ``None``, the "no
+            slice" default, so it reads the whole feature tensor.
 
         Returns
         -------
diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py
index 3d1fe62..5adf21e 100644
--- a/bdpy/dataform/features.py
+++ b/bdpy/dataform/features.py
@@ -165,11 +165,13 @@ def get(
         feature_slice: slice, int, Ellipsis or tuple of those, optional
             Index applied to the feature axes (axes 1 and up), as produced by
             ``numpy.s_[...]``. Basic forward indexing only: slices with a
-            positive step, integers, and at most one ``Ellipsis``. Fancy
-            indexing, a negative step, booleans and ``numpy.newaxis`` raise
-            ``ValueError``. With chunked HDF5 storage this is a genuine
-            partial read; with the legacy layout the files are loaded in full
-            and then sliced.
+            positive step and integer bounds, integers, and at most one
+            ``Ellipsis``. Fancy indexing, a negative step, a non-integer slice
+            bound, booleans, and ``numpy.newaxis`` *inside a tuple* raise
+            ``ValueError``. A bare ``numpy.newaxis`` is ``None``, the "no
+            slice" default, so it reads the whole feature tensor. With chunked
+            HDF5 storage this is a genuine partial read; with the legacy layout
+            the files are loaded in full and then sliced.
 
         Returns
         -------
diff --git a/tests/dataform/test_feature_hdf5.py b/tests/dataform/test_feature_hdf5.py
index dab9738..8ba70da 100644
--- a/tests/dataform/test_feature_hdf5.py
+++ b/tests/dataform/test_feature_hdf5.py
@@ -12,13 +12,12 @@
 
 from bdpy.dataform import Features
 from bdpy.dataform._feature_store import (
-    MatFeatureStore,
-    SUPPORTED_FORMAT_VERSION,
     FORMAT_ATTR,
     FORMAT_NAME,
     FORMAT_VERSION_ATTR,
     SUPPORTED_FORMAT_VERSION,
     HDF5FeatureStore,
+    MatFeatureStore,
 )
 from bdpy.dataform.feature_hdf5 import (
     FeatureWriter,

From 2763672978facfb772c946493e467f7ac36c1a3b Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Tue, 29 Sep 2026 13:35:37 +0000
Subject: [PATCH 09/11] perf(dataform): read a .mat feature axis once per
 iteration, not per slab

FeatureStore.iter_chunks called read() for every slab, and MatFeatureStore.read
loads every selected stimulus file. A k-slab feature-axis iteration therefore
re-read the whole layer k times: 1088 s for conv3_4 (100, 256, 56, 56) against
5 s for one read plus in-memory slicing. The legacy layout has no file boundary
on the feature axes, so the per-slab read bought nothing.

Split the feature-axis branch of the loop into _iter_feature_slabs, which gets
the slab slices the base class has already computed, and override it in
MatFeatureStore to read the selection once and cut the slabs out of it. Axis
and size validation, the slab boundaries and the sample-axis branch all stay in
one place, so the two backends cannot drift apart on semantics. The slices are
handed over as a generator rather than a list, so iterating a large axis at
size 1 does not first materialize one slice object per slab.

Sample-axis iteration is deliberately untouched: there each slab reads only its
own files, which is the one case where the legacy layout can stream a layer
larger than memory. Nothing is kept on the store, so this is a per-call
strategy and not a cache. The slabs are views into the selection, exactly as
they were when read separately, so holding a slab still holds the array.

Iterating an empty axis now returns before reading anything, matching the
previous behaviour where the loop body simply never ran.
---
 bdpy/dataform/_feature_store.py      |  96 +++++++++++++++++++---
 bdpy/dataform/features.py            |  15 +++-
 docs/dataform_features.md            |  14 +++-
 tests/dataform/test_feature_store.py | 115 ++++++++++++++++++++++++++-
 4 files changed, 218 insertions(+), 22 deletions(-)

diff --git a/bdpy/dataform/_feature_store.py b/bdpy/dataform/_feature_store.py
index 9ebddc5..58f283c 100644
--- a/bdpy/dataform/_feature_store.py
+++ b/bdpy/dataform/_feature_store.py
@@ -278,6 +278,15 @@ def iter_chunks(
         ------
         tuple of (slice, numpy.ndarray)
             The slice applied to `axis`, and the corresponding slab.
+
+        Notes
+        -----
+        How much is resident at once depends on the backend. Along the sample
+        axis every backend reads one slab at a time. Along a feature axis only
+        a backend with real partial reads can do so: the legacy ``.mat``
+        layout has no file boundary there, so it reads the selection once and
+        yields views into it, and the blocks keep that array alive for as long
+        as they are held.
         """
         selected_shape = self._selected_shape(layer, labels)
         if axis < 0:
@@ -295,18 +304,47 @@ def iter_chunks(
         if size < 1:
             raise ValueError("size must be positive, got {}".format(size))
 
-        for start in range(0, length, size):
-            sl = slice(start, min(start + size, length))
-            if axis == 0:
-                block_labels = (
-                    list(self.labels) if labels is None else list(labels)
-                )[sl]
-                yield sl, self.read(layer, block_labels)
-            else:
-                # axis >= 1 indexes the feature axes, where axis 1 of the array
-                # is entry 0 of the feature-slice tuple.
-                indexers = (slice(None),) * (axis - 1) + (sl,)
-                yield sl, self.read(layer, labels, indexers)
+        if length == 0:
+            # Nothing to iterate. Returning here rather than falling into the
+            # loop keeps a backend that reads the selection up front from
+            # reading anything the per-slab path would not have read either.
+            return
+
+        # A generator, not a list: a materialized one would cost an object per
+        # slab, which for a large axis at size 1 is millions of them.
+        slices = (
+            slice(start, min(start + size, length))
+            for start in range(0, length, size)
+        )
+        if axis == 0:
+            all_labels = list(self.labels) if labels is None else list(labels)
+            for sl in slices:
+                yield sl, self.read(layer, all_labels[sl])
+        else:
+            yield from self._iter_feature_slabs(layer, labels, axis, slices)
+
+    def _iter_feature_slabs(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]],
+        axis: int,
+        slices: Iterator[slice],
+    ) -> Iterator[Tuple[slice, np.ndarray]]:
+        """Yield the feature-axis slabs of `layer` named by `slices`.
+
+        This is the seam for a backend that cannot read a feature axis
+        partially. `axis` is already normalized and validated to a feature axis
+        (``>= 1``), and `slices` partitions that axis, so an implementation
+        decides only *how* a slab is obtained -- never which slabs there are.
+
+        The default reads each slab on its own, which is what a backend with
+        real partial reads wants.
+        """
+        for sl in slices:
+            # axis >= 1 indexes the feature axes, where axis 1 of the array is
+            # entry 0 of the feature-slice tuple.
+            indexers = (slice(None),) * (axis - 1) + (sl,)
+            yield sl, self.read(layer, labels, indexers)
 
     def _selected_shape(
         self,
@@ -341,7 +379,10 @@ class MatFeatureStore(FeatureStore):
     under `key` with a leading sample axis. There is no file boundary on the
     feature axes, so `feature_slice` is applied after the full stimulus files
     have been loaded and concatenated -- the read cost is the same as before,
-    and the slice only saves the caller from doing it themselves.
+    and the slice only saves the caller from doing it themselves. For the same
+    reason :meth:`iter_chunks` along a feature axis reads the selection once
+    and cuts the slabs out of it in memory; only sample-axis iteration is
+    bounded by the slab.
 
     Parameters
     ----------
@@ -415,6 +456,35 @@ def read(
             features = features[(slice(None), *indexers)]
         return features
 
+    def _iter_feature_slabs(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]],
+        axis: int,
+        slices: Iterator[slice],
+    ) -> Iterator[Tuple[slice, np.ndarray]]:
+        """Read the selection once, then cut every slab out of it in memory.
+
+        This layout has no file boundary on the feature axes, so one slab
+        already costs a read of every selected stimulus file. Reading per slab
+        therefore re-read the whole selection once per slab, making a k-slab
+        iteration k times the cost of a single read.
+
+        Nothing is kept on the store, so this is a per-call strategy rather
+        than a cache. The selection is a plain array held for the duration of
+        the iteration, and the yielded slabs are views into it -- exactly as
+        they were when each slab was read separately -- so holding a slab keeps
+        the whole selection alive. Sample-axis iteration is deliberately left
+        to the base class, where each slab reads only its own files and a layer
+        larger than memory can still be processed end to end.
+        """
+        features = self.read(layer, labels)
+        for sl in slices:
+            # `features` still carries its sample axis, so axis N of the array
+            # is entry N here -- unlike the feature-slice tuple the base class
+            # builds, whose entry 0 is axis 1.
+            yield sl, features[(slice(None),) * axis + (sl,)]
+
     def _collect_layers(self) -> List[str]:
         return sorted(
             d
diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py
index 5adf21e..d65bff0 100644
--- a/bdpy/dataform/features.py
+++ b/bdpy/dataform/features.py
@@ -246,10 +246,17 @@ def iter_chunks(
 
         Notes
         -----
-        The one-slab-at-a-time guarantee holds when the requested labels all
-        live in a single feature directory, which is the usual case. When they
-        are spread across several `dpath` entries, no single store can stream
-        them, so the selection is read in full and then split.
+        The one-slab-at-a-time guarantee is a property of chunked HDF5 storage
+        and of `axis=0` on any layout. The legacy `.mat` layout has no file
+        boundary on the feature axes, so iterating one reads the selection once
+        -- the same peak as a single :meth:`get` -- and yields views into it,
+        which keep that array alive for as long as a block is held. It is still
+        read once rather than once per slab.
+
+        The guarantee also assumes the requested labels all live in a single
+        feature directory, which is the usual case. When they are spread across
+        several `dpath` entries, no single store can stream them, so the
+        selection is read in full and then split.
 
         Yields
         ------
diff --git a/docs/dataform_features.md b/docs/dataform_features.md
index ffc12ec..cccb82a 100644
--- a/docs/dataform_features.md
+++ b/docs/dataform_features.md
@@ -140,10 +140,16 @@ for sl, block in features.iter_chunks('conv5', axis=1):
 ```
 
 `axis=0` iterates over stimuli instead of features. `iter_chunks` takes no
-`feature_slice` -- slice each `block` as it comes out instead. The
-one-slab-at-a-time guarantee holds when the requested labels all live in one
-feature directory; when they are spread across several `dpath` entries no single
-store can stream them, so the selection is read in full and then split.
+`feature_slice` -- slice each `block` as it comes out instead.
+
+Only one slab is resident at a time for chunked HDF5, and for `axis=0` on either
+layout. On the legacy `.mat` layout there is no file boundary on the feature
+axes, so iterating one reads the selection once -- the same peak as a single
+`get()` -- and yields views into it, which keep that array alive for as long as a
+block is held. It is still read once rather than once per slab, which is what
+makes `iter_chunks` usable on that layout at all. The same full-read fallback
+applies when the requested labels are spread across several `dpath` entries,
+since no single store can stream them.
 
 ### Writing
 
diff --git a/tests/dataform/test_feature_store.py b/tests/dataform/test_feature_store.py
index c46624f..9b34585 100644
--- a/tests/dataform/test_feature_store.py
+++ b/tests/dataform/test_feature_store.py
@@ -117,9 +117,122 @@ def test_statistic_matches(self):
                     )
 
 
+class TestMatIterChunksReadsOnce(_BackendPair):
+    """The legacy backend must read a feature axis once, not once per slab.
+
+    ``MatFeatureStore.read`` loads every selected stimulus file, so reading a
+    slab at a time made a k-slab iteration cost k full reads of the layer
+    (measured: 1088 s against 5 s for one read plus in-memory slicing).
+    """
+
+    def setUp(self):
+        super().setUp()
+        self.store = MatFeatureStore(self.matdir)
+        self.reads = []
+        original = self.store.read
+
+        def spy(*args, **kwargs):
+            self.reads.append((args, kwargs))
+            return original(*args, **kwargs)
+
+        self.store.read = spy
+
+    def test_feature_axis_reads_the_selection_once(self):
+        slab_counts = {}
+        for size in (1, 3, 4, 24, 100, None):
+            with self.subTest(size=size):
+                self.reads.clear()
+                blocks = list(self.store.iter_chunks('conv5', axis=1, size=size))
+                self.assertEqual(len(self.reads), 1)
+                slab_counts[size] = len(blocks)
+        # Not vacuous: the slab count really does vary, and the worst case is
+        # the 24 reads this test exists to prevent.
+        self.assertEqual(slab_counts[1], 24)
+        self.assertEqual(slab_counts[100], 1)
+        self.assertGreater(len(set(slab_counts.values())), 1)
+
+    def test_sample_axis_still_reads_one_slab_at_a_time(self):
+        # Preloading here would destroy the one guarantee the legacy backend
+        # can make: a layer larger than memory streams along the sample axis.
+        blocks = list(self.store.iter_chunks('conv5', axis=0, size=4))
+        self.assertEqual(len(blocks), 3)
+        self.assertEqual(len(self.reads), 3)
+        for (args, _), (sl, _) in zip(self.reads, blocks):
+            self.assertEqual(list(args[1]), LABELS[sl])
+
+    def test_nothing_is_cached_between_calls(self):
+        for call in (1, 2):
+            with self.subTest(call=call):
+                self.reads.clear()
+                list(self.store.iter_chunks('conv5', axis=1, size=4))
+                self.assertEqual(len(self.reads), 1)
+
+    def test_blocks_match_ground_truth_and_the_other_backend(self):
+        cases = [
+            ('conv5', None, 1),
+            ('conv5', None, 2),
+            ('conv5', None, 3),
+            ('conv5', None, -1),
+            ('conv5', None, -3),
+            ('fc8', None, 1),
+            ('conv5', ['img0009', 'img0001'], 1),
+            ('conv5', ['img0009', 'img0001', 'img0009'], 2),
+            ('conv5', ['img0003'], 3),
+        ]
+        for layer, label, axis in cases:
+            for size in (1, 3, 1000, None):
+                with self.subTest(layer=layer, label=label, axis=axis, size=size):
+                    rows = (
+                        slice(None) if label is None
+                        else [LABELS.index(s) for s in label]
+                    )
+                    expected = self.stacked[layer][rows]
+                    norm = axis if axis >= 0 else axis + expected.ndim
+                    from_mat = list(self.from_mat.iter_chunks(
+                        layer, label=label, axis=axis, size=size))
+                    for sl, block in from_mat:
+                        assert_array_equal(
+                            block, expected[(slice(None),) * norm + (sl,)]
+                        )
+                    # The two backends need not agree on slab boundaries -- the
+                    # HDF5 default follows the on-disk chunk extent -- so the
+                    # reassembled layer is what has to match.
+                    from_h5 = list(self.from_h5.iter_chunks(
+                        layer, label=label, axis=axis, size=size))
+                    assert_array_equal(
+                        np.concatenate([b for _, b in from_mat], axis=norm),
+                        np.concatenate([b for _, b in from_h5], axis=norm),
+                    )
+
+    def test_writing_to_one_block_does_not_disturb_another(self):
+        # Slabs are views into one array now, but they partition the axis, so
+        # they cannot alias each other.
+        blocks = list(self.from_mat.iter_chunks('conv5', axis=1, size=4))
+        blocks[0][1][:] = 0
+        assert_array_equal(blocks[1][1], self.stacked['conv5'][:, 4:8])
+
+    def test_validation_precedes_any_read(self):
+        cases = (
+            ({'axis': 9}, 'axis 9 is out of range for selected shape (11, 24, 5, 5)'),
+            ({'axis': -9}, 'axis -5 is out of range for selected shape (11, 24, 5, 5)'),
+            ({'size': 0}, 'size must be positive, got 0'),
+        )
+        for kwargs, message in cases:
+            with self.subTest(**kwargs):
+                self.reads.clear()
+                with self.assertRaises(ValueError) as ctx:
+                    list(self.store.iter_chunks('conv5', **kwargs))
+                self.assertEqual(str(ctx.exception), message)
+                # Reading the selection up front must not outrun validation.
+                self.assertEqual(len(self.reads), 0)
+                with self.assertRaises(ValueError) as from_h5:
+                    list(HDF5FeatureStore(self.h5dir).iter_chunks('conv5', **kwargs))
+                self.assertEqual(str(from_h5.exception), message)
+
+
 class TestIterChunks(_BackendPair):
     def test_reassembles_to_a_full_read(self):
-        for axis in (0, 1):
+        for axis in (0, 1, 2, 3):
             for features in (self.from_mat, self.from_h5):
                 with self.subTest(axis=axis, backend=type(features).__name__):
                     blocks = list(features.iter_chunks('conv5', axis=axis, size=4))

From 47c79098dc513524c71b8fb51367ad5a4d49e12e Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Tue, 6 Oct 2026 05:10:58 +0000
Subject: [PATCH 10/11] refactor(dataform): separate iter_chunks index
 arithmetic from I/O

The arithmetic behind iter_chunks lived inline next to the reads, and was
written more than once: axis normalisation in both Features.iter_chunks
and FeatureStore.iter_chunks (with different messages), slab
partitioning in the store and again in the multi-directory fallback, and
the "apply sl on axis" index three ways. feature_slice validation was
left to each backend's read(), and the .mat backend ran it only after
loading every file. Testing any of this meant going through files, which
is why the iteration test was a 36-case product over both backends.

Move the arithmetic into small pure functions -- _normalize_axis,
_check_size, _partition, _axis_index, _default_iter_size and
_iter_array_slabs -- and use them on every path; the
FeatureStore._default_iter_size method they replace is removed. Make
FeatureStore.read a template method that validates feature_slice and
then calls the backend's _read, so feature_slice validation happens
once, before the backend reads anything, on every backend.

Two behaviour changes:
- With size=None, the multi-directory fallback now splits by the same
  byte budget as a store without on-disk chunking, instead of yielding
  the whole axis as one slab. The data is unchanged; statistic() reads
  through get() and is unaffected.
- An out-of-range axis is reported as the caller passed it (axis -9 for
  a 4-D selection) rather than after normalisation (it used to say -5).

Tests: the helpers are checked in memory (TestIterationHelpers),
including the table of rejected feature_slice specs against the
validator. I/O tests keep only what each backend owns: cutting the right
slab (one feature-axis case per axis with a partial last slab and
repeated labels, replacing the product), rejecting a bad feature_slice
before its _read runs, and not overriding read(). The fallback's default
size is checked with the budget patched small enough to need several
slabs.

TestIterChunks tests that the new ones subsume are dropped:
test_yielded_slices_address_the_right_block and test_restricts_to_labels
(the per-axis ground-truth test),
test_slices_cover_the_axis_without_overlap (_partition) and
test_rejects_bad_axis (test_validation_precedes_selection_read);
test_reassembles_to_a_full_read keeps only the sample axis, which no
other test iterates on both backends.
---
 bdpy/dataform/_feature_store.py      | 154 ++++++++++-----
 bdpy/dataform/features.py            |  47 +++--
 tests/dataform/test_feature_store.py | 278 ++++++++++++++++++---------
 3 files changed, 310 insertions(+), 169 deletions(-)

diff --git a/bdpy/dataform/_feature_store.py b/bdpy/dataform/_feature_store.py
index 58f283c..d57c7ff 100644
--- a/bdpy/dataform/_feature_store.py
+++ b/bdpy/dataform/_feature_store.py
@@ -173,6 +173,77 @@ def _validate_feature_slice(indexers: Tuple) -> None:
         )
 
 
+def _normalize_axis(axis: int, shape: Tuple[int, ...]) -> int:
+    """Return `axis` as a non-negative axis of an array of `shape`.
+
+    Raises
+    ------
+    ValueError
+        If `axis` is out of range for `shape`.
+    """
+    normalized = axis + len(shape) if axis < 0 else axis
+    if not 0 <= normalized < len(shape):
+        raise ValueError(
+            "axis {} is out of range for selected shape {}".format(axis, shape)
+        )
+    return normalized
+
+
+def _check_size(size: int) -> None:
+    """Reject a slab size below 1."""
+    if size < 1:
+        raise ValueError("size must be positive, got {}".format(size))
+
+
+def _partition(length: int, size: int) -> Iterator[slice]:
+    """Split ``range(length)`` into consecutive slices of at most `size`.
+
+    The slices cover the range exactly once, in order; only the last may be
+    shorter. ``length == 0`` yields nothing. `size` must already be positive.
+
+    A generator rather than a list: a materialized one would cost an object per
+    slab, which for a large axis at size 1 is millions of them.
+    """
+    for start in range(0, length, size):
+        yield slice(start, min(start + size, length))
+
+
+def _axis_index(axis: int, sl: slice) -> Tuple[slice, ...]:
+    """Index tuple that applies `sl` to `axis` and leaves earlier axes whole."""
+    return (slice(None),) * axis + (sl,)
+
+
+def _default_iter_size(
+    shape: Tuple[int, ...], axis: int, itemsize: int, extent: Optional[int]
+) -> int:
+    """Default slab size along `axis` of an array of `shape`.
+
+    The on-disk chunk `extent` along the axis when there is one, so iteration is
+    chunk-aligned; otherwise as many elements as fit in the byte budget. Always
+    at least 1 and at most the axis length (when that is positive).
+    """
+    if extent is not None:
+        return max(1, min(extent, shape[axis]))
+    other = 1
+    for i, n in enumerate(shape):
+        if i != axis:
+            other *= int(n)
+    per_element = max(1, other * itemsize)
+    # Read the module constant at call time so it can be patched in tests.
+    return max(1, min(shape[axis], DEFAULT_TARGET_CHUNK_BYTES // per_element))
+
+
+def _iter_array_slabs(
+    array: np.ndarray, axis: int, slices: Iterator[slice]
+) -> Iterator[Tuple[slice, np.ndarray]]:
+    """Yield ``(sl, view)`` for each of `slices` applied to `axis` of `array`.
+
+    The slabs are views, so holding one keeps `array` alive.
+    """
+    for sl in slices:
+        yield sl, array[_axis_index(axis, sl)]
+
+
 class FeatureStore(ABC):
     """Interface to a collection of DNN features on disk.
 
@@ -199,7 +270,6 @@ def shape(self, layer: str) -> Tuple[int, ...]:
     def dtype(self, layer: str) -> np.dtype:
         """Dtype of the features in `layer`."""
 
-    @abstractmethod
     def read(
         self,
         layer: str,
@@ -208,6 +278,10 @@ def read(
     ) -> np.ndarray:
         """Read features from `layer`.
 
+        `feature_slice` is validated here, before any I/O, so every backend
+        accepts and rejects exactly the same specs. Backends implement
+        :meth:`_read` and must not override this method.
+
         Parameters
         ----------
         layer : str
@@ -230,6 +304,16 @@ def read(
         numpy.ndarray
             Array of shape ``(n_labels, *sliced_feature_shape)``.
         """
+        return self._read(layer, labels, _normalize_feature_slice(feature_slice))
+
+    @abstractmethod
+    def _read(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]],
+        indexers: Tuple,
+    ) -> np.ndarray:
+        """Backend read: `indexers` is already validated, ``()`` for no slice."""
 
     def chunk_extent(self, layer: str, axis: int) -> Optional[int]:
         """On-disk chunk extent of `layer` along `axis`, if the store has one.
@@ -289,20 +373,17 @@ def iter_chunks(
         as they are held.
         """
         selected_shape = self._selected_shape(layer, labels)
-        if axis < 0:
-            axis += len(selected_shape)
-        if not 0 <= axis < len(selected_shape):
-            raise ValueError(
-                "axis {} is out of range for selected shape {}".format(
-                    axis, selected_shape
-                )
-            )
+        axis = _normalize_axis(axis, selected_shape)
 
         length = selected_shape[axis]
         if size is None:
-            size = self._default_iter_size(layer, axis, selected_shape)
-        if size < 1:
-            raise ValueError("size must be positive, got {}".format(size))
+            size = _default_iter_size(
+                selected_shape,
+                axis,
+                np.dtype(self.dtype(layer)).itemsize,
+                self.chunk_extent(layer, axis),
+            )
+        _check_size(size)
 
         if length == 0:
             # Nothing to iterate. Returning here rather than falling into the
@@ -310,12 +391,7 @@ def iter_chunks(
             # reading anything the per-slab path would not have read either.
             return
 
-        # A generator, not a list: a materialized one would cost an object per
-        # slab, which for a large axis at size 1 is millions of them.
-        slices = (
-            slice(start, min(start + size, length))
-            for start in range(0, length, size)
-        )
+        slices = _partition(length, size)
         if axis == 0:
             all_labels = list(self.labels) if labels is None else list(labels)
             for sl in slices:
@@ -342,9 +418,9 @@ def _iter_feature_slabs(
         """
         for sl in slices:
             # axis >= 1 indexes the feature axes, where axis 1 of the array is
-            # entry 0 of the feature-slice tuple.
-            indexers = (slice(None),) * (axis - 1) + (sl,)
-            yield sl, self.read(layer, labels, indexers)
+            # entry 0 of the feature-slice tuple. Go through read(), the single
+            # entry point that validates.
+            yield sl, self.read(layer, labels, _axis_index(axis - 1, sl))
 
     def _selected_shape(
         self,
@@ -356,20 +432,6 @@ def _selected_shape(
         n_samples = full[0] if labels is None else len(labels)
         return (n_samples, *full[1:])
 
-    def _default_iter_size(
-        self, layer: str, axis: int, selected_shape: Tuple[int, ...]
-    ) -> int:
-        extent = self.chunk_extent(layer, axis)
-        if extent is not None:
-            return max(1, min(extent, selected_shape[axis]))
-        # No on-disk chunking to align with: fall back to the byte budget.
-        itemsize = np.dtype(self.dtype(layer)).itemsize
-        other = 1
-        for i, n in enumerate(selected_shape):
-            if i != axis:
-                other *= int(n)
-        per_element = max(1, other * itemsize)
-        return max(1, min(selected_shape[axis], DEFAULT_TARGET_CHUNK_BYTES // per_element))
 
 
 class MatFeatureStore(FeatureStore):
@@ -432,11 +494,11 @@ def dtype(self, layer: str) -> np.dtype:
         sample = _load_array_with_key(self._key, self.path(layer, self._labels[0]))
         return sample.dtype
 
-    def read(
+    def _read(
         self,
         layer: str,
-        labels: Optional[Sequence[str]] = None,
-        feature_slice: FeatureSlice = None,
+        labels: Optional[Sequence[str]],
+        indexers: Tuple,
     ) -> np.ndarray:
         if labels is None:
             labels = self._labels
@@ -451,7 +513,6 @@ def read(
                 arrays = pool.map(load, paths)
         features = np.concatenate(arrays, axis=0)
 
-        indexers = _normalize_feature_slice(feature_slice)
         if indexers:
             features = features[(slice(None), *indexers)]
         return features
@@ -478,12 +539,7 @@ def _iter_feature_slabs(
         to the base class, where each slab reads only its own files and a layer
         larger than memory can still be processed end to end.
         """
-        features = self.read(layer, labels)
-        for sl in slices:
-            # `features` still carries its sample axis, so axis N of the array
-            # is entry N here -- unlike the feature-slice tuple the base class
-            # builds, whose entry 0 is axis 1.
-            yield sl, features[(slice(None),) * axis + (sl,)]
+        yield from _iter_array_slabs(self.read(layer, labels), axis, slices)
 
     def _collect_layers(self) -> List[str]:
         return sorted(
@@ -562,14 +618,12 @@ def chunk_extent(self, layer: str, axis: int) -> Optional[int]:
             return None
         return int(chunks[axis])
 
-    def read(
+    def _read(
         self,
         layer: str,
-        labels: Optional[Sequence[str]] = None,
-        feature_slice: FeatureSlice = None,
+        labels: Optional[Sequence[str]],
+        indexers: Tuple,
     ) -> np.ndarray:
-        indexers = _normalize_feature_slice(feature_slice)
-
         with self._open(layer) as f:
             dset = f[FEATURES_DATASET]
 
diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py
index d65bff0..610dafc 100644
--- a/bdpy/dataform/features.py
+++ b/bdpy/dataform/features.py
@@ -26,8 +26,13 @@
     FeatureStore,
     HDF5FeatureStore,
     MatFeatureStore,
+    _check_size,
+    _default_iter_size,
     _determine_num_parallel,
+    _iter_array_slabs,
     _load_array_with_key,
+    _normalize_axis,
+    _partition,
     detect_format,
 )
 
@@ -256,7 +261,8 @@ def iter_chunks(
         The guarantee also assumes the requested labels all live in a single
         feature directory, which is the usual case. When they are spread across
         several `dpath` entries, no single store can stream them, so the
-        selection is read in full and then split.
+        selection is read in full and then split; without `size`, the slabs
+        then follow the byte budget, as on a store without on-disk chunking.
 
         Yields
         ------
@@ -288,22 +294,13 @@ def iter_chunks(
         else:
             labels = list(label)
 
-        # Validate here rather than in the store, so that both dispatch paths
-        # below reject the same inputs the same way. The multi-store fallback
-        # does not go through FeatureStore.iter_chunks and would otherwise
-        # silently yield nothing for size < 1, and raise IndexError rather than
-        # ValueError for an out-of-range axis.
-        ndim = len(self.shape(layer))
-        if axis < 0:
-            axis += ndim
-        if not 0 <= axis < ndim:
-            raise ValueError(
-                'axis {} is out of range for a {}-dimensional selection'.format(
-                    axis, ndim
-                )
-            )
-        if size is not None and size < 1:
-            raise ValueError('size must be positive, got {}'.format(size))
+        # Validate up front with the same helpers the stores use, so that both
+        # dispatch paths below reject the same inputs the same way, and the
+        # fallback rejects them before reading the selection.
+        selected_shape = (len(labels), *self.shape(layer)[1:])
+        axis = _normalize_axis(axis, selected_shape)
+        if size is not None:
+            _check_size(size)
 
         store = self.__store_for(labels)
         if store is not None:
@@ -311,16 +308,16 @@ def iter_chunks(
             return
 
         # Labels span several directories, so no single store can stream them.
-        # Fall back to slicing a full read, which still yields the same blocks.
+        # Fall back to slicing a full read; the default size follows the same
+        # byte budget as a store without on-disk chunking.
         features = self.__read(layer, labels)
-        length = features.shape[axis]
         if size is None:
-            size = length
-        for start in range(0, length, size):
-            sl = slice(start, min(start + size, length))
-            index: List[Any] = [slice(None)] * features.ndim
-            index[axis] = sl
-            yield sl, features[tuple(index)]
+            size = _default_iter_size(
+                features.shape, axis, features.dtype.itemsize, extent=None
+            )
+        yield from _iter_array_slabs(
+            features, axis, _partition(features.shape[axis], size)
+        )
 
     def statistic(self, statistic: str = 'mean', layer: Optional[str] = None):
         # NOTE: return type is ambiguous. currently, it is inferred as Unkown | Any
diff --git a/tests/dataform/test_feature_store.py b/tests/dataform/test_feature_store.py
index 9b34585..41bba40 100644
--- a/tests/dataform/test_feature_store.py
+++ b/tests/dataform/test_feature_store.py
@@ -10,15 +10,25 @@
 import tempfile
 import unittest
 import warnings
+from unittest import mock
 
 import h5py
 import numpy as np
 from numpy.testing import assert_array_equal
 
 from bdpy.dataform import Features, convert_features_to_hdf5, save_features
+from bdpy.dataform import _feature_store
 from bdpy.dataform._feature_store import (
+    FeatureStore,
     HDF5FeatureStore,
     MatFeatureStore,
+    _axis_index,
+    _check_size,
+    _default_iter_size,
+    _iter_array_slabs,
+    _normalize_axis,
+    _normalize_feature_slice,
+    _partition,
     detect_format,
 )
 
@@ -46,6 +56,103 @@
 ]
 
 
+# feature_slice specs outside basic forward indexing. Validation happens once,
+# in FeatureStore.read, so the table is checked against the validator itself;
+# TestFeatureSliceValidation checks that every backend goes through it.
+REJECTED_FEATURE_SLICES = [
+    ('negative step', np.s_[::-1]),
+    ('reversed with bounds', np.s_[16:4:-1]),
+    ('zero step', slice(None, None, 0)),
+    ('list', np.s_[[3, 1, 7]]),
+    ('ndarray', np.array([1, 2])),
+    ('range', range(3)),
+    ('bool', True),
+    ('numpy bool', np.True_),
+    ('newaxis in tuple', (np.newaxis, slice(None))),
+    ('two ellipses', (Ellipsis, Ellipsis)),
+    ('float start', slice(1.5, 3)),
+    ('float stop', slice(1, 3.5)),
+    ('float step', slice(None, None, 2.0)),
+    ('bool start', slice(True, 3)),
+]
+
+
+class TestIterationHelpers(unittest.TestCase):
+    """The index arithmetic behind iter_chunks, on plain in-memory values."""
+
+    def test_normalize_axis(self):
+        shape = (11, 24, 5, 5)
+        for axis, expected in ((0, 0), (3, 3), (-1, 3), (-4, 0)):
+            with self.subTest(axis=axis):
+                self.assertEqual(_normalize_axis(axis, shape), expected)
+
+    def test_normalize_axis_rejects_out_of_range(self):
+        shape = (11, 24, 5, 5)
+        for axis, message in (
+            (4, 'axis 4 is out of range for selected shape (11, 24, 5, 5)'),
+            (-5, 'axis -5 is out of range for selected shape (11, 24, 5, 5)'),
+        ):
+            with self.subTest(axis=axis):
+                with self.assertRaises(ValueError) as ctx:
+                    _normalize_axis(axis, shape)
+                self.assertEqual(str(ctx.exception), message)
+
+    def test_check_size(self):
+        _check_size(1)
+        for bad in (0, -1):
+            with self.subTest(size=bad):
+                with self.assertRaises(ValueError):
+                    _check_size(bad)
+
+    def test_partition(self):
+        cases = [
+            ('exact', 6, 3, [slice(0, 3), slice(3, 6)]),
+            ('partial last', 7, 3, [slice(0, 3), slice(3, 6), slice(6, 7)]),
+            ('size above length', 4, 10, [slice(0, 4)]),
+            ('size one', 3, 1, [slice(0, 1), slice(1, 2), slice(2, 3)]),
+            ('empty', 0, 3, []),
+        ]
+        for name, length, size, expected in cases:
+            with self.subTest(case=name):
+                self.assertEqual(list(_partition(length, size)), expected)
+
+    def test_axis_index(self):
+        sl = slice(2, 4)
+        self.assertEqual(_axis_index(0, sl), (sl,))
+        self.assertEqual(_axis_index(2, sl), (slice(None), slice(None), sl))
+
+    def test_default_size_follows_the_chunk_extent(self):
+        shape = (3, 10, 4)
+        self.assertEqual(_default_iter_size(shape, 1, 8, extent=4), 4)
+        # Never longer than the axis, never below one.
+        self.assertEqual(_default_iter_size(shape, 1, 8, extent=50), 10)
+        self.assertEqual(_default_iter_size(shape, 1, 8, extent=0), 1)
+
+    def test_default_size_without_extent_follows_the_byte_budget(self):
+        # One element along axis 1 spans 3 * 4 items of 8 bytes = 96 bytes.
+        shape = (3, 10, 4)
+        with mock.patch.object(_feature_store, 'DEFAULT_TARGET_CHUNK_BYTES', 300):
+            self.assertEqual(_default_iter_size(shape, 1, 8, extent=None), 3)
+        with mock.patch.object(_feature_store, 'DEFAULT_TARGET_CHUNK_BYTES', 10):
+            self.assertEqual(_default_iter_size(shape, 1, 8, extent=None), 1)
+        with mock.patch.object(_feature_store, 'DEFAULT_TARGET_CHUNK_BYTES', 10 ** 6):
+            self.assertEqual(_default_iter_size(shape, 1, 8, extent=None), 10)
+
+    def test_iter_array_slabs(self):
+        array = np.arange(2 * 7 * 3).reshape(2, 7, 3)
+        slabs = list(_iter_array_slabs(array, 1, _partition(7, 3)))
+        self.assertEqual([sl for sl, _ in slabs], list(_partition(7, 3)))
+        for sl, block in slabs:
+            assert_array_equal(block, array[:, sl])
+        assert_array_equal(np.concatenate([b for _, b in slabs], axis=1), array)
+
+    def test_rejected_feature_slices(self):
+        for name, spec in REJECTED_FEATURE_SLICES:
+            with self.subTest(index=name):
+                with self.assertRaises(ValueError):
+                    _normalize_feature_slice(spec)
+
+
 class _BackendPair(unittest.TestCase):
     """Builds the same features as a .mat tree and as chunked HDF5."""
 
@@ -167,43 +274,6 @@ def test_nothing_is_cached_between_calls(self):
                 list(self.store.iter_chunks('conv5', axis=1, size=4))
                 self.assertEqual(len(self.reads), 1)
 
-    def test_blocks_match_ground_truth_and_the_other_backend(self):
-        cases = [
-            ('conv5', None, 1),
-            ('conv5', None, 2),
-            ('conv5', None, 3),
-            ('conv5', None, -1),
-            ('conv5', None, -3),
-            ('fc8', None, 1),
-            ('conv5', ['img0009', 'img0001'], 1),
-            ('conv5', ['img0009', 'img0001', 'img0009'], 2),
-            ('conv5', ['img0003'], 3),
-        ]
-        for layer, label, axis in cases:
-            for size in (1, 3, 1000, None):
-                with self.subTest(layer=layer, label=label, axis=axis, size=size):
-                    rows = (
-                        slice(None) if label is None
-                        else [LABELS.index(s) for s in label]
-                    )
-                    expected = self.stacked[layer][rows]
-                    norm = axis if axis >= 0 else axis + expected.ndim
-                    from_mat = list(self.from_mat.iter_chunks(
-                        layer, label=label, axis=axis, size=size))
-                    for sl, block in from_mat:
-                        assert_array_equal(
-                            block, expected[(slice(None),) * norm + (sl,)]
-                        )
-                    # The two backends need not agree on slab boundaries -- the
-                    # HDF5 default follows the on-disk chunk extent -- so the
-                    # reassembled layer is what has to match.
-                    from_h5 = list(self.from_h5.iter_chunks(
-                        layer, label=label, axis=axis, size=size))
-                    assert_array_equal(
-                        np.concatenate([b for _, b in from_mat], axis=norm),
-                        np.concatenate([b for _, b in from_h5], axis=norm),
-                    )
-
     def test_writing_to_one_block_does_not_disturb_another(self):
         # Slabs are views into one array now, but they partition the axis, so
         # they cannot alias each other.
@@ -211,10 +281,10 @@ def test_writing_to_one_block_does_not_disturb_another(self):
         blocks[0][1][:] = 0
         assert_array_equal(blocks[1][1], self.stacked['conv5'][:, 4:8])
 
-    def test_validation_precedes_any_read(self):
+    def test_validation_precedes_selection_read(self):
         cases = (
             ({'axis': 9}, 'axis 9 is out of range for selected shape (11, 24, 5, 5)'),
-            ({'axis': -9}, 'axis -5 is out of range for selected shape (11, 24, 5, 5)'),
+            ({'axis': -9}, 'axis -9 is out of range for selected shape (11, 24, 5, 5)'),
             ({'size': 0}, 'size must be positive, got 0'),
         )
         for kwargs, message in cases:
@@ -231,25 +301,15 @@ def test_validation_precedes_any_read(self):
 
 
 class TestIterChunks(_BackendPair):
-    def test_reassembles_to_a_full_read(self):
-        for axis in (0, 1, 2, 3):
-            for features in (self.from_mat, self.from_h5):
-                with self.subTest(axis=axis, backend=type(features).__name__):
-                    blocks = list(features.iter_chunks('conv5', axis=axis, size=4))
-                    joined = np.concatenate([b for _, b in blocks], axis=axis)
-                    assert_array_equal(joined, self.stacked['conv5'])
-
-    def test_yielded_slices_address_the_right_block(self):
-        full = self.from_h5.get('conv5')
-        for sl, block in self.from_h5.iter_chunks('conv5', axis=1, size=5):
-            assert_array_equal(block, full[:, sl])
-
-    def test_slices_cover_the_axis_without_overlap(self):
-        length = self.from_h5.shape('conv5')[1]
-        covered = []
-        for sl, _ in self.from_h5.iter_chunks('conv5', axis=1, size=5):
-            covered.extend(range(*sl.indices(length)))
-        self.assertEqual(covered, list(range(length)))
+    def test_sample_axis_reassembles_to_a_full_read(self):
+        # The sample axis is iterated by label, not through a feature slice, so
+        # it is a separate path from the one test_feature_axis_blocks_match_ground_truth
+        # covers.
+        for features in (self.from_mat, self.from_h5):
+            with self.subTest(store=type(features._Features__stores[0]).__name__):
+                blocks = list(features.iter_chunks('conv5', axis=0, size=4))
+                joined = np.concatenate([b for _, b in blocks], axis=0)
+                assert_array_equal(joined, self.stacked['conv5'])
 
     def test_blocks_can_be_sliced_by_the_caller(self):
         # iter_chunks deliberately takes no feature_slice; slicing the blocks as
@@ -259,12 +319,6 @@ def test_blocks_can_be_sliced_by_the_caller(self):
             wanted[:, sl] = block[:, :, 1:2]
         assert_array_equal(wanted, self.stacked['conv5'][:, :, 1:2])
 
-    def test_restricts_to_labels(self):
-        labels = ['img0009', 'img0001']
-        blocks = list(self.from_h5.iter_chunks('conv5', label=labels, axis=1, size=5))
-        joined = np.concatenate([b for _, b in blocks], axis=1)
-        assert_array_equal(joined, self.stacked['conv5'][[9, 1]])
-
     def test_default_size_is_chunk_aligned(self):
         store = HDF5FeatureStore(self.h5dir)
         extent = store.chunk_extent('conv5', 1)
@@ -300,6 +354,55 @@ def test_multi_directory_validates_like_a_single_one(self):
                 with self.assertRaises(ValueError):
                     list(self.from_h5.iter_chunks('conv5', **kwargs))
 
+    def test_feature_axis_blocks_match_ground_truth(self):
+        # What each backend owns is cutting a feature-axis slab out of what it
+        # reads; the axis and size arithmetic is checked in
+        # TestIterationHelpers. Each size leaves a partial last slab, and the
+        # labels repeat and are out of order.
+        labels = ['img0009', 'img0001', 'img0009']
+        cases = [('conv5', 1, 5), ('conv5', 2, 2), ('conv5', 3, 2), ('fc8', 1, 7)]
+        for features in (self.from_mat, self.from_h5):
+            for layer, axis, size in cases:
+                expected = self.stacked[layer][[LABELS.index(s) for s in labels]]
+                with self.subTest(store=type(features._Features__stores[0]).__name__,
+                                  layer=layer, axis=axis):
+                    blocks = list(features.iter_chunks(
+                        layer, label=labels, axis=axis, size=size))
+                    self.assertEqual(
+                        [sl for sl, _ in blocks],
+                        list(_partition(expected.shape[axis], size)),
+                    )
+                    for sl, block in blocks:
+                        assert_array_equal(block, expected[_axis_index(axis, sl)])
+
+    def test_multi_directory_default_size_follows_the_byte_budget(self):
+        # Spanning two directories takes the in-memory fallback, which used to
+        # yield the whole axis as one slab when size was None. It now splits
+        # the way a store without on-disk chunking does. The budget is shrunk
+        # so that the small fixture spans several slabs.
+        other_dir = tempfile.TemporaryDirectory()
+        self.addCleanup(other_dir.cleanup)
+        other = other_dir.name
+        rng = np.random.default_rng(1)
+        for layer, shape in zip(LAYERS, SHAPES):
+            data = rng.random((1, *shape[1:]))
+            save_features(os.path.join(other, layer + '.h5'), data, ['other0000'])
+            if layer == 'conv5':
+                other_data = data
+        spread = Features([self.matdir, other])
+        both = [LABELS[0], 'other0000']
+        expected = np.concatenate([self.stacked['conv5'][:1], other_data])
+
+        with mock.patch.object(_feature_store, 'DEFAULT_TARGET_CHUNK_BYTES', 1000):
+            from_spread = list(spread.iter_chunks('conv5', label=both, axis=1))
+            # The .mat store has no chunk extent, so it uses the same budget.
+            from_store = list(self.from_mat.iter_chunks('conv5', label=LABELS[:2], axis=1))
+
+        self.assertGreater(len(from_spread), 1)
+        self.assertEqual([sl for sl, _ in from_spread], [sl for sl, _ in from_store])
+        for sl, block in from_spread:
+            assert_array_equal(block, expected[:, sl])
+
     def test_negative_axis_is_accepted(self):
         by_negative = list(self.from_h5.iter_chunks('conv5', axis=-3, size=4))
         by_positive = list(self.from_h5.iter_chunks('conv5', axis=1, size=4))
@@ -308,10 +411,6 @@ def test_negative_axis_is_accepted(self):
             self.assertEqual(sl_a, sl_b)
             assert_array_equal(a, b)
 
-    def test_rejects_bad_axis(self):
-        with self.assertRaises(ValueError):
-            list(self.from_h5.iter_chunks('conv5', axis=9))
-
 
 class TestPartialReads(_BackendPair):
     def test_slice_is_pushed_down_to_h5py(self):
@@ -387,23 +486,6 @@ class TestFeatureSliceValidation(_BackendPair):
         ('open slice', np.s_[:]),
     ]
 
-    rejected = [
-        ('negative step', np.s_[::-1]),
-        ('reversed with bounds', np.s_[16:4:-1]),
-        ('zero step', slice(None, None, 0)),
-        ('list', np.s_[[3, 1, 7]]),
-        ('ndarray', np.array([1, 2])),
-        ('range', range(3)),
-        ('bool', True),
-        ('numpy bool', np.True_),
-        ('newaxis in tuple', (np.newaxis, slice(None))),
-        ('two ellipses', (Ellipsis, Ellipsis)),
-        ('float start', slice(1.5, 3)),
-        ('float stop', slice(1, 3.5)),
-        ('float step', slice(None, None, 2.0)),
-        ('bool start', slice(True, 3)),
-    ]
-
     def test_accepted_agree_across_backends_and_numpy(self):
         for name, spec in self.accepted:
             with self.subTest(index=name):
@@ -414,15 +496,23 @@ def test_accepted_agree_across_backends_and_numpy(self):
                 assert_array_equal(from_mat, expected)
                 assert_array_equal(from_h5, expected)
 
-    def test_rejected_by_both_backends(self):
-        for name, spec in self.rejected:
-            with self.subTest(index=name):
-                for backend, features in (
-                    ('mat', self.from_mat), ('hdf5', self.from_h5),
-                ):
-                    with self.subTest(backend=backend):
-                        with self.assertRaises(ValueError):
-                            features.get('conv5', feature_slice=spec)
+    def test_backends_do_not_override_read(self):
+        # FeatureStore.read is where feature_slice is validated; a backend that
+        # overrode it could skip validation, which no spy on _read would see.
+        for store_cls in (MatFeatureStore, HDF5FeatureStore):
+            with self.subTest(store=store_cls.__name__):
+                self.assertIs(store_cls.read, FeatureStore.read)
+
+    def test_rejected_by_both_backends_before_backend_read(self):
+        # The full table of rejected specs is checked in TestIterationHelpers;
+        # here each backend only has to route through the shared validation,
+        # and reject before the backend's _read runs.
+        for store in (MatFeatureStore(self.matdir), HDF5FeatureStore(self.h5dir)):
+            with self.subTest(store=type(store).__name__):
+                with mock.patch.object(store, '_read') as backend_read:
+                    with self.assertRaises(ValueError):
+                        store.read('conv5', feature_slice=np.s_[::-1])
+                backend_read.assert_not_called()
 
     def test_bool_is_not_read_as_an_integer(self):
         # isinstance(True, int) is True in Python, so a naive integer check

From a9331c7538f9e3d3eeb1ddcbfbc04b892d646e7e Mon Sep 17 00:00:00 2001
From: Kenya Otsuka <kenya.otsuka.0720+github@gmail.com>
Date: Tue, 6 Oct 2026 11:37:09 +0000
Subject: [PATCH 11/11] fix(dataform): never clobber a file published while a
 writer is staging

overwrite=False was enforced only when the writer was created; close()
then published with os.replace, which replaces whatever is at the target.
If another process (or save_features) published the same path while this
writer was still staging, it was silently overwritten.

Without overwrite, publish with os.link, which fails if the target
exists, so checking for it and publishing become one atomic step; the
scratch file is removed afterwards. With overwrite, keep os.replace. On a
filesystem that cannot hard-link (EPERM, ENOTSUP/EOPNOTSUPP, ENOSYS),
fall back to checking right before the rename -- a short window remains
there, as documented, but no longer the whole duration of the write.

When close() finds the target taken, the writer discards its scratch
file, is left aborted, and raises FileExistsError. convert_features_to_hdf5
treats that as its existing skip-if-exists case: a layer another process
published meanwhile is skipped, and the next layer is converted.
---
 bdpy/dataform/feature_hdf5.py       | 125 +++++++++++++++++++++-------
 tests/dataform/test_feature_hdf5.py |  66 +++++++++++++++
 2 files changed, 160 insertions(+), 31 deletions(-)

diff --git a/bdpy/dataform/feature_hdf5.py b/bdpy/dataform/feature_hdf5.py
index 3d01560..02fb8af 100644
--- a/bdpy/dataform/feature_hdf5.py
+++ b/bdpy/dataform/feature_hdf5.py
@@ -43,6 +43,7 @@
     convert_features_to_hdf5('features_mat', 'features_h5')
 """
 
+import errno
 import os
 import uuid
 import weakref
@@ -111,27 +112,68 @@ def _staging_path(path: str) -> str:
 def _prepare_target(path: str, overwrite: bool) -> str:
     """Refuse to clobber, then name a scratch file next to `path`.
 
-    Writing goes to a scratch file in the *same directory*, so the final
-    :func:`os.replace` is an atomic rename on the same filesystem: either the
-    finished file appears at `path`, or nothing does. A half-written file must
-    never be left where a reader -- or the converter's skip-if-exists check --
-    would take it for a complete one.
+    Writing goes to a scratch file in the *same directory*, so publishing it
+    (see :func:`_promote`) is a single atomic step on the same filesystem:
+    either the finished file appears at `path`, or nothing does. A half-written
+    file must never be left where a reader -- or the converter's skip-if-exists
+    check -- would take it for a complete one.
+
+    This check fails early, before any work is done; :func:`_promote` repeats
+    it when publishing, since another writer may publish `path` meanwhile.
 
     The caller opens the returned path with mode ``"x"``, which both creates it
     exclusively and lets HDF5 apply the process umask, so published files keep
     the permissions they had before staging existed.
     """
     if os.path.exists(path) and not overwrite:
-        raise FileExistsError(
-            "{} already exists. Pass overwrite=True to replace it.".format(path)
-        )
+        raise _exists_error(path)
     _makedirs_for(path)
     return _staging_path(path)
 
 
-def _promote(tmp_path: str, path: str) -> None:
-    """Move the finished scratch file into place (atomic on one filesystem)."""
-    os.replace(tmp_path, path)
+def _exists_error(path: str) -> FileExistsError:
+    return FileExistsError(
+        "{} already exists. Pass overwrite=True to replace it.".format(path)
+    )
+
+
+#: errno values with which os.link reports that the filesystem cannot hard-link.
+_NO_HARDLINK_ERRNOS = frozenset(
+    {errno.EPERM, errno.ENOTSUP, errno.EOPNOTSUPP, errno.ENOSYS}
+)
+
+
+def _promote(tmp_path: str, path: str, overwrite: bool) -> None:
+    """Publish the finished scratch file at `path`.
+
+    With `overwrite`, an atomic rename replaces whatever is there. Without it,
+    the file is published with :func:`os.link`, which fails if `path` exists:
+    checking for a file and publishing ours are one atomic step, so a file that
+    another writer published while this one was staging is never clobbered.
+
+    On a filesystem that cannot hard-link, this falls back to checking for
+    `path` and then renaming, which leaves a short window between the two.
+
+    Raises
+    ------
+    FileExistsError
+        If `path` exists and `overwrite` is False.
+    """
+    if overwrite:
+        os.replace(tmp_path, path)
+        return
+    try:
+        os.link(tmp_path, path)
+    except FileExistsError:
+        raise _exists_error(path) from None
+    except OSError as exc:
+        if exc.errno not in _NO_HARDLINK_ERRNOS:
+            raise
+        if os.path.exists(path):
+            raise _exists_error(path) from None
+        os.replace(tmp_path, path)
+    else:
+        _discard(tmp_path)
 
 
 def _discard(tmp_path: str) -> None:
@@ -157,8 +199,9 @@ def save_features(
     Parameters
     ----------
     path : str
-        Output file. An existing file is kept unless `overwrite` is True. The
-        write is atomic: on failure nothing is left at `path`.
+        Output file. An existing file -- including one another writer
+        publishes meanwhile -- is kept unless `overwrite` is True. The write is
+        atomic: on failure nothing is left at `path`.
     features : numpy.ndarray
         Feature array of shape ``(n_samples, *feature_shape)``.
     labels : sequence of str
@@ -182,7 +225,7 @@ def save_features(
         If `labels` does not have one entry per sample, or `features` has fewer
         than two axes.
     FileExistsError
-        If `path` exists and `overwrite` is False.
+        If `path` exists or is published meanwhile, and `overwrite` is False.
     """
     features = np.asarray(features)
     if features.ndim < 2:
@@ -227,8 +270,11 @@ class FeatureWriter:
     Parameters
     ----------
     path : str
-        Output file. An existing file is kept unless `overwrite` is True, and
-        nothing is written there until :meth:`close` succeeds.
+        Output file. An existing file -- including one another writer
+        publishes meanwhile -- is kept unless `overwrite` is True, and nothing
+        is written there until :meth:`close` succeeds. On a filesystem without
+        hard links, a short window remains in which a file published at the
+        same moment can still be replaced.
     feature_shape : sequence of int
         Shape of a single sample's features, without the sample axis.
     dtype : numpy.dtype
@@ -247,7 +293,8 @@ class FeatureWriter:
     Raises
     ------
     FileExistsError
-        If `path` exists and `overwrite` is False.
+        If `path` exists, or (from :meth:`close`) was published meanwhile,
+        and `overwrite` is False.
 
     Examples
     --------
@@ -286,6 +333,7 @@ def __init__(
         # Fail before doing any work if the target is occupied, then write to
         # a scratch file so `path` stays untouched until close() succeeds.
         self._path = path
+        self._overwrite = overwrite
         self._aborted = False
         self._tmp_path: Optional[str] = _prepare_target(path, overwrite)
         self._file: Optional[h5py.File] = h5py.File(self._tmp_path, "x")
@@ -411,6 +459,8 @@ def close(self) -> None:
         ------
         RuntimeError
             If the writer was already aborted.
+        FileExistsError
+            If the target was published meanwhile; the writer is then aborted.
         """
         if self._aborted:
             raise RuntimeError(
@@ -422,7 +472,11 @@ def close(self) -> None:
             self._file.close()
             self._file = None
         if self._tmp_path is not None:
-            _promote(self._tmp_path, self._path)
+            try:
+                _promote(self._tmp_path, self._path, self._overwrite)
+            except FileExistsError:
+                self.abort()
+                raise
             self._tmp_path = None
             self._finalizer.detach()  # the scratch file is now the output
 
@@ -494,7 +548,8 @@ def convert_features_to_hdf5(
         Variable name inside the legacy files (default: ``'feat'``).
     overwrite : bool, optional
         Overwrite an existing output file instead of skipping it
-        (default: False).
+        (default: False). Without it, a layer that another process publishes
+        while this one is converting it is skipped as well.
     batch_size : int, optional
         Number of stimulus files read per batch. Must be positive.
     target_chunk_bytes : int, optional
@@ -539,18 +594,26 @@ def convert_features_to_hdf5(
         # partial file instead of publishing it. A published partial file would
         # be taken for a finished one by the skip check above, and the layer
         # would stay silently truncated across re-runs.
-        with FeatureWriter(
-            out_path,
-            feature_shape=full_shape[1:],
-            dtype=store.dtype(layer),
-            n_samples=full_shape[0],
-            target_chunk_bytes=target_chunk_bytes,
-            compression=compression,
-            overwrite=overwrite,
-        ) as writer:
-            for start in range(0, len(all_labels), batch_size):
-                batch = all_labels[start:start + batch_size]
-                writer.extend(store.read(layer, batch), batch)
+        try:
+            with FeatureWriter(
+                out_path,
+                feature_shape=full_shape[1:],
+                dtype=store.dtype(layer),
+                n_samples=full_shape[0],
+                target_chunk_bytes=target_chunk_bytes,
+                compression=compression,
+                overwrite=overwrite,
+            ) as writer:
+                for start in range(0, len(all_labels), batch_size):
+                    batch = all_labels[start:start + batch_size]
+                    writer.extend(store.read(layer, batch), batch)
+        except FileExistsError:
+            # Only reachable without overwrite: another process published this
+            # layer while we were converting it. Treat it like the check above;
+            # the writer has already discarded our scratch file.
+            if verbose:
+                print("{} already exists. Skipped.".format(out_path))
+            continue
 
         if verbose:
             print("Saved {} ({} samples).".format(out_path, full_shape[0]))
diff --git a/tests/dataform/test_feature_hdf5.py b/tests/dataform/test_feature_hdf5.py
index 8ba70da..8a30b1c 100644
--- a/tests/dataform/test_feature_hdf5.py
+++ b/tests/dataform/test_feature_hdf5.py
@@ -1,3 +1,4 @@
+import errno
 import gc
 import glob
 import os
@@ -5,6 +6,7 @@
 import tempfile
 import unittest
 import warnings
+from unittest import mock
 
 import h5py
 import numpy as np
@@ -28,6 +30,10 @@
 from .test_features import prepare_mat_features
 
 
+def _scratch_files(directory):
+    return [n for n in os.listdir(directory) if n.endswith('.partial')]
+
+
 def _current_umask():
     mask = os.umask(0)
     os.umask(mask)
@@ -207,6 +213,27 @@ def test_existing_file_is_skipped_unless_overwrite(self):
         convert_features_to_hdf5(self.matdir, self.h5dir, overwrite=True)
         self.assertGreater(os.path.getmtime(marker), mtime - 100)
 
+    def test_layer_published_by_another_process_meanwhile_is_skipped(self):
+        # Two converters on the same destination: the other one publishes
+        # conv5 while this one is still reading it. This one must keep the
+        # other's file, discard its own, and go on to the next layer.
+        other = self.stacked['conv5'][::-1].copy()
+        out_path = os.path.join(self.h5dir, 'conv5.h5')
+        original_read = MatFeatureStore.read
+
+        def read_then_race(store, layer, labels=None, feature_slice=None):
+            if layer == 'conv5' and not os.path.exists(out_path):
+                save_features(out_path, other, self.labels)
+            return original_read(store, layer, labels, feature_slice)
+
+        with mock.patch.object(MatFeatureStore, 'read', read_then_race):
+            convert_features_to_hdf5(self.matdir, self.h5dir)
+
+        converted = Features(self.h5dir)
+        assert_array_equal(converted.get('conv5'), other)
+        assert_array_equal(converted.get('fc8'), self.stacked['fc8'])
+        self.assertEqual(_scratch_files(self.h5dir), [])
+
     def test_small_batches_produce_the_same_file(self):
         # Batching is an implementation detail; it must not affect the result.
         convert_features_to_hdf5(self.matdir, self.h5dir, batch_size=3)
@@ -510,6 +537,45 @@ def test_dropped_writer_does_not_leak_a_scratch_file(self):
         gc.collect()
         self.assertEqual(self._leftovers(self.tmpdir.name), ['mat'])
 
+    # --- another writer publishes first ---------------------------------
+
+    def test_close_does_not_clobber_a_file_published_meanwhile(self):
+        # The overwrite check in __init__ is not enough on its own: another
+        # writer can publish the same path while this one is still staging.
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.extend(self.data, self.labels)
+        save_features(self.path, self.data * 2, self.labels)
+        with self.assertRaises(FileExistsError):
+            writer.close()
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data * 2
+        )
+        self.assertEqual(_scratch_files(self.tmpdir.name), [])
+        # The failed publish leaves the writer aborted, not half-closed.
+        with self.assertRaises(RuntimeError):
+            writer.close()
+
+    def test_publishes_without_hard_links(self):
+        # A filesystem that cannot hard-link falls back to check-then-rename.
+        with mock.patch('os.link', side_effect=OSError(errno.EPERM, 'no links')):
+            save_features(self.path, self.data, self.labels)
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+        self.assertEqual(_scratch_files(self.tmpdir.name), [])
+
+    def test_without_hard_links_a_file_published_meanwhile_is_kept(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.extend(self.data, self.labels)
+        save_features(self.path, self.data * 2, self.labels)
+        with mock.patch('os.link', side_effect=OSError(errno.EPERM, 'no links')):
+            with self.assertRaises(FileExistsError):
+                writer.close()
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data * 2
+        )
+        self.assertEqual(_scratch_files(self.tmpdir.name), [])
+
     def test_close_after_abort_is_refused(self):
         # Silently doing nothing would let a caller believe it published.
         writer = FeatureWriter(self.path, (8,), np.float32)