DSP.jl
DSP.jl copied to clipboard
xcorr errors with integer arguments
julia> using DSP
julia> xcorr(rand(Int, 5), rand(Int, 5))
ERROR: InexactError: trunc(Int64, 7.28283193350032e36)
Stacktrace:
[1] trunc at ./float.jl:693 [inlined]
[2] round at ./float.jl:359 [inlined]
[3] _broadcast_getindex_evalf at ./broadcast.jl:574 [inlined]
[4] _broadcast_getindex at ./broadcast.jl:557 [inlined]
[5] getindex at ./broadcast.jl:507 [inlined]
[6] macro expansion at ./broadcast.jl:838 [inlined]
[7] macro expansion at ./simdloop.jl:73 [inlined]
[8] copyto! at ./broadcast.jl:837 [inlined]
[9] copyto! at ./broadcast.jl:792 [inlined]
[10] copy at ./broadcast.jl:768 [inlined]
[11] materialize at ./broadcast.jl:748 [inlined]
[12] conv(::Array{Int64,1}, ::Array{Int64,1}) at /home/glynch/.julia/packages/DSP/8zBXI/src/dspbase.jl:141
[13] xcorr(::Array{Int64,1}, ::Array{Int64,1}) at /home/glynch/.julia/packages/DSP/8zBXI/src/dspbase.jl:201
[14] top-level scope at none:0
I managed to narrow down this issue. When cross-correlation or convolution (xcorr or conv) intake integer vectors with values larger than half the bit size of the system (i.e. larger than typemax(Int32) on a system where Int defaults to Int64) the calculation can create final float values larger than typemax(Int64). These floats cannot be rounded and an error is thrown. Typical users probably won't run into this as it's not likely you would be using these functions on extremely large values.
Nonetheless, there are 2 solutions I can see.
- Don't convert the output of conv or xcorr to Integers. Do users expect integers back from a convolution when they supply integer vectors?
- Limit the integer input to values that are the typemax of the size below system size (i.e. Int32 on 64 bit systems and maybe Int16 on 32 bit systems?)
Either would be relatively easy to implement. Just need to know which is preferred and I can create a PR.
Using floats would work in more scenarios, though you are right that most users will probably not run into this issue. Limiting the inputs would have the benefit of not changing the interface, but it wouldn't be guaranteed to prevent the error:
julia> a = fill(typemax(Int32), 5); conv(a, a)
ERROR: InexactError: trunc(Int64, 1.3835058042397262e19)
I'm a new maintainer so I don't know what the policy on breaking changes is, but I think not converting the floating point result to integers would make the most sense.
What do other people think?
I tend to prefer keeping integers if supplied. We could, however, consider widening the type similar to what sum
does, e.g.:
julia> a = fill(typemax(Int32), 5); sum(a)
10737418235
Relatedly, we could provide conv!(dest, a, b)
which would allow the caller to implicitly choose the output type by providing a suitable dest
array.
Changing the Integer handling of conv to conv(u::StridedVector{T}, v::StridedVector{T}) where {T<:Integer} = round.(Int128, conv(float(u), float(v)))
from conv(u::StridedVector{T}, v::StridedVector{T}) where {T<:Integer} = round.(Int, conv(float(u), float(v)))
Atleast solves the problem for input of typemax(Int32) but just pushes the problem off to a larger type. Now a = fill(typemax(Int64),5); conv(a,a)
gives the same error. I'm not sure how to solve that problem without resorting to BigInt. Or we could limit input to Int32 and return Int128, preventing any errors (I think).