 # Unexpected Tullio behavior

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()
``````

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.

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.

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]`

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.

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

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

Yes, my problem can be expressed that way.

1 Like