incorporate code suggestions

Author: lkdvos

docs/make.jl

@@ -50,6 +50,9 @@ mathengine = MathJax3(
     )
 )
 
+# docstrings don't need `using TensorKit`
+DocMeta.setdocmeta!(TensorKit, :DocTestSetup, :(using TensorKit); recursive = true)
+
 makedocs(;
     modules = [TensorKit, TensorKitSectors],
     sitename = "TensorKit.jl",

docs/src/lib/fusiontrees.md

@@ -25,7 +25,7 @@ split
 TensorKit.join
 merge
 elementary_trace
-planar_trace(f::FusionTree{I, N}, q::Index2Tuple{N₃, N₃}) where {I, N, N₃}
+planar_trace(f::FusionTree, q::Index2Tuple)
 artin_braid
 braid(f::FusionTree{I, N}, p::IndexTuple{N}, levels::IndexTuple{N}) where {I, N}
 permute(f::FusionTree{I, N}, p::IndexTuple{N}) where {I, N}

docs/src/man/fusiontrees.md

@@ -124,7 +124,7 @@ For this, we provide an interface
 
 [`braid(f::FusionTree{I, N}, p::IndexTuple{N}, levels::IndexTuple{N})`](@ref braid(::FusionTree{I, N}, ::IndexTuple{N}, ::IndexTuple{N}) where {I, N})
 
-where the braid is specified as a permutation, such that the new sector at position `i` was originally at position `permutation[i]`, and where every uncoupled sector is also assigned a level or depth.
+where the braid is specified as a permutation, such that the new sector at position `i` was originally at position `p[i]`, and where every uncoupled sector is also assigned a level or depth.
 The permutation is decomposed into swaps between neighbouring sectors, and when two sectors are swapped, their respective level will determine whether the left sector is braided over or under its right neighbor.
 This interface does not allow to specify the most general braid, and in particular will never wind one line around another, but can be used as a more general building block for arbitrary braids than the elementary Artin generators.
 A graphical example makes this probably more clear, i.e for `levels = (1, 2, 3, 4, 5)` and `permutation = (5, 3, 1, 4, 2)`, the corresponding braid is given by

src/fusiontrees/basic_manipulations.jl

@@ -3,8 +3,7 @@
 # -> rewrite generic fusion tree in basis of fusion trees in standard form
 # -> only depend on Fsymbol
 """
-    split(f::FusionTree{I, N}, M::Int)
-    -> (::FusionTree{I, M}, ::FusionTree{I, N-M+1})
+    split(f::FusionTree{I, N}, M::Int) -> (::FusionTree{I, M}, ::FusionTree{I, N - M + 1})
 
 Split a fusion tree into two. The first tree has as uncoupled sectors the first `M`
 uncoupled sectors of the input tree `f`, whereas its coupled sector corresponds to the
@@ -15,6 +14,18 @@ operation is the inverse of `join` in the sense that if `f == join(split(f, M)..
 holds for all `M` between `0` and `N`, where `N` is the number of uncoupled sectors of `f`.
 
 See also [`join`](@ref) and [`insertat`](@ref).
+
+## Examples
+
+```jldoctest
+julia> f = FusionTree{Z2Irrep}((1, 1, 0), 0, (false, false, false));
+
+julia> f₁, f₂ = TensorKit.split(f, 2)
+(FusionTree{Irrep[ℤ₂]}((1, 1), 0, (false, false), ()), FusionTree{Irrep[ℤ₂]}((0, 0), 0, (false, false), ()))
+
+julia> TensorKit.join(f₁, f₂) == f
+true
+```
 """
 @inline function split(f::FusionTree{I, N}, M::Int) where {I, N} # inline helps with constant propagation of M
     0 <= M <= N ||
@@ -52,6 +63,17 @@ operation is the inverse of split, in the sense that `f == join(split(f, M)...)`
 holds for all `M` between `0` and `N`, where `N` is the number of uncoupled sectors of `f`.
 
 See also [`split`](@ref) and [`insertat`](@ref).
+
+## Examples
+
+```jldoctest
+julia> f₁ = FusionTree{Z2Irrep}((1, 1), 0, (false, false));
+
+julia> f₂ = FusionTree{Z2Irrep}((0, 0), 0, (false, false));
+
+julia> f = TensorKit.join(f₁, f₂)
+FusionTree{Irrep[ℤ₂]}((1, 1, 0), 0, (false, false, false), (0,))
+```
 """
 @inline function join(f₁::FusionTree{I, N₁}, f₂::FusionTree{I, N₂}) where {I, N₁, N₂}
     (f₁.coupled == f₂.uncoupled[1] && !f₂.isdual[1]) ||
@@ -574,4 +596,3 @@ function flip((f₁, f₂)::FusionTreePair{I, N₁, N₂}, ind; inv::Bool = fals
     end
     return SingletonDict((f₁′, f₂′) => factor)
 end
-

src/fusiontrees/duality_manipulations.jl

@@ -573,37 +573,39 @@ CacheStyle(::typeof(fstranspose), k::FSBTransposeKey{I}) where {I} =
 # -> composite manipulations that depend on the duality (rigidity) and pivotal structure
 # -> planar manipulations that do not require braiding, everything is in Fsymbol (A/Bsymbol)
 
-function planar_trace(
-        (f₁, f₂)::FusionTreePair{I}, (p1, p2)::Index2Tuple{N₁, N₂}, (q1, q2)::Index2Tuple{N₃, N₃}
-    ) where {I, N₁, N₂, N₃}
-    N = N₁ + N₂ + 2N₃
-    @assert length(f₁) + length(f₂) == N
-    if N₃ == 0
-        return transpose((f₁, f₂), (p1, p2))
+function planar_trace((f₁, f₂)::FusionTreePair, (p₁, p₂)::Index2Tuple, (q₁, q₂)::Index2Tuple)
+    length(q₁) == length(q₂) ||
+        throw(ArgumentError(lazy"trace index tuples q₁ and q₂ must have equal length, got $(length(q₁)) and $(length(q₂))"))
+    I = sectortype(f₁)
+    N = length(p₁) + length(p₂) + 2 * length(q₁)
+    length(f₁) + length(f₂) == N ||
+        throw(ArgumentError(lazy"fusion tree pair has $(length(f₁) + length(f₂)) indices, but permutation expects $N = $(length(p₁)) + $(length(p₂)) + 2×$(length(q₁))"))
+    if isempty(q₁)
+        return transpose((f₁, f₂), (p₁, p₂))
     end
 
     linearindex = (
         ntuple(identity, Val(length(f₁)))...,
         reverse(length(f₁) .+ ntuple(identity, Val(length(f₂))))...,
     )
 
-    q1′ = TupleTools.getindices(linearindex, q1)
-    q2′ = TupleTools.getindices(linearindex, q2)
-    p1′, p2′ = let q′ = (q1′..., q2′...)
+    q₁′ = TupleTools.getindices(linearindex, q₁)
+    q₂′ = TupleTools.getindices(linearindex, q₂)
+    p₁′, p₂′ = let q′ = (q₁′..., q₂′...)
         (
-            map(l -> l - count(l .> q′), TupleTools.getindices(linearindex, p1)),
-            map(l -> l - count(l .> q′), TupleTools.getindices(linearindex, p2)),
+            map(l -> l - count(l .> q′), TupleTools.getindices(linearindex, p₁)),
+            map(l -> l - count(l .> q′), TupleTools.getindices(linearindex, p₂)),
         )
     end
 
     T = fusionscalartype(I)
-    F₁ = fusiontreetype(I, N₁)
-    F₂ = fusiontreetype(I, N₂)
+    F₁ = fusiontreetype(I, length(p₁))
+    F₂ = fusiontreetype(I, length(p₂))
     newtrees = FusionTreeDict{Tuple{F₁, F₂}, T}()
     if FusionStyle(I) isa UniqueFusion
         (f₁′, f₂′), coeff′ = repartition((f₁, f₂), N)
-        for (f₁′′, coeff′′) in planar_trace(f₁′, (q1′, q2′))
-            (f12′′′, coeff′′′) = transpose((f₁′′, f₂′), (p1′, p2′))
+        for (f₁′′, coeff′′) in planar_trace(f₁′, (q₁′, q₂′))
+            (f12′′′, coeff′′′) = transpose((f₁′′, f₂′), (p₁′, p₂′))
             coeff = coeff′ * coeff′′ * coeff′′′
             iszero(coeff) || (newtrees[f12′′′] = get(newtrees, f12′′′, zero(coeff)) + coeff)
         end
@@ -612,9 +614,9 @@ function planar_trace(
         src = FusionTreeBlock([(f₁, f₂)])
         dst, U = repartition(src, N)
         for ((f₁′, f₂′), coeff′) in zip(fusiontrees(dst), U)
-            for (f₁′′, coeff′′) in planar_trace(f₁′, (q1′, q2′))
+            for (f₁′′, coeff′′) in planar_trace(f₁′, (q₁′, q₂′))
                 src′ = FusionTreeBlock([(f₁′′, f₂′)])
-                dst′, U′ = transpose(src′, (p1′, p2′))
+                dst′, U′ = transpose(src′, (p₁′, p₂′))
                 for (f12′′′, coeff′′′) in zip(fusiontrees(dst′), U′)
                     coeff = coeff′ * coeff′′ * coeff′′′
                     iszero(coeff) || (newtrees[f12′′′] = get(newtrees, f12′′′, zero(coeff)) + coeff)
@@ -626,20 +628,23 @@ function planar_trace(
 end
 
 """
-    planar_trace(f::FusionTree{I,N}, (q1, q2)::Index2Tuple{N₃,N₃}) where {I,N,N₃}
-        -> <:AbstractDict{FusionTree{I,N-2*N₃}, <:Number}
+    planar_trace(f::FusionTree, (q₁, q₂)::Index2Tuple)
+        -> <:AbstractDict{<:FusionTree, <: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.
+Perform a planar trace of the uncoupled indices of the fusion tree `f` at `q₁` with those at `q₂`,
+where `q₁[i]` is connected to `q₂[i]` for all `i`. The result is returned as a dictionary of output
+trees and corresponding coefficients.
 """
-function planar_trace(f::FusionTree{I, N}, (q1, q2)::Index2Tuple{N₃, N₃}) where {I, N, N₃}
+function planar_trace(f::FusionTree, (q₁, q₂)::Index2Tuple)
+    length(q₁) == length(q₂) ||
+        throw(ArgumentError(lazy"trace index tuples q₁ and q₂ must have equal length, got $(length(q₁)) and $(length(q₂))"))
+    I = sectortype(f)
     T = fusionscalartype(I)
-    F = fusiontreetype(I, N - 2 * N₃)
+    F = fusiontreetype(I, length(f) - 2 * length(q₁))
     newtrees = FusionTreeDict{F, T}()
-    N₃ === 0 && return push!(newtrees, f => one(T))
+    isempty(q₁) && return push!(newtrees, f => one(T))
 
-    for (i, j) in zip(q1, q2)
+    for (i, j) in zip(q₁, q₂)
         (f.uncoupled[i] == dual(f.uncoupled[j]) && f.isdual[i] != f.isdual[j]) ||
             return newtrees
     end
@@ -649,29 +654,30 @@ function planar_trace(f::FusionTree{I, N}, (q1, q2)::Index2Tuple{N₃, N₃}) wh
     # tracing away neighbouring pairs.
     k = 1
     local i, j
-    while k <= N₃
-        if mod1(q1[k] + 1, N) == q2[k]
-            i = q1[k]
-            j = q2[k]
+    while k <= length(q₁)
+        if mod1(q₁[k] + 1, length(f)) == q₂[k]
+            i = q₁[k]
+            j = q₂[k]
             break
-        elseif mod1(q2[k] + 1, N) == q1[k]
-            i = q2[k]
-            j = q1[k]
+        elseif mod1(q₂[k] + 1, length(f)) == q₁[k]
+            i = q₂[k]
+            j = q₁[k]
             break
         else
             k += 1
         end
     end
-    k > N₃ && throw(ArgumentError("Not a planar trace"))
+    k > length(q₁) &&
+        throw(ArgumentError(lazy"indices $q₁ and $q₂ do not form a valid planar trace on a fusion tree with $(length(f)) legs: no neighboring pair found among the remaining trace indices"))
 
-    q1′ = let i = i, j = j
-        map(l -> (l - (l > i) - (l > j)), TupleTools.deleteat(q1, k))
+    q₁′ = let i = i, j = j
+        map(l -> (l - (l > i) - (l > j)), TupleTools.deleteat(q₁, k))
     end
-    q2′ = let i = i, j = j
-        map(l -> (l - (l > i) - (l > j)), TupleTools.deleteat(q2, k))
+    q₂′ = let i = i, j = j
+        map(l -> (l - (l > i) - (l > j)), TupleTools.deleteat(q₂, 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′, (q₁′, q₂′))
             coeff = coeff′ * coeff′′
             if !iszero(coeff)
                 newtrees[f′′] = get(newtrees, f′′, zero(coeff)) + coeff

src/spaces/homspace.jl

@@ -141,7 +141,7 @@ fusiontrees(W::HomSpace) = fusionblockstructure(W).fusiontreelist
 """
     fusionblocks(W::HomSpace)
 
-Return the fusionblocks corresponding to all valid fusion channels of a given `HomSpace`,
+Return the [`FusionTreeBlock`](@ref)s corresponding to all valid fusion channels of a given `HomSpace`,
 grouped by their uncoupled charges.
 """
 function fusionblocks(W::HomSpace)

src/tensors/tensoroperations.jl

@@ -190,23 +190,20 @@ TO.tensorcost(t::AbstractTensorMap, i::Int) = dim(space(t, i))
 """
     trace_permute!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap,
                    (p₁, p₂)::Index2Tuple, (q₁, q₂)::Index2Tuple,
-                   α::Number, β::Number, backend=TO.DefaultBackend())
+                   α::Number, β::Number, backend = TO.DefaultBackend())
 
 Return the updated `tdst`, which is the result of adding `α * tsrc` to `tdst` after permuting
 the indices of `tsrc` according to `(p₁, p₂)` and furthermore tracing the indices in `q₁` and `q₂`.
 """
 function trace_permute!(
         tdst::AbstractTensorMap,
-        tsrc::AbstractTensorMap,
-        (p₁, p₂)::Index2Tuple,
-        (q₁, q₂)::Index2Tuple,
-        α::Number,
-        β::Number,
-        backend = TO.DefaultBackend()
+        tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple, (q₁, q₂)::Index2Tuple,
+        α::Number, β::Number, backend = TO.DefaultBackend()
     )
     # some input checks
     S = check_spacetype(tdst, tsrc)
-    if !(BraidingStyle(sectortype(S)) isa SymmetricBraiding)
+    I = sectortype(S)
+    if !(BraidingStyle(I) isa SymmetricBraiding)
         throw(SectorMismatch("only tensors with symmetric braiding rules can be contracted; try `@planar` instead"))
     end
     (N₃ = length(q₁)) == length(q₂) ||
@@ -223,60 +220,61 @@ function trace_permute!(
                     q₁ = $(q₁), q₂ = $(q₂)"))
     end
 
-    I = sectortype(S)
     if I === Trivial
-        cod = codomain(tsrc)
-        dom = domain(tsrc)
-        n = length(cod)
         TO.tensortrace!(tdst[], tsrc[], (p₁, p₂), (q₁, q₂), false, α, β, backend)
-        return tdst
+    else
+        _trace_permute!(FusionStyle(I), tdst, tsrc, (p₁, p₂), (q₁, q₂), α, β, backend)
     end
 
-    cod = codomain(tsrc)
-    dom = domain(tsrc)
-    n = length(cod)
+    return tdst
+end
+
+function _trace_permute!(::UniqueFusion, tdst, tsrc, (p₁, p₂), (q₁, q₂), α, β, backend)
     scale!(tdst, β)
-    r₁ = (p₁..., q₁...)
-    r₂ = (p₂..., q₂...)
-
-    # TODO: Is it worth to add multithreading support?
-    if FusionStyle(I) isa UniqueFusion
-        for (f₁, f₂) in fusiontrees(tsrc)
-            (f₁′, f₂′), coeff = permute((f₁, f₂), (r₁, r₂))
-            f₁′′, g₁ = split(f₁′, N₁)
-            f₂′′, g₂ = split(f₂′, N₂)
-            g₁ == g₂ || continue
-            coeff *= dim(g₁.coupled) / dim(g₁.uncoupled[1])
-            for i in 2:length(g₁.uncoupled)
-                if !(g₁.isdual[i])
-                    coeff *= twist(g₁.uncoupled[i])
-                end
+    r₁, r₂ = (p₁..., q₁...), (p₂..., q₂...)
+
+    for (f₁, f₂) in fusiontrees(tsrc)
+        (f₁′, f₂′), coeff = permute((f₁, f₂), (r₁, r₂))
+        f₁′′, g₁ = split(f₁′, N₁)
+        f₂′′, g₂ = split(f₂′, N₂)
+        g₁ == g₂ || continue
+        coeff *= dim(g₁.coupled) / dim(g₁.uncoupled[1])
+        for i in 2:length(g₁.uncoupled)
+            if !(g₁.isdual[i])
+                coeff *= twist(g₁.uncoupled[i])
             end
-            C = tdst[f₁′′, f₂′′]
-            A = tsrc[f₁, f₂]
-            α′ = α * coeff
-            TO.tensortrace!(C, A, (p₁, p₂), (q₁, q₂), false, α′, One(), backend)
         end
-    else
-        for src in fusionblocks(tsrc)
-            dst, U = permute(src, (r₁, r₂))
-            for (i, (f₁, f₂)) in enumerate(fusiontrees(src))
-                for (j, (f₁′, f₂′)) in enumerate(fusiontrees(dst))
-                    coeff = U[j, i]
-                    f₁′′, g₁ = split(f₁′, N₁)
-                    f₂′′, g₂ = split(f₂′, N₂)
-                    g₁ == g₂ || continue
-                    coeff *= dim(g₁.coupled) / dim(g₁.uncoupled[1])
-                    for i in 2:length(g₁.uncoupled)
-                        if !(g₁.isdual[i])
-                            coeff *= twist(g₁.uncoupled[i])
-                        end
+        C = tdst[f₁′′, f₂′′]
+        A = tsrc[f₁, f₂]
+        α′ = α * coeff
+        TO.tensortrace!(C, A, (p₁, p₂), (q₁, q₂), false, α′, One(), backend)
+    end
+
+    return tdst
+end
+
+function _trace_permute!(::FusionStyle, tdst, tsrc, (p₁, p₂), (q₁, q₂), α, β, backend)
+    scale!(tdst, β)
+    r₁, r₂ = (p₁..., q₁...), (p₂..., q₂...)
+
+    for src in fusionblocks(tsrc)
+        dst, U = permute(src, (r₁, r₂))
+        for (i, (f₁, f₂)) in enumerate(fusiontrees(src))
+            for (j, (f₁′, f₂′)) in enumerate(fusiontrees(dst))
+                coeff = U[j, i]
+                f₁′′, g₁ = split(f₁′, N₁)
+                f₂′′, g₂ = split(f₂′, N₂)
+                g₁ == g₂ || continue
+                coeff *= dim(g₁.coupled) / dim(g₁.uncoupled[1])
+                for i in 2:length(g₁.uncoupled)
+                    if !(g₁.isdual[i])
+                        coeff *= twist(g₁.uncoupled[i])
                     end
-                    C = tdst[f₁′′, f₂′′]
-                    A = tsrc[f₁, f₂]
-                    α′ = α * coeff
-                    TO.tensortrace!(C, A, (p₁, p₂), (q₁, q₂), false, α′, One(), backend)
                 end
+                C = tdst[f₁′′, f₂′′]
+                A = tsrc[f₁, f₂]
+                α′ = α * coeff
+                TO.tensortrace!(C, A, (p₁, p₂), (q₁, q₂), false, α′, One(), backend)
             end
         end
     end