julia icon indicating copy to clipboard operation
julia copied to clipboard

maximum(abs, ::Array{Complex{Int}}; dims=...) fails with an internal error

Open eschnett opened this issue 1 year ago • 4 comments

This code fails:

julia> A = Complex{Int}.(rand(Complex{Int8}, 2, 2, 2));
julia> maximum(abs, A; dims=(1, 2))

The reason is that Julia wants to store the result of abs(::Complex{T}) as type T. This fails for integer T since abs returns a real number (not an integer) when applied to a complex number.

The error is probably in the line

_realtype(::Union{typeof(abs),typeof(abs2)}, T) = _realtype(T)

in reducedim.jl. This lines doesn't take into account that abs may return a type that's different from its input.

The error is

ERROR: InexactError: Int64(137.0109484676316)
Stacktrace:
  [1] Int64
    @ ./float.jl:912 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./array.jl:1024 [inlined]
  [4] setindex!
    @ ./multidimensional.jl:698 [inlined]
  [5] _mapreducedim!(f::typeof(abs), op::typeof(max), R::Array{Int64, 3}, A::Array{Complex{Int64}, 3})
    @ Base ./reducedim.jl:311
  [6] mapreducedim!(f::Function, op::Function, R::Array{Int64, 3}, A::Array{Complex{Int64}, 3})
    @ Base ./reducedim.jl:324
  [7] _mapreduce_dim(f::Function, op::Function, ::Base._InitialValue, A::Array{Complex{Int64}, 3}, dims::Tuple{Int64, Int64})
    @ Base ./reducedim.jl:371
  [8] mapreduce
    @ ./reducedim.jl:357 [inlined]
  [9] _maximum
    @ ./reducedim.jl:1039 [inlined]
 [10] #maximum#837
    @ ./reducedim.jl:1011 [inlined]
 [11] top-level scope
    @ REPL[64]:1

As a work-around one can write maximum(abs ∘ identity, ...).

I am using Julia 1.10.4.

eschnett avatar Jun 24 '24 19:06 eschnett

Sounds very similar to https://github.com/JuliaLang/julia/pull/53461 (which was accumulate and cumsum, misusing reduce_first), for someone that wants to tackle this

vtjnash avatar Jun 24 '24 20:06 vtjnash

Sounds more like we need

_realtype(::typeof(abs2), T) = _realtype(T)
_realtype(::typeof(abs), T) = float(_realtype(T))

?

dkarrasch avatar Jun 24 '24 21:06 dkarrasch

Not quite, abs(::T<:Real) returns T. We need to distinguish between abs(::T) and abs(::Complex{T}).

Maybe this?

_realtype(::typeof(abs2), T) = _realtype(T)
_realtype(::typeof(abs), T) = T
_realtype(::typeof(abs), Complex{T}) = float(T)

eschnett avatar Jun 24 '24 22:06 eschnett

Indeed, I was confusing abs and norm. norm(1) returns 1.0.

dkarrasch avatar Jun 25 '24 07:06 dkarrasch