### Install and Use Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Installs the Tullio package and demonstrates basic index notation for summation and broadcasting.

```julia
using Pkg; Pkg.add("Tullio")
using Tullio, Test
M = rand(1:20, 3, 7)

@tullio S[1,c] := M[r,c]  # sum over r ∈ 1:3, for each c ∈ 1:7
@test S == sum(M, dims=1) 

@tullio Q[ρ,c] := M[ρ,c] + sqrt(S[1,c])  # loop over ρ & c, no sum -- broadcasting
@test Q ≈ M .+ sqrt.(S)
```

--------------------------------

### Product Reduction with @tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Perform a product reduction over specified indices. This example calculates a product over 'k' and assigns it to 'Z'.

```julia
@tullio (*) Z[j] := X[ind[k],j] * exp(-Y[k])
```

--------------------------------

### Tullio Keyword Options Example

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Illustrates various keyword options for the @tullio macro, including threading, AVX vectorization, gradient calculation, and array assignment methods.

```julia
@tullio threads=true fastmath=true avx=true tensor=true cuda=256 grad=Base verbose=false A[i,j] := ...
```

--------------------------------

### Compute gradients for matrix multiplication

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

This example demonstrates computing gradients for a matrix multiplication defined using `@tullio`. It shows how to obtain gradients with respect to both input matrices.

```julia
using Tullio, Zygote

# Matrix multiply gradient
mm(x, y) = @tullio z[i, j] := x[i, k] * y[k, j]
x1 = rand(3, 4)
y1 = rand(4, 5)
dx, dy = gradient(sum ∘ mm, x1, y1)
# dx ≈ ones(3,5) * y1'
# dy ≈ x1' * ones(3,5)
```

--------------------------------

### Tullio with OffsetArrays

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Example demonstrating the use of Tullio in conjunction with OffsetArrays for handling arrays with non-standard indexing.

```julia
using Tullio, OffsetArrays
```

--------------------------------

### Compute min/max gradients

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

This example shows how to compute gradients for reduction operations like `max`. Note that at ties, a subgradient is returned.

```julia
using Tullio, Zygote

# Min/max gradients (subgradient at ties)
f1(x) = @tullio (max) z = x[i]
gradient(f1, [1, 2, 3, 4])  # returns ([0, 0, 0, 1],)
```

--------------------------------

### Matrix multiplication with gradient support

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

This example shows how to perform matrix multiplication using Tullio's `@tensor` macro and obtain gradients using Zygote.jl. The `@tensor` macro automatically provides gradient definitions.

```julia
using Tullio, TensorOperations, Zygote

# With gradient support
gradient(sum, Tullio.@tensor(C[i,j] := A[i,k] * B[k,j]))
```

--------------------------------

### Product over gathered elements

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

This example demonstrates a product reduction over gathered elements, combining gather operations with a reduction across a specific dimension.

```julia
using Tullio

ind = vcat(1:3, 1:3)
V = 1:6
@tullio (*) Z[j] := M[ind[k], j] * exp(-V[k])
# Product over k dimension
```

--------------------------------

### Compute gradient for log function with division

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

This example demonstrates computing gradients for a function involving `log` and division, defined using `@tullio`. Gradients are computed for both input matrices.

```julia
using Tullio, Zygote

# Log function with gradient
flog(x, y) = @tullio z[i] := log(x[i, j]) / y[j, i]
r_x, r_y = rand(2, 3), rand(3, 2)
gradient(sum ∘ flog, r_x, r_y)
```

--------------------------------

### Tullio.@einsum Macro Usage

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Example of using the Tullio.@einsum macro, a variant designed to allow running Einsum.jl's tests.

```julia
Tullio.@einsum
```

--------------------------------

### Broadcasting with Reduction using @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Perform broadcasting combined with a reduction operation. This example sums over the `r` dimension for each `c`.

```julia
# Broadcasting with reduction - sum over y dimension
M = rand(3, 7)
@tullio S[1, c] := M[r, c]  # sum over r for each c
# Equivalent to: S == sum(M, dims=1)
```

--------------------------------

### Scalar Reduction Example

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Illustrates a Tullio.jl macro for a scalar reduction operation, where the result is a single value computed from a tensor contraction.

```julia
@tullio s := A[i,j] * log(B[j,i])
```

--------------------------------

### Custom initialization for reductions

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Specify a custom initial value for reduction operations using the `init` keyword argument. This is particularly useful for operations like `max` or `min` where a specific starting value is required.

```julia
using Tullio

# Custom initialization for reductions
@tullio (max) m := A[i]  init=0.0
```

--------------------------------

### Basic Array Creation with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Use `@tullio` with `:=` to create new arrays from index expressions. This example demonstrates squaring elements of a range.

