TensorCast.jl

This package lets you work with multi-dimensional arrays in index notation, by defining a few macros which translate this to broadcasting, permuting, and reducing operations.

The first is @cast, which deals both with "casting" into new shapes (including going to and from an array-of-arrays) and with broadcasting:

@cast A[row][col] := B[row, col]        # slice a matrix B into rows, also @cast A[r] := B[r,:]

@cast C[(i,j), (k,ℓ)] := D.x[i,j,k,ℓ]   # reshape a 4-tensor D.x to give a matrix

@cast E[φ,γ] = F[φ]^2 * exp(G[γ])       # broadcast E .= F.^2 .* exp.(G') into existing E

@cast _[i] := isodd(i) ? log(i) : V[i]  # broadcast a function of the index values

@cast T[x,y,n] := outer(M[:,n])[x,y]    # generalised mapslices, vector -> matrix function

Second, @reduce takes sums (or other reductions) over the indicated directions. Among such sums is matrix multiplication, which can be done more efficiently using @matmul instead:

@reduce K[_,b] := prod(a,c) L.field[a,b,c]           # product over dims=(1,3), drop dims=3

@reduce S[i] = sum(n) -P[i,n] * log(P[i,n]/Q[n])     # sum!(S, @. -P*log(P/Q')) into exising S

@matmul M[i,j] := sum(k,k′) U[i,k,k′] * V[(k,k′),j]  # matrix multiplication, plus reshape

The same notation with @cast applies a function accepting the dims keyword, without reducing:

@cast W[i,j,c,n] := cumsum(c) X[c,i,j,n]^2           # permute, broadcast, cumsum(; dims=3)

All of these are converted into array commands like reshape and permutedims and eachslice, plus a broadcasting expression if needed, and sum / sum!, or * / mul!. This package just provides a convenient notation.

From version 0.4, it relies on TransmuteDims.jl to handle re-ordering of dimensions, and LazyStack.jl to handle slices. It should also now work with OffsetArrays.jl:

using OffsetArrays
@cast R[n,c] := n^2 + rand(3)[c]  (n in -5:5)        # arbitrary indexing

And it can be used with some packages which modify broadcasting:

using Strided, LoopVectorization, LazyArrays
@cast @strided E[φ,γ] = F[φ]^2 * exp(G[γ])           # multi-threaded
@reduce @turbo S[i] := sum(n) -P[i,n] * log(P[i,n])  # SIMD-enhanced
@reduce @lazy M[i,j] := sum(k) U[i,k] * V[j,k]       # non-materialised

It should work automatically with most array types. This includes GPU arrays such as those from CUDA.jl, whose broadcasting is executed on the device.

Installation

using Pkg; Pkg.add("TensorCast")

The current version requires Julia 1.6 or later. There are a few pages of documentation.

Elsewhere

Similar notation is also used by some other packages, although all of them use an implicit sum over repeated indices. TensorOperations.jl performs Einstein-convention contractions and traces:

@tensor A[i] := B[i,j] * C[j,k] * D[k]      # matrix multiplication, A = B * C * D
@tensor D[i] := 2 * E[i] + F[i,k,k]         # partial trace of F only, Dᵢ = 2Eᵢ + Σⱼ Fᵢⱼⱼ

More general contractions are allowed by OMEinsum.jl, but only one term:

@ein Z[i,j,ξ] := X[i,k,ξ] * Y[j,k,ξ]        # batched matrix multiplication
Z = ein" ikξ,jkξ -> ijξ "(X,Y)              # numpy-style notation

Einsum.jl and Tullio.jl allow arbitrary (element-wise) functions:

@einsum S[i] := -P[i,n] * log(P[i,n]/Q[n])  # sum over n, for each i (also with @reduce above)
@einsum G[i] := 2 * E[i] + F[i,k,k]         # the sum includes everyting:  Gᵢ = Σⱼ (2Eᵢ + Fᵢⱼⱼ)
@tullio Z[i,j] := abs(A[i+x, j+y] * K[x,y]) # convolution, summing over x and y

Notice that @einsum and @tullio sum the entire right hand side, like @reduce does, while @tensor sums individual terms.

These produce very different code for actually doing what you request: The macros @tensor and @ein work out a sequence of basic tensor operations (like contraction and traces), while @einsum and @tullio write the necessary set of nested loops directly (plus optimisations). This package's macros @cast, @reduce and @matmul instead write everything in terms of whole-array operations (like reshape, permutedims and broadcasting).

For those who speak Python, @cast and @reduce allow similar operations to einshape or einops (minus the cool video, but plus broadcasting). In the tests, this file translates many examples. Python's einsum is closer OMEinsum.@ein and TensorOperations.@tensor, and this package's @matmul.

About

This was a holiday project to learn a bit of metaprogramming, originally TensorSlice.jl. But it suffered a little scope creep.