symmetry: port the OME-Zarr symmetry test (readOmeZarr) to Python

.github/workflows/test-symmetry.yml

@@ -97,11 +97,14 @@ jobs:
         with:
           python-version: "3.12"
 
-      - name: Install NDR-python (with dev extras)
+      # formats carries zarr/numcodecs (OME-Zarr) and the other format
+      # libraries the symmetry tests read with; under strict mode a missing
+      # one is a failure, not a skip.
+      - name: Install NDR-python (with dev and formats extras)
         working-directory: NDR-python
         run: |
           python -m pip install --upgrade pip
-          pip install -e ".[dev]"
+          pip install -e ".[dev,formats]"
 
       # -- Stage 1: MATLAB makeArtifacts ------------------------------------
       - name: "Stage 1: MATLAB makeArtifacts"
@@ -178,21 +181,10 @@ jobs:
             import matlab.unittest.TestSuite
             import matlab.unittest.TestRunner
 
+            % Every MATLAB symmetry test has a Python half, so every case here,
+            % pythonArtifacts included, must run (VH-Lab/NDR-python#36).
             suite = TestSuite.fromPackage("ndr.symmetry.readArtifacts", "IncludingSubpackages", true);
 
-            % KNOWN GAPS: MATLAB symmetry tests with no Python half, whose
-            % pythonArtifacts case therefore cannot run. Deselected by name and
-            % listed in the log, never left to come back Incomplete. Remove an
-            % entry when its Python half lands.
-            %   readOmeZarr -- NDR-python has no OME-Zarr symmetry test yet.
-            noPythonHalf = ["readOmeZarr"];
-            gap = contains({suite.Name}, "SourceType=pythonArtifacts") & ...
-                contains({suite.Name}, noPythonHalf);
-            for k = find(gap)
-                fprintf("Deselected (no Python half yet): %s\n", suite(k).Name);
-            end
-            suite = suite(~gap);
-
             fprintf("\n=== MATLAB readArtifacts: %d test(s) ===\n", numel(suite));
             assert(~isempty(suite), "No MATLAB readArtifacts tests found.")
 

tests/_omezarr_fixture.py

@@ -6,9 +6,14 @@
 ``ferret-lightsheet/formats/OurUsualZarrFormat.md``: one shared level-0
 array plus separate ``mean`` and ``max`` downsample chains. Without
 ``with_chunks`` only metadata files are written; with ``with_chunks`` a
-seeded uint16 volume is written through the ``zarr`` library so
-chunks exist on disk and their ground-truth downsamples come back for
-comparison.
+seeded uint16 volume is written so chunks exist on disk and their
+ground-truth downsamples come back for comparison.
+
+The chunks are uncompressed by default. ``codec="blosc-zstd"`` instead writes
+every array through zarr-python with Blosc+Zstd, the codec the lab's own data
+uses and the one ``makeExampleFixture.m`` writes. The OME-Zarr symmetry test
+needs that: MATLAB's reader hand-rolls the Blosc/Zstd path, so only a fixture
+compressed by zarr-python itself checks it against real zarr-python output.
 """
 
 from __future__ import annotations
@@ -122,13 +127,47 @@ def _write_chunks(array_dir: Path, vol: np.ndarray, chunks: tuple[int, ...]) ->
                     (array_dir / key).write_bytes(chunk.tobytes(order="C"))
 
 
+def _write_zarr_array(array_dir: Path, vol: np.ndarray, chunks: tuple[int, ...]) -> None:
+    """Write ``vol`` as a zarr v2 array compressed with Blosc+Zstd, via zarr-python.
+
+    Replaces the ``.zarray`` written by :func:`_write_zarray` with zarr-python's
+    own, so the metadata and the chunks are both exactly what zarr-python
+    produces. Same codec settings as ``makeExampleFixture.m``: Blosc container,
+    Zstd, level 5, byte shuffle.
+    """
+    import zarr
+    from numcodecs import Blosc
+
+    kwargs = {
+        "mode": "w",
+        "shape": vol.shape,
+        "chunks": chunks,
+        "dtype": "<u2",
+        "compressor": Blosc(cname="zstd", clevel=5, shuffle=Blosc.SHUFFLE),
+        "fill_value": 0,
+    }
+    try:
+        arr = zarr.open_array(str(array_dir), zarr_format=2, **kwargs)
+    except TypeError:  # zarr 2.x: no zarr_format keyword, and v2 is the default
+        arr = zarr.open_array(str(array_dir), **kwargs)
+    arr[...] = vol
+
+
 def make_example_fixture(
     parent_dir: str | Path,
     *,
     with_chunks: bool = False,
     seed: int = 20260906,
+    codec: str = "none",
 ) -> tuple[Path, GroundTruth]:
-    """Write a dual-pyramid OME-Zarr store and return its root and ground truth."""
+    """Write a dual-pyramid OME-Zarr store and return its root and ground truth.
+
+    ``codec`` is ``"none"`` (uncompressed chunks written by hand, the default)
+    or ``"blosc-zstd"`` (arrays written by zarr-python with Blosc+Zstd). It
+    only matters with ``with_chunks``.
+    """
+    if codec not in ("none", "blosc-zstd"):
+        raise ValueError(f'codec must be "none" or "blosc-zstd", not {codec!r}.')
     parent = Path(parent_dir)
     parent.mkdir(parents=True, exist_ok=True)
     fixture_dir = parent / "example.zarr"
@@ -205,11 +244,12 @@ def make_example_fixture(
     mean2 = _downsample_block(mean1, "mean")
     max2 = _downsample_block(max1, "max")
 
-    _write_chunks(fixture_dir / "0", lvl0, chunks0)
-    _write_chunks(fixture_dir / "mean" / "1", mean1, chunks1)
-    _write_chunks(fixture_dir / "max" / "1", max1, chunks1)
-    _write_chunks(fixture_dir / "mean" / "2", mean2, chunks2)
-    _write_chunks(fixture_dir / "max" / "2", max2, chunks2)
+    write = _write_zarr_array if codec == "blosc-zstd" else _write_chunks
+    write(fixture_dir / "0", lvl0, chunks0)
+    write(fixture_dir / "mean" / "1", mean1, chunks1)
+    write(fixture_dir / "max" / "1", max1, chunks1)
+    write(fixture_dir / "mean" / "2", mean2, chunks2)
+    write(fixture_dir / "max" / "2", max2, chunks2)
 
     gt.level0 = lvl0
     gt.mean_level1 = mean1

tests/symmetry/conftest.py

@@ -101,3 +101,21 @@ def pytest_runtest_makereport(item, call):
         "comparison to make this green."
     )
     return report
+
+
+def matlab_sha256(array) -> str:
+    """SHA-256 of an array's bytes in MATLAB's ``typecast(x(:), 'uint8')`` order.
+
+    Both ports' readers index arrays the same way (MATLAB ``x(i,j,k)`` is
+    numpy ``x[i-1, j-1, k-1]``), but ``x(:)`` linearizes in column-major
+    (Fortran) order, and MATLAB hosts are little-endian. So this flattens in
+    Fortran order and encodes little-endian before hashing, which gives the
+    digest a MATLAB symmetry test computes for the same array.
+    """
+    import hashlib
+
+    import numpy as np
+
+    a = np.asarray(array)
+    a = a.astype(a.dtype.newbyteorder("<"), copy=False)
+    return hashlib.sha256(a.ravel(order="F").tobytes()).hexdigest()

tests/symmetry/make_artifacts/format/test_read_omezarr.py

@@ -0,0 +1,94 @@
+"""Generate symmetry artifacts for the OME-Zarr format layer (readArray).
+
+Zarr is lossless, so if MATLAB and Python decode the same chunks to
+different values, one of them is wrong. This writes a dual-pyramid OME-Zarr
+store with zarr-python itself (zarr v2, Blosc+Zstd: the codec the lab's data
+uses), then records the SHA-256 of what ndr.format.omezarr.readArray returns
+for a fixed set of regions. MATLAB's readArtifacts/readOmeZarr reads the same
+store with its hand-rolled Blosc/Zstd reader and must reproduce every digest.
+
+The regions and the artifact layout (example.zarr, metadata.json,
+checksums.json) are those of the MATLAB half, so either port can read the
+other's. Digests are computed over the bytes in MATLAB's typecast(x(:)) order
+(see conftest.matlab_sha256).
+
+Mirrors tools/tests/+ndr/+symmetry/+makeArtifacts/+format/readOmeZarr.m
+"""
+
+import json
+import shutil
+
+import numpy as np
+import pytest
+
+from tests.symmetry.conftest import PYTHON_ARTIFACTS, matlab_sha256
+
+ARTIFACT_DIR = PYTHON_ARTIFACTS / "format" / "readOmeZarr" / "testReadOmeZarrArtifacts"
+
+SEED = 20260906
+
+# (name, pyramid, level, region) -- the MATLAB half's regionSpecs. A region is
+# [start, stop] per axis (c, z, y, x), 1-based and inclusive; None is the
+# whole array. Chosen to cover a whole level, one chunk, a chunk straddle,
+# and the x-axis edge chunk (10 voxels in chunks of 4).
+REGIONS = [
+    ("whole_mean_1", "mean", 1, None),
+    ("whole_max_1", "max", 1, None),
+    ("whole_mean_2", "mean", 2, None),
+    ("whole_max_2", "max", 2, None),
+    ("inside_one_chunk", "mean", 1, [[1, 1, 1, 1], [1, 2, 3, 3]]),
+    ("straddle", "mean", 1, [[1, 3, 3, 3], [1, 6, 6, 6]]),
+    ("reach_edge_x", "mean", 1, [[1, 1, 1, 7], [1, 8, 8, 10]]),
+    ("inside_edge_x", "mean", 1, [[1, 1, 1, 9], [1, 4, 4, 10]]),
+]
+
+
+class TestReadOmeZarr:
+    @pytest.fixture(autouse=True)
+    def _setup(self):
+        pytest.importorskip("zarr")
+        pytest.importorskip("numcodecs")
+
+    def test_read_ome_zarr_artifacts(self):
+        from ndr.format.omezarr.readArray import readArray
+        from tests._omezarr_fixture import make_example_fixture
+
+        if ARTIFACT_DIR.exists():
+            shutil.rmtree(ARTIFACT_DIR)
+        ARTIFACT_DIR.mkdir(parents=True)
+
+        fixture_dir, gt = make_example_fixture(
+            ARTIFACT_DIR, with_chunks=True, seed=SEED, codec="blosc-zstd"
+        )
+
+        # The reader must return the ground truth before its digests are
+        # worth publishing: a digest of wrong data would only test agreement.
+        np.testing.assert_array_equal(readArray(str(fixture_dir), "mean", 1), gt.level0)
+        np.testing.assert_array_equal(readArray(str(fixture_dir), "mean", 2), gt.mean_level1)
+        np.testing.assert_array_equal(readArray(str(fixture_dir), "max", 2), gt.max_level1)
+
+        metadata = {
+            "fixtureName": fixture_dir.name,
+            "level0Shape": list(gt.shape0),
+            "level0ChunkShape": list(gt.chunk_shape0),
+            "dtype": gt.dtype,
+            "seed": SEED,
+        }
+        (ARTIFACT_DIR / "metadata.json").write_text(json.dumps(metadata, indent=2))
+
+        regions = []
+        for name, pyramid, level, region in REGIONS:
+            if region is None:
+                got = readArray(str(fixture_dir), pyramid, level)
+            else:
+                got = readArray(str(fixture_dir), pyramid, level, Region=np.array(region))
+            regions.append(
+                {
+                    "name": name,
+                    "pyramid": pyramid,
+                    "level": level,
+                    "region": region if region is not None else [],
+                    "sha256": matlab_sha256(got),
+                }
+            )
+        (ARTIFACT_DIR / "checksums.json").write_text(json.dumps({"regions": regions}, indent=2))

tests/symmetry/read_artifacts/format/test_read_omezarr.py

@@ -0,0 +1,73 @@
+"""Read and verify symmetry artifacts for the OME-Zarr format layer (readArray).
+
+Parameterized over both artifact sources. Each source ships its own
+example.zarr inside its artifact directory, plus the SHA-256 of what its own
+readArray returned for each region. This reads that store with the Python
+readArray and checks every digest: a mismatch means the two ports decode the
+same lossless bytes differently, so one of them is wrong.
+
+Skips when the artifact directory for a source does not exist; under
+NDR_SYMMETRY_CHECK_STRICT (CI) that skip fails, see tests/symmetry/conftest.py.
+
+Mirrors tools/tests/+ndr/+symmetry/+readArtifacts/+format/readOmeZarr.m
+"""
+
+import json
+
+import numpy as np
+import pytest
+
+from tests.symmetry.conftest import SOURCE_TYPES, SYMMETRY_BASE, matlab_sha256
+
+
+@pytest.fixture(params=SOURCE_TYPES)
+def source_type(request):
+    return request.param
+
+
+def _artifact_dir(source_type):
+    return SYMMETRY_BASE / source_type / "format" / "readOmeZarr" / "testReadOmeZarrArtifacts"
+
+
+class TestReadOmeZarr:
+    def test_checksums(self, source_type):
+        pytest.importorskip("zarr")
+        from ndr.format.omezarr.readArray import readArray
+
+        artifact_dir = _artifact_dir(source_type)
+        if not artifact_dir.exists():
+            pytest.skip(f"No artifacts from {source_type}")
+
+        metadata = json.loads((artifact_dir / "metadata.json").read_text())
+        fixture_dir = artifact_dir / metadata["fixtureName"]
+        assert fixture_dir.is_dir(), f"fixture {fixture_dir} missing"
+
+        # The point is the Blosc+Zstd decode path; an uncompressed fixture
+        # would compare nothing interesting.
+        zarray = json.loads((fixture_dir / "0" / ".zarray").read_text())
+        assert zarray["compressor"]["id"] == "blosc"
+        assert zarray["compressor"]["cname"] == "zstd"
+
+        entries = json.loads((artifact_dir / "checksums.json").read_text())["regions"]
+        if isinstance(entries, dict):  # MATLAB jsonencode of a one-element array
+            entries = [entries]
+        assert entries, f"no regions in {source_type} checksums.json"
+
+        mismatches = []
+        for e in entries:
+            region = e["region"]
+            if region is None or len(region) == 0:
+                got = readArray(str(fixture_dir), e["pyramid"], int(e["level"]))
+            else:
+                got = readArray(
+                    str(fixture_dir), e["pyramid"], int(e["level"]), Region=np.array(region)
+                )
+            actual = matlab_sha256(got)
+            if actual != e["sha256"].lower():
+                mismatches.append(f'{e["name"]} ({e["pyramid"]} level {e["level"]})')
+
+        assert not mismatches, (
+            f"SHA-256 mismatch against {source_type} for: {', '.join(mismatches)}. "
+            "The two ports decode the same lossless bytes differently, so one "
+            "of the readers is wrong."
+        )

tests/test_omezarr_reader.py

@@ -188,3 +188,47 @@ def test_maps_common_dtypes(self):
     def test_rejects_unsupported(self):
         with pytest.raises(ValueError):
             ndr_reader_omezarr.zarrClass("<c16")
+
+
+class TestBloscZstdFixture:
+    """The ``codec="blosc-zstd"`` fixture the OME-Zarr symmetry test writes.
+
+    It must really be compressed, with the codec MATLAB's reader decodes, and
+    still read back as the same ground truth as the uncompressed fixture.
+    """
+
+    def test_arrays_are_blosc_zstd(self, tmp_path):
+        import json
+
+        root, _gt = make_example_fixture(tmp_path, with_chunks=True, codec="blosc-zstd")
+        zstd_frames = 0
+        for array in ("0", "mean/1", "mean/2", "max/1", "max/2"):
+            meta = json.loads((root / array / ".zarray").read_text())
+            assert meta["zarr_format"] == 2
+            assert meta["compressor"]["id"] == "blosc"
+            assert meta["compressor"]["cname"] == "zstd"
+            assert meta["order"] == "C" and not meta["filters"]
+            for chunk_file in (root / array).glob("*.*.*.*"):
+                frame = chunk_file.read_bytes()
+                flags = frame[2]
+                assert frame[0] == 2, f"{chunk_file} is not a Blosc v1 frame"
+                # Blosc stores incompressible chunks verbatim (MEMCPYED, 0x02),
+                # and the noise in the full level-0 chunks is incompressible.
+                # Count the frames that really went through Zstd (codec bits
+                # 4 << 5, MEMCPYED clear): without some, the Zstd decode path
+                # this fixture exists to exercise would never run.
+                if flags >> 5 == 4 and not flags & 0x02:
+                    zstd_frames += 1
+        assert zstd_frames > 0, "no chunk is genuinely Zstd-compressed"
+
+    def test_reads_back_the_ground_truth(self, tmp_path):
+        from ndr.format.omezarr.readArray import readArray
+
+        root, gt = make_example_fixture(tmp_path, with_chunks=True, codec="blosc-zstd")
+        np.testing.assert_array_equal(readArray(str(root), "mean", 1), gt.level0)
+        np.testing.assert_array_equal(readArray(str(root), "mean", 2), gt.mean_level1)
+        np.testing.assert_array_equal(readArray(str(root), "max", 3), gt.max_level2)
+
+    def test_rejects_unknown_codec(self, tmp_path):
+        with pytest.raises(ValueError, match="codec"):
+            make_example_fixture(tmp_path, with_chunks=True, codec="lz4")