Unexpected Tullio behavior

General Usage Performance
,
1

In the program below, I had expected a, a2 to produce the same result, but they don’t. It appears that Tullio (which I like!) does not deal with vectors of matrices in the way that I had expected.

The output is:

[ Info: loop 1
true
[ Info: loop 2
false

Two questions:

  1. What is the best way of running a loop like the one indicated?
  2. What other gotchas should I be aware of?

using Tullio, LoopVectorization


function test()
    J = [2; 3];  I = 4; R= 5; Jmax = maximum(J)
    x = Vector{ Matrix{Float64} }(undef, 2);
    a = fill(typemax(Float64), I, Jmax);  a2 = fill(typemax(Float64),I, Jmax);  
    for m = 1:2
        x[m] = randn(I,J[m])
    end
    for m = 1:2
        @info "loop $m"
        @tullio a[i,j] = x[$m][i,j+0]              (j in 1:J[m])
        xm = x[m]
        @tullio a2[i,j] = xm[i,j+0]              (j in 1:J[m])
        println( a[:,1:J[m]] ≈ a2[:,1:J[m]] )
#~         println(a);  println(a2)
    end
end

test()
2

I am by no means an expert on Tullio but like you am a satisfied user.

I think this particular use-case is quite strange for Tullio. You have only one loop, which is just summing numbers. The only “difficult” thing about it is the indexing. I’d suggest just writing a small function which does this for you, except maybe it already exists (?):

a[i,j] = sum(x[m],dims=2)

I think Tullio is best for heavy matrix algebra where there are several loops to consider.

3

Thanks @tbeason. My actual program is much more complicated and demanding than this. I’ve just tried to create a small working example that demonstrates the behavior that confuses me.

4

I think your indexing is wrong somehow. Especially on a2. If you are just summing columns, a2 should just be a vector, ie a2[i] and not a2[i,j]

5

No, I’m not summing columns. This is just a simple example to demonstrate the issue, but I can create an example involving sums that produces the same problem.

6

This is a gotcha I didn’t think of. Here’s what’s happening, first: the expected range for j is 1:3, but in fact it runs 1:2, so doesn’t fill the array:

julia> m = 2;

julia> axes(x[m])
(Base.OneTo(4), Base.OneTo(3))

julia> J[m]
3

julia> @tullio a[i,j] = x[$m][i,j+0]   (j in 1:J[m])  verbose=true
┌ Warning: LoopVectorization failed 
...
┌ Info: left index ranges
│   i = Base.OneTo(4)
└   j = 1:2
4×3 Array{Float64,2}:
 -0.696945  1.8923    Inf
  0.199954  0.645115  Inf
  1.30176   0.318919  Inf
 -0.669092  0.93308   Inf

These ranges are calculated like so (printed by verbose=2), and axes(first(x), 2) assumes that any array of arrays has all slices the same size:

│        local 𝒶𝓍j = intersect(axes(a, 2), axes(first(x), 2) .- 0, var"≪1:J[m]≫")
│        local 𝒶𝓍i = axes(a, 1)
│        axes(first(x), 1) == axes(a, 1) || throw("range of index i must agree")

Changing that would be pretty hard, as the range of k in some A[i][k] would have to be worked out inside the loop over i, and it really doesn’t think like that. The ranges of all indices are fixed, and then a kernel function is called on (small parts of) the rectangular iteration space. (It should really check though that sizes are the same. If this x were reversed, it could instead end up reading out of bounds.)

But here there is no iteration over m, it’s fixed but the macro doesn’t realise. What ought to work, but doesn’t, is “interpolating” the whole x[m]. It tries, but doesn’t quite unwrap the $ signs correctly. I think that’s a bug:

julia> @tullio a[i,j] = $(x[m])[i,j+0]   (j in 1:J[m])
ERROR: syntax: "$" expression outside quote around /Users/me/.julia/dev/Tullio/src/macro.jl:1001

If fixed that would just be a shorthand way to write your a2 method. Which should also be faster, and LoopVectorization-friendly. I hope that your real problem can be expressed this way?

2 Likes
7

Thanks @mcabbott as usual: your support has been awesome!

Yes, my problem can be expressed that way.

1 Like