diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md index aa3d1b297..8f88c8ad6 100644 --- a/docs/src/Changelog.md +++ b/docs/src/Changelog.md @@ -86,6 +86,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - Fusion tree block manipulations of (named) product sectors are computed as Kronecker products of the manipulations of their parts, splitting off the `UniqueFusion` sectors, such that they share the caches of the individual sectors ([#592](https://github.com/QuantumKitHub/TensorKit.jl/pull/592)) - The TensorOperations ChainRules rules no longer copy or retain `C` when `β = Zero()`, as is the case for every tensor allocated by `@tensor` ([#555](https://github.com/QuantumKitHub/TensorKit.jl/pull/555)) - The Enzyme `norm` rule reuses the cached primal and the `tr` rule avoids an unnecessary copy ([#554](https://github.com/QuantumKitHub/TensorKit.jl/pull/554)) +- Fusion tree block manipulations no longer construct and index the canonical destination block, but keep the trees of the result in the order in which they are produced; the buffers of the intermediate blocks of a braid are recycled, and blocks with a single tree are braided as their tree pair (internal) ## [0.17.2](https://github.com/QuantumKitHub/TensorKit.jl/compare/v0.17.1...v0.17.2) - 2026-09-20 diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl index b737f354f..b7fe74ce6 100644 --- a/src/fusiontrees/braiding_manipulations.jl +++ b/src/fusiontrees/braiding_manipulations.jl @@ -14,10 +14,12 @@ 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`. + +For sectors without `UniqueFusion`, the braid of `f` has to be a single fusion tree, which is the +case if `f` is the only fusion tree with its uncoupled and coupled sectors; this throws otherwise. """ 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] @@ -53,7 +55,7 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N} if i == 1 c = N > 2 ? inner[1] : coupled′ - R = oftype(oneT, (inv ? conj(Rsymbol(b, a, c)) : Rsymbol(a, b, c))) + R = oftype(oneT, (inv ? conj(only(Rsymbol(b, a, c))) : only(Rsymbol(a, b, c)))) f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices) return f′ => R end @@ -64,13 +66,13 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N} a = inner_extended[i - 1] c = inner_extended[i] e = inner_extended[i + 1] - c′ = only(a ⊗ d) + c′ = only(_artin_channels(a, b, d, e)) coeff = oftype( oneT, if inv - conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′) + conj(only(Rsymbol(d, c, e)) * only(Fsymbol(d, a, b, e, c′, c))) * only(Rsymbol(d, a, c′)) else - Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′)) + only(Rsymbol(c, d, e)) * conj(only(Fsymbol(d, a, b, e, c′, c)) * only(Rsymbol(a, d, c′))) end ) inner′ = TupleTools.setindex(inner, c′, i - 1) @@ -78,9 +80,19 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N} return f′ => coeff end +# the inner lines `c′` of `a ⊗ d → c′ → c′ ⊗ b → e`, i.e. `intersect(a ⊗ d, e ⊗ conj(b))`, +# without constructing the second fusion product +_artin_channels(a, b, d, e) = Iterators.filter(c′ -> !iszero(Nsymbol(c′, b, e)), a ⊗ d) + function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where {I, N} + n = length(src) # braiding a complete block keeps its size + return artin_braid!(Matrix{sectorscalartype(I)}(undef, n, n), src, i; inv) # zeroed by `artin_braid!` +end +# `U` is overwritten with the braid matrix of `src`, a complete block +function artin_braid!(U::AbstractMatrix, 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")) - isempty(src) && return src => zeros(sectorscalartype(I), 0, 0) + size(U) == (length(src), length(src)) || throw(DimensionMismatch(lazy"braid matrix of size $(size(U)) for a block of $(length(src)) trees")) + isempty(src) && return src => U uncoupled = src.uncoupled[1] a, b = uncoupled[i], uncoupled[i + 1] @@ -91,12 +103,10 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where 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)) + rows = TreeRows(fusiontreetype(I, N, 0), length(src)) oneT = one(sectorscalartype(I)) - - indexmap = treeindex_map(dst) - U = zeros(sectorscalartype(I), length(dst), length(src)) + fill!(U, zero(oneT)) if isunit(a) || isunit(b) # braiding with trivial sector: simple and always possible for (col, (f, f₂)) in enumerate(fusiontrees(src)) @@ -113,10 +123,10 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where vertices′ = TupleTools.setindex(vertices′, vertices[i - 1], i) end f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′) - row = indexmap[treeindex_data((f′, f₂))] + row = row!(rows, (f′, f₂)) @inbounds U[row, col] = oneT end - return dst => U + return FusionTreeBlock(rows.trees) => U end BraidingStyle(I) isa NoBraiding && @@ -132,7 +142,7 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where 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₂))] + row = row!(rows, (f′, f₂)) @inbounds U[row, col] = R else # GenericFusion μ = vertices[1] @@ -142,7 +152,7 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where iszero(R) && continue vertices′ = TupleTools.setindex(vertices, ν, 1) f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices′) - row = indexmap[treeindex_data((f′, f₂))] + row = row!(rows, (f′, f₂)) @inbounds U[row, col] = R end end @@ -155,38 +165,39 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where c = inner_extended[i] e = inner_extended[i + 1] if FusionStyle(I) isa MultiplicityFreeFusion - for c′ in intersect(a ⊗ d, e ⊗ conj(b)) + R = inv ? Rsymbol(d, c, e) : Rsymbol(c, d, e) + for c′ in _artin_channels(a, b, d, e) coeff = if inv - conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′) + conj(R * 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′)) + R * 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₂))] + row = row!(rows, (f′, f₂)) @inbounds U[row, col] = coeff end else # GenericFusion - for c′ in intersect(a ⊗ d, e ⊗ conj(b)) - Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e) + Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e) + μ = vertices[i - 1] + ν = vertices[i] + for c′ in _artin_channels(a, b, 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) + inner′ = TupleTools.setindex(inner, c′, i - 1) + for σ in axes(Rmat2, 1) + vertices_σ = TupleTools.setindex(vertices, σ, i - 1) + for λ in axes(Fmat, 2) coeff = zero(oneT) - for ρ in 1:Nsymbol(d, c, e), κ in 1:Nsymbol(d, a, c′) + for ρ in axes(Rmat1, 2), κ in axes(Rmat2, 2) 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) + vertices′ = TupleTools.setindex(vertices_σ, λ, i) f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′) - row = indexmap[treeindex_data((f′, f₂))] + row = row!(rows, (f′, f₂)) @inbounds U[row, col] = coeff end end @@ -194,7 +205,7 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where end end - return dst => U + return FusionTreeBlock(rows.trees) => U end # braid fusion tree @@ -216,9 +227,8 @@ 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(lazy"not a valid permutation: $p")) - @assert FusionStyle(I) isa UniqueFusion braid_style = BraidingStyle(I) - if braid_style isa SymmetricBraiding # this assumes Fsymbols are 1! + if FusionStyle(I) isa UniqueFusion && braid_style isa SymmetricBraiding # Fsymbols are 1 coeff = one(sectorscalartype(I)) for i in 1:N for j in 1:(i - 1) @@ -310,42 +320,52 @@ end @cached function fsbraid(key::K)::_fsdicttype(K) where {I, N₁, N₂, K <: Union{FSPBraidKey{I, N₁, N₂}, FSBBraidKey{I, N₁, N₂}}} if K <: FSPBraidKey - ((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 + f, p, levels = key + return _braid_single(f::FusionTreePair{I}, p, levels) else src, (p1, p2), (l1, l2) = key + src = src::FusionTreeBlock{I} + + if length(src) == 1 # then so is every intermediate block: the braid is a product of scalars + f′, coeff = _braid_single(only(fusiontrees(src)), (p1, p2), (l1, l2)) + return FusionTreeBlock([f′]) => fill(coeff, 1, 1) + end p = linearizepermutation(p1, p2, numout(src), numin(src)) levels = (l1..., reverse(l2)...) - dst, U = repartition(src, numind(src)) + dst, U = repartition(src, Val(numind(src))) + # the swaps keep the size of the block, so that the matrices of the previous step are recycled + U_tmp = Matrix{sectorscalartype(I)}(undef, length(dst), length(dst)) + U_spare = nothing braid_style = BraidingStyle(I) for s in permutation2swaps(p) _check_levels(braid_style, levels, s) inv = levels[s] > levels[s + 1] - dst, U_tmp = artin_braid(dst, s; inv) - U = U_tmp * U + dst, _ = artin_braid!(U_tmp, dst, s; inv) + U_spare, U = U, U_spare isa typeof(U_tmp) ? mul!(U_spare, U_tmp, 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::FusionTreeBlock{I}, N₁) - U = U_tmp * U - return dst => U - end + N₂ == 0 && return dst => U + dst, U_tmp = repartition(dst::FusionTreeBlock{I}, Val(N₁)) + return dst => U_tmp * U end end +# also for a tree pair that is the only one with its uncoupled sectors (`only` throws otherwise) +function _braid_single((f₁, f₂)::FusionTreePair, (p1, p2)::Index2Tuple{N₁, N₂}, (l1, l2)::Index2Tuple) where {N₁, N₂} + 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 + CacheStyle(::typeof(fsbraid), k::Union{FSPBraidKey{I}, FSBBraidKey{I}}) where {I} = FusionStyle(I) isa UniqueFusion ? NoCache() : GlobalLRUCache() diff --git a/src/fusiontrees/duality_manipulations.jl b/src/fusiontrees/duality_manipulations.jl index 9e613d50e..171ef2ba7 100644 --- a/src/fusiontrees/duality_manipulations.jl +++ b/src/fusiontrees/duality_manipulations.jl @@ -53,16 +53,15 @@ function _bendright_treepair((f₁, f₂)::FusionTreePair) return (a, b, c), (f₁′, f₂′) end +# also for a final splitting vertex without multiplicity (`only` throws otherwise) function bendright((f₁, f₂)::FusionTreePair) - I = sectortype((f₁, f₂)) N₁ = numout((f₁, f₂)) - @assert FusionStyle(I) === UniqueFusion() (a, b, c), (f₁′, f₂′) = _bendright_treepair((f₁, f₂)) # compute the coefficient coeff₀ = sqrtdim(c) * invsqrtdim(a) f₁.isdual[N₁] && (coeff₀ *= conj(frobenius_schur_phase(dual(b)))) - coeff = coeff₀ * Bsymbol(a, b, c) + coeff = coeff₀ * only(Bsymbol(a, b, c)) return (f₁′, f₂′) => coeff end @@ -71,18 +70,9 @@ function bendright(src::FusionTreeBlock) N₁ = numout(src) N₂ = numin(src) @assert N₁ > 0 - uncoupled_dst = ( - TupleTools.front(src.uncoupled[1]), - (src.uncoupled[2]..., dual(src.uncoupled[1][N₁])), - ) - isdual_dst = ( - TupleTools.front(src.isdual[1]), - (src.isdual[2]..., !(src.isdual[1][N₁])), - ) - - dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src)) - indexmap = treeindex_map(dst) - U = zeros(fusionscalartype(I), length(dst), length(src)) + + rows = TreeRows(fusiontreetype(I, N₁ - 1, N₂ + 1), length(src)) + U = zeros(fusionscalartype(I), length(src), length(src)) for (col, (f₁, f₂)) in enumerate(fusiontrees(src)) (a, b, c), (f₁′, f₂′) = _bendright_treepair((f₁, f₂)) @@ -92,7 +82,7 @@ function bendright(src::FusionTreeBlock) end if FusionStyle(I) isa MultiplicityFreeFusion coeff = coeff₀ * Bsymbol(a, b, c) - row = indexmap[treeindex_data((f₁′, f₂′))] + row = row!(rows, (f₁′, f₂′)) @inbounds U[row, col] = coeff else Bmat = Bsymbol(a, b, c) @@ -106,12 +96,12 @@ function bendright(src::FusionTreeBlock) iszero(coeff) && continue vertices₂ = N₂ > 0 ? (f₂.vertices..., ν) : () f₂′ = FusionTree(uncoupled₂, coupled₂, isdual₂, inner₂, vertices₂) - row = indexmap[treeindex_data((f₁′, f₂′))] + row = row!(rows, (f₁′, f₂′)) @inbounds U[row, col] = coeff end end end - return dst => U + return FusionTreeBlock(rows.trees) => U end @doc """ @@ -138,7 +128,6 @@ 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 @@ -150,18 +139,9 @@ function bendleft(src::FusionTreeBlock) N₁ = numout(src) N₂ = numin(src) @assert N₂ > 0 - uncoupled_dst = ( - (src.uncoupled[1]..., dual(src.uncoupled[2][N₂])), - TupleTools.front(src.uncoupled[2]), - ) - isdual_dst = ( - (src.isdual[1]..., !(src.isdual[2][N₂])), - TupleTools.front(src.isdual[2]), - ) - - dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src)) - indexmap = treeindex_map(dst) - U = zeros(fusionscalartype(I), length(dst), length(src)) + + rows = TreeRows(fusiontreetype(I, N₁ + 1, N₂ - 1), length(src)) + U = zeros(fusionscalartype(I), length(src), length(src)) for (col, (f₁, f₂)) in enumerate(fusiontrees(src)) (a, b, c), (f₂′, f₁′) = _bendright_treepair((f₂, f₁)) @@ -171,7 +151,7 @@ function bendleft(src::FusionTreeBlock) end if FusionStyle(I) isa MultiplicityFreeFusion coeff = coeff₀ * Bsymbol(a, b, c) - row = indexmap[treeindex_data((f₁′, f₂′))] + row = row!(rows, (f₁′, f₂′)) @inbounds U[row, col] = conj(coeff) else Bmat = Bsymbol(a, b, c) @@ -185,12 +165,12 @@ function bendleft(src::FusionTreeBlock) iszero(coeff) && continue vertices₁ = N₁ > 0 ? (f₁.vertices..., ν) : () f₁′ = FusionTree(uncoupled₁, coupled₁, isdual₁, inner₁, vertices₁) - row = indexmap[treeindex_data((f₁′, f₂′))] + row = row!(rows, (f₁′, f₂′)) @inbounds U[row, col] = conj(coeff) end end end - return dst => U + return FusionTreeBlock(rows.trees) => U end @doc """ @@ -236,17 +216,11 @@ function foldright((f₁, f₂)::FusionTreePair) 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) + rows = TreeRows(fusiontreetype(I, N₁ - 1, numin(src) + 1), length(src)) f₁, f₂ = first(fusiontrees(src)) a::I = f₁.uncoupled[1] @@ -260,7 +234,7 @@ function foldright(src::FusionTreeBlock) c::I = f₁.coupled cache₃ = Dict((b, c) => Asymbol(a, b, c)) - U = zeros(eltype(coeff₁), length(dst), length(src)) + U = zeros(eltype(coeff₁), length(src), length(src)) for (col, (f₁, f₂)) in enumerate(fusiontrees(src)) f₁′s, coeffs₁ = get!(cache₁, f₁) do multi_Fmove(f₁) @@ -280,12 +254,12 @@ function foldright(src::FusionTreeBlock) if isduala coeff *= κₐ end - row = indexmap[treeindex_data((f₁′, f₂′))] + row = row!(rows, (f₁′, f₂′)) @inbounds U[row, col] += coeff end end end - return dst => U + return FusionTreeBlock(rows.trees) => U end @@ -319,21 +293,12 @@ function foldleft((f₁, f₂)::FusionTreePair) end 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) + rows = TreeRows(fusiontreetype(I, N₂ + 1, N₁ - 1), length(src)) f₁, f₂ = first(fusiontrees(src)) a::I = f₂.uncoupled[1] @@ -347,7 +312,7 @@ function foldleft(src::FusionTreeBlock) c::I = f₂.coupled cache₃ = Dict((b, c) => Asymbol(a, b, c)) - U = zeros(eltype(coeff₂), length(dst), length(src)) + U = zeros(eltype(coeff₂), length(src), length(src)) for (col, (f₁, f₂)) in enumerate(fusiontrees(src)) f₂′s, coeffs₂ = get!(cache₂, f₂) do multi_Fmove(f₂) @@ -367,12 +332,12 @@ function foldleft(src::FusionTreeBlock) if isduala coeff *= conj(κₐ) end - row = indexmap[treeindex_data((f₁′, f₂′))] + row = row!(rows, (f₁′, f₂′)) @inbounds U[row, col] += coeff end end end - return dst => U + return FusionTreeBlock(rows.trees) => U end # clockwise cyclic permutation while preserving (N₁, N₂): foldright & bendleft @@ -482,7 +447,7 @@ function _repartition_body(N) if N == 0 ex = quote T = fusionscalartype(sectortype(src)) - if FusionStyle(src) === UniqueFusion() + if src isa FusionTreePair return src => one(T) else U = copyto!(zeros(T, length(src), length(src)), LinearAlgebra.I) diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index 4ea3ab915..dcc63ae42 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -117,6 +117,10 @@ const FusionTreePair{I, N₁, N₂} = Tuple{FusionTree{I, N₁}, FusionTree{I, N Collection of fusion-splitting tree pairs that share the same uncoupled sectors. Mostly internal structure to speed up non-`UniqueFusion` fusiontree manipulation computations. +Manipulating a block requires it to be complete, i.e. to contain all the fusion-splitting tree +pairs with its uncoupled sectors and duality flags, as `FusionTreeBlock{I}(uncoupled, isdual)` +does. The resulting block holds its trees in the order in which they are produced, rather than in +the canonical order of `FusionTreeBlock{I}(uncoupled, isdual)`. """ struct FusionTreeBlock{I, N₁, N₂, F <: FusionTreePair{I, N₁, N₂}} trees::Vector{F} @@ -278,9 +282,36 @@ function treeindex_data((f₁, f₂)) 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))) + +# The trees of the destination block of a block manipulation, in the order in which they are +# produced, and the row of each of them in the transformation matrix. The source block has to be +# complete, such that the destination block has as many trees (`n`). +struct TreeRows{F, D} + trees::Vector{F} + index::Union{D, Nothing} # `nothing`: few trees, the rows are found by a linear search + n::Int + # no default outer constructor, of which `D` would be unbound if `index === nothing` + TreeRows{F, D}(trees::Vector{F}, index, n::Int) where {F, D} = new{F, D}(trees, index, n) +end +function TreeRows(::Type{F}, n::Int) where {F} + trees = sizehint!(Vector{F}(), n) + K = Base.promote_op(treeindex_data, F) + index = n > 32 ? sizehint!(Dict{K, Int}(), n) : nothing + return TreeRows{F, Dict{K, Int}}(trees, index, n) +end +function row!(rows::TreeRows, f) + key = treeindex_data(f) + i = if isnothing(rows.index) + findfirst(f′ -> treeindex_data(f′) == key, rows.trees) + else + get(rows.index, key, nothing) + end + isnothing(i) || return i + length(rows.trees) < rows.n || + throw(ArgumentError("manipulating a block requires it to be complete")) + push!(rows.trees, f) + isnothing(rows.index) || (rows.index[key] = length(rows.trees)) + return length(rows.trees) end diff --git a/test/setup.jl b/test/setup.jl index dca3a6c4e..d44366ae9 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -9,12 +9,12 @@ export sectorlist, fast_sectorlist export default_spacelist, factorization_spacelist, ad_spacelist export VIBM, VIBMRepA4 export test_ad_rrule -export _isunitary, _isone +export _isunitary, _isone, _inorder using Random using Test: @test using TensorKit -using TensorKit: ℙ, PlanarTrivial +using TensorKit: ℙ, PlanarTrivial, FusionTreeBlock using TensorKitSectors using TensorOperations: IndexTuple, Index2Tuple using Base.Iterators: take, product @@ -149,6 +149,15 @@ _isunitary(x::Number; kwargs...) = isapprox(x * x', one(x); kwargs...) _isunitary(x; kwargs...) = isunitary(x; kwargs...) _isone(x; kwargs...) = isapprox(x, one(x); kwargs...) +# the rows of the transformation matrix `U` onto the trees of `dst`, in the order of the trees of +# `ref`, as block manipulations keep their trees in the order in which they are produced +_inorder(U::Number, dst, ref) = (@test dst == ref; U) +function _inorder(U::AbstractMatrix, dst::FusionTreeBlock, ref::FusionTreeBlock) + rows = indexin(fusiontrees(ref), fusiontrees(dst)) + @test length(dst) == length(ref) && !any(isnothing, rows) + return U[something.(rows, 1), :] +end + uniquefusionsectorlist = ( Z2Irrep, Z3Irrep, Z4Irrep, Z3Irrep ⊠ Z4Irrep, U1Irrep, diff --git a/test/symmetries/doubletree.jl b/test/symmetries/doubletree.jl index 0af5db19c..1d51b6730 100644 --- a/test/symmetries/doubletree.jl +++ b/test/symmetries/doubletree.jl @@ -57,15 +57,13 @@ using TensorKitSectors # single bend dst, U = @constinferred TK.bendright(src) dst2, U2 = @constinferred TK.bendleft(dst) - @test src == dst2 - @test _isone(U2 * U) + @test _isone(_inorder(U2 * U, dst2, src)) # double bend dst1, U1 = @constinferred TK.bendleft(src) dst2, U2 = @constinferred TK.bendleft(dst1) dst3, U3 = @constinferred TK.bendright(dst2) dst4, U4 = @constinferred TK.bendright(dst3) - @test src == dst4 - @test _isone(U4 * U3 * U2 * U1) + @test _isone(_inorder(U4 * U3 * U2 * U1, dst4, src)) if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) all_inds = (ntuple(identity, numout(src))..., reverse(ntuple(i -> i + numout(src), numin(src)))...) @@ -88,15 +86,13 @@ using TensorKitSectors # single bend dst, U = @constinferred TK.foldleft(src) dst2, U2 = @constinferred TK.foldright(dst) - @test src == dst2 - @test _isone(U2 * U) + @test _isone(_inorder(U2 * U, dst2, src)) # double bend dst1, U1 = @constinferred TK.foldright(src) dst2, U2 = @constinferred TK.foldright(dst1) dst3, U3 = @constinferred TK.foldleft(dst2) dst4, U4 = @constinferred TK.foldleft(dst3) - @test src == dst4 - @test _isone(U4 * U3 * U2 * U1) + @test _isone(_inorder(U4 * U3 * U2 * U1, dst4, src)) if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) all_inds = TupleTools.circshift((ntuple(identity, numout(src))..., reverse(ntuple(i -> i + numout(src), numin(src)))...), -2) @@ -121,7 +117,7 @@ using TensorKitSectors # @test _isunitary(U) dst′, U′ = repartition(dst, N) - @test _isone(U * U′) + @test _isone(U * _inorder(U′, dst′, src)) if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) all_inds = (ntuple(identity, numout(src))..., reverse(ntuple(i -> i + numout(src), numin(src)))...) @@ -152,11 +148,11 @@ using TensorKitSectors dst, U = @constinferred transpose(src, (p1, p2)) dst′, U′ = @constinferred transpose(dst, (ip1, ip2)) - @test _isone(U * U′) + @test _isone(U * _inorder(U′, dst′, src)) if BraidingStyle(I) isa Bosonic dst″, U″ = permute(src, (p1, p2)) - @test U″ ≈ U + @test _inorder(U″, dst″, dst) ≈ U end if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) @@ -217,7 +213,7 @@ using TensorKitSectors else dst′, U′ = @constinferred TensorKit.braid(dst, (ip1, ip2), (il1, il2)) end - @test _isone(U * U′) + @test _isone(U * _inorder(U′, dst′, src)) # check fusiontensor compatibility if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) @@ -235,6 +231,31 @@ using TensorKitSectors end end + FusionStyle(I) isa UniqueFusion || !(BraidingStyle(I) isa HasBraiding) || + @testset "Double fusion tree: braiding a block with a single tree" begin + a, b = random_fusion(I, Val(2)) + cs = filter(c -> Nsymbol(a, b, c) == 1, collect(a ⊗ b)) + if !isempty(cs) + src1 = FusionTreeBlock{I}(((a, b), (rand(cs),)), ((rand(Bool), rand(Bool)), (rand(Bool),))) + @test length(src1) == 1 + for n in 0:3 + p = (randperm(3)...,) + p1, p2 = p[1:n], p[(n + 1):3] + ip = invperm(p) + ilevels = TupleTools.getindices((1, 2, 3), p) + dst, U = @constinferred TensorKit.braid(src1, (p1, p2), ((1, 2), (3,))) + @test length(dst) == 1 && size(U) == (1, 1) + dst′, U′ = @constinferred TensorKit.braid(dst, (ip[1:2], ip[3:3]), (ilevels[1:n], ilevels[(n + 1):3])) + @test dst′ == src1 + @test _isone(U * U′) + if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) + A1 = fusiontensor(only(fusiontrees(src1))) + @test permutedims(A1, (p1..., p2...)) ≈ only(U) * fusiontensor(only(fusiontrees(dst))) + end + end + end + end + UnitStyle(I) isa SimpleUnit && @testset "Double fusion tree: blocks of a HomSpace" begin V = Vect[I](Dict(s => 1 for s in (out..., incoming))) W = V ⊗ V' ← V ⊗ V' # never empty: the unit couples both sides