diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md index b91e0e68b..a6ab14c3b 100644 --- a/docs/src/Changelog.md +++ b/docs/src/Changelog.md @@ -43,6 +43,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - Recoupling matrices are stored such that they can leverage BLAS, also for mixed complex tensor with real recoupling coefficients. - Constructing a `FusionTreeBlock` only takes the fusion product of the shorter side, and `fusionblocks` groups the cached fusion trees of the space instead of constructing a block for every combination of sectors (internal) +- 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)) ## [0.17.2](https://github.com/QuantumKitHub/TensorKit.jl/compare/v0.17.1...v0.17.2) - 2026-09-20 diff --git a/ext/TensorKitEnzymeExt/utility.jl b/ext/TensorKitEnzymeExt/utility.jl index 0994bd850..bb9a21db0 100644 --- a/ext/TensorKitEnzymeExt/utility.jl +++ b/ext/TensorKitEnzymeExt/utility.jl @@ -147,6 +147,7 @@ end @inline EnzymeRules.inactive(::typeof(TensorKit.artin_braid), ::Any, ::Any) = nothing @inline EnzymeRules.inactive(::typeof(TensorKit.fstranspose), ::Any) = nothing @inline EnzymeRules.inactive(::typeof(TensorKit.fstranspose), ::Any, ::Any) = nothing +@inline EnzymeRules.inactive(::typeof(TensorKit.product_transform), ::Any...) = nothing @inline EnzymeRules.inactive(::typeof(TensorKit.treetransposer), ::Any...) = nothing @inline EnzymeRules.inactive(::typeof(TensorKit.treebraider), ::Any...) = nothing @inline EnzymeRules.inactive(::typeof(TensorKit.insertleftunit), ::HomSpace, ::Any) = nothing diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl index 368e326a7..0ad6353ed 100644 --- a/src/fusiontrees/braiding_manipulations.jl +++ b/src/fusiontrees/braiding_manipulations.jl @@ -286,6 +286,13 @@ function braid(src::Union{FusionTreePair, FusionTreeBlock}, p::Index2Tuple, leve @assert TupleTools.isperm((p[1]..., p[2]...)) return fsbraid((src, p, levels)) end +function braid(src::ProductSectorBlock, p::Index2Tuple, levels::Index2Tuple) + isempty(src) && return fsbraid((src, p, levels)) + # the return type is asserted since inference gives up on the recursion on older Julia versions + return product_transform(sectorscalartype(sectortype(src)), src) do x + return braid(x, p, levels) + end::_fsdicttype(typeof((src, p, levels))) +end const FSPBraidKey{I, N₁, N₂} = Tuple{FusionTreePair{I}, Index2Tuple{N₁, N₂}, Index2Tuple} const FSBBraidKey{I, N₁, N₂} = Tuple{FusionTreeBlock{I}, Index2Tuple{N₁, N₂}, Index2Tuple} diff --git a/src/fusiontrees/duality_manipulations.jl b/src/fusiontrees/duality_manipulations.jl index e399f87ab..cae44aca3 100644 --- a/src/fusiontrees/duality_manipulations.jl +++ b/src/fusiontrees/duality_manipulations.jl @@ -526,6 +526,13 @@ function Base.transpose(src::Union{FusionTreePair, FusionTreeBlock}, p::Index2Tu iscyclicpermutation(p′) || throw(ArgumentError("invalid cyclic or planar permutation p = $p")) return fstranspose((src, p)) end +function Base.transpose(src::ProductSectorBlock, p::Index2Tuple) + isempty(src) && return fstranspose((src, p)) + # the return type is asserted since inference gives up on the recursion on older Julia versions + return product_transform(fusionscalartype(sectortype(src)), src) do x + return transpose(x, p) + end::_fsdicttype(typeof((src, p))) +end const FSPTransposeKey{I, N₁, N₂} = Tuple{FusionTreePair{I}, Index2Tuple{N₁, N₂}} const FSBTransposeKey{I, N₁, N₂} = Tuple{FusionTreeBlock{I}, Index2Tuple{N₁, N₂}} diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index 40a2961b6..ea8334462 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -127,6 +127,13 @@ function FusionTreeBlock{I}( isdual::Tuple{NTuple{N₁, Bool}, NTuple{N₂, Bool}}; sizehint::Int = 0 ) where {I <: Sector, N₁, N₂} + return fusiontreeblock(I, uncoupled, isdual; sizehint) +end +# construct the block by enumerating its coupled sectors and fusion trees +function fusiontreeblock( + ::Type{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) @@ -379,6 +386,9 @@ end # Fusion tree iterators include("iterator.jl") +# Fusion trees of product sectors +include("producttrees.jl") + # Manipulate fusion trees include("basic_manipulations.jl") include("duality_manipulations.jl") diff --git a/src/fusiontrees/producttrees.jl b/src/fusiontrees/producttrees.jl new file mode 100644 index 000000000..0fba6aec4 --- /dev/null +++ b/src/fusiontrees/producttrees.jl @@ -0,0 +1,162 @@ +# Product sectors +# --------------- +# Fusion tree manipulations of a product sector are computed from those of two parts of the product, +# as the Kronecker product of their recoupling matrices. + +# products of at least two sectors, also with names +const _MultiSectorTuple = Tuple{Sector, Sector, Vararg{Sector}} +const ProductSectorType = Union{ + ProductSector{<:_MultiSectorTuple}, NamedSector{<:NamedTuple{<:Any, <:_MultiSectorTuple}}, +} +const ProductSectorBlock = FusionTreeBlock{<:ProductSectorType} + +# The two parts of a product are its non-`UniqueFusion` and its `UniqueFusion` sectors if it has both, +# and otherwise its first and its remaining sectors, on which we recurse. +Base.@assume_effects :foldable function _partpositions(::Type{ProductSector{T}}) where {T} + positions = ntuple(identity, fieldcount(T)) + isunique = map(i -> FusionStyle(fieldtype(T, i)) isa UniqueFusion, positions) + any(isunique) && !all(isunique) && + return Base.filter(i -> !isunique[i], positions), Base.filter(i -> isunique[i], positions) + return (first(positions),), Base.tail(positions) +end +# parts with a single sector are not wrapped in a `ProductSector`, to share the caches of that sector +_sectorproduct(sectors::Tuple{Sector}) = only(sectors) +_sectorproduct(sectors::Tuple) = ProductSector(sectors) +_sectortuple(s::ProductSector) = s.sectors +_sectortuple(s::Sector) = (s,) + +function deconstructsector(s::ProductSector) + positions₁, positions₂ = _partpositions(typeof(s)) + return _sectorproduct(TupleTools.getindices(s.sectors, positions₁)), + _sectorproduct(TupleTools.getindices(s.sectors, positions₂)) +end +Base.@assume_effects :foldable function deconstructsectortype(::Type{ProductSector{T}}) where {T} + return map(_partpositions(ProductSector{T})) do positions + return length(positions) == 1 ? fieldtype(T, only(positions)) : + ProductSector{Tuple{map(i -> fieldtype(T, i), positions)...}} + end +end +function reconstructsector(::Type{ProductSector{T}}, s₁::Sector, s₂::Sector) where {T} + sectors = (_sectortuple(s₁)..., _sectortuple(s₂)...) + positions = TupleTools.invperm(TupleTools.vcat(_partpositions(ProductSector{T})...)) + return ProductSector{T}(TupleTools.getindices(sectors, positions)) +end + +# named products are split and joined as the corresponding `ProductSector` +_productsectortype(::Type{NamedSector{NamedTuple{K, T}}}) where {K, T} = ProductSector{T} +deconstructsector(s::NamedSector) = deconstructsector(ProductSector(s)) +deconstructsectortype(::Type{I}) where {I <: NamedSector} = deconstructsectortype(_productsectortype(I)) +reconstructsector(::Type{I}, s₁::Sector, s₂::Sector) where {I <: NamedSector} = + I(reconstructsector(_productsectortype(I), s₁, s₂)) + +# the vertex labels of a product are the linear indices of the labels of its parts, with the second +# part running fastest, such that a `UniqueFusion` second part leaves the labels of the first unchanged +_vertexlines(f::FusionTree) = (f.uncoupled[1], f.innerlines..., f.coupled) +_vertexmultiplicities(f::FusionTree) = + ntuple(i -> Nsymbol(_vertexlines(f)[i], f.uncoupled[i + 1], _vertexlines(f)[i + 1]), length(f) - 1) + +_maptree(proj, f::FusionTree, vertices) = FusionTree( + map(proj, f.uncoupled), proj(f.coupled), f.isdual, map(proj, f.innerlines), vertices +) +function deconstructtree(f::FusionTree) + vertices₁ = vertices₂ = f.vertices + if FusionStyle(last(deconstructsectortype(sectortype(f)))) isa UniqueFusion + vertices₂ = map(Returns(1), f.vertices) + elseif FusionStyle(f) isa GenericFusion && length(f) > 1 + ns = _vertexmultiplicities(_maptree(last ∘ deconstructsector, f, map(Returns(1), f.vertices))) + vertices₁₂ = map((μ, n) -> (divrem(μ - 1, n) .+ 1), f.vertices, ns) + vertices₁ = first.(vertices₁₂) + vertices₂ = last.(vertices₁₂) + end + return _maptree(first ∘ deconstructsector, f, vertices₁), + _maptree(last ∘ deconstructsector, f, vertices₂) +end +function reconstructtree(::Type{I}, f₁::FusionTree, f₂::FusionTree) where {I} + vertices = if FusionStyle(I) isa GenericFusion && length(f₁) > 1 && !(FusionStyle(f₂) isa UniqueFusion) + map((μ, ν, n) -> ν + (μ - 1) * n, f₁.vertices, f₂.vertices, _vertexmultiplicities(f₂)) + else + f₁.vertices + end + join(a, b) = reconstructsector(I, a, b) + return fusiontreetype(I, length(f₁))( + map(join, f₁.uncoupled, f₂.uncoupled), join(f₁.coupled, f₂.coupled), f₁.isdual, + map(join, f₁.innerlines, f₂.innerlines), vertices + ) +end + +# the canonical order of a product block is that of the Kronecker product of the blocks of its +# parts, with the first part running fastest +function productblock(::Type{I}, block₁::FusionTreeBlock, block₂::FusionTreeBlock) where {I} + trees = fusiontreetype(I, numout(block₁), numin(block₁))[] + sizehint!(trees, length(block₁) * length(block₂)) + for (g₁, g₂) in fusiontrees(block₂) + for (f₁, f₂) in fusiontrees(block₁) + push!(trees, (reconstructtree(I, f₁, g₁), reconstructtree(I, f₂, g₂))) + end + end + return FusionTreeBlock(trees) +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 <: ProductSectorType, N₁, N₂} + # blocks of `UniqueFusion` sectors are enumerated directly + FusionStyle(I) isa UniqueFusion && return fusiontreeblock(I, uncoupled, isdual; sizehint) + + # the second part holds the `UniqueFusion` sectors of a mixed product, which are cheaper to construct + I₁, I₂ = deconstructsectortype(I) + block₂ = FusionTreeBlock{I₂}(map.(last ∘ deconstructsector, uncoupled), isdual) + isempty(block₂) && return FusionTreeBlock(Vector{fusiontreetype(I, N₁, N₂)}()) + + block₁ = FusionTreeBlock{I₁}(map.(first ∘ deconstructsector, uncoupled), isdual) + return productblock(I, block₁, block₂) +end + +# the distinct tree pairs, and the position of every tree pair among them +function _distincttrees(trees) + distinct = unique(trees) + return distinct, indexin(trees, distinct) +end +# the tree pairs of the two parts of a tree pair +_deconstructtrees((f₁, f₂)) = map(tuple, deconstructtree(f₁), deconstructtree(f₂)) + +# apply `transform` to the block `src`, where every tree pair of a `UniqueFusion` block maps to a +# single tree pair +function _unique_transform(transform::F, ::Type{T}, src::FusionTreeBlock) where {F, T} + FusionStyle(src) isa UniqueFusion || return transform(src) + dst = map(transform, fusiontrees(src)) + return FusionTreeBlock(map(first, dst)) => LinearAlgebra.diagm(convert.(T, last.(dst))) +end + +# apply `transform` to the block `src` through its parts, with the recoupling matrix with scalar type `T` +function product_transform(transform::F, ::Type{T}, src::FusionTreeBlock{I}) where {F, T, I} + if FusionStyle(I) isa UniqueFusion + return _unique_transform(transform, T, src) + elseif FusionStyle(last(deconstructsectortype(I))) isa UniqueFusion && numind(src) > 0 + # the second part holds the `UniqueFusion` sectors of a mixed product, with a single tree pair + (g₁, g₂), coeff = transform(last(_deconstructtrees(first(fusiontrees(src))))) + dst₁, U₁ = transform(FusionTreeBlock(map(first ∘ _deconstructtrees, fusiontrees(src)))) + dst = FusionTreeBlock( + map(((f₁, f₂),) -> (reconstructtree(I, f₁, g₁), reconstructtree(I, f₂, g₂)), fusiontrees(dst₁)) + ) + U = Matrix{T}(undef, size(U₁)) + U .= coeff .* U₁ + return dst => U + end + src₁, cols₁ = _distincttrees(map(first ∘ _deconstructtrees, fusiontrees(src))) + src₂, cols₂ = _distincttrees(map(last ∘ _deconstructtrees, fusiontrees(src))) + dst₁, U₁ = _unique_transform(transform, T, FusionTreeBlock(src₁)) + dst₂, U₂ = _unique_transform(transform, T, FusionTreeBlock(src₂)) + + dst = productblock(I, dst₁, dst₂) + U = Matrix{T}(undef, length(dst), length(src)) + for (col, i, j) in zip(axes(U, 2), cols₁, cols₂) + row = 0 + for j₂ in axes(U₂, 1), j₁ in axes(U₁, 1) + U[row += 1, col] = U₁[j₁, i] * U₂[j₂, j] + end + end + return dst => U +end diff --git a/test/setup.jl b/test/setup.jl index 2030228a5..ef8f1896f 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -158,7 +158,8 @@ uniquefusionsectorlist = ( ) simplefusionsectorlist = ( CU1Irrep, SU2Irrep, FibonacciAnyon, IsingAnyon, - FermionParity ⊠ U1Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, + @NamedSector{f::FermionParity, s::SU2Irrep, c::U1Irrep}, # named Hubbard-like + FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, ZNElement{6, 3} ⊠ SU2Irrep, # complex F symbols, anyonic braiding Z3Element{1} ⊠ FibonacciAnyon ⊠ FibonacciAnyon, # complex F symbols, no braiding ) diff --git a/test/symmetries/doubletree.jl b/test/symmetries/doubletree.jl index 180e19bdc..38d0aaea0 100644 --- a/test/symmetries/doubletree.jl +++ b/test/symmetries/doubletree.jl @@ -12,7 +12,7 @@ using TensorKitSectors @timedtestset "Double fusion trees for $(TensorKit.type_repr(I))" verbose = true for I in (fast_tests ? fast_sectorlist : sectorlist) Istr = TensorKit.type_repr(I) - N = I <: ProductSector ? 3 : 4 + N = I <: Union{ProductSector, NamedSector} ? 3 : 4 if UnitStyle(I) isa SimpleUnit out = random_fusion(I, Val(N))