Data array is flattened in numcodecs which reduced the compression ratio that ZFP can provide on multi-dimension arrays

[halehawk]

Minimal, reproducible code sample, a copy-pastable example if possible

In ensure_contiguous_ndarray:
if arr.flags.c_contiguous or arr.flags.f_contiguous:
        # can flatten without copy
        arr = arr.reshape(-1, order='A')

Problem description

ZFP is a lossy compressor that expects to get better compression ratio when compressing multi-dimension arrays. Currently in the zfpy.encode, ensure_contiguous_ndarray function flattens the array and feeds to ZFP compressor. In this way, ZFP compressor only can get less expected compression ratio. So I tried to add a switch in ensure_contiguous_ndarray, if the codec is ZFP, I skipped the flattening. However, this only works on C-order array, it doesn't work on F-order array because ZFP decoder only returns C-order array. When the decoder returns C-order multi-dimension array for an original F-order multi-dimension array, ZFP compressor cannot pass the F-order array encode and decode test.

Version and installation information

Please provide the following:

• 0.9.0 of numcodecs.__version__
• Python 3.7, 3.8
• Operating system (Linux/Windows)
• pip install

[halehawk]

@lindstro @jakirkham please add your comments.

[lindstro]

I believe there are multiple independent issues here.

First, arrays should not be flattened when used with zfp. I see no need to even make arrays contiguous as zfp supports strided compression.

Second, it is true that zfp returns decompressed arrays in C order (the default layout in NumPy), regardless of the original layout. zfp does not capture the original layout, nor is it clear that it should. For instance, if the array being compressed is not contiguous (think array-of-struct storage), would we really want to return a decompressed array full of holes? We did consider supporting strided decompression with zfpy.decompress_numpy(), but as I recall there were some issues with this. That said, I believe one could use zfpy._decompress() to decompress into an existing strided array passed as out argument.

Third, it is an issue if zfp's tests break if given an original array to compress in Fortran order, and I'd be surprised if this were not handled properly. @halehawk, can you provide a simple example that illustrates such a bug? If this is indeed the case, then we need to fix the zfp implementation/tests, not numcodecs.

[rabernat]

Sounds like we could fix the flattening pretty easily with a new keyword argument here flatten=True

numcodecs/numcodecs/compat.py

Line 53 in a81856e

def ensure_contiguous_ndarray(buf, max_buffer_size=None):

All the existing codecs would work the same with the default flatten=True, but we could call it from zfp here

numcodecs/numcodecs/zfpy.py

Line 58 in a81856e

buf = ensure_contiguous_ndarray(buf)

with flatten=False

[halehawk]

But using Ryan's way, it cannot pass the test at test_zfpy.py if the array is https://github.com/zarr-developers/numcodecs/blob/a81856efbf8dfd8ea925acaac01235a096d285e4/numcodecs/tests/test_zfpy.py#L41 np.asfortranarray(np.random.normal(loc=1000, scale=1, size=(5, 10, 20))), because zfp decoder returns a C-order multidimensional array.

…

On Tue, Feb 8, 2022 at 3:59 PM Ryan Abernathey ***@***.***> wrote: Sounds like we could fix the flattening pretty easily with a new keyword argument here flatten=True https://github.com/zarr-developers/numcodecs/blob/a81856efbf8dfd8ea925acaac01235a096d285e4/numcodecs/compat.py#L53 All the existing codecs would work the same with the default flatten=True, but we could call it from zfp here https://github.com/zarr-developers/numcodecs/blob/a81856efbf8dfd8ea925acaac01235a096d285e4/numcodecs/zfpy.py#L58 with flatten=False — Reply to this email directly, view it on GitHub <#303 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACAPEFDTBEMOQUEQCUVP5K3U2GN47ANCNFSM5N3H636A> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[rabernat]

So why don't we just disallow F order for the zfp codec? I am becoming more and more convinced that it is unnecessary to support it. See discussion in zarr-developers/community#41.

[halehawk]

@rabernat I agree with you it is great if we can skip Fortran array test here. I tested xarray.open_dataset with backend netcdf4, it always opens a netcdf file as a c-order data array even the netcdf file was written to f-order originally. I didn't check with hdf5 or scipy backends. If anyone can confirm that hdf5 and scipy are the same cases, then I think it is all right to skip the Fortran array test in the test_zfpy.py.

[rabernat]

HDF5 and NetCDF simply do not support F order on disk at all. I think it's fine for a particular codec to not support F-order memory layout. The user can always transform their data if need be if they happen to be creating F order. We just need to raise an informative error in the event of seeing F order inputs to the codec.

[lindstro]

