Add comprehensive documentation for planar tensor contractions (@planar and @plansor)

docs/src/lib/tensors.md

@@ -212,6 +212,17 @@ contract!
 ⊗(::AbstractTensorMap, ::AbstractTensorMap)
 ```
 
+### Planar tensor contractions
+
+For fermionic and anyonic tensors, planar contractions are needed to correctly handle
+non-trivial braiding statistics. See the manual section on
+[Fermionic tensor contractions](@ref) and [Anyonic tensor contractions](@ref) for
+detailed information and examples.
+
+The `@planar` macro performs planar tensor contractions that respect non-trivial braiding,
+while `@plansor` automatically dispatches between standard and planar contractions based
+on the `BraidingStyle` of the tensors.
+
 ## `TensorMap` factorizations
 
 The factorisation methods are powered by [MatrixAlgebraKit.jl](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl)

docs/src/man/tensors.md

@@ -1265,8 +1265,158 @@ Some examples:
 
 ## Fermionic tensor contractions
 
-TODO
+When working with fermionic tensors, i.e., tensors with `BraidingStyle(sectortype(tensor)) isa Fermionic`,
+the non-trivial braiding statistics of fermionic particles must be accounted for in tensor
+contractions. Unlike bosonic tensors where the order of operations doesn't matter (up to 
+trivial permutations), fermionic tensors require careful tracking of how tensor factors are
+reordered and how indices are permuted.
+
+### The `@planar` macro
+
+For fermionic tensors, use the `@planar` macro instead of `@tensor`. The `@planar` macro
+correctly handles the signs that arise from reordering fermionic indices:
+
+```julia
+# Example with fermionic tensors
+V = Vect[FermionParity](0 => 2, 1 => 2)  # Fermionic vector space
+A = randn(V ← V)
+B = randn(V ← V')
+C = randn(V' ← V)
+
+# Use @planar for fermionic contractions
+@planar D[i; j] := A[i; j] + B[i; k] * C[k; j]
+```
+
+The key difference from `@tensor` is that `@planar` uses planar operations (`planaradd!`,
+`planarcontract!`, `planartrace!`) that respect the fermionic braiding when permuting indices.
+These operations automatically insert the appropriate signs (``-1`` raised to the appropriate
+power) when indices with odd fermion parity are exchanged.
+
+### Braiding tensors
+
+The braiding of fermionic indices can be made explicit using braiding tensors, denoted by `τ`.
+A braiding tensor `τ[a b; c d]` represents the operation that permutes indices, inserting
+the appropriate fermionic signs. The braiding tensor must always have exactly 2 input and 2 output indices.
+
+The braiding tensor `τ` is special: TensorKit.jl automatically constructs it based on the
+spaces of the indices it connects. You don't need to define it yourself; simply use the
+symbol `τ` in your `@planar` expressions. For examples of using braiding tensors, see the
+test suite in `test/planar.jl`.
+
+### Key differences from bosonic contractions
+
+1. **Order matters**: In `@tensor`, the order of factors doesn't affect the result (modulo
+   permutations). In `@planar`, reordering factors can introduce minus signs.
+
+2. **Index permutations**: When indices are permuted, `@planar` automatically accounts for
+   the fermionic statistics, inserting signs when odd-parity indices are exchanged.
+
+3. **Performance**: `@planar` has some overhead compared to `@tensor` due to the additional
+   bookkeeping required for braiding. For bosonic tensors, stick to `@tensor`.
+
+For more details on the mathematical framework underlying fermionic tensor contractions,
+see the paper [arXiv:2404.14611](https://arxiv.org/abs/2404.14611), which describes the
+implementation of fermionic tensor networks in TensorKit.jl.
 
 ## Anyonic tensor contractions
 
-TODO
+Anyonic tensors arise in the study of quantum systems with anyonic excitations, such as
+fractional quantum Hall systems or topological quantum computing. In these systems, the
+braiding of anyonic particles is described by a non-trivial braiding (or R-matrix), and
+the tensor contractions must respect this braiding structure.
+
+### Using `@planar` for anyonic tensors
+
+Just like fermionic tensors, anyonic tensors require the use of `@planar` instead of `@tensor`:
+
+```julia
+# Example with Fibonacci anyons
+V = Vect[FibonacciAnyon](:I => 2, :τ => 2)
+A = randn(V ← V)
+B = randn(V ← V')
+C = randn(V' ← V)
+
+# Use @planar for anyonic contractions
+@planar D[i; j] := A[i; j] + B[i; k] * C[k; j]
+```
+
+For anyonic systems, the braiding is even more general than for fermions: instead of just
+signs, the braiding tensors can have non-trivial matrix elements that depend on the
+fusion-splitting structure of the anyons. The `@planar` macro handles all of this
+automatically through the braiding tensor `τ`.
+
+### Anyonic braiding tensors
+
+The braiding tensors for anyons encode the non-trivial R-matrices (or F-matrices for the
+fusion structure) of the anyonic theory. Similar to the fermionic case, these can be used
+explicitly:
+
+```julia
+# Explicit braiding in anyonic contractions
+V = Vect[FibonacciAnyon](:I => 2, :τ => 2)
+A = randn(V ← V ⊗ V)
+B = randn(V ⊗ V ← V)
+# Using braiding tensor to permute indices
+@planar C[i; j] := A[i; k l] * τ[k l; m n] * B[m n; j]
+```
+
+The matrix elements of `τ` are determined by the specific anyonic theory (e.g., Fibonacci,
+Ising, etc.) and are computed automatically by TensorKit.jl based on the sector type.
+
+### Differences from bosonic and fermionic contractions
+
+1. **General braiding**: Anyonic braiding is the most general case, with fermionic being
+   a special case (ℤ₂-graded braiding) and bosonic being the trivial case.
+
+2. **Fusion rules**: Anyonic tensors have more complex fusion rules than fermions or bosons,
+   which affects how tensor contractions are computed.
+
+3. **Computational cost**: Anyonic tensor contractions can be more expensive than fermionic
+   or bosonic contractions due to the more complex braiding structure.
+
+## Using `@plansor` for generic code
+
+When writing generic code that should work for bosonic, fermionic, and anyonic tensors,
+use the `@plansor` macro. This macro automatically dispatches to the appropriate
+implementation based on the `BraidingStyle` of the tensors:
+
+```julia
+function generic_contraction(A::AbstractTensorMap, B::AbstractTensorMap)
+    @plansor C[i; j] := A[i; k] * B[k; j]
+    return C
+end
+
+# Works correctly for all three cases:
+# 1. Bosonic tensors: dispatches to @tensor (fast)
+V_bosonic = ℂ^2
+A_bosonic = randn(V_bosonic ← V_bosonic)
+B_bosonic = randn(V_bosonic ← V_bosonic)
+C_bosonic = generic_contraction(A_bosonic, B_bosonic)
+
+# 2. Fermionic tensors: dispatches to @planar (correct signs)
+V_fermionic = Vect[FermionParity](0 => 2, 1 => 2)
+A_fermionic = randn(V_fermionic ← V_fermionic)
+B_fermionic = randn(V_fermionic ← V_fermionic)
+C_fermionic = generic_contraction(A_fermionic, B_fermionic)
+
+# 3. Anyonic tensors: dispatches to @planar (correct braiding)
+V_anyonic = Vect[FibonacciAnyon](:I => 2, :τ => 2)
+A_anyonic = randn(V_anyonic ← V_anyonic)
+B_anyonic = randn(V_anyonic ← V_anyonic)
+C_anyonic = generic_contraction(A_anyonic, B_anyonic)
+```
+
+The `@plansor` macro incurs a small runtime overhead for bosonic tensors due to the dispatch,
+but ensures correctness for all braiding types.
+
+## Summary: `@tensor` vs `@planar` vs `@plansor`
+
+| Macro | Use case | Braiding handling | Performance |
+|-------|----------|-------------------|-------------|
+| `@tensor` | Bosonic tensors (no symmetries or abelian symmetries) | Ignores braiding (assumes trivial) | Fastest |
+| `@planar` | Fermionic or anyonic tensors | Correctly handles non-trivial braiding | Slower due to braiding |
+| `@plansor` | Generic code (unknown braiding type) | Runtime dispatch based on `BraidingStyle` | Small dispatch overhead |
+
+**Recommendation**: Use `@tensor` for bosonic tensors in performance-critical code,
+`@planar` for fermionic/anyonic tensors, and `@plansor` for generic library code or when
+the braiding type is not known at compile time.

src/planar/macros.jl

@@ -1,4 +1,71 @@
 # NEW MACROS: @planar and @plansor
+"""
+    @planar [kwargs...] expr
+
+Perform a planar tensor contraction, i.e., a tensor contraction that respects the 
+non-trivial braiding of fermionic or anyonic tensor indices. This macro is the 
+planar counterpart of the `@tensor` macro from TensorOperations.jl.
+
+The `@planar` macro should be used when working with tensors that have fermionic or 
+anyonic symmetries (i.e., `BraidingStyle(sectortype(tensor)) isa Fermionic` or 
+`BraidingStyle(sectortype(tensor)) isa Anyonic`). For such tensors, the order of 
+tensor contractions and index permutations matters due to non-trivial braiding, and 
+the `@planar` macro correctly handles the braiding tensors `τ` that encode these 
+braidings.
+
+# Syntax
+
+The syntax is identical to `@tensor`, with the addition that braiding tensors can be 
+explicitly specified using `τ[i j; k l]` to denote a braiding between indices.
+
+```julia
+@planar C[i; j] := A[i; k] * B[k; j]
+@planar C[i; j] := A[i; k] * τ[k l; m n] * B[m; n] * τ[n o; j l]
+```
+
+# Keyword arguments
+
+The same keyword arguments as `@tensor` are supported:
+- `opt=true` or `opt=optexpr`: optimize the contraction order
+- `backend=backend`: specify the backend for the contraction
+- `allocator=allocator`: specify the allocator for temporary tensors
+- `order=(...)`: manually specify the contraction order
+- `contractcheck=true`: add runtime checks for contraction validity
+
+# Key differences from `@tensor`
+
+While `@tensor` assumes trivial braiding (appropriate for bosonic tensors), `@planar` 
+correctly handles non-trivial braiding statistics. Specifically:
+
+1. `@tensor` ignores the order of tensor factors in a product (they commute)
+2. `@planar` respects the order and handles explicit braiding tensors `τ`
+3. `@planar` uses planar operations (`planaradd!`, `planarcontract!`, etc.) that 
+   account for braiding when permuting indices
+
+For bosonic tensors (trivial braiding), `@planar` reduces to the same operations as 
+`@tensor`, but with potential performance overhead. Use `@plansor` to automatically 
+dispatch to the appropriate implementation based on the braiding style.
+
+# Example
+
+```julia
+# Fermionic tensors
+V = Vect[FermionParity](0 => 2, 1 => 2)
+A = randn(V ← V)
+B = randn(V ← V')
+C = randn(V' ← V)
+
+# Correct handling of fermionic statistics
+@planar D[i; j] := A[i; j] + B[i; k] * C[k; j]
+
+# With explicit braiding tensor
+A2 = randn(V ← V ⊗ V)
+B2 = randn(V ⊗ V ← V)
+@planar C2[i; j] := A2[i; k l] * τ[k l; m n] * B2[m n; j]
+```
+
+See also: [`@plansor`](@ref), [`@tensor`](@ref)
+"""
 macro planar(args::Vararg{Expr})
     isempty(args) && throw(ArgumentError("No arguments passed to `planar`"))
 
@@ -91,6 +158,76 @@ function planarparser(planarexpr, kwargs...)
     return parser
 end
 
+"""
+    @plansor [kwargs...] expr
+
+Automatically dispatch between `@tensor` and `@planar` based on the braiding style 
+of the tensors involved in the contraction.
+
+This macro provides a unified interface for tensor contractions that works correctly 
+for both bosonic tensors (trivial braiding) and fermionic/anyonic tensors (non-trivial 
+braiding). At runtime, it checks `BraidingStyle(sectortype(tensor))` and:
+
+- If `Bosonic`: uses the standard `@tensor` implementation (efficient, ignores braiding)
+- If `Fermionic` or `Anyonic`: uses the `@planar` implementation (handles braiding correctly)
+
+# Syntax
+
+The syntax is the same as `@tensor` and `@planar`:
+
+```julia
+@plansor C[i; j] := A[i; k] * B[k; j]
+@plansor C[i; j] := A[i; k] * τ[k l; m n] * B[m; n] * τ[n o; j l]
+```
+
+# Keyword arguments
+
+The same keyword arguments as `@tensor` and `@planar` are supported:
+- `opt=true` or `opt=optexpr`: optimize the contraction order
+- `backend=backend`: specify the backend for the contraction
+- `allocator=allocator`: specify the allocator for temporary tensors
+- `order=(...)`: manually specify the contraction order
+- `contractcheck=true`: add runtime checks for contraction validity
+
+# When to use `@plansor` vs `@tensor` vs `@planar`
+
+- Use `@tensor` when you know your tensors have trivial braiding (bosonic, or no symmetries)
+  and want maximum performance
+- Use `@planar` when you know your tensors have non-trivial braiding (fermionic, anyonic)
+  and need correct handling of braiding
+- Use `@plansor` when writing generic code that should work for both cases, or when you're 
+  unsure about the braiding properties of your tensors
+
+# Performance considerations
+
+For bosonic tensors, `@plansor` has a small runtime overhead compared to `@tensor` due to 
+the runtime dispatch. For performance-critical code with known bosonic tensors, prefer 
+`@tensor`. For fermionic/anyonic tensors, `@plansor` and `@planar` have identical performance.
+
+# Example
+
+```julia
+# Works correctly for both bosonic and fermionic tensors
+function my_contraction(A::AbstractTensorMap, B::AbstractTensorMap)
+    @plansor C[i; j] := A[i; k] * B[k; j]
+    return C
+end
+
+# Bosonic case (uses @tensor internally)
+V_bosonic = ℂ^2
+A_bosonic = randn(V_bosonic ← V_bosonic)
+B_bosonic = randn(V_bosonic ← V_bosonic)
+C_bosonic = my_contraction(A_bosonic, B_bosonic)
+
+# Fermionic case (uses @planar internally)
+V_fermionic = Vect[FermionParity](0 => 2, 1 => 2)
+A_fermionic = randn(V_fermionic ← V_fermionic)
+B_fermionic = randn(V_fermionic ← V_fermionic)
+C_fermionic = my_contraction(A_fermionic, B_fermionic)
+```
+
+See also: [`@planar`](@ref), [`@tensor`](@ref)
+"""
 macro plansor(args::Vararg{Expr})
     isempty(args) && throw(ArgumentError("No arguments passed to `planar`"))