# Unexpected Tullio behavior
**URL:** https://discourse.julialang.org/t/unexpected-tullio-behavior/49371
**Category:** Performance
**Tags:** macros, tullio
**Created:** [October 31, 2020, 3:24pm UTC](https://discourse.julialang.org/t/unexpected-tullio-behavior/49371 "2020-10-31T15:24:53Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)
#### Post date: [October 31, 2020, 3:24pm UTC](https://discourse.julialang.org/t/unexpected-tullio-behavior/49371/1 "2020-10-31T15:24:53Z")
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?
```julia
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()
```
---
### Author:  [@tbeason](https://discourse.julialang.org/u/tbeason)
#### Post date: [October 31, 2020, 3:39pm UTC](https://discourse.julialang.org/t/unexpected-tullio-behavior/49371/2 "2020-10-31T15:39:16Z")
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 (?):
```julia
a[i,j] = sum(x[m],dims=2)
```
I think Tullio is best for heavy matrix algebra where there are several loops to consider.
---
### Author:  [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)
#### Post date: [October 31, 2020, 3:41pm UTC](https://discourse.julialang.org/t/unexpected-tullio-behavior/49371/3 "2020-10-31T15:41:32Z")
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.
---
### Author:  [@tbeason](https://discourse.julialang.org/u/tbeason)
#### Post date: [October 31, 2020, 3:49pm UTC](https://discourse.julialang.org/t/unexpected-tullio-behavior/49371/4 "2020-10-31T15:49:28Z")
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]`
---
### Author:  [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)
#### Post date: [October 31, 2020, 3:51pm UTC](https://discourse.julialang.org/t/unexpected-tullio-behavior/49371/5 "2020-10-31T15:51:25Z")
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.
---
### Author:  [@mcabbott](https://discourse.julialang.org/u/mcabbott)
#### Post date: [October 31, 2020, 4:38pm UTC](https://discourse.julialang.org/t/unexpected-tullio-behavior/49371/6 "2020-10-31T16:38:41Z")
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
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:
```julia
│ 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
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?
---
### Author:  [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)
#### Post date: [October 31, 2020, 5:17pm UTC](https://discourse.julialang.org/t/unexpected-tullio-behavior/49371/7 "2020-10-31T17:17:35Z")
Thanks @mcabbott as usual: your support has been awesome!
Yes, my problem can be expressed that way.