Preserve FixedScaleOffset array dimensions

docs/release.md

@@ -14,6 +14,12 @@
 
 ## Unreleased
 
+### Bug fixes
+
+* Preserve multidimensional array shapes in `FixedScaleOffset` and preserve logical
+  coordinates when copying equally shaped arrays across different memory orders.
+  By {user}`shixi-li <shixi-li>`, {issue}`852`
+
 (release_0.17.0)=
 
 * Allow building against a system Zlib. Note, this is only required if blosc is still

src/numcodecs/compat.py

@@ -185,13 +185,21 @@ def ndarray_copy(src, dst) -> NDArrayLike:
     src = ensure_ndarray_like(src)
     dst = ensure_ndarray_like(dst)
 
-    # flatten source array
-    src = src.reshape(-1, order="A")
-
     # ensure same data type
     if dst.dtype != object:
+        # a different itemsize requires a contiguous last axis for the dtype view
+        if src.dtype.itemsize != dst.dtype.itemsize:
+            src = src.reshape(-1, order="A")
         src = src.view(dst.dtype)
 
+    # preserve logical coordinates when equally shaped arrays use different memory orders
+    if src.shape == dst.shape:
+        np.copyto(dst, src)
+        return dst
+
+    # flatten source array
+    src = src.reshape(-1, order="A")
+
     # reshape source to match destination
     if src.shape != dst.shape:
         if dst.flags.f_contiguous:

src/numcodecs/fixedscaleoffset.py

@@ -84,30 +84,34 @@ def encode(self, buf):
         # normalise input
         arr = ensure_ndarray(buf).view(self.dtype)
 
-        # flatten to simplify implementation
-        arr = arr.reshape(-1, order='A')
+        # preserve the historical one-element shape for scalar inputs
+        if arr.ndim == 0:
+            arr = arr.reshape(-1)
 
         # compute scale offset
         enc = (arr - self.offset) * self.scale
 
         # round to nearest integer
         enc = np.around(enc)
 
-        # convert dtype
-        return enc.astype(self.astype, copy=False)
+        # preserve the legacy byte order, including for non-contiguous inputs
+        order = "F" if arr.flags.f_contiguous else "C"
+        return enc.astype(self.astype, order=order, copy=False)
 
     def decode(self, buf, out=None):
         # interpret buffer as numpy array
         enc = ensure_ndarray(buf).view(self.astype)
 
-        # flatten to simplify implementation
-        enc = enc.reshape(-1, order='A')
+        # preserve the historical one-element shape for scalar inputs
+        if enc.ndim == 0:
+            enc = enc.reshape(-1)
 
         # decode scale offset
         dec = (enc / self.scale) + self.offset
 
-        # convert dtype
-        dec = dec.astype(self.dtype, copy=False)
+        # preserve the legacy byte order when copying into a differently shaped output
+        order = "F" if enc.flags.f_contiguous else "C"
+        dec = dec.astype(self.dtype, order=order, copy=False)
 
         # handle output
         return ndarray_copy(dec, out)

tests/test_compat.py

@@ -4,7 +4,7 @@
 import numpy as np
 import pytest
 
-from numcodecs.compat import ensure_bytes, ensure_contiguous_ndarray, ensure_text
+from numcodecs.compat import ensure_bytes, ensure_contiguous_ndarray, ensure_text, ndarray_copy
 
 
 def test_ensure_text():
@@ -109,3 +109,43 @@ def test_ensure_contiguous_ndarray_max_buffer_size():
         for buf in buffers:
             with pytest.raises(ValueError):
                 ensure_contiguous_ndarray(buf, max_buffer_size=max_buffer_size)
+
+
+@pytest.mark.parametrize(("source_order", "destination_order"), [("C", "F"), ("F", "C")])
+@pytest.mark.parametrize("destination_type", ["ndarray", "memoryview"])
+def test_ndarray_copy_same_shape_preserves_logical_coordinates(
+    source_order, destination_order, destination_type
+):
+    shape = (2, 3, 4)
+    source = np.arange(np.prod(shape), dtype="<f4").reshape(shape, order=source_order)
+    destination_array = np.empty(shape, dtype="<f4", order=destination_order)
+    destination = (
+        memoryview(destination_array) if destination_type == "memoryview" else destination_array
+    )
+
+    ndarray_copy(source, destination)
+
+    np.testing.assert_array_equal(destination_array, source)
+
+
+@pytest.mark.parametrize("source_layout", ["F", "sliced"])
+@pytest.mark.parametrize("destination_order", ["C", "F"])
+@pytest.mark.parametrize("destination_type", ["ndarray", "memoryview"])
+def test_ndarray_copy_different_itemsize_preserves_bytes(
+    source_layout, destination_order, destination_type
+):
+    if source_layout == "F":
+        source = np.arange(12, dtype="<u2").reshape(3, 4, order="F")
+    else:
+        source = np.arange(48, dtype="<u2").reshape(6, 8)[::2, ::2]
+    destination_array = np.empty((3, 8), dtype="u1", order=destination_order)
+    destination = (
+        memoryview(destination_array) if destination_type == "memoryview" else destination_array
+    )
+
+    result = ndarray_copy(source, destination)
+
+    assert result.shape == destination_array.shape
+    assert result.dtype == destination_array.dtype
+    assert np.shares_memory(result, destination_array)
+    assert destination_array.tobytes(order="A") == source.tobytes(order="A")