HDF5 and NetCDF simply do not support F order on disk at all.

I partially agree. HDF5 imposes C order and (to my knowledge) does not support strided storage, where strides are recorded in the file. I believe it does not support strides on reads either, e.g., to force F order. But you can always store F order data in HDF5 by simply reversing the order of dimensions.

zfpy actually examines the NumPy array's strides and passes them on to zfp. Because zfp uses F order, zfpy reverses the dimensions in anticipation that in most cases, the NumPy order will be the default C order, which promotes good data locality during compression.

Now, the original memory layout is not recorded by zfp, just like how this information is absent from HDF5. Therefore, if upon decompression the user wants to recover the array in the same order as it was compressed in, that kind of information would have to be stored separately and the data would have to be permuted.

That said, it is not clear to me why it's necessary to ensure that compressed and decompressed data appear in the same order. Both arrays will agree on the dimensions of the array, just not how it is organized in memory. Is there a compelling reason to guarantee that both arrays have the same layout? After all, numcodecs forces both contiguous and flat storage, so the original memory layout is already lost, no?

[jakirkham]

Ryan's suggestion seems reasonable.

Forget whether we discussed this, but does zfp work with N-D data? If not, what happens if it is fed an array with N greater than what is supported by zfp.

Also how would this work for decompression? Does zfp return an array of the right dimensionality?

[lindstro]

Only when 1 <= N <= 4. :-) For higher-dimensional arrays, you need to consolidate some dimensions (ideally the ones that are least correlated). This could for instance be done using a reshape operation, though you'd probably want to keep track of the dimensions of the original array.

We're planning to support higher-dimensional arrays in future versions of zfp. It is unlikely that the compression scheme would support more than four dimensions (it's unusual to have data that is correlated in that many dimensions), but the API would allow you to describe higher-dimensional arrays and to identify the smoothest dimensions along which to compress (may be less than four).

The current zfp C API does not accommodate N > 4; the zfpy module will raise an exception if asked to compress a NumPy array with N > 4.

[jakirkham]

Thanks for the info Peter :)

Follow up question, when the data is decompressed, will it return the expected N-D array or is it 1-D (IOW up to the user to reshape)?

Side note: In microscopy it is common to have 5-D data (TZYXC or time + 3D spatial + channel).

[lindstro]

Thanks for the info Peter :)

Follow up question, when the data is decompressed, will it return the expected N-D array or is it 1-D (IOW up to the user to reshape)?

zfpy will return a decompressed N-D NumPy array of the same shape as the array that was compressed. However, the way zfpy currently interfaces with numcodecs, it is handed a 1D reshaped array to compress, which is the issue at hand. In this case, zfpy will also return a 1D decompressed array. Moreover, compression is severely handicapped by not being able to exploit higher-dimensional correlations.

Side note: In microscopy it is common to have 5-D data (TZYXC or time + 3D spatial + channel).

Right, though I suspect there's little correlation across channels. In radiation transport, we deal with 6D data (3 spatial + 2 angular + energy group) that is correlated in all six dimensions, but there are diminishing returns in compressing more than three dimensions (at least as far as zfp is concerned). I agree that zfp should handle such data sets, but probably in the API layer and not in the compression algorithm itself.

[halehawk]

Peter, I saw your zfpy pip package has Python 3.9 and 3.10 version, does your 3.9 macos version work on mac? Did you merge my PR for building that pip package?

…

On Thu, Feb 10, 2022 at 2:14 PM Peter Lindstrom ***@***.***> wrote: Thanks for the info Peter :) Follow up question, when the data is decompressed, will it return the expected N-D array or is it 1-D (IOW up to the user to reshape)? zfpy will return a decompressed N-D NumPy array of the same shape as the array that was compressed. However, the way zfpy currently interfaces with numcodecs, it is handed a 1D reshaped array to compress, which is the issue at hand. In this case, zfpy will also return a 1D decompressed array. Moreover, compression is severely handicapped by not being able to exploit higher-dimensional correlations. Side note: In microscopy it is common to have 5-D data (TZYXC or time + 3D spatial + channel). Right, though I suspect there's little correlation across channels. In radiation transport, we deal with 6D data (3 spatial + 2 angular + energy group) that is correlated in all six dimensions, but there are diminishing returns in compressing more than three dimensions (at least as far as zfp is concerned). I agree that zfp should handle such data sets, but probably in the API layer and not in the compression algorithm itself. — Reply to this email directly, view it on GitHub <#303 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACAPEFFUNFYA5JQ4XBUIILLU2QTFBANCNFSM5N3H636A> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you modified the open/close state.Message ID: ***@***.***>