Factorize fusion tree block manipulations of product sectors

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
 

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

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}

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₂}}

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")

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

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
 )

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))