Fiddling with sumlog
#1
Conversation
| @inline function _sumlog_op((sig1, ex1), (sig2, ex2)) | ||
| sig = sig1 * sig2 | ||
| # sig = ifelse(sig2<0, sig2, sig1 * sig2) |
There was a problem hiding this comment.
Here I was wondering whether, instead of xj < 0 && throw(... at every step, you could just ensure that sig ends up negative so that the final result throws. But it was slower.
| @inline _sumlog(::Type{T}, x) where {T} = sum(log, x) | ||
| # The exported `exponent(x)` checks for `NaN` etc, this function doesn't, which is fine as `sig` keeps track. | ||
| _exponent(x::Base.IEEEFloat) = Base.Math._exponent_finite_nonzero(x) | ||
| Base.@assume_effects :nothrow _exponent(x::AbstractFloat) = Int(exponent(x)) # e.g. for BigFloat |
There was a problem hiding this comment.
This Base.@assume_effects I got from a comment in Base, but it doesn't do what I wanted here.
exponent(::BigFloat) is a C call, with the error check we'd like to skip in the same function.
| if isnothing(iter) | ||
| T = Base._return_type(first, Tuple{typeof(x)}) | ||
| return T <: Number ? zero(float(T)) : 0.0 |
There was a problem hiding this comment.
For the empty iterator case I think you can't avoid relying on inference.
There was a problem hiding this comment.
Does this still infer ok? In another library (maybe it was MappedArrays?) I've seen this written something like
infer() = Base._return_type(first, Tuple{typeof(x)})
T = infer()There was a problem hiding this comment.
I have never seen this, I don't know why it would be different, but maybe there is some black magic.
There was a problem hiding this comment.
Found it:
JuliaArrays/MappedArrays.jl@0e1b57f
| float_xj = float(xj) | ||
| significand(float_xj), _exponent(float_xj) | ||
| end | ||
| return log(sig) + IrrationalConstants.logtwo * T(ex) |
There was a problem hiding this comment.
Previously, logtwo * ex was giving Float64 always.
| float_xj = float(xj) | ||
| significand(float_xj), _exponent(float_xj) | ||
| end | ||
| return log(sig) + IrrationalConstants.logtwo * T(ex) |
| if isnothing(iter) | ||
| T = Base._return_type(first, Tuple{typeof(x)}) | ||
| return T <: Number ? zero(float(T)) : 0.0 |
There was a problem hiding this comment.
Does this still infer ok? In another library (maybe it was MappedArrays?) I've seen this written something like
infer() = Base._return_type(first, Tuple{typeof(x)})
T = infer()| @test sumlog([1,2,3]) isa Float64 | ||
| @test sumlog([1,2,3]) ≈ sum(log, [1,2,3]) | ||
| @test sumlog([1 2; 3 4]; dims=1) ≈ sum(log, [1 2; 3 4]; dims=1) | ||
| @test sumlog(Int(x) for x in Float64[1,2,3]) ≈ sum(log, [1,2,3]) | ||
| end | ||
| end |
|
|
||
| Does not accept a `dims` keyword. | ||
| """ | ||
| sumlog(f, x) = sumlog(Iterators.map(f, x)) |
There was a problem hiding this comment.
Should we include a sumlog(f, x; dims=:) version for arrays? Seems like an easy update to include f inside the function and have it default to identity
There was a problem hiding this comment.
I wondered a little. But the way that's currently set up, it uses float(eltype(x)) to know whether it can run this algorithm. If f isn't type-stable, then perhaps things get awkward... surely it can be done but not 1 line.
Here if you want it is my attempt to make
sumexpwork with dims, and on arbitrary iterators. Surely not free of mistakes yet, but...Looks like I touched every line, just about, in fiddling with this.
Float16 doesn't work so well yet, tests fail. Maybe it ought to accumulate in Float32 or something.
BigFloat has issues with NaN etc, tests skipped.