InexactError when the intermediate reduction type differs from the final type

[disberd]

I came into an InexactError error when using Tullio to perform reduction and then abs2 over complex vectors.

The MWE might be easier to explain the issue:

using FFTW, Tullio
x = rand(120) .- .5
y = rand(120) .- .5
u = rand(500) .- .5
v = rand(500) .- .5
a = rand(120)
k = 2π/.01

function tfunc(u,v,a,k,x,y)
    @tullio vals[i] := exp(im * (k * (x[j]*u[i] + y[j]*v[i]) + a[j])) |> abs2
end

tfunc(u,v,a,k,x,y)

Which results in

ERROR: LoadError: InexactError: Float64(-13.863000778992978 + 2.707549484306247im)

In my example the issue could easily be resolved by performing the abs2 directly on the vals array returned by the @tullio macro. I briefly discussed this with @mcabbott on slack and this appears to only happen when threading is enabled in @tullio.

[mcabbott]

Thanks for the note. This might be tricky to fix, will think a bit.

The other type inference issue which needs a look is this -- ideally it would realise that bool + bool gives an Int:

b = randn(2,10) .> 0
@tullio s[i] := b[i,j]  # InexactError: Bool(4)

[mcabbott]

This came up in https://discourse.julialang.org/t/fast-4d-argmax/58566

One possibility would be that, when the macro sees a finaliser (such as abs2), it always removes the reduction indices from consideration by the tiling code. This would allow threading on other indices, and should be safe.

Another would be that, during the setup, it could test that the type of the intermediate reduction matches the final type, and throw an error. (If threading is enabled.) This would at least avoid the confusing situation where the above case works for small arrays, but fails for large ones:

julia> x = rand(5) .- .5; y = copy(x); u = copy(x); v = copy(x); a = copy(x);

julia> tfunc(u,v,a,k,x,y)
5-element Vector{Float64}:
  0.9609872326367754
  6.200966967475577
  2.8161379721911195
 12.243066458238065
  7.800953389121679