tests/test_fixedscaleoffset.py

@@ -1,8 +1,9 @@
 import itertools
+from typing import Literal
 
 import numpy as np
 import pytest
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_array_almost_equal, assert_array_equal
 
 from numcodecs.fixedscaleoffset import FixedScaleOffset
 from tests.common import (
@@ -55,9 +56,181 @@ def test_encode(offset: float, scale: float, expected: list[int]):
     assert np.dtype(astype) == actual.dtype
 
 
+@pytest.mark.parametrize("order", ["C", "F"])
+@pytest.mark.parametrize("input_type", ["ndarray", "memoryview"])
+def test_preserves_dimensions(order: Literal["C", "F"], input_type):
+    shape = (8, 4, 5)
+    arr = 200.0 + np.arange(np.prod(shape), dtype="<f4").reshape(shape, order=order) * 0.001
+    codec = FixedScaleOffset(
+        offset=200.0,
+        scale=1000.0,
+        dtype="<f4",
+        astype="<u2",
+    )
+
+    encoded = codec.encode(memoryview(arr) if input_type == "memoryview" else arr)
+    decoded = codec.decode(memoryview(encoded) if input_type == "memoryview" else encoded)
+
+    assert encoded.shape == shape
+    assert encoded.dtype == codec.astype
+    assert decoded.shape == shape
+    assert decoded.dtype == arr.dtype
+    assert_array_almost_equal(arr, decoded, decimal=3)
+
+    serialized = encoded.tobytes(order="A")
+    decoded_serialized = codec.decode(serialized)
+    assert decoded_serialized.shape == (arr.size,)
+    assert_array_almost_equal(arr.reshape(-1, order="A"), decoded_serialized, decimal=3)
+
+    out = np.empty_like(arr)
+    decoded_out = codec.decode(serialized, out=out)
+    assert decoded_out is out
+    assert_array_almost_equal(arr, decoded_out, decimal=3)
+
+    legacy_encoded = np.around((arr.reshape(-1, order="A") - codec.offset) * codec.scale).astype(
+        codec.astype,
+        copy=False,
+    )
+    assert encoded.tobytes(order="A") == legacy_encoded.tobytes(order="A")
+
+
+@pytest.mark.parametrize(("input_order", "out_order"), [("C", "F"), ("F", "C")])
+def test_decode_out_preserves_logical_coordinates(
+    input_order: Literal["C", "F"], out_order: Literal["C", "F"]
+):
+    shape = (2, 3, 4)
+    arr = np.arange(np.prod(shape), dtype="<f4").reshape(shape, order=input_order)
+    codec = FixedScaleOffset(
+        offset=0,
+        scale=1,
+        dtype="<f4",
+        astype="<i4",
+    )
+    encoded = codec.encode(arr)
+    out = np.empty(shape, dtype="<f4", order=out_order)
+
+    decoded = codec.decode(encoded, out=out)
+
+    assert decoded is out
+    assert decoded.shape == shape
+    assert decoded.dtype == arr.dtype
+    assert_array_equal(decoded, arr)
+
+
+def test_scalar_preserves_legacy_shape():
+    arr = np.array(3.5, dtype="<f4")
+    codec = FixedScaleOffset(offset=0, scale=10, dtype="<f4")
+
+    encoded = codec.encode(arr)
+    decoded = codec.decode(encoded)
+
+    assert encoded.shape == (1,)
+    assert decoded.shape == (1,)
+    assert_array_equal(decoded, arr.reshape(-1))
+
+
+@pytest.mark.parametrize("out_shape", [None, (), (1,), (1, 1)])
+def test_decode_scalar_preserves_legacy_shape(out_shape):
+    codec = FixedScaleOffset(offset=1, scale=10, dtype="<f4", astype="<i2")
+    encoded = np.array(25, dtype=codec.astype)
+    out = None if out_shape is None else np.empty(out_shape, dtype=codec.dtype)
+
+    decoded = codec.decode(encoded, out=out)
+
+    expected_shape = (1,) if out is None else out.shape
+    assert decoded.shape == expected_shape
+    assert decoded.dtype == codec.dtype
+    assert_array_equal(decoded, np.full(expected_shape, 3.5, dtype=codec.dtype))
+    if out is not None:
+        assert decoded is out
+
+
+@pytest.mark.parametrize("shape", [(0,), (0, 3), (2, 0, 3)])
+@pytest.mark.parametrize("order", ["C", "F"])
+def test_empty_array_preserves_dimensions(shape, order):
+    codec = FixedScaleOffset(offset=0, scale=10, dtype="<f4", astype="<i2")
+    arr = np.empty(shape, dtype=codec.dtype, order=order)
+
+    encoded = codec.encode(arr)
+    decoded = codec.decode(encoded)
+    out = np.empty_like(arr)
+
+    assert encoded.shape == shape
+    assert encoded.dtype == codec.astype
+    assert decoded.shape == shape
+    assert decoded.dtype == arr.dtype
+    assert_array_equal(decoded, arr)
+    assert codec.decode(encoded, out=out) is out
+    assert_array_equal(out, arr)
+    assert codec.decode(encoded.tobytes()).shape == (0,)
+
+
+@pytest.mark.parametrize(("input_order", "out_order"), [("C", "F"), ("F", "C")])
+def test_decode_different_output_shape(input_order, out_order):
+    codec = FixedScaleOffset(offset=0, scale=1, dtype="<f4", astype="<i2")
+    arr = np.arange(24, dtype=codec.dtype).reshape(2, 3, 4, order=input_order)
+    out = np.empty((4, 6), dtype=codec.dtype, order=out_order)
+
+    decoded = codec.decode(codec.encode(arr), out=out)
+
+    assert decoded is out
+    assert_array_equal(decoded, np.arange(24, dtype=codec.dtype).reshape(4, 6, order=out_order))
+
+
+@pytest.mark.parametrize("order", ["C", "F"])
+def test_decode_memoryview_to_bytearray(order):
+    codec = FixedScaleOffset(offset=0, scale=1, dtype="<f4", astype="<i2")
+    arr = np.arange(6, dtype=codec.dtype).reshape(2, 3, order=order)
+    encoded = codec.encode(arr)
+    out = bytearray(arr.nbytes)
+
+    decoded = codec.decode(memoryview(encoded), out=out)
+
+    assert decoded.shape == (arr.nbytes,)
+    assert decoded.dtype == np.dtype("u1")
+    assert np.shares_memory(decoded, np.asarray(memoryview(out)))
+    assert bytes(out) == arr.tobytes(order="A")
+
+
+@pytest.mark.parametrize("order", ["C", "F"])
+def test_noncontiguous_input_preserves_legacy_bytes(order):
+    codec = FixedScaleOffset(offset=0, scale=10, dtype="<f4", astype="<i2")
+    arr = np.arange(24, dtype=codec.dtype).reshape(4, 6, order=order)[::2, :]
+    encoded = codec.encode(arr)
+    legacy_encoded = np.around((arr.reshape(-1, order="A") - codec.offset) * codec.scale).astype(
+        codec.astype
+    )
+
+    assert encoded.shape == arr.shape
+    assert encoded.dtype == codec.astype
+    assert encoded.tobytes(order="A") == legacy_encoded.tobytes()
+    assert_array_equal(codec.decode(encoded), arr)
+    assert_array_equal(codec.decode(encoded.tobytes(order="A")), arr.reshape(-1, order="A"))
+    out = np.empty(arr.shape, dtype=codec.dtype)
+    assert codec.decode(encoded.tobytes(order="A"), out=out) is out
+    assert_array_equal(out, arr)
+
+
+@pytest.mark.parametrize("order", ["C", "F"])
+def test_decode_noncontiguous_input_preserves_legacy_order(order):
+    codec = FixedScaleOffset(offset=1, scale=2, dtype="<f4", astype="<i2")
+    encoded = np.arange(24, dtype=codec.astype).reshape(4, 6, order=order)[::2, :]
+    expected = (encoded / codec.scale + codec.offset).astype(codec.dtype)
+
+    decoded = codec.decode(encoded)
+
+    assert decoded.shape == encoded.shape
+    assert decoded.dtype == codec.dtype
+    assert_array_equal(decoded, expected)
+    out = np.empty(encoded.size, dtype=codec.dtype)
+    assert codec.decode(memoryview(encoded), out=out) is out
+    assert_array_equal(out, expected.reshape(-1, order="C"))
+
+
 def test_config():
     codec = FixedScaleOffset(dtype='<f8', astype='<i4', scale=10, offset=100)
     check_config(codec)
+    assert "keep_dimensions" not in codec.get_config()
 
 
 def test_repr():