Not an issue, just a general question

[Jutho]

Dear Tim,

I have some questions regarding Cartesian. They might be too specific to discuss in the Julia forum, but also of sufficiently wide interest not to discuss them via private communication via email. So I hope you don't mind I abuse the Github issue functionality of your package for this purpose.

1. I assume most of Cartesian.jl package has now become redundant and is replaced by the functionality in Base.Cartesian? This is then used extensively to define the functions in "base/multidimensional.jl". So my questions will actually about Base.Cartesian instead, which is clearly a fascinating piece of programming.

2. I would like to understand better both the motivation and the inner workings of some of the implementations in multidimensional.jl, which relies on the macros and functions in cartesian.jl. Do you have some more advanced documentation about this? As somebody with very little metaprogramming experience, I have a hard time deciphering how all of this works exactly.

3. In particular, taking permutedims! as an example:
   a) Let us first the motivation. In principle, this could be written using a single for loop that uses linear indexing over the destination (this is done anyway) and for every index computes the corresponding indices or linear index for the source array, by doing operations on a vector of strides and counts. If I understand correctly, the different macros basically store the strides in individual local variables stride_1 up to stride_{N+1}, at which point the function body becomes dependent on N. Is the motivation for this efficiency, i.e. that it is much quicker to act on local variables then on the items of an array? (I have only a very basic understanding of how computer architecture and function calls work).

b) While I think I have some understanding of what @ngenerate is doing, I will not try to guess but rather would like to ask if you could briefly explain what it is doing exactly with building a cache dictionary? How does this relate to the @nsplat which is also defined in cartesian.jl

c) I think I understand what @nLoops is doing, and even how it is doing it (which is very nice by the way). Would there be a huge computation overhead if you just try to use the same mechanism (i.e. building a variable number of for loops by building an expr that is recursively grown) but would not have the advanced @ngenerate framework etc that automatically expands these expressions again and directly build the corresponding function? In that case you would probably be building the expr object in every function call?

3. I would like to apply the machinery of cartesian.jl in the context of tensor contractions and tensor traces, where there would be up to three independent N values that would need to be dispatched on and would require a variable number of for loops (i.e a loop over N1 indices that need to be contracted, a loop over N2 uncontracted indices of the first tensor and a loop over N3 uncontracted indices of the second tensor). Is this possible with the current macros or would this require a further extension of the functionality in cartesian.jl ?

Thanks,

Jutho

[timholy]

Do you have some more advanced documentation about this?

Not yet about the implementations in multidimensional; the only documentation for Cartesian is in the README. I plan to move this to "developer documentation" when that gets created (or create it myself).

Is the motivation for this efficiency, i.e. that it is much quicker to act on local variables then on the items of an array?

Yes, accessing/modifying items of an array requires a lookup to/from memory; stack variables save that.

briefly explain what it is doing exactly with building a cache dictionary

I wish it didn't have to, but it's currently necessary. The idea is (1) you have to generate a different version for each dimensionality, (2) code generation is slow, and (3) if you don't want to pre-generate versions for arbitrary dimensions (e.g., it probably doesn't make sense to pre-generate one for 27-dimensional arrays), then you need a mechanism to generate them on the fly if necessary but cache the result so it can be reused later. JuliaLang/julia#5395, if it gets implemented, will allow us to use Julia's internal compiler-aware method caches instead.

Would there be a huge computation overhead if you just try to use the same mechanism ... but would not have the advanced @ngenerate framework etc that automatically expands these expressions again and directly build the corresponding function?

That's fine. For my own use I often write code like the first example in this section.

I would like to apply the machinery of cartesian.jl in the context of tensor contractions and tensor traces, where there would be up to three independent N values that would need to be dispatched on and would require a variable number of for loops (i.e a loop over N1 indices that need to be contracted, a loop over N2 uncontracted indices of the first tensor and a loop over N3 uncontracted indices of the second tensor). Is this possible with the current macros or would this require a further extension of the functionality in cartesian.jl ?

Your best model would be in base/broadcast.jl, where there is a loop over both dimensionality and narray. One of the differences between the version in the Cartesian package and what's in base is that in base the variables are called i_1, i_2, ... instead of i1, i2. The reason for the change was to allow multi-dimensional iteration, e.g., iter_i_j.

[timholy]

There's another good example, by the way, in JuliaLang/julia#5857. I'm not sure, but it may be closer to your application.