Simple multifusion example for testing

src/TensorKitSectors.jl

@@ -21,6 +21,7 @@ export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep
 export ProductSector, TimeReversed
 export FermionParity, FermionNumber, FermionSpin
 export PlanarTrivial, FibonacciAnyon, IsingAnyon
+export IsingBimod, CatType, 𝒞, 𝒟, ℳ, ℳᵒᵖ
 
 # unicode exports
 # ---------------
@@ -50,6 +51,7 @@ include("irreps.jl")    # irreps of symmetry groups, with bosonic braiding
 include("product.jl")   # direct product of different sectors
 include("fermions.jl")  # irreps with defined fermionparity and fermionic braiding
 include("anyons.jl")    # non-group sectors
+include("multifusion.jl") # multifusion example, namely Rep Z2 ⊕ Rep Z2 ≅ Ising
 
 # precompile
 # ----------

src/multifusion.jl

@@ -0,0 +1,186 @@
+# Rep Z2 ⊕ Rep Z2 ≅ Ising is worked out here
+# 𝒞 = 𝒟 = RepZ2 ≅ {1, ψ}, while ℳ = Vec ≅ {σ}
+# this is mainly meant for testing within TensorKit without relying on MultiTensorKit
+#------------------------------------------------------------------------------#
+
+@enum CatType 𝒞 = 1 ℳ = 3 ℳᵒᵖ = 2 𝒟 = 4
+# possible TODO: get rid of CatType and use Int instead -> prevents the need to export 𝒞 etc
+function Base.getindex(a::CatType, label::Int)
+    return IsingBimod(a, label)
+end
+
+struct IsingBimod <: Sector
+    type::CatType
+    label::Int
+    function IsingBimod(type::CatType, label::Int)
+        if type == 𝒞
+            0 ≤ label ≤ 1 ||
+                throw(ArgumentError("Invalid 𝒞 label for Ising bimodule: $(label)"))
+        elseif type == ℳ
+            label == 0 ||
+                throw(ArgumentError("Invalid ℳ label for Ising bimodule: $(label)"))
+        elseif type == ℳᵒᵖ
+            label == 0 ||
+                throw(ArgumentError("Invalid ℳᵒᵖ label for Ising bimodule: $(label)"))
+        elseif type == 𝒟
+            0 ≤ label ≤ 1 ||
+                throw(ArgumentError("Invalid 𝒟 label for Ising bimodule: $(label)"))
+        end
+        return new(type, label)
+    end
+end
+
+isC(x::IsingBimod) = x.type == 𝒞
+isM(x::IsingBimod) = x.type == ℳ
+isMop(x::IsingBimod) = x.type == ℳᵒᵖ
+isD(x::IsingBimod) = x.type == 𝒟
+
+function ismodulecategory(a::IsingBimod)
+    return a.type in (ℳ, ℳᵒᵖ)
+end
+
+const all_isingbimod_objects = (IsingBimod(𝒞, 0), IsingBimod(𝒞, 1), IsingBimod(ℳᵒᵖ, 0),
+                                IsingBimod(ℳ, 0), IsingBimod(𝒟, 0), IsingBimod(𝒟, 1))
+
+Base.IteratorSize(::Type{<:SectorValues{IsingBimod}}) = Base.SizeUnknown()
+Base.iterate(::SectorValues{IsingBimod}, i=1) = iterate(all_isingbimod_objects, i)
+Base.length(::SectorValues{IsingBimod}) = length(all_isingbimod_objects)
+
+⊗(a::IsingBimod, b::IsingBimod) = IsingBimodIterator(a, b)
+
+struct IsingBimodIterator
+    a::IsingBimod
+    b::IsingBimod
+end
+
+Base.IteratorSize(::Type{IsingBimodIterator}) = Base.SizeUnknown()
+Base.IteratorEltype(::Type{IsingBimodIterator}) = Base.HasEltype()
+Base.eltype(::Type{IsingBimodIterator}) = IsingBimod
+
+function Base.iterate(iter::IsingBimodIterator, state=0)
+    if isC(iter.a) && isC(iter.b) # 𝒞 × 𝒞 -> 𝒞
+        return state == 0 ? (IsingBimod(𝒞, mod(iter.a.label + iter.b.label, 2)), 1) :
+               nothing
+    end
+
+    if isD(iter.a) && isD(iter.b) # 𝒟 × 𝒟 -> 𝒟
+        return state == 0 ? (IsingBimod(𝒟, mod(iter.a.label + iter.b.label, 2)), 1) :
+               nothing
+    end
+
+    if isM(iter.a) && isMop(iter.b) # ℳ × ℳop -> 𝒞
+        return state < 2 ? (IsingBimod(𝒞, state), state + 1) : nothing
+    end
+
+    if isMop(iter.a) && isM(iter.b) # ℳop × ℳ -> 𝒟
+        return state < 2 ? (IsingBimod(𝒟, state), state + 1) : nothing
+    end
+
+    if isC(iter.a) && isM(iter.b) # 𝒞 × ℳ -> ℳ
+        return state == 0 ? (iter.b, 1) : nothing
+    end
+
+    if isMop(iter.a) && isC(iter.b) # ℳop x 𝒞 -> ℳop
+        return state == 0 ? (iter.a, 1) : nothing
+    end
+
+    if isM(iter.a) && isD(iter.b) # ℳ x 𝒟 -> ℳ
+        return state == 0 ? (iter.a, 1) : nothing
+    end
+
+    if isD(iter.a) && isMop(iter.b) # 𝒟 x ℳop -> ℳop
+        return state == 0 ? (iter.b, 1) : nothing
+    end
+
+    return nothing
+end
+
+function Base.convert(::Type{IsingAnyon}, a::IsingBimod) # identify RepZ2 ⊕ RepZ2 ≅ Ising
+    ismodulecategory(a) && return IsingAnyon(:σ)
+    return IsingAnyon(a.label == 0 ? :I : :ψ)
+end
+
+FusionStyle(::Type{IsingBimod}) = SimpleFusion() # no multiplicities
+BraidingStyle(::Type{IsingBimod}) = NoBraiding() # because of module categories
+
+function Nsymbol(a::IsingBimod, b::IsingBimod, c::IsingBimod)
+    # if a and b can fuse, then so can dual(a) and c, and c and dual(b)
+    # only needs to be explicitly checked when CatTypes differ or when there's a module category involved
+    if a.type != b.type || a.type != c.type || b.type != c.type ||
+       any(ismodulecategory, (a, b, c))
+        c ∈ a ⊗ b && dual(b) ∈ dual(c) ⊗ a && dual(a) ∈ b ⊗ dual(c) ||
+            throw(ArgumentError("invalid fusion channel"))
+    end
+    return Nsymbol(convert(IsingAnyon, a), convert(IsingAnyon, b), convert(IsingAnyon, c))
+end
+
+function Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I<:IsingBimod}
+    Nsymbol(a, b, e) && Nsymbol(e, c, d) &&
+        Nsymbol(b, c, f) && Nsymbol(a, f, d) || return 0.0
+    return Fsymbol(convert(IsingAnyon, a), convert(IsingAnyon, b), convert(IsingAnyon, c),
+                   convert(IsingAnyon, d), convert(IsingAnyon, e), convert(IsingAnyon, f))
+end
+
+function Base.conj(a::IsingBimod) # ℳ ↔ ℳop when conjugating elements within these
+    if isC(a) || isD(a) # self-conjugate within RepZ2
+        return a
+    elseif isM(a)
+        return IsingBimod(ℳᵒᵖ, a.label)
+    else
+        return IsingBimod(ℳ, a.label)
+    end
+end
+
+function rightone(a::IsingBimod)
+    if isC(a) || isD(a)
+        return IsingBimod(a.type, 0)
+    elseif isM(a) # ℳ as right-𝒟 module category
+        return IsingBimod(𝒟, 0)
+    else
+        return IsingBimod(𝒞, 0)
+    end
+end
+
+function leftone(a::IsingBimod)
+    if isC(a) || isD(a)
+        return IsingBimod(a.type, 0)
+    elseif isM(a) # ℳ as left-𝒞 module category
+        return IsingBimod(𝒞, 0)
+    else
+        return IsingBimod(𝒟, 0)
+    end
+end
+
+function Base.one(a::IsingBimod)
+    if isC(a) || isD(a)
+        return IsingBimod(a.type, 0)
+    else
+        throw(DomainError("unit of module category $(a.type) doesn't exist"))
+    end
+end
+
+Base.one(::Type{IsingBimod}) = throw(ArgumentError("one of Type IsingBimod doesn't exist"))
+
+Base.isreal(::Type{IsingBimod}) = false
+
+function Base.show(io::IO, a::IsingBimod)
+    if isC(a)
+        print(io, "𝒞[$(a.label)]")
+    elseif isM(a)
+        print(io, "ℳ[$(a.label)]")
+    elseif isMop(a)
+        print(io, "ℳᵒᵖ[$(a.label)]")
+    elseif isD(a)
+        print(io, "𝒟[$(a.label)]")
+    end
+end
+
+function Base.isless(a::IsingBimod, b::IsingBimod)
+    return isless((a.type, a.label), (b.type, b.label)) # order 𝒞 < ℳᵒᵖ < ℳ < 𝒟, and then within each cat according to label
+end
+
+function Base.hash(a::IsingBimod, h::UInt)
+    return hash(a.label, hash(a.type, h))
+end
+
+dim(a::IsingBimod) = dim(convert(IsingAnyon, a))