```julia
using Tullio

# Create array from index expression
@tullio A[i] := (1:10)[i]^2
# A == [1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
```

--------------------------------

### Compute simple gradients with Zygote

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Tullio integrates with Zygote.jl for automatic differentiation. This example shows how to compute the gradient of a simple sum operation with respect to its input vector.

```julia
using Tullio, Zygote

# Simple gradient
x = rand(3)
gradient(x -> sum(@tullio y[i] := 2*x[i]), x)
# Returns ([2, 2, 2],)
```

--------------------------------

### Batched Matrix Multiplication with Unusual Index Layout

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

This example shows batched matrix multiplication with a non-standard index layout, summing over `j` to avoid explicit `permutedims`.

```julia
# Batched matrix multiplication with unusual index layout
A = randn(20, 30, 500)
B = randn(500, 40, 30)
@tullio C[i, k, b] := A[i, j, b] * B[b, k, j]
# Sum over j, avoiding expensive permutedims
```

--------------------------------

### Reduction and Assignment with @tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Perform reductions (e.g., summation) over specific indices and assign the result to a target array. This example sums over 'y' and writes the result into 'S'.

```julia
@tullio S[x] = P[x,y] * log(Q[x,y] / R[y])
```

--------------------------------

### Tullio.@tensor Macro Usage

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Example of using the Tullio.@tensor macro, which leverages TensorOperations.jl for expression evaluation and provides gradient definitions.

```julia
Tullio.@tensor
```

--------------------------------

### Compute gradient for scalar output function

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

This example shows how to compute the gradient of a function that produces a scalar output using `@tullio`. The gradient is computed with respect to the input vector.

```julia
using Tullio, Zygote

# Scalar output gradient
s2(x) = @tullio s := exp(x[i]) / x[j]
gradient(s2, rand(10))
```

--------------------------------

### Function Application Outside Summation with @tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Use pipe operators '|>' or '<|' to apply functions to the result *after* the summation is complete. This example calculates the log-sum-exp.

```julia
@tullio lse[j] := log <| exp(mat[i,j])
```

--------------------------------

### Accumulation with @tullio and +=

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Perform in-place accumulation using `+=` with `@tullio`. This example sums over the `j` dimension and adds the result to array `B`.

```julia
# Accumulation with +=
B = copy(A)
D .= 3
@tullio B[i] += D[i, j]  # sum over j and add to B
```

--------------------------------

### In-place Array Update with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Use `@tullio` with `=` for in-place updates of an existing array. This example adds 100 to each element based on indices.

```julia
# In-place update
D = similar(A, 10, 10)
@tullio D[i, j] = A[i] + 100  # write into existing array
```

--------------------------------

### Tullio Expression with Explicit Index Ranges and Dual Numbers

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

An example of a more complex Tullio expression where output index ranges must be explicitly provided. It also demonstrates the use of `grad=Dual` for automatic differentiation.

```julia
@tullio out[x, y] := @inbounds(begin  # sum over k
        a,b = off[k]
        mat[mod(x+a), mod(y+b)]
    end) (x in axes(mat,1), y in axes(mat,2)) grad=Dual nograd=off
```

--------------------------------

### Tullio Cost Model and Fast Operations

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Explains the Tullio cost model for operations and shows how to use the 'threads=false' option to avoid allocation for fast operations, ensuring zero memory allocation after warmup.

```julia
# Cost model influences threading decisions
# Expensive operations (log, exp) have higher cost
Tullio.COSTS  # Dict(:+ => 0, :log => 10, :exp => 10, ...)

# Fast operations avoid allocation
m!(C, A, B) = @tullio C[i, k] = A[i, j] * B[j, k]  threads=false
C, A, B = rand(4, 4), rand(4, 4), rand(4, 4)
@allocated m!(C, A, B)  # 0 bytes after warmup
```

--------------------------------

### Matrix Multiplication with Gradients using Tullio and Tracker

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates how to compute matrix multiplication using Tullio and then calculate its gradients with respect to the input matrices using the Tracker library.

```julia
using Tullio
mul(A, B) = @tullio C[i,k] := A[i,j] * B[j,k]
```

```julia
A = rand(3,40); B = rand(40,500);
```

```julia
A * B ≈ mul(A, B) # true
```

```julia
using Tracker # or Zygote
ΔA = Tracker.gradient((A,B) -> sum(mul(A, B)), A, B)[1]
```

```julia
ΔA ≈ ones(3,500) * B' # true
```

--------------------------------

### FFT Implementation with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Shows how to implement a Discrete Fourier Transform using Tullio's index notation and complex exponentials.

