IsingBimodule for testing multifusion category within TensorKit

src/TensorKit.jl

@@ -15,7 +15,8 @@ export BraidingStyle, SymmetricBraiding, Bosonic, Fermionic, Anyonic, NoBraiding
 export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep
 export ProductSector
 export FermionParity, FermionNumber, FermionSpin
-export FibonacciAnyon, IsingAnyon
+export FibonacciAnyon, IsingAnyon, IsingBimodule
+export leftone, rightone
 
 export VectorSpace, Field, ElementarySpace # abstract vector spaces
 export InnerProductStyle, NoInnerProduct, HasInnerProduct, EuclideanInnerProduct
@@ -32,7 +33,7 @@ export SpaceMismatch, SectorMismatch, IndexError # error types
 
 # general vector space methods
 export space, field, dual, dim, reduceddim, dims, fuse, flip, isdual, oplus,
-       insertleftunit, insertrightunit, removeunit
+       leftoneunit, rightoneunit, insertleftunit, insertrightunit, removeunit
 
 # partial order for vector spaces
 export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic
@@ -41,7 +42,6 @@ export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic
 export sectortype, sectors, hassector, Nsymbol, Fsymbol, Rsymbol, Bsymbol,
        frobeniusschur, twist, otimes, sectorscalartype
 export fusiontrees, braid, permute, transpose
-export ZNSpace, SU2Irrep, U1Irrep, CU1Irrep
 # other fusion tree manipulations, should not be exported:
 # export insertat, split, merge, repartition, artin_braid,
 #        bendleft, bendright, foldleft, foldright, cycleclockwise, cycleanticlockwise
@@ -229,6 +229,10 @@ include("planar/planaroperations.jl")
 # deprecations: to be removed in version 1.0 or sooner
 include("auxiliary/deprecate.jl")
 
+# Additional methods for IsingBimodule Sector
+# ----------------------------------------
+include("spaces/multifusionspace.jl")
+
 # Extensions
 # ----------
 function __init__()

src/fusiontrees/manipulations.jl

@@ -786,13 +786,13 @@ function elementary_trace(f::FusionTree{I,N}, i) where {I<:Sector,N}
         vertices_ = TupleTools.front(f.vertices)
         f_ = FusionTree(uncoupled_, coupled_, isdual_, inner_, vertices_)
         fs = FusionTree((b,), b, (!f.isdual[1],), (), ())
-        for (f_′, coeff) in merge(fs, f_, unit, 1) # coloring gets reversed here, should be the other unit
+        for (f_′, coeff) in merge(fs, f_, unit, 1)
             f_′.innerlines[1] == unit || continue
             uncoupled′ = Base.tail(Base.tail(f_′.uncoupled))
             isdual′ = Base.tail(Base.tail(f_′.isdual))
             inner′ = N <= 4 ? () : Base.tail(Base.tail(f_′.innerlines))
             vertices′ = N <= 3 ? () : Base.tail(Base.tail(f_′.vertices))
-            f′ = FusionTree(uncoupled′, unit, isdual′, inner′, vertices′) # and this one?
+            f′ = FusionTree(uncoupled′, unit, isdual′, inner′, vertices′)
             coeff *= sqrtdim(b)
             if !(f.isdual[N])
                 coeff *= conj(frobeniusschur(b))
@@ -835,7 +835,7 @@ function artin_braid(f::FusionTree{I,N}, i; inv::Bool=false) where {I<:Sector,N}
     vertices = f.vertices
     oneT = one(sectorscalartype(I))
 
-    if isone(uncoupled[i]) || isone(uncoupled[i + 1])
+    if isone(a) || isone(b)
         # braiding with trivial sector: simple and always possible
         inner′ = inner
         vertices′ = vertices

src/spaces/gradedspace.jl

@@ -132,6 +132,23 @@ function Base.axes(V::GradedSpace{I}, c::I) where {I<:Sector}
 end
 
 Base.oneunit(S::Type{<:GradedSpace{I}}) where {I<:Sector} = S(one(I) => 1)
