Finishing touches to a working MultiTensorKit code

Project.toml

@@ -5,15 +5,19 @@ version = "0.1.0"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+BlockTensorKit = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f"
+TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 
 [compat]
 Aqua = "0.8.9"
 Artifacts = "1.10, 1"
 JSON3 = "1.14.1"
 SafeTestsets = "0.1"
-TensorKitSectors = "0.1.2"
+TensorKitSectors = "0.1.4"
 Test = "1.10"
 TestExtras = "0.3"
 julia = "1.10"

src/MultiTensorKit.jl

@@ -3,9 +3,19 @@ module MultiTensorKit
 export BimoduleSector, A4Object
 
 using JSON3
+using DelimitedFiles
 using Artifacts
 using TensorKitSectors
 
+using TupleTools
+import TupleTools: insertafter
+
+using BlockTensorKit
+import BlockTensorKit: SumSpace
+
+using TensorKit
+import TensorKit: hasblock, dim
+
 include("bimodulesector.jl")
 
 end

src/bimodulesector.jl

@@ -2,6 +2,11 @@ struct BimoduleSector{Name} <: Sector
     i::Int
     j::Int
     label::Int
+    function BimoduleSector{:A4}(i::Int, j::Int, label::Int)
+        i <= 12 && j <= 12 || throw(DomainError("object outside the matrix A4"))
+        return label <= _numlabels(BimoduleSector{:A4}, i, j) ? new{:A4}(i, j, label) :
+               throw(DomainError("label outside category A4($i, $j)"))
+    end
 end
 BimoduleSector{Name}(data::NTuple{3,Int}) where {Name} = BimoduleSector{Name}(data...)
 const A4Object = BimoduleSector{:A4}
@@ -16,23 +21,27 @@ function Base.convert(::Type{BimoduleSector{Name}}, d::NTuple{3,Int}) where {Nam
     return BimoduleSector{Name}(d...)
 end
 
-Base.IteratorSize(::Type{SectorValues{<:BimoduleSector}}) = Base.SizeUnknown()
+Base.IteratorSize(::Type{<:SectorValues{<:BimoduleSector}}) = Base.SizeUnknown()
 
 # TODO: generalize?
-# TODO: numlabels?
 function Base.iterate(iter::SectorValues{A4Object}, (I, label)=(1, 1))
     I > 12 * 12 && return nothing
     i, j = CartesianIndices((12, 12))[I].I
-    maxlabel = numlabels(A4Object, i, j)
+    maxlabel = _numlabels(A4Object, i, j)
     return if label > maxlabel
         iterate(iter, (I + 1, 1))
     else
         A4Object(i, j, label), (I, label + 1)
     end
 end
 
+function Base.length(::SectorValues{A4Object})
+    return sum(_numlabels(A4Object, i, j) for i in 1:12, j in 1:12)
+end
+
 TensorKitSectors.FusionStyle(::Type{A4Object}) = GenericFusion()
 TensorKitSectors.BraidingStyle(::Type{A4Object}) = NoBraiding()
+TensorKitSectors.sectorscalartype(::Type{A4Object}) = ComplexF64
 
 function TensorKitSectors.:⊗(a::A4Object, b::A4Object)
     @assert a.j == b.i
@@ -42,8 +51,8 @@ function TensorKitSectors.:⊗(a::A4Object, b::A4Object)
                     if (a_l == a.label && b_l == b.label)]
 end
 
-function _numlabels(::Type{A4Object}, i, j)
-    return Ncache = _get_Ncache(A4Object)
+function _numlabels(::Type{T}, i, j) where {T<:BimoduleSector}
+    return length(_get_dual_cache(T)[2][i, j])
 end
 
 # Data from files
@@ -82,8 +91,7 @@ function TensorKitSectors.Nsymbol(a::I, b::I, c::I) where {I<:A4Object}
     return get(_get_Ncache(I)[i, j, k], (a.label, b.label, c.label), 0)
 end
 
-# TODO: can we define dual for modules?
-const Dualcache = IdDict{Type{<:BimoduleSector},Vector{Tuple{Int,Vector{Int}}}}()
+const Dualcache = IdDict{Type{<:BimoduleSector},Tuple{Vector{Int64},Matrix{Vector{Int64}}}}()
 
 function _get_dual_cache(::Type{T}) where {T<:BimoduleSector}
     global Dualcache
@@ -95,74 +103,122 @@ end
 function extract_dual(::Type{A4Object})
     N = _get_Ncache(A4Object)
     ncats = size(N, 1)
-    return map(1:ncats) do i
+    Is = zeros(Int, ncats)
+
+    map(1:ncats) do i
         Niii = N[i, i, i]
-        nobj = maximum(first, keys(Niii))
-
-        # find identity object:
-        # I x I -> I, a x I -> a, I x a -> a
-        I = findfirst(1:nobj) do u
-            get(Niii, (u, u, u), 0) == 1 || return false
-            for j in 1:nobj
-                get(Niii, (j, u, j), 0) == 1 || return false
-                get(Niii, (u, j, j), 0) == 1 || return false
+        nobji = maximum(first, keys(N[i, i, i]))
+        # want to return a leftone and rightone for each entry in multifusion cat
+        # leftone/rightone needs to at least be the unit object within a fusion cat
+        Is[i] = findfirst(1:nobji) do a
+            get(Niii, (a, a, a), 0) == 1 || return false # I x I -> I
+            for othera in 1:nobji
+                get(Niii, (othera, a, othera), 0) == 1 || return false # a x I -> a
+                get(Niii, (a, othera, othera), 0) == 1 || return false # I x a -> a
+            end
+
+            # check leftone
+            map(1:ncats) do j
+                nobjj = maximum(first, keys(N[j, j, j]))
+                for b in 1:nobjj
+                    get(N[i, j, j], (a, b, b), 0) == 1 || return false # I = leftone(b)
+                end
+            end
+
+            # check rightone
+            map(1:ncats) do k
+                nobjk = maximum(first, keys(N[k, k, k]))
+                for c in 1:nobjk
+                    get(N[k, i, k], (c, a, c), 0) == 1 || return false # I = rightone(c)
+                end
             end
             return true
         end
+    end
 
-        # find duals
-        # a x abar -> I
-        duals = map(1:nobj) do j
-            return findfirst(1:nobj) do k
-                return get(Niii, (j, k, I), 0) == 1
+    allduals = Matrix{Vector{Int}}(undef, ncats, ncats) # ncats square matrix of vectors
+    for i in 1:ncats
+        nobji = maximum(first, keys(N[i, i, i]))
+        for j in 1:ncats
+            allduals[i, j] = Int[]
+
+            nobjj = maximum(first, keys(N[j, j, j]))
+            # the nested vectors contain the duals of the objects in 𝒞_ij, which are in C_ji 
+            Niji = N[i, j, i] # 𝒞_ij x 𝒞_ji -> C_ii
+            Njij = N[j, i, j] # 𝒞_ji x 𝒞_ij -> C_jj
+            for i_ob in 1:nobji, j_ob in 1:nobjj
+                get(Niji, (i_ob, j_ob, Is[i]), 0) == 1 || continue # leftone(c_ij) ∈ c_ij x c_ji
+                get(Njij, (j_ob, i_ob, Is[j]), 0) == 1 || continue # rightone(c_ij) ∈ c_ji x c_ij
+                push!(allduals[i, j], j_ob)
             end
         end
-
-        return I, duals
     end
+    return Is, allduals
 end
 
 function Base.one(a::BimoduleSector)
-    a.i == a.j || error("don't know how to define one for modules")
-    return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[a.i][1])
+    a.i == a.j || error("unit object for module categories is ill-defined")
+    return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[1][a.i])
+end
+
+function TensorKitSectors.leftone(a::BimoduleSector)
+    return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[1][a.i])
+end
+
+function TensorKitSectors.rightone(a::BimoduleSector)
+    return A4Object(a.j, a.j, _get_dual_cache(typeof(a))[1][a.j])
 end
 
 function Base.conj(a::BimoduleSector)
-    a.i == a.j || error("don't know how to define dual for modules")
-    return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[a.i][2][a.label])
+    return A4Object(a.j, a.i, _get_dual_cache(typeof(a))[2][a.i, a.j][a.label])
 end
 
 function extract_Fsymbol(::Type{A4Object})
-    return mapreduce((x, y) -> cat(x, y; dims=1), 1:12) do i
-        filename = joinpath(artifact_path, "A4", "Fsymbol_$i.json")
+    result = Dict{NTuple{4,Int},Dict{NTuple{6,Int},Array{ComplexF64,4}}}()
+    for i in 1:12
+        filename = joinpath(artifact_path, "A4", "Fsymbol_$i.txt")
         @assert isfile(filename) "cannot find $filename"
-        json_string = read(filename, String)
-        Farray_part = copy(JSON3.read(json_string))
-        return map(enumerate(Farray_part[Symbol(i)])) do (I, x)
-            j, k, l = Tuple(CartesianIndices((12, 12, 12))[I])
-            y = Dict{NTuple{6,Int},Array{ComplexF64,4}}()
-            for (key, v) in x
-                a, b, c, d, e, f = parse.(Int, split(string(key)[2:(end - 1)], ", "))
+        Farray_part = readdlm(filename)
+        for ((i, j, k, l), colordict) in convert_Fs(Farray_part)
+            result[(i, j, k, l)] = Dict{NTuple{6,Int},Array{ComplexF64,4}}()
+            for ((a, b, c, d, e, f), Fvals) in colordict
                 a_ob, b_ob, c_ob, d_ob, e_ob, f_ob = A4Object.(((i, j, a), (j, k, b),
                                                                 (k, l, c), (i, l, d),
                                                                 (i, k, e), (j, l, f)))
-                result = Array{ComplexF64,4}(undef,
-                                             (Nsymbol(a_ob, b_ob, e_ob),
-                                              Nsymbol(e_ob, c_ob, d_ob),
-                                              Nsymbol(b_ob, c_ob, f_ob),
-                                              Nsymbol(a_ob, f_ob, d_ob)))
-                map!(result, reshape(v, size(result))) do cmplxdict
-                    return complex(cmplxdict[:re], cmplxdict[:im])
+                result[(i, j, k, l)][(a, b, c, d, e, f)] = zeros(ComplexF64,
+                                                                 Nsymbol(a_ob, b_ob, e_ob),
+                                                                 Nsymbol(e_ob, c_ob, d_ob),
+                                                                 Nsymbol(b_ob, c_ob, f_ob),
+                                                                 Nsymbol(a_ob, f_ob, d_ob))
+                for (I, v) in Fvals
+                    result[(i, j, k, l)][(a, b, c, d, e, f)][I] = v
                 end
-
-                y[(a, b, c, d, e, f)] = result
             end
         end
     end
+    return result
+end
+
+function convert_Fs(Farray_part::Matrix{Float64}) # Farray_part is a matrix with 16 columns
+    data_dict = Dict{NTuple{4,Int},
+                     Dict{NTuple{6,Int},Vector{Pair{CartesianIndex{4},ComplexF64}}}}()
+    # want to make a Dict with keys (i,j,k,l) and vals 
+    # a Dict with keys (a,b,c,d,e,f) and vals 
+    # a pair of (mu, nu, rho, sigma) and the F value
+    for row in eachrow(Farray_part)
+        i, j, k, l, a, b, c, d, e, f, mu, nu, rho, sigma = Int.(@view(row[1:14]))
+        v = complex(row[15], row[16])
+        colordict = get!(data_dict, (i, j, k, l),
+                         Dict{NTuple{6,Int},Vector{Pair{CartesianIndex{4},ComplexF64}}}())
+        Fdict = get!(colordict, (a, b, c, d, e, f),
+                     Vector{Pair{CartesianIndex{4},ComplexF64}}())
+        push!(Fdict, CartesianIndex(mu, nu, rho, sigma) => v)
+    end
+    return data_dict
 end
 
 const Fcache = IdDict{Type{<:BimoduleSector},
-                      Array{Dict{NTuple{6,Int},Array{ComplexF64,4}},4}}()
+                      Dict{NTuple{4,Int64},Dict{NTuple{6,Int64},Array{ComplexF64,4}}}}()
 
 function _get_Fcache(::Type{T}) where {T<:BimoduleSector}
     global Fcache
@@ -173,16 +229,95 @@ end
 
 function TensorKitSectors.Fsymbol(a::I, b::I, c::I, d::I, e::I,
                                   f::I) where {I<:A4Object}
-    # TODO: should this error or return 0?
-    (a.j == b.i && b.j == c.i && a.i == d.i && c.j == d.j &&
-     a.i == e.i && b.j == e.j && f.i == a.j && f.j == c.j) ||
-        throw(ArgumentError("invalid fusion channel"))
+    # required to keep track of multiplicities where F-move is partially unallowed
+    # also deals with invalid fusion channels
+    Nabe = Nsymbol(a, b, e)
+    Necd = Nsymbol(e, c, d)
+    Nbcf = Nsymbol(b, c, f)
+    Nafd = Nsymbol(a, f, d)
+
+    Nabe > 0 && Necd > 0 && Nbcf > 0 && Nafd > 0 ||
+        return zeros(sectorscalartype(I), Nabe, Necd, Nbcf, Nafd)
 
     i, j, k, l = a.i, a.j, b.j, c.j
-    return get(_get_Fcache(I)[i, j, k, l],
-               (a.label, b.label, c.label, d.label, e.label, f.label)) do
-        return zeros(sectorscalartype(A4Object),
-                     (Nsymbol(a, b, e), Nsymbol(e, c, d), Nsymbol(b, c, f),
-                      Nsymbol(a, f, d)))
+    colordict = _get_Fcache(I)[i, j, k, l]
+    return colordict[(a.label, b.label, c.label, d.label, e.label, f.label)]
+end
+
+# interface with TensorKit where necessary
+#-----------------------------------------
+
+function TensorKit.blocksectors(W::TensorMapSpace{S,N₁,N₂}) where
+         {S<:Union{Vect[A4Object],
+                   SumSpace{Vect[A4Object]}},N₁,N₂}
+    codom = codomain(W)
+    dom = domain(W)
+    if N₁ == 0 && N₂ == 0 # 0x0-dimensional TensorMap is just a scalar, return all units
+        # this is a problem in full contractions where the coloring outside is 𝒞
+        return NTuple{12,A4Object}(one(A4Object(i, i, 1)) for i in 1:12) # have to return all units b/c no info on W in this case
+    elseif N₁ == 0
+        @assert N₂ != 0 "one of Type A4Object doesn't exist"
+        return filter!(isone, collect(blocksectors(dom)))
+    elseif N₂ == 0
+        @assert N₁ != 0 "one of Type A4Object 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
+
+# TODO: definition for zero of GradedSpace?
+
+function dim(V::GradedSpace{<:BimoduleSector})
+    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
+
+# limited oneunit 
+function Base.oneunit(S::GradedSpace{<:BimoduleSector})
+    allequal(a.i for a in sectors(S)) && allequal(a.j for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S are not all equal"))
+    first(sectors(S)).i == first(sectors(S)).j ||
+        throw(ArgumentError("sectors of $S are non-diagonal"))
+    sector = one(first(sectors(S)))
+    return ℂ[A4Object](sector => 1)
+end
+
+function Base.oneunit(S::SumSpace{<:GradedSpace{<:BimoduleSector}})
+    allequal(a.i for a in sectors(S)) && allequal(a.j for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S are not all equal"))
+    first(sectors(S)).i == first(sectors(S)).j ||
+        throw(ArgumentError("sectors of $S are non-diagonal"))
+    sector = one(first(sectors(S)))
+    return SumSpace(ℂ[A4Object](sector => 1))
+end
+
+# maybe from the homspace
+function TensorKit.insertrightunit(P::ProductSpace{V,N}, ::Val{i}=Val(length(P));
+                                   conj::Bool=false,
+                                   dual::Bool=false) where {i,V<:GradedSpace{<:I},N} where {I<:BimoduleSector}
+    #possible change to rightone of correct space for N = 0
+    u = N > 0 ? oneunit(P[1]) : 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
+
+function TensorKit.insertleftunit(P::ProductSpace{V,N}, ::Val{i}=Val(length(P) + 1);
+                                  conj::Bool=false,
+                                  dual::Bool=false) where {i,V<:GradedSpace{<:I},N} where {I<:BimoduleSector}
+    u = N > 0 ? oneunit(P[1]) : 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