```julia
using FFTW # Functions of the indices are OK:
S = [0,1,0,0, 0,0,0,0]
fft(S) ≈ @tullio F[k] := S[x] * exp(-im*pi/8 * (k-1) * x)  (k ∈ axes(S,1))
```

--------------------------------

### Standard Matrix Multiplication with Einsum

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Demonstrates standard matrix multiplication using the @einsum macro for tensor contraction.

```julia
using Tullio

# Standard matrix multiply syntax
@einsum C[i, j] := A[i, k] * B[k, j]
```

--------------------------------

### Calculating the Trace with Einsum

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Illustrates how to compute the trace of a matrix using the @einsum macro by contracting diagonal elements.

```julia
@einsum tr = A[i, i]
```

--------------------------------

### Einsum.jl Macro Expansion

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Illustrates the explicit for-loop expansion of the Einsum.jl macro, showing how it generates manual loops for computation.

```julia
@einsum C[i,j] := A[i,k] * B[k,j]         # Einsum.jl
```

--------------------------------

### Asymmetric Padding with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Apply asymmetric padding using `pad(i, before, after)` to specify different padding amounts on each side.

```julia
# Asymmetric padding
A = [i^2 for i in 1:10]
@tullio G[i+_] := A[pad(i, 2, 3)]  # 2 before, 3 after
# G == [0, 0, 1, 4, 9, ..., 81, 100, 0, 0, 0]
```

--------------------------------

### Custom Reduction Functions with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Illustrates defining and using custom reduction functions (e.g., product, custom min) with Tullio's notation.

```julia
# Reduction over any function:
@tullio (*) P[i] := A[i+k]  (k in 0:2) # product
@tullio (max) X[i,_] := D[i,j]         # maximum(D, dims=2), almost

min1(x,y) = ifelse(first(x) < first(y), x, y); # findmin(D, dims=1), almost:
@tullio (min1) Ts[j+_] := (D[i,j], (i,j))  init=(typemax(Int), (0,0))
```

--------------------------------

### Multi-line Expressions with Begin-End Blocks

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Shows how to use begin...end blocks within @tullio for complex operations, intermediate assignments, and custom index handling.

```julia
using Tullio, OffsetArrays

A = rand(10, 10)
offsets = [(a, b) for a in -2:2 for b in -2:2 if a >= b]

# Stencil with tuple unpacking
@tullio out[x, y] := begin
    a, b = offsets[k]
    i = clamp(x+a, extrema(axes(A, 1))...)
    @inbounds A[i, clamp(y+b), 1] * 10
end
```

```julia
# Convolution with cyclic boundaries
mat = zeros(10, 10, 1)
mat[2, 2] = 101
kern = @tullio _[i, j] := 1/(1+i^2+j^2)  (i in -3:3, j in -3:3)

@tullio out[x, y, c] := begin
    xi = mod(x+i, axes(mat, 1))
    @inbounds trunc(Int, mat[xi, mod(y+j), c] * kern[i, j])
end (x in 1:10, y in 1:10)
```

```julia
# Applying vector of functions
fs = [sin, cos, tan]
xs = randn(3, 100)
@tullio ys[r, c] := (fs[r])(xs[r, c])
```

--------------------------------

### Reduction with Custom Initialization using @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Specify a custom initial value for reductions using the `init=` keyword argument.

```julia
# Maximum with custom initialization
@tullio (max) m := A[i]  init=200
# m == 200 (init larger than all elements)
```

--------------------------------

### Broadcast Reductions with Tullio for Efficiency

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Shows how Tullio can optimize broadcast reductions, avoiding large allocations and speeding up operations like log. Compares Tullio with standard Julia broadcasting and reductions.

```julia
sum_opp(X, Y=X) = @tullio s := X[i,j] * log(Y[j,i])
```

```julia
sum_part(X, Y=X) = @tullio S[i] := X[i,j] * log(Y[j,i])
```

```julia
X = rand(1000,1000);
```

```julia
@btime sum_opp($X)
```

```julia
@btime sum($X .* log.(transpose($X)))
```

```julia
@btime sum_part($X)'
```

```julia
@btime sum($X .* log.(transpose($X)), dims=2)
```

--------------------------------

### Convolution Implementations with Tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Illustrates different convolution implementations using Tullio, including those with padding, and compares their performance against standard Julia libraries like Flux and DSP.

```julia
conv1(x,k) = @tullio y[i+_, j+_] := x[i+a, j+b] * k[a,b]
```

```julia
conv2(x,k) = @tullio y[i+_, j+_] := x[2i-a, 2j-b] * k[a,b]
```

```julia
conv3(x,k) = @tullio y[i+_, j+_] := x[pad(i-a,3), pad(j-b,3)] * k[a,b]
```

