diff --git a/docs/Project.toml b/docs/Project.toml index 695415bc4..bb471e442 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,7 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f" [compat] diff --git a/docs/make.jl b/docs/make.jl index 0c2e0a35a..398fde0f7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,7 @@ using Documenter using Random using TensorKit, TensorKitSectors +using TensorKit: FusionTreePair, Index2Tuple pages = ["Home" => "index.md", "Manual" => ["man/intro.md", "man/tutorial.md", "man/categories.md", diff --git a/docs/src/lib/sectors.md b/docs/src/lib/sectors.md index 4029b68cf..09e725b0e 100644 --- a/docs/src/lib/sectors.md +++ b/docs/src/lib/sectors.md @@ -90,7 +90,7 @@ insertat split merge elementary_trace -planar_trace(f::FusionTree{I,N}, q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {I<:Sector,N,N₃} +planar_trace(f::FusionTree{I,N}, q::Index2Tuple{N₃,N₃}) where {I<:Sector,N,N₃} artin_braid braid(f::FusionTree{I,N}, levels::NTuple{N,Int}, p::NTuple{N,Int}) where {I<:Sector,N} permute(f::FusionTree{I,N}, p::NTuple{N,Int}) where {I<:Sector,N} @@ -113,8 +113,8 @@ Finally, these are used to define large manipulations of fusion-splitting tree p are then used in the index manipulation of `AbstractTensorMap` objects. The following methods defined on fusion splitting tree pairs have an associated definition for tensors. ```@docs -repartition(::FusionTree{I,N₁}, ::FusionTree{I,N₂}, ::Int) where {I<:Sector,N₁,N₂} -transpose(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂} -braid(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple, ::IndexTuple, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂} -permute(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂} +repartition(::FusionTreePair, ::Int) +transpose(::FusionTreePair{I}, ::Index2Tuple{N₁,N₂}) where {I,N₁,N₂} +braid(::FusionTreePair, ::Index2Tuple, ::Index2Tuple) +permute(::FusionTreePair, ::Index2Tuple) ``` diff --git a/docs/src/man/sectors.md b/docs/src/man/sectors.md index 30590ef00..768e3f608 100644 --- a/docs/src/man/sectors.md +++ b/docs/src/man/sectors.md @@ -1155,7 +1155,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 @@ -1180,7 +1180,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 @@ -1206,7 +1206,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 6f7467e37..1bcd39135 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -198,6 +198,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/fusiontrees/fusiontreeblocks.jl b/src/fusiontrees/fusiontreeblocks.jl new file mode 100644 index 000000000..ecd06739b --- /dev/null +++ b/src/fusiontrees/fusiontreeblocks.jl @@ -0,0 +1,677 @@ +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 + +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 + +Base.propertynames(::FusionTreeBlock, private::Bool=false) = (:trees, :uncoupled, :isdual) + +sectortype(::Type{<:FusionTreeBlock{I}}) where {I} = I +numout(fs::FusionTreeBlock) = numout(typeof(fs)) +numout(::Type{<:FusionTreeBlock{I,N₁}}) where {I,N₁} = N₁ +numin(fs::FusionTreeBlock) = numin(typeof(fs)) +numin(::Type{<:FusionTreeBlock{I,N₁,N₂}}) where {I,N₁,N₂} = N₂ +numind(fs::FusionTreeBlock) = numind(typeof(fs)) +numind(::Type{T}) where {T<:FusionTreeBlock} = numin(T) + numout(T) + +fusiontrees(block::FusionTreeBlock) = block.trees +Base.length(block::FusionTreeBlock) = length(fusiontrees(block)) + +# 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 + +# Manipulations +# ------------- +function transformation_matrix(transform, dst::FusionTreeBlock{I}, + src::FusionTreeBlock{I}) where {I} + U = zeros(sectorscalartype(I), length(dst), length(src)) + indexmap = treeindex_map(dst) + for (col, f) in enumerate(fusiontrees(src)) + for (f′, c) in transform(f) + row = indexmap[f′] + U[row, col] = c + end + end + return U +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 ? leftone(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(frobeniusschur(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[(f₁′, f₂′)] + @inbounds U[row, col] = coeff + end + end + end + + return dst, U +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 ? leftone(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(frobeniusschur(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 + +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) + 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)) + # 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(frobeniusschur(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 = (leftone(c1),) # or rightone(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 + +# !! 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(frobeniusschur(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 = (leftone(c1),) # or rightone(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 + +function cycleclockwise(src::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::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 + +@inline function repartition(src::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 +``` +=# +@generated function repartition(src::FusionTreeBlock, ::Val{N}) where {N} + return _repartition_body(numout(src) - N) +end +function _repartition_body(N) + if N == 0 + ex = quote + T = sectorscalartype(sectortype(src)) + U = copyto!(zeros(T, length(src), length(src)), LinearAlgebra.I) + return src, U + 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 + +function Base.transpose(src::FusionTreeBlock, p::Index2Tuple{N₁,N₂}) where {N₁,N₂} + N = N₁ + N₂ + @assert numind(src) == N + p′ = linearizepermutation(p..., numout(src), numin(src)) + @assert iscyclicpermutation(p′) + return _fstranspose((src, p)) +end + +const _FSTransposeKey{I,N₁,N₂} = Tuple{<:FusionTreeBlock{I},Index2Tuple{N₁,N₂}} + +@cached function _fstranspose(key::_FSTransposeKey{I,N₁,N₂})::Tuple{FusionTreeBlock{I,N₁, + N₂}, + Matrix{sectorscalartype(I)}} where {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 + +function CacheStyle(::typeof(_fstranspose), k::_FSTransposeKey{I}) where {I} + if FusionStyle(I) == UniqueFusion() + return NoCache() + else + return GlobalLRUCache() + end +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′ = rightone(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 + +function braid(src::FusionTreeBlock{I,N,0}, p::NTuple{N,Int}, + levels::NTuple{N,Int}) where {I,N} + TupleTools.isperm(p) || throw(ArgumentError("not a valid permutation: $p")) + + if FusionStyle(I) isa UniqueFusion && BraidingStyle(I) isa SymmetricBraiding + uncoupled′ = TupleTools._permute(src.uncoupled[1], p) + isdual′ = TupleTools._permute(src.isdual[1], p) + dst = FusionTreeBlock{I}(uncoupled′, isdual′; sizehint=length(src)) + U = transformation_matrix(dst, src) do (f₁, f₂) + return ((f₁′, f₂) => c for (f₁′, c) in braid(f₁, p, levels)) + end + else + dst, U = repartition(src, N) # TODO: can we avoid this? + for s in permutation2swaps(p) + inv = levels[s] > levels[s + 1] + dst, U_tmp = artin_braid(dst, s; inv) + U = U_tmp * U + end + end + return dst, U +end + +function braid(src::FusionTreeBlock{I}, p::Index2Tuple{N₁,N₂}, + levels::Index2Tuple) where {I,N₁,N₂} + @assert numind(src) == N₁ + N₂ + @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 _FSBraidKey{I,N₁,N₂} = Tuple{<:FusionTreeBlock{I},Index2Tuple{N₁,N₂},Index2Tuple} + +@cached function _fsbraid(key::_FSBraidKey{I,N₁,N₂})::Tuple{FusionTreeBlock{I,N₁,N₂, + fusiontreetype(I, + N₁, + N₂)}, + Matrix{sectorscalartype(I)}} where {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)) + + if FusionStyle(I) isa UniqueFusion && BraidingStyle(I) isa SymmetricBraiding + uncoupled′ = TupleTools._permute(dst.uncoupled[1], p) + isdual′ = TupleTools._permute(dst.isdual[1], p) + + dst′ = FusionTreeBlock{I}(uncoupled′, isdual′) + U_tmp = transformation_matrix(dst′, dst) do (f₁, f₂) + return ((f₁′, f₂) => c for (f₁, c) in braid(f₁, p, levels)) + end + dst = dst′ + U = U_tmp * U + else + for s in permutation2swaps(p) + inv = levels[s] > levels[s + 1] + dst, U_tmp = artin_braid(dst, s; inv) + U = U_tmp * U + end + end + + if N₂ == 0 + return dst, U + else + dst, U_tmp = repartition(dst, N₁) + U = U_tmp * U + return dst, U + end +end + +function CacheStyle(::typeof(_fsbraid), k::_FSBraidKey{I}) where {I} + if FusionStyle(I) isa UniqueFusion + return NoCache() + else + return GlobalLRUCache() + end +end + +function permute(src::FusionTreeBlock{I}, p::Index2Tuple) where {I} + @assert BraidingStyle(I) 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/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index 3e694e523..64daff4f7 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -92,6 +92,8 @@ function FusionTree(uncoupled::NTuple{N,I}, coupled::I, end FusionTree(uncoupled::Tuple{I,Vararg{I}}) where {I<:Sector} = FusionTree(uncoupled, one(I)) +const FusionTreePair{I,N₁,N₂} = Tuple{FusionTree{I,N₁},FusionTree{I,N₂}} + # Properties sectortype(::Type{<:FusionTree{I}}) where {I<:Sector} = I FusionStyle(::Type{<:FusionTree{I}}) where {I<:Sector} = FusionStyle(I) @@ -134,7 +136,7 @@ end Base.:(==)(f₁::FusionTree, f₂::FusionTree) = false # 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} if N === 0 FusionTree{I,0,0,0} elseif N === 1 @@ -143,6 +145,10 @@ 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 # converting to actual array function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I,0}) where {I} @@ -199,8 +205,7 @@ 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} +function Base.convert(A::Type{<:AbstractArray}, (f₁, f₂)::FusionTreePair{I}) where {I} F₁ = convert(A, f₁) F₂ = convert(A, f₂) sz1 = size(F₁) @@ -224,11 +229,12 @@ 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") +include("fusiontreeblocks.jl") + +# Manipulate fusion trees +include("manipulations.jl") # auxiliary routines # _abelianinner: generate the inner indices for given outer indices in the abelian case diff --git a/src/fusiontrees/manipulations.jl b/src/fusiontrees/manipulations.jl index f19bb06ed..ed9869760 100644 --- a/src/fusiontrees/manipulations.jl +++ b/src/fusiontrees/manipulations.jl @@ -244,8 +244,7 @@ end # -> 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₂} +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] @@ -274,19 +273,18 @@ function flip(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}, i::Int; 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₂} +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)) + (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₂} +function bendright((f₁, f₂)::FusionTreePair{I,N₁,N₂}) where {I,N₁,N₂} # map final splitting vertex (a, b)<-c to fusion vertex a<-(c, dual(b)) @assert N₁ > 0 c = f₁.coupled @@ -332,15 +330,14 @@ function bendright(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {I< end end # change to N₁ + 1, N₂ - 1 -function bendleft(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I} +function bendleft((f₁, f₂)::FusionTreePair{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₁)) + 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₂} +function foldright((f₁, f₂)::FusionTreePair{I,N₁,N₂}) where {I,N₁,N₂} # map first splitting vertex (a, b)<-c to fusion vertex b<-(dual(a), c) @assert N₁ > 0 a = f₁.uncoupled[1] @@ -359,7 +356,6 @@ function foldright(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {I< 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 = (leftone(c1),) # or rightone(a) @@ -372,6 +368,7 @@ function foldright(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {I< 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 @@ -380,7 +377,7 @@ function foldright(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {I< 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₂) + for (fr, coeff2) in fr_coeffs coeff = factor * coeff1 * conj(coeff2) if (@isdefined newtrees) newtrees[(fl, fr)] = get(newtrees, (fl, fr), zero(coeff)) + @@ -396,11 +393,10 @@ function foldright(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {I< end end # change to N₁ + 1, N₂ - 1 -function foldleft(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I} +function foldleft((f₁, f₂)::FusionTreePair{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₁)) + for ((f₂′, f₁′), coeff) in foldright((f₂, f₁))) end # COMPOSITE DUALITY MANIPULATIONS PART 1: Repartition and transpose @@ -423,11 +419,11 @@ function iscyclicpermutation(v1, v2) end # clockwise cyclic permutation while preserving (N₁, N₂): foldright & bendleft -function cycleclockwise(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I<:Sector} +function cycleclockwise((f₁, f₂)::FusionTreePair{I}) where {I} local newtrees if length(f₁) > 0 - for ((f1a, f2a), coeffa) in foldright(f₁, f₂) - for ((f1b, f2b), coeffb) in bendleft(f1a, f2a) + 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 @@ -437,8 +433,8 @@ function cycleclockwise(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I<:Sect end end else - for ((f1a, f2a), coeffa) in bendleft(f₁, f₂) - for ((f1b, f2b), coeffb) in foldright(f1a, f2a) + 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 @@ -452,11 +448,11 @@ function cycleclockwise(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I<:Sect end # anticlockwise cyclic permutation while preserving (N₁, N₂): foldleft & bendright -function cycleanticlockwise(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I<:Sector} +function cycleanticlockwise((f₁, f₂)::FusionTreePair{I}) where {I} local newtrees if length(f₂) > 0 - for ((f1a, f2a), coeffa) in foldleft(f₁, f₂) - for ((f1b, f2b), coeffb) in bendright(f1a, f2a) + 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 @@ -466,8 +462,8 @@ function cycleanticlockwise(f₁::FusionTree{I}, f₂::FusionTree{I}) where {I<: end end else - for ((f1a, f2a), coeffa) in bendright(f₁, f₂) - for ((f1b, f2b), coeffb) in foldleft(f1a, f2a) + 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 @@ -482,8 +478,8 @@ 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} + 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 @@ -492,35 +488,32 @@ outgoing (`f₁`) and incoming sectors (`f₂`) respectively (with identical cou 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₂} +@inline function repartition((f₁, f₂)::FusionTreePair, N::Int) f₁.coupled == f₂.coupled || throw(SectorMismatch()) - @assert 0 <= N <= N₁ + N₂ - return _recursive_repartition(f₁, f₂, Val(N)) + @assert 0 <= N <= length(f₁) + length(f₂) + 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} +function _recursive_repartition((f₁, f₂)::FusionTreePair{I,N₁,N₂}, + ::Val{N}) where {I,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) + FF = Tuple{F₁,F₂} 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)) + 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) + newtrees = fusiontreedict(I){FF,T}((f₁′′, f₂′′) => coeff1 * coeff2) end end end @@ -529,9 +522,8 @@ function _recursive_repartition(f₁::FusionTree{I,N₁}, 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} + 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 @@ -540,17 +532,15 @@ outgoing (`t1`) and incoming sectors (`t2`) respectively (with identical coupled 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₂} +function Base.transpose((f₁, f₂)::FusionTreePair{I}, p::Index2Tuple{N₁,N₂}) where {I,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)) + p′ = linearizepermutation(p..., length(f₁), length(f₂)) + @assert iscyclicpermutation(p′) + return fstranspose(((f₁, f₂), p)) end -const FSTransposeKey{I<:Sector,N₁,N₂} = Tuple{<:FusionTree{I},<:FusionTree{I}, - IndexTuple{N₁},IndexTuple{N₂}} +const FSTransposeKey{I,N₁,N₂} = Tuple{<:FusionTreePair{I},Index2Tuple{N₁,N₂}} function _fsdicttype(I, N₁, N₂) F₁ = fusiontreetype(I, N₁) @@ -560,13 +550,11 @@ function _fsdicttype(I, N₁, N₂) end @cached function fstranspose(key::FSTransposeKey{I,N₁,N₂})::_fsdicttype(I, N₁, - N₂) where {I<:Sector, - N₁, - N₂} - f₁, f₂, p1, p2 = key + N₂) where {I,N₁,N₂} + (f₁, f₂), (p1, p2) = key N = N₁ + N₂ p = linearizepermutation(p1, p2, length(f₁), length(f₂)) - newtrees = repartition(f₁, f₂, N₁) + newtrees = repartition((f₁, f₂), N₁) length(p) == 0 && return newtrees i1 = findfirst(==(1), p) @assert i1 !== nothing @@ -575,7 +563,7 @@ end while 1 < i1 <= Nhalf local newtrees′ for ((f1a, f2a), coeffa) in newtrees - for ((f1b, f2b), coeffb) in cycleanticlockwise(f1a, f2a) + for ((f1b, f2b), coeffb) in cycleanticlockwise((f1a, f2a)) coeff = coeffa * coeffb if (@isdefined newtrees′) newtrees′[(f1b, f2b)] = get(newtrees′, (f1b, f2b), zero(coeff)) + coeff @@ -590,7 +578,7 @@ end while Nhalf < i1 local newtrees′ for ((f1a, f2a), coeffa) in newtrees - for ((f1b, f2b), coeffb) in cycleclockwise(f1a, f2a) + for ((f1b, f2b), coeffb) in cycleclockwise((f1a, f2a)) coeff = coeffa * coeffb if (@isdefined newtrees′) newtrees′[(f1b, f2b)] = get(newtrees′, (f1b, f2b), zero(coeff)) + coeff @@ -605,7 +593,7 @@ end return newtrees end -function CacheStyle(::typeof(fstranspose), k::FSTransposeKey{I}) where {I<:Sector} +function CacheStyle(::typeof(fstranspose), k::FSTransposeKey{I}) where {I} if FusionStyle(I) isa UniqueFusion return NoCache() else @@ -618,13 +606,12 @@ end # -> 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₃} +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) + return transpose((f₁, f₂), (p1, p2)) end linearindex = (ntuple(identity, Val(length(f₁)))..., @@ -641,9 +628,9 @@ function planar_trace(f₁::FusionTree{I}, f₂::FusionTree{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′) + 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 @@ -655,15 +642,14 @@ function planar_trace(f₁::FusionTree{I}, f₂::FusionTree{I}, end """ - planar_trace(f::FusionTree{I,N}, q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {I<:Sector,N,N₃} + 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::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {I<:Sector,N,N₃} +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}() @@ -697,7 +683,7 @@ function planar_trace(f::FusionTree{I,N}, 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′) + for (f′′, coeff′′) in planar_trace(f′, (q1′, q2′)) coeff = coeff′ * coeff′′ if !iszero(coeff) newtrees[f′′] = get(newtrees, f′′, zero(coeff)) + coeff @@ -709,13 +695,13 @@ 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} + 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<:Sector,N} +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 || isone(f.coupled) || @@ -820,7 +806,7 @@ applying `artin_braid(f′, i; inv = true)` to all the outputs `f′` of 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} +function artin_braid(f::FusionTree{I,N}, 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 = f.uncoupled @@ -898,7 +884,8 @@ function artin_braid(f::FusionTree{I,N}, i; inv::Bool=false) where {I<:Sector,N} return fusiontreedict(I)(f′ => coeff) elseif FusionStyle(I) isa SimpleFusion local newtrees - for c′ in intersect(a ⊗ d, e ⊗ conj(b)) + 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)) * @@ -919,7 +906,8 @@ function artin_braid(f::FusionTree{I,N}, i; inv::Bool=false) where {I<:Sector,N} return newtrees else # GenericFusion local newtrees - for c′ in intersect(a ⊗ d, e ⊗ conj(b)) + 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) @@ -950,7 +938,7 @@ end # braid fusion tree """ - braid(f::FusionTree{<:Sector, N}, levels::NTuple{N, Int}, p::NTuple{N, Int}) + 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 @@ -963,9 +951,7 @@ that if `i` and `j` cross, ``τ_{i,j}`` is applied if `levels[i] < levels[j]` an ``τ_{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} +function braid(f::FusionTree{I,N}, p::NTuple{N,Int}, levels::NTuple{N,Int}) where {I,N} TupleTools.isperm(p) || throw(ArgumentError("not a valid permutation: $p")) if FusionStyle(I) isa UniqueFusion && BraidingStyle(I) isa SymmetricBraiding coeff = one(sectorscalartype(I)) @@ -1012,17 +998,15 @@ Perform a permutation of the uncoupled indices of the fusion tree `f` and return 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} +function permute(f::FusionTree{I,N}, p::NTuple{N,Int}) where {I,N} @assert BraidingStyle(I) isa SymmetricBraiding - return braid(f, ntuple(identity, Val(N)), p) + return braid(f, p, ntuple(identity, Val(N))) 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} + 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 @@ -1036,27 +1020,23 @@ respectively, which determines how indices braid. In particular, if `i` and `j` 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₂ +function braid((f₁, f₂)::FusionTreePair, (p1, p2)::Index2Tuple, + (levels1, levels2)::Index2Tuple) + @assert length(f₁) + length(f₂) == length(p1) + length(p2) @assert length(f₁) == length(levels1) && length(f₂) == length(levels2) @assert TupleTools.isperm((p1..., p2...)) - return fsbraid((f₁, f₂, levels1, levels2, p1, p2)) + return fsbraid(((f₁, f₂), (p1, p2), (levels1, levels2))) end -const FSBraidKey{I<:Sector,N₁,N₂} = Tuple{<:FusionTree{I},<:FusionTree{I}, - IndexTuple,IndexTuple, - IndexTuple{N₁},IndexTuple{N₂}} +const FSBraidKey{I,N₁,N₂} = Tuple{<:FusionTreePair{I},Index2Tuple{N₁,N₂},Index2Tuple} -@cached function fsbraid(key::FSBraidKey{I,N₁,N₂})::_fsdicttype(I, N₁, - N₂) where {I<:Sector,N₁,N₂} - (f₁, f₂, l1, l2, p1, p2) = key +@cached function fsbraid(key::FSBraidKey{I,N₁,N₂})::_fsdicttype(I, N₁, N₂) where {I,N₁,N₂} + ((f₁, f₂), (p1, p2), (l1, l2)) = 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₁) + for ((f, f0), coeff1) in repartition((f₁, f₂), N₁ + N₂) + for (f′, coeff2) in braid(f, p, levels) + for ((f₁′, f₂′), coeff3) in repartition((f′, f0), N₁) if @isdefined newtrees newtrees[(f₁′, f₂′)] = get(newtrees, (f₁′, f₂′), zero(coeff3)) + coeff1 * coeff2 * coeff3 @@ -1078,9 +1058,8 @@ function CacheStyle(::typeof(fsbraid), k::FSBraidKey{I}) where {I<:Sector} 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} + 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 @@ -1089,10 +1068,9 @@ outgoing (`t1`) and incoming sectors (`t2`) respectively (with identical coupled 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 +function permute((f₁, f₂)::FusionTreePair, (p1, p2)::Index2Tuple) + @assert BraidingStyle(sectortype(f₁)) isa SymmetricBraiding levels1 = ntuple(identity, length(f₁)) levels2 = length(f₁) .+ ntuple(identity, length(f₂)) - return braid(f₁, f₂, levels1, levels2, p1, p2) + return braid((f₁, f₂), (p1, p2), (levels1, levels2)) end diff --git a/src/planar/planaroperations.jl b/src/planar/planaroperations.jl index 15fce65d0..8d689834f 100644 --- a/src/planar/planaroperations.jl +++ b/src/planar/planaroperations.jl @@ -95,7 +95,7 @@ function planartrace!(C::AbstractTensorMap, 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, β′, backend, allocator) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index a4a184e5f..3889f7551 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -230,7 +230,7 @@ function planarcontract!(C::AbstractTensorMap, 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′′ @@ -281,7 +281,7 @@ function planarcontract!(C::AbstractTensorMap, 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 88d0c3b25..7b0ea9cf2 100644 --- a/src/tensors/diagonal.jl +++ b/src/tensors/diagonal.jl @@ -203,7 +203,7 @@ function permute(d::DiagonalTensorMap, (p₁, p₂)::Index2Tuple{1,1}; d′ = typeof(d)(undef, dual(d.domain)) for (c, b) in blocks(d) f = only(fusiontrees(codomain(d), c)) - ((f′, _), coeff) = only(permute(f, f, p₁, p₂)) + ((f′, _), coeff) = only(permute((f, f), (p₁, p₂))) c′ = f′.coupled scale!(block(d′, c′), b, coeff) end diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl index c4e2b9228..6cda0c584 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′ @@ -558,7 +558,7 @@ function _add_abelian_kernel_threaded!(tdst, tsrc, p, transformer, α, β, backe end function _add_abelian_block!(tdst, tsrc, p, transformer, f₁, f₂, α, β, backend...) - (f₁′, f₂′), coeff = first(transformer(f₁, f₂)) + (f₁′, f₂′), coeff = first(transformer((f₁, f₂))) @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, β, backend...) return nothing @@ -618,7 +618,7 @@ function _add_general_kernel_nonthreaded!(tdst, tsrc, p, transformer, α, β, ba tdst = scale!(tdst, β) end for (f₁, f₂) in fusiontrees(tsrc) - for ((f₁′, f₂′), coeff) in transformer(f₁, f₂) + for ((f₁′, f₂′), coeff) in transformer((f₁, f₂)) @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, One(), backend...) end @@ -683,7 +683,7 @@ 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₂) + for ((f₁′, f₂′), coeff) in fusiontreetransform((f₁, f₂)) @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, One(), backend...) end diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index f29bdf809..8a3c3d1e9 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -71,8 +71,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 6110093b6..5910c3322 100644 --- a/src/tensors/tensoroperations.jl +++ b/src/tensors/tensoroperations.jl @@ -199,7 +199,7 @@ function trace_permute!(tdst::AbstractTensorMap, r₁ = (p₁..., q₁...) r₂ = (p₂..., q₂...) for (f₁, f₂) in fusiontrees(tsrc) - for ((f₁′, f₂′), coeff) in permute(f₁, f₂, r₁, r₂) + for ((f₁′, f₂′), coeff) in permute((f₁, f₂), (r₁, r₂)) f₁′′, g₁ = split(f₁′, N₁) f₂′′, g₂ = split(f₂′, N₂) g₁ == g₂ || continue diff --git a/src/tensors/treetransformers.jl b/src/tensors/treetransformers.jl index 046af5edc..8c9601ad3 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 = only(transform((f₁, f₂))) j = structure_dst.fusiontreeindices[(f₃, f₄)] stridestructure_dst = structure_dst.fusiontreestructure[j] stridestructure_src = structure_src.fusiontreestructure[i] @@ -62,62 +62,99 @@ 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) + N₁ = numout(Vsrc) + N₂ = numin(Vsrc) - # 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] + isdual_src = (map(isdual, codomain(Vsrc).spaces), map(isdual, domain(Vsrc).spaces)) - uncoupled_dst = TupleTools.getindices(uncoupled, (p[1]..., p[2]...)) - inds_dst = findall(==(uncoupled_dst), uncoupleds_dst) + data = Vector{_GenericTransformerData{T,N}}() - fusiontrees_outer_dst = structure_dst.fusiontreelist[inds_dst] + nthreads = get_num_manipulation_threads() + if nthreads > 1 + fusiontreeblocks = Vector{FusionTreeBlock{I,N₁,N₂,fusiontreetype(I, N₁, N₂)}}() + for cod_uncoupled_src in sectors(codomain(Vsrc)), + dom_uncoupled_src in sectors(domain(Vsrc)) - 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 + fs_src = FusionTreeBlock{I}((cod_uncoupled_src, dom_uncoupled_src), isdual_src) + trees_src = fusiontrees(fs_src) + if !isempty(trees_src) + push!(fusiontreeblocks, fs_src) 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)) + nblocks = length(fusiontreeblocks) + + resize!(data, nblocks) + 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 = fusiontreeblocks[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 + transformer = GenericTreeTransformer{T,N}(data) + else + isdual_src = (map(isdual, codomain(Vsrc).spaces), map(isdual, domain(Vsrc).spaces)) + for cod_uncoupled_src in sectors(codomain(Vsrc)), + dom_uncoupled_src in sectors(domain(Vsrc)) + + fs_src = FusionTreeBlock{I}((cod_uncoupled_src, dom_uncoupled_src), isdual_src) + trees_src = fusiontrees(fs_src) + isempty(trees_src) && continue + + fs_dst, U = transform(fs_src) + matrix = copy(transpose(U)) # TODO: should we avoid this + + 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) + + push!(data, (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) Δt = Base.time() - t₀ @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) return transformer @@ -166,14 +203,14 @@ 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) 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 @@ -181,14 +218,14 @@ 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) end @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/fusiontrees.jl b/test/fusiontrees.jl index d758f9a83..f94e52c6e 100644 --- a/test/fusiontrees.jl +++ b/test/fusiontrees.jl @@ -89,13 +89,13 @@ ti = time() @test c′ == one(c′) return t′ end - braid_i_to_1 = braid(f1, levels, (i, (1:(i - 1))..., ((i + 1):N)...)) + 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, levels, p) + for (t′, coeff′) in braid(t, p, levels) trees3[t′] = get(trees3, t′, zero(coeff′)) + coeff * coeff′ end end @@ -142,7 +142,7 @@ ti = time() @test bf ≈ bf′ atol = 1e-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) @@ -151,7 +151,7 @@ ti = time() end @test bf2 ≈ bf2′ atol = 1e-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) @@ -160,7 +160,7 @@ ti = time() end @test bf2 ≈ bf2′ atol = 1e-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 @@ -170,7 +170,7 @@ ti = time() 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 @@ -180,7 +180,7 @@ ti = time() 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 @@ -190,7 +190,7 @@ ti = time() 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 @@ -200,7 +200,7 @@ ti = time() 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 @@ -273,11 +273,11 @@ ti = time() ip = invperm(p) levels = ntuple(identity, N) - d = @constinferred braid(f, levels, p) + d = @constinferred braid(f, p, levels) d2 = Dict{typeof(f),valtype(d)}() levels2 = p for (f2, coeff) in d - for (f1, coeff2) in braid(f2, levels2, ip) + for (f1, coeff2) in braid(f2, ip, levels2) d2[f1] = get(d2, f1, zero(coeff)) + coeff2 * coeff end end @@ -334,7 +334,7 @@ ti = time() perm = ((N .+ (1:N))..., (1:N)...) levels = ntuple(identity, 2 * N) for (t, coeff) in trees1 - for (t′, coeff′) in braid(t, levels, perm) + for (t′, coeff′) in braid(t, perm, levels) trees3[t′] = get(trees3, t′, zero(valtype(trees3))) + coeff * coeff′ end end @@ -381,12 +381,12 @@ ti = time() @testset "Double fusion tree $Istr: repartioning" begin for n in 0:(2 * N) - d = @constinferred TK.repartition(f1, f2, $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) + for ((f1′′, f2′′), coeff2) in TK.repartition((f1′, f2′), N) d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff2 * coeff end end @@ -432,12 +432,12 @@ ti = time() ip = invperm(p) ip1, ip2 = ip[1:N], ip[(N + 1):(2N)] - d = @constinferred TensorKit.permute(f1, f2, p1, p2) + d = @constinferred TensorKit.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′ = TensorKit.permute(f1′, f2′, ip1, ip2) + d′ = TensorKit.permute((f1′, f2′), (ip1, ip2)) for ((f1′′, f2′′), coeff2) in d′ d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff2 * coeff @@ -490,12 +490,12 @@ ti = time() 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) + 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) + d′ = transpose((f1′, f2′), (ip1, ip2)) for ((f1′′, f2′′), coeff2) in d′ d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff2 * coeff end @@ -509,7 +509,7 @@ ti = time() end if BraidingStyle(I) isa Bosonic - d3 = permute(f1, f2, p1, p2) + 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))) @@ -546,14 +546,15 @@ ti = time() end end @testset "Double fusion tree $Istr: planar trace" begin - d1 = transpose(f1, f1, (N + 1, 1:N..., ((2N):-1:(N + 3))...), (N + 2,)) + d1 = transpose((f1, f1), ((N + 1, 1:N..., ((2N):-1:(N + 3))...), (N + 2,))) f1front, = TK.split(f1, N - 1) T = typeof(Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1, 1, 1, 1]) 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,)) + TK.planar_trace((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 end