Docs additions related to constructing a given symmetric TensorMap

.github/workflows/Documentation.yml

@@ -8,6 +8,7 @@ on:
       - 'release-'
     tags: '*'
   pull_request:
+  workflow_dispatch:
 
 jobs:
   build:

docs/make.jl

@@ -12,4 +12,4 @@ makedocs(; modules=[TensorKit],
                 "Library" => ["lib/sectors.md", "lib/spaces.md", "lib/tensors.md"],
                 "Index" => ["index/index.md"]])
 
-deploydocs(; repo="github.com/Jutho/TensorKit.jl.git")
+deploydocs(; repo="github.com/Jutho/TensorKit.jl.git", push_preview=true)

docs/src/lib/sectors.md

@@ -18,7 +18,10 @@ SU2Irrep
 CU1Irrep
 ProductSector
 FermionParity
+FermionNumber
+FermionSpin
 FibonacciAnyon
+IsingAnyon
 FusionTree
 ```
 
@@ -42,6 +45,8 @@ Base.isreal(::Type{<:Sector})
 TensorKit.vertex_labeltype
 TensorKit.vertex_ind2label
 ⊠(::Sector, ::Sector)
+fusiontrees(uncoupled::NTuple{N,I}, coupled::I,
+                     isdual::NTuple{N,Bool}) where {N,I<:Sector}
 ```
 
 ## Methods for manipulating fusion trees

docs/src/lib/spaces.md

@@ -16,6 +16,7 @@ ComplexSpace
 GradedSpace
 CompositeSpace
 ProductSpace
+HomSpace
 ```
 
 ## Useful constants
@@ -38,7 +39,9 @@ field
 sectortype
 sectors
 hassector
-dim
+dim(::VectorSpace)
+dim(::ElementarySpace, ::Sector)
+dim(P::ProductSpace{<:ElementarySpace,N}, sector::NTuple{N,<:Sector}) where {N}
 dims
 blocksectors(::ProductSpace)
 blocksectors(::HomSpace)

docs/src/lib/tensors.md

@@ -24,23 +24,38 @@ TrivialTensorMap
 TrivialTensor
 ```
 
-## Specific `TensorMap` constructors
+## `TensorMap` constructors
 
+### General constructors
+
+A general `TensorMap` can be constructed by specifying its data, codmain and domain in one
+of the following ways:
 ```@docs
-id
-isomorphism
-unitary
-isometry
+TensorMap(::AbstractDict{<:Sector,<:DenseMatrix}, ::ProductSpace{S,N₁},
+                   ::ProductSpace{S,N₂}) where {S<:IndexSpace,N₁,N₂}
+TensorMap(::Any, ::Type{T}, codom::ProductSpace{S},
+                   dom::ProductSpace{S}) where {S<:IndexSpace,T<:Number}
+TensorMap(::DenseArray, ::ProductSpace{S,N₁}, ::ProductSpace{S,N₂};
+                   tol) where {S<:IndexSpace,N₁,N₂}
 ```
 
-Additionally, several special-purpose methods exist to generate data according to specific distributions:
-
+Several special-purpose methods exist to generate data according to specific distributions:
 ```@docs
 randuniform
 randnormal
 randisometry
 ```
 
+### Specific constructors
+
+Additionally, several special-purpose constructors exist to generate data according to specific distributions:
+```@docs
+id
+isomorphism
+unitary
+isometry
+```
+
 ## Accessing properties and data
 
 The following methods exist to obtain type information:
@@ -71,10 +86,28 @@ blocksectors(::AbstractTensorMap)
 blockdim(::AbstractTensorMap, ::Sector)
 block
 blocks