```julia
x100 = rand(100,100); k7 = randn(7,7);
```

```julia
@btime conv1($x100, $k7);
```

```julia
@btime conv2($x100, $k7);
```

```julia
@btime conv3($x100, $k7);
```

```julia
using Flux
x104 = reshape(x100,(100,100,1,1)); k74 = reshape(k7,(7,7,1,1));
```

```julia
conv1(x100, k7) ≈ @btime CrossCor($k74, false)($x104)
```

```julia
conv2(x100, k7) ≈ @btime Conv($k74, false, stride=2)($x104)
```

```julia
conv3(x100, k7) ≈ @btime Conv($k74, false, pad=3)($x104)
```

```julia
using DSP
@btime DSP.conv($x100, $k7);
```

--------------------------------

### In-place Product Accumulation with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Perform in-place product accumulation using `*=` with `@tullio`, calculating cumulative products.

```julia
# In-place product with *=
C = copy(A)
@tullio (*) C[k] = ifelse(i <= k, A[i], 1)
# C == cumprod(A)
```

--------------------------------

### Zero-Padded Indexing with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Utilize `pad(i)` for zero-padded boundary conditions, introducing zeros outside the original array bounds.

```julia
# Zero-padded boundary
@tullio E[i] := B[pad(i)]  (i in -2:7)
# E has zeros outside original range
```

--------------------------------

### Tullio Performance Tuning Constants

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Demonstrates how to check and adjust Tullio's BLOCK and TILE constants for threading and cache optimization. These constants influence threading decisions based on array size and memory usage.

```julia
using Tullio

# Check current settings
Tullio.BLOCK[]  # default 2^18, threshold for threading
Tullio.TILE[]   # default 512 bytes, tile size for cache

# Adjust for testing/small arrays
Tullio.BLOCK[] = 32
Tullio.TILE[] = 32
```

--------------------------------

### Reduction over min/max with Tullio Gradients

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates gradient calculation for a reduction operation (maximum) using Tullio. Note that complete reductions to a scalar may not work on GPU and gradients are only calculated for arrays.

```julia
Tracker.gradient(x -> (@tullio (max) res := x[i]^3), [1,2,3,-2,-1,3])[1]
```

--------------------------------

### GPU Accelerated Matrix Multiplication with Tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Shows how to use Tullio for matrix multiplication on the GPU by leveraging CUDA and KernelAbstractions. Verifies correctness and gradient calculation on the GPU.

```julia
using CUDA, KernelAbstractions # Now defined with a GPU version:
mul(A, B) = @tullio C[i,k] := A[i,j] * B[j,k]
```

```julia
cu(A * B) ≈ mul(cu(A), cu(B)) # true
```

```julia
cu(ΔA) ≈ Tracker.gradient((A,B) -> sum(mul(A, B)), cu(A), cu(B))[1] # true
```

--------------------------------

### Dimension Names and Pretty Printing with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Illustrates using NamedDims and AxisKeys with Tullio for dimension naming and improved output formatting.

```julia
using NamedDims, AxisKeys # Dimension names, plus pretty printing:
@tullio M[row=i, col=j, z=k] := A[i+j-1]  (j in 1:15, k in 1:2)
@tullio S[i] := M[col=j-i, z=k, row=i+1] # sum over j,k
```

--------------------------------

### In-place Updates with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates modifying an existing matrix in-place using Tullio's assignment syntax.

```julia
dbl!(M, S) = @tullio M[r,c] = 2 * S[1,c]  # write into existing matrix, M .= 2 .* S
dbl!(M, S)
@test all(M[r,c] == 2*S[1,c] for r ∈ 1:3, c ∈ 1:7)
```

--------------------------------

### Outer Product Calculation with Einsum

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Demonstrates the calculation of an outer product between two vectors using the @einsum macro.

```julia
v = rand(5)
@einsum O[i, j] := v[i] * v[j]
```

--------------------------------

### Reduction Operations with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Shows how to perform reduction operations like maximum using Tullio's index notation.

```julia
@tullio (max) X[i] := abs2(T[j,i,δ])  # reduce using max, over j and δ
@test X == dropdims(maximum(abs2, T, dims=(1,3)), dims=(1,3))
```

--------------------------------

### In-place Matrix Multiplication with Einsum

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Shows how to perform in-place matrix multiplication using the @einsum macro with the += operator.

```julia
C = zeros(3, 3)
@einsum C[i, j] += A[i, k] * B[k, j]
```

--------------------------------

### Finalizers and Reductions with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates the use of finalizer operators like '<|' and '|>' for post-computation operations and various reduction functions.