test/runtests.jl

@@ -75,6 +75,101 @@ end
     @test size(Rsymbol(a, b, c)) == (0, 0)
 end
 
+@testset "IsingBimod sector" begin
+    I = IsingBimod
+    Istr = TensorKitSectors.type_repr(I)
+    @testset "Basic type properties" begin
+        @test eval(Meta.parse(sprint(show, I))) == I
+        @test eval(Meta.parse(TensorKitSectors.type_repr(I))) == I
+    end
+
+    M = IsingBimod(CatType(3), 0)
+    Mop = IsingBimod(CatType(2), 0)
+    C0 = IsingBimod(CatType(1), 0)
+    C1 = IsingBimod(CatType(1), 1)
+    D0 = IsingBimod(CatType(4), 0)
+    D1 = IsingBimod(CatType(4), 1)
+    C = rand([C0, C1])
+    D = rand([D0, D1])
+    s = rand((M, Mop, C, D))
+
+    @testset "Basic properties" begin
+        @test @constinferred(one(C1)) == @constinferred(leftone(C1)) ==
+              @constinferred(rightone(C1))
+        @test one(D1) == leftone(D1) == rightone(D1)
+        @test one(C1) == leftone(M) == rightone(Mop)
+        @test one(D1) == rightone(M) == leftone(Mop)
+
+        @test eval(Meta.parse(sprint(show, s))) == s
+        @test @constinferred(hash(s)) == hash(deepcopy(s))
+        @constinferred dual(s)
+        @test dual(dual(s)) == s
+        @constinferred dim(s)
+        @constinferred frobeniusschur(s)
+        @constinferred convert(IsingAnyon, s)
+
+        @constinferred Bsymbol(C, C, C)
+        @constinferred Fsymbol(D, D, D, D, D, D)
+    end
+
+    @testset "Sector $Istr: Value iterator" begin
+        @test eltype(values(I)) == I
+        @test_throws ArgumentError one(I)
+        sprev = C0 # first in SectorValues
+        for (i, s) in enumerate(values(I))
+            @test !isless(s, sprev) # confirm compatibility with sort order
+            @test s == @constinferred (values(I)[i])
+            @test findindex(values(I), s) == i
+            sprev = s
+            i >= 10 && break
+        end
+        @test C0 == first(values(I))
+        @test (@constinferred findindex(values(I), C0)) == 1
+        for s in smallset(I)
+            @test (@constinferred values(I)[findindex(values(I), s)]) == s
+        end
+    end
+
+    @testset "Fusion rules and F-symbols" begin
+        argerr = ArgumentError("invalid fusion channel")
+        # forbidden fusions
+        for obs in [(C, D), (D, C), (M, M), (Mop, Mop), (D, M), (M, C), (Mop, D), (C, Mop)]
+            @test isempty(⊗(obs...))
+            @test_throws argerr Nsymbol(obs..., s)
+        end
+
+        # allowed fusions
+        for obs in [(C, C), (D, D), (M, Mop), (Mop, M), (C, M), (Mop, C), (M, D), (D, Mop)]
+            @test !isempty(⊗(obs...))
+        end
+
+        @test Nsymbol(C, C, one(C)) == Nsymbol(D, D, one(D)) == 1
+        @test Nsymbol(C, M, M) == Nsymbol(Mop, C, Mop) == 1
+        @test Nsymbol(M, D, M) == Nsymbol(D, Mop, Mop) == 1
+
+        @test_throws argerr Nsymbol(M, Mop, D)
+        @test_throws argerr Nsymbol(Mop, M, C)
+        @test Nsymbol(M, Mop, C) == Nsymbol(Mop, M, D) == 1
+
+        # non-trivial F-symbol checks
+        @test Fsymbol(Mop, M, D1, D0, D0, M) == 0 # ℳᵒᵖ x ℳ x 𝒟 → ℳ allowed, but 𝒟-labels mismatch
+        @test Fsymbol(Mop, M, D1, D0, D1, M) == 1
+        @test Fsymbol(M, Mop, C1, C0, C0, Mop) == 0
+        @test Fsymbol(M, Mop, C1, C0, C1, Mop) == 1
+
+        @test Fsymbol(C, M, D, M, M, M) == (C.label * D.label == 0 ? 1 : -1) # 𝒞 x ℳ x 𝒟 → ℳ allowed
+        @test_throws argerr Fsymbol(M, Mop, M, Mop, C, D) == 0 # IsingAnyon conversion would give non-zero
+        @test_throws argerr Fsymbol(Mop, M, Mop, M, D, C) == 0
+
+        @test Fsymbol(M, Mop, M, M, C, D) ==
+              (C.label * D.label == 0 ? inv(sqrt(2)) : -inv(sqrt(2))) # ℳ x ℳᵒᵖ x ℳ → ℳ allowed
+        @test Fsymbol(Mop, M, Mop, Mop, D, C) ==
+              (C.label * D.label == 0 ? inv(sqrt(2)) : -inv(sqrt(2))) # ℳᵒᵖ x ℳ x ℳᵒᵖ → ℳᵒᵖ allowed
+
+        @test_throws argerr Fsymbol(M, Mop, M, Mop, C, D)
+    end
+end
+
 @testset "Aqua" begin
     using Aqua: Aqua
     Aqua.test_all(TensorKitSectors)