-fusiontrees
+fusiontrees(::AbstractTensorMap)
 hasblock
 ```
 
+For `TensorMap`s with `Trivial` `sectortype`, the data can be directly accessed and
+manipulated in a straightforward way:
+```@docs
+Base.getindex(t::TrivialTensorMap)
+Base.getindex(t::TrivialTensorMap, indices::Vararg{Int})
+Base.setindex!(t::TrivialTensorMap, ::Any, indices::Vararg{Int})
+```
+
+For general `TensorMap`s, this can be done using custom `getindex` and `setindex!` methods:
+```@docs
+Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I},
+              sectors::Tuple{Vararg{I}}) where {N₁,N₂,I<:Sector}
+Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I},
+              f₁::FusionTree{I,N₁},
+              f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector}
+Base.setindex!(::TensorMap{<:IndexSpace,N₁,N₂,I}, ::Any, ::FusionTree{I,N₁}, ::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector}
+```
+
 ## `TensorMap` operations
 
 The operations that can be performed on a `TensorMap` can be organized into *index

docs/src/man/tensors.md

@@ -224,7 +224,7 @@ via `t[f₁,f₂]`, and is returned as a `StridedArray` of size
 implementation does not distinguish between `FusionStyle isa UniqueFusion` or
 `FusionStyle isa MultipleFusion`, in the former case the fusion tree is completely
 characterized by the uncoupled sectors, and so the subblocks can also be accessed as
-`t[(a1, …, aN₁), (b1, …, bN₂)]`. When there is no symmetry at all, i.e.
+`t[(a1, …, aN₁, b1, …, bN₂)]`. When there is no symmetry at all, i.e.
 `sectortype(t) == Trivial`, `t[]` returns the raw tensor data as a `StridedArray` of size
 `(dim(V1), …, dim(VN₁), dim(W1), …, dim(WN₂))`, whereas `block(t, Trivial())` returns the
 same data as a `DenseMatrix` of size `(dim(V1) * … * dim(VN₁), dim(W1) * … * dim(WN₂))`.
@@ -313,13 +313,28 @@ e.g. for `CartesianSpace` or `ComplexSpace`. Then the `data` array is just resha
 matrix form and referred to as such in the resulting `TensorMap` instance. When `spacetype`
 is `GradedSpace`, the `TensorMap` constructor will try to reconstruct the tensor data such
 that the resulting tensor `t` satisfies `data == convert(Array, t)`. This might not be
-possible, if the data does not respect the symmetry structure. Let's sketch this with a
-simple example
+possible, if the data does not respect the symmetry structure. This procedure can be
+sketched using a simple physical example, namely the SWAP gate on two qubits,
+```math
+\begin{align*}
+\mathrm{SWAP}: \mathbb{C}^2 \otimes \mathbb{C}^2 & \to \mathbb{C}^2 \otimes \mathbb{C}^2\\
+|i\rangle \otimes |j\rangle &\mapsto |j\rangle \otimes |i\rangle.
+\end{align*}
+```
+This operator can be rewritten in terms of the familiar Heisenberg exchange interaction
+``\vec{S}_i \cdot \vec{S}_j`` as
+```math
+\mathrm{SWAP} = 2 \vec{S}_i \cdot \vec{S}_j + \frac{1}{2} 𝟙,
+```
+where ``\vec{S} = (S^x, S^y, S^z)`` and the spin-1/2 generators of SU₂ ``S^k`` are defined
+defined in terms of the ``2 \times 2`` Pauli matrices ``\sigma^k`` as
+``S^k = \frac{1}{2}\sigma^k``. The SWAP gate can be realized as a rank-4 `TensorMap` in the
+following way:
 ```@repl tensors
+# encode the matrix elements of the swap gate into a rank-4 array, where the first two
+# indices correspond to the codomain and the last two indices correspond to the domain
 data = zeros(2,2,2,2)
-# encode the operator (σ_x * σ_x + σ_y * σ_y + σ_z * σ_z)/2
-# that is, the swap gate, which maps the last two indices on the first two in reversed order
-# also known as Heisenberg interaction between two spin 1/2 particles
+# the swap gate then maps the last two indices on the first two in reversed order
 data[1,1,1,1] = data[2,2,2,2] = data[1,2,2,1] = data[2,1,1,2] = 1
 V1 = ℂ^2 # generic qubit hilbert space
 t1 = TensorMap(data, V1 ⊗ V1, V1 ⊗ V1)
@@ -329,12 +344,12 @@ V3 = U1Space(1/2=>1,-1/2=>1) # restricted space that only uses the `σ_z` rotati
 t3 = TensorMap(data, V3 ⊗ V3, V3 ⊗ V3)
 for (c,b) in blocks(t3)
     println("Data for block $c :")