```julia
# Finalisers <| or |> are applied after sum (the two are equivalent):
@tullio N2[j] := sqrt <| M[i,j]^2     # N2 ≈ map(norm, eachcol(M))
@tullio n3[_] := A[i]^3  |> (_)^(1/3) # n3[1] ≈ norm(A,3), with _ anon. func.
```

--------------------------------

### Control gradient computation with @tullio options

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Control gradient computation using `grad=false` to disable it entirely, or `grad=Dual` to use ForwardDiff.jl. The default is `grad=Base` which uses Zygote.jl.

```julia
using Tullio

# Disable gradient computation
@tullio A[i] := (1:10)[i]^2  grad=false

# Use ForwardDiff for gradients
@tullio A[i] := (1:10)[i]^2  grad=Dual
```

--------------------------------

### Index Shifts with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates how to handle index shifts in array operations using Tullio's notation, including automatic shifts.

```julia
# Shifts -- range of i calculated in terms of that given for j:
@tullio M[i,j] := A[i+j-1]  (j in 1:15)  # i in 1:7
@tullio M[i+_,j] := A[i+j]  (j in 1:15)  # i in 0:6, automatic shift "i+_"
```

--------------------------------

### Manual For-Loop Implementation

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Provides a direct, manual implementation of matrix multiplication using nested for-loops, similar to what Einsum.jl might expand to.

```julia
T = promote_type(eltype(A), eltype(B))
C = Array{T}(undef, size(A,1), size(B,2))
@inbounds for j in 1:size(B,2)
    for i in 1:size(A,1)
        acc = zero(T)
        for k in 1:size(A,2)
            acc += A[i,k] * B[k,j]
        end
        C[i,j] = acc
    end
end
```

--------------------------------

### Tullio.jl Performance Comparison

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Compares the performance of Tullio.jl with OpenBLAS for matrix multiplication, highlighting its speed for tensor contractions.

```julia
using Tullio, LoopVectorization, NNlib, BenchmarkTools
```

--------------------------------

### Tensor Operations Macro Comparison

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Compares the syntax of different Julia packages for tensor operations, highlighting their similarities and differences.

```julia
@tensor C[i,j] := A[i,k] * B[k,j]         # TensorOperations.jl
```

```julia
@ein C[i,j] := A[i,k] * B[k,j]            # OMEinsum.jl
```

```julia
@matmul C[i,j] := sum(k) A[i,k] * B[k,j]  # TensorCast.jl
```

--------------------------------

### Concatenation and Broadcasting with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates concatenation of arrays and addition of a constant using Tullio's index notation.

```julia
R = [rand(Int8, 3, 4) for δ in 1:5]

@tullio T[j,i,δ] := R[δ][i,j] + 10im  # three nested loops -- concatenation
@test T == permutedims(cat(R...; dims=3), (2,1,3)) .+ 10im
```

--------------------------------

### Index by Values with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates indexing into an array using values from another array (K) within Tullio's notation.

```julia
# Index by the values in K
@tullio D[i,j] := A[2K[j]+i] ÷ K[j] # extrema(K)==(-1,2) implies i ∈ 3:17
```

--------------------------------

### Complex Matrix Expression with Tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates a more complex matrix expression involving subtractions and transposes, comparing a Tullio implementation with a sequential approach.

```julia
M, Σ = randn(100,100), randn(100,100);
```

```julia
@tullio R4[i, j] := (M[μ, i] - M[μ,j])' * Σ[μ,ν] * (M[ν, i] - M[ν, j]);
```

```julia
begin
  S = M' * Σ * M  # two N^3 operations, instead of one N^4
  @tullio R3[i,j] := S[i,i] + S[j,j] - S[i,j] - S[j,i]
end;
```

```julia
R3 ≈ R4
```

--------------------------------

### Wrapped and Padded Arrays with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Illustrates handling wrapped (modulo) and padded array access using Tullio's notation.

```julia
# Wrapped & padded:
@tullio M[i,j] := A[mod(i+j)]  (j in 1:15, i in 1:15)   # wraps around, mod(i+j, axes(A,1))
@tullio M[i,j] := A[clamp(i+j)]  (j in 1:15, i in 1:15) # instead repeats "100"
@tullio M[i+_,j] := A[pad(i+j, 3)]  (j in 1:15)         # fills with zeros
```

--------------------------------

### Automatic Differentiation with Tullio and Zygote

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates automatic differentiation using Tullio's `grad=Dual` option and Zygote.jl. It defines a `rowmap` function that computes gradients and then uses Zygote to compute the gradient of a composed function.

```julia
using Zygote, ForwardDiff
rowmap(fs, xs) = @tullio ys[r,c] := (fs[r])(xs[r,c]) grad=Dual nograd=fs
Zygote.gradient(sum∘rowmap, fs, ones(3,2))
```

