More multi line confusion. Inconsistent local scoping: soft or hard??

[Emmanuel-R8]

A lot of examples with inconsistent results. But bottom line, the behaviour breaks my assumption that macros define a hard local scope (they are hygienic).

What is even weirder is that the macro behaves a mix of soft and hard scope. I am utterly confused.

(Note: I posted the same on Discourse: https://discourse.julialang.org/t/tullio-multiline-syntax-unknown-symbols/77616/5 )

using Tullio

A = rand(4, 4, 5, 5)
B = rand(4, 4)

# OK
@tullio C1[k, l] := A[i, j, k, l] * B[i, j]
C1

# Return a sum instead of 4 x 4 matrix. Intended behaviour?
@tullio C7 := A[i, j, k, l] * B[i, j]

# OK
@tullio C2[k, l] := begin
    A[i, j, k, l] * B[i, j]
end
C1 ≈ C2

# Return a sum instead of 4 x 4 matrix. Intended behaviour?
@tullio C6 := begin
    A[i, j, k, l] * B[i, j]
end


# Breaks because base Julia would expect the matrix to be already existing
@tullio C3[k, l] := begin
    TMP[k, l] = A[i, j, k, l] * B[i, j]
end

# Therefore this works
@tullio C3[k, l] := begin
    TMP = A[i, j, k, l] * B[i, j]
end
C1 ≈ C3

# TMP was just used, but defined in a local scope. Not exported.
# This doesn't work
@tullio C4[k, l] := begin
    TMP[k, l] = A[i, j, k, l] * B[i, j]
    TMP
end
C1 ≈ C4

TMP = zeros(5, 5)
# Breaks with error:
# ERROR: LoadError: MethodError: no method matching zero(::Type{Matrix{Float64}})
@tullio C5[k, l] := begin
    TMP[k, l] = A[i, j, k, l] * B[i, j]
    TMP
end
C1 ≈ C5

TMP = zeros(5, 5)
# But this works
@tullio C5[k, l] := begin
    TMP[k, l] = A[i, j, k, l] * B[i, j]
    TMP[k, l]
end
C1 ≈ C5

# expected a 2-array
@tullio C8 := begin
    TMP[k, l] = A[i, j, k, l] * B[i, j]
    TMP[k, l]
end

# expected a 2-array
@tullio C9 := begin
    TMP[k, l] = A[i, j, k, l] * B[i, j]
    TMP
end


##########################################

A = rand(4, 4, 5, 5)
B = rand(5, 5, 6, 6)
C = rand(6, 6)

@tullio D1[i, j, m, n] := A[i, j, k, l] * B[k, l, m, n]
@tullio E1[i, j] := D1[i, j, m, n] * C[m, n]

# Breaks D2 undefined. Exoected with soft local scope. Not within macro that should define the array.
@tullio F2 := begin
    D2[i, j, m, n] = A[i, j, k, l] * B[k, l, m, n]
    E2[i, j] = D[i, j, m, n] * C[m, n]
end

# Breaks D3 undefined. No idea why
@tullio F3 := begin
    D3 = A[i, j, k, l] * B[k, l, m, n]
    D3[i, j, m, n] * C[m, n]
end

# Bounds error
D4 = zeros(4, 4, 6, 6)
@tullio F4 := begin
    D4 = A[i, j, k, l] * B[k, l, m, n]
    E4 = D4[i, j, m, n] * C[m, n]
end

# Bounds error
D8 = zeros(4, 4, 6, 6)
@tullio F8[i, j] := begin
    D8 = A[i, j, k, l] * B[k, l, m, n]
    E8 = D8[i, j, m, n] * C[m, n]
end

# Bounds error
D9 = zeros(4, 4, 6, 6)
E9 = zeros(4, 4)
@tullio F9[i, j] := begin
    D9 = A[i, j, k, l] * B[k, l, m, n]
    E9 = D9[i, j, m, n] * C[m, n]
end


# Breaks: returns a sum instead of 2 x 2 matrix
D5 = zeros(4, 4, 6, 6)
@tullio F5 := begin
    D5[i, j, m, n] = A[i, j, k, l] * B[k, l, m, n]
    E5 = D5[i, j, m, n] * C[m, n]
end

# Breaks: E6 not defined
D6 = zeros(4, 4, 6, 6)
@tullio F6 := begin
    D6[i, j, m, n] = A[i, j, k, l] * B[k, l, m, n]
    E6[i, j] = D6[i, j, m, n] * C[m, n]
end

# Works! D7 is modified as soft local scope
# D7 needs to be defined outside of scope and captured
# E7 cannot be already defined. See example 9
D7 = zeros(4, 4, 6, 6)
@tullio F7[i, j] := begin
    D7[i, j, m, n] = A[i, j, k, l] * B[k, l, m, n]
    E7 = D7[i, j, m, n] * C[m, n]
end
D7

# Breaks D10 undefined
@tullio F10[i, j] := begin
    D10 = zeros(4, 4, 6, 6)
    D10[i, j, m, n] = A[i, j, k, l] * B[k, l, m, n]
    E10 = D10[i, j, m, n] * C[m, n]
end

[mcabbott]

I strongly urge you not to write into an array within the block, like left[i] := begin TMP[k, l] = .... There is no guarantee that this won't produce garbage, by writing in parallel to the same location. Maybe someone has a safe use for it, but understand that this is a way of hacking this package to do things it doesn't really support, and requires you to think through whether its multi-threading and blocked-access can hurt you.

The basic use is that the right hand side should be a pure function of some arrays of numbers, which is evaluated at every value of every visible index.

It won't really handle arrays of arrays, except maybe arrays of StaticArrays. That's what the zero(::Matrix) errors are saying.

Macros don't have to define a scope, but for Tullio, everything is evaluated within a function, which does. So soft scope weirdness shouldn't apply.

I'm not sure what most of these examples are hoping to achieve. If there are particular formulae you want to encode, then perhaps that's easier than explaining why these don't work?

Or if you want to explore what it does, the "Internals" section writes out an example. Which I hope should make the output of verbose=2 comprehensible; for many of these you will be able to see what it's trying to do and why it fails.

[Emmanuel-R8]

Thanks.

Looking only at the C examples. @tullio C1[k, l] := A[i, j, k, l] * B[i, j] is the one-liner version. All the other C examples aim at replicating that with a begin ... end block.

The first version is easy and works:

@tullio C2[k, l] := begin
    A[i, j, k, l] * B[i, j]
end

After that I explored different versions that seem plausible alternatives, and when it doesn't work understand why and improve the docs. The code I have in a project is about 10 lines of Tullio that I would like to have in a single block to minimise allocations and hopefully better optimisation.

I renumbered the examples.

---

Having read your post on Discourse, the following is clear, although I wasn't certain that indexes that are not repeated should be summed over. But at long as this is the consistent behaviour, all good. Those examples all work and return a scalar.

@tullio C3 := A[i, j, k, l] * B[i, j]

@tullio C4 := begin
    TMP4 = A[i, j, k, l] * B[i, j]
    TMP4
end

@tullio C5 := begin
    A[i, j, k, l] * B[i, j]
end

---

But then this should not work, but C6 is instead a 2x2 matrix!

@tullio C6[k, l] := begin
    TMP6 = A[i, j, k, l] * B[i, j]
end

The allocation to TMP should be a scalar, and should throw an error give C2 is an array.

---

This seems reasonable. First line creates the proper 2x2 matrix, the second returns TMP as an object to be assigned to C9, the same way this is valid.

A = rand(2, 2)
B = A

However, if TMP7 is not already allocated, the error is TMP7 not allocated. If TMP7 is already allocated, the error is no method matching zero(::Type{Matrix{Float64}}). You mention this error occurs when dealing with array of arrays. I don't see where that happens.

@tullio C7 := begin
    TMP7[k, l] = A[i, j, k, l] * B[i, j]
    TMP7
end

If the problem is that Tullio doesn't guess the dimensions of TMP7, the following does not help. TMP8 is not defined.

@tullio C8[k, l] := begin
    TMP8[k, l] = A[i, j, k, l] * B[i, j]
    TMP8
end

---

If the code inside the block is pure Julia, I would have expected this to work:

@tullio C9[k, l] := begin
    TMP9[k, l] = A[i, j, k, l] * B[i, j]
end

Why? Because in Julia an allocation statement returns the allocated value.

function test()
    return a = 2
end

test()
# returns 2

---

Since the errors sometimes refer to temporary matrices not being allocated. For example, this breaks with TMP10 is not defined (same as C8)

@tullio C10[k, l] := begin
    TMP10[k, l] = A[i, j, k, l] * B[i, j]
    TMP10[k, l]
end

Let's have a look when defining.

TMP11 = zeros(5, 5)
@tullio C11[k, l] := begin
    TMP11[k, l] = A[i, j, k, l] * B[i, j]
    TMP11
end

This breaks. Error is no method matching zero(::Type{Matrix{Float64}}.

---

This one is what I mean about soft vs hard scope. Same as C10, but adding indexes to TMP.

TMP11 = zeros(5, 5)
@tullio C11[k, l] := begin
    TMP11[k, l] = A[i, j, k, l] * B[i, j]
    TMP11[k, l]
end

It works giving the same result as the one-liner. BUT:

• TMP11 is modified (not filled with 0.0 anymore). This is soft scope.
• Even stranger, TMP11 is completely different from C11 and seems completely random

[mcabbott]

although I wasn't certain that indexes that are not repeated should be summed over

This is perhaps weird if you expect more Einstein-ish behaviour. It just sums the entire RHS over all free indices. See also #133

but C6 is instead a 2x2 matrix!

@tullio C6[k, l] := begin
    TMP6 = A[i, j, k, l] * B[i, j]
end

This is just like the difference between writing return y and return TMP = y in an ordinary function. The assignment is to a new local scalar variable, but the returned value is identical. C6[k, l] := means it allocates a matrix, and calls this inner block N^2 times.

TMP7 is already allocated

@tullio C7 := begin
    TMP7[k, l] = A[i, j, k, l] * B[i, j]
    TMP7
end

This writes into TMP7 at every iteration. If it decides to multi-thread, it will write different values simultaneously. But then it returns the whole array TMP7, at every i,j,k,l, and tries to sum those to make C7. But (to multi-thread this sum) it has to know the neutral element, and it doesn't. If you put say first(TMP7) there then I think it would run, but the result will still be junk.

This is soft scope.

No, it's not. You can always mutate the contents of a global variable, setindex!(x, ...). What scope controls is whether assignment x = ... makes a new local or not.

stranger, TMP11 is completely different from C11 and seems completely random

There are 4 loops over i, j, k, l, and in the innermost one you write into TMP7. Which element remains depends on what iteration happens last, which for small arrays will be the last, and for larger arrays will be a data race. The result of the innermost expression is what is summed (over i,j) to make C11[k, l], but this happens outside that.

[Emmanuel-R8]

I have read carefully your answers and come up with a list of principles to keep in mind. Please have look to see where I am off. I must be off somewhere since in the examples below, I still don't understand the errors. Note, I have the impression that the errors have little to do with the mistakes in the code. That's macros...

STATEMENT 0:
Summation over all indexes in the RHS unless in the indexes are in the LHS. We don't need to have indexes repeated to have summation.

STATEMENT 1:
Using a block equivalent to using the LHS of @tullio LHS := and evaluating verbatim

LHS = begin
    line 1
    line 2
    line 3
end

The value return by the block is line 3.

STATEMENT 2:
The names of the indexes in LHS have to be consistent with that of the last line. (LHS and line 3). For example if [i, j] in LHS, line 3 would also have [i, j].

STATEMENT 3:
If we have an array A, writing A alone on the last line is equivalent to A[i, j, k, ...]. Therefore full sum over A. (Exception to statement 0)

STATEMENT 4:
Evaluation of lines, example with the final one.
Indices of left and right are compared to see where to sum in line 3. Any index not in LHS will be summed over.
Line 3 is evaluated, then the dimensions of the result (line 3) is used to check the dimensions of LHS.

STATEMENT 4 bis:
However, what are indexes on the left and indices on the right ? If the final line is in the form A = B, the evaluation is LHS := (A = B) which should be equivalent to LHS := B
I don't know the exact impact mismatch dimensions between A, B and C. For example C[k, l] := (A = B[i, j, k, l]). Which dimensions are compared, which one raise an error? (see 6 and 6_2)

STATEMENT 5:
Julia macros are not hygienic by default. Same behaviour as function unless local var is used. Tullio allows for environment capture. Julia is Common Lisp by default, and not Scheme. (I thought the opposite...)

---

Here is the list of tests. I understand more than before, but I still don't understand 6, 6_2, 7, 9, 10. I know you explained 7, but I still don't get it.

# OK. Statement 0
@tullio C2[k, l] := begin
    A[i, j, k, l] * B[i, j]
end

# OK. Statement 0
@tullio C3 := A[i, j, k, l] * B[i, j]

# OK. Statement 0
@tullio C4 := begin
    TMP4 = A[i, j, k, l] * B[i, j]
    TMP4
end

# OK. Statement 0
@tullio C5 := begin
    A[i, j, k, l] * B[i, j]
end

# The result is a 5x5 matrix. 
# Problem of Statement 4 bis.
# TMP6 = A[i, j, k, l] * B[i, j] alone means summing over everything and getting a scalar. But this conflicts with `C6[k, l]`
# This suggests that TMP6 is ignored.
@tullio C6[k, l] := begin
    TMP6 = A[i, j, k, l] * B[i, j]
end

# BUT this breaks with TMP6 not defined. TMP6 clearly has an impact negating the previous conclusion.
@tullio C6_2[k, l] := begin
    TMP6[k, l] = A[i, j, k, l] * B[i, j]
end

# Error TMP7 not defined. 
# Statement 3: TMP7 means sum(TMP7) which is a scalar. 
# C7 is also a scalar. 
# I don't understand
@tullio C7 := begin
    TMP7[k, l] = A[i, j, k, l] * B[i, j]
    TMP7
end

# Statement 3: TMP8 means sum(TMP8) which is a scalar.
# Breaks because of Statement 4
@tullio C8[k, l] := begin
    TMP8[k, l] = A[i, j, k, l] * B[i, j]
    TMP8
end

# TMP9 not defined. What has TMP9 to do with the error. Dimensions are correct. 
@tullio C9[k, l] := begin
    TMP9[k, l] = A[i, j, k, l] * B[i, j]
end

# TMP10 not defined. Ditto. Why?
@tullio C10[k, l] := begin
    TMP10[k, l] = A[i, j, k, l] * B[i, j]
    TMP10[k, l]
end

# Unhygienic macro. Can capture TMP11
TMP11 = zeros(5, 5)
@tullio C11[k, l] := begin
    TMP11[k, l] = A[i, j, k, l] * B[i, j]
    TMP11
end

# All further examples make sense.
TMP11 = zeros(5, 5)
@tullio C11[k, l] := begin
    TMP11[k, l] = A[i, j, k, l] * B[i, j]
    TMP11[k, l]
end

TMP12 = zeros(5, 5)
@tullio C12[k, l] := begin
    TMP12[k, l] = A[i, j, k, l] * B[i, j]
end

TMP13 = zeros(5, 5)
@tullio C13[k, l] := begin
    TMP13[k, l] = A[i, j, k, l] * B[i, j]
    TMP13
end