-    b |> disp
+    disp(b)
     println()
 end
 ```
-Hence, we recognize that the Heisenberg interaction has eigenvalue ``-1`` in the coupled
-spin zero sector (`SUIrrep(0)`), and eigenvalue ``+1`` in the coupled spin 1 sector
+Hence, we recognize that the exchange interaction has eigenvalue ``-1`` in the coupled spin
+zero sector (`SU2Irrep(0)`), and eigenvalue ``+1`` in the coupled spin 1 sector
 (`SU2Irrep(1)`). Using `Irrep[U₁]` instead, we observe that both coupled charge
 `U1Irrep(+1)` and `U1Irrep(-1)` have eigenvalue ``+1``. The coupled charge `U1Irrep(0)`
 sector is two-dimensional, and has an eigenvalue ``+1`` and an eigenvalue ``-1``.
@@ -357,6 +372,68 @@ axes(P, (SU2Irrep(1), SU2Irrep(0), SU2Irrep(2)))
 Note that the length of the range is the degeneracy dimension of that sector, times the
 dimension of the internal representation space, i.e. the quantum dimension of that sector.
 
+### Assigning block data after initialization
+
+In order to avoid having to know the internal structure of each representation space to
+properly construct the full `data` array, it is often simpler to assign the block data
+directly after initializing an all zero `TensorMap` with the correct spaces. While this may
+seem more difficult at first sight since it requires knowing the exact entries associated to
+each valid combination of domain uncoupled sectors, coupled sector and codomain uncoupled
+sectors, this is often a far more natural procedure in practice.
+
+A first option is to directly set the full matrix block for each coupled sector in the
+`TensorMap`. For the example with U₁ symmetry, this can be done as
+```@repl tensors
+t4 = TensorMap(zeros, V3 ⊗ V3, V3 ⊗ V3);
+block(t4, U1Irrep(0)) .= [1 0; 0 1];
+block(t4, U1Irrep(1)) .= [1;;];
+block(t4, U1Irrep(-1)) .= [1;;];
+for (c, b) in blocks(t4)
+    println("Data for block $c :")
+    disp(b)
+    println()
+end
+```
+While this indeed does not require considering the internal structure of the representation
+spaces, it still requires knowing the precise row and column indices corresponding to each
+set of uncoupled sectors in the codmain and domain respectively to correctly assign the
+nonzero entries in each block.
+
+Perhaps the most natural way of constructing a particular `TensorMap` is to directly assign
+the data slices for each splitting - fusion tree pair using the `fusiontrees(::TensorMap)`
+method. This returns an iterator over all tuples `(f₁, f₂)` of splitting - fusion tree pairs
+corresponding to all ways in which the set of domain uncoupled sectors can fuse to a coupled
+sector and split back into the set of codomain uncoupled sectors. By directly setting the
+corresponding data slice `t[f₁, f₂]` of size
+`(dims(codomain(t), f₁.uncoupled)..., dims(domain(t), f₂.uncoupled)...)`, we can construct
+all the block data without worrying about the internal ordering of row and column indices in
+each block. In addition, the corresponding value of each fusion tree slice is often directly
+informed by the object we are trying to construct in the first place. For example, in order
+to construct the Heisenberg exchange interaction on two spin-1/2 particles ``i`` and ``j``
+as an SU₂ symmetric `TensorMap`, we can make use of the observation that
+```math
+\vec{S}_i \cdot \vec{S}_j = \frac{1}{2} \left( \left( \vec{S}_i \cdot \vec{S}_j \right)^2 - \vec{S}_i^2 - \vec{S}_j^2 \right).
+```
+Recalling some basic group theory, we know that the
+[quadratic Casimir of SU₂](https://en.wikipedia.org/wiki/Representation_theory_of_SU(2)#The_Casimir_element),
+``\vec{S}^2``, has a well-defined eigenvalue ``j(j+1)`` on every irrep of spin ``j``. From
+the above expressions, we can therefore directly read off the eigenvalues of the SWAP gate
+in terms of this Casimir eigenvalue on the domain uncoupled sectors and the coupled sector.
+This gives us exactly the prescription we need to assign the data slice corresponding to
+each splitting - fusion tree pair:
+```@repl tensors
+C(s::SU2Irrep) = s.j * (s.j + 1)
+t5 = TensorMap(zeros, V2 ⊗ V2, V2 ⊗ V2);
+for (f₁, f₂) in fusiontrees(t5)
+    t5[f₁, f₂] .= C(f₂.coupled) - C(f₂.uncoupled[1]) - C(f₂.uncoupled[2]) + 1/2
+end
+for (c, b) in blocks(t5)
+    println("Data for block $c :")
+    disp(b)
+    println()
+end
+```
+
 ### Constructing similar tensors
 
 A third way to construct a `TensorMap` instance is to use `Base.similar`, i.e.

src/fusiontrees/fusiontrees.jl

@@ -4,13 +4,19 @@
     struct FusionTree{I, N, M, L, T}
 
 Represents a fusion tree of sectors of type `I<:Sector`, fusing (or splitting) `N` uncoupled
-sectors to a coupled sector. (It actually represents a splitting tree, but fusion tree
-is a more common term). The `isdual` field indicates whether an isomorphism is present
-(if the corresponding value is true) or not. The field `uncoupled` contains the sectors
-coming out of the splitting trees, before the possible 𝑍 isomorphism. This fusion tree
-has `M=max(0, N-2)` inner lines. Furthermore, for `FusionStyle(I) isa GenericFusion`,
-the `L=max(0, N-1)` corresponding vertices carry a label of type `T`. If `FusionStyle(I)
-isa MultiplicityFreeFusion, `T = Nothing`.
+sectors to a coupled sector. It actually represents a splitting tree, but fusion tree
+is a more common term.
+
+## Fields
+- `uncoupled::NTuple{N,I}`: the uncoupled sectors coming out of the splitting tree, before
+  the possible 𝑍 isomorphism (see `isdual`).
+- `coupled::I`: the coupled sector.
+- `isdual::NTuple{N,Bool}`: indicates whether a 𝑍 isomorphism is present (`true`) or not
+  (`false`) for each uncoupled sector.
+- `innerlines::NTuple{M,I}`: the labels of the M=max(0, N-2)` inner lines of the splitting
+  tree.
+- `vertices::NTuple{L,T}`: the `L=max(0, N-1)` labels of type `T` of the vertices of the
+  splitting tree. If `FusionStyle(I) isa MultiplicityFreeFusion`, then `T = Nothing`.
 """
 struct FusionTree{I<:Sector,N,M,L,T}
     uncoupled::NTuple{N,I}

src/fusiontrees/iterator.jl

@@ -1,6 +1,15 @@
 # FusionTreeIterator:
 # iterate over fusion trees for fixed coupled and uncoupled sector labels
 #==============================================================================#
+
+"""
+    fusiontrees(uncoupled::NTuple{N,I}[,
+        coupled::I=one(I)[, isdual::NTuple{N,Bool}=ntuple(n -> false, length(uncoupled))]])
+        where {N,I<:Sector} -> FusionTreeIterator{I,N}
+
+Return an iterator over all fusion trees with a given coupled sector label `coupled` and
+uncoupled sector labels and isomorphisms `uncoupled` and `isdual` respectively.
+"""
 function fusiontrees(uncoupled::NTuple{N,I}, coupled::I,
                      isdual::NTuple{N,Bool}) where {N,I<:Sector}
     return FusionTreeIterator{I,N}(uncoupled, coupled, isdual)

src/sectors/anyons.jl

@@ -35,10 +35,13 @@ Fsymbol(::Vararg{PlanarTrivial,6}) = 1
     struct FibonacciAnyon <: Sector
     FibonacciAnyon(s::Symbol)
 
-Represents the anyons (isomorphism classes of simple objects) of the Fibonacci fusion
-category. It can take two values, corresponding to the trivial sector
-`FibonacciAnyon(:I)` and the non-trivial sector `FibonacciAnyon(:τ)` with fusion rules
-``τ ⊗ τ = 1 ⊕ τ``.
+Represents the anyons of the Fibonacci modular fusion category. It can take two values,
+corresponding to the trivial sector `FibonacciAnyon(:I)` and the non-trivial sector
+`FibonacciAnyon(:τ)` with fusion rules ``τ ⊗ τ = 1 ⊕ τ``.
+
+## Fields
+- `isone::Bool`: indicates whether the sector corresponds the to trivial anyon `:I`
+  (`true`), or the non-trivial anyon `:τ` (`false`).
 """
 struct FibonacciAnyon <: Sector
     isone::Bool
@@ -149,10 +152,13 @@ Base.isless(a::FibonacciAnyon, b::FibonacciAnyon) = isless(!a.isone, !b.isone)
     struct IsingAnyon <: Sector
     IsingAnyon(s::Symbol)
 
-Represents the anyons (isomorphism classes of simple objects) of the Ising fusion category.
-It can take three values, corresponding to the trivial sector `IsingAnyon(:I)` and the
-non-trivial sectors `IsingAnyon(:σ)` and `IsingAnyon(:ψ)`, with fusion rules
-``ψ ⊗ ψ = 1``, ``σ ⊗ ψ = σ``, and ``σ ⊗ σ = 1 ⊕ ψ``.
+Represents the anyons of the Ising modular fusion category. It can take three values,
+corresponding to the trivial sector `IsingAnyon(:I)` and the non-trivial sectors
+`IsingAnyon(:σ)` and `IsingAnyon(:ψ)`, with fusion rules ``ψ ⊗ ψ = 1``, ``σ ⊗ ψ = σ``, and
+``σ ⊗ σ = 1 ⊕ ψ``.
+
+## Fields
+- `s::Symbol`: the label of the represented anyon, which can be `:I`, `:σ`, or `:ψ`.
 """
 struct IsingAnyon <: Sector
     s::Symbol

src/sectors/fermions.jl

@@ -4,7 +4,10 @@
 Represents sectors with fermion parity. The fermion parity is a ℤ₂ quantum number that
 yields an additional sign when two odd fermions are exchanged.
 
-See also: `FermionNumber`, `FermionSpin`
+## Fields
+- `isodd::Bool`: indicates whether the fermion parity is odd (`true`) or even (`false`).
+
+See also: [`FermionNumber`](@ref), [`FermionSpin`](@ref)
 """
 struct FermionParity <: Sector
     isodd::Bool
@@ -66,6 +69,15 @@ Base.isless(a::FermionParity, b::FermionParity) = isless(a.isodd, b.isodd)
 # Common fermionic combinations
 # -----------------------------
 
+"""
+    const FermionNumber = U1Irrep ⊠ FermionParity
+    FermionNumber(a::Int)
+
+Represents the fermion number as the direct product of a ``U₁`` irrep `a` and a fermion
+parity, with the restriction that the fermion parity is odd if and only if `a` is odd.
+
+See also: [`U1Irrep`](@ref), [`FermionParity`](@ref)
+"""
 const FermionNumber = U1Irrep ⊠ FermionParity
 const fU₁ = FermionNumber
 FermionNumber(a::Int) = U1Irrep(a) ⊠ FermionParity(isodd(a))
@@ -74,9 +86,18 @@ type_repr(::Type{FermionNumber}) = "FermionNumber"
 # convenience default converter -> allows Vect[FermionNumber](1 => 1)
 Base.convert(::Type{FermionNumber}, a::Int) = FermionNumber(a)
 
+"""
+    const FermionSpin = SU2Irrep ⊠ FermionParity
+    FermionSpin(j::Real)
+
+Represents the fermion spin as the direct product of a ``SU₂`` irrep `j` and a fermion
+parity, with the restriction that the fermion parity is odd if `2 * j` is odd.
+
+See also: [`SU2Irrep`](@ref), [`FermionParity`](@ref)
+"""
 const FermionSpin = SU2Irrep ⊠ FermionParity
 const fSU₂ = FermionSpin
-FermionSpin(a::Real) = (s = SU2Irrep(a);
+FermionSpin(j::Real) = (s = SU2Irrep(j);
                         s ⊠ FermionParity(isodd(twice(s.j))))
 type_repr(::Type{FermionSpin}) = "FermionSpin"