+
+"""
+    leftoneunit(S::GradedSpace{I}) where {I<:Sector} -> GradedSpace{I}
+
+Return the corresponding vector space of type `GradedSpace{I}` that represents the trivial
+one-dimensional space consisting of the left unit of the objects in  `Sector` `I`.
+"""
+leftoneunit(S::GradedSpace{I}) where {I<:Sector} = oneunit(typeof(S))
+
+"""
+    rightoneunit(S::GradedSpace{I}) where {I<:Sector} -> GradedSpace{I}
+
+Return the corresponding vector space of type `GradedSpace{I}` that represents the trivial
+one-dimensional space consisting of the right unit of the objects in `Sector` `I`.
+"""
+rightoneunit(S::GradedSpace{I}) where {I<:Sector} = oneunit(typeof(S))
+
 Base.zero(S::Type{<:GradedSpace{I}}) where {I<:Sector} = S(one(I) => 0)
 
 # TODO: the following methods can probably be implemented more efficiently for

src/spaces/multifusionspace.jl

@@ -0,0 +1,111 @@
+# additional interface to deal with IsingBimodule Sector
+#------------------------------------------------------------------------------
+
+function dim(V::Vect[IsingBimodule])
+    T = Base.promote_op(*, Int, real(sectorscalartype(sectortype(V))))
+    return reduce(+, dim(V, c) * dim(c) for c in sectors(V);
+                  init=zero(T))
+end
+
+function scalar(t::AbstractTensorMap{T,Vect[IsingBimodule],0,0}) where {T}
+    Bs = collect(blocks(t))
+    inds = findall(!iszero ∘ last, Bs)
+    isempty(inds) && return zero(scalartype(t))
+    return only(last(Bs[only(inds)]))
+end
+
+# no custom fuse: we choose to return empty graded space when fusion is forbidden
+
+function Base.oneunit(S::Vect[IsingBimodule])
+    !isempty(sectors(S)) || throw(ArgumentError("Cannot determine type of empty space"))
+    allequal(a.row for a in sectors(S)) && allequal(a.col for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S are not all equal"))
+    first(sectors(S)).row == first(sectors(S)).col ||
+        throw(ArgumentError("sectors of $S are non-diagonal"))
+    sector = one(first(sectors(S)))
+    return spacetype(S)(sector => 1)
+end
+
+Base.zero(S::Type{Vect[IsingBimodule]}) = Vect[IsingBimodule]()
+
+function blocksectors(W::TensorMapSpace{Vect[IsingBimodule],N₁,N₂}) where {N₁,N₂}
+    codom = codomain(W)
+    dom = domain(W)
+    if N₁ == 0 && N₂ == 0
+        return [IsingBimodule(1, 1, 0), IsingBimodule(2, 2, 0)]
+    elseif N₁ == 0
+        @assert N₂ != 0 "one of Type IsingBimodule doesn't exist"
+        return filter!(isone, collect(blocksectors(dom))) # is this what we want? doesn't allow M/Mop to end at empty space
+    elseif N₂ == 0
+        @assert N₁ != 0 "one of Type IsingBimodule doesn't exist"
+        return filter!(isone, collect(blocksectors(codom)))
+    elseif N₂ <= N₁ # keep intersection
+        return filter!(c -> hasblock(codom, c), collect(blocksectors(dom)))
+    else
+        return filter!(c -> hasblock(dom, c), collect(blocksectors(codom)))
+    end
+end
+
+function rightoneunit(S::Vect[IsingBimodule])
+    !isempty(sectors(S)) || throw(ArgumentError("Cannot determine type of empty space"))
+    allequal(a.col for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S do not have the same rightone"))
+
+    sector = rightone(first(sectors(S)))
+    return spacetype(S)(sector => 1)
+end
+
+function leftoneunit(S::Vect[IsingBimodule])
+    !isempty(sectors(S)) || throw(ArgumentError("Cannot determine type of empty space"))
+    allequal(a.row for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S do not have the same leftone"))
+
+    sector = leftone(first(sectors(S)))
+    return spacetype(S)(sector => 1)
+end
+
+function insertrightunit(P::ProductSpace{Vect[IsingBimodule],N}, ::Val{i};
+                         conj::Bool=false,
+                         dual::Bool=false) where {i,N}
+    i > N && error("cannot insert a sensible right unit onto $P at index $(i+1)")
+    # possible change to rightone of correct space for N = 0
+    u = N > 0 ? rightoneunit(P[i]) : error("no unit object in $P")
+    if dual
+        u = TensorKit.dual(u)
+    end
+    if conj
+        u = TensorKit.conj(u)
+    end
+    return ProductSpace(TupleTools.insertafter(P.spaces, i, (u,)))
+end
+
+# TODO?: overwrite defaults at level of HomSpace and TensorMap?
+function insertleftunit(P::ProductSpace{Vect[IsingBimodule],N}, ::Val{i}; # want no defaults?
+                        conj::Bool=false,
+                        dual::Bool=false) where {i,N}
+    i > N && error("cannot insert a sensible left unit onto $P at index $i") # do we want this to error in the diagonal case?
+    u = N > 0 ? leftoneunit(P[i]) : error("no unit object in $P")
+    if dual
+        u = TensorKit.dual(u)
+    end
+    if conj
+        u = TensorKit.conj(u)
+    end
+    return ProductSpace(TupleTools.insertafter(P.spaces, i - 1, (u,)))
+end
+
+# is this even necessary? can let it error at fusiontrees.jl:93 from the one(IsingBimodule) call
+# but these errors are maybe more informative
+function FusionTree(uncoupled::Tuple{IsingBimodule,Vararg{IsingBimodule}})
+    coupled = collect(⊗(uncoupled...))
+    if length(coupled) == 0 # illegal fusion somewhere
+        throw(ArgumentError("Forbidden fusion with uncoupled sectors $uncoupled"))
+    else # allowed fusions require inner lines
+        error("fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`")
+    end
+end
+
+# this one might also be overkill, since `FusionTreeIterator`s don't check whether the fusion is allowed
+function fusiontrees(uncoupled::Tuple{IsingBimodule,Vararg{IsingBimodule}})
+    return throw(ArgumentError("coupled sector must be provided for IsingBimodule fusion"))
+end

src/tensors/linalg.jl

@@ -526,10 +526,17 @@ function ⊗(t1::AbstractTensorMap, t2::AbstractTensorMap)
         throw(SpaceMismatch("spacetype(t1) ≠ spacetype(t2)"))
     cod1, cod2 = codomain(t1), codomain(t2)
     dom1, dom2 = domain(t1), domain(t2)
-    cod = cod1 ⊗ cod2
-    dom = dom1 ⊗ dom2
+    p12 = ((codomainind(t1)..., (codomainind(t2) .+ numind(t1))...),
+           (domainind(t1)..., (domainind(t2) .+ numind(t1))...))
+
     T = promote_type(scalartype(t1), scalartype(t2))
-    t = zerovector!(similar(t1, T, cod ← dom))
+    TC = promote_type(sectorscalartype(sectortype(t1)), T)
+    t = TO.tensoralloc_contract(TC,
+                                t1, ((codomainind(t1)..., domainind(t1)...), ()), false,
+                                t2, ((), (codomainind(t2)..., domainind(t2)...)), false,
+                                p12, Val(false))
+
+    zerovector!(t)
     for (f1l, f1r) in fusiontrees(t1)
         for (f2l, f2r) in fusiontrees(t2)
             c1 = f1l.coupled # = f1r.coupled

test/fusiontrees.jl

@@ -211,7 +211,7 @@ ti = time()
             end
         end
     end
-    @testset "Fusion tree $Istr: elementy artin braid" begin
+    @testset "Fusion tree $Istr: elementary artin braid" begin
         N = length(out)
         isdual = ntuple(n -> rand(Bool), N)
         for in in ⊗(out...)
@@ -268,7 +268,7 @@ ti = time()
         end
     end
     @testset "Fusion tree $Istr: braiding and permuting" begin
-        f = rand(collect(fusiontrees(out, in, isdual)))
+        f = rand(collect(it))
         p = tuple(randperm(N)...)
         ip = invperm(p)
 
@@ -379,7 +379,7 @@ ti = time()
     f1 = rand(collect(fusiontrees(out, incoming, ntuple(n -> rand(Bool), N))))
     f2 = rand(collect(fusiontrees(out[randperm(N)], incoming, ntuple(n -> rand(Bool), N))))
 
-    @testset "Double fusion tree $Istr: repartioning" begin
+    @testset "Double fusion tree $Istr: repartitioning" begin
         for n in 0:(2 * N)
             d = @constinferred TK.repartition(f1, f2, $n)
             @test dim(incoming) ≈
@@ -548,7 +548,7 @@ ti = time()
     @testset "Double fusion tree $Istr: planar trace" begin
         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])
+        T = sectorscalartype(I)
         d2 = Dict{typeof((f1front, f1front)),T}()
         for ((f1′, f2′), coeff′) in d1
             for ((f1′′, f2′′), coeff′′) in