--------------------------------

### Enable verbose output for debugging

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Use `verbose=true` or `verbose=<level>` to enable detailed output during `@tullio` macro expansion, which is helpful for debugging performance issues or understanding the generated code.

```julia
using Tullio

# Verbose output for debugging
@tullio A[i] := (1:10)[i]^2  verbose=true
@tullio A[i] := (1:10)[i]^2  verbose=2  # maximum verbosity
```

--------------------------------

### Product Reduction with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Perform a product reduction over an array using the `(*)` operator with `@tullio`.

```julia
using Tullio

A = [i^2 for i in 1:10]

# Product reduction
@tullio (*) P[_] := float(A[i])
# P == [prod(A)]
```

--------------------------------

### Control vectorization with @tullio options

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Control LoopVectorization.jl usage with `avx=false` to disable it or `avx=<unroll_factor>` to set the unroll factor for AVX instructions.

```julia
using Tullio

# Disable LoopVectorization
@tullio A[i] := (1:10)[i]^2  avx=false

# Set AVX unroll factor
@tullio A[i] := (1:10)[i]^2  avx=4
```

--------------------------------

### Tullio One-Based Indexing Adjustment

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Shows how to use `_` in the left-hand side of a @tullio assignment to automatically insert the necessary shift for one-based indexing when creating a new array.

```julia
A[i+_] := 
```

--------------------------------

### Chained Matrix Multiplication Comparison

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Compares the performance of chained matrix multiplication using standard Julia operators versus a Tullio implementation. Tullio can be slower if it doesn't recognize a more efficient algorithm.

```julia
M1, M2, M3 = randn(30,30), randn(30,30), randn(30,30);
```

```julia
@btime $M1 * $M2 * $M3;
```

```julia
@btime @tullio M4[i,l] := $M1[i,j] * $M2[j,k] * $M3[k,l];
```

--------------------------------

### Array Creation Functions in Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates that functions creating arrays (like rand) are evaluated once per operation within Tullio's notation.

```julia
# Functions which create arrays are evaluated once:
@tullio R[i,j] := abs.((rand(Int8, 5)[i], rand(Int8, 5)[j]))
```

--------------------------------

### Gradient Kernel with ForwardDiff Dual Numbers

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates an alternative gradient kernel implementation using ForwardDiff's Dual numbers, suitable for cases without finalizers.

```julia
function ∇act!(::Type, ΔC, ΔA, ΔB, C, A, B, ax_i, ax_j, ax_k, keep)
    eps1 = ForwardDiff.Dual(0, (1,0))
    eps2 = ForwardDiff.Dual(0, (0,1))
    for k in ax_k, i in ax_i, j in ax_j
        res = (A[i,k] + eps1) * (B[k,j] + eps2)
        ΔA[i,k] += ForwardDiff.partials(res, 1) * ΔC[i,j]
        ΔB[k,j] += ForwardDiff.partials(res, 2) * ΔC[i,j]
    end
end
```

--------------------------------

### Convolution with OffsetArrays and Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Shows how to perform convolution using OffsetArrays and Tullio's index notation.

```julia
using OffsetArrays # Convolve a filter:
K = OffsetArray([1,-1,2,-1,1], -2:2)
@tullio C[i] := A[i+j] * K[j]    # j ∈ -2:2 implies i ∈ 3:19
```

--------------------------------

### Batched Matrix Multiplication with Tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Performs batched matrix multiplication where the batch index is first in the second matrix. Compares performance against NNlib.

```julia
bmm_rev(A, B) = @tullio C[i,k,b] := A[i,j,b] * B[b,k,j]  # (sum over j)
```

```julia
A = randn(20,30,500); B = randn(500,40,30);
```

```julia
bmm_rev(A, B) ≈ NNlib.batched_mul(A, permutedims(B, (3,2,1)))  # true
```

```julia
@btime bmm_rev($A, $B);
```

```julia
@btime NNlib.batched_mul($A, permutedims($B, (3,2,1)));
```

--------------------------------

### 2D Padded Convolution with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Demonstrates a 2D convolution using padded array access for handling boundaries.

```julia
# 2D padded convolution
mat = rand(10, 10)
kern = randn(7, 7)
@tullio H[i+_, j+_] := mat[pad(i, 3), pad(j, 3)] * kern[i, j]
```

--------------------------------

### Tullio Function to Create and Execute Operation

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Illustrates the helper function generated by Tullio.jl that sets up the necessary axes, determines the return type, and calls the core computation function.

