diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 72f40322..b1483f5f 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -8,6 +8,7 @@ on: - 'release-' tags: '*' pull_request: + workflow_dispatch: jobs: build: diff --git a/docs/make.jl b/docs/make.jl index d1cec878..5238a42a 100644 --- a/docs/make.jl +++ b/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) diff --git a/docs/src/lib/sectors.md b/docs/src/lib/sectors.md index 4f45b387..a117f171 100644 --- a/docs/src/lib/sectors.md +++ b/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 diff --git a/docs/src/lib/spaces.md b/docs/src/lib/spaces.md index 87781715..99246af8 100644 --- a/docs/src/lib/spaces.md +++ b/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) diff --git a/docs/src/lib/tensors.md b/docs/src/lib/tensors.md index c460ebb9..fe89f491 100644 --- a/docs/src/lib/tensors.md +++ b/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 diff --git a/docs/src/man/tensors.md b/docs/src/man/tensors.md index 897a7711..4d3e3d4f 100644 --- a/docs/src/man/tensors.md +++ b/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. diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index 2bbf6064..568ba95a 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/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} diff --git a/src/fusiontrees/iterator.jl b/src/fusiontrees/iterator.jl index e14ef562..92e8865a 100644 --- a/src/fusiontrees/iterator.jl +++ b/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) diff --git a/src/sectors/anyons.jl b/src/sectors/anyons.jl index 168179e0..f32adf32 100644 --- a/src/sectors/anyons.jl +++ b/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 diff --git a/src/sectors/fermions.jl b/src/sectors/fermions.jl index a28bee74..a1d4b522 100644 --- a/src/sectors/fermions.jl +++ b/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" diff --git a/src/sectors/irreps.jl b/src/sectors/irreps.jl index f55b0e8b..0eb3ad03 100644 --- a/src/sectors/irreps.jl +++ b/src/sectors/irreps.jl @@ -23,7 +23,8 @@ struct IrrepTable end const Irrep A constant of a singleton type used as `Irrep[G]` with `G<:Group` a type of group, to -construct or obtain a concrete subtype of `AbstractIrrep{G}` that implements the data structure used to represent irreducible representations of the group `G`. +construct or obtain a concrete subtype of `AbstractIrrep{G}` that implements the data +structure used to represent irreducible representations of the group `G`. """ const Irrep = IrrepTable() @@ -66,14 +67,17 @@ end # ZNIrrep: irreps of Z_N are labelled by integers mod N; do we ever want N > 64? """ + struct ZNIrrep{N} <: AbstractIrrep{ℤ{N}} ZNIrrep{N}(n::Integer) Irrep[ℤ{N}](n::Integer) -Represents irreps of the group ``ℤ_N`` for some value of `N<64`. (We need 2*(N-1) <= 127 -in order for a ⊗ b to work correctly.) For `N` equals `2`, `3` or `4`, `ℤ{N}` can be -replaced by `ℤ₂`, `ℤ₃`, `ℤ₄`, whereas `Parity` is a synonym for `Irrep{ℤ₂}`. An arbitrary -`Integer` `n` can be provided to the constructor, but only the value `mod(n, N)` is -relevant. +Represents irreps of the group ``ℤ_N`` for some value of `N<64`. (We need 2*(N-1) <= 127 in +order for a ⊗ b to work correctly.) For `N` equals `2`, `3` or `4`, `ℤ{N}` can be replaced +by `ℤ₂`, `ℤ₃`, `ℤ₄`. An arbitrary `Integer` `n` can be provided to the constructor, but only +the value `mod(n, N)` is relevant. + +## Fields +- `n::Int8`: the integer label of the irrep, modulo `N`. """ struct ZNIrrep{N} <: AbstractIrrep{ℤ{N}} n::Int8 @@ -107,14 +111,18 @@ Base.isless(c1::ZNIrrep{N}, c2::ZNIrrep{N}) where {N} = isless(c1.n, c2.n) # U1Irrep: irreps of U1 are labelled by integers """ - U1Irrep(j::Real) - Irrep[U₁](j::Real) + struct U1Irrep <: AbstractIrrep{U₁} + U1Irrep(charge::Real) + Irrep[U₁](charge::Real) Represents irreps of the group ``U₁``. The irrep is labelled by a charge, which should be an integer for a linear representation. However, it is often useful to allow half integers to represent irreps of ``U₁`` subgroups of ``SU₂``, such as the Sz of spin-1/2 system. Hence, the charge is stored as a `HalfInt` from the package HalfIntegers.jl, but can be entered as arbitrary `Real`. The sequence of the charges is: 0, 1/2, -1/2, 1, -1, ... + +## Fields +- `charge::HalfInt`: the label of the irrep, which can be any half integer. """ struct U1Irrep <: AbstractIrrep{U₁} charge::HalfInt @@ -153,12 +161,16 @@ function Base.show(io::IO, ::SU2IrrepException) end """ + struct SU2Irrep <: AbstractIrrep{SU₂} SU2Irrep(j::Real) Irrep[SU₂](j::Real) Represents irreps of the group ``SU₂``. The irrep is labelled by a half integer `j` which can be entered as an abitrary `Real`, but is stored as a `HalfInt` from the HalfIntegers.jl package. + +## Fields +- `j::HalfInt`: the label of the irrep, which can be any non-negative half integer. """ struct SU2Irrep <: AbstractIrrep{SU₂} j::HalfInt @@ -219,13 +231,18 @@ Base.isless(s1::SU2Irrep, s2::SU2Irrep) = isless(s1.j, s2.j) # U₁ ⋊ C (U₁ and charge conjugation) """ + struct CU1Irrep <: AbstractIrrep{CU₁} CU1Irrep(j, s = ifelse(j>zero(j), 2, 0)) Irrep[CU₁](j, s = ifelse(j>zero(j), 2, 0)) Represents irreps of the group ``U₁ ⋊ C`` (``U₁`` and charge conjugation or reflection), -which is also known as just `O₂`. The irrep is labelled by a positive half integer `j` (the -``U₁`` charge) and an integer `s` indicating the behaviour under charge conjugation. They -take values: +which is also known as just `O₂`. + +## Fields +- `j::HalfInt`: the value of the ``U₁`` charge. +- `s::Int`: the representation of charge conjugation. + +They can take values: * if `j == 0`, `s = 0` (trivial charge conjugation) or `s = 1` (non-trivial charge conjugation) * if `j > 0`, `s = 2` (two-dimensional representation) diff --git a/src/spaces/vectorspaces.jl b/src/spaces/vectorspaces.jl index 14595580..d8e764f0 100644 --- a/src/spaces/vectorspaces.jl +++ b/src/spaces/vectorspaces.jl @@ -58,12 +58,11 @@ Return the vector space associated to object `a`. """ function space end -""" +@doc """ dim(V::VectorSpace) -> Int Return the total dimension of the vector space `V` as an Int. -""" -function dim end +""" dim(::VectorSpace) """ dual(V::VectorSpace) -> VectorSpace @@ -106,6 +105,12 @@ const IndexSpace = ElementarySpace field(V::ElementarySpace) = field(typeof(V)) # field(::Type{<:ElementarySpace{𝕜}}) where {𝕜} = 𝕜 +@doc """ + dim(V::ElementarySpace, s::Sector) -> Int + +Return the degeneracy dimension corresponding to the sector `s` of the vector space `V`. +""" dim(::ElementarySpace, ::Sector) + """ oneunit(V::S) where {S<:ElementarySpace} -> S diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl index 604608ae..55195e4f 100644 --- a/src/tensors/abstracttensor.jl +++ b/src/tensors/abstracttensor.jl @@ -185,15 +185,16 @@ function adjointtensorindices(t::AbstractTensorMap, p::Index2Tuple) end @doc """ - blocks(t::AbstractTensorMap) + blocks(t::AbstractTensorMap) -> SectorDict{<:Sector,<:DenseMatrix} -Return an iterator over all blocks of a tensor, i.e. all coupled sectors and their corresponding blocks. +Return an iterator over all blocks of a tensor, i.e. all coupled sectors and their +corresponding blocks. See also [`block`](@ref), [`blocksectors`](@ref), [`blockdim`](@ref) and [`hasblock`](@ref). """ blocks @doc """ - block(t::AbstractTensorMap, c::Sector) + block(t::AbstractTensorMap, c::Sector) -> DenseMatrix Return the block of a tensor corresponding to a coupled sector `c`. @@ -223,7 +224,7 @@ Return the dimensions of the block of a tensor corresponding to a coupled sector fusiontrees(t::AbstractTensorMap) Return an iterator over all splitting - fusion tree pairs of a tensor. -""" fusiontrees +""" fusiontrees(::AbstractTensorMap) # Equality and approximality #---------------------------- diff --git a/src/tensors/tensor.jl b/src/tensors/tensor.jl index 2ec6a3b5..e93d792a 100644 --- a/src/tensors/tensor.jl +++ b/src/tensors/tensor.jl @@ -114,6 +114,26 @@ dim(t::TensorMap) = mapreduce(x -> length(x[2]), +, blocks(t); init=0) # General TensorMap constructors #-------------------------------- # constructor starting from block data +""" + TensorMap(data::AbstractDict{<:Sector,<:DenseMatrix}, codomain::ProductSpace{S,N₁}, + domain::ProductSpace{S,N₂}) where {S<:ElementarySpace,N₁,N₂} + TensorMap(data, codomain ← domain) + TensorMap(data, domain → codomain) + +Construct a `TensorMap` by explicitly specifying its block data. + +## Arguments +- `data::AbstractDict{<:Sector,<:DenseMatrix}`: dictionary containing the block data for + each coupled sector `c` as a `DenseMatrix` of size + `(blockdim(codomain, c), blockdim(domain, c))`. +- `codomain::ProductSpace{S,N₁}`: the codomain as a `ProductSpace` of `N₁` spaces of type + `S<:ElementarySpace`. +- `domain::ProductSpace{S,N₂}`: the domain as a `ProductSpace` of `N₂` spaces of type + `S<:ElementarySpace`. + +Alternatively, the domain and codomain can be specified by passing a [`HomSpace`](@ref) +using the syntax `codomain ← domain` or `domain → codomain`. +""" function TensorMap(data::AbstractDict{<:Sector,<:DenseMatrix}, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}) where {S<:IndexSpace,N₁,N₂} I = sectortype(S) @@ -199,6 +219,32 @@ function _buildblockstructure(P::ProductSpace{S,N}, blocksectors) where {S<:Inde return treeranges, blockdims end +""" + TensorMap([f, eltype,] codomain::ProductSpace{S,N₁}, domain::ProductSpace{S,N₂}) + where {S<:ElementarySpace,N₁,N₂} + TensorMap([f, eltype,], codomain ← domain) + TensorMap([f, eltype,], domain → codomain) + +Construct a `TensorMap` from a general callable that produces block data for each coupled +sector. + +## Arguments +- `f`: callable object that returns a `DenseMatrix`, or `UndefInitializer`. +- `eltype::Type{<:Number}`: element type of the data. +- `codomain::ProductSpace{S,N₁}`: the codomain as a `ProductSpace` of `N₁` spaces of type + `S<:ElementarySpace`. +- `domain::ProductSpace{S,N₂}`: the domain as a `ProductSpace` of `N₂` spaces of type + `S<:ElementarySpace`. + +If `eltype` is left unspecified, `f` should support the calling syntax `f(::Tuple{Int,Int})` +such that `f((m, n))` returns a `DenseMatrix` with `size(f((m, n))) == (m, n)`. If `eltype` is +specified, `f` is instead called as `f(eltype, (m, n))`. In the case where `f` is left +unspecified or `undef` is passed explicitly, a `TensorMap` with uninitialized data is +generated. + +Alternatively, the domain and codomain can be specified by passing a [`HomSpace`](@ref) +using the syntax `codomain ← domain` or `domain → codomain`. +""" function TensorMap(f, ::Type{T}, codom::ProductSpace{S}, dom::ProductSpace{S}) where {S<:IndexSpace,T<:Number} return TensorMap(d -> f(T, d), codom, dom) @@ -263,6 +309,40 @@ Tensor(T::Type{<:Number}, P::TensorSpace{S}) where {S<:IndexSpace} = TensorMap(T Tensor(P::TensorSpace{S}) where {S<:IndexSpace} = TensorMap(P, one(P)) # constructor starting from a dense array +""" + TensorMap(data::DenseArray, codomain::ProductSpace{S,N₁}, domain::ProductSpace{S,N₂}; + tol=sqrt(eps(real(float(eltype(data)))))) where {S<:ElementarySpace,N₁,N₂} + TensorMap(data, codomain ← domain) + TensorMap(data, domain → codomain) + +Construct a `TensorMap` from a plain multidimensional array. + +## Arguments +- `data::DenseArray`: tensor data as a plain array. +- `codomain::ProductSpace{S,N₁}`: the codomain as a `ProductSpace` of `N₁` spaces of type + `S<:ElementarySpace`. +- `domain::ProductSpace{S,N₂}`: the domain as a `ProductSpace` of `N₂` spaces of type + `S<:ElementarySpace`. +- `tol=sqrt(eps(real(float(eltype(data)))))::Float64`: + +Here, `data` can be specified in two ways. It can either be a `DenseArray` of rank `N₁ + N₂` +whose size matches that of the domain and codomain spaces, +`size(data) == (dims(codomain)..., dims(domain)...)`, or a `DenseMatrix` where +`size(data) == (dim(codomain), dim(domain))`. The `TensorMap` constructor will then +reconstruct the tensor data such that the resulting tensor `t` satisfies +`data == convert(Array, t)`. For the case where `sectortype(S) == Trivial`, the `data` array +is simply reshaped into matrix form and referred to as such in the resulting `TensorMap` +instance. When `S<:GradedSpace`, the `data` array has to be compatible with how how each +sector in every space `V` is assigned to an index range within `1:dim(V)`. + +Alternatively, the domain and codomain can be specified by passing a [`HomSpace`](@ref) +using the syntax `codomain ← domain` or `domain → codomain`. + +!!! note + This constructor only works for `sectortype` values for which conversion to a plain + array is possible, and only in the case where the `data` actually respects the specified + symmetry structure. +""" function TensorMap(data::DenseArray, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}; tol=sqrt(eps(real(float(eltype(data)))))) where {S<:IndexSpace,N₁,N₂} (d1, d2) = (dim(codom), dim(dom)) @@ -466,6 +546,21 @@ blocks(t::TensorMap) = t.data fusiontrees(t::TrivialTensorMap) = ((nothing, nothing),) fusiontrees(t::TensorMap) = TensorKeyIterator(t.rowr, t.colr) +""" + Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I}, + sectors::NTuple{N₁+N₂,I}) where {N₁,N₂,I<:Sector} + -> StridedViews.StridedView + t[sectors] + +Return a view into the data slice of `t` corresponding to the splitting - fusion tree pair +with combined uncoupled charges `sectors`. In particular, if `sectors == (s1..., s2...)` +where `s1` and `s2` correspond to the coupled charges in the codomain and domain +respectively, then a `StridedViews.StridedView` of size +`(dims(codomain(t), s1)..., dims(domain(t), s2))` is returned. + +This method is only available for the case where `FusionStyle(I) isa UniqueFusion`, +since it assumes a uniquely defined coupled charge. +""" @inline function Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I}, sectors::Tuple{Vararg{I}}) where {N₁,N₂,I<:Sector} FusionStyle(I) isa UniqueFusion || @@ -488,6 +583,20 @@ end return t[map(sectortype(t), sectors)] end +""" + Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I}, + f₁::FusionTree{I,N₁}, + f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + -> StridedViews.StridedView + t[f₁, f₂] + +Return a view into the data slice of `t` corresponding to the splitting - fusion tree pair +`(f₁, f₂)`. In particular, if `f₁.coupled == f₂.coupled == c`, then a +`StridedViews.StridedView` of size +`(dims(codomain(t), f₁.uncoupled)..., dims(domain(t), f₂.uncoupled))` is returned which +represents the slice of `block(t, c)` whose row indices correspond to `f₁.uncoupled` and +column indices correspond to `f₂.uncoupled`. +""" @inline function Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I}, f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} @@ -502,6 +611,21 @@ end return sreshape(StridedView(t.data[c])[t.rowr[c][f₁], t.colr[c][f₂]], d) end end + +""" + Base.setindex!(t::TensorMap{<:IndexSpace,N₁,N₂,I}, + v, + f₁::FusionTree{I,N₁}, + f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + t[f₁, f₂] = v + +Copies `v` into the data slice of `t` corresponding to the splitting - fusion tree pair +`(f₁, f₂)`. Here, `v` can be any object that can be copied into a `StridedViews.StridedView` +of size `(dims(codomain(t), f₁.uncoupled)..., dims(domain(t), f₂.uncoupled))` using +`Base.copy!`. + +See also [`Base.getindex(::TensorMap{<:IndexSpace,N₁,N₂,I<:Sector}, ::FusionTree{I<:Sector,N₁}, ::FusionTree{I<:Sector,N₂})`](@ref) +""" @propagate_inbounds function Base.setindex!(t::TensorMap{<:IndexSpace,N₁,N₂,I}, v, f₁::FusionTree{I,N₁}, @@ -510,6 +634,13 @@ end end # For a tensor with trivial symmetry, allow no argument indexing +""" + Base.getindex(t::TrivialTensorMap) + t[] + +Return a view into the data of `t` as a `StridedViews.StridedView` of size +`(dims(codomain(t))..., dims(domain(t))...)`. +""" @inline function Base.getindex(t::TrivialTensorMap) return sreshape(StridedView(t.data), (dims(codomain(t))..., dims(domain(t))...)) end @@ -520,12 +651,25 @@ end @inline Base.setindex!(t::TrivialTensorMap, v, ::Nothing, ::Nothing) = setindex!(t, v) # For a tensor with trivial symmetry, allow direct indexing +""" + Base.getindex(t::TrivialTensorMap, indices::Vararg{Int}) + t[indices] + +Return a view into the data slice of `t` corresponding to `indices`, by slicing the +`StridedViews.StridedView` into the full data array. +""" @inline function Base.getindex(t::TrivialTensorMap, indices::Vararg{Int}) data = t[] @boundscheck checkbounds(data, indices...) @inbounds v = data[indices...] return v end +""" + Base.setindex!(t::TrivialTensorMap, v, indices::Vararg{Int}) + t[indices] = v + +Assigns `v` to the data slice of `t` corresponding to `indices`. +""" @inline function Base.setindex!(t::TrivialTensorMap, v, indices::Vararg{Int}) data = t[] @boundscheck checkbounds(data, indices...)