Reductions in memory with int64 data type#

In this benchmark we will perform a reduction whose operands will be arrays of type int64. We will compare ironArray’s performance with Zarr, HDF5, Numpy and TileDB. The data used will be the same for each library.

Let’s go:

[1]:
%load_ext memprofiler

import numpy as np
import iarray as ia
import os
import zarr
import dask.array as da
import h5py
import hdf5plugin as h5plugin
import tiledb

In case the reader would like to run this benchmark with a different data type it is enough to change the following variable value:

[2]:
dtype = np.int64

ironArray#

We will create an array with synthetic data with the arange function.

[3]:
%%time

shape = [1000, 100, 100, 100]
stop = 1
for i in range(len(shape)):
    stop *= shape[i]

ia_array_i64 = ia.arange(shape=shape, start=-10, stop=stop - 10, dtype=dtype)

print("cratio: ", round(ia_array_i64.cratio, 2))
chunks = ia_array_i64.chunks
print(ia_array_i64.info)
cratio:  26.3
type   : IArray
shape  : (1000, 100, 100, 100)
chunks : (128, 16, 32, 32)
blocks : (32, 4, 8, 16)
cratio : 26.30

CPU times: user 32.9 s, sys: 1.35 s, total: 34.2 s
Wall time: 20.4 s
[4]:
%%mprof_run 1.iarray::mean_memory
ia_reduc = ia.mean(ia_array_i64, axis=(3, 2, 0)).data
memprofiler: used 38.27 MiB RAM (peak of 38.27 MiB) in 0.8405 s, total RAM usage 832.23 MiB

Zarr#

[5]:
array_zarr = zarr.array(ia_array_i64.data)
print(array_zarr.info)
array_zdask = da.from_zarr(array_zarr)
Type               : zarr.core.Array
Data type          : int64
Shape              : (1000, 100, 100, 100)
Chunk shape        : (125, 13, 13, 25)
Order              : C
Read-only          : False
Compressor         : Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0)
Store type         : builtins.dict
No. bytes          : 8000000000 (7.5G)
No. bytes stored   : 160456467 (153.0M)
Storage ratio      : 49.9
Chunks initialized : 2048/2048

[6]:
%%mprof_run 2.zarr::mean_memory
zdask_reduc = da.mean(array_zdask, axis=(3, 2, 0))
zarr_reduc = zarr.create(shape=ia_reduc.shape, dtype=ia_reduc.dtype)
da.to_zarr(zdask_reduc, zarr_reduc)
memprofiler: used 955.86 MiB RAM (peak of 956.18 MiB) in 1.3548 s, total RAM usage 1951.97 MiB
[7]:
np.testing.assert_array_equal(zarr_reduc, ia_reduc)
del array_zdask
del zdask_reduc
del zarr_reduc

HDF5#

[8]:
h5_urlpath = "mean.hdf5"

if not os.path.exists(h5_urlpath):
    with h5py.File(h5_urlpath, "w") as f:
        h5_array = f.create_dataset("h5_mean", ia_array_i64.shape, dtype=ia_array_i64.dtype, **h5plugin.Blosc())
        ia_array_i64.copyto(h5_array)


h5_file = h5py.File(h5_urlpath, "r", driver='core', backing_store=False)
array_h5 = h5_file['h5_mean']


array_h5dask = da.from_array(array_h5)

h5_reduc_urlpath = "reduc.hdf5"
ia.remove_urlpath(h5_reduc_urlpath)
[9]:
%%mprof_run 3.hdf5::mean_memory

h5dask_reduc = da.mean(array_h5dask, axis=(3, 2, 0))

f = h5py.File(h5_reduc_urlpath, "w", driver='core', backing_store=False)
h5_reduc = f.create_dataset(name=h5_reduc_urlpath, shape=ia_reduc.shape, dtype=ia_reduc.dtype, **h5plugin.Blosc())
da.to_hdf5(h5_reduc_urlpath, '/x', h5dask_reduc)
memprofiler: used 0.26 MiB RAM (peak of 2106.66 MiB) in 2.8027 s, total RAM usage 2258.49 MiB
[10]:
del h5dask_reduc
del array_h5dask
del h5_reduc

TileDB#

[11]:
# TileDB does not come with defaults for chunk shape (tiles), so using ironArray's one
adom = tiledb.Domain(
        tiledb.Dim(name="rows", domain=(0, shape[0] - 1), dtype=ia_array_i64.dtype, tile=chunks[0]),
        tiledb.Dim(name="cols", domain=(0, shape[1] - 1), dtype=ia_array_i64.dtype, tile=chunks[1]),
        tiledb.Dim(name="else", domain=(0, shape[2] - 1), dtype=ia_array_i64.dtype, tile=chunks[2]),
        tiledb.Dim(name="else2", domain=(0, shape[3] - 1), dtype=ia_array_i64.dtype, tile=chunks[3]),
    )

filters = tiledb.FilterList([tiledb.ByteShuffleFilter(), tiledb.LZ4Filter(5)])
aschema = tiledb.ArraySchema(
    domain=adom,  sparse=False, attrs=[tiledb.Attr(name="a", dtype=ia_array_i64.dtype, filters=filters)]
)

# Create the (empty) array on disk.
tiledb_name = "mem://mean.tiledb"
if not os.path.exists(tiledb_name):
    tiledb.DenseArray.create(tiledb_name, aschema)
    with tiledb.DenseArray(tiledb_name, mode="w") as A_opt:
        A_opt[:] = ia_array_i64.data

dask_tiledb = da.from_tiledb(tiledb_name, attribute='a')
[12]:
%%mprof_run 4.tiledb::mean_memory
res_dask_opt = da.mean(dask_tiledb, axis=(3,2,0))
res_dask_opt.to_tiledb("mem://res_tiledb")
memprofiler: used 795.09 MiB RAM (peak of 833.54 MiB) in 11.8555 s, total RAM usage 2803.95 MiB
[13]:
del res_dask_opt
del dask_tiledb
tiledb.remove(tiledb_name)
tiledb.remove("mem://res_tiledb")

NumPy#

Let’s do the same computation with NumPy. First we will get the NumPy array from the ironArray one:

[14]:
%%time

np_array = ia_array_i64.data
CPU times: user 12.8 s, sys: 1.53 s, total: 14.3 s
Wall time: 5.17 s

and now the actual reduction:

[15]:
%%mprof_run 5.numpy::mean_memory
np_reduc = np.mean(np_array, axis=(3, 2, 0))
memprofiler: used 0.01 MiB RAM (peak of 0.01 MiB) in 0.5867 s, total RAM usage 10116.65 MiB
[16]:
np.testing.assert_array_equal(ia_reduc, np_reduc)
del np_array
del np_reduc

Results#

Finally, we will do a plot of the memory and time consumption with each different library.

[17]:
%mprof_plot .*::mean