```julia
function make(A, B)
    ax_i = axes(A,1)
    ax_j = axes(B,2)
    ax_k = axes(A,2) # and check this is == axes(B,1)
    rhs(A,B,i,j,k) = tanh(A[i,k] * B[k,j])
    T = Core.Compiler.return_type(rhs, eltype.((A,B,1,1,1))) # plus a fallback
    C = similar(A, T, (ax_i, ax_j))
    Tullio.threader(act!, Array{T}, C, (A,B), (ax_i,ax_j), (ax_k,), +, 64^3)
    return C
end
```

--------------------------------

### Control threading with @tullio options

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Control multithreading behavior using `threads=false` to disable threading or `threads=<threshold>` to set a minimum number of iterations before threading is used.

```julia
using Tullio

# Disable threading
@tullio A[i] := (1:10)[i]^2  threads=false

# Set threading threshold (number of iterations)
@tullio A[i] := (1:10)[i]^2  threads=64^3
```

--------------------------------

### Summation and Array Creation with @tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Use this macro to perform summations over specified indices and create new arrays. It infers ranges for indices not present on the left-hand side.

```julia
@tullio M[x,y,c] := N[x+i, y+j,c] * K[i,j]
```

--------------------------------

### Modulo Indexing for Cyclic Boundaries

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Use `mod(i)` within `@tullio` for cyclic or periodic boundary conditions, wrapping indices around.

```julia
using Tullio

B = 1:5

# Cyclic/periodic boundary - wraps around
@tullio C[i] := B[mod(i)]  (i in 1:10)  avx=false
# C == [1, 2, 3, 4, 5, 1, 2, 3, 4, 5]
```

--------------------------------

### Tullio Padding Indices

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates the use of `pad` for extending index ranges in @tullio expressions, inserting zeros (or a specified padding value) outside the original array bounds. Note that this implementation is not currently very fast.

```julia
A[pad(i,3)]
```

```julia
pad=NaN
```

--------------------------------

### Logical AND Reduction with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Perform a logical AND reduction over boolean conditions using the `(&)` operator with `@tullio`.

```julia
# Logical AND
@tullio (&) p := A[i] > 0
# p == true (all positive)
```

--------------------------------

### Tullio Including Other Constants

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Shows the preferred method for including other constants in @tullio expressions, using the `$` prefix. It also notes that no gradient is calculated for such constants.

```julia
A[i] := $d * B[i]
```

--------------------------------

### Set global defaults for @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Configure global default settings for `@tullio` macro options such as threading, vectorization, and gradient computation using keyword arguments directly in the macro call.

```julia
using Tullio

# Set global defaults
@tullio threads=false avx=true grad=Base
```

--------------------------------

### Custom Padding Value with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Specify a custom padding value using the `pad=` keyword argument with `pad(i)`.

```julia
# Pad with specific value
@tullio F[i] := B[pad(i, 3)]  pad=100  # pad with 100 instead of 0
```

--------------------------------

### Tullio Eval for Automatic Differentiation

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Shows how Tullio.jl's `Eval` structure is used to enable automatic differentiation by providing forward and reverse pass functions.

```julia
C = Tullio.Eval(make, ∇make)(A, B)
```

--------------------------------

### Tullio Macro with Function Application

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Demonstrates a Tullio.jl macro call that includes a function application (tanh) after the core tensor operation.

```julia
@tullio C[i,j] := tanh <| A[i,k] * B[k,j]
```

--------------------------------

### Generic In-place Update with ^= and @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Use the generic in-place update operator `^=` with a reduction function like `(max)` for flexible updates.

```julia
# Generic in-place update with ^=
C = copy(A)
@tullio (max) C[i] ^= 5i
# C == max.(5:5:50, A)
```

--------------------------------

### Multi-dimensional Reduction with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Perform reductions across multiple dimensions using `@tullio`, as shown here for finding the maximum over `j` and `k`.

```julia
# Multi-dimensional max
Q = rand(1:1000, 4, 5, 6)
@tullio (max) R[i] := Q[i, j, k]  # max over j, k
# R == vec(maximum(Q, dims=(2,3)))
```

--------------------------------

### Clamped Indexing for Edge Repetition

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Employ `clamp(i)` in `@tullio` to repeat edge values for clamped boundary conditions.

```julia
# Clamped boundary - repeats edge values
@tullio D[i] := B[clamp(i)]  (i in 1:10)  avx=false
# D == [1, 2, 3, 4, 5, 5, 5, 5, 5, 5]
```

--------------------------------

### Einsum.jl compatibility with @einsum

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Use the `@einsum` macro for compatibility with Einsum.jl syntax. This macro automatically disables threading, AVX, and gradients for performance.

```julia
using Tullio: @einsum

A = rand(3, 3)
B = rand(3, 3)

# Example usage of @einsum would go here, but is not provided in the source.
```

