Merge pull request #23 from milankl/mk/round

faster and more abstract round

Project.toml

@@ -1,7 +1,7 @@
 name = "BitInformation"
 uuid = "de688a37-743e-4ac2-a6f0-bd62414d1aa7"
 authors = ["Milan <[email protected]> and contributors"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

README.md

@@ -4,17 +4,20 @@
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4774191.svg)](https://doi.org/10.5281/zenodo.4774191)
 
 BitInformation.jl is a package for the analysis of bitwise information in Julia Arrays.
-Based on counting the occurrences of bits in floats (or generally any bittype) across various
-dimensions of an array, this package provides functions to calculate quantities like the 
-bitwise real information content, the mutual information, the redundancy or preserved information
-between arrays.
+Based on counting the occurrences of bits in floats (Float16/32/64 or generally any bittype)
+across various dimensions of an array, this package provides functions to calculate quantities
+like the bitwise real information content, the mutual information, the redundancy or preserved
+information between arrays.
 
 BitInformation.jl also implements various rounding modes (round,shave,set_one, etc.)
-efficiently with bitwise operations. Furthermore, transormations like XOR-delta, bittranspose,
+efficiently with bitwise operations. `round(x,i)` implements IEEE's round to nearest tie to even
+for any float retaining `i` mantissa bits. Furthermore, transormations like XOR-delta, bittranspose,
 or signed_exponent are implemented.
 
-If you'd like to propose changes, or contribute in any form raise an issue, create a pull request
-or contact me. Contributions are highly appreciated!
+If you'd like to propose changes, or contribute in any form raise an issue, create a 
+[pull request](https://github.com/milankl/BitInformation.jl/pulls)
+or raise an [issue](https://github.com/milankl/BitInformation.jl/issues).
+Contributions are highly appreciated!
 
 ## Functionality
 
@@ -33,3 +36,9 @@ where `]` opens the package manager. The latest version is automatically install
 
 This project is funded by the [Copernicus Programme](https://www.copernicus.eu/en/copernicus-services/atmosphere) through the [ECMWF summer of weather code 2020 and 2021](https://esowc.ecmwf.int/)
 
+## Reference
+
+If you use this package, please cite the following publication
+
+> M Klöwer, M Razinger, JJ Dominguez, PD Düben and TN Palmer, 2021. *Compressing atmospheric data into its real information content*. **Nature Computational Science** 1, 713–724. [10.1038/s43588-021-00156-2](https://doi.org/10.1038/s43588-021-00156-2)
+

src/BitInformation.jl

@@ -14,7 +14,7 @@ module BitInformation
     include("bitstring.jl")
     include("bittranspose.jl")
     include("shave_set_groom.jl")
-    include("round.jl")
+    include("round_nearest.jl")
     include("xor_delta.jl")
     include("signed_exponent.jl")
     include("bit_information.jl")

src/round_nearest.jl

@@ -0,0 +1,97 @@
+
+"""Shift integer to push the mantissa in the right position. Used to determine
+round up or down in the tie case. `keepbits` is the number of mantissa bits to
+be kept (i.e. not zero-ed) after rounding."""
+function get_shift(::Type{T},keepbits::Integer) where {T<:Base.IEEEFloat}
+    return Base.significand_bits(T) - keepbits
+end
+
+"""Returns for a Float-type `T` and `keepbits`, the number of mantissa bits to be
+kept/non-zeroed after rounding, half of the unit in the last place as unsigned integer.
+Used in round (nearest) to add ulp/2 just before round down to achieve round nearest.
+Technically ulp/2 here is just smaller than ulp/2 which rounds down the ties. For 
+a tie round up +1 is added in `round(T,keepbits)`."""
+function get_ulp_half(  ::Type{T},
+                        keepbits::Integer
+                        ) where {T<:Base.IEEEFloat}
+    # convert to signed for arithmetic bitshift
+    # trick to push in 0s from the left and 1s from the right
+    return ~unsigned(signed(~Base.significand_mask(T)) >> (keepbits+1))
+end
+
+"""Returns a mask that's 1 for all bits that are kept after rounding and 0 for the
+discarded trailing bits. E.g.
+```
+julia> get_keep_mask(Float16,5)
+0xffe0
+```."""
+function get_keep_mask( ::Type{T},
+                        keepbits::Integer
+                        ) where {T<:Base.IEEEFloat}
+    # convert to signed for arithmetic bitshift
+    # trick to push in 1s from the left and 0s from the right
+    return unsigned(signed(~Base.significand_mask(T)) >> keepbits)
+end
+
+"""IEEE's round to nearest tie to even for Float16/32/64."""
+function Base.round(x::T,               # Float to be rounded
+                    ulp_half::UIntT,    # obtain from get_ulp_half,
+                    shift::Integer,     # get_shift, and
+                    keepmask::UIntT     # get_keep_mask
+                    ) where {T<:Base.IEEEFloat,UIntT<:Unsigned}
+    ui = reinterpret(UIntT,x)                       # bitpattern as uint
+    ui += ulp_half + ((ui >> shift) & one(UIntT))   # add ulp/2 with tie to even
+    return reinterpret(T,ui & keepmask)             # set trailing bits to zero
+end
+
+"""Scalar version of `round(::Float,keepbits)` that first obtains
+`shift, ulp_half, keep_mask` and then rounds."""
+function Base.round(x::T,
+                    keepbits::Integer
+                    ) where {T<:Base.IEEEFloat}
+    return round(x,get_ulp_half(T,keepbits),get_shift(T,keepbits),get_keep_mask(T,keepbits))
+end
+
+"""IEEE's round to nearest tie to even for a float array `X` in-place. Calculates from `keepbits`
+only once the variables `ulp_half`, `shift` and `keep_mask` and loops over every element of the
+array."""
+function round!(X::AbstractArray{T},            # any array with element type T
+                keepbits::Integer               # how many mantissa bits to keep
+                ) where {T<:Base.IEEEFloat}     # constrain element type to Float16/32/64
+
+    ulp_half = get_ulp_half(T,keepbits)         # half of unit in last place (ulp)
+    shift = get_shift(T,keepbits)               # integer used for bitshift 
+    keep_mask = get_keep_mask(T,keepbits)       # mask to zero trailing mantissa bits
+
+    @inbounds for i in eachindex(X)             # apply rounding to each element
+        X[i] = round(X[i],ulp_half,shift,keep_mask)
+    end
+
+    return X
+end
+
+"""IEEE's round to nearest tie to even for a float array `X` which returns a rounded copy of `X`."""
+function Base.round(X::AbstractArray{T},        # any array with element type T
+                    keepbits::Integer           # mantissa bits to keep
+                    ) where {T<:Base.IEEEFloat} # constrain element type to Float32/64
+
+    Xcopy = copy(X)                             # copy array to avoid in-place changes
+    round!(Xcopy,keepbits)                      # in-place round the copied array
+    return Xcopy
+end
+
+"""Checks a given `mantissabit` of `x` for eveness. 1=odd, 0=even. Mantissa bits 
+are positive for the mantissa (`mantissabit = 1` is the first mantissa bit), `mantissa = 0`
+is the last exponent bit, and negative for the other exponent bits."""
+function Base.iseven(x::T,
+                    mantissabit::Integer
+                    ) where {T<:Base.IEEEFloat}
+
+    mask = Base.sign_mask(T) >> (Base.exponent_bits(T) + mantissabit)
+    return 0x0 == reinterpret(typeof(mask),x) & mask
+end
+
+"""Checks a given `mantissabit` of `x` for oddness. 1=odd, 0=even. Mantissa bits 
+are positive for the mantissa (`mantissabit = 1` is the first mantissa bit), `mantissa = 0`
+is the last exponent bit, and negative for the other exponent bits."""
+Base.isodd(x::T,mantissabit::Integer) where {T<:Base.IEEEFloat} = ~iseven(x,mantissabit)

src/shave_set_groom.jl

@@ -133,3 +133,17 @@ function groom(X::AbstractArray{Float64},nsb::Integer)
 
     return Y
 end
+
+kouzround(x::Union{Float32,Float64},nsb::Integer) = shave(2x-shave(x,nsb),nsb)
+
+function kouzround(x::AbstractArray{Float32},nsb::Integer)
+    y = similar(x)
+    mask = ~mask32(nsb)
+    for i in eachindex(x)
+        y[i] = shave(2x[i]-shave(x[i],mask),mask)
+    end
+    return y
+end
+
+"""Number of significant bits `nsb` given the number of significant digits `nsd`."""
+nsb(nsd::Integer) = Integer(ceil(log(10)/log(2)*nsd))