diff --git a/docs/make.jl b/docs/make.jl index 057981577..55234c2c9 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,7 @@ using Documenter using Random using TensorKit +using TensorKit: FusionTreePair, Index2Tuple using TensorKit.TensorKitSectors using TensorKit.MatrixAlgebraKit using DocumenterInterLinks diff --git a/docs/src/man/sectors.md b/docs/src/man/sectors.md index 73c1630a0..b2fea3a6c 100644 --- a/docs/src/man/sectors.md +++ b/docs/src/man/sectors.md @@ -1159,7 +1159,7 @@ the splitting tree. The `FusionTree` interface to duality and line bending is given by -[`repartition(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, N::Int)`](@ref repartition) +[`repartition(f1::FusionTreePair{I,N₁,N₂}, N::Int)`](@ref repartition) which takes a splitting tree `f1` with `N₁` outgoing sectors, a fusion tree `f2` with `N₂` incoming sectors, and applies line bending such that the resulting splitting and fusion @@ -1184,7 +1184,7 @@ With this basic function, we can now perform arbitrary combinations of braids or permutations with line bendings, to completely reshuffle where sectors appear. The interface provided for this is given by -[`braid(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, levels1::NTuple{N₁,Int}, levels2::NTuple{N₂,Int}, p1::NTuple{N₁′,Int}, p2::NTuple{N₂′,Int})`](@ref braid(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple, ::IndexTuple, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂}) +[`braid((f₁, f₂)::FusionTreePair, (p1, p2)::Index2Tuple, (levels1, levels2)::Index2Tuple)`](@ref braid(::TensorKit.FusionTreePair, ::Index2Tuple, ::Index2Tuple)) where we now have splitting tree `f1` with `N₁` outgoing sectors, a fusion tree `f2` with `N₂` incoming sectors, `levels1` and `levels2` assign a level or depth to the corresponding @@ -1211,7 +1211,7 @@ As before, there is a simplified interface for the case where `BraidingStyle(I) isa SymmetricBraiding` and the levels are not needed. This is simply given by -[`permute(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, p1::NTuple{N₁′,Int}, p2::NTuple{N₂′,Int})`](@ref permute(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂}) +[`permute((f₁, f₂)::FusionTreePair, (p1, p2)::Index2Tuple)`](@ref permute(::FusionTreePair, ::Index2Tuple)) The `braid` and `permute` routines for double fusion trees will be the main access point for corresponding manipulations on tensors. As a consequence, results from this routine are diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 74d23c8ce..1b62a95ab 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -119,7 +119,7 @@ using ScopedValues using TensorKitSectors import TensorKitSectors: dim, BraidingStyle, FusionStyle, ⊠, ⊗ -import TensorKitSectors: dual, type_repr +import TensorKitSectors: dual, type_repr, fusiontensor import TensorKitSectors: twist using Base: @boundscheck, @propagate_inbounds, @constprop, @@ -210,6 +210,19 @@ function set_num_transformer_threads(n::Int) return TRANSFORMER_THREADS[] = n end +const TREEMANIPULATION_THREADS = Ref(1) + +get_num_manipulation_threads() = TREEMANIPULATION_THREADS[] + +function set_num_manipulation_threads(n::Int) + N = Base.Threads.nthreads() + if n > N + n = N + Strided._set_num_threads_warn(n) + end + return TREEMANIPULATION_THREADS[] = n +end + # Definitions and methods for tensors #------------------------------------- # general definitions diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 4b97603ca..23436bf1d 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -27,6 +27,20 @@ function permutation2swaps(perm) return swaps end +# one-argument version: check whether `p` is a cyclic permutation (of `1:length(p)`) +function iscyclicpermutation(p) + N = length(p) + @inbounds for i in 1:N + p[mod1(i + 1, N)] == mod1(p[i] + 1, N) || return false + end + return true +end +# two-argument version: check whether `v1` is a cyclic permutation of `v2` +function iscyclicpermutation(v1, v2) + length(v1) == length(v2) || return false + return iscyclicpermutation(indexin(v1, v2)) +end + _kron(A, B, C, D...) = _kron(_kron(A, B), C, D...) function _kron(A, B) sA = size(A) diff --git a/src/fusiontrees/basic_manipulations.jl b/src/fusiontrees/basic_manipulations.jl new file mode 100644 index 000000000..3f2ca74a7 --- /dev/null +++ b/src/fusiontrees/basic_manipulations.jl @@ -0,0 +1,273 @@ +# BASIC MANIPULATIONS: +#---------------------------------------------- +# -> rewrite generic fusion tree in basis of fusion trees in standard form +# -> only depend on Fsymbol + +""" + insertat(f::FusionTree{I, N₁}, i::Int, f₂::FusionTree{I, N₂}) + -> <:AbstractDict{<:FusionTree{I, N₁+N₂-1}, <:Number} + +Attach a fusion tree `f₂` to the uncoupled leg `i` of the fusion tree `f₁` and bring it +into a linear combination of fusion trees in standard form. This requires that +`f₂.coupled == f₁.uncoupled[i]` and `f₁.isdual[i] == false`. +""" +function insertat(f₁::FusionTree{I}, i::Int, f₂::FusionTree{I, 0}) where {I} + # this actually removes uncoupled line i, which should be trivial + (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) || + throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])")) + coeff = one(sectorscalartype(I)) + + uncoupled = TupleTools.deleteat(f₁.uncoupled, i) + coupled = f₁.coupled + isdual = TupleTools.deleteat(f₁.isdual, i) + if length(uncoupled) <= 2 + inner = () + else + inner = TupleTools.deleteat(f₁.innerlines, max(1, i - 2)) + end + if length(uncoupled) <= 1 + vertices = () + else + vertices = TupleTools.deleteat(f₁.vertices, max(1, i - 1)) + end + f = FusionTree(uncoupled, coupled, isdual, inner, vertices) + return fusiontreedict(I)(f => coeff) +end +function insertat(f₁::FusionTree{I}, i, f₂::FusionTree{I, 1}) where {I} + # identity operation + (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) || + throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])")) + coeff = one(sectorscalartype(I)) + isdual′ = TupleTools.setindex(f₁.isdual, f₂.isdual[1], i) + f = FusionTree{I}(f₁.uncoupled, f₁.coupled, isdual′, f₁.innerlines, f₁.vertices) + return fusiontreedict(I)(f => coeff) +end +function insertat(f₁::FusionTree{I}, i, f₂::FusionTree{I, 2}) where {I} + # elementary building block, + (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) || + throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])")) + uncoupled = f₁.uncoupled + coupled = f₁.coupled + inner = f₁.innerlines + b, c = f₂.uncoupled + isdual = f₁.isdual + isdualb, isdualc = f₂.isdual + if i == 1 + uncoupled′ = (b, c, tail(uncoupled)...) + isdual′ = (isdualb, isdualc, tail(isdual)...) + inner′ = (uncoupled[1], inner...) + vertices′ = (f₂.vertices..., f₁.vertices...) + coeff = one(sectorscalartype(I)) + f′ = FusionTree(uncoupled′, coupled, isdual′, inner′, vertices′) + return fusiontreedict(I)(f′ => coeff) + end + uncoupled′ = TupleTools.insertafter(TupleTools.setindex(uncoupled, b, i), i, (c,)) + isdual′ = TupleTools.insertafter(TupleTools.setindex(isdual, isdualb, i), i, (isdualc,)) + inner_extended = (uncoupled[1], inner..., coupled) + a = inner_extended[i - 1] + d = inner_extended[i] + e′ = uncoupled[i] + if FusionStyle(I) isa MultiplicityFreeFusion + local newtrees + for e in a ⊗ b + coeff = conj(Fsymbol(a, b, c, d, e, e′)) + iszero(coeff) && continue + inner′ = TupleTools.insertafter(inner, i - 2, (e,)) + f′ = FusionTree(uncoupled′, coupled, isdual′, inner′) + if @isdefined newtrees + push!(newtrees, f′ => coeff) + else + newtrees = fusiontreedict(I)(f′ => coeff) + end + end + return newtrees + else + local newtrees + κ = f₂.vertices[1] + λ = f₁.vertices[i - 1] + for e in a ⊗ b + inner′ = TupleTools.insertafter(inner, i - 2, (e,)) + Fmat = Fsymbol(a, b, c, d, e, e′) + for μ in axes(Fmat, 1), ν in axes(Fmat, 2) + coeff = conj(Fmat[μ, ν, κ, λ]) + iszero(coeff) && continue + vertices′ = TupleTools.setindex(f₁.vertices, ν, i - 1) + vertices′ = TupleTools.insertafter(vertices′, i - 2, (μ,)) + f′ = FusionTree(uncoupled′, coupled, isdual′, inner′, vertices′) + if @isdefined newtrees + push!(newtrees, f′ => coeff) + else + newtrees = fusiontreedict(I)(f′ => coeff) + end + end + end + return newtrees + end +end +function insertat(f₁::FusionTree{I, N₁}, i, f₂::FusionTree{I, N₂}) where {I, N₁, N₂} + F = fusiontreetype(I, N₁ + N₂ - 1) + (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) || + throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])")) + T = sectorscalartype(I) + coeff = one(T) + if length(f₁) == 1 + return fusiontreedict(I){F, T}(f₂ => coeff) + end + if i == 1 + uncoupled = (f₂.uncoupled..., tail(f₁.uncoupled)...) + isdual = (f₂.isdual..., tail(f₁.isdual)...) + inner = (f₂.innerlines..., f₂.coupled, f₁.innerlines...) + vertices = (f₂.vertices..., f₁.vertices...) + coupled = f₁.coupled + f′ = FusionTree(uncoupled, coupled, isdual, inner, vertices) + return fusiontreedict(I){F, T}(f′ => coeff) + else # recursive definition + N2 = length(f₂) + f₂′, f₂′′ = split(f₂, N2 - 1) + local newtrees::fusiontreedict(I){F, T} + for (f, coeff) in insertat(f₁, i, f₂′′) + for (f′, coeff′) in insertat(f, i, f₂′) + if @isdefined newtrees + coeff′′ = coeff * coeff′ + newtrees[f′] = get(newtrees, f′, zero(coeff′′)) + coeff′′ + else + newtrees = fusiontreedict(I){F, T}(f′ => coeff * coeff′) + end + end + end + return newtrees + end +end + +""" + split(f::FusionTree{I, N}, M::Int) + -> (::FusionTree{I, M}, ::FusionTree{I, N-M+1}) + +Split a fusion tree into two. The first tree has as uncoupled sectors the first `M` +uncoupled sectors of the input tree `f`, whereas its coupled sector corresponds to the +internal sector between uncoupled sectors `M` and `M+1` of the original tree `f`. The +second tree has as first uncoupled sector that same internal sector of `f`, followed by +remaining `N-M` uncoupled sectors of `f`. It couples to the same sector as `f`. This +operation is the inverse of `insertat` in the sense that if +`f₁, f₂ = split(t, M) ⇒ f == insertat(f₂, 1, f₁)`. +""" +@inline function split(f::FusionTree{I, N}, M::Int) where {I, N} + if M > N || M < 0 + throw(ArgumentError("M should be between 0 and N = $N")) + elseif M === N + (f, FusionTree{I}((f.coupled,), f.coupled, (false,), (), ())) + elseif M === 1 + isdual1 = (f.isdual[1],) + isdual2 = TupleTools.setindex(f.isdual, false, 1) + f₁ = FusionTree{I}((f.uncoupled[1],), f.uncoupled[1], isdual1, (), ()) + f₂ = FusionTree{I}(f.uncoupled, f.coupled, isdual2, f.innerlines, f.vertices) + return f₁, f₂ + elseif M === 0 + u = leftunit(f.uncoupled[1]) + f₁ = FusionTree{I}((), u, (), ()) + uncoupled2 = (u, f.uncoupled...) + coupled2 = f.coupled + isdual2 = (false, f.isdual...) + innerlines2 = N >= 2 ? (f.uncoupled[1], f.innerlines...) : () + if FusionStyle(I) isa GenericFusion + vertices2 = (1, f.vertices...) + return f₁, FusionTree{I}(uncoupled2, coupled2, isdual2, innerlines2, vertices2) + else + return f₁, FusionTree{I}(uncoupled2, coupled2, isdual2, innerlines2) + end + else + uncoupled1 = ntuple(n -> f.uncoupled[n], M) + isdual1 = ntuple(n -> f.isdual[n], M) + innerlines1 = ntuple(n -> f.innerlines[n], max(0, M - 2)) + coupled1 = f.innerlines[M - 1] + vertices1 = ntuple(n -> f.vertices[n], M - 1) + + uncoupled2 = ntuple(N - M + 1) do n + return n == 1 ? f.innerlines[M - 1] : f.uncoupled[M + n - 1] + end + isdual2 = ntuple(N - M + 1) do n + return n == 1 ? false : f.isdual[M + n - 1] + end + innerlines2 = ntuple(n -> f.innerlines[M - 1 + n], N - M - 1) + coupled2 = f.coupled + vertices2 = ntuple(n -> f.vertices[M - 1 + n], N - M) + + f₁ = FusionTree{I}(uncoupled1, coupled1, isdual1, innerlines1, vertices1) + f₂ = FusionTree{I}(uncoupled2, coupled2, isdual2, innerlines2, vertices2) + return f₁, f₂ + end +end + +""" + merge(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, c::I, μ = 1) + -> <:AbstractDict{<:FusionTree{I, N₁+N₂}, <:Number} + +Merge two fusion trees together to a linear combination of fusion trees whose uncoupled +sectors are those of `f₁` followed by those of `f₂`, and where the two coupled sectors of +`f₁` and `f₂` are further fused to `c`. In case of +`FusionStyle(I) == GenericFusion()`, also a degeneracy label `μ` for the fusion of +the coupled sectors of `f₁` and `f₂` to `c` needs to be specified. +""" +function merge(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, c::I) where {I, N₁, N₂} + if FusionStyle(I) isa GenericFusion + throw(ArgumentError("vertex label for merging required")) + end + return merge(f₁, f₂, c, 1) +end +function merge(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, c::I, μ) where {I, N₁, N₂} + if !(c in f₁.coupled ⊗ f₂.coupled) + throw(SectorMismatch("cannot fuse sectors $(f₁.coupled) and $(f₂.coupled) to $c")) + end + if μ > Nsymbol(f₁.coupled, f₂.coupled, c) + throw(ArgumentError("invalid fusion vertex label $μ")) + end + f₀ = FusionTree{I}((f₁.coupled, f₂.coupled), c, (false, false), (), (μ,)) + f, coeff = first(insertat(f₀, 1, f₁)) # takes fast path, single output + @assert coeff == one(coeff) + return insertat(f, N₁ + 1, f₂) +end +function merge(f₁::FusionTree{I, 0}, f₂::FusionTree{I, 0}, c::I, μ) where {I} + Nsymbol(f₁.coupled, f₂.coupled, c) == μ == 1 || + throw(SectorMismatch("cannot fuse sectors $(f₁.coupled) and $(f₂.coupled) to $c")) + return fusiontreedict(I)(f₁ => Fsymbol(c, c, c, c, c, c)[1, 1, 1, 1]) +end + +# flip a duality flag of a fusion tree +function flip((f₁, f₂)::FusionTreePair{I, N₁, N₂}, i::Int; inv::Bool = false) where {I, N₁, N₂} + @assert 0 < i ≤ N₁ + N₂ + if i ≤ N₁ + a = f₁.uncoupled[i] + χₐ = frobenius_schur_phase(a) + θₐ = twist(a) + if !inv + factor = f₁.isdual[i] ? χₐ * θₐ : one(θₐ) + else + factor = f₁.isdual[i] ? one(θₐ) : χₐ * conj(θₐ) + end + isdual′ = TupleTools.setindex(f₁.isdual, !f₁.isdual[i], i) + f₁′ = FusionTree{I}(f₁.uncoupled, f₁.coupled, isdual′, f₁.innerlines, f₁.vertices) + return SingletonDict((f₁′, f₂) => factor) + else + i -= N₁ + a = f₂.uncoupled[i] + χₐ = frobenius_schur_phase(a) + θₐ = twist(a) + if !inv + factor = f₂.isdual[i] ? χₐ * one(θₐ) : θₐ + else + factor = f₂.isdual[i] ? conj(θₐ) : χₐ * one(θₐ) + end + isdual′ = TupleTools.setindex(f₂.isdual, !f₂.isdual[i], i) + f₂′ = FusionTree{I}(f₂.uncoupled, f₂.coupled, isdual′, f₂.innerlines, f₂.vertices) + return SingletonDict((f₁, f₂′) => factor) + end +end +function flip((f₁, f₂)::FusionTreePair{I, N₁, N₂}, ind; inv::Bool = false) where {I, N₁, N₂} + f₁′, f₂′ = f₁, f₂ + factor = one(sectorscalartype(I)) + for i in ind + (f₁′, f₂′), s = only(flip((f₁′, f₂′), i; inv)) + factor *= s + end + return SingletonDict((f₁′, f₂′) => factor) +end diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl new file mode 100644 index 000000000..a75b4e268 --- /dev/null +++ b/src/fusiontrees/braiding_manipulations.jl @@ -0,0 +1,376 @@ +# BRAIDING MANIPULATIONS: +#----------------------------------------------- +# -> manipulations that depend on a braiding +# -> requires both Fsymbol and Rsymbol +""" + artin_braid(f::FusionTree, i; inv::Bool = false) -> <:AbstractDict{typeof(f), <:Number} + +Perform an elementary braid (Artin generator) of neighbouring uncoupled indices `i` and +`i+1` on a fusion tree `f`, and returns the result as a dictionary of output trees and +corresponding coefficients. + +The keyword `inv` determines whether index `i` will braid above or below index `i+1`, i.e. +applying `artin_braid(f′, i; inv = true)` to all the outputs `f′` of +`artin_braid(f, i; inv = false)` and collecting the results should yield a single fusion +tree with non-zero coefficient, namely `f` with coefficient `1`. This keyword has no effect +if `BraidingStyle(sectortype(f)) isa SymmetricBraiding`. +""" +function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N} + 1 <= i < N || throw(ArgumentError(lazy"Cannot swap outputs i=$i and i+1 out of only $N outputs")) + @assert FusionStyle(I) === UniqueFusion() + + uncoupled = f.uncoupled + a, b = uncoupled[i], uncoupled[i + 1] + uncoupled′ = TupleTools.setindex(uncoupled, b, i) + uncoupled′ = TupleTools.setindex(uncoupled′, a, i + 1) + coupled′ = f.coupled + isdual′ = TupleTools.setindex(f.isdual, f.isdual[i], i + 1) + isdual′ = TupleTools.setindex(isdual′, f.isdual[i + 1], i) + inner = f.innerlines + inner_extended = (uncoupled[1], inner..., coupled′) + vertices = f.vertices + oneT = one(sectorscalartype(I)) + + if isunit(a) || isunit(b) + # braiding with trivial sector: simple and always possible + inner′ = inner + vertices′ = vertices + if i > 1 # we also need to alter innerlines and vertices + inner′ = TupleTools.setindex( + inner, + inner_extended[isunit(a) ? (i + 1) : (i - 1)], + i - 1 + ) + vertices′ = TupleTools.setindex(vertices′, vertices[i], i - 1) + vertices′ = TupleTools.setindex(vertices′, vertices[i - 1], i) + end + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′) + return f′ => oneT + end + + BraidingStyle(I) isa NoBraiding && + throw(SectorMismatch("Cannot braid sectors $(uncoupled[i]) and $(uncoupled[i + 1])")) + + if i == 1 + c = N > 2 ? inner[1] : coupled′ + R = oftype(oneT, (inv ? conj(Rsymbol(b, a, c)) : Rsymbol(a, b, c))) + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices) + return f′ => R + end + + # case i > 1: other naming convention + b = uncoupled[i] + d = uncoupled[i + 1] + a = inner_extended[i - 1] + c = inner_extended[i] + e = inner_extended[i + 1] + c′ = first(a ⊗ d) + coeff = oftype( + oneT, + if inv + conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′) + else + Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′)) + end + ) + inner′ = TupleTools.setindex(inner, c′, i - 1) + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′) + return f′ => coeff +end + +function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where {I, N} + 1 <= i < N || throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs")) + + uncoupled = src.uncoupled[1] + a, b = uncoupled[i], uncoupled[i + 1] + uncoupled′ = TupleTools.setindex(uncoupled, b, i) + uncoupled′ = TupleTools.setindex(uncoupled′, a, i + 1) + coupled′ = rightunit(src.uncoupled[1][N]) + + isdual = src.isdual[1] + isdual′ = TupleTools.setindex(isdual, isdual[i], i + 1) + isdual′ = TupleTools.setindex(isdual′, isdual[i + 1], i) + dst = FusionTreeBlock{I}((uncoupled′, ()), (isdual′, ()); sizehint = length(src)) + + oneT = one(sectorscalartype(I)) + + indexmap = treeindex_map(dst) + U = zeros(sectorscalartype(I), length(dst), length(src)) + + if isone(a) || isone(b) # braiding with trivial sector: simple and always possible + for (col, (f, f₂)) in enumerate(fusiontrees(src)) + inner = f.innerlines + inner_extended = (uncoupled[1], inner..., coupled′) + vertices = f.vertices + inner′ = inner + vertices′ = vertices + if i > 1 # we also need to alter innerlines and vertices + inner′ = TupleTools.setindex( + inner, inner_extended[isone(a) ? (i + 1) : (i - 1)], i - 1 + ) + vertices′ = TupleTools.setindex(vertices′, vertices[i], i - 1) + vertices′ = TupleTools.setindex(vertices′, vertices[i - 1], i) + end + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′) + row = indexmap[treeindex_data((f′, f₂))] + @inbounds U[row, col] = oneT + end + return dst => U + end + + BraidingStyle(I) isa NoBraiding && + throw(SectorMismatch(lazy"Cannot braid sectors $a and $b")) + + for (col, (f, f₂)) in enumerate(fusiontrees(src)) + inner = f.innerlines + inner_extended = (uncoupled[1], inner..., coupled′) + vertices = f.vertices + + if i == 1 + c = N > 2 ? inner[1] : coupled′ + if FusionStyle(I) isa MultiplicityFreeFusion + R = oftype(oneT, (inv ? conj(Rsymbol(b, a, c)) : Rsymbol(a, b, c))) + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices) + row = indexmap[treeindex_data((f′, f₂))] + @inbounds U[row, col] = R + else # GenericFusion + μ = vertices[1] + Rmat = inv ? Rsymbol(b, a, c)' : Rsymbol(a, b, c) + for ν in axes(Rmat, 2) + R = oftype(oneT, Rmat[μ, ν]) + iszero(R) && continue + vertices′ = TupleTools.setindex(vertices, ν, 1) + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices′) + row = indexmap[treeindex_data((f′, f₂))] + @inbounds U[row, col] = R + end + end + continue + end + # case i > 1: other naming convention + b = uncoupled[i] + d = uncoupled[i + 1] + a = inner_extended[i - 1] + c = inner_extended[i] + e = inner_extended[i + 1] + if FusionStyle(I) isa UniqueFusion + c′ = first(a ⊗ d) + coeff = oftype( + oneT, + if inv + conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′) + else + Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′)) + end + ) + inner′ = TupleTools.setindex(inner, c′, i - 1) + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′) + row = indexmap[treeindex_data((f′, f₂))] + @inbounds U[row, col] = coeff + elseif FusionStyle(I) isa SimpleFusion + cs = collect(I, intersect(a ⊗ d, e ⊗ conj(b))) + for c′ in cs + coeff = oftype( + oneT, + if inv + conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′) + else + Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′)) + end + ) + iszero(coeff) && continue + inner′ = TupleTools.setindex(inner, c′, i - 1) + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′) + row = indexmap[treeindex_data((f′, f₂))] + @inbounds U[row, col] = coeff + end + else # GenericFusion + cs = collect(I, intersect(a ⊗ d, e ⊗ conj(b))) + for c′ in cs + Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e) + Rmat2 = inv ? Rsymbol(d, a, c′)' : Rsymbol(a, d, c′) + Fmat = Fsymbol(d, a, b, e, c′, c) + μ = vertices[i - 1] + ν = vertices[i] + for σ in 1:Nsymbol(a, d, c′) + for λ in 1:Nsymbol(c′, b, e) + coeff = zero(oneT) + for ρ in 1:Nsymbol(d, c, e), κ in 1:Nsymbol(d, a, c′) + coeff += Rmat1[ν, ρ] * conj(Fmat[κ, λ, μ, ρ]) * + conj(Rmat2[σ, κ]) + end + iszero(coeff) && continue + vertices′ = TupleTools.setindex(vertices, σ, i - 1) + vertices′ = TupleTools.setindex(vertices′, λ, i) + inner′ = TupleTools.setindex(inner, c′, i - 1) + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′) + row = indexmap[treeindex_data((f′, f₂))] + @inbounds U[row, col] = coeff + end + end + end + end + end + + return dst => U +end + +# braid fusion tree +""" + braid(f::FusionTree{<:Sector, N}, p::NTuple{N, Int}, levels::NTuple{N, Int}) + -> <:AbstractDict{typeof(t), <:Number} + +Perform a braiding of the uncoupled indices of the fusion tree `f` and return the result as +a `<:AbstractDict` of output trees and corresponding coefficients. The braiding is +determined by specifying that the new sector at position `k` corresponds to the sector that +was originally at the position `i = p[k]`, and assigning to every index `i` of the original +fusion tree a distinct level or depth `levels[i]`. This permutation is then decomposed into +elementary swaps between neighbouring indices, where the swaps are applied as braids such +that if `i` and `j` cross, ``τ_{i,j}`` is applied if `levels[i] < levels[j]` and +``τ_{j,i}^{-1}`` if `levels[i] > levels[j]`. This does not allow to encode the most general +braid, but a general braid can be obtained by combining such operations. +""" +braid(f::FusionTree{I, N}, p::IndexTuple{N}, levels::IndexTuple{N}) where {I, N} = + braid(f, (p, ()), (levels, ())) +function braid(f::FusionTree{I, N}, (p, _)::Index2Tuple{N, 0}, (levels, _)::Index2Tuple{N, 0}) where {I, N} + TupleTools.isperm(p) || throw(ArgumentError("not a valid permutation: $p")) + @assert FusionStyle(I) isa UniqueFusion + if BraidingStyle(I) isa SymmetricBraiding # this assumes Fsymbols are 1! + coeff = one(sectorscalartype(I)) + for i in 1:N + for j in 1:(i - 1) + if p[j] > p[i] + a, b = f.uncoupled[p[j]], f.uncoupled[p[i]] + coeff *= Rsymbol(a, b, first(a ⊗ b)) + end + end + end + uncoupled′ = TupleTools._permute(f.uncoupled, p) + coupled′ = f.coupled + isdual′ = TupleTools._permute(f.isdual, p) + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′) + return f′ => coeff + else + T = sectorscalartype(I) + c = one(T) + for s in permutation2swaps(p) + inv = levels[s] > levels[s + 1] + f, c′ = artin_braid(f, s; inv) + c *= c′ + l = levels[s] + levels = TupleTools.setindex(levels, levels[s + 1], s) + levels = TupleTools.setindex(levels, l, s + 1) + end + + return f => c + end +end + +# permute fusion tree +""" + permute(f::FusionTree, p::NTuple{N, Int}) -> <:AbstractDict{typeof(t), <:Number} + +Perform a permutation of the uncoupled indices of the fusion tree `f` and returns the result +as a `<:AbstractDict` of output trees and corresponding coefficients; this requires that +`BraidingStyle(sectortype(f)) isa SymmetricBraiding`. +""" +function permute(f::FusionTree{I, N}, p::NTuple{N, Int}) where {I, N} + @assert BraidingStyle(I) isa SymmetricBraiding + return braid(f, p, ntuple(identity, Val(N))) +end + +# braid double fusion tree +""" + braid((f₁, f₂)::FusionTreePair, (p1, p2)::Index2Tuple, (levels1, levels2)::Index2Tuple) + -> <:AbstractDict{<:FusionTreePair{I, N₁, N₂}}, <:Number} + +Input is a fusion-splitting tree pair that describes the fusion of a set of incoming +uncoupled sectors to a set of outgoing uncoupled sectors, represented using the splitting +tree `f₁` and fusion tree `f₂`, such that the incoming sectors `f₂.uncoupled` are fused to +`f₁.coupled == f₂.coupled` and then to the outgoing sectors `f₁.uncoupled`. Compute new +trees and corresponding coefficients obtained from repartitioning and braiding the tree such +that sectors `p1` become outgoing and sectors `p2` become incoming. The uncoupled indices in +splitting tree `f₁` and fusion tree `f₂` have levels (or depths) `levels1` and `levels2` +respectively, which determines how indices braid. In particular, if `i` and `j` cross, +``τ_{i,j}`` is applied if `levels[i] < levels[j]` and ``τ_{j,i}^{-1}`` if `levels[i] > +levels[j]`. This does not allow to encode the most general braid, but a general braid can +be obtained by combining such operations. +""" +function braid(src::Union{FusionTreePair, FusionTreeBlock}, p::Index2Tuple, levels::Index2Tuple) + @assert numind(src) == length(p[1]) + length(p[2]) + @assert numout(src) == length(levels[1]) && numin(src) == length(levels[2]) + @assert TupleTools.isperm((p[1]..., p[2]...)) + return fsbraid((src, p, levels)) +end + +const FSPBraidKey{I, N₁, N₂} = Tuple{FusionTreePair{I}, Index2Tuple{N₁, N₂}, Index2Tuple} +const FSBBraidKey{I, N₁, N₂} = Tuple{FusionTreeBlock{I}, Index2Tuple{N₁, N₂}, Index2Tuple} + +Base.@assume_effects :foldable function _fsdicttype(::Type{T}) where {I, N₁, N₂, T <: FSPBraidKey{I, N₁, N₂}} + E = sectorscalartype(I) + return Pair{fusiontreetype(I, N₁, N₂), E} +end +Base.@assume_effects :foldable function _fsdicttype(::Type{T}) where {I, N₁, N₂, T <: FSBBraidKey{I, N₁, N₂}} + F₁ = fusiontreetype(I, N₁) + F₂ = fusiontreetype(I, N₂) + E = sectorscalartype(I) + return Pair{FusionTreeBlock{I, N₁, N₂, Tuple{F₁, F₂}}, Matrix{E}} +end + +@cached function fsbraid(key::K)::_fsdicttype(K) where {I, N₁, N₂, K <: FSPBraidKey{I, N₁, N₂}} + ((f₁, f₂), (p1, p2), (l1, l2)) = key + p = linearizepermutation(p1, p2, length(f₁), length(f₂)) + levels = (l1..., reverse(l2)...) + (f, f0), coeff1 = repartition((f₁, f₂), N₁ + N₂) + f′, coeff2 = braid(f, p, levels) + (f₁′, f₂′), coeff3 = repartition((f′, f0), N₁) + return (f₁′, f₂′) => coeff1 * coeff2 * coeff3 +end +@cached function fsbraid(key::K)::_fsdicttype(K) where {I, N₁, N₂, K <: FSBBraidKey{I, N₁, N₂}} + src, (p1, p2), (l1, l2) = key + + p = linearizepermutation(p1, p2, numout(src), numin(src)) + levels = (l1..., reverse(l2)...) + + dst, U = repartition(src, numind(src)) + + for s in permutation2swaps(p) + inv = levels[s] > levels[s + 1] + dst, U_tmp = artin_braid(dst, s; inv) + U = U_tmp * U + l = levels[s] + levels = TupleTools.setindex(levels, levels[s + 1], s) + levels = TupleTools.setindex(levels, l, s + 1) + end + + if N₂ == 0 + return dst => U + else + dst, U_tmp = repartition(dst, N₁) + U = U_tmp * U + return dst => U + end +end + +CacheStyle(::typeof(fsbraid), k::FSPBraidKey{I}) where {I} = + FusionStyle(I) isa UniqueFusion ? NoCache() : GlobalLRUCache() +CacheStyle(::typeof(fsbraid), k::FSBBraidKey{I}) where {I} = + FusionStyle(I) isa UniqueFusion ? NoCache() : GlobalLRUCache() + +""" + permute((f₁, f₂)::FusionTreePair, (p1, p2)::Index2Tuple) + -> <:AbstractDict{<:FusionTreePair{I, N₁, N₂}}, <:Number} + +Input is a double fusion tree that describes the fusion of a set of incoming uncoupled +sectors to a set of outgoing uncoupled sectors, represented using the individual trees of +outgoing (`t1`) and incoming sectors (`t2`) respectively (with identical coupled sector +`t1.coupled == t2.coupled`). Computes new trees and corresponding coefficients obtained from +repartitioning and permuting the tree such that sectors `p1` become outgoing and sectors +`p2` become incoming. +""" +function permute(src::Union{FusionTreePair, FusionTreeBlock}, p::Index2Tuple) + @assert BraidingStyle(src) isa SymmetricBraiding + levels1 = ntuple(identity, numout(src)) + levels2 = numout(src) .+ ntuple(identity, numin(src)) + return braid(src, p, (levels1, levels2)) +end diff --git a/src/fusiontrees/duality_manipulations.jl b/src/fusiontrees/duality_manipulations.jl new file mode 100644 index 000000000..7d06a2b54 --- /dev/null +++ b/src/fusiontrees/duality_manipulations.jl @@ -0,0 +1,782 @@ +# ELEMENTARY DUALITY MANIPULATIONS: A- and B-moves +#--------------------------------------------------------- +# -> elementary manipulations that depend on the duality (rigidity) and pivotal structure +# -> planar manipulations that do not require braiding, everything is in Fsymbol (A/Bsymbol) +# -> B-move (bendleft, bendright) is simple in standard basis +# -> A-move (foldleft, foldright) is complicated, needs to be reexpressed in standard form + +@doc """ + bendright((f₁, f₂)::FusionTreePair) -> (f₃, f₄) => coeff + bendright(src::FusionTreeBlock) -> dst => coeffs + +Map the final splitting vertex `a ⊗ b ← c` of `src` to a fusion vertex `a ← c ⊗ dual(b)` in `dst`. +For `FusionStyle(src) === UniqueFusion()`, both `src` and `dst` are simple `FusionTreePair`s, and the +transformation consists of a single coefficient `coeff`. +For generic `FusionStyle`s, the input and output consist of `FusionTreeBlock`s that bundle together +all trees with the same uncoupled charges, and `coeffs` now forms a transformation matrix. + +``` + ╰─┬─╯ | | | ╰─┬─╯ | | | + ╰─┬─╯ | | ╰─┬─╯ | | + ╰ ⋯ ┬╯ | ╰ ⋯ ┬╯ | + | | → ╰─┬─╯ + ╭ ⋯ ┴╮ | ╭ ⋯ ╯ + ╭─┴─╮ | | ╭─┴─╮ + ╭─┴─╮ | ╰─╯ ╭─┴─╮ | +``` + +See also [`bendleft`](@ref). +""" bendright + +function bendright((f₁, f₂)::FusionTreePair) + I = sectortype((f₁, f₂)) + @assert FusionStyle(I) === UniqueFusion() + N₁, N₂ = numout((f₁, f₂)), numin((f₁, f₂)) + @assert N₁ > 0 + c = f₁.coupled + a = N₁ == 1 ? leftunit(f₁.uncoupled[1]) : (N₁ == 2 ? f₁.uncoupled[1] : f₁.innerlines[end]) + b = f₁.uncoupled[N₁] + + # construct the new fusiontree pair + uncoupled1 = TupleTools.front(f₁.uncoupled) + isdual1 = TupleTools.front(f₁.isdual) + inner1 = N₁ > 2 ? TupleTools.front(f₁.innerlines) : () + vertices1 = N₁ > 1 ? TupleTools.front(f₁.vertices) : () + f₁′ = FusionTree{I}(uncoupled1, a, isdual1, inner1, vertices1) + + uncoupled2 = (f₂.uncoupled..., dual(b)) + isdual2 = (f₂.isdual..., !(f₁.isdual[N₁])) + inner2 = N₂ > 1 ? (f₂.innerlines..., c) : () + vertices2 = N₂ > 0 ? (f₂.vertices..., 1) : () + f₂′ = FusionTree{I}(uncoupled2, a, isdual2, inner2, vertices2) + + # compute the coefficient + coeff₀ = sqrtdim(c) * invsqrtdim(a) + f₁.isdual[N₁] && (coeff₀ *= conj(frobenius_schur_phase(dual(b)))) + coeff = coeff₀ * Bsymbol(a, b, c) + + return (f₁′, f₂′) => coeff +end +function bendright(src::FusionTreeBlock) + uncoupled_dst = ( + TupleTools.front(src.uncoupled[1]), + (src.uncoupled[2]..., dual(src.uncoupled[1][end])), + ) + isdual_dst = ( + TupleTools.front(src.isdual[1]), + (src.isdual[2]..., !(src.isdual[1][end])), + ) + I = sectortype(src) + N₁ = numout(src) + N₂ = numin(src) + @assert N₁ > 0 + + dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src)) + indexmap = treeindex_map(dst) + U = zeros(sectorscalartype(I), length(dst), length(src)) + + for (col, (f₁, f₂)) in enumerate(fusiontrees(src)) + c = f₁.coupled + a = N₁ == 1 ? leftunit(f₁.uncoupled[1]) : + (N₁ == 2 ? f₁.uncoupled[1] : f₁.innerlines[end]) + b = f₁.uncoupled[N₁] + + uncoupled1 = TupleTools.front(f₁.uncoupled) + isdual1 = TupleTools.front(f₁.isdual) + inner1 = N₁ > 2 ? TupleTools.front(f₁.innerlines) : () + vertices1 = N₁ > 1 ? TupleTools.front(f₁.vertices) : () + f₁′ = FusionTree(uncoupled1, a, isdual1, inner1, vertices1) + + uncoupled2 = (f₂.uncoupled..., dual(b)) + isdual2 = (f₂.isdual..., !(f₁.isdual[N₁])) + inner2 = N₂ > 1 ? (f₂.innerlines..., c) : () + + coeff₀ = sqrtdim(c) * invsqrtdim(a) + if f₁.isdual[N₁] + coeff₀ *= conj(frobenius_schur_phase(dual(b))) + end + if FusionStyle(I) isa MultiplicityFreeFusion + coeff = coeff₀ * Bsymbol(a, b, c) + vertices2 = N₂ > 0 ? (f₂.vertices..., 1) : () + f₂′ = FusionTree(uncoupled2, a, isdual2, inner2, vertices2) + row = indexmap[treeindex_data((f₁′, f₂′))] + @inbounds U[row, col] = coeff + else + Bmat = Bsymbol(a, b, c) + μ = N₁ > 1 ? f₁.vertices[end] : 1 + for ν in axes(Bmat, 2) + coeff = coeff₀ * Bmat[μ, ν] + iszero(coeff) && continue + vertices2 = N₂ > 0 ? (f₂.vertices..., ν) : () + f₂′ = FusionTree(uncoupled2, a, isdual2, inner2, vertices2) + row = indexmap[treeindex_data((f₁′, f₂′))] + @inbounds U[row, col] = coeff + end + end + end + + return dst => U +end + +@doc """ + bendleft((f₁, f₂)::FusionTreePair) -> (f₃, f₄) => coeff + bendleft(src::FusionTreeBlock) -> dst => coeffs + +Map the final fusion vertex `a ← c ⊗ dual(b)` of `src` to a splitting vertex `a ⊗ b ← c` in `dst`. +For `FusionStyle(src) === UniqueFusion()`, both `src` and `dst` are simple `FusionTreePair`s, and the +transformation consists of a single coefficient `coeff`. +For generic `FusionStyle`s, the input and output consist of `FusionTreeBlock`s that bundle together +all trees with the same uncoupled charges, and `coeffs` now forms a transformation matrix. + +``` + ╰─┬─╯ | ╭─╮ ╰─┬─╯ | + ╰─┬─╯ | | ╰─┬─╯ + ╰ ⋯ ┬╯ | ╰ ⋯ ╮ + | | → ╭─┴─╮ + ╭ ⋯ ┴╮ | ╭ ⋯ ┴╮ | + ╭─┴─╮ | | ╭─┴─╮ | | + ╭─┴─╮ | | | ╭─┴─╮ | | | +``` + +See also [`bendright`](@ref). +""" bendleft + +function bendleft((f₁, f₂)::FusionTreePair) + @assert FusionStyle((f₁, f₂)) === UniqueFusion() + (f₂′, f₁′), coeff = bendright((f₂, f₁)) + return (f₁′, f₂′) => conj(coeff) +end + +# !! note that this is more or less a copy of bendright through +# (f1, f2) => conj(coeff) for ((f2, f1), coeff) in bendleft(src) +function bendleft(src::FusionTreeBlock) + uncoupled_dst = ( + (src.uncoupled[1]..., dual(src.uncoupled[2][end])), + TupleTools.front(src.uncoupled[2]), + ) + isdual_dst = ( + (src.isdual[1]..., !(src.isdual[2][end])), + TupleTools.front(src.isdual[2]), + ) + I = sectortype(src) + N₁ = numin(src) + N₂ = numout(src) + @assert N₁ > 0 + + dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src)) + indexmap = treeindex_map(dst) + U = zeros(sectorscalartype(I), length(dst), length(src)) + + for (col, (f₂, f₁)) in enumerate(fusiontrees(src)) + c = f₁.coupled + a = N₁ == 1 ? leftunit(f₁.uncoupled[1]) : + (N₁ == 2 ? f₁.uncoupled[1] : f₁.innerlines[end]) + b = f₁.uncoupled[N₁] + + uncoupled1 = TupleTools.front(f₁.uncoupled) + isdual1 = TupleTools.front(f₁.isdual) + inner1 = N₁ > 2 ? TupleTools.front(f₁.innerlines) : () + vertices1 = N₁ > 1 ? TupleTools.front(f₁.vertices) : () + f₁′ = FusionTree(uncoupled1, a, isdual1, inner1, vertices1) + + uncoupled2 = (f₂.uncoupled..., dual(b)) + isdual2 = (f₂.isdual..., !(f₁.isdual[N₁])) + inner2 = N₂ > 1 ? (f₂.innerlines..., c) : () + + coeff₀ = sqrtdim(c) * invsqrtdim(a) + if f₁.isdual[N₁] + coeff₀ *= conj(frobenius_schur_phase(dual(b))) + end + if FusionStyle(I) isa MultiplicityFreeFusion + coeff = coeff₀ * Bsymbol(a, b, c) + vertices2 = N₂ > 0 ? (f₂.vertices..., 1) : () + f₂′ = FusionTree(uncoupled2, a, isdual2, inner2, vertices2) + row = indexmap[treeindex_data((f₂′, f₁′))] + @inbounds U[row, col] = conj(coeff) + else + Bmat = Bsymbol(a, b, c) + μ = N₁ > 1 ? f₁.vertices[end] : 1 + for ν in axes(Bmat, 2) + coeff = coeff₀ * Bmat[μ, ν] + iszero(coeff) && continue + vertices2 = N₂ > 0 ? (f₂.vertices..., ν) : () + f₂′ = FusionTree(uncoupled2, a, isdual2, inner2, vertices2) + row = indexmap[treeindex_data((f₂′, f₁′))] + @inbounds U[row, col] = conj(coeff) + end + end + end + + return dst => U +end + + +@doc """ + foldright((f₁, f₂)::FusionTreePair) -> (f₃, f₄) => coeff + foldright(src::FusionTreeBlock) -> dst => coeffs + +Map the first splitting vertex `a ⊗ b ← c` of `src` to a fusion vertex `a ← c ⊗ dual(b)` in `dst`. +For `FusionStyle(src) === UniqueFusion()`, both `src` and `dst` are simple `FusionTreePair`s, and the +transformation consists of a single coefficient `coeff`. +For generic `FusionStyle`s, the input and output consist of `FusionTreeBlock`s that bundle together +all trees with the same uncoupled charges, and `coeffs` now forms a transformation matrix. + +``` + | ╰─┬─╯ | | ╰─┬─╯ | | | + | ╰─┬─╯ | ╰─┬─╯ | | + | ╰ ⋯ ┬╯ ╰─┬─╯ | + | | → ╰ ⋯ ┬╯ + | ╭ ⋯ ┴╮ | + | ╭─┴─╮ | ╭─ ⋯ ┴╮ + ╰───┴─╮ | | ╭─┴─╮ | +``` + +See also [`foldleft`](@ref). +""" foldright + +function foldright((f₁, f₂)::FusionTreePair) + I = sectortype((f₁, f₂)) + @assert FusionStyle(I) === UniqueFusion() + @assert length(f₁) > 0 + + # compute new trees + a = f₁.uncoupled[1] + isduala = f₁.isdual[1] + c1 = dual(a) + c2 = f₁.coupled + c = first(c1 ⊗ c2) + fl = FusionTree{I}(Base.tail(f₁.uncoupled), c, Base.tail(f₁.isdual)) + fr = FusionTree{I}((c1, f₂.uncoupled...), c, (!isduala, f₂.isdual...)) + + # compute new coefficients + factor = sqrtdim(a) + isduala || (factor *= conj(frobenius_schur_phase(a))) + + return (fl, fr) => factor +end + +function foldright(src::FusionTreeBlock) + uncoupled_dst = ( + Base.tail(src.uncoupled[1]), + (dual(first(src.uncoupled[1])), src.uncoupled[2]...), + ) + isdual_dst = (Base.tail(src.isdual[1]), (!first(src.isdual[1]), src.isdual[2]...)) + I = sectortype(src) + N₁ = numout(src) + @assert N₁ > 0 + + dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src)) + indexmap = treeindex_map(dst) + U = zeros(sectorscalartype(I), length(dst), length(src)) + + for (col, (f₁, f₂)) in enumerate(fusiontrees(src)) + # map first splitting vertex (a, b)<-c to fusion vertex b<-(dual(a), c) + a = f₁.uncoupled[1] + isduala = f₁.isdual[1] + factor = sqrtdim(a) + if !isduala + factor *= conj(frobenius_schur_phase(a)) + end + c1 = dual(a) + c2 = f₁.coupled + uncoupled = Base.tail(f₁.uncoupled) + isdual = Base.tail(f₁.isdual) + if FusionStyle(I) isa UniqueFusion + c = first(c1 ⊗ c2) + fl = FusionTree{I}(Base.tail(f₁.uncoupled), c, Base.tail(f₁.isdual)) + fr = FusionTree{I}((c1, f₂.uncoupled...), c, (!isduala, f₂.isdual...)) + row = indexmap[treeindex_data((fl, fr))] + @inbounds U[row, col] = factor + else + if N₁ == 1 + cset = (leftunit(c1),) # or rightunit(a) + elseif N₁ == 2 + cset = (f₁.uncoupled[2],) + else + cset = ⊗(Base.tail(f₁.uncoupled)...) + end + for c in c1 ⊗ c2 + c ∈ cset || continue + for μ in 1:Nsymbol(c1, c2, c) + fc = FusionTree((c1, c2), c, (!isduala, false), (), (μ,)) + frs_coeffs = insertat(fc, 2, f₂) + for (fl′, coeff1) in insertat(fc, 2, f₁) + N₁ > 1 && !isone(fl′.innerlines[1]) && continue + coupled = fl′.coupled + uncoupled = Base.tail(Base.tail(fl′.uncoupled)) + isdual = Base.tail(Base.tail(fl′.isdual)) + inner = N₁ <= 3 ? () : Base.tail(Base.tail(fl′.innerlines)) + vertices = N₁ <= 2 ? () : Base.tail(Base.tail(fl′.vertices)) + fl = FusionTree{I}(uncoupled, coupled, isdual, inner, vertices) + for (fr, coeff2) in frs_coeffs + coeff = factor * coeff1 * conj(coeff2) + row = indexmap[treeindex_data((fl, fr))] + @inbounds U[row, col] = coeff + end + end + end + end + end + end + + return dst => U +end + +@doc """ + foldleft((f₁, f₂)::FusionTreePair) -> (f₃, f₄) => coeff + foldleft(src::FusionTreeBlock) -> dst => coeffs + +Map the first fusion vertex `a ← c ⊗ dual(b)` of `src` to a splitting vertex `a ⊗ b ← c` in `dst`. +For `FusionStyle(src) === UniqueFusion()`, both `src` and `dst` are simple `FusionTreePair`s, and the +transformation consists of a single coefficient `coeff`. +For generic `FusionStyle`s, the input and output consist of `FusionTreeBlock`s that bundle together +all trees with the same uncoupled charges, and `coeffs` now forms a transformation matrix. + +``` + ╭───┬─╯ | | ╰─┬─╯ | + | ╰─┬─╯ | ╰ ⋯ ┬╯ + | ╰ ⋯ ┬╯ | + | | → ╭ ⋯ ┴╮ + | ╭ ⋯ ┴╮ ╭─┴─╮ | + | ╭─┴─╮ | ╭─┴─╮ | | + | ╭─┴─╮ | | ╭─┴─╮ | | | +``` + +See also [`foldright`](@ref). +""" foldleft + +function foldleft((f₁, f₂)::FusionTreePair) + @assert FusionStyle((f₁, f₂)) === UniqueFusion() + (f₂′, f₁′), coeff = foldright((f₂, f₁)) + return (f₁′, f₂′) => conj(coeff) +end + +# !! note that this is more or less a copy of foldright through +# (f1, f2) => conj(coeff) for ((f2, f1), coeff) in foldright(src) +function foldleft(src::FusionTreeBlock) + uncoupled_dst = ( + (dual(first(src.uncoupled[2])), src.uncoupled[1]...), + Base.tail(src.uncoupled[2]), + ) + isdual_dst = ( + (!first(src.isdual[2]), src.isdual[1]...), + Base.tail(src.isdual[2]), + ) + I = sectortype(src) + N₁ = numin(src) + N₂ = numout(src) + @assert N₁ > 0 + + dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src)) + indexmap = treeindex_map(dst) + U = zeros(sectorscalartype(I), length(dst), length(src)) + + for (col, (f₂, f₁)) in enumerate(fusiontrees(src)) + # map first splitting vertex (a, b)<-c to fusion vertex b<-(dual(a), c) + a = f₁.uncoupled[1] + isduala = f₁.isdual[1] + factor = sqrtdim(a) + if !isduala + factor *= conj(frobenius_schur_phase(a)) + end + c1 = dual(a) + c2 = f₁.coupled + uncoupled = Base.tail(f₁.uncoupled) + isdual = Base.tail(f₁.isdual) + if FusionStyle(I) isa UniqueFusion + c = first(c1 ⊗ c2) + fl = FusionTree{I}(Base.tail(f₁.uncoupled), c, Base.tail(f₁.isdual)) + fr = FusionTree{I}((c1, f₂.uncoupled...), c, (!isduala, f₂.isdual...)) + row = indexmap[treeindex_data((fr, fl))] + @inbounds U[row, col] = conj(factor) + else + if N₁ == 1 + cset = (leftunit(c1),) # or rightunit(a) + elseif N₁ == 2 + cset = (f₁.uncoupled[2],) + else + cset = ⊗(Base.tail(f₁.uncoupled)...) + end + for c in c1 ⊗ c2 + c ∈ cset || continue + for μ in 1:Nsymbol(c1, c2, c) + fc = FusionTree((c1, c2), c, (!isduala, false), (), (μ,)) + fr_coeffs = insertat(fc, 2, f₂) + for (fl′, coeff1) in insertat(fc, 2, f₁) + N₁ > 1 && !isone(fl′.innerlines[1]) && continue + coupled = fl′.coupled + uncoupled = Base.tail(Base.tail(fl′.uncoupled)) + isdual = Base.tail(Base.tail(fl′.isdual)) + inner = N₁ <= 3 ? () : Base.tail(Base.tail(fl′.innerlines)) + vertices = N₁ <= 2 ? () : Base.tail(Base.tail(fl′.vertices)) + fl = FusionTree{I}(uncoupled, coupled, isdual, inner, vertices) + for (fr, coeff2) in fr_coeffs + coeff = factor * coeff1 * conj(coeff2) + row = indexmap[treeindex_data((fr, fl))] + @inbounds U[row, col] = conj(coeff) + end + end + end + end + end + end + return dst => U +end + +# clockwise cyclic permutation while preserving (N₁, N₂): foldright & bendleft +# anticlockwise cyclic permutation while preserving (N₁, N₂): foldleft & bendright +# These are utility functions that preserve the type of the input/output trees, +# and are therefore used to craft type-stable transpose implementations. + +function cycleclockwise(src::Union{FusionTreePair, FusionTreeBlock}) + if numout(src) > 0 + tmp, U₁ = foldright(src) + dst, U₂ = bendleft(tmp) + else + tmp, U₁ = bendleft(src) + dst, U₂ = foldright(tmp) + end + return dst => U₂ * U₁ +end +function cycleanticlockwise(src::Union{FusionTreePair, FusionTreeBlock}) + if numin(src) > 0 + tmp, U₁ = foldleft(src) + dst, U₂ = bendright(tmp) + else + tmp, U₁ = bendright(src) + dst, U₂ = foldleft(tmp) + end + return dst => U₂ * U₁ +end + +# COMPOSITE DUALITY MANIPULATIONS PART 1: Repartition and transpose +#------------------------------------------------------------------- +# -> composite manipulations that depend on the duality (rigidity) and pivotal structure +# -> planar manipulations that do not require braiding, everything is in Fsymbol (A/Bsymbol) +# -> transpose expressed as cyclic permutation + +# repartition double fusion tree +""" + repartition((f₁, f₂)::FusionTreePair{I, N₁, N₂}, N::Int) where {I, N₁, N₂} + -> <:AbstractDict{<:FusionTreePair{I, N, N₁+N₂-N}}, <:Number} + +Input is a double fusion tree that describes the fusion of a set of incoming uncoupled +sectors to a set of outgoing uncoupled sectors, represented using the individual trees of +outgoing (`f₁`) and incoming sectors (`f₂`) respectively (with identical coupled sector +`f₁.coupled == f₂.coupled`). Computes new trees and corresponding coefficients obtained from +repartitioning the tree by bending incoming to outgoing sectors (or vice versa) in order to +have `N` outgoing sectors. +""" +@inline function repartition(src::Union{FusionTreePair, FusionTreeBlock}, N::Int) + @assert 0 <= N <= numind(src) + return repartition(src, Val(N)) +end + +#= +Using a generated function here to ensure type stability by unrolling the loops: +```julia +dst, U = bendleft/right(src) + +# repeat the following 2 lines N - 1 times +dst, Utmp = bendleft/right(dst) +U = Utmp * U + +return dst, U +``` +=# +function _repartition_body(N) + if N == 0 + ex = quote + T = sectorscalartype(sectortype(src)) + if FusionStyle(src) === UniqueFusion() + return src => one(T) + else + U = copyto!(zeros(T, length(src), length(src)), LinearAlgebra.I) + return src, U + end + end + else + f = N < 0 ? bendleft : bendright + ex_rep = Expr(:block) + for _ in 1:(abs(N) - 1) + push!(ex_rep.args, :((dst, Utmp) = $f(dst))) + push!(ex_rep.args, :(U = Utmp * U)) + end + ex = quote + dst, U = $f(src) + $ex_rep + return dst => U + end + end + return ex +end +@generated function repartition(src::Union{FusionTreePair, FusionTreeBlock}, ::Val{N}) where {N} + return _repartition_body(numout(src) - N) +end + +""" + transpose((f₁, f₂)::FusionTreePair{I}, p::Index2Tuple{N₁, N₂}) where {I, N₁, N₂} + -> <:AbstractDict{<:FusionTreePair{I, N₁, N₂}}, <:Number} + +Input is a double fusion tree that describes the fusion of a set of incoming uncoupled +sectors to a set of outgoing uncoupled sectors, represented using the individual trees of +outgoing (`t1`) and incoming sectors (`t2`) respectively (with identical coupled sector +`t1.coupled == t2.coupled`). Computes new trees and corresponding coefficients obtained from +repartitioning and permuting the tree such that sectors `p1` become outgoing and sectors +`p2` become incoming. +""" +function Base.transpose(src::Union{FusionTreePair, FusionTreeBlock}, p::Index2Tuple) + numind(src) == length(p[1]) + length(p[2]) || throw(ArgumentError("invalid permutation p")) + p′ = linearizepermutation(p..., numout(src), numin(src)) + iscyclicpermutation(p′) || throw(ArgumentError("invalid permutation p")) + return fstranspose((src, p)) +end + +const FSPTransposeKey{I, N₁, N₂} = Tuple{FusionTreePair{I}, Index2Tuple{N₁, N₂}} +const FSBTransposeKey{I, N₁, N₂} = Tuple{FusionTreeBlock{I}, Index2Tuple{N₁, N₂}} + +Base.@assume_effects :foldable function _fsdicttype(::Type{T}) where {I, N₁, N₂, T <: FSPTransposeKey{I, N₁, N₂}} + E = sectorscalartype(I) + return Pair{fusiontreetype(I, N₁, N₂), E} +end +Base.@assume_effects :foldable function _fsdicttype(::Type{T}) where {I, N₁, N₂, T <: FSBTransposeKey{I, N₁, N₂}} + F₁ = fusiontreetype(I, N₁) + F₂ = fusiontreetype(I, N₂) + E = sectorscalartype(I) + return Pair{FusionTreeBlock{I, N₁, N₂, Tuple{F₁, F₂}}, Matrix{E}} +end + +@cached function fstranspose(key::K)::_fsdicttype(K) where {I, N₁, N₂, K <: Union{FSPTransposeKey{I, N₁, N₂}, FSBTransposeKey{I, N₁, N₂}}} + src, (p1, p2) = key + + N = N₁ + N₂ + p = linearizepermutation(p1, p2, numout(src), numin(src)) + + dst, U = repartition(src, N₁) + length(p) == 0 && return dst => U + i1 = findfirst(==(1), p)::Int + i1 == 1 && return dst => U + + Nhalf = N >> 1 + while 1 < i1 ≤ Nhalf + dst, U_tmp = cycleanticlockwise(dst) + U = U_tmp * U + i1 -= 1 + end + while Nhalf < i1 + dst, U_tmp = cycleclockwise(dst) + U = U_tmp * U + i1 = mod1(i1 + 1, N) + end + + return dst => U +end + +CacheStyle(::typeof(fstranspose), k::FSPTransposeKey{I}) where {I} = + FusionStyle(I) isa UniqueFusion ? NoCache() : GlobalLRUCache() +CacheStyle(::typeof(fstranspose), k::FSBTransposeKey{I}) where {I} = + FusionStyle(I) isa UniqueFusion ? NoCache() : GlobalLRUCache() + +# COMPOSITE DUALITY MANIPULATIONS PART 2: Planar traces +#------------------------------------------------------------------- +# -> composite manipulations that depend on the duality (rigidity) and pivotal structure +# -> planar manipulations that do not require braiding, everything is in Fsymbol (A/Bsymbol) + +function planar_trace( + (f₁, f₂)::FusionTreePair{I}, (p1, p2)::Index2Tuple{N₁, N₂}, (q1, q2)::Index2Tuple{N₃, N₃} + ) where {I, N₁, N₂, N₃} + N = N₁ + N₂ + 2N₃ + @assert length(f₁) + length(f₂) == N + if N₃ == 0 + return transpose((f₁, f₂), (p1, p2)) + end + + linearindex = ( + ntuple(identity, Val(length(f₁)))..., + reverse(length(f₁) .+ ntuple(identity, Val(length(f₂))))..., + ) + + q1′ = TupleTools.getindices(linearindex, q1) + q2′ = TupleTools.getindices(linearindex, q2) + p1′, p2′ = let q′ = (q1′..., q2′...) + ( + map(l -> l - count(l .> q′), TupleTools.getindices(linearindex, p1)), + map(l -> l - count(l .> q′), TupleTools.getindices(linearindex, p2)), + ) + end + + T = sectorscalartype(I) + F₁ = fusiontreetype(I, N₁) + F₂ = fusiontreetype(I, N₂) + newtrees = FusionTreeDict{Tuple{F₁, F₂}, T}() + if FusionStyle(I) isa UniqueFusion + (f₁′, f₂′), coeff′ = repartition((f₁, f₂), N) + for (f₁′′, coeff′′) in planar_trace(f₁′, (q1′, q2′)) + (f12′′′, coeff′′′) = transpose((f₁′′, f₂′), (p1′, p2′)) + coeff = coeff′ * coeff′′ * coeff′′′ + iszero(coeff) || (newtrees[f12′′′] = get(newtrees, f12′′′, zero(coeff)) + coeff) + end + else + # TODO: this is a bit of a hack to fix the traces for now + src = FusionTreeBlock([(f₁, f₂)]) + dst, U = repartition(src, N) + for ((f₁′, f₂′), coeff′) in zip(fusiontrees(dst), U) + for (f₁′′, coeff′′) in planar_trace(f₁′, (q1′, q2′)) + src′ = FusionTreeBlock([(f₁′′, f₂′)]) + dst′, U′ = transpose(src′, (p1′, p2′)) + for (f12′′′, coeff′′′) in zip(fusiontrees(dst′), U′) + coeff = coeff′ * coeff′′ * coeff′′′ + iszero(coeff) || (newtrees[f12′′′] = get(newtrees, f12′′′, zero(coeff)) + coeff) + end + end + end + end + return newtrees +end + +""" + planar_trace(f::FusionTree{I,N}, (q1, q2)::Index2Tuple{N₃,N₃}) where {I,N,N₃} + -> <:AbstractDict{FusionTree{I,N-2*N₃}, <:Number} + +Perform a planar trace of the uncoupled indices of the fusion tree `f` at `q1` with those at +`q2`, where `q1[i]` is connected to `q2[i]` for all `i`. The result is returned as a dictionary +of output trees and corresponding coefficients. +""" +function planar_trace(f::FusionTree{I, N}, (q1, q2)::Index2Tuple{N₃, N₃}) where {I, N, N₃} + T = sectorscalartype(I) + F = fusiontreetype(I, N - 2 * N₃) + newtrees = FusionTreeDict{F, T}() + N₃ === 0 && return push!(newtrees, f => one(T)) + + for (i, j) in zip(q1, q2) + (f.uncoupled[i] == dual(f.uncoupled[j]) && f.isdual[i] != f.isdual[j]) || + return newtrees + end + k = 1 + local i, j + while k <= N₃ + if mod1(q1[k] + 1, N) == q2[k] + i = q1[k] + j = q2[k] + break + elseif mod1(q2[k] + 1, N) == q1[k] + i = q2[k] + j = q1[k] + break + else + k += 1 + end + end + k > N₃ && throw(ArgumentError("Not a planar trace")) + + q1′ = let i = i, j = j + map(l -> (l - (l > i) - (l > j)), TupleTools.deleteat(q1, k)) + end + q2′ = let i = i, j = j + map(l -> (l - (l > i) - (l > j)), TupleTools.deleteat(q2, k)) + end + for (f′, coeff′) in elementary_trace(f, i) + for (f′′, coeff′′) in planar_trace(f′, (q1′, q2′)) + coeff = coeff′ * coeff′′ + if !iszero(coeff) + newtrees[f′′] = get(newtrees, f′′, zero(coeff)) + coeff + end + end + end + return newtrees +end + +# trace two neighbouring indices of a single fusion tree +""" + elementary_trace(f::FusionTree{I,N}, i) where {I,N} -> <:AbstractDict{FusionTree{I,N-2}, <:Number} + +Perform an elementary trace of neighbouring uncoupled indices `i` and +`i+1` on a fusion tree `f`, and returns the result as a dictionary of output trees and +corresponding coefficients. +""" +function elementary_trace(f::FusionTree{I, N}, i) where {I, N} + (N > 1 && 1 <= i <= N) || + throw(ArgumentError("Cannot trace outputs i=$i and i+1 out of only $N outputs")) + i < N || isunit(f.coupled) || + throw(ArgumentError("Cannot trace outputs i=$N and 1 of fusion tree that couples to non-trivial sector")) + + T = sectorscalartype(I) + F = fusiontreetype(I, N - 2) + newtrees = FusionTreeDict{F, T}() + + j = mod1(i + 1, N) + b = f.uncoupled[i] + b′ = f.uncoupled[j] + # if trace is zero, return empty dict + (b == dual(b′) && f.isdual[i] != f.isdual[j]) || return newtrees + if i < N + inner_extended = (leftunit(f.uncoupled[1]), f.uncoupled[1], f.innerlines..., f.coupled) + a = inner_extended[i] + d = inner_extended[i + 2] + a == d || return newtrees + uncoupled′ = TupleTools.deleteat(TupleTools.deleteat(f.uncoupled, i + 1), i) + isdual′ = TupleTools.deleteat(TupleTools.deleteat(f.isdual, i + 1), i) + coupled′ = f.coupled + if N <= 4 + inner′ = () + else + inner′ = i <= 2 ? Base.tail(Base.tail(f.innerlines)) : + TupleTools.deleteat(TupleTools.deleteat(f.innerlines, i - 1), i - 2) + end + if N <= 3 + vertices′ = () + else + vertices′ = i <= 2 ? Base.tail(Base.tail(f.vertices)) : + TupleTools.deleteat(TupleTools.deleteat(f.vertices, i), i - 1) + end + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′) + coeff = sqrtdim(b) + if i > 1 + c = f.innerlines[i - 1] + if FusionStyle(I) isa MultiplicityFreeFusion + coeff *= Fsymbol(a, b, dual(b), a, c, rightunit(a)) + else + μ = f.vertices[i - 1] + ν = f.vertices[i] + coeff *= Fsymbol(a, b, dual(b), a, c, rightunit(a))[μ, ν, 1, 1] + end + end + if f.isdual[i] + coeff *= frobenius_schur_phase(b) + end + push!(newtrees, f′ => coeff) + return newtrees + else # i == N + unit = leftunit(b) + if N == 2 + f′ = FusionTree{I}((), unit, (), (), ()) + coeff = sqrtdim(b) + if !(f.isdual[N]) + coeff *= conj(frobenius_schur_phase(b)) + end + push!(newtrees, f′ => coeff) + return newtrees + end + uncoupled_ = TupleTools.front(f.uncoupled) + inner_ = TupleTools.front(f.innerlines) + coupled_ = f.innerlines[end] + isdual_ = TupleTools.front(f.isdual) + vertices_ = TupleTools.front(f.vertices) + f_ = FusionTree(uncoupled_, coupled_, isdual_, inner_, vertices_) + fs = FusionTree((b,), b, (!f.isdual[1],), (), ()) + for (f_′, coeff) in merge(fs, f_, unit, 1) + f_′.innerlines[1] == unit || continue + uncoupled′ = Base.tail(Base.tail(f_′.uncoupled)) + isdual′ = Base.tail(Base.tail(f_′.isdual)) + inner′ = N <= 4 ? () : Base.tail(Base.tail(f_′.innerlines)) + vertices′ = N <= 3 ? () : Base.tail(Base.tail(f_′.vertices)) + f′ = FusionTree(uncoupled′, unit, isdual′, inner′, vertices′) + coeff *= sqrtdim(b) + if !(f.isdual[N]) + coeff *= conj(frobenius_schur_phase(b)) + end + newtrees[f′] = get(newtrees, f′, zero(coeff)) + coeff + end + return newtrees + end +end diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index 1665cdb27..895a5d619 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -105,15 +105,98 @@ function FusionTree( end FusionTree(uncoupled::Tuple{I, Vararg{I}}) where {I <: Sector} = FusionTree(uncoupled, unit(I)) +""" + FusionTreePair{I, N₁, N₂} + +Type alias for a fusion-splitting tree pair of sectortype `I`, with `N₁` splitting legs and +`N₂` fusion legs. +""" +const FusionTreePair{I, N₁, N₂} = Tuple{FusionTree{I, N₁}, FusionTree{I, N₂}} + +""" + FusionTreeBlock{I, N₁, N₂, F <: FusionTreePair{I, N₁, N₂}} + +Collection of fusion-splitting tree pairs that share the same uncoupled sectors. +Mostly internal structure to speed up non-`UniqueFusion` fusiontree manipulation computations. +""" +struct FusionTreeBlock{I, N₁, N₂, F <: FusionTreePair{I, N₁, N₂}} + trees::Vector{F} +end + +function FusionTreeBlock{I}( + uncoupled::Tuple{NTuple{N₁, I}, NTuple{N₂, I}}, + isdual::Tuple{NTuple{N₁, Bool}, NTuple{N₂, Bool}}; + sizehint::Int = 0 + ) where {I <: Sector, N₁, N₂} + F = fusiontreetype(I, N₁, N₂) + trees = Vector{F}(undef, 0) + sizehint > 0 && sizehint!(trees, sizehint) + + if N₁ == N₂ == 0 + return FusionTreeBlock(trees) + elseif N₁ == 0 + cs = sort!(collect(filter(isone, ⊗(uncoupled[2]...)))) + elseif N₂ == 0 + cs = sort!(collect(filter(isone, ⊗(uncoupled[1]...)))) + else + cs = sort!(collect(intersect(⊗(uncoupled[1]...), ⊗(uncoupled[2]...)))) + end + + for c in cs + for f₁ in fusiontrees(uncoupled[1], c, isdual[1]), + f₂ in fusiontrees(uncoupled[2], c, isdual[2]) + + push!(trees, (f₁, f₂)) + end + end + return FusionTreeBlock(trees) +end + # Properties +# ---------- sectortype(::Type{<:FusionTree{I}}) where {I <: Sector} = I -FusionStyle(::Type{<:FusionTree{I}}) where {I <: Sector} = FusionStyle(I) -BraidingStyle(::Type{<:FusionTree{I}}) where {I <: Sector} = BraidingStyle(I) -Base.length(::Type{<:FusionTree{<:Sector, N}}) where {N} = N +sectortype(::Type{<:FusionTreePair{I}}) where {I <: Sector} = I +sectortype(::Type{<:FusionTreeBlock{I}}) where {I} = I + +FusionStyle(f::Union{FusionTree, FusionTreePair, FusionTreeBlock}) = + FusionStyle(typeof(f)) +FusionStyle(::Type{F}) where {F <: Union{FusionTree, FusionTreePair, FusionTreeBlock}} = + FusionStyle(sectortype(F)) + +BraidingStyle(f::Union{FusionTree, FusionTreePair, FusionTreeBlock}) = + BraidingStyle(typeof(f)) +BraidingStyle(::Type{F}) where {F <: Union{FusionTree, FusionTreePair, FusionTreeBlock}} = + BraidingStyle(sectortype(F)) -FusionStyle(f::FusionTree) = FusionStyle(typeof(f)) -BraidingStyle(f::FusionTree) = BraidingStyle(typeof(f)) Base.length(f::FusionTree) = length(typeof(f)) +Base.length(::Type{<:FusionTree{<:Sector, N}}) where {N} = N +# Note: cannot define the following since FusionTreePair is a const for a Tuple +# Base.length(::Type{<:FusionTreePair{<:Sector, N₁, N₂}}) where {N₁, N₂} = N₁ + N₂ +# Base.length(f::FusionTreePair) = length(typeof(f)) +Base.length(block::FusionTreeBlock) = length(fusiontrees(block)) + +numout(fs::Union{FusionTreePair, FusionTreeBlock}) = numout(typeof(fs)) +numout(::Type{<:FusionTreePair{I, N₁}}) where {I, N₁} = N₁ +numout(::Type{<:FusionTreeBlock{I, N₁}}) where {I, N₁} = N₁ +numin(fs::Union{FusionTreePair, FusionTreeBlock}) = numin(typeof(fs)) +numin(::Type{<:FusionTreePair{I, N₁, N₂}}) where {I, N₁, N₂} = N₂ +numin(::Type{<:FusionTreeBlock{I, N₁, N₂}}) where {I, N₁, N₂} = N₂ +numind(fs::Union{FusionTreePair, FusionTreeBlock}) = numind(typeof(fs)) +numind(::Type{T}) where {T <: Union{FusionTreePair, FusionTreeBlock}} = numin(T) + numout(T) + +Base.propertynames(::FusionTreeBlock, private::Bool = false) = (:trees, :uncoupled, :isdual) +Base.@constprop :aggressive function Base.getproperty(block::FusionTreeBlock, prop::Symbol) + if prop === :uncoupled + f₁, f₂ = first(block.trees) + return f₁.uncoupled, f₂.uncoupled + elseif prop === :isdual + f₁, f₂ = first(block.trees) + return f₁.isdual, f₂.isdual + else + return getfield(block, prop) + end +end +fusiontrees(block::FusionTreeBlock) = block.trees # Hashing, important for using fusion trees as key in a dictionary function Base.hash(f::FusionTree{I}, h::UInt) where {I} @@ -146,8 +229,26 @@ function Base.:(==)(f₁::FusionTree{I, N}, f₂::FusionTree{I, N}) where {I <: end Base.:(==)(f₁::FusionTree, f₂::FusionTree) = false +# Within one block, all values of uncoupled and isdual are equal, so avoid hashing these +function treeindex_data((f₁, f₂)) + I = sectortype(f₁) + if FusionStyle(I) isa GenericFusion + return (f₁.coupled, f₁.innerlines..., f₂.innerlines...), + (f₁.vertices..., f₂.vertices...) + elseif FusionStyle(I) isa MultipleFusion + return (f₁.coupled, f₁.innerlines..., f₂.innerlines...) + else # there should be only a single element anyways + return false + end +end +function treeindex_map(fs::FusionTreeBlock) + I = sectortype(fs) + return fusiontreedict(I)(treeindex_data(f) => ind for (ind, f) in enumerate(fusiontrees(fs))) +end + + # Facilitate getting correct fusion tree types -function fusiontreetype(::Type{I}, N::Int) where {I <: Sector} +Base.@assume_effects :foldable function fusiontreetype(::Type{I}, N::Int) where {I <: Sector} return if N === 0 FusionTree{I, 0, 0, 0} elseif N === 1 @@ -156,52 +257,57 @@ function fusiontreetype(::Type{I}, N::Int) where {I <: Sector} FusionTree{I, N, N - 2, N - 1} end end +Base.@assume_effects :foldable function fusiontreetype(::Type{I}, N₁::Int, N₂::Int) where {I <: Sector} + return Tuple{fusiontreetype(I, N₁), fusiontreetype(I, N₂)} +end + +fusiontreedict(I) = FusionStyle(I) isa UniqueFusion ? SingletonDict : FusionTreeDict # converting to actual array -function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, 0}) where {I} - X = convert(A, fusiontensor(unit(I), unit(I), unit(I)))[1, 1, :] - return X -end -function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, 1}) where {I} +Base.convert(A::Type{<:AbstractArray}, f::FusionTree) = convert(A, fusiontensor(f)) +# TODO: is this piracy? +Base.convert(A::Type{<:AbstractArray}, (f₁, f₂)::FusionTreePair) = + convert(A, fusiontensor((f₁, f₂))) + +fusiontensor(::FusionTree{I, 0}) where {I} = fusiontensor(unit(I), unit(I), unit(I))[1, 1, :] +function fusiontensor(f::FusionTree{I, 1}) where {I} c = f.coupled if f.isdual[1] sqrtdc = sqrtdim(c) - Zcbartranspose = sqrtdc * convert(A, fusiontensor(dual(c), c, unit(c)))[:, :, 1, 1] + Zcbartranspose = sqrtdc * fusiontensor(dual(c), c, unit(c))[:, :, 1, 1] X = conj!(Zcbartranspose) # we want Zcbar^† else - X = convert(A, fusiontensor(c, unit(c), c))[:, 1, :, 1, 1] + X = fusiontensor(c, unit(c), c)[:, 1, :, 1, 1] end return X end - -function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, 2}) where {I} +function fusiontensor(f::FusionTree{I, 2}) where {I} a, b = f.uncoupled isduala, isdualb = f.isdual c = f.coupled μ = (FusionStyle(I) isa GenericFusion) ? f.vertices[1] : 1 - C = convert(A, fusiontensor(a, b, c))[:, :, :, μ] + C = fusiontensor(a, b, c)[:, :, :, μ] X = C if isduala - Za = convert(A, FusionTree((a,), a, (isduala,), ())) + Za = fusiontensor(FusionTree((a,), a, (isduala,), ())) @tensor X[a′, b, c] := Za[a′, a] * X[a, b, c] end if isdualb - Zb = convert(A, FusionTree((b,), b, (isdualb,), ())) + Zb = fusiontensor(FusionTree((b,), b, (isdualb,), ())) @tensor X[a, b′, c] := Zb[b′, b] * X[a, b, c] end return X end - -function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, N}) where {I, N} +function fusiontensor(f::FusionTree{I, N}) where {I, N} tailout = (f.innerlines[1], TupleTools.tail2(f.uncoupled)...) isdualout = (false, TupleTools.tail2(f.isdual)...) ftail = FusionTree(tailout, f.coupled, isdualout, Base.tail(f.innerlines), Base.tail(f.vertices)) - Ctail = convert(A, ftail) + Ctail = fusiontensor(ftail) f₁ = FusionTree( (f.uncoupled[1], f.uncoupled[2]), f.innerlines[1], (f.isdual[1], f.isdual[2]), (), (f.vertices[1],) ) - C1 = convert(A, f₁) + C1 = fusiontensor(f₁) dtail = size(Ctail) d1 = size(C1) X = similar(C1, (d1[1], d1[2], Base.tail(dtail)...)) @@ -214,22 +320,19 @@ function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, N}) where {I, N ) end -# TODO: is this piracy? -function Base.convert( - A::Type{<:AbstractArray}, (f₁, f₂)::Tuple{FusionTree{I}, FusionTree{I}} - ) where {I} - F₁ = convert(A, f₁) - F₂ = convert(A, f₂) +function fusiontensor((f₁, f₂)::FusionTreePair) + F₁ = fusiontensor(f₁) + F₂ = fusiontensor(f₂) sz1 = size(F₁) sz2 = size(F₂) d1 = TupleTools.front(sz1) d2 = TupleTools.front(sz2) - return reshape( reshape(F₁, TupleTools.prod(d1), sz1[end]) * reshape(F₂, TupleTools.prod(d2), sz2[end])', (d1..., d2...) ) end +fusiontensor(src::FusionTreeBlock) = sum(fusiontensor, fusiontrees(src)) # Show methods function Base.show(io::IO, t::FusionTree{I}) where {I <: Sector} @@ -247,12 +350,14 @@ function Base.show(io::IO, t::FusionTree{I}) where {I <: Sector} end end -# Manipulate fusion trees -include("manipulations.jl") - # Fusion tree iterators include("iterator.jl") +# Manipulate fusion trees +include("basic_manipulations.jl") +include("duality_manipulations.jl") +include("braiding_manipulations.jl") + # auxiliary routines # _abelianinner: generate the inner indices for given outer indices in the abelian case _abelianinner(outer::Tuple{}) = () diff --git a/src/fusiontrees/manipulations.jl b/src/fusiontrees/manipulations.jl deleted file mode 100644 index 1564b1b67..000000000 --- a/src/fusiontrees/manipulations.jl +++ /dev/null @@ -1,1108 +0,0 @@ -fusiontreedict(I) = FusionStyle(I) isa UniqueFusion ? SingletonDict : FusionTreeDict - -# BASIC MANIPULATIONS: -#---------------------------------------------- -# -> rewrite generic fusion tree in basis of fusion trees in standard form -# -> only depend on Fsymbol - -""" - insertat(f::FusionTree{I, N₁}, i::Int, f₂::FusionTree{I, N₂}) - -> <:AbstractDict{<:FusionTree{I, N₁+N₂-1}, <:Number} - -Attach a fusion tree `f₂` to the uncoupled leg `i` of the fusion tree `f₁` and bring it -into a linear combination of fusion trees in standard form. This requires that -`f₂.coupled == f₁.uncoupled[i]` and `f₁.isdual[i] == false`. -""" -function insertat(f₁::FusionTree{I}, i::Int, f₂::FusionTree{I, 0}) where {I} - # this actually removes uncoupled line i, which should be trivial - (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) || - throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])")) - coeff = one(sectorscalartype(I)) - - uncoupled = TupleTools.deleteat(f₁.uncoupled, i) - coupled = f₁.coupled - isdual = TupleTools.deleteat(f₁.isdual, i) - if length(uncoupled) <= 2 - inner = () - else - inner = TupleTools.deleteat(f₁.innerlines, max(1, i - 2)) - end - if length(uncoupled) <= 1 - vertices = () - else - vertices = TupleTools.deleteat(f₁.vertices, max(1, i - 1)) - end - f = FusionTree(uncoupled, coupled, isdual, inner, vertices) - return fusiontreedict(I)(f => coeff) -end -function insertat(f₁::FusionTree{I}, i, f₂::FusionTree{I, 1}) where {I} - # identity operation - (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) || - throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])")) - coeff = one(sectorscalartype(I)) - isdual′ = TupleTools.setindex(f₁.isdual, f₂.isdual[1], i) - f = FusionTree{I}(f₁.uncoupled, f₁.coupled, isdual′, f₁.innerlines, f₁.vertices) - return fusiontreedict(I)(f => coeff) -end -function insertat(f₁::FusionTree{I}, i, f₂::FusionTree{I, 2}) where {I} - # elementary building block, - (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) || - throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])")) - uncoupled = f₁.uncoupled - coupled = f₁.coupled - inner = f₁.innerlines - b, c = f₂.uncoupled - isdual = f₁.isdual - isdualb, isdualc = f₂.isdual - if i == 1 - uncoupled′ = (b, c, tail(uncoupled)...) - isdual′ = (isdualb, isdualc, tail(isdual)...) - inner′ = (uncoupled[1], inner...) - vertices′ = (f₂.vertices..., f₁.vertices...) - coeff = one(sectorscalartype(I)) - f′ = FusionTree(uncoupled′, coupled, isdual′, inner′, vertices′) - return fusiontreedict(I)(f′ => coeff) - end - uncoupled′ = TupleTools.insertafter(TupleTools.setindex(uncoupled, b, i), i, (c,)) - isdual′ = TupleTools.insertafter(TupleTools.setindex(isdual, isdualb, i), i, (isdualc,)) - inner_extended = (uncoupled[1], inner..., coupled) - a = inner_extended[i - 1] - d = inner_extended[i] - e′ = uncoupled[i] - if FusionStyle(I) isa MultiplicityFreeFusion - local newtrees - for e in a ⊗ b - coeff = conj(Fsymbol(a, b, c, d, e, e′)) - iszero(coeff) && continue - inner′ = TupleTools.insertafter(inner, i - 2, (e,)) - f′ = FusionTree(uncoupled′, coupled, isdual′, inner′) - if @isdefined newtrees - push!(newtrees, f′ => coeff) - else - newtrees = fusiontreedict(I)(f′ => coeff) - end - end - return newtrees - else - local newtrees - κ = f₂.vertices[1] - λ = f₁.vertices[i - 1] - for e in a ⊗ b - inner′ = TupleTools.insertafter(inner, i - 2, (e,)) - Fmat = Fsymbol(a, b, c, d, e, e′) - for μ in axes(Fmat, 1), ν in axes(Fmat, 2) - coeff = conj(Fmat[μ, ν, κ, λ]) - iszero(coeff) && continue - vertices′ = TupleTools.setindex(f₁.vertices, ν, i - 1) - vertices′ = TupleTools.insertafter(vertices′, i - 2, (μ,)) - f′ = FusionTree(uncoupled′, coupled, isdual′, inner′, vertices′) - if @isdefined newtrees - push!(newtrees, f′ => coeff) - else - newtrees = fusiontreedict(I)(f′ => coeff) - end - end - end - return newtrees - end -end -function insertat(f₁::FusionTree{I, N₁}, i, f₂::FusionTree{I, N₂}) where {I, N₁, N₂} - F = fusiontreetype(I, N₁ + N₂ - 1) - (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) || - throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])")) - T = sectorscalartype(I) - coeff = one(T) - if length(f₁) == 1 - return fusiontreedict(I){F, T}(f₂ => coeff) - end - if i == 1 - uncoupled = (f₂.uncoupled..., tail(f₁.uncoupled)...) - isdual = (f₂.isdual..., tail(f₁.isdual)...) - inner = (f₂.innerlines..., f₂.coupled, f₁.innerlines...) - vertices = (f₂.vertices..., f₁.vertices...) - coupled = f₁.coupled - f′ = FusionTree(uncoupled, coupled, isdual, inner, vertices) - return fusiontreedict(I){F, T}(f′ => coeff) - else # recursive definition - N2 = length(f₂) - f₂′, f₂′′ = split(f₂, N2 - 1) - local newtrees::fusiontreedict(I){F, T} - for (f, coeff) in insertat(f₁, i, f₂′′) - for (f′, coeff′) in insertat(f, i, f₂′) - if @isdefined newtrees - coeff′′ = coeff * coeff′ - newtrees[f′] = get(newtrees, f′, zero(coeff′′)) + coeff′′ - else - newtrees = fusiontreedict(I){F, T}(f′ => coeff * coeff′) - end - end - end - return newtrees - end -end - -""" - split(f::FusionTree{I, N}, M::Int) - -> (::FusionTree{I, M}, ::FusionTree{I, N-M+1}) - -Split a fusion tree into two. The first tree has as uncoupled sectors the first `M` -uncoupled sectors of the input tree `f`, whereas its coupled sector corresponds to the -internal sector between uncoupled sectors `M` and `M+1` of the original tree `f`. The -second tree has as first uncoupled sector that same internal sector of `f`, followed by -remaining `N-M` uncoupled sectors of `f`. It couples to the same sector as `f`. This -operation is the inverse of `insertat` in the sense that if -`f₁, f₂ = split(t, M) ⇒ f == insertat(f₂, 1, f₁)`. -""" -@inline function split(f::FusionTree{I, N}, M::Int) where {I, N} - if M > N || M < 0 - throw(ArgumentError("M should be between 0 and N = $N")) - elseif M === N - (f, FusionTree{I}((f.coupled,), f.coupled, (false,), (), ())) - elseif M === 1 - isdual1 = (f.isdual[1],) - isdual2 = TupleTools.setindex(f.isdual, false, 1) - f₁ = FusionTree{I}((f.uncoupled[1],), f.uncoupled[1], isdual1, (), ()) - f₂ = FusionTree{I}(f.uncoupled, f.coupled, isdual2, f.innerlines, f.vertices) - return f₁, f₂ - elseif M === 0 - u = leftunit(f.uncoupled[1]) - f₁ = FusionTree{I}((), u, (), ()) - uncoupled2 = (u, f.uncoupled...) - coupled2 = f.coupled - isdual2 = (false, f.isdual...) - innerlines2 = N >= 2 ? (f.uncoupled[1], f.innerlines...) : () - if FusionStyle(I) isa GenericFusion - vertices2 = (1, f.vertices...) - return f₁, FusionTree{I}(uncoupled2, coupled2, isdual2, innerlines2, vertices2) - else - return f₁, FusionTree{I}(uncoupled2, coupled2, isdual2, innerlines2) - end - else - uncoupled1 = ntuple(n -> f.uncoupled[n], M) - isdual1 = ntuple(n -> f.isdual[n], M) - innerlines1 = ntuple(n -> f.innerlines[n], max(0, M - 2)) - coupled1 = f.innerlines[M - 1] - vertices1 = ntuple(n -> f.vertices[n], M - 1) - - uncoupled2 = ntuple(N - M + 1) do n - return n == 1 ? f.innerlines[M - 1] : f.uncoupled[M + n - 1] - end - isdual2 = ntuple(N - M + 1) do n - return n == 1 ? false : f.isdual[M + n - 1] - end - innerlines2 = ntuple(n -> f.innerlines[M - 1 + n], N - M - 1) - coupled2 = f.coupled - vertices2 = ntuple(n -> f.vertices[M - 1 + n], N - M) - - f₁ = FusionTree{I}(uncoupled1, coupled1, isdual1, innerlines1, vertices1) - f₂ = FusionTree{I}(uncoupled2, coupled2, isdual2, innerlines2, vertices2) - return f₁, f₂ - end -end - -""" - merge(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, c::I, μ = 1) - -> <:AbstractDict{<:FusionTree{I, N₁+N₂}, <:Number} - -Merge two fusion trees together to a linear combination of fusion trees whose uncoupled -sectors are those of `f₁` followed by those of `f₂`, and where the two coupled sectors of -`f₁` and `f₂` are further fused to `c`. In case of -`FusionStyle(I) == GenericFusion()`, also a degeneracy label `μ` for the fusion of -the coupled sectors of `f₁` and `f₂` to `c` needs to be specified. -""" -function merge(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, c::I) where {I, N₁, N₂} - if FusionStyle(I) isa GenericFusion - throw(ArgumentError("vertex label for merging required")) - end - return merge(f₁, f₂, c, 1) -end -function merge(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, c::I, μ) where {I, N₁, N₂} - if !(c in f₁.coupled ⊗ f₂.coupled) - throw(SectorMismatch("cannot fuse sectors $(f₁.coupled) and $(f₂.coupled) to $c")) - end - if μ > Nsymbol(f₁.coupled, f₂.coupled, c) - throw(ArgumentError("invalid fusion vertex label $μ")) - end - f₀ = FusionTree{I}((f₁.coupled, f₂.coupled), c, (false, false), (), (μ,)) - f, coeff = first(insertat(f₀, 1, f₁)) # takes fast path, single output - @assert coeff == one(coeff) - return insertat(f, N₁ + 1, f₂) -end -function merge(f₁::FusionTree{I, 0}, f₂::FusionTree{I, 0}, c::I, μ) where {I} - Nsymbol(f₁.coupled, f₂.coupled, c) == μ == 1 || - throw(SectorMismatch("cannot fuse sectors $(f₁.coupled) and $(f₂.coupled) to $c")) - return fusiontreedict(I)(f₁ => Fsymbol(c, c, c, c, c, c)[1, 1, 1, 1]) -end - -# ELEMENTARY DUALITY MANIPULATIONS: A- and B-moves -#--------------------------------------------------------- -# -> elementary manipulations that depend on the duality (rigidity) and pivotal structure -# -> planar manipulations that do not require braiding, everything is in Fsymbol (A/Bsymbol) -# -> B-move (bendleft, bendright) is simple in standard basis -# -> A-move (foldleft, foldright) is complicated, needs to be reexpressed in standard form - -# flip a duality flag of a fusion tree -function flip(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, i::Int; inv::Bool = false) where {I <: Sector, N₁, N₂} - @assert 0 < i ≤ N₁ + N₂ - if i ≤ N₁ - a = f₁.uncoupled[i] - χₐ = frobenius_schur_phase(a) - θₐ = twist(a) - if !inv - factor = f₁.isdual[i] ? χₐ * θₐ : one(θₐ) - else - factor = f₁.isdual[i] ? one(θₐ) : χₐ * conj(θₐ) - end - isdual′ = TupleTools.setindex(f₁.isdual, !f₁.isdual[i], i) - f₁′ = FusionTree{I}(f₁.uncoupled, f₁.coupled, isdual′, f₁.innerlines, f₁.vertices) - return SingletonDict((f₁′, f₂) => factor) - else - i -= N₁ - a = f₂.uncoupled[i] - χₐ = frobenius_schur_phase(a) - θₐ = twist(a) - if !inv - factor = f₂.isdual[i] ? χₐ * one(θₐ) : θₐ - else - factor = f₂.isdual[i] ? conj(θₐ) : χₐ * one(θₐ) - end - isdual′ = TupleTools.setindex(f₂.isdual, !f₂.isdual[i], i) - f₂′ = FusionTree{I}(f₂.uncoupled, f₂.coupled, isdual′, f₂.innerlines, f₂.vertices) - return SingletonDict((f₁, f₂′) => factor) - end -end -function flip(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, ind; inv::Bool = false) where {I <: Sector, N₁, N₂} - f₁′, f₂′ = f₁, f₂ - factor = one(sectorscalartype(I)) - for i in ind - (f₁′, f₂′), s = only(flip(f₁′, f₂′, i; inv)) - factor *= s - end - return SingletonDict((f₁′, f₂′) => factor) -end - -# change to N₁ - 1, N₂ + 1 -function bendright(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}) where {I <: Sector, N₁, N₂} - # map final splitting vertex (a, b)<-c to fusion vertex a<-(c, dual(b)) - @assert N₁ > 0 - c = f₁.coupled - a = N₁ == 1 ? leftunit(f₁.uncoupled[1]) : - (N₁ == 2 ? f₁.uncoupled[1] : f₁.innerlines[end]) - b = f₁.uncoupled[N₁] - - uncoupled1 = TupleTools.front(f₁.uncoupled) - isdual1 = TupleTools.front(f₁.isdual) - inner1 = N₁ > 2 ? TupleTools.front(f₁.innerlines) : () - vertices1 = N₁ > 1 ? TupleTools.front(f₁.vertices) : () - f₁′ = FusionTree(uncoupled1, a, isdual1, inner1, vertices1) - - uncoupled2 = (f₂.uncoupled..., dual(b)) - isdual2 = (f₂.isdual..., !(f₁.isdual[N₁])) - inner2 = N₂ > 1 ? (f₂.innerlines..., c) : () - - coeff₀ = sqrtdim(c) * invsqrtdim(a) - if f₁.isdual[N₁] - coeff₀ *= conj(frobenius_schur_phase(dual(b))) - end - if FusionStyle(I) isa MultiplicityFreeFusion - coeff = coeff₀ * Bsymbol(a, b, c) - vertices2 = N₂ > 0 ? (f₂.vertices..., 1) : () - f₂′ = FusionTree(uncoupled2, a, isdual2, inner2, vertices2) - return SingletonDict((f₁′, f₂′) => coeff) - else - local newtrees - Bmat = Bsymbol(a, b, c) - μ = N₁ > 1 ? f₁.vertices[end] : 1 - for ν in axes(Bmat, 2) - coeff = coeff₀ * Bmat[μ, ν] - iszero(coeff) && continue - vertices2 = N₂ > 0 ? (f₂.vertices..., ν) : () - f₂′ = FusionTree(uncoupled2, a, isdual2, inner2, vertices2) - if @isdefined newtrees - push!(newtrees, (f₁′, f₂′) => coeff) - else - newtrees = FusionTreeDict((f₁′, f₂′) => coeff) - end - end - return newtrees - end -end -# change to N₁ + 1, N₂ - 1 -function bendleft(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I} - # map final fusion vertex c<-(a, b) to splitting vertex (c, dual(b))<-a - return fusiontreedict(I)( - (f₁′, f₂′) => conj(coeff) for ((f₂′, f₁′), coeff) in bendright(f₂, f₁) - ) -end - -# change to N₁ - 1, N₂ + 1 -function foldright(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}) where {I <: Sector, N₁, N₂} - # map first splitting vertex (a, b)<-c to fusion vertex b<-(dual(a), c) - @assert N₁ > 0 - a = f₁.uncoupled[1] - isduala = f₁.isdual[1] - factor = sqrtdim(a) - if !isduala - factor *= conj(frobenius_schur_phase(a)) - end - c1 = dual(a) - c2 = f₁.coupled - uncoupled = Base.tail(f₁.uncoupled) - isdual = Base.tail(f₁.isdual) - if FusionStyle(I) isa UniqueFusion - c = first(c1 ⊗ c2) - fl = FusionTree{I}(Base.tail(f₁.uncoupled), c, Base.tail(f₁.isdual)) - fr = FusionTree{I}((c1, f₂.uncoupled...), c, (!isduala, f₂.isdual...)) - return fusiontreedict(I)((fl, fr) => factor) - else - hasmultiplicities = FusionStyle(a) isa GenericFusion - local newtrees - if N₁ == 1 - cset = (leftunit(c1),) # or rightunit(a) - elseif N₁ == 2 - cset = (f₁.uncoupled[2],) - else - cset = ⊗(Base.tail(f₁.uncoupled)...) - end - for c in c1 ⊗ c2 - c ∈ cset || continue - for μ in 1:Nsymbol(c1, c2, c) - fc = FusionTree((c1, c2), c, (!isduala, false), (), (μ,)) - for (fl′, coeff1) in insertat(fc, 2, f₁) - N₁ > 1 && !isunit(fl′.innerlines[1]) && continue - coupled = fl′.coupled - uncoupled = Base.tail(Base.tail(fl′.uncoupled)) - isdual = Base.tail(Base.tail(fl′.isdual)) - inner = N₁ <= 3 ? () : Base.tail(Base.tail(fl′.innerlines)) - vertices = N₁ <= 2 ? () : Base.tail(Base.tail(fl′.vertices)) - fl = FusionTree{I}(uncoupled, coupled, isdual, inner, vertices) - for (fr, coeff2) in insertat(fc, 2, f₂) - coeff = factor * coeff1 * conj(coeff2) - if (@isdefined newtrees) - newtrees[(fl, fr)] = get(newtrees, (fl, fr), zero(coeff)) + - coeff - else - newtrees = fusiontreedict(I)((fl, fr) => coeff) - end - end - end - end - end - return newtrees - end -end -# change to N₁ + 1, N₂ - 1 -function foldleft(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I} - # map first fusion vertex c<-(a, b) to splitting vertex (dual(a), c)<-b - return fusiontreedict(I)( - (f₁′, f₂′) => conj(coeff) for ((f₂′, f₁′), coeff) in foldright(f₂, f₁) - ) -end - -# COMPOSITE DUALITY MANIPULATIONS PART 1: Repartition and transpose -#------------------------------------------------------------------- -# -> composite manipulations that depend on the duality (rigidity) and pivotal structure -# -> planar manipulations that do not require braiding, everything is in Fsymbol (A/Bsymbol) -# -> transpose expressed as cyclic permutation -# one-argument version: check whether `p` is a cyclic permutation (of `1:length(p)`) -function iscyclicpermutation(p) - N = length(p) - @inbounds for i in 1:N - p[mod1(i + 1, N)] == mod1(p[i] + 1, N) || return false - end - return true -end -# two-argument version: check whether `v1` is a cyclic permutation of `v2` -function iscyclicpermutation(v1, v2) - length(v1) == length(v2) || return false - return iscyclicpermutation(indexin(v1, v2)) -end - -# clockwise cyclic permutation while preserving (N₁, N₂): foldright & bendleft -function cycleclockwise(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I <: Sector} - local newtrees - if length(f₁) > 0 - for ((f1a, f2a), coeffa) in foldright(f₁, f₂) - for ((f1b, f2b), coeffb) in bendleft(f1a, f2a) - coeff = coeffa * coeffb - if (@isdefined newtrees) - newtrees[(f1b, f2b)] = get(newtrees, (f1b, f2b), zero(coeff)) + coeff - else - newtrees = fusiontreedict(I)((f1b, f2b) => coeff) - end - end - end - else - for ((f1a, f2a), coeffa) in bendleft(f₁, f₂) - for ((f1b, f2b), coeffb) in foldright(f1a, f2a) - coeff = coeffa * coeffb - if (@isdefined newtrees) - newtrees[(f1b, f2b)] = get(newtrees, (f1b, f2b), zero(coeff)) + coeff - else - newtrees = fusiontreedict(I)((f1b, f2b) => coeff) - end - end - end - end - return newtrees -end - -# anticlockwise cyclic permutation while preserving (N₁, N₂): foldleft & bendright -function cycleanticlockwise(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I <: Sector} - local newtrees - if length(f₂) > 0 - for ((f1a, f2a), coeffa) in foldleft(f₁, f₂) - for ((f1b, f2b), coeffb) in bendright(f1a, f2a) - coeff = coeffa * coeffb - if (@isdefined newtrees) - newtrees[(f1b, f2b)] = get(newtrees, (f1b, f2b), zero(coeff)) + coeff - else - newtrees = fusiontreedict(I)((f1b, f2b) => coeff) - end - end - end - else - for ((f1a, f2a), coeffa) in bendright(f₁, f₂) - for ((f1b, f2b), coeffb) in foldleft(f1a, f2a) - coeff = coeffa * coeffb - if (@isdefined newtrees) - newtrees[(f1b, f2b)] = get(newtrees, (f1b, f2b), zero(coeff)) + coeff - else - newtrees = fusiontreedict(I)((f1b, f2b) => coeff) - end - end - end - end - return newtrees -end - -# repartition double fusion tree -""" - repartition(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, N::Int) where {I, N₁, N₂} - -> <:AbstractDict{Tuple{FusionTree{I, N}, FusionTree{I, N₁+N₂-N}}, <:Number} - -Input is a double fusion tree that describes the fusion of a set of incoming uncoupled -sectors to a set of outgoing uncoupled sectors, represented using the individual trees of -outgoing (`f₁`) and incoming sectors (`f₂`) respectively (with identical coupled sector -`f₁.coupled == f₂.coupled`). Computes new trees and corresponding coefficients obtained from -repartitioning the tree by bending incoming to outgoing sectors (or vice versa) in order to -have `N` outgoing sectors. -""" -@inline function repartition( - f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, N::Int - ) where {I <: Sector, N₁, N₂} - f₁.coupled == f₂.coupled || throw(SectorMismatch()) - @assert 0 <= N <= N₁ + N₂ - return _recursive_repartition(f₁, f₂, Val(N)) -end - -function _recursive_repartition( - f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}, ::Val{N} - ) where {I <: Sector, N₁, N₂, N} - # recursive definition is only way to get correct number of loops for - # GenericFusion, but is too complex for type inference to handle, so we - # precompute the parameters of the return type - F₁ = fusiontreetype(I, N) - F₂ = fusiontreetype(I, N₁ + N₂ - N) - T = sectorscalartype(I) - coeff = one(T) - if N == N₁ - return fusiontreedict(I){Tuple{F₁, F₂}, T}((f₁, f₂) => coeff) - else - local newtrees::fusiontreedict(I){Tuple{F₁, F₂}, T} - for ((f₁′, f₂′), coeff1) in (N < N₁ ? bendright(f₁, f₂) : bendleft(f₁, f₂)) - for ((f₁′′, f₂′′), coeff2) in _recursive_repartition(f₁′, f₂′, Val(N)) - if (@isdefined newtrees) - push!(newtrees, (f₁′′, f₂′′) => coeff1 * coeff2) - else - newtrees = fusiontreedict(I){Tuple{F₁, F₂}, T}( - (f₁′′, f₂′′) => coeff1 * coeff2 - ) - end - end - end - return newtrees - end -end - -""" - transpose(f₁::FusionTree{I}, f₂::FusionTree{I}, - p1::NTuple{N₁, Int}, p2::NTuple{N₂, Int}) where {I, N₁, N₂} - -> <:AbstractDict{Tuple{FusionTree{I, N₁}, FusionTree{I, N₂}}, <:Number} - -Input is a double fusion tree that describes the fusion of a set of incoming uncoupled -sectors to a set of outgoing uncoupled sectors, represented using the individual trees of -outgoing (`t1`) and incoming sectors (`t2`) respectively (with identical coupled sector -`t1.coupled == t2.coupled`). Computes new trees and corresponding coefficients obtained from -repartitioning and permuting the tree such that sectors `p1` become outgoing and sectors -`p2` become incoming. -""" -function Base.transpose( - f₁::FusionTree{I}, f₂::FusionTree{I}, p1::IndexTuple{N₁}, p2::IndexTuple{N₂} - ) where {I <: Sector, N₁, N₂} - N = N₁ + N₂ - @assert length(f₁) + length(f₂) == N - p = linearizepermutation(p1, p2, length(f₁), length(f₂)) - @assert iscyclicpermutation(p) - return fstranspose((f₁, f₂, p1, p2)) -end - -const FSTransposeKey{I <: Sector, N₁, N₂} = Tuple{ - <:FusionTree{I}, <:FusionTree{I}, - IndexTuple{N₁}, IndexTuple{N₂}, -} - -function _fsdicttype(I, N₁, N₂) - F₁ = fusiontreetype(I, N₁) - F₂ = fusiontreetype(I, N₂) - T = sectorscalartype(I) - return fusiontreedict(I){Tuple{F₁, F₂}, T} -end - -@cached function fstranspose(key::FSTransposeKey{I, N₁, N₂})::_fsdicttype(I, N₁, N₂) where {I <: Sector, N₁, N₂} - f₁, f₂, p1, p2 = key - N = N₁ + N₂ - p = linearizepermutation(p1, p2, length(f₁), length(f₂)) - newtrees = repartition(f₁, f₂, N₁) - length(p) == 0 && return newtrees - i1 = findfirst(==(1), p) - @assert i1 !== nothing - i1 == 1 && return newtrees - Nhalf = N >> 1 - while 1 < i1 <= Nhalf - local newtrees′ - for ((f1a, f2a), coeffa) in newtrees - for ((f1b, f2b), coeffb) in cycleanticlockwise(f1a, f2a) - coeff = coeffa * coeffb - if (@isdefined newtrees′) - newtrees′[(f1b, f2b)] = get(newtrees′, (f1b, f2b), zero(coeff)) + coeff - else - newtrees′ = fusiontreedict(I)((f1b, f2b) => coeff) - end - end - end - newtrees = newtrees′ - i1 -= 1 - end - while Nhalf < i1 - local newtrees′ - for ((f1a, f2a), coeffa) in newtrees - for ((f1b, f2b), coeffb) in cycleclockwise(f1a, f2a) - coeff = coeffa * coeffb - if (@isdefined newtrees′) - newtrees′[(f1b, f2b)] = get(newtrees′, (f1b, f2b), zero(coeff)) + coeff - else - newtrees′ = fusiontreedict(I)((f1b, f2b) => coeff) - end - end - end - newtrees = newtrees′ - i1 = mod1(i1 + 1, N) - end - return newtrees -end - -function CacheStyle(::typeof(fstranspose), k::FSTransposeKey{I}) where {I <: Sector} - if FusionStyle(I) isa UniqueFusion - return NoCache() - else - return GlobalLRUCache() - end -end - -# COMPOSITE DUALITY MANIPULATIONS PART 2: Planar traces -#------------------------------------------------------------------- -# -> composite manipulations that depend on the duality (rigidity) and pivotal structure -# -> planar manipulations that do not require braiding, everything is in Fsymbol (A/Bsymbol) - -function planar_trace( - f₁::FusionTree{I}, f₂::FusionTree{I}, - p1::IndexTuple{N₁}, p2::IndexTuple{N₂}, - q1::IndexTuple{N₃}, q2::IndexTuple{N₃} - ) where {I <: Sector, N₁, N₂, N₃} - N = N₁ + N₂ + 2N₃ - @assert length(f₁) + length(f₂) == N - if N₃ == 0 - return transpose(f₁, f₂, p1, p2) - end - - linearindex = ( - ntuple(identity, Val(length(f₁)))..., - reverse(length(f₁) .+ ntuple(identity, Val(length(f₂))))..., - ) - - q1′ = TupleTools.getindices(linearindex, q1) - q2′ = TupleTools.getindices(linearindex, q2) - p1′, p2′ = let q′ = (q1′..., q2′...) - ( - map(l -> l - count(l .> q′), TupleTools.getindices(linearindex, p1)), - map(l -> l - count(l .> q′), TupleTools.getindices(linearindex, p2)), - ) - end - - T = sectorscalartype(I) - F₁ = fusiontreetype(I, N₁) - F₂ = fusiontreetype(I, N₂) - newtrees = FusionTreeDict{Tuple{F₁, F₂}, T}() - for ((f₁′, f₂′), coeff′) in repartition(f₁, f₂, N) - for (f₁′′, coeff′′) in planar_trace(f₁′, q1′, q2′) - for (f12′′′, coeff′′′) in transpose(f₁′′, f₂′, p1′, p2′) - coeff = coeff′ * coeff′′ * coeff′′′ - if !iszero(coeff) - newtrees[f12′′′] = get(newtrees, f12′′′, zero(coeff)) + coeff - end - end - end - end - return newtrees -end - -""" - planar_trace(f::FusionTree{I,N}, q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {I<:Sector,N,N₃} - -> <:AbstractDict{FusionTree{I,N-2*N₃}, <:Number} - -Perform a planar trace of the uncoupled indices of the fusion tree `f` at `q1` with those at -`q2`, where `q1[i]` is connected to `q2[i]` for all `i`. The result is returned as a dictionary -of output trees and corresponding coefficients. -""" -function planar_trace( - f::FusionTree{I, N}, q1::IndexTuple{N₃}, q2::IndexTuple{N₃} - ) where {I <: Sector, N, N₃} - T = sectorscalartype(I) - F = fusiontreetype(I, N - 2 * N₃) - newtrees = FusionTreeDict{F, T}() - N₃ === 0 && return push!(newtrees, f => one(T)) - - for (i, j) in zip(q1, q2) - (f.uncoupled[i] == dual(f.uncoupled[j]) && f.isdual[i] != f.isdual[j]) || - return newtrees - end - k = 1 - local i, j - while k <= N₃ - if mod1(q1[k] + 1, N) == q2[k] - i = q1[k] - j = q2[k] - break - elseif mod1(q2[k] + 1, N) == q1[k] - i = q2[k] - j = q1[k] - break - else - k += 1 - end - end - k > N₃ && throw(ArgumentError("Not a planar trace")) - - q1′ = let i = i, j = j - map(l -> (l - (l > i) - (l > j)), TupleTools.deleteat(q1, k)) - end - q2′ = let i = i, j = j - map(l -> (l - (l > i) - (l > j)), TupleTools.deleteat(q2, k)) - end - for (f′, coeff′) in elementary_trace(f, i) - for (f′′, coeff′′) in planar_trace(f′, q1′, q2′) - coeff = coeff′ * coeff′′ - if !iszero(coeff) - newtrees[f′′] = get(newtrees, f′′, zero(coeff)) + coeff - end - end - end - return newtrees -end - -# trace two neighbouring indices of a single fusion tree -""" - elementary_trace(f::FusionTree{I,N}, i) where {I<:Sector,N} -> <:AbstractDict{FusionTree{I,N-2}, <:Number} - -Perform an elementary trace of neighbouring uncoupled indices `i` and -`i+1` on a fusion tree `f`, and returns the result as a dictionary of output trees and -corresponding coefficients. -""" -function elementary_trace(f::FusionTree{I, N}, i) where {I <: Sector, N} - (N > 1 && 1 <= i <= N) || - throw(ArgumentError("Cannot trace outputs i=$i and i+1 out of only $N outputs")) - i < N || isunit(f.coupled) || - throw(ArgumentError("Cannot trace outputs i=$N and 1 of fusion tree that couples to non-trivial sector")) - - T = sectorscalartype(I) - F = fusiontreetype(I, N - 2) - newtrees = FusionTreeDict{F, T}() - - j = mod1(i + 1, N) - b = f.uncoupled[i] - b′ = f.uncoupled[j] - # if trace is zero, return empty dict - (b == dual(b′) && f.isdual[i] != f.isdual[j]) || return newtrees - if i < N - inner_extended = (leftunit(f.uncoupled[1]), f.uncoupled[1], f.innerlines..., f.coupled) - a = inner_extended[i] - d = inner_extended[i + 2] - a == d || return newtrees - uncoupled′ = TupleTools.deleteat(TupleTools.deleteat(f.uncoupled, i + 1), i) - isdual′ = TupleTools.deleteat(TupleTools.deleteat(f.isdual, i + 1), i) - coupled′ = f.coupled - if N <= 4 - inner′ = () - else - inner′ = i <= 2 ? Base.tail(Base.tail(f.innerlines)) : - TupleTools.deleteat(TupleTools.deleteat(f.innerlines, i - 1), i - 2) - end - if N <= 3 - vertices′ = () - else - vertices′ = i <= 2 ? Base.tail(Base.tail(f.vertices)) : - TupleTools.deleteat(TupleTools.deleteat(f.vertices, i), i - 1) - end - f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′) - coeff = sqrtdim(b) - if i > 1 - c = f.innerlines[i - 1] - if FusionStyle(I) isa MultiplicityFreeFusion - coeff *= Fsymbol(a, b, dual(b), a, c, rightunit(a)) - else - μ = f.vertices[i - 1] - ν = f.vertices[i] - coeff *= Fsymbol(a, b, dual(b), a, c, rightunit(a))[μ, ν, 1, 1] - end - end - if f.isdual[i] - coeff *= frobenius_schur_phase(b) - end - push!(newtrees, f′ => coeff) - return newtrees - else # i == N - unit = leftunit(b) - if N == 2 - f′ = FusionTree{I}((), unit, (), (), ()) - coeff = sqrtdim(b) - if !(f.isdual[N]) - coeff *= conj(frobenius_schur_phase(b)) - end - push!(newtrees, f′ => coeff) - return newtrees - end - uncoupled_ = TupleTools.front(f.uncoupled) - inner_ = TupleTools.front(f.innerlines) - coupled_ = f.innerlines[end] - isdual_ = TupleTools.front(f.isdual) - vertices_ = TupleTools.front(f.vertices) - f_ = FusionTree(uncoupled_, coupled_, isdual_, inner_, vertices_) - fs = FusionTree((b,), b, (!f.isdual[1],), (), ()) - for (f_′, coeff) in merge(fs, f_, unit, 1) - f_′.innerlines[1] == unit || continue - uncoupled′ = Base.tail(Base.tail(f_′.uncoupled)) - isdual′ = Base.tail(Base.tail(f_′.isdual)) - inner′ = N <= 4 ? () : Base.tail(Base.tail(f_′.innerlines)) - vertices′ = N <= 3 ? () : Base.tail(Base.tail(f_′.vertices)) - f′ = FusionTree(uncoupled′, unit, isdual′, inner′, vertices′) - coeff *= sqrtdim(b) - if !(f.isdual[N]) - coeff *= conj(frobenius_schur_phase(b)) - end - newtrees[f′] = get(newtrees, f′, zero(coeff)) + coeff - end - return newtrees - end -end - -# BRAIDING MANIPULATIONS: -#----------------------------------------------- -# -> manipulations that depend on a braiding -# -> requires both Fsymbol and Rsymbol -""" - artin_braid(f::FusionTree, i; inv::Bool = false) -> <:AbstractDict{typeof(f), <:Number} - -Perform an elementary braid (Artin generator) of neighbouring uncoupled indices `i` and -`i+1` on a fusion tree `f`, and returns the result as a dictionary of output trees and -corresponding coefficients. - -The keyword `inv` determines whether index `i` will braid above or below index `i+1`, i.e. -applying `artin_braid(f′, i; inv = true)` to all the outputs `f′` of -`artin_braid(f, i; inv = false)` and collecting the results should yield a single fusion -tree with non-zero coefficient, namely `f` with coefficient `1`. This keyword has no effect -if `BraidingStyle(sectortype(f)) isa SymmetricBraiding`. -""" -function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I <: Sector, N} - 1 <= i < N || - throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs")) - uncoupled = f.uncoupled - a, b = uncoupled[i], uncoupled[i + 1] - uncoupled′ = TupleTools.setindex(uncoupled, b, i) - uncoupled′ = TupleTools.setindex(uncoupled′, a, i + 1) - coupled′ = f.coupled - isdual′ = TupleTools.setindex(f.isdual, f.isdual[i], i + 1) - isdual′ = TupleTools.setindex(isdual′, f.isdual[i + 1], i) - inner = f.innerlines - inner_extended = (uncoupled[1], inner..., coupled′) - vertices = f.vertices - oneT = one(sectorscalartype(I)) - - if isunit(a) || isunit(b) - # braiding with trivial sector: simple and always possible - inner′ = inner - vertices′ = vertices - if i > 1 # we also need to alter innerlines and vertices - inner′ = TupleTools.setindex( - inner, - inner_extended[isunit(a) ? (i + 1) : (i - 1)], - i - 1 - ) - vertices′ = TupleTools.setindex(vertices′, vertices[i], i - 1) - vertices′ = TupleTools.setindex(vertices′, vertices[i - 1], i) - end - f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′) - return fusiontreedict(I)(f′ => oneT) - end - - BraidingStyle(I) isa NoBraiding && - throw(SectorMismatch("Cannot braid sectors $(uncoupled[i]) and $(uncoupled[i + 1])")) - - if i == 1 - c = N > 2 ? inner[1] : coupled′ - if FusionStyle(I) isa MultiplicityFreeFusion - R = oftype(oneT, (inv ? conj(Rsymbol(b, a, c)) : Rsymbol(a, b, c))) - f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices) - return fusiontreedict(I)(f′ => R) - else # GenericFusion - μ = vertices[1] - Rmat = inv ? Rsymbol(b, a, c)' : Rsymbol(a, b, c) - local newtrees - for ν in axes(Rmat, 2) - R = oftype(oneT, Rmat[μ, ν]) - iszero(R) && continue - vertices′ = TupleTools.setindex(vertices, ν, 1) - f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices′) - if (@isdefined newtrees) - push!(newtrees, f′ => R) - else - newtrees = fusiontreedict(I)(f′ => R) - end - end - return newtrees - end - end - # case i > 1: other naming convention - b = uncoupled[i] - d = uncoupled[i + 1] - a = inner_extended[i - 1] - c = inner_extended[i] - e = inner_extended[i + 1] - if FusionStyle(I) isa UniqueFusion - c′ = first(a ⊗ d) - coeff = oftype( - oneT, - if inv - conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′) - else - Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′)) - end - ) - inner′ = TupleTools.setindex(inner, c′, i - 1) - f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′) - return fusiontreedict(I)(f′ => coeff) - elseif FusionStyle(I) isa SimpleFusion - local newtrees - for c′ in intersect(a ⊗ d, e ⊗ dual(b)) - coeff = oftype( - oneT, - if inv - conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′) - else - Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′)) - end - ) - iszero(coeff) && continue - inner′ = TupleTools.setindex(inner, c′, i - 1) - f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′) - if (@isdefined newtrees) - push!(newtrees, f′ => coeff) - else - newtrees = fusiontreedict(I)(f′ => coeff) - end - end - return newtrees - else # GenericFusion - local newtrees - for c′ in intersect(a ⊗ d, e ⊗ dual(b)) - Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e) - Rmat2 = inv ? Rsymbol(d, a, c′)' : Rsymbol(a, d, c′) - Fmat = Fsymbol(d, a, b, e, c′, c) - μ = vertices[i - 1] - ν = vertices[i] - for σ in 1:Nsymbol(a, d, c′) - for λ in 1:Nsymbol(c′, b, e) - coeff = zero(oneT) - for ρ in 1:Nsymbol(d, c, e), κ in 1:Nsymbol(d, a, c′) - coeff += Rmat1[ν, ρ] * conj(Fmat[κ, λ, μ, ρ]) * conj(Rmat2[σ, κ]) - end - iszero(coeff) && continue - vertices′ = TupleTools.setindex(vertices, σ, i - 1) - vertices′ = TupleTools.setindex(vertices′, λ, i) - inner′ = TupleTools.setindex(inner, c′, i - 1) - f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′) - if (@isdefined newtrees) - push!(newtrees, f′ => coeff) - else - newtrees = fusiontreedict(I)(f′ => coeff) - end - end - end - end - return newtrees - end -end - -# braid fusion tree -""" - braid(f::FusionTree{<:Sector, N}, levels::NTuple{N, Int}, p::NTuple{N, Int}) - -> <:AbstractDict{typeof(t), <:Number} - -Perform a braiding of the uncoupled indices of the fusion tree `f` and return the result as -a `<:AbstractDict` of output trees and corresponding coefficients. The braiding is -determined by specifying that the new sector at position `k` corresponds to the sector that -was originally at the position `i = p[k]`, and assigning to every index `i` of the original -fusion tree a distinct level or depth `levels[i]`. This permutation is then decomposed into -elementary swaps between neighbouring indices, where the swaps are applied as braids such -that if `i` and `j` cross, ``τ_{i,j}`` is applied if `levels[i] < levels[j]` and -``τ_{j,i}^{-1}`` if `levels[i] > levels[j]`. This does not allow to encode the most general -braid, but a general braid can be obtained by combining such operations. -""" -function braid( - f::FusionTree{I, N}, levels::NTuple{N, Int}, p::NTuple{N, Int} - ) where {I <: Sector, N} - TupleTools.isperm(p) || throw(ArgumentError("not a valid permutation: $p")) - if FusionStyle(I) isa UniqueFusion && BraidingStyle(I) isa SymmetricBraiding - coeff = one(sectorscalartype(I)) - for i in 1:N - for j in 1:(i - 1) - if p[j] > p[i] - a, b = f.uncoupled[p[j]], f.uncoupled[p[i]] - coeff *= Rsymbol(a, b, first(a ⊗ b)) - end - end - end - uncoupled′ = TupleTools._permute(f.uncoupled, p) - coupled′ = f.coupled - isdual′ = TupleTools._permute(f.isdual, p) - f′ = FusionTree{I}(uncoupled′, coupled′, isdual′) - return fusiontreedict(I)(f′ => coeff) - else - T = sectorscalartype(I) - coeff = one(T) - trees = FusionTreeDict(f => coeff) - newtrees = empty(trees) - for s in permutation2swaps(p) - inv = levels[s] > levels[s + 1] - for (f, c) in trees - for (f′, c′) in artin_braid(f, s; inv) - newtrees[f′] = get(newtrees, f′, zero(coeff)) + c * c′ - end - end - l = levels[s] - levels = TupleTools.setindex(levels, levels[s + 1], s) - levels = TupleTools.setindex(levels, l, s + 1) - trees, newtrees = newtrees, trees - empty!(newtrees) - end - return trees - end -end - -# permute fusion tree -""" - permute(f::FusionTree, p::NTuple{N, Int}) -> <:AbstractDict{typeof(t), <:Number} - -Perform a permutation of the uncoupled indices of the fusion tree `f` and returns the result -as a `<:AbstractDict` of output trees and corresponding coefficients; this requires that -`BraidingStyle(sectortype(f)) isa SymmetricBraiding`. -""" -function permute(f::FusionTree{I, N}, p::NTuple{N, Int}) where {I <: Sector, N} - @assert BraidingStyle(I) isa SymmetricBraiding - return braid(f, ntuple(identity, Val(N)), p) -end - -# braid double fusion tree -""" - braid(f₁::FusionTree{I}, f₂::FusionTree{I}, - levels1::IndexTuple, levels2::IndexTuple, - p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {I<:Sector, N₁, N₂} - -> <:AbstractDict{Tuple{FusionTree{I, N₁}, FusionTree{I, N₂}}, <:Number} - -Input is a fusion-splitting tree pair that describes the fusion of a set of incoming -uncoupled sectors to a set of outgoing uncoupled sectors, represented using the splitting -tree `f₁` and fusion tree `f₂`, such that the incoming sectors `f₂.uncoupled` are fused to -`f₁.coupled == f₂.coupled` and then to the outgoing sectors `f₁.uncoupled`. Compute new -trees and corresponding coefficients obtained from repartitioning and braiding the tree such -that sectors `p1` become outgoing and sectors `p2` become incoming. The uncoupled indices in -splitting tree `f₁` and fusion tree `f₂` have levels (or depths) `levels1` and `levels2` -respectively, which determines how indices braid. In particular, if `i` and `j` cross, -``τ_{i,j}`` is applied if `levels[i] < levels[j]` and ``τ_{j,i}^{-1}`` if `levels[i] > -levels[j]`. This does not allow to encode the most general braid, but a general braid can -be obtained by combining such operations. -""" -function braid( - f₁::FusionTree{I}, f₂::FusionTree{I}, - levels1::IndexTuple, levels2::IndexTuple, - p1::IndexTuple{N₁}, p2::IndexTuple{N₂} - ) where {I <: Sector, N₁, N₂} - @assert length(f₁) + length(f₂) == N₁ + N₂ - @assert length(f₁) == length(levels1) && length(f₂) == length(levels2) - @assert TupleTools.isperm((p1..., p2...)) - return fsbraid((f₁, f₂, levels1, levels2, p1, p2)) -end -const FSBraidKey{I <: Sector, N₁, N₂} = Tuple{ - <:FusionTree{I}, <:FusionTree{I}, - IndexTuple, IndexTuple, - IndexTuple{N₁}, IndexTuple{N₂}, -} - -@cached function fsbraid(key::FSBraidKey{I, N₁, N₂})::_fsdicttype(I, N₁, N₂) where {I <: Sector, N₁, N₂} - (f₁, f₂, l1, l2, p1, p2) = key - p = linearizepermutation(p1, p2, length(f₁), length(f₂)) - levels = (l1..., reverse(l2)...) - local newtrees - for ((f, f0), coeff1) in repartition(f₁, f₂, N₁ + N₂) - for (f′, coeff2) in braid(f, levels, p) - for ((f₁′, f₂′), coeff3) in repartition(f′, f0, N₁) - if @isdefined newtrees - newtrees[(f₁′, f₂′)] = get(newtrees, (f₁′, f₂′), zero(coeff3)) + - coeff1 * coeff2 * coeff3 - else - newtrees = fusiontreedict(I)((f₁′, f₂′) => coeff1 * coeff2 * coeff3) - end - end - end - end - return newtrees -end - -function CacheStyle(::typeof(fsbraid), k::FSBraidKey{I}) where {I <: Sector} - if FusionStyle(I) isa UniqueFusion - return NoCache() - else - return GlobalLRUCache() - end -end - -""" - permute(f₁::FusionTree{I}, f₂::FusionTree{I}, - p1::NTuple{N₁, Int}, p2::NTuple{N₂, Int}) where {I, N₁, N₂} - -> <:AbstractDict{Tuple{FusionTree{I, N₁}, FusionTree{I, N₂}}, <:Number} - -Input is a double fusion tree that describes the fusion of a set of incoming uncoupled -sectors to a set of outgoing uncoupled sectors, represented using the individual trees of -outgoing (`t1`) and incoming sectors (`t2`) respectively (with identical coupled sector -`t1.coupled == t2.coupled`). Computes new trees and corresponding coefficients obtained from -repartitioning and permuting the tree such that sectors `p1` become outgoing and sectors -`p2` become incoming. -""" -function permute( - f₁::FusionTree{I}, f₂::FusionTree{I}, - p1::IndexTuple{N₁}, p2::IndexTuple{N₂} - ) where {I <: Sector, N₁, N₂} - @assert BraidingStyle(I) isa SymmetricBraiding - levels1 = ntuple(identity, length(f₁)) - levels2 = length(f₁) .+ ntuple(identity, length(f₂)) - return braid(f₁, f₂, levels1, levels2, p1, p2) -end diff --git a/src/planar/planaroperations.jl b/src/planar/planaroperations.jl index a06f4431d..b1cb7d438 100644 --- a/src/planar/planaroperations.jl +++ b/src/planar/planaroperations.jl @@ -99,9 +99,11 @@ function planartrace!( end β′ = One() for (f₁, f₂) in fusiontrees(A) - for ((f₁′, f₂′), coeff) in planar_trace(f₁, f₂, p₁, p₂, q₁, q₂) + for ((f₁′, f₂′), coeff) in planar_trace((f₁, f₂), (p₁, p₂), (q₁, q₂)) TO.tensortrace!( - C[f₁′, f₂′], A[f₁, f₂], (p₁, p₂), (q₁, q₂), false, α * coeff, β′, + C[f₁′, f₂′], + A[f₁, f₂], (p₁, p₂), (q₁, q₂), false, + α * coeff, β′, backend, allocator ) end diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index 27df834c4..7a70ee3f1 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -134,6 +134,25 @@ Return the fusiontrees corresponding to all valid fusion channels of a given `Ho """ fusiontrees(W::HomSpace) = fusionblockstructure(W).fusiontreelist +""" + fusionblocks(W::HomSpace) + +Return the fusionblocks corresponding to all valid fusion channels of a given `HomSpace`, +grouped by their uncoupled charges. +""" +function fusionblocks(W::HomSpace) + I = sectortype(W) + N₁, N₂ = numout(W), numin(W) + isdual_src = (map(isdual, codomain(W).spaces), map(isdual, domain(W).spaces)) + fblocks = Vector{FusionTreeBlock{I, N₁, N₂, fusiontreetype(I, N₁, N₂)}}() + for cod_uncoupled_src in sectors(codomain(W)), dom_uncoupled_src in sectors(domain(W)) + fs_src = FusionTreeBlock{I}((cod_uncoupled_src, dom_uncoupled_src), isdual_src) + trees_src = fusiontrees(fs_src) + isempty(trees_src) || push!(fblocks, fs_src) + end + return fblocks +end + # Operations on HomSpaces # ----------------------- """ diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl index 393bd905e..b5b613177 100644 --- a/src/tensors/abstracttensor.jl +++ b/src/tensors/abstracttensor.jl @@ -269,6 +269,8 @@ end return (f₁, f₂) end +fusionblocks(t::AbstractTensorMap) = fusionblocks(space(t)) + # tensor data: block access #--------------------------- @doc """ diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index 9c92c67a1..f9e78a965 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -239,7 +239,7 @@ function planarcontract!( inv_braid = τ_levels[cindA[1]] > τ_levels[cindA[2]] for (f₁, f₂) in fusiontrees(B) local newtrees - for ((f₁′, f₂′), coeff′) in transpose(f₁, f₂, cindB, oindB) + for ((f₁′, f₂′), coeff′) in transpose((f₁, f₂), (cindB, oindB)) for (f₁′′, coeff′′) in artin_braid(f₁′, 1; inv = inv_braid) f12 = (f₁′′, f₂′) coeff = coeff′ * coeff′′ @@ -294,7 +294,7 @@ function planarcontract!( for (f₁, f₂) in fusiontrees(A) local newtrees - for ((f₁′, f₂′), coeff′) in transpose(f₁, f₂, oindA, cindA) + for ((f₁′, f₂′), coeff′) in transpose((f₁, f₂), (oindA, cindA)) for (f₂′′, coeff′′) in artin_braid(f₂′, 1; inv = inv_braid) f12 = (f₁′, f₂′′) coeff = coeff′ * conj(coeff′′) diff --git a/src/tensors/diagonal.jl b/src/tensors/diagonal.jl index 904e40cba..9f71272d5 100644 --- a/src/tensors/diagonal.jl +++ b/src/tensors/diagonal.jl @@ -194,23 +194,31 @@ end function permute( d::DiagonalTensorMap, (p₁, p₂)::Index2Tuple{1, 1}; copy::Bool = false ) - if p₁ === (1,) && p₂ === (2,) - return copy ? Base.copy(d) : d - elseif p₁ === (2,) && p₂ === (1,) # transpose - if has_shared_permute(d, (p₁, p₂)) # tranpose for bosonic sectors - return DiagonalTensorMap(copy ? Base.copy(d.data) : d.data, dual(d.domain)) - end - d′ = typeof(d)(undef, dual(d.domain)) + # special cases handled first + p₁ === (1,) && p₂ === (2,) && return copy ? Base.copy(d) : d + (p₁, p₂) === ((2,), (1,)) || # has to be transpose + throw(ArgumentError("invalid permutation $((p₁, p₂)) for tensor in space $(space(d))")) + has_shared_permute(d, (p₁, p₂)) && # tranpose for bosonic sectors + return DiagonalTensorMap(copy ? Base.copy(d.data) : d.data, dual(d.domain)) + + d′ = typeof(d)(undef, dual(d.domain)) + if FusionStyle(sectortype(d)) === UniqueFusion() for (c, b) in blocks(d) f = only(fusiontrees(codomain(d), c)) - ((f′, _), coeff) = only(permute(f, f, p₁, p₂)) + (f′, _), coeff = permute((f, f), (p₁, p₂)) c′ = f′.coupled scale!(block(d′, c′), b, coeff) end - return d′ else - throw(ArgumentError("invalid permutation $((p₁, p₂)) for tensor in space $(space(d))")) + for src in fusionblocks(d) + dst, U = permute(src, (p₁, p₂)) + c = only(fusiontrees(src))[1].coupled + c′ = only(fusiontrees(dst))[1].coupled + scale!(block(d′, c′), block(d, c), only(U)) + end end + + return d′ end # VectorInterface diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl index 408de09f0..92f77af14 100644 --- a/src/tensors/indexmanipulations.jl +++ b/src/tensors/indexmanipulations.jl @@ -17,7 +17,7 @@ function flip(t::AbstractTensorMap, I; inv::Bool = false) P = flip(space(t), I) t′ = similar(t, P) for (f₁, f₂) in fusiontrees(t) - (f₁′, f₂′), factor = only(flip(f₁, f₂, I; inv)) + (f₁′, f₂′), factor = only(flip((f₁, f₂), I; inv)) scale!(t′[f₁′, f₂′], t[f₁, f₂], factor) end return t′ @@ -482,22 +482,13 @@ function add_transform!( else I = sectortype(tdst) if I === Trivial - _add_trivial_kernel!(tdst, tsrc, p, transformer, α, β, backend...) - elseif FusionStyle(I) === UniqueFusion() - if use_threaded_transform(tdst, transformer) - _add_abelian_kernel_threaded!(tdst, tsrc, p, transformer, α, β, backend...) - else - _add_abelian_kernel_nonthreaded!( - tdst, tsrc, p, transformer, α, β, backend... - ) - end + add_trivial_kernel!(tdst, tsrc, p, transformer, α, β, backend...) else + style = FusionStyle(I) if use_threaded_transform(tdst, transformer) - _add_general_kernel_threaded!(tdst, tsrc, p, transformer, α, β, backend...) + add_kernel_threaded!(style, tdst, tsrc, p, transformer, α, β, backend...) else - _add_general_kernel_nonthreaded!( - tdst, tsrc, p, transformer, α, β, backend... - ) + add_kernel_nonthreaded!(style, tdst, tsrc, p, transformer, α, β, backend...) end end end @@ -514,93 +505,131 @@ end # Trivial implementations # ----------------------- -function _add_trivial_kernel!(tdst, tsrc, p, transformer, α, β, backend...) +function add_trivial_kernel!(tdst, tsrc, p, transformer, α, β, backend...) TO.tensoradd!(tdst[], tsrc[], p, false, α, β, backend...) return nothing end -# Abelian implementations -# ----------------------- -function _add_abelian_kernel_nonthreaded!( - tdst, tsrc, p, transformer::AbelianTreeTransformer, α, β, backend... +# Non-threaded implementations +# ---------------------------- +function add_kernel_nonthreaded!( + ::UniqueFusion, tdst, tsrc, p, transformer, α, β, backend... + ) + for (f₁, f₂) in fusiontrees(tsrc) + _add_transform_single!(tdst, tsrc, p, (f₁, f₂), transformer, α, β, backend...) + end + return nothing +end +function add_kernel_nonthreaded!( + ::UniqueFusion, tdst, tsrc, p, transformer::AbelianTreeTransformer, α, β, backend... ) for subtransformer in transformer.data _add_transform_single!(tdst, tsrc, p, subtransformer, α, β, backend...) end return nothing end +function add_kernel_nonthreaded!(::FusionStyle, tdst, tsrc, p, transformer, α, β, backend...) + # preallocate buffers + buffers = allocate_buffers(tdst, tsrc, transformer) -function _add_abelian_kernel_threaded!( - tdst, tsrc, p, transformer::AbelianTreeTransformer, - α, β, backend...; + for src in fusionblocks(tsrc) + if length(src) == 1 + _add_transform_single!(tdst, tsrc, p, src, transformer, α, β, backend...) + else + _add_transform_multi!(tdst, tsrc, p, src, transformer, α, β, backend...) + end + end + return nothing +end +# specialization in the case of TensorMap +function add_kernel_nonthreaded!( + ::FusionStyle, tdst, tsrc, p, transformer::GenericTreeTransformer, α, β, backend... + ) + # preallocate buffers + buffers = allocate_buffers(tdst, tsrc, transformer) + + for subtransformer in transformer.data + # Special case without intermediate buffers whenever there is only a single block + if length(subtransformer[1]) == 1 + _add_transform_single!(tdst, tsrc, p, subtransformer, α, β, backend...) + else + _add_transform_multi!(tdst, tsrc, p, subtransformer, buffers, α, β, backend...) + end + end + return nothing +end + +# Threaded implementations +# ------------------------ +function add_kernel_threaded!( + ::UniqueFusion, tdst, tsrc, p, transformer, α, β, backend...; ntasks::Int = get_num_transformer_threads() ) - nblocks = length(transformer.data) + trees = fusiontrees(tsrc) + nblocks = length(trees) counter = Threads.Atomic{Int}(1) Threads.@sync for _ in 1:min(ntasks, nblocks) Threads.@spawn begin while true local_counter = Threads.atomic_add!(counter, 1) local_counter > nblocks && break - @inbounds subtransformer = transformer.data[local_counter] - _add_transform_single!(tdst, tsrc, p, subtransformer, α, β, backend...) + @inbounds (f₁, f₂) = trees[local_counter] + _add_transform_single!(tdst, tsrc, p, (f₁, f₂), transformer, α, β, backend...) end end end return nothing end - -function _add_transform_single!( - tdst, tsrc, p, (coeff, struct_dst, struct_src)::_AbelianTransformerData, - α, β, backend... +function add_kernel_threaded!( + ::UniqueFusion, tdst, tsrc, p, transformer::AbelianTreeTransformer, α, β, backend...; + ntasks::Int = get_num_transformer_threads() ) - subblock_dst = StridedView(tdst.data, struct_dst...) - subblock_src = StridedView(tsrc.data, struct_src...) - TO.tensoradd!(subblock_dst, subblock_src, p, false, α * coeff, β, backend...) - return nothing -end - -function _add_abelian_kernel_nonthreaded!(tdst, tsrc, p, transformer, α, β, backend...) - for (f₁, f₂) in fusiontrees(tsrc) - _add_abelian_block!(tdst, tsrc, p, transformer, f₁, f₂, α, β, backend...) - end - return nothing -end - -function _add_abelian_kernel_threaded!(tdst, tsrc, p, transformer, α, β, backend...) - Threads.@sync for (f₁, f₂) in fusiontrees(tsrc) - Threads.@spawn _add_abelian_block!(tdst, tsrc, p, transformer, f₁, f₂, α, β, backend...) + nblocks = length(transformer.data) + counter = Threads.Atomic{Int}(1) + Threads.@sync for _ in 1:min(ntasks, nblocks) + Threads.@spawn begin + while true + local_counter = Threads.atomic_add!(counter, 1) + local_counter > nblocks && break + @inbounds subtransformer = transformer.data[local_counter] + _add_transform_single!(tdst, tsrc, p, subtransformer, α, β, backend...) + end + end end return nothing end -function _add_abelian_block!(tdst, tsrc, p, transformer, f₁, f₂, α, β, backend...) - (f₁′, f₂′), coeff = first(transformer(f₁, f₂)) - @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, β, backend...) - return nothing -end - -# Non-abelian implementations -# --------------------------- -function _add_general_kernel_nonthreaded!( - tdst, tsrc, p, transformer::GenericTreeTransformer, α, β, backend... +function add_kernel_threaded!( + ::FusionStyle, tdst, tsrc, p, transformer, α, β, backend...; + ntasks::Int = get_num_transformer_threads() ) - # preallocate buffers - buffers = allocate_buffers(tdst, tsrc, transformer) + allblocks = fusionblocks(tsrc) + nblocks = length(allblocks) - for subtransformer in transformer.data - # Special case without intermediate buffers whenever there is only a single block - if length(subtransformer[1]) == 1 - _add_transform_single!(tdst, tsrc, p, subtransformer, α, β, backend...) - else - _add_transform_multi!(tdst, tsrc, p, subtransformer, buffers, α, β, backend...) + counter = Threads.Atomic{Int}(1) + Threads.@sync for _ in 1:min(ntasks, nblocks) + Threads.@spawn begin + # preallocate buffers for each task + buffers = allocate_buffers(tdst, tsrc, transformer) + + while true + local_counter = Threads.atomic_add!(counter, 1) + local_counter > nblocks && break + @inbounds src = allblocks[local_counter] + if length(src) == 1 + _add_transform_single!(tdst, tsrc, p, src, transformer, α, β, backend...) + else + _add_transform_multi!(tdst, tsrc, p, src, transformer, buffers, α, β, backend...) + end + end end end + return nothing end - -function _add_general_kernel_threaded!( - tdst, tsrc, p, transformer::GenericTreeTransformer, α, β, backend...; +# specialization in the case of TensorMap +function add_kernel_threaded!( + ::FusionStyle, tdst, tsrc, p, transformer::GenericTreeTransformer, α, β, backend...; ntasks::Int = get_num_transformer_threads() ) nblocks = length(transformer.data) @@ -627,20 +656,30 @@ function _add_general_kernel_threaded!( return nothing end -function _add_general_kernel_nonthreaded!(tdst, tsrc, p, transformer, α, β, backend...) - if iszero(β) - tdst = zerovector!(tdst) - elseif !isone(β) - tdst = scale!(tdst, β) - end - for (f₁, f₂) in fusiontrees(tsrc) - for ((f₁′, f₂′), coeff) in transformer(f₁, f₂) - @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, One(), backend...) - end - end +# Auxiliary methods +# ----------------- +function _add_transform_single!(tdst, tsrc, p, (f₁, f₂)::FusionTreePair, transformer, α, β, backend...) + (f₁′, f₂′), coeff = transformer((f₁, f₂)) + @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, β, backend...) + return nothing +end +function _add_transform_single!(tdst, tsrc, p, src::FusionTreeBlock, transformer, α, β, backend...) + dst, U = transformer(src) + f₁, f₂ = only(fusiontrees(src)) + f₁′, f₂′ = only(fusiontrees(dst)) + coeff = only(U) + @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, β, backend...) + return nothing +end +function _add_transform_single!( + tdst, tsrc, p, (coeff, struct_dst, struct_src)::_AbelianTransformerData, + α, β, backend... + ) + subblock_dst = StridedView(tdst.data, struct_dst...) + subblock_src = StridedView(tsrc.data, struct_src...) + TO.tensoradd!(subblock_dst, subblock_src, p, false, α * coeff, β, backend...) return nothing end - function _add_transform_single!( tdst, tsrc, p, (basistransform, structs_dst, structs_src)::_GenericTransformerData, α, β, backend... @@ -652,6 +691,36 @@ function _add_transform_single!( return nothing end +function _add_transform_multi!(tdst, tsrc, p, src::FusionTreeBlock, transformer, (buffer1, buffer2), α, β, backend...) + dst, U = transformer(src) + rows, cols = size(U) + sz_src = size(tsrc[first(fusiontrees(src))...]) + blocksize = prod(sz_src) + matsize = ( + prod(TupleTools.getindices(sz_src, codomainind(tsrc))), + prod(TupleTools.getindices(sz_src, domainind(tsrc))), + ) + + # Filling up a buffer with contiguous data + buffer_src = StridedView(buffers2, (blocksize, cols), (1, blocksize), 0) + for (i, (f₁, f₂)) in enumerate(fusiontrees(src)) + subblock_src = sreshape(tsrc[f₁, f₂], matsize) + _copyto!(buffer_src[:, i], subblock_src) + end + + # Resummation into a second buffer using BLAS + buffer_dst = StridedView(buffer1, (blocksize, rows), (1, blocksize), 0) + mul!(buffer_dst, buffer_src, U, α, Zero()) + + # Filling up the output + for (i, (f₃, f₄)) in enumerate(fusiontrees(dst)) + subblock_dst = tdst[f₃, f₄] + bufblock_dst = sreshape(buffer_dst[:, i], sz_src) + TO.tensoradd!(subblock_dst, bufblock_dst, p, false, One(), β, backend...) + end + + return nothing +end function _add_transform_multi!( tdst, tsrc, p, (basistransform, (sz_dst, structs_dst), (sz_src, structs_src)), (buffer1, buffer2), α, β, backend... @@ -683,25 +752,3 @@ function _add_transform_multi!( return nothing end - -function _add_general_kernel_threaded!(tdst, tsrc, p, transformer, α, β, backend...) - if iszero(β) - tdst = zerovector!(tdst) - elseif !isone(β) - tdst = scale!(tdst, β) - end - Threads.@sync for s₁ in sectors(codomain(tsrc)), s₂ in sectors(domain(tsrc)) - Threads.@spawn _add_nonabelian_sector!(tdst, tsrc, p, transformer, s₁, s₂, α, backend...) - end - return nothing -end - -function _add_nonabelian_sector!(tdst, tsrc, p, fusiontreetransform, s₁, s₂, α, backend...) - for (f₁, f₂) in fusiontrees(tsrc) - (f₁.uncoupled == s₁ && f₂.uncoupled == s₂) || continue - for ((f₁′, f₂′), coeff) in fusiontreetransform(f₁, f₂) - @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, One(), backend...) - end - end - return nothing -end diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index d09609edc..592a40a5b 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -73,8 +73,6 @@ end Construct the identity endomorphism on space `V`, i.e. return a `t::TensorMap` with `domain(t) == codomain(t) == V`, where either `scalartype(t) = T` if `T` is a `Number` type or `storagetype(t) = T` if `T` is a `DenseVector` type. - -See also [`one!`](@ref). """ id, id! id(V::TensorSpace) = id(Float64, V) diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl index 7bd0ab96b..091b59cd4 100644 --- a/src/tensors/tensoroperations.jl +++ b/src/tensors/tensoroperations.jl @@ -200,38 +200,63 @@ function trace_permute!( end I = sectortype(S) - # TODO: is it worth treating UniqueFusion separately? Is it worth to add multithreading support? if I === Trivial cod = codomain(tsrc) dom = domain(tsrc) n = length(cod) TO.tensortrace!(tdst[], tsrc[], (p₁, p₂), (q₁, q₂), false, α, β, backend) - # elseif FusionStyle(I) isa UniqueFusion - else - cod = codomain(tsrc) - dom = domain(tsrc) - n = length(cod) - scale!(tdst, β) - r₁ = (p₁..., q₁...) - r₂ = (p₂..., q₂...) + return tdst + end + + cod = codomain(tsrc) + dom = domain(tsrc) + n = length(cod) + scale!(tdst, β) + r₁ = (p₁..., q₁...) + r₂ = (p₂..., q₂...) + + # TODO: Is it worth to add multithreading support? + if FusionStyle(I) isa UniqueFusion for (f₁, f₂) in fusiontrees(tsrc) - for ((f₁′, f₂′), coeff) in permute(f₁, f₂, r₁, r₂) - f₁′′, g₁ = split(f₁′, N₁) - f₂′′, g₂ = split(f₂′, N₂) - g₁ == g₂ || continue - coeff *= dim(g₁.coupled) / dim(g₁.uncoupled[1]) - for i in 2:length(g₁.uncoupled) - if !(g₁.isdual[i]) - coeff *= twist(g₁.uncoupled[i]) + (f₁′, f₂′), coeff = permute((f₁, f₂), (r₁, r₂)) + f₁′′, g₁ = split(f₁′, N₁) + f₂′′, g₂ = split(f₂′, N₂) + g₁ == g₂ || continue + coeff *= dim(g₁.coupled) / dim(g₁.uncoupled[1]) + for i in 2:length(g₁.uncoupled) + if !(g₁.isdual[i]) + coeff *= twist(g₁.uncoupled[i]) + end + end + C = tdst[f₁′′, f₂′′] + A = tsrc[f₁, f₂] + α′ = α * coeff + TO.tensortrace!(C, A, (p₁, p₂), (q₁, q₂), false, α′, One(), backend) + end + else + for src in fusionblocks(tsrc) + dst, U = permute(src, (r₁, r₂)) + for (i, (f₁, f₂)) in enumerate(fusiontrees(src)) + for (j, (f₁′, f₂′)) in enumerate(fusiontrees(dst)) + coeff = U[j, i] + f₁′′, g₁ = split(f₁′, N₁) + f₂′′, g₂ = split(f₂′, N₂) + g₁ == g₂ || continue + coeff *= dim(g₁.coupled) / dim(g₁.uncoupled[1]) + for i in 2:length(g₁.uncoupled) + if !(g₁.isdual[i]) + coeff *= twist(g₁.uncoupled[i]) + end end + C = tdst[f₁′′, f₂′′] + A = tsrc[f₁, f₂] + α′ = α * coeff + TO.tensortrace!(C, A, (p₁, p₂), (q₁, q₂), false, α′, One(), backend) end - C = tdst[f₁′′, f₂′′] - A = tsrc[f₁, f₂] - α′ = α * coeff - TO.tensortrace!(C, A, (p₁, p₂), (q₁, q₂), false, α′, One(), backend) end end end + return tdst end diff --git a/src/tensors/treetransformers.jl b/src/tensors/treetransformers.jl index 36cd3926d..b7a636de0 100644 --- a/src/tensors/treetransformers.jl +++ b/src/tensors/treetransformers.jl @@ -26,7 +26,7 @@ function AbelianTreeTransformer(transform, p, Vdst, Vsrc) for i in 1:L f₁, f₂ = structure_src.fusiontreelist[i] - (f₃, f₄), coeff = only(transform(f₁, f₂)) + (f₃, f₄), coeff = transform((f₁, f₂)) j = structure_dst.fusiontreeindices[(f₃, f₄)] stridestructure_dst = structure_dst.fusiontreestructure[j] stridestructure_src = structure_src.fusiontreestructure[i] @@ -62,55 +62,72 @@ function GenericTreeTransformer(transform, p, Vdst, Vsrc) fusionstructure_dst = structure_dst.fusiontreestructure structure_src = fusionblockstructure(Vsrc) fusionstructure_src = structure_src.fusiontreestructure - I = sectortype(Vsrc) - - uncoupleds_src = map(structure_src.fusiontreelist) do (f₁, f₂) - return TupleTools.vcat(f₁.uncoupled, dual.(f₂.uncoupled)) - end - uncoupleds_src_unique = unique(uncoupleds_src) - - uncoupleds_dst = map(structure_dst.fusiontreelist) do (f₁, f₂) - return TupleTools.vcat(f₁.uncoupled, dual.(f₂.uncoupled)) - end + I = sectortype(Vsrc) T = sectorscalartype(I) N = numind(Vdst) - L = length(uncoupleds_src_unique) - data = Vector{_GenericTransformerData{T, N}}(undef, L) - - # TODO: this can be multithreaded - for (i, uncoupled) in enumerate(uncoupleds_src_unique) - inds_src = findall(==(uncoupled), uncoupleds_src) - fusiontrees_outer_src = structure_src.fusiontreelist[inds_src] - - uncoupled_dst = TupleTools.getindices(uncoupled, (p[1]..., p[2]...)) - inds_dst = findall(==(uncoupled_dst), uncoupleds_dst) - - fusiontrees_outer_dst = structure_dst.fusiontreelist[inds_dst] - - matrix = zeros(sectorscalartype(I), length(inds_dst), length(inds_src)) - for (row, (f₁, f₂)) in enumerate(fusiontrees_outer_src) - for ((f₃, f₄), coeff) in transform(f₁, f₂) - col = findfirst(==((f₃, f₄)), fusiontrees_outer_dst)::Int - matrix[row, col] = coeff + N₁ = numout(Vsrc) + N₂ = numin(Vsrc) + + fblocks = fusionblocks(Vsrc) + nblocks = length(fblocks) + data = Vector{_GenericTransformerData{T, N}}(undef, nblocks) + + nthreads = get_num_manipulation_threads() + if nthreads > 1 + counter = Threads.Atomic{Int}(1) + Threads.@sync for _ in 1:min(nthreads, nblocks) + Threads.@spawn begin + while true + local_counter = Threads.atomic_add!(counter, 1) + local_counter > nblocks && break + fs_src = fblocks[local_counter] + fs_dst, U = transform(fs_src) + matrix = copy(transpose(U)) # TODO: should we avoid this + + trees_src = fusiontrees(fs_src) + inds_src = map(Base.Fix1(getindex, structure_src.fusiontreeindices), trees_src) + trees_dst = fusiontrees(fs_dst) + inds_dst = map(Base.Fix1(getindex, structure_dst.fusiontreeindices), trees_dst) + + # size is shared between blocks, so repack: + # from [(sz, strides, offset), ...] to (sz, [(strides, offset), ...]) + sz_src, newstructs_src = repack_transformer_structure( + fusionstructure_src, inds_src + ) + sz_dst, newstructs_dst = repack_transformer_structure( + fusionstructure_dst, inds_dst + ) + + data[local_counter] = (matrix, (sz_dst, newstructs_dst), (sz_src, newstructs_src)) + end end end - - # size is shared between blocks, so repack: - # from [(sz, strides, offset), ...] to (sz, [(strides, offset), ...]) - sz_src, newstructs_src = repack_transformer_structure(fusionstructure_src, inds_src) - sz_dst, newstructs_dst = repack_transformer_structure(fusionstructure_dst, inds_dst) - - @debug( - "Created recoupling block for uncoupled: $uncoupled", - sz = size(matrix), sparsity = count(!iszero, matrix) / length(matrix) - ) - - data[i] = (matrix, (sz_dst, newstructs_dst), (sz_src, newstructs_src)) + transformer = GenericTreeTransformer{T, N}(data) + else + for (i, fs_src) in enumerate(fblocks) + fs_dst, U = transform(fs_src) + matrix = copy(transpose(U)) # TODO: should we avoid this + + trees_src = fusiontrees(fs_src) + inds_src = map(Base.Fix1(getindex, structure_src.fusiontreeindices), trees_src) + trees_dst = fusiontrees(fs_dst) + inds_dst = map(Base.Fix1(getindex, structure_dst.fusiontreeindices), trees_dst) + + # size is shared between blocks, so repack: + # from [(sz, strides, offset), ...] to (sz, [(strides, offset), ...]) + sz_src, newstructs_src = repack_transformer_structure( + fusionstructure_src, inds_src + ) + sz_dst, newstructs_dst = repack_transformer_structure( + fusionstructure_dst, inds_dst + ) + + data[i] = matrix, (sz_dst, newstructs_dst), (sz_src, newstructs_src) + end + transformer = GenericTreeTransformer{T, N}(data) end - transformer = GenericTreeTransformer{T, N}(data) - # sort by (approximate) weight to facilitate multi-threading strategies sort!(transformer) @@ -118,9 +135,9 @@ function GenericTreeTransformer(transform, p, Vdst, Vsrc) @debug( "TreeTransformer for $Vsrc to $Vdst via $p", - nblocks = length(data), - sz_median = size(data[cld(end, 2)][1], 1), - sz_max = size(data[1][1], 1), + nblocks = length(transformer.data), + sz_median = size(transformer.data[cld(end, 2)][1], 1), + sz_max = size(transformer.data[1][1], 1), Δt ) @@ -145,6 +162,13 @@ function allocate_buffers( sz = buffersize(transformer) return similar(tdst.data, sz), similar(tsrc.data, sz) end +function allocate_buffers( + tdst::AbstractTensorMap, tsrc::AbstractTensorMap, transformer + ) + # be pessimistic and assume the worst for now + sz = dim(space(tsrc)) + return similar(storagetype(tdst), sz), similar(storagetype(tsrc), sz) +end function treetransformertype(Vdst, Vsrc) I = sectortype(Vdst) @@ -171,7 +195,7 @@ end # braid is special because it has levels function treebraider(::AbstractTensorMap, ::AbstractTensorMap, p::Index2Tuple, levels) - return fusiontreetransform(f1, f2) = braid(f1, f2, levels..., p...) + return fusiontreetransform(f) = braid(f, p, levels) end function treebraider(tdst::TensorMap, tsrc::TensorMap, p::Index2Tuple, levels) return treebraider(space(tdst), space(tsrc), p, levels) @@ -179,7 +203,7 @@ end @cached function treebraider( Vdst::TensorMapSpace, Vsrc::TensorMapSpace, p::Index2Tuple, levels )::treetransformertype(Vdst, Vsrc) - fusiontreebraider(f1, f2) = braid(f1, f2, levels..., p...) + fusiontreebraider(f) = braid(f, p, levels) return TreeTransformer(fusiontreebraider, p, Vdst, Vsrc) end @@ -187,7 +211,7 @@ for (transform, treetransformer) in ((:permute, :treepermuter), (:transpose, :treetransposer)) @eval begin function $treetransformer(::AbstractTensorMap, ::AbstractTensorMap, p::Index2Tuple) - return fusiontreetransform(f1, f2) = $transform(f1, f2, p...) + return fusiontreetransform(f) = $transform(f, p) end function $treetransformer(tdst::TensorMap, tsrc::TensorMap, p::Index2Tuple) return $treetransformer(space(tdst), space(tsrc), p) @@ -195,7 +219,7 @@ for (transform, treetransformer) in @cached function $treetransformer( Vdst::TensorMapSpace, Vsrc::TensorMapSpace, p::Index2Tuple )::treetransformertype(Vdst, Vsrc) - fusiontreetransform(f1, f2) = $transform(f1, f2, p...) + fusiontreetransform(f) = $transform(f, p) return TreeTransformer(fusiontreetransform, p, Vdst, Vsrc) end end diff --git a/test/symmetries/fusiontrees.jl b/test/symmetries/fusiontrees.jl index c001f4e26..e577839f9 100644 --- a/test/symmetries/fusiontrees.jl +++ b/test/symmetries/fusiontrees.jl @@ -1,12 +1,18 @@ using Test, TestExtras using TensorKit +using TensorKit: FusionTreeBlock import TensorKit as TK using Random: randperm using TensorOperations +using MatrixAlgebraKit: isunitary # TODO: remove this once type_repr works for all included types using TensorKitSectors +_isunitary(x::Number; kwargs...) = isapprox(x * x', one(x); kwargs...) +_isunitary(x; kwargs...) = isunitary(x; kwargs...) +_isone(x; kwargs...) = isapprox(x, one(x); kwargs...) + @isdefined(TestSetup) || include("../setup.jl") using .TestSetup @@ -18,7 +24,7 @@ using .TestSetup in = rand(collect(⊗(out...))) numtrees = length(fusiontrees(out, in, isdual)) @test numtrees == count(n -> true, fusiontrees(out, in, isdual)) - while !(0 < numtrees < 30) + while !(0 < numtrees < 30) && !(one(in) in ⊗(out...)) out = ntuple(n -> randsector(I), N) in = rand(collect(⊗(out...))) numtrees = length(fusiontrees(out, in, isdual)) @@ -99,20 +105,21 @@ using .TestSetup @test c′ == one(c′) return t′ end - braid_i_to_1 = braid(f1, levels, (i, (1:(i - 1))..., ((i + 1):N)...)) - trees2 = Dict(_reinsert_partial_tree(t, f2) => c for (t, c) in braid_i_to_1) - trees3 = empty(trees2) - p = (((N + 1):(N + i - 1))..., (1:N)..., ((N + i):(2N - 1))...) - levels = ((i:(N + i - 1))..., (1:(i - 1))..., ((i + N):(2N - 1))...) - for (t, coeff) in trees2 - for (t′, coeff′) in braid(t, levels, p) - trees3[t′] = get(trees3, t′, zero(coeff′)) + coeff * coeff′ - end - end - for (t, coeff) in trees3 - coeff′ = get(trees, t, zero(coeff)) - @test isapprox(coeff′, coeff; atol = 1.0e-12, rtol = 1.0e-12) - end + # TODO: restore this? + # braid_i_to_1 = braid(f1, (i, (1:(i - 1))..., ((i + 1):N)...), levels) + # trees2 = Dict(_reinsert_partial_tree(t, f2) => c for (t, c) in braid_i_to_1) + # trees3 = empty(trees2) + # p = (((N + 1):(N + i - 1))..., (1:N)..., ((N + i):(2N - 1))...) + # levels = ((i:(N + i - 1))..., (1:(i - 1))..., ((i + N):(2N - 1))...) + # for (t, coeff) in trees2 + # for (t′, coeff′) in braid(t, p, levels) + # trees3[t′] = get(trees3, t′, zero(coeff′)) + coeff * coeff′ + # end + # end + # for (t, coeff) in trees3 + # coeff′ = get(trees, t, zero(coeff)) + # @test isapprox(coeff′, coeff; atol = 1.0e-12, rtol = 1.0e-12) + # end if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) Af1 = convert(Array, f1) @@ -154,7 +161,7 @@ using .TestSetup @test bf ≈ bf′ atol = 1.0e-12 end - d2 = @constinferred TK.planar_trace(f, (1, 3), (2, 4)) + d2 = @constinferred TK.planar_trace(f, ((1, 3), (2, 4))) oind2 = (5, 6, 7) bf2 = tensortrace(af, (:a, :a, :b, :b, :c, :d, :e)) bf2′ = zero(bf2) @@ -163,7 +170,7 @@ using .TestSetup end @test bf2 ≈ bf2′ atol = 1.0e-12 - d2 = @constinferred TK.planar_trace(f, (5, 6), (2, 1)) + d2 = @constinferred TK.planar_trace(f, ((5, 6), (2, 1))) oind2 = (3, 4, 7) bf2 = tensortrace(af, (:a, :b, :c, :d, :b, :a, :e)) bf2′ = zero(bf2) @@ -172,7 +179,7 @@ using .TestSetup end @test bf2 ≈ bf2′ atol = 1.0e-12 - d2 = @constinferred TK.planar_trace(f, (1, 4), (6, 3)) + d2 = @constinferred TK.planar_trace(f, ((1, 4), (6, 3))) bf2 = tensortrace(af, (:a, :b, :c, :c, :d, :a, :e)) bf2′ = zero(bf2) for (f2′, coeff) in d2 @@ -182,7 +189,7 @@ using .TestSetup q1 = (1, 3, 5) q2 = (2, 4, 6) - d3 = @constinferred TK.planar_trace(f, q1, q2) + d3 = @constinferred TK.planar_trace(f, (q1, q2)) bf3 = tensortrace(af, (:a, :a, :b, :b, :c, :c, :d)) bf3′ = zero(bf3) for (f3′, coeff) in d3 @@ -192,7 +199,7 @@ using .TestSetup q1 = (1, 3, 5) q2 = (6, 2, 4) - d3 = @constinferred TK.planar_trace(f, q1, q2) + d3 = @constinferred TK.planar_trace(f, (q1, q2)) bf3 = tensortrace(af, (:a, :b, :b, :c, :c, :a, :d)) bf3′ = zero(bf3) for (f3′, coeff) in d3 @@ -202,7 +209,7 @@ using .TestSetup q1 = (1, 2, 3) q2 = (6, 5, 4) - d3 = @constinferred TK.planar_trace(f, q1, q2) + d3 = @constinferred TK.planar_trace(f, (q1, q2)) bf3 = tensortrace(af, (:a, :b, :c, :c, :b, :a, :d)) bf3′ = zero(bf3) for (f3′, coeff) in d3 @@ -212,7 +219,7 @@ using .TestSetup q1 = (1, 2, 4) q2 = (6, 3, 5) - d3 = @constinferred TK.planar_trace(f, q1, q2) + d3 = @constinferred TK.planar_trace(f, (q1, q2)) bf3 = tensortrace(af, (:a, :b, :b, :c, :c, :a, :d)) bf3′ = zero(bf3) for (f3′, coeff) in d3 @@ -226,89 +233,86 @@ using .TestSetup @testset "Fusion tree $Istr: elementary artin braid" begin N = length(out) isdual = ntuple(n -> rand(Bool), N) - for in in ⊗(out...) - for i in 1:(N - 1) - for f in fusiontrees(out, in, isdual) - d1 = @constinferred TK.artin_braid(f, i) - @test norm(values(d1)) ≈ 1 - d2 = empty(d1) - for (f1, coeff1) in d1 - for (f2, coeff2) in TK.artin_braid(f1, i; inv = true) - d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 - end - end - for (f2, coeff2) in d2 - if f2 == f - @test coeff2 ≈ 1 - else - @test isapprox(coeff2, 0; atol = 1.0e-12, rtol = 1.0e-12) - end - end + if FusionStyle(I) isa UniqueFusion + for in in ⊗(out...) + src = only(fusiontrees(out, in, isdual)) + + for i in 1:(N - 1) + dst, U = @constinferred TK.artin_braid(src, i) + @test _isunitary(U) + dst′, U′ = @constinferred TK.artin_braid(dst, i; inv = true) + @test U' ≈ U′ end end - end - - f = rand(collect(it)) - d1 = TK.artin_braid(f, 2) - d2 = empty(d1) - for (f1, coeff1) in d1 - for (f2, coeff2) in TK.artin_braid(f1, 3) - d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 - end - end - d1 = d2 - d2 = empty(d1) - for (f1, coeff1) in d1 - for (f2, coeff2) in TK.artin_braid(f1, 3; inv = true) - d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 - end - end - d1 = d2 - d2 = empty(d1) - for (f1, coeff1) in d1 - for (f2, coeff2) in TK.artin_braid(f1, 2; inv = true) - d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 + else + src = FusionTreeBlock{I}((out, ()), (isdual, ())) + length(src) > 0 && for i in 1:(N - 1) + dst, U = @constinferred TK.artin_braid(src, i) + @test _isunitary(U) + dst′, U′ = @constinferred TK.artin_braid(dst, i; inv = true) + @test U' ≈ U′ end end - d1 = d2 - for (f1, coeff1) in d1 - if f1 == f - @test coeff1 ≈ 1 - else - @test isapprox(coeff1, 0; atol = 1.0e-12, rtol = 1.0e-12) + + # Test double braid-unbraid + if FusionStyle(I) isa UniqueFusion + in = rand(collect(⊗(out...))) + src = only(fusiontrees(out, in, isdual)) + dst, U = TK.artin_braid(src, 2) + dst′, U′ = TK.artin_braid(dst, 3) + dst″, U″ = TK.artin_braid(dst′, 3; inv = true) + dst‴, U‴ = TK.artin_braid(dst″, 2; inv = true) + @test U * U′ * U″ * U‴ ≈ 1 + else + src = FusionTreeBlock{I}((out, ()), (isdual, ())) + if length(src) > 0 + dst, U = TK.artin_braid(src, 2) + dst′, U′ = TK.artin_braid(dst, 3) + dst″, U″ = TK.artin_braid(dst′, 3; inv = true) + dst‴, U‴ = TK.artin_braid(dst″, 2; inv = true) + @test U * U′ * U″ * U‴ ≈ LinearAlgebra.I end end end + @testset "Fusion tree $Istr: braiding and permuting" begin - f = rand(collect(it)) p = tuple(randperm(N)...) ip = invperm(p) - levels = ntuple(identity, N) - d = @constinferred braid(f, levels, p) - d2 = Dict{typeof(f), valtype(d)}() levels2 = p - for (f2, coeff) in d - for (f1, coeff2) in braid(f2, levels2, ip) - d2[f1] = get(d2, f1, zero(coeff)) + coeff2 * coeff - end - end - for (f1, coeff2) in d2 - if f1 == f - @test coeff2 ≈ 1 - else - @test isapprox(coeff2, 0; atol = 1.0e-12, rtol = 1.0e-12) - end + + if FusionStyle(I) isa UniqueFusion + f = only(fusiontrees(out, in, isdual)) + else + f = FusionTreeBlock{I}((out, ()), (isdual, ())) end + dst, U = @constinferred braid(f, (p, ()), (levels, ())) + @test _isunitary(U) + dst′, U′ = braid(dst, (ip, ()), (levels2, ())) + @test U' ≈ U′ + + if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) - Af = convert(Array, f) - Afp = permutedims(Af, (p..., N + 1)) - Afp2 = zero(Afp) - for (f1, coeff) in d - Afp2 .+= coeff .* convert(Array, f1) + if FusionStyle(I) isa UniqueFusion + A = permutedims(fusiontensor(f), (p..., N + 1)) + A′ = U * fusiontensor(dst) + @test A ≈ A′ + else + A = map(x -> permutedims(fusiontensor(x[1]), (p..., N + 1)), fusiontrees(f)) + A′ = map(fusiontensor ∘ first, fusiontrees(dst)) + for (i, Ai) in enumerate(A) + Aj = sum(U[i, :] .* A′) + @test Ai ≈ Aj + end end - @test Afp ≈ Afp2 + # Af = convert(Array, f) + # Afp = permutedims(Af, (p..., N + 1)) + # Afp2 = zero(Afp) + # for (f1, coeff) in d + # Afp2 .+= coeff .* convert(Array, f1) + # end + # @test Afp ≈ Afp2 end end @@ -347,11 +351,21 @@ using .TestSetup end perm = ((N .+ (1:N))..., (1:N)...) levels = ntuple(identity, 2 * N) + for (t, coeff) in trees1 - for (t′, coeff′) in braid(t, levels, perm) + if FusionStyle(I) isa UniqueFusion + t′, coeff′ = braid(t, perm, levels) trees3[t′] = get(trees3, t′, zero(valtype(trees3))) + coeff * coeff′ + else + src = FusionTreeBlock([(t, FusionTree{I}((), t.coupled, (), ()))]) + dst, U = braid(src, (perm, ()), (levels, ())) + for ((t′, _), coeff′) in zip(fusiontrees(dst), U) + trees3[t′] = get(trees3, t′, zero(valtype(trees3))) + coeff * coeff′ + + end end end + for (t, coeff) in trees3 coeff′ = get(trees2, t, zero(coeff)) @test isapprox(coeff, coeff′; atol = 1.0e-12, rtol = 1.0e-12) @@ -393,58 +407,47 @@ using .TestSetup numtrees = count(n -> true, fusiontrees((out..., map(dual, out)...))) end incoming = rand(collect(⊗(out...))) - f1 = rand(collect(fusiontrees(out, incoming, ntuple(n -> rand(Bool), N)))) - f2 = rand(collect(fusiontrees(out[randperm(N)], incoming, ntuple(n -> rand(Bool), N)))) + + if FusionStyle(I) isa UniqueFusion + f1 = rand(collect(fusiontrees(out, incoming, ntuple(n -> rand(Bool), N)))) + f2 = rand(collect(fusiontrees(out[randperm(N)], incoming, ntuple(n -> rand(Bool), N)))) + src = (f1, f2) + if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) + A = fusiontensor(src) + end + else + src = FusionTreeBlock{I}((out, out), (ntuple(n -> rand(Bool), N), ntuple(n -> rand(Bool), N))) + if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) + A = map(fusiontensor, fusiontrees(src)) + end + end @testset "Double fusion tree $Istr: repartitioning" begin for n in 0:(2 * N) - d = @constinferred TK.repartition(f1, f2, $n) - @test dim(incoming) ≈ - sum(abs2(coef) * dim(f1.coupled) for ((f1, f2), coef) in d) - d2 = Dict{typeof((f1, f2)), valtype(d)}() - for ((f1′, f2′), coeff) in d - for ((f1′′, f2′′), coeff2) in TK.repartition(f1′, f2′, N) - d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff2 * coeff - end - end - for ((f1′, f2′), coeff2) in d2 - if f1 == f1′ && f2 == f2′ - @test coeff2 ≈ 1 - else - @test isapprox(coeff2, 0; atol = 1.0e-12, rtol = 1.0e-12) - end - end + dst, U = @constinferred TK.repartition(src, $n) + # @test _isunitary(U) + + dst′, U′ = repartition(dst, N) + @test _isone(U * U′) + if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) - Af1 = convert(Array, f1) - Af2 = permutedims(convert(Array, f2), [N:-1:1; N + 1]) - sz1 = size(Af1) - sz2 = size(Af2) - d1 = prod(sz1[1:(end - 1)]) - d2 = prod(sz2[1:(end - 1)]) - dc = sz1[end] - A = reshape( - reshape(Af1, (d1, dc)) * reshape(Af2, (d2, dc))', - (sz1[1:(end - 1)]..., sz2[1:(end - 1)]...) - ) - A2 = zero(A) - for ((f1′, f2′), coeff) in d - Af1′ = convert(Array, f1′) - Af2′ = permutedims(convert(Array, f2′), [(2N - n):-1:1; 2N - n + 1]) - sz1′ = size(Af1′) - sz2′ = size(Af2′) - d1′ = prod(sz1′[1:(end - 1)]) - d2′ = prod(sz2′[1:(end - 1)]) - dc′ = sz1′[end] - A2 += coeff * - reshape( - reshape(Af1′, (d1′, dc′)) * reshape(Af2′, (d2′, dc′))', - (sz1′[1:(end - 1)]..., sz2′[1:(end - 1)]...) - ) + all_inds = (ntuple(identity, numout(src))..., reverse(ntuple(i -> i + numout(src), numin(src)))...) + p₁ = ntuple(i -> all_inds[i], numout(dst)) + p₂ = reverse(ntuple(i -> all_inds[i + numout(dst)], numin(dst))) + if FusionStyle(I) isa UniqueFusion + @test permutedims(A, (p₁..., p₂...)) ≈ U * fusiontensor(dst) + else + A′ = map(Base.Fix2(permutedims, (p₁..., p₂...)), A) + A″ = map(fusiontensor, fusiontrees(dst)) + for (i, Ai) in enumerate(A′) + Aj = sum(U[:, i] .* A″) + @test Ai ≈ Aj + end end - @test A ≈ A2 end end end + @testset "Double fusion tree $Istr: permutation" begin if BraidingStyle(I) isa SymmetricBraiding for n in 0:(2N) @@ -453,54 +456,52 @@ using .TestSetup ip = invperm(p) ip1, ip2 = ip[1:N], ip[(N + 1):(2N)] - d = @constinferred TK.permute(f1, f2, p1, p2) - @test dim(incoming) ≈ - sum(abs2(coef) * dim(f1.coupled) for ((f1, f2), coef) in d) - d2 = Dict{typeof((f1, f2)), valtype(d)}() - for ((f1′, f2′), coeff) in d - d′ = TK.permute(f1′, f2′, ip1, ip2) - for ((f1′′, f2′′), coeff2) in d′ - d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + - coeff2 * coeff - end - end - for ((f1′, f2′), coeff2) in d2 - if f1 == f1′ && f2 == f2′ - @test coeff2 ≈ 1 - else - @test abs(coeff2) < 1.0e-12 - end - end + dst, U = @constinferred TensorKit.permute(src, (p1, p2)) + # @test _isunitary(U) + dst′, U′ = @constinferred TensorKit.permute(dst, (ip1, ip2)) + # @test U' ≈ U′ + @test _isone(U * U′) if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) - Af1 = convert(Array, f1) - Af2 = convert(Array, f2) - sz1 = size(Af1) - sz2 = size(Af2) - d1 = prod(sz1[1:(end - 1)]) - d2 = prod(sz2[1:(end - 1)]) - dc = sz1[end] - A = reshape( - reshape(Af1, (d1, dc)) * reshape(Af2, (d2, dc))', - (sz1[1:(end - 1)]..., sz2[1:(end - 1)]...) - ) - Ap = permutedims(A, (p1..., p2...)) - A2 = zero(Ap) - for ((f1′, f2′), coeff) in d - Af1′ = convert(Array, f1′) - Af2′ = convert(Array, f2′) - sz1′ = size(Af1′) - sz2′ = size(Af2′) - d1′ = prod(sz1′[1:(end - 1)]) - d2′ = prod(sz2′[1:(end - 1)]) - dc′ = sz1′[end] - A2 += coeff * reshape( - reshape(Af1′, (d1′, dc′)) * - reshape(Af2′, (d2′, dc′))', - (sz1′[1:(end - 1)]..., sz2′[1:(end - 1)]...) - ) + if FusionStyle(I) isa UniqueFusion + @test permutedims(A, (p1..., p2...)) ≈ U * fusiontensor(dst) + else + A′ = map(Base.Fix2(permutedims, (p1..., p2...)), A) + A″ = map(fusiontensor, fusiontrees(dst)) + for (i, Ai) in enumerate(A′) + Aj = sum(U[:, i] .* A″) + @test Ai ≈ Aj + end end - @test Ap ≈ A2 + # + # Af1 = convert(Array, f1) + # Af2 = convert(Array, f2) + # sz1 = size(Af1) + # sz2 = size(Af2) + # d1 = prod(sz1[1:(end - 1)]) + # d2 = prod(sz2[1:(end - 1)]) + # dc = sz1[end] + # A = reshape( + # reshape(Af1, (d1, dc)) * reshape(Af2, (d2, dc))', + # (sz1[1:(end - 1)]..., sz2[1:(end - 1)]...) + # ) + # Ap = permutedims(A, (p1..., p2...)) + # A2 = zero(Ap) + # for ((f1′, f2′), coeff) in d + # Af1′ = convert(Array, f1′) + # Af2′ = convert(Array, f2′) + # sz1′ = size(Af1′) + # sz2′ = size(Af2′) + # d1′ = prod(sz1′[1:(end - 1)]) + # d2′ = prod(sz2′[1:(end - 1)]) + # dc′ = sz1′[end] + # A2 += coeff * reshape( + # reshape(Af1′, (d1′, dc′)) * + # reshape(Af2′, (d2′, dc′))', + # (sz1′[1:(end - 1)]..., sz2′[1:(end - 1)]...) + # ) + # end + # @test Ap ≈ A2 end end end @@ -515,75 +516,81 @@ using .TestSetup ip′ = tuple(getindex.(Ref(vcat(1:n, (2N):-1:(n + 1))), ip)...) ip1, ip2 = ip′[1:N], ip′[(2N):-1:(N + 1)] - d = @constinferred transpose(f1, f2, p1, p2) - @test dim(incoming) ≈ - sum(abs2(coef) * dim(f1.coupled) for ((f1, f2), coef) in d) - d2 = Dict{typeof((f1, f2)), valtype(d)}() - for ((f1′, f2′), coeff) in d - d′ = transpose(f1′, f2′, ip1, ip2) - for ((f1′′, f2′′), coeff2) in d′ - d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff2 * coeff - end - end - for ((f1′, f2′), coeff2) in d2 - if f1 == f1′ && f2 == f2′ - @test coeff2 ≈ 1 - else - @test abs(coeff2) < 1.0e-12 - end - end + dst, U = @constinferred transpose(src, (p1, p2)) + # @test _isunitary(U) + dst′, U′ = @constinferred transpose(dst, (ip1, ip2)) + # @test U' ≈ U′ + @test _isone(U * U′) if BraidingStyle(I) isa Bosonic - d3 = permute(f1, f2, p1, p2) - for (f1′, f2′) in union(keys(d), keys(d3)) - coeff1 = get(d, (f1′, f2′), zero(valtype(d))) - coeff3 = get(d3, (f1′, f2′), zero(valtype(d3))) - @test isapprox(coeff1, coeff3; atol = 1.0e-12) - end + dst″, U″ = permute(src, (p1, p2)) + @test U″ ≈ U end if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) - Af1 = convert(Array, f1) - Af2 = convert(Array, f2) - sz1 = size(Af1) - sz2 = size(Af2) - d1 = prod(sz1[1:(end - 1)]) - d2 = prod(sz2[1:(end - 1)]) - dc = sz1[end] - A = reshape( - reshape(Af1, (d1, dc)) * reshape(Af2, (d2, dc))', - (sz1[1:(end - 1)]..., sz2[1:(end - 1)]...) - ) - Ap = permutedims(A, (p1..., p2...)) - A2 = zero(Ap) - for ((f1′, f2′), coeff) in d - Af1′ = convert(Array, f1′) - Af2′ = convert(Array, f2′) - sz1′ = size(Af1′) - sz2′ = size(Af2′) - d1′ = prod(sz1′[1:(end - 1)]) - d2′ = prod(sz2′[1:(end - 1)]) - dc′ = sz1′[end] - A2 += coeff * reshape( - reshape(Af1′, (d1′, dc′)) * - reshape(Af2′, (d2′, dc′))', - (sz1′[1:(end - 1)]..., sz2′[1:(end - 1)]...) - ) + if FusionStyle(I) isa UniqueFusion + @test permutedims(A, (p1..., p2...)) ≈ U * fusiontensor(dst) + else + A′ = map(Base.Fix2(permutedims, (p1..., p2...)), A) + A″ = map(fusiontensor, fusiontrees(dst)) + for (i, Ai) in enumerate(A′) + Aj = sum(U[:, i] .* A″) + @test Ai ≈ Aj + end end - @test Ap ≈ A2 + + # Af1 = convert(Array, f1) + # Af2 = convert(Array, f2) + # sz1 = size(Af1) + # sz2 = size(Af2) + # d1 = prod(sz1[1:(end - 1)]) + # d2 = prod(sz2[1:(end - 1)]) + # dc = sz1[end] + # A = reshape( + # reshape(Af1, (d1, dc)) * reshape(Af2, (d2, dc))', + # (sz1[1:(end - 1)]..., sz2[1:(end - 1)]...) + # ) + # Ap = permutedims(A, (p1..., p2...)) + # A2 = zero(Ap) + # for ((f1′, f2′), coeff) in d + # Af1′ = convert(Array, f1′) + # Af2′ = convert(Array, f2′) + # sz1′ = size(Af1′) + # sz2′ = size(Af2′) + # d1′ = prod(sz1′[1:(end - 1)]) + # d2′ = prod(sz2′[1:(end - 1)]) + # dc′ = sz1′[end] + # A2 += coeff * reshape( + # reshape(Af1′, (d1′, dc′)) * + # reshape(Af2′, (d2′, dc′))', + # (sz1′[1:(end - 1)]..., sz2′[1:(end - 1)]...) + # ) + # end + # @test Ap ≈ A2 end end end - @testset "Double fusion tree $Istr: planar trace" begin - d1 = transpose(f1, f1, (N + 1, 1:N..., ((2N):-1:(N + 3))...), (N + 2,)) + false && @testset "Double fusion tree $Istr: planar trace" begin + + + if FusionStyle(I) isa UniqueFusion + f1, f1 = src + dst, U = transpose((f1, f1), ((N + 1, 1:N..., ((2N):-1:(N + 3))...), (N + 2,))) + d1 = zip((dst,), (U,)) + else + f1, f1 = first(fusiontrees(src)) + dst, U = transpose(src, ((N + 1, 1:N..., ((2N):-1:(N + 3))...), (N + 2,))) + @show size(U) + d1 = zip(fusiontrees(dst), U[:, 1]) + end + f1front, = TK.split(f1, N - 1) T = TensorKitSectors._Fscalartype(I) d2 = Dict{typeof((f1front, f1front)), T}() for ((f1′, f2′), coeff′) in d1 for ((f1′′, f2′′), coeff′′) in TK.planar_trace( - f1′, f2′, (2:N...,), (1, ((2N):-1:(N + 3))...), (N + 1,), - (N + 2,) + (f1′, f2′), ((2:N...,), (1, ((2N):-1:(N + 3))...)), ((N + 1,), (N + 2,)) ) coeff = coeff′ * coeff′′ d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff