Generalising functions to support GenericUnit

Project.toml

@@ -43,7 +43,7 @@ Random = "1"
 SafeTestsets = "0.1"
 ScopedValues = "1.3.0"
 Strided = "2"
-TensorKitSectors = "=0.3.0, 0.3.2"
+TensorKitSectors = "0.3.3"
 TensorOperations = "5.1"
 Test = "1"
 TestExtras = "0.2,0.3"

src/TensorKit.jl

@@ -11,11 +11,11 @@ module TensorKit
 export Sector, AbstractIrrep, Irrep
 export FusionStyle, UniqueFusion, MultipleFusion, MultiplicityFreeFusion, SimpleFusion, GenericFusion
 export UnitStyle, SimpleUnit, GenericUnit
-export BraidingStyle, SymmetricBraiding, Bosonic, Fermionic, Anyonic, NoBraiding
+export BraidingStyle, SymmetricBraiding, Bosonic, Fermionic, Anyonic, NoBraiding, HasBraiding
 export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep
 export ProductSector
 export FermionParity, FermionNumber, FermionSpin
-export FibonacciAnyon, IsingAnyon
+export FibonacciAnyon, IsingAnyon, IsingBimodule
 
 # Export common vector space, fusion tree and tensor types
 export VectorSpace, Field, ElementarySpace # abstract vector spaces
@@ -34,7 +34,10 @@ export SpaceMismatch, SectorMismatch, IndexError # error types
 # Export general vector space methods
 export space, field, dual, dim, reduceddim, dims, fuse, flip, isdual
 export unitspace, zerospace, oplus, ominus
+export leftunitspace, rightunitspace, isunitspace
 export insertleftunit, insertrightunit, removeunit
+
+# partial order for vector spaces
 export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic
 
 # Reexport methods for sectors and properties thereof

src/auxiliary/auxiliary.jl

@@ -80,3 +80,9 @@ function _copyto!(A::StridedView{<:Any, 1}, B::StridedView{<:Any, 2})
 
     return A
 end
+
+@static if VERSION < v"1.11" # TODO: remove once support for v1.10 is dropped
+    _allequal(f, xs) = allequal(Base.Generator(f, xs))
+else
+    _allequal(f, xs) = allequal(f, xs)
+end

src/fusiontrees/fusiontrees.jl

@@ -92,7 +92,7 @@ function FusionTree{I}(
         uncoupled::NTuple{N}, coupled = unit(I), isdual = ntuple(Returns(false), N)
     ) where {I <: Sector, N}
     FusionStyle(I) isa UniqueFusion ||
-        error("fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`")
+        throw(ArgumentError("fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`"))
     return FusionTree{I}(
         map(s -> convert(I, s), uncoupled), convert(I, coupled), isdual,
         _abelianinner(map(s -> convert(I, s), (uncoupled..., dual(coupled))))

src/fusiontrees/iterator.jl

@@ -17,6 +17,9 @@ function fusiontrees(uncoupled::Tuple{Vararg{I}}, coupled::I) where {I <: Sector
     return fusiontrees(uncoupled, coupled, isdual)
 end
 function fusiontrees(uncoupled::Tuple{I, Vararg{I}}) where {I <: Sector}
+    UnitStyle(I) isa GenericUnit && throw(
+        ArgumentError("Must specify coupled sector when UnitStyle is GenericUnit.")
+    )
     coupled = unit(I)
     isdual = ntuple(n -> false, length(uncoupled))
     return fusiontrees(uncoupled, coupled, isdual)

src/spaces/gradedspace.jl

@@ -90,10 +90,8 @@ field(::Type{<:GradedSpace}) = ℂ
 InnerProductStyle(::Type{<:GradedSpace}) = EuclideanInnerProduct()
 
 function dim(V::GradedSpace)
-    return reduce(
-        +, dim(V, c) * dim(c) for c in sectors(V);
-        init = zero(dim(unit(sectortype(V))))
-    )
+    init = 0 * dim(first(allunits(sectortype(V))))
+    return sum(c -> dim(c) * dim(V, c), sectors(V); init = init)
 end
 function dim(V::GradedSpace{I, <:AbstractDict}, c::I) where {I <: Sector}
     return get(V.dims, isdual(V) ? dual(c) : c, 0)
@@ -123,8 +121,11 @@ function flip(V::GradedSpace{I}) where {I <: Sector}
     end
 end
 
-unitspace(S::Type{<:GradedSpace{I}}) where {I <: Sector} = S(unit(I) => 1)
-zerospace(S::Type{<:GradedSpace{I}}) where {I <: Sector} = S(unit(I) => 0)
+function unitspace(S::Type{<:GradedSpace{I}}) where {I <: Sector}
+    return S(unit => 1 for unit in allunits(I))
+end
+zerospace(S::Type{<:GradedSpace}) = S()
+
 # TODO: the following methods can probably be implemented more efficiently for
 # `FiniteGradedSpace`, but we don't expect them to be used often in hot loops, so
 # these generic definitions (which are still quite efficient) are good for now.

src/spaces/homspace.jl

@@ -93,11 +93,16 @@ function blocksectors(W::HomSpace)
     N₁ = length(codom)
     N₂ = length(dom)
     I = sectortype(W)
-    # TODO: is sort! still necessary now that blocksectors of ProductSpace is sorted?
-    if N₂ <= N₁
-        return sort!(filter!(c -> hasblock(codom, c), blocksectors(dom)))
+    if N₁ == N₂ == 0
+        return allunits(I)
+    elseif N₁ == 0
+        return filter!(isunit, collect(blocksectors(dom))) # module space cannot end in empty space
+    elseif N₂ == 0
+        return filter!(isunit, collect(blocksectors(codom)))
+    elseif N₂ <= N₁
+        return filter!(c -> hasblock(codom, c), collect(blocksectors(dom)))
     else
-        return sort!(filter!(c -> hasblock(dom, c), blocksectors(codom)))
+        return filter!(c -> hasblock(dom, c), collect(blocksectors(codom)))
     end
 end
 
@@ -227,7 +232,7 @@ function TensorOperations.tensorcontract(
 end
 
 """
-    insertleftunit(W::HomSpace, i=numind(W) + 1; conj=false, dual=false)
+    insertleftunit(W::HomSpace, i = numind(W) + 1; conj = false, dual = false)
 
 Insert a trivial vector space, isomorphic to the underlying field, at position `i`,
 which can be specified as an `Int` or as `Val(i)` for improved type stability.
@@ -237,7 +242,8 @@ See also [`insertrightunit`](@ref insertrightunit(::HomSpace, ::Val{i}) where {i
 [`removeunit`](@ref removeunit(::HomSpace, ::Val{i}) where {i}).
 """
 function insertleftunit(
-        W::HomSpace, ::Val{i} = Val(numind(W) + 1); conj::Bool = false, dual::Bool = false
+        W::HomSpace, ::Val{i} = Val(numind(W) + 1);
+        conj::Bool = false, dual::Bool = false
     ) where {i}
     if i ≤ numout(W)
         return insertleftunit(codomain(W), Val(i); conj, dual) ← domain(W)
@@ -247,7 +253,7 @@ function insertleftunit(
 end
 
 """
-    insertrightunit(W::HomSpace, i=numind(W); conj=false, dual=false)
+    insertrightunit(W::HomSpace, i = numind(W); conj = false, dual = false)
 
 Insert a trivial vector space, isomorphic to the underlying field, after position `i`,
 which can be specified as an `Int` or as `Val(i)` for improved type stability.
@@ -257,7 +263,8 @@ See also [`insertleftunit`](@ref insertleftunit(::HomSpace, ::Val{i}) where {i})
 [`removeunit`](@ref removeunit(::HomSpace, ::Val{i}) where {i}).
 """
 function insertrightunit(
-        W::HomSpace, ::Val{i} = Val(numind(W)); conj::Bool = false, dual::Bool = false
+        W::HomSpace, ::Val{i} = Val(numind(W));
+        conj::Bool = false, dual::Bool = false
     ) where {i}
     if i ≤ numout(W)
         return insertrightunit(codomain(W), Val(i); conj, dual) ← domain(W)

src/spaces/productspace.jl

@@ -1,6 +1,6 @@
 """
-    struct ProductSpace{S<:ElementarySpace, N} <: CompositeSpace{S}
-    ProductSpace(spaces::NTuple{N, S}) where {S<:ElementarySpace, N}
+    struct ProductSpace{S <: ElementarySpace, N} <: CompositeSpace{S}
+    ProductSpace(spaces::NTuple{N, S}) where {S <: ElementarySpace, N}
 
 A `ProductSpace` is a tensor product space of `N` vector spaces of type `S <: ElementarySpace`.
 Only tensor products between [`ElementarySpace`](@ref) objects of the same type are allowed.
@@ -87,7 +87,7 @@ end
 
 # more specific methods
 """
-    sectors(P::ProductSpace{S, N}) where {S<:ElementarySpace}
+    sectors(P::ProductSpace{S, N}) where {S <: ElementarySpace}
 
 Return an iterator over all possible combinations of sectors (represented as an
 `NTuple{N, sectortype(S)}`) that can appear within the tensor product space `P`.
@@ -151,7 +151,7 @@ function blocksectors(P::ProductSpace{S, N}) where {S, N}
     end
     bs = Vector{I}()
     if N == 0
-        push!(bs, unit(I))
+        append!(bs, allunits(I))
     elseif N == 1
         for s in sectors(P)
             push!(bs, first(s))
@@ -196,7 +196,7 @@ hasblock(P::ProductSpace, c::Sector) = !isempty(fusiontrees(P, c))
     blockdim(P::ProductSpace, c::Sector)
 
 Return the total dimension of a coupled sector `c` in the product space, by summing over
-all `dim(P, s)` for all tuples of sectors `s::NTuple{N, <:Sector}` that can fuse to  `c`,
+all `dim(P, s)` for all tuples of sectors `s::NTuple{N, Sector}` that can fuse to  `c`,
 counted with the correct multiplicity (i.e. number of ways in which `s` can fuse to `c`).
 
 See also [`hasblock`](@ref) and [`blocksectors`](@ref).
@@ -228,8 +228,8 @@ end
 
 # unit element with respect to the monoidal structure of taking tensor products
 """
-    one(::S) where {S<:ElementarySpace} -> ProductSpace{S, 0}
-    one(::ProductSpace{S}) where {S<:ElementarySpace} -> ProductSpace{S, 0}
+    one(::S) where {S <: ElementarySpace} -> ProductSpace{S, 0}
+    one(::ProductSpace{S}) where {S <: ElementarySpace} -> ProductSpace{S, 0}
 
 Return a tensor product of zero spaces of type `S`, i.e. this is the unit object under the
 tensor product operation, such that `V ⊗ one(V) == V`.
@@ -248,7 +248,7 @@ fuse(P::ProductSpace{S, 0}) where {S <: ElementarySpace} = unitspace(S)
 fuse(P::ProductSpace{S}) where {S <: ElementarySpace} = fuse(P.spaces...)
 
 """
-    insertleftunit(P::ProductSpace, i::Int=length(P) + 1; conj=false, dual=false)
+    insertleftunit(P::ProductSpace, i::Int = length(P) + 1; conj = false, dual = false)
 
 Insert a trivial vector space, isomorphic to the underlying field, at position `i`,
 which can be specified as an `Int` or as `Val(i)` for improved type stability.
@@ -260,7 +260,18 @@ function insertleftunit(
         P::ProductSpace, ::Val{i} = Val(length(P) + 1);
         conj::Bool = false, dual::Bool = false
     ) where {i}
-    u = unitspace(spacetype(P))
+    N = length(P)
+    I = sectortype(P)
+    if UnitStyle(I) isa SimpleUnit
+        u = unitspace(spacetype(P))
+    else
+        N > 0 || throw(ArgumentError("cannot insert a sensible unit space in the empty product space"))
+        if i == N + 1
+            u = rightunitspace(P[N])
+        else
+            u = leftunitspace(P[i])
+        end
+    end
     if dual
         u = TensorKit.dual(u)
     end
@@ -271,7 +282,7 @@ function insertleftunit(
 end
 
 """
-    insertrightunit(P::ProductSpace, i=lenght(P); conj=false, dual=false)
+    insertrightunit(P::ProductSpace, i = length(P); conj = false, dual = false)
 
 Insert a trivial vector space, isomorphic to the underlying field, after position `i`,
 which can be specified as an `Int` or as `Val(i)` for improved type stability.
@@ -283,7 +294,14 @@ function insertrightunit(
         P::ProductSpace, ::Val{i} = Val(length(P));
         conj::Bool = false, dual::Bool = false
     ) where {i}
-    u = unitspace(spacetype(P))
+    N = length(P)
+    I = sectortype(P)
+    if UnitStyle(I) isa SimpleUnit
+        u = unitspace(spacetype(P))
+    else
+        N > 0 || throw(ArgumentError("cannot insert a sensible unit space in the empty product space"))
+        u = rightunitspace(P[i])
+    end
     if dual
         u = TensorKit.dual(u)
     end
@@ -305,7 +323,8 @@ and [`insertrightunit`](@ref insertrightunit(::ProductSpace, ::Val{i}) where {i}
 """
 function removeunit(P::ProductSpace, ::Val{i}) where {i}
     1 ≤ i ≤ length(P) || _boundserror(P, i)
-    isisomorphic(P[i], unitspace(P[i])) || _nontrivialspaceerror(P, i)
+    I = sectortype(P)
+    isunitspace(P[i]) || _nontrivialspaceerror(P, i)
     return ProductSpace{spacetype(P)}(TupleTools.deleteat(P.spaces, i))
 end
 

src/spaces/vectorspaces.jl

@@ -41,8 +41,8 @@ generally, objects in linear monoidal categories.
 abstract type VectorSpace end
 
 """
-    field(a) -> Type{𝔽<:Field}
-    field(::Type{T}) -> Type{𝔽<:Field}
+    field(a) -> Type{𝔽 <: Field}
+    field(::Type{T}) -> Type{𝔽 <: Field}
 
 Return the type of field over which object `a` (e.g. a vector space or a tensor) is defined.
 This also works in type domain.
@@ -130,10 +130,13 @@ Always returns `false` for spaces where `V == conj(V)`, i.e. vector spaces over
 """ isconj(::ElementarySpace)
 
 """
-    unitspace(V::S) where {S<:ElementarySpace} -> S
+    unitspace(V::S) where {S <: ElementarySpace} -> S
 
 Return the corresponding vector space of type `S` that represents the trivial
 one-dimensional space, i.e. the space that is isomorphic to the corresponding field.
+For vector spaces where `I = sectortype(S)` has a semi-simple unit structure
+(`UnitStyle(I) == GenericUnit()`), this returns a multi-dimensional space corresponding to all unit sectors:
+`dim(unitspace(V), s) == 1` for all `s in allunits(I)`. 
 
 !!! note
     `unitspace(V)`is different from `one(V)`. The latter returns the empty product space
@@ -144,7 +147,7 @@ Base.oneunit(V::ElementarySpace) = unitspace(V)
 Base.oneunit(::Type{V}) where {V <: ElementarySpace} = unitspace(V)
 
 """
-    zerospace(V::S) where {S<:ElementarySpace} -> S
+    zerospace(V::S) where {S <: ElementarySpace} -> S
 
 Return the corresponding vector space of type `S` that represents the zero-dimensional or empty space.
 This is the zero element of the direct sum of vector spaces.
@@ -154,6 +157,68 @@ zerospace(V::ElementarySpace) = zerospace(typeof(V))
 Base.zero(V::ElementarySpace) = zerospace(V)
 Base.zero(::Type{V}) where {V <: ElementarySpace} = zerospace(V)
 
+"""
+    leftunitspace(V::S) where {S <: ElementarySpace} -> S
+
+Return the corresponding vector space of type `S` that represents the trivial
+one-dimensional space, i.e. the space that is isomorphic to the corresponding field. For vector spaces 
+of type `GradedSpace{I}`, this one-dimensional space contains the unique left unit of the objects in `Sector` `I` present
+in the vector space.
+"""
+function leftunitspace(V::ElementarySpace)
+    I = sectortype(V)
+    if UnitStyle(I) isa SimpleUnit
+        return unitspace(typeof(V))
+    else
+        !isempty(sectors(V)) || throw(ArgumentError("Cannot determine the left unit of an empty space"))
+        _allequal(leftunit, sectors(V)) ||
+            throw(ArgumentError("sectors of $V do not have the same left unit"))
+
+        sector = leftunit(first(sectors(V)))
+        return spacetype(V)(sector => 1)
+    end
+end
+
+"""
+    rightunitspace(V::S) where {S <: ElementarySpace} -> S
+
+Return the corresponding vector space of type `ElementarySpace` that represents the trivial
+one-dimensional space, i.e. the space that is isomorphic to the corresponding field. For vector spaces 
+of type `GradedSpace{I}`, this corresponds to the right unit of the objects in `Sector` `I` present
+in the vector space.
+"""
+function rightunitspace(V::ElementarySpace)
+    I = sectortype(V)
+    if UnitStyle(I) isa SimpleUnit
+        return unitspace(typeof(V))
+    else
+        !isempty(sectors(V)) || throw(ArgumentError("Cannot determine the right unit of an empty space"))
+        _allequal(rightunit, sectors(V)) ||
+            throw(ArgumentError("sectors of $V do not have the same right unit"))
+
+        sector = rightunit(first(sectors(V)))
+        return spacetype(V)(sector => 1)
+    end
+end
+
+"""
+    isunitspace(V::S) where {S <: ElementarySpace} -> Bool
+
+Return whether the elementary space `V` is a unit space, i.e. is isomorphic to the
+trivial one-dimensional space. For vector spaces of type `GradedSpace{I}` where `Sector` `I` has a
+semi-simple unit structure, this returns `true` if `V` is isomorphic to either the left, right or
+semi-simple unit space.
+"""
+function isunitspace(V::ElementarySpace)
+    I = sectortype(V)
+    return if isa(UnitStyle(I), SimpleUnit)
+        isisomorphic(V, unitspace(V))
+    else
+        (dim(V) == 0 || !all(isunit, sectors(V))) && return false
+        return true
+    end
+end
+
 """
     ⊕(V₁::S, V₂::S, V₃::S...) where {S<:ElementarySpace} -> S
     oplus(V₁::S, V₂::S, V₃::S...) where {S<:ElementarySpace} -> S

src/tensors/diagonal.jl

@@ -306,7 +306,7 @@ function LinearAlgebra.pinv(d::DiagonalTensorMap; kwargs...)
     if iszero(atol)
         rtol = get(kwargs, :rtol, zero(real(T)))
     else
-        rtol = sqrt(eps(real(float(oneunit(T))))) * length(d.data)
+        rtol = sqrt(eps(real(float(one(T))))) * length(d.data)
     end
     pdata = let tol = max(atol, rtol * maximum(abs, d.data))
         map(x -> abs(x) < tol ? zero(x) : pinv(x), d.data)

src/tensors/indexmanipulations.jl

@@ -267,7 +267,7 @@ function has_shared_twist(t, inds)
     if BraidingStyle(I) == NoBraiding()
         for i in inds
             cs = sectors(space(t, i))
-            all(isone, cs) || throw(SectorMismatch(lazy"Cannot twist sectors $cs"))
+            all(isunit, cs) || throw(SectorMismatch(lazy"Cannot twist sectors $cs"))
         end
         return true
     elseif BraidingStyle(I) == Bosonic()

src/tensors/linalg.jl

@@ -296,7 +296,7 @@ function LinearAlgebra.rank(
         t::AbstractTensorMap;
         atol::Real = 0, rtol::Real = atol > 0 ? 0 : _default_rtol(t)
     )
-    r = dim(one(sectortype(t))) * 0
+    r = 0 * dim(first(allunits(sectortype(t))))
     dim(t) == 0 && return r
     S = LinearAlgebra.svdvals(t)
     tol = max(atol, rtol * maximum(first, values(S)))