--------------------------------

### Downsampling with Tullio.jl

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Illustrates downsampling an array by averaging adjacent elements using Tullio's index notation.

```julia
using Tullio
A = [abs2(i - 11) for i in 1:21]

# Downsample -- range of i is that allowed by both terms: 
@tullio B[i] := (A[2i] + A[2i+1])/2  # 1:10 == intersect(1:10, 0:10)
```

--------------------------------

### Stencil Operation with Tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Implements a stencil operation using the @tullio macro. It defines offsets and uses `clamp` to keep indices within bounds, demonstrating how to access existing arrays for computation.

```julia
offsets = [(a,b) for a in -2:2 for b in -2:2 if a>=b] # vector of tuples

@tullio out[x,y,1] = begin
        a,b = offsets[k]
        i = clamp(x+a, extrema(axes(mat,1))...)
        # j = clamp(y+b, extrema(axes(mat,2))...) # can be written clamp(y+b)
        @inbounds mat[i, clamp(y+b), 1] * 10
    end # ranges of x,y read from out[x,y,1]
```

--------------------------------

### Apply functions after reduction with pipe operators

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Use `|>` or `<|` to apply functions like `sqrt`, `log`, or `tanh` after a reduction operation. This is useful for computing norms or numerically stable operations like log-sum-exp. The underscore `_` can be used for anonymous functions.

```julia
using Tullio

A = [i^2 for i in 1:10]
M = rand(10, 10)

# Column norms - sqrt applied after sum
@tullio N[j] := sqrt <| M[i, j]^2
# Equivalent to: N ≈ map(norm, eachcol(M))

# Alternative syntax with |>
@tullio N2[j] := M[i, j]^2 |> sqrt
# Same as above

# Log-sum-exp (numerically stable softmax helper)
mat = randn(5, 10)
@tullio lse[j] := log <| exp(mat[i, j])
# lse ≈ vec(log.(sum(exp.(mat), dims=1)))

# Using underscore for anonymous function
@tullio n3[_] := A[i]^3 |> (_)^(1/3)
# n3[1] ≈ norm(A, 3)

# Finaliser with array indexing (non-reduced index)
B = rand(10, 10)
@tullio (max) E[i] := float(B[i, j]) |> atan(_, A[i])
# E ≈ vec(atan.(maximum(B, dims=2), A))

# Tanh activation after matrix multiply
W = rand(3, 3)
x = rand(3, 4)
@tullio y[i, k] := W[i, j] * x[j, k] |> tanh  avx=false
```

--------------------------------

### Logical OR Reduction with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Execute a logical OR reduction over boolean conditions using the `(|)` operator with `@tullio`.

```julia
# Logical OR
@tullio (|) q := A[i] > 50  avx=false
# q == true (some > 50)
```

--------------------------------

### Index Shifts and Stencils with @tullio

Source: https://context7.com/mcabbott/tullio.jl/llms.txt

Implement operations involving shifted indices, such as convolutions or stencils. Tullio automatically infers index ranges.

```julia
using Tullio, OffsetArrays

A = [i^2 for i in 1:21]

# Simple index shift - range automatically inferred
@tullio B[i] := A[2i+1] + A[i]  # i in 1:10
```

```julia
# Convolution with kernel
K = OffsetArray([1, -1, 2, -1, 1], -2:2)
@tullio C[i] := A[i+j] * K[j]  # j in -2:2 implies i in 3:19
```

```julia
# Stencil operation with explicit ranges
@tullio M[i, j] := A[i+j-1]  (j in 1:15)  # i automatically computed
```

```julia
# Downsampling
@tullio D[i] := (A[2i] + A[2i+1]) / 2  # i in 1:10
```

```julia
# 2D convolution
N = rand(100, 100, 3)
kernel = randn(7, 7)
@tullio out[x, y, c] := N[x+i-1, y+j-1, c] * kernel[i, j]
```

```julia
# Magic shift to force 1-based output
@tullio Z[i+_] := A[2i+10]  # axes(Z,1) === Base.OneTo(5)
```

--------------------------------

### Tullio Index Range Specification

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Shows how to explicitly specify index ranges when they cannot be inferred by the @tullio macro, using the syntax `A[i] := ... (i in range)`.

```julia
A[i] := i^2  (i in 1:10)
```

--------------------------------

### Applying a Vector of Functions with Tullio

Source: https://github.com/mcabbott/tullio.jl/blob/master/README.md

Applies a vector of functions to elements of a matrix using the @tullio macro. This showcases how to use functions stored in an array within a Tullio expression.

```julia
fs = [sin, cos, tan]
xs = randn(3,100)
@tullio ys[r,c] := (fs[r])(xs[r,c])
```