Bitround in python or julia

[aaronspring]

is xr_bitround the bitround from xarray as discussed here? For performance I'd suggest to use round from BitInformation. I've never actually compared the speed, but recalculation of constants is avoided reaching 16GB/s at which point it's memory-bound (and therefore likely an upper bound). Also gives us more flexibility in case we want to change something on the rounding front.

Originally posted by @milankl in #21 (comment)

[milankl]

For comparison, BitInformation.jl on my MacBook reaches

julia> using BitInformation,BenchmarkTools
julia> A = rand(Float32,10000,10000);    # 400MB array, only fits into RAM not in caches
julia> B = rand(Float64,10000,10000);    # 800MB array
julia> @btime round!($A,7);
  24.991 ms (0 allocations: 0 bytes)

julia> @btime round!($B,7);
  50.709 ms (0 allocations: 0 bytes)

julia> sizeof(A)/1000^2/25e-3
16000.0

julia> sizeof(A)/1000^2/50e-3
16000.0

so about 16GB/s for Float32/Float64 with in-place rounding without any memory allocation and all data has to be pulled from RAM and does not fit into low level caches (which would allow much higher speeds actually). With memory allocations (round! is in-place, round creates a copy of the array and rounds that one in-place) this drops to less than 2GB/s for both Float32 and Float64

julia> @btime round($A,7);
  215.357 ms (2 allocations: 381.47 MiB)

julia> @btime round($B,7);
  458.840 ms (2 allocations: 762.94 MiB)

julia> sizeof(A)/1000^2/215e-3
1860.46511627907

julia> sizeof(B)/1000^2/459e-3
1742.9193899782135

[aaronspring]

how does the julia version handle lazy or out-of-memory data? Imagine I'd want to bitround 1 TB, which is larger than memory often, that should work with xarray and dask

[aaronspring]

can confirm that bitround.round is faster than xr_bitround despite the julia-python wrapping

[milankl]

how does the julia version handle lazy or out-of-memory data? Imagine I'd want to bitround 1 TB, which is larger than memory often, that should work with xarray and dask

At the moment the round function is written for ::AbstractArray in BitInformation.jl which means it'll accept any type that's a subtype of AbstractArray. The function itself does not hardcode how to access elements of such an array, i.e. whether you load everything into ram or keep things lazy isn't prescribed in BitInformation.round. I don't know how the interfaces work between python arrays that support lazy reading and julia, but from the Julia side this should automatically supported. Rounding is also a trivial in parallel, so you can go for any chunking you like.

[milankl]

can confirm that bitround.round is faster than xr_bitround despite the julia-python wrapping

Could you provide some timings, I'm curious.

[milankl]

Also note that the zarr bitround @rabernat's PR zarr-developers/numcodecs#299 only contains float32 at the moment. BitInformation.round supports Float64/32/16 (and any other float type if you define it as <:Base.IEEEFloat and add methods for Base.significand_mask(::MyFloatType) and Base.significand_bits(::MyFloatType))

[aaronspring]

we should definitely use jl_bitround as default

[rabernat]

I plan to add support for the other dtypes, but just have not had time to work on the numcodecs PR recently. @aaronspring - if you want to help out, I would welcome your help on the PR. What is xr_bitround?

[aaronspring]

xr_bitround is just wrapping your PR, which works on arrays, here on xr.DataArrays and xr.Datasets, so nothing really new.
probably out of this project some further developments arise for your PR

[aaronspring]

for chunked data with dask, the julia wrapper doesnt work, see #47