Compatibility with MultiTensorKit (#247)

* This PR extends the `one` functionality to a `leftone` and `rightone` in order to support multifusion categories, and updates the fusion tree manipulations accordingly.
* This PR also includes a number of other important fixes, in the docs, in `eltype(::BlockIterator)`, in the scalar type determination in tensor operations, etc.

---------

Co-authored-by: Lukas Devos <[email protected]>

Author: borisdevos

docs/src/man/categories.md

@@ -229,7 +229,7 @@ morphism from ``I`` to ``V``. To map morphisms from ``\mathrm{Hom}(W,V)`` to ele
 ``V ⊗ W^*``, i.e. morphisms in ``\mathrm{Hom}(I, V ⊗ W^*)``, we use another morphism
 ``\mathrm{Hom}(I, W ⊗ W^*)`` which can be considered as the inverse of the evaluation map.
 
-Hence, duality in a monoidal category is defined via an *exact paring*, i.e. two families
+Hence, duality in a monoidal category is defined via an *exact pairing*, i.e. two families
 of non-degenerate morphisms, the evaluation (or co-unit) ``ϵ_V: {}^{∨}V ⊗ V → I`` and the
 coevaluation (or unit) ``η_V: I → V ⊗ {}^{∨}V`` which satisfy the "snake rules":
 

docs/src/man/sectors.md

@@ -773,7 +773,7 @@ groups.
 
 Other methods for `ElementarySpace`, such as [`dual`](@ref), [`fuse`](@ref) and
 [`flip`](@ref) also work. In fact, `GradedSpace` is the reason `flip` exists, cause
-in this case it is different then `dual`. The existence of flip originates from the
+in this case it is different than `dual`. The existence of flip originates from the
 non-trivial isomorphism between ``R_{\overline{a}}`` and ``R_{a}^*``, i.e. the
 representation space of the dual ``\overline{a}`` of sector ``a`` and the dual of the
 representation space of sector ``a``. In order for `flip(V)` to be isomorphic to `V`, it is
@@ -898,7 +898,7 @@ for the specific case ``N_1=4`` and ``N_2=3``. We can separate this tree into th
 part ``(b_1⊗b_2)⊗b_3 → c`` and the splitting part ``c→(((a_1⊗a_2)⊗a_3)⊗a_4)``. Given that
 the fusion tree can be considered to be the adjoint of a corresponding splitting tree
 ``c→(b_1⊗b_2)⊗b_3``, we now first consider splitting trees in isolation. A splitting tree
-which goes from one coupled sectors ``c`` to ``N`` uncoupled sectors ``a_1``, ``a_2``, …,
+which goes from one coupled sector ``c`` to ``N`` uncoupled sectors ``a_1``, ``a_2``, …,
 ``a_N`` needs ``N-2`` additional internal sector labels ``e_1``, …, ``e_{N-2}``, and, if
 `FusionStyle(I) isa GenericFusion`, ``N-1`` additional multiplicity labels ``μ_1``,
 …, ``μ_{N-1}``. We henceforth refer to them as vertex labels, as they are associated with
@@ -912,7 +912,7 @@ the orthogonality condition
 which now forces all internal lines ``e_k`` and vertex labels ``μ_l`` to be the same.
 
 There is one subtle remark that we have so far ignored. Within the specific subtypes of
-`Sector`, we do not explicitly distinguish between ``R_a^*`` (simply denoted as ``a`^*``
+`Sector`, we do not explicitly distinguish between ``R_a^*`` (simply denoted as ``a^*``
 and graphically depicted as an upgoing arrow ``a``) and ``R_{\bar{a}}`` (simply denoted as
 ``\bar{a}`` and depicted with a downgoing arrow), i.e. between the dual space of ``R_a`` on
 which the conjugated irrep acts, or the irrep ``\bar{a}`` to which the complex conjugate of

docs/src/man/spaces.md

@@ -297,8 +297,8 @@ corresponding spaces, but in general none of those will be canonical.
 There are also a number of convenience functions to create isomorphic spaces. The function
 `fuse(V1, V2, ...)` or `fuse(V1 ⊗ V2 ⊗ ...)` returns an elementary space that is isomorphic
 to `V1 ⊗ V2 ⊗ ...`. The function `flip(V::ElementarySpace)` returns a space that is
-isomorphic to `V` but has `isdual(flip(V)) == isdual(V')`, i.e. if `V` is a normal space
-than `flip(V)` is a dual space. `flip(V)` is different from `dual(V)` in the case of
+isomorphic to `V` but has `isdual(flip(V)) == isdual(V')`, i.e., if `V` is a normal space,
+then `flip(V)` is a dual space. `flip(V)` is different from `dual(V)` in the case of
 [`GradedSpace`](@ref). It is useful to flip a tensor index from a ket to a bra (or
 vice versa), by contracting that index with a unitary map from `V1` to `flip(V1)`. We refer
 to the reference on [vector space methods](@ref s_spacemethods) for further information.

src/fusiontrees/iterator.jl

@@ -157,7 +157,7 @@ function _fusiontree_iterate(uncoupledsectors::NTuple{N},
         nextout = iterate(outiterN, outstateN)
         nextout === nothing && return nothing
         b, outstateN = nextout
-        vertexiterN = c ⊗ dual(b)
+        vertexiterN = coupled ⊗ dual(b)
         nextline = iterate(vertexiterN)
     end
     a, vertexstateN = nextline

src/fusiontrees/manipulations.jl

@@ -165,8 +165,9 @@ operation is the inverse of `insertat` in the sense that if
         f₂ = FusionTree{I}(f.uncoupled, f.coupled, isdual2, f.innerlines, f.vertices)
         return f₁, f₂
     elseif M === 0
-        f₁ = FusionTree{I}((), one(I), (), ())
-        uncoupled2 = (one(I), f.uncoupled...)
+        u = leftone(f.uncoupled[1])
+        f₁ = FusionTree{I}((), u, (), ())
+        uncoupled2 = (u, f.uncoupled...)
         coupled2 = f.coupled
         isdual2 = (false, f.isdual...)
         innerlines2 = N >= 2 ? (f.uncoupled[1], f.innerlines...) : ()
@@ -230,7 +231,7 @@ function merge(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂},
     return insertat(f, N₁ + 1, f₂)
 end
 function merge(f₁::FusionTree{I,0}, f₂::FusionTree{I,0}, c::I, μ) where {I}
-    isone(c) ||
+    Nsymbol(f₁.coupled, f₂.coupled, c) == μ == 1 ||
         throw(SectorMismatch("cannot fuse sectors $(f₁.coupled) and $(f₂.coupled) to $c"))
     return fusiontreedict(I)(f₁ => Fsymbol(c, c, c, c, c, c)[1, 1, 1, 1])
 end
@@ -289,7 +290,8 @@ function bendright(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {I<
     # map final splitting vertex (a, b)<-c to fusion vertex a<-(c, dual(b))
     @assert N₁ > 0
     c = f₁.coupled
-    a = N₁ == 1 ? one(I) : (N₁ == 2 ? f₁.uncoupled[1] : f₁.innerlines[end])
+    a = N₁ == 1 ? leftone(f₁.uncoupled[1]) :
+        (N₁ == 2 ? f₁.uncoupled[1] : f₁.innerlines[end])
     b = f₁.uncoupled[N₁]
 
     uncoupled1 = TupleTools.front(f₁.uncoupled)
@@ -360,7 +362,7 @@ function foldright(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {I<
         hasmultiplicities = FusionStyle(a) isa GenericFusion
         local newtrees
         if N₁ == 1
-            cset = (one(c1),)
+            cset = (leftone(c1),) # or rightone(a)
         elseif N₁ == 2
             cset = (f₁.uncoupled[2],)
         else
@@ -729,7 +731,8 @@ function elementary_trace(f::FusionTree{I,N}, i) where {I<:Sector,N}
     # if trace is zero, return empty dict
     (b == dual(b′) && f.isdual[i] != f.isdual[j]) || return newtrees
     if i < N
-        inner_extended = (one(I), f.uncoupled[1], f.innerlines..., f.coupled)
+        inner_extended = (leftone(f.uncoupled[1]), f.uncoupled[1], f.innerlines...,
+                          f.coupled)
         a = inner_extended[i]
         d = inner_extended[i + 2]
         a == d || return newtrees
@@ -753,11 +756,11 @@ function elementary_trace(f::FusionTree{I,N}, i) where {I<:Sector,N}
         if i > 1
             c = f.innerlines[i - 1]
             if FusionStyle(I) isa MultiplicityFreeFusion
-                coeff *= Fsymbol(a, b, dual(b), a, c, one(I))
+                coeff *= Fsymbol(a, b, dual(b), a, c, rightone(a))
             else
                 μ = f.vertices[i - 1]
                 ν = f.vertices[i]
-                coeff *= Fsymbol(a, b, dual(b), a, c, one(I))[μ, ν, 1, 1]
+                coeff *= Fsymbol(a, b, dual(b), a, c, rightone(a))[μ, ν, 1, 1]
             end
         end
         if f.isdual[i]
@@ -766,8 +769,9 @@ function elementary_trace(f::FusionTree{I,N}, i) where {I<:Sector,N}
         push!(newtrees, f′ => coeff)
         return newtrees
     else # i == N
+        unit = leftone(b)
         if N == 2
-            f′ = FusionTree{I}((), one(I), (), (), ())
+            f′ = FusionTree{I}((), unit, (), (), ())
             coeff = sqrtdim(b)
             if !(f.isdual[N])
                 coeff *= conj(frobeniusschur(b))
@@ -778,18 +782,17 @@ function elementary_trace(f::FusionTree{I,N}, i) where {I<:Sector,N}
         uncoupled_ = TupleTools.front(f.uncoupled)
         inner_ = TupleTools.front(f.innerlines)
         coupled_ = f.innerlines[end]
-        @assert coupled_ == dual(b)
         isdual_ = TupleTools.front(f.isdual)
         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_, one(I), 1)
-            f_′.innerlines[1] == one(I) || continue
+        for (f_′, coeff) in merge(fs, f_, unit, 1) # coloring gets reversed here, should be the other unit
+            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′, one(I), isdual′, inner′, vertices′)
+            f′ = FusionTree(uncoupled′, unit, isdual′, inner′, vertices′) # and this one?
             coeff *= sqrtdim(b)
             if !(f.isdual[N])
                 coeff *= conj(frobeniusschur(b))

src/spaces/gradedspace.jl

@@ -12,7 +12,7 @@ isomorphism classes of simple objects of a unitary and pivotal (pre-)fusion cate
 
 Here `dims` represents the degeneracy or multiplicity of every sector.
 
-The data structure `D` of `dims` will depend on the result `Base.IteratorElsize(values(I))`;
+The data structure `D` of `dims` will depend on the result `Base.IteratorSize(values(I))`;
 if the result is of type `HasLength` or `HasShape`, `dims` will be stored in a
 `NTuple{N,Int}` with `N = length(values(I))`. This requires that a sector `s::I` can be
 transformed into an index via `s == getindex(values(I), i)` and
@@ -189,16 +189,16 @@ end
 
 function Base.show(io::IO, V::GradedSpace{I}) where {I<:Sector}
     print(io, type_repr(typeof(V)), "(")
-    seperator = ""
+    separator = ""
     comma = ", "
     io2 = IOContext(io, :typeinfo => I)
     for c in sectors(V)
         if isdual(V)
-            print(io2, seperator, dual(c), "=>", dim(V, c))
+            print(io2, separator, dual(c), "=>", dim(V, c))
         else
-            print(io2, seperator, c, "=>", dim(V, c))
+            print(io2, separator, c, "=>", dim(V, c))
         end
-        seperator = comma
+        separator = comma
     end
     print(io, ")")
     V.dual && print(io, "'")

src/tensors/blockiterator.jl

@@ -10,6 +10,5 @@ end
 
 Base.IteratorSize(::BlockIterator) = Base.HasLength()
 Base.IteratorEltype(::BlockIterator) = Base.HasEltype()
-Base.eltype(::Type{<:BlockIterator{T}}) where {T} = blocktype(T)
+Base.eltype(::Type{<:BlockIterator{T}}) where {T} = Pair{sectortype(T),blocktype(T)}
 Base.length(iter::BlockIterator) = length(iter.structure)
-Base.isdone(iter::BlockIterator, state...) = Base.isdone(iter.structure, state...)

src/tensors/tensoroperations.jl

@@ -52,7 +52,8 @@ end
 
 function TO.tensoradd_type(TC, A::AbstractTensorMap, ::Index2Tuple{N₁,N₂},
                            ::Bool) where {N₁,N₂}
-    M = similarstoragetype(A, TC)
+    I = sectortype(A)
+    M = similarstoragetype(A, sectorscalartype(I) <: Real ? TC : complex(TC))
     return tensormaptype(spacetype(A), N₁, N₂, M)
 end
 
@@ -113,10 +114,11 @@ function TO.tensorcontract_type(TC,
                                 A::AbstractTensorMap, ::Index2Tuple, ::Bool,
                                 B::AbstractTensorMap, ::Index2Tuple, ::Bool,
                                 ::Index2Tuple{N₁,N₂}) where {N₁,N₂}
-    M = similarstoragetype(A, TC)
-    M == similarstoragetype(B, TC) ||
-        throw(ArgumentError("incompatible storage types:\n$(M) ≠ $(similarstoragetype(B, TC))"))
     spacetype(A) == spacetype(B) || throw(SpaceMismatch("incompatible space types"))
+    I = sectortype(A)
+    M = similarstoragetype(A, sectorscalartype(I) <: Real ? TC : complex(TC))
+    MB = similarstoragetype(B, sectorscalartype(I) <: Real ? TC : complex(TC))
+    M == MB || throw(ArgumentError("incompatible storage types:\n$(M) ≠ $(MB)"))
     return tensormaptype(spacetype(A), N₁, N₂, M)
 end
 

test/diagonal.jl

@@ -29,7 +29,8 @@ diagspacelist = ((ℂ^4)', ℂ[Z2Irrep](0 => 2, 1 => 3),
             next = @constinferred Nothing iterate(bs, state)
             b2 = @constinferred block(t, first(blocksectors(t)))
             @test b1 == b2
-            @test eltype(bs) === typeof(b1) === TensorKit.blocktype(t)
+            @test eltype(bs) === Pair{typeof(c),typeof(b1)}
+            @test typeof(b1) === TensorKit.blocktype(t)
             # basic linear algebra
             @test isa(@constinferred(norm(t)), real(T))
             @test norm(t)^2 ≈ dot(t, t)

test/tensors.jl

@@ -47,7 +47,9 @@ for V in spacelist
                 next = @constinferred Nothing iterate(bs, state)
                 b2 = @constinferred block(t, first(blocksectors(t)))
                 @test b1 == b2
-                @test eltype(bs) === typeof(b1) === TensorKit.blocktype(t)
+                @test eltype(bs) === Pair{typeof(c),typeof(b1)}
+                @test typeof(b1) === TensorKit.blocktype(t)
+                @test typeof(c) === sectortype(t)
             end
         end
         @timedtestset "Tensor Dict conversion" begin
@@ -107,7 +109,9 @@ for V in spacelist
                 next = @constinferred Nothing iterate(bs, state)
                 b2 = @constinferred block(t', first(blocksectors(t')))
                 @test b1 == b2
-                @test eltype(bs) === typeof(b1) === TensorKit.blocktype(t')
+                @test eltype(bs) === Pair{typeof(c),typeof(b1)}
+                @test typeof(b1) === TensorKit.blocktype(t')
+                @test typeof(c) === sectortype(t)
                 # linear algebra
                 @test isa(@constinferred(norm(t)), real(T))
                 @test norm(t)^2 ≈ dot(t, t)