From b3fe4310a79a67e67c374cc1efa2fca225291f58 Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Wed, 20 Mar 2024 09:08:52 +0100
Subject: [PATCH 01/21] Add .vscode to gitignore

---
 .gitignore | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index ae7f128c..5df92310 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,4 +5,5 @@
 *-old
 __pycache__
 .ipynb*
-Manifest.toml
\ No newline at end of file
+Manifest.toml
+.vscode
\ No newline at end of file

From 440dccee5d44260b0337cf78719bd588838c36db Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Wed, 20 Mar 2024 09:33:19 +0100
Subject: [PATCH 02/21] Add changelog.md

---
 Changelog.md | 15 +++++++++++++++
 1 file changed, 15 insertions(+)
 create mode 100644 Changelog.md

diff --git a/Changelog.md b/Changelog.md
new file mode 100644
index 00000000..295d6a18
--- /dev/null
+++ b/Changelog.md
@@ -0,0 +1,15 @@
+# Planned changes for v1.0
+
+Features that are planned to be implemented before the release of v1.0.0, in no particular order.
+
+- [ ] Separate `Sectors` module
+- [ ] Make `TrivialTensorMap` and `TensorMap` be the same
+- [ ] Simplify `TensorMap` type to hide `rowr` and `colr`
+- [ ] Make `AdjointTensorMap` generic
+- [ ] Rewrite planar operations in order to be AD-compatible
+- [ ] Fix rrules for fermionic tensors
+- [ ] Fix GPU support
+- [ ] Proper threading support
+- [ ] Rewrite documentation
+
+# Changelog

From 79fcfa51df4a9492e05a1934e965dfbfab2ab2a3 Mon Sep 17 00:00:00 2001
From: Lukas <37111893+lkdvos@users.noreply.github.com>
Date: Thu, 21 Mar 2024 09:53:20 +0100
Subject: [PATCH 03/21] Change type signature:
 `AbstractTensorMap{<:Number,<:IndexSpace,N1,N2}` (#1)

* Add scalartype parameter to `AbstractTensorMap`

* Add scalartype in definition of TensorMap

* Add scalartype to AdjointTensorMap

* Add scalartype to BraidingTensor

* Implement changes for indexmanipulations

* Implement scalartype for linalg

* implement changes for tensoroperations

* implement changes for planaroperations

* implement changes for ChainRules

* small fixes

* Fix typo

* Fix missing change

* Formatter

* Apply scalartype changes to tests

* Fix some ambiguities

* De-specialize `trace_permute`

* Despecialize `_contract!`

* Remove unused type-parameter

* Despecialize `planaradd!` and `planartrace!`

* Small fixes

* Change braidingtensor syntax to resolve ambiguities

* despecialize some linalg methods

* formatter

* Add changelog entry

* Fix ambiguity on Julia1.6
---
 Changelog.md                      |   7 ++
 ext/TensorKitChainRulesCoreExt.jl |   8 +-
 src/planar/planaroperations.jl    |  63 +++++++------
 src/tensors/abstracttensor.jl     |  66 +++++++-------
 src/tensors/adjoint.jl            |  30 +++----
 src/tensors/braidingtensor.jl     | 128 ++++++++++++--------------
 src/tensors/indexmanipulations.jl | 100 ++++++++++-----------
 src/tensors/linalg.jl             |  18 ++--
 src/tensors/tensor.jl             | 143 ++++++++++++++++--------------
 src/tensors/tensoroperations.jl   | 111 +++++++++++++----------
 test/planar.jl                    |   4 +-
 11 files changed, 352 insertions(+), 326 deletions(-)

diff --git a/Changelog.md b/Changelog.md
index 295d6a18..12a5874a 100644
--- a/Changelog.md
+++ b/Changelog.md
@@ -13,3 +13,10 @@ Features that are planned to be implemented before the release of v1.0.0, in no
 - [ ] Rewrite documentation
 
 # Changelog
+
+### `AbstractTensorMap{E,S,N₁,N₂}`
+
+This adds the scalar type as a parameter to the `AbstractTensorMap` type. This is useful in
+the contexts where different types of tensors are used (`AdjointTensor`, `BraidingTensor`,
+...), which still have the same scalartype. Additionally, this removes many specializations
+for methods in order to reduce ambiguity errors.
diff --git a/ext/TensorKitChainRulesCoreExt.jl b/ext/TensorKitChainRulesCoreExt.jl
index b36976fd..cb38d009 100644
--- a/ext/TensorKitChainRulesCoreExt.jl
+++ b/ext/TensorKitChainRulesCoreExt.jl
@@ -22,8 +22,8 @@ function _repartition(p::Union{IndexTuple,Index2Tuple}, ::Index2Tuple{N₁}) whe
     return _repartition(p, N₁)
 end
 function _repartition(p::Union{IndexTuple,Index2Tuple},
-                      ::AbstractTensorMap{<:Any,N₁}) where {N₁}
-    return _repartition(p, N₁)
+                      t::AbstractTensorMap)
+    return _repartition(p, numout(t))
 end
 
 TensorKit.block(t::ZeroTangent, c::Sector) = t
@@ -59,8 +59,8 @@ function ChainRulesCore.rrule(::typeof(Base.copy), t::AbstractTensorMap)
 end
 
 ChainRulesCore.ProjectTo(::T) where {T<:AbstractTensorMap} = ProjectTo{T}()
-function (::ProjectTo{T1})(x::T2) where {S,N1,N2,T1<:AbstractTensorMap{S,N1,N2},
-                                         T2<:AbstractTensorMap{S,N1,N2}}
+function (::ProjectTo{T1})(x::T2) where {S,N1,N2,T1<:AbstractTensorMap{<:Any,S,N1,N2},
+                                         T2<:AbstractTensorMap{<:Any,S,N1,N2}}
     T1 === T2 && return x
     y = similar(x, scalartype(T1))
     for (c, b) in blocks(y)
diff --git a/src/planar/planaroperations.jl b/src/planar/planaroperations.jl
index ee23bcd7..1d763631 100644
--- a/src/planar/planaroperations.jl
+++ b/src/planar/planaroperations.jl
@@ -1,34 +1,40 @@
 # planar versions of tensor operations add!, trace! and contract!
-function planaradd!(C::AbstractTensorMap{S,N₁,N₂},
-                    A::AbstractTensorMap{S},
-                    p::Index2Tuple{N₁,N₂},
-                    α::Number,
-                    β::Number,
-                    backend::Backend...) where {S,N₁,N₂}
+function planaradd!(C::AbstractTensorMap,
+                    A::AbstractTensorMap,
+                    p::Index2Tuple,
+                    α::Number, β::Number,
+                    backend::Backend...)
     return add_transpose!(C, A, p, α, β, backend...)
 end
 
-function planartrace!(C::AbstractTensorMap{S,N₁,N₂},
-                      A::AbstractTensorMap{S},
-                      p::Index2Tuple{N₁,N₂},
-                      q::Index2Tuple{N₃,N₃},
+function planartrace!(C::AbstractTensorMap,
+                      A::AbstractTensorMap,
+                      (p₁, p₂)::Index2Tuple,
+                      (q₁, q₂)::Index2Tuple,
                       α::Number,
                       β::Number,
-                      backend::Backend...) where {S,N₁,N₂,N₃}
+                      backend::Backend...)
+    (S = spacetype(C)) == spacetype(A) ||
+        throw(SpaceMismatch("incompatible spacetypes"))
     if BraidingStyle(sectortype(S)) == Bosonic()
         return trace_permute!(C, A, p, q, α, β, backend...)
     end
+    (N₃ = length(q₁)) == length(q₂) ||
+        throw(IndexError("number of trace indices does not match"))
+    N₁, N₂ = length(p₁), length(p₂)
 
     @boundscheck begin
-        all(i -> space(A, p[1][i]) == space(C, i), 1:N₁) ||
-            throw(SpaceMismatch("trace: A = $(codomain(A))←$(domain(A)),
-                    C = $(codomain(C))←$(domain(C)), p1 = $(p1), p2 = $(p2)"))
-        all(i -> space(A, p[2][i]) == space(C, N₁ + i), 1:N₂) ||
-            throw(SpaceMismatch("trace: A = $(codomain(A))←$(domain(A)),
-                    C = $(codomain(C))←$(domain(C)), p1 = $(p1), p2 = $(p2)"))
-        all(i -> space(A, q[1][i]) == dual(space(A, q[2][i])), 1:N₃) ||
-            throw(SpaceMismatch("trace: A = $(codomain(A))←$(domain(A)),
-                    q1 = $(q1), q2 = $(q2)"))
+        numout(C) == N₁ || throw(IndexError("number of output indices does not match"))
+        numin(C) == N₂ || throw(IndexError("number of input indices does not match"))
+        all(i -> space(A, p₁[i]) == space(C, i), 1:N₁) ||
+            throw(SpaceMismatch("trace: A = $(space(A)),
+                    C = $(space(C)), p₁ = $(p₁), p₂ = $(p₂)"))
+        all(i -> space(A, p₂[i]) == space(C, N₁ + i), 1:N₂) ||
+            throw(SpaceMismatch("trace: A = $(space(A)),
+                    C = $(space(C)), p₁ = $(p₁), p₂ = $(p₂)"))
+        all(i -> space(A, q₁[i]) == dual(space(A, q₂[i])), 1:N₃) ||
+            throw(SpaceMismatch("trace: A = $(space(A)),
+                    q1 = $(q₁), q2 = $(q₂)"))
     end
 
     if iszero(β)
@@ -37,23 +43,24 @@ function planartrace!(C::AbstractTensorMap{S,N₁,N₂},
         rmul!(C, β)
     end
     for (f₁, f₂) in fusiontrees(A)
-        for ((f₁′, f₂′), coeff) in planar_trace(f₁, f₂, p..., q...)
-            TO.tensortrace!(C[f₁′, f₂′], p, A[f₁, f₂], q, :N, α * coeff, true, backend...)
+        for ((f₁′, f₂′), coeff) in planar_trace(f₁, f₂, p₁, p₂, q₁, q₂)
+            TO.tensortrace!(C[f₁′, f₂′], (p₁, p₂), A[f₁, f₂], (q₁, q₂), :N, α * coeff, true,
+                            backend...)
         end
     end
     return C
 end
 
-function planarcontract!(C::AbstractTensorMap{S,N₁,N₂},
-                         A::AbstractTensorMap{S},
+function planarcontract!(C::AbstractTensorMap,
+                         A::AbstractTensorMap,
                          pA::Index2Tuple,
-                         B::AbstractTensorMap{S},
+                         B::AbstractTensorMap,
                          pB::Index2Tuple,
-                         pAB::Index2Tuple{N₁,N₂},
+                         pAB::Index2Tuple,
                          α::Number,
                          β::Number,
-                         backend::Backend...) where {S,N₁,N₂}
-    if BraidingStyle(sectortype(S)) == Bosonic()
+                         backend::Backend...)
+    if BraidingStyle(sectortype(C)) == Bosonic()
         return contract!(C, A, pA, B, pB, pAB, α, β, backend...)
     end
 
diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl
index 55195e4f..447538b4 100644
--- a/src/tensors/abstracttensor.jl
+++ b/src/tensors/abstracttensor.jl
@@ -3,45 +3,45 @@
 # Abstract Tensor type
 #----------------------
 """
-    abstract type AbstractTensorMap{S<:IndexSpace, N₁, N₂} end
+    abstract type AbstractTensorMap{E<:Number, S<:IndexSpace, N₁, N₂} end
 
-Abstract supertype of all tensor maps, i.e. linear maps between tensor products
-of vector spaces of type `S<:IndexSpace`. An `AbstractTensorMap` maps from
-an input space of type `ProductSpace{S, N₂}` to an output space of type
-`ProductSpace{S, N₁}`.
+Abstract supertype of all tensor maps, i.e. linear maps between tensor products of vector
+spaces of type `S<:IndexSpace`, with element type `E`. An `AbstractTensorMap` maps from an
+input space of type `ProductSpace{S, N₂}` to an output space of type `ProductSpace{S, N₁}`.
 """
-abstract type AbstractTensorMap{S<:IndexSpace,N₁,N₂} end
+abstract type AbstractTensorMap{E<:Number,S<:IndexSpace,N₁,N₂} end
+
 """
-    AbstractTensor{S<:IndexSpace, N} = AbstractTensorMap{S, N, 0}
+    AbstractTensor{E,S,N} = AbstractTensorMap{E,S,N,0}
 
-Abstract supertype of all tensors, i.e. elements in the tensor product space
-of type `ProductSpace{S, N}`, built from elementary spaces of type `S<:IndexSpace`.
+Abstract supertype of all tensors, i.e. elements in the tensor product space of type
+`ProductSpace{S, N}`, with element type `E`.
 
-An `AbstractTensor{S, N}` is actually a special case `AbstractTensorMap{S, N, 0}`,
-i.e. a tensor map with only a non-trivial output space.
+An `AbstractTensor{E, S, N}` is actually a special case `AbstractTensorMap{E, S, N, 0}`,
+i.e. a tensor map with only non-trivial output spaces.
 """
-const AbstractTensor{S<:IndexSpace,N} = AbstractTensorMap{S,N,0}
+const AbstractTensor{E,S,N} = AbstractTensorMap{E,S,N,0}
 
 # tensor characteristics
 #------------------------
-Base.eltype(::Union{T,Type{T}}) where {T<:AbstractTensorMap} = scalartype(T)
+Base.eltype(::Type{<:AbstractTensorMap{E}}) where {E} = E
 
 """
     spacetype(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Type{S<:IndexSpace}
 
 Return the type of the elementary space `S` of a tensor.
 """
-spacetype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = S
+spacetype(::Type{<:AbstractTensorMap{E,S}}) where {E,S<:IndexSpace} = S
 
 """
     sectortype(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Type{I<:Sector}
 
 Return the type of sector `I` of a tensor.
 """
-sectortype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = sectortype(S)
+sectortype(::Type{T}) where {T<:AbstractTensorMap} = sectortype(spacetype(T))
 
-function InnerProductStyle(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace}
-    return InnerProductStyle(S)
+function InnerProductStyle(::Type{T}) where {T<:AbstractTensorMap}
+    return InnerProductStyle(spacetype(T))
 end
 
 """
@@ -49,7 +49,7 @@ end
 
 Return the type of field `𝕂` of a tensor.
 """
-field(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = field(S)
+field(::Type{T}) where {T<:AbstractTensorMap} = field(spacetype(T))
 
 """
     numout(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Int
@@ -58,7 +58,7 @@ Return the number of output spaces of a tensor. This is equivalent to the number
 
 See also [`numin`](@ref) and [`numind`](@ref).
 """
-numout(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₁
+numout(::Type{<:AbstractTensorMap{E,S,N₁}}) where {E,S,N₁} = N₁
 
 """
     numin(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Int
@@ -67,7 +67,7 @@ Return the number of input spaces of a tensor. This is equivalent to the number
 
 See also [`numout`](@ref) and [`numind`](@ref).
 """
-numin(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₂
+numin(::Type{<:AbstractTensorMap{E,S,N₁,N₂}}) where {E,S,N₁,N₂} = N₂
 
 """
     numind(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Int
@@ -77,7 +77,7 @@ total number of spaces in the domain and codomain of that tensor.
 
 See also [`numout`](@ref) and [`numin`](@ref).
 """
-numind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₁ + N₂
+numind(::Type{T}) where {T<:AbstractTensorMap} = numin(T) + numout(T)
 
 function similarstoragetype(TT::Type{<:AbstractTensorMap}, ::Type{T}) where {T}
     return Core.Compiler.return_type(similar, Tuple{storagetype(TT),Type{T}})
@@ -97,7 +97,7 @@ similarstoragetype(t::AbstractTensorMap, T) = similarstoragetype(typeof(t), T)
 const order = numind
 
 @doc """
-    codomain(t::AbstractTensorMap{S,N₁,N₂}, [i::Int]) -> ProductSpace{S,N₁}
+    codomain(t::AbstractTensorMap{E,S,N₁,N₂}, [i::Int]) -> ProductSpace{S,N₁}
 
 Return the codomain of a tensor, i.e. the product space of the output spaces. If `i` is
 specified, return the `i`-th output space. Implementations should provide `codomain(t)`.
@@ -109,7 +109,7 @@ codomain(t::AbstractTensorMap, i) = codomain(t)[i]
 target(t::AbstractTensorMap) = codomain(t) # categorical terminology
 
 @doc """
-    domain(t::AbstractTensorMap{S,N₁,N₂}, [i::Int]) -> ProductSpace{S,N₂}
+    domain(t::AbstractTensorMap{E,S,N₁,N₂}, [i::Int]) -> ProductSpace{S,N₂}
 
 Return the domain of a tensor, i.e. the product space of the input spaces. If `i` is
 specified, return the `i`-th input space. Implementations should provide `domain(t)`.
@@ -121,7 +121,7 @@ domain(t::AbstractTensorMap, i) = domain(t)[i]
 source(t::AbstractTensorMap) = domain(t) # categorical terminology
 
 """
-    space(t::AbstractTensorMap{S,N₁,N₂}, [i::Int]) -> HomSpace{S,N₁,N₂}
+    space(t::AbstractTensorMap{E,S,N₁,N₂}, [i::Int]) -> HomSpace{S,N₁,N₂}
 
 The index information of a tensor, i.e. the `HomSpace` of its domain and codomain. If `i` is specified, return the `i`-th index space.
 """
@@ -143,8 +143,8 @@ Return all indices of the codomain of a tensor.
 
 See also [`domainind`](@ref) and [`allind`](@ref).
 """
-function codomainind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂}
-    return ntuple(n -> n, N₁)
+function codomainind(::Type{T}) where {T<:AbstractTensorMap}
+    return ntuple(identity, numout(T))
 end
 codomainind(t::AbstractTensorMap) = codomainind(typeof(t))
 
@@ -155,8 +155,8 @@ Return all indices of the domain of a tensor.
 
 See also [`codomainind`](@ref) and [`allind`](@ref).
 """
-function domainind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂}
-    return ntuple(n -> N₁ + n, N₂)
+function domainind(::Type{T}) where {T<:AbstractTensorMap}
+    return ntuple(n -> numout(T) + n, numin(T))
 end
 domainind(t::AbstractTensorMap) = domainind(typeof(t))
 
@@ -167,13 +167,13 @@ Return all indices of a tensor, i.e. the indices of its domain and codomain.
 
 See also [`codomainind`](@ref) and [`domainind`](@ref).
 """
-function allind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂}
-    return ntuple(n -> n, N₁ + N₂)
+function allind(::Type{T}) where {T<:AbstractTensorMap}
+    return ntuple(n -> n, numind(T))
 end
 allind(t::AbstractTensorMap) = allind(typeof(t))
 
-function adjointtensorindex(::AbstractTensorMap{<:IndexSpace,N₁,N₂}, i) where {N₁,N₂}
-    return ifelse(i <= N₁, N₂ + i, i - N₁)
+function adjointtensorindex(t::AbstractTensorMap, i)
+    return ifelse(i <= numout(t), numin(t) + i, i - numout(t))
 end
 
 function adjointtensorindices(t::AbstractTensorMap, indices::IndexTuple)
@@ -258,7 +258,7 @@ end
 # Conversion to Array:
 #----------------------
 # probably not optimized for speed, only for checking purposes
-function Base.convert(::Type{Array}, t::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂}
+function Base.convert(::Type{Array}, t::AbstractTensorMap)
     I = sectortype(t)
     if I === Trivial
         convert(Array, t[])
diff --git a/src/tensors/adjoint.jl b/src/tensors/adjoint.jl
index 43ad110d..073361bf 100644
--- a/src/tensors/adjoint.jl
+++ b/src/tensors/adjoint.jl
@@ -1,19 +1,19 @@
 # AdjointTensorMap: lazy adjoint
 #==========================================================#
 """
-    struct AdjointTensorMap{S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{S, N₁, N₂}
+    struct AdjointTensorMap{E, S, N₁, N₂, ...} <: AbstractTensorMap{E, S, N₁, N₂}
 
 Specific subtype of [`AbstractTensorMap`](@ref) that is a lazy wrapper for representing the
 adjoint of an instance of [`TensorMap`](@ref).
 """
-struct AdjointTensorMap{S<:IndexSpace,N₁,N₂,I<:Sector,A,F₁,F₂} <:
-       AbstractTensorMap{S,N₁,N₂}
-    parent::TensorMap{S,N₂,N₁,I,A,F₂,F₁}
+struct AdjointTensorMap{E,S,N₁,N₂,I,A,F₁,F₂} <:
+       AbstractTensorMap{E,S,N₁,N₂}
+    parent::TensorMap{E,S,N₂,N₁,I,A,F₂,F₁}
 end
 
 #! format: off
-const AdjointTrivialTensorMap{S<:IndexSpace,N₁,N₂,A<:DenseMatrix} =
-    AdjointTensorMap{S,N₁,N₂,Trivial,A,Nothing,Nothing}
+const AdjointTrivialTensorMap{E,S,N₁,N₂,A<:DenseMatrix} =
+    AdjointTensorMap{E,S,N₁,N₂,Trivial,A,Nothing,Nothing}
 #! format: on
 
 # Constructor: construct from taking adjoint of a tensor
@@ -29,8 +29,8 @@ domain(t::AdjointTensorMap) = codomain(t.parent)
 blocksectors(t::AdjointTensorMap) = blocksectors(t.parent)
 
 #! format: off
-storagetype(::Type{<:AdjointTensorMap{<:IndexSpace,N₁,N₂,Trivial,A}}) where {N₁,N₂,A<:DenseMatrix} = A
-storagetype(::Type{<:AdjointTensorMap{<:IndexSpace,N₁,N₂,I,<:SectorDict{I,A}}}) where {N₁, N₂,I<:Sector,A<:DenseMatrix} = A
+storagetype(::Type{<:AdjointTrivialTensorMap{E,S,N₁,N₂,A}}) where {E,S,N₁,N₂,A<:DenseMatrix} = A
+storagetype(::Type{<:AdjointTensorMap{E,S,N₁,N₂,I,<:SectorDict{I,A}}}) where {E,S,N₁,N₂,I<:Sector,A<:DenseMatrix} = A
 #! format: on
 
 dim(t::AdjointTensorMap) = dim(t.parent)
@@ -44,8 +44,8 @@ blocks(t::AdjointTensorMap) = (c => b' for (c, b) in blocks(t.parent))
 fusiontrees(::AdjointTrivialTensorMap) = ((nothing, nothing),)
 fusiontrees(t::AdjointTensorMap) = TensorKeyIterator(t.parent.colr, t.parent.rowr)
 
-function Base.getindex(t::AdjointTensorMap{S,N₁,N₂,I},
-                       f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {S,N₁,N₂,I}
+function Base.getindex(t::AdjointTensorMap{E,S,N₁,N₂,I},
+                       f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂,I}
     c = f₁.coupled
     @boundscheck begin
         c == f₂.coupled || throw(SectorMismatch())
@@ -55,9 +55,9 @@ function Base.getindex(t::AdjointTensorMap{S,N₁,N₂,I},
                                                    t.parent.colr[c][f₁]])',
                     (dims(codomain(t), f₁.uncoupled)..., dims(domain(t), f₂.uncoupled)...))
 end
-@propagate_inbounds function Base.setindex!(t::AdjointTensorMap{S,N₁,N₂}, v,
+@propagate_inbounds function Base.setindex!(t::AdjointTensorMap{E,S,N₁,N₂,I}, v,
                                             f₁::FusionTree{I,N₁},
-                                            f₂::FusionTree{I,N₂}) where {S,N₁,N₂,I}
+                                            f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂,I}
     return copy!(getindex(t, f₁, f₂), v)
 end
 
@@ -90,16 +90,16 @@ end
 function Base.summary(t::AdjointTensorMap)
     return print("AdjointTensorMap(", codomain(t), " ← ", domain(t), ")")
 end
-function Base.show(io::IO, t::AdjointTensorMap{S}) where {S<:IndexSpace}
+function Base.show(io::IO, t::AdjointTensorMap)
     if get(io, :compact, false)
         print(io, "AdjointTensorMap(", codomain(t), " ← ", domain(t), ")")
         return
     end
     println(io, "AdjointTensorMap(", codomain(t), " ← ", domain(t), "):")
-    if sectortype(S) == Trivial
+    if sectortype(t) === Trivial
         Base.print_array(io, t[])
         println(io)
-    elseif FusionStyle(sectortype(S)) isa UniqueFusion
+    elseif FusionStyle(sectortype(t)) isa UniqueFusion
         for (f₁, f₂) in fusiontrees(t)
             println(io, "* Data for sector ", f₁.uncoupled, " ← ", f₂.uncoupled, ":")
             Base.print_array(io, t[f₁, f₂])
diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl
index 4b8489d6..f5924bf0 100644
--- a/src/tensors/braidingtensor.jl
+++ b/src/tensors/braidingtensor.jl
@@ -2,7 +2,7 @@
 # special (2,2) tensor that implements a standard braiding operation
 #====================================================================#
 """
-    struct BraidingTensor{S<:IndexSpace} <: AbstractTensorMap{S, 2, 2}
+    struct BraidingTensor{E,S<:IndexSpace} <: AbstractTensorMap{E, S, 2, 2}
     BraidingTensor(V1::S, V2::S, adjoint::Bool=false) where {S<:IndexSpace}
 
 Specific subtype of [`AbstractTensorMap`](@ref) for representing the braiding tensor that
@@ -11,12 +11,11 @@ braids the first input over the second input; its inverse can be obtained as the
 It holds that `domain(BraidingTensor(V1, V2)) == V1 ⊗ V2` and
 `codomain(BraidingTensor(V1, V2)) == V2 ⊗ V1`.
 """
-struct BraidingTensor{S<:IndexSpace,A} <: AbstractTensorMap{S,2,2}
+struct BraidingTensor{E,S} <: AbstractTensorMap{E,S,2,2}
     V1::S
     V2::S
     adjoint::Bool
-    function BraidingTensor{S,A}(V1::S, V2::S,
-                                 adjoint::Bool=false) where {S<:IndexSpace,A<:DenseMatrix}
+    function BraidingTensor{E,S}(V1::S, V2::S, adjoint::Bool=false) where {E,S<:IndexSpace}
         for a in sectors(V1)
             for b in sectors(V2)
                 for c in (a ⊗ b)
@@ -25,30 +24,31 @@ struct BraidingTensor{S<:IndexSpace,A} <: AbstractTensorMap{S,2,2}
                 end
             end
         end
-        return new{S,A}(V1, V2, adjoint)
+        return new{E,S}(V1, V2, adjoint)
         # partial construction: only construct rowr and colr when needed
     end
 end
 function BraidingTensor(V1::S, V2::S, adjoint::Bool=false) where {S<:IndexSpace}
     if BraidingStyle(sectortype(S)) isa SymmetricBraiding
-        return BraidingTensor{S,Matrix{Float64}}(V1, V2, adjoint)
+        return BraidingTensor{Float64,S}(V1, V2, adjoint)
     else
-        return BraidingTensor{S,Matrix{ComplexF64}}(V1, V2, adjoint)
+        return BraidingTensor{ComplexF64,S}(V1, V2, adjoint)
     end
 end
-function BraidingTensor(V::HomSpace{S}, adjoint::Bool=false) where {S<:IndexSpace}
+function BraidingTensor(V::HomSpace, adjoint::Bool=false)
     domain(V) == reverse(codomain(V)) ||
         throw(SpaceMismatch("Cannot define a braiding on $V"))
     return BraidingTensor(V[1], V[2], adjoint)
 end
-function Base.adjoint(b::BraidingTensor{S,A}) where {S<:IndexSpace,A<:DenseMatrix}
-    return BraidingTensor{S,A}(b.V1, b.V2, !b.adjoint)
+function Base.adjoint(b::BraidingTensor{E,S}) where {E,S}
+    return BraidingTensor{E,S}(b.V1, b.V2, !b.adjoint)
 end
 
 domain(b::BraidingTensor) = b.adjoint ? b.V2 ⊗ b.V1 : b.V1 ⊗ b.V2
 codomain(b::BraidingTensor) = b.adjoint ? b.V1 ⊗ b.V2 : b.V2 ⊗ b.V1
 
-storagetype(::Type{BraidingTensor{S,A}}) where {S<:IndexSpace,A<:DenseMatrix} = A
+# TODO: check if this can be removed/is correct
+storagetype(::Type{BraidingTensor{E,S}}) where {E,S} = Matrix{E}
 
 blocksectors(b::BraidingTensor) = blocksectors(b.V1 ⊗ b.V2)
 hasblock(b::BraidingTensor, s::Sector) = s ∈ blocksectors(b)
@@ -87,8 +87,8 @@ function fusiontrees(b::BraidingTensor)
     return TensorKeyIterator(rowr, colr)
 end
 
-function Base.getindex(b::BraidingTensor{S}) where {S}
-    sectortype(S) == Trivial || throw(SectorMismatch())
+function Base.getindex(b::BraidingTensor)
+    sectortype(b) === Trivial || throw(SectorMismatch())
     (V1, V2) = domain(b)
     d = (dim(V2), dim(V1), dim(V1), dim(V2))
     return sreshape(StridedView(block(b, Trivial())), d)
@@ -192,13 +192,13 @@ blocks(b::BraidingTensor) = blocks(TensorMap(b))
 # Index manipulations
 # -------------------
 has_shared_permute(t::BraidingTensor, args...) = false
-function add_transform!(tdst::AbstractTensorMap{S,N₁,N₂},
-                        tsrc::BraidingTensor{S},
-                        (p₁, p₂)::Index2Tuple{N₁,N₂},
+function add_transform!(tdst::AbstractTensorMap,
+                        tsrc::BraidingTensor,
+                        (p₁, p₂)::Index2Tuple,
                         fusiontreetransform,
                         α::Number,
                         β::Number,
-                        backend::Backend...) where {S,N₁,N₂}
+                        backend::Backend...)
     return add_transform!(tdst, copy(tsrc), (p₁, p₂), fusiontreetransform, α, β, backend...)
 end
 
@@ -209,30 +209,37 @@ end
 # TensorOperations
 # ----------------
 # TODO: implement specialized methods
-function TO.tensoradd!(C::AbstractTensorMap{S,N₁,N₂}, pC::Index2Tuple{N₁,N₂},
-                       A::BraidingTensor{S}, conjA::Symbol, α::Number, β::Number,
-                       backend::Backend...) where {S,N₁,N₂}
+
+function TO.tensoradd!(C::AbstractTensorMap, pC::Index2Tuple,
+                       A::BraidingTensor, conjA::Symbol, α::Number, β::Number,
+                       backend::Backend...)
     return TO.tensoradd!(C, pC, copy(A), conjA, α, β, backend...)
 end
 
 # Planar operations
 # -----------------
-function planaradd!(C::AbstractTensorMap{S,N₁,N₂},
-                    A::BraidingTensor{S},
-                    p::Index2Tuple{N₁,N₂},
+# TODO: implement specialized methods
+
+function planaradd!(C::AbstractTensorMap,
+                    A::BraidingTensor, p::Index2Tuple,
                     α::Number, β::Number,
-                    backend::Backend...) where {S,N₁,N₂}
+                    backend::Backend...)
     return planaradd!(C, copy(A), p, α, β, backend...)
 end
 
-function planarcontract!(C::AbstractTensorMap{S,N₁,N₂},
-                         A::BraidingTensor{S},
-                         (oindA, cindA)::Index2Tuple{2,2},
-                         B::AbstractTensorMap{S},
-                         (cindB, oindB)::Index2Tuple{2,N₃},
-                         (p1, p2)::Index2Tuple{N₁,N₂},
+function planarcontract!(C::AbstractTensorMap,
+                         A::BraidingTensor,
+                         (oindA, cindA)::Index2Tuple,
+                         B::AbstractTensorMap,
+                         (cindB, oindB)::Index2Tuple,
+                         (p1, p2)::Index2Tuple,
                          α::Number, β::Number,
-                         backend::Backend...) where {S,N₁,N₂,N₃}
+                         backend::Backend...)
+    # special case only defined for contracting 2 indices
+    length(oindA) == length(cindA) == 2 ||
+        return planarcontract!(C, copy(A), (oindA, cindA), B, (cindB, oindB), (p1, p2),
+                               α, β, backend...)
+
     codA, domA = codomainind(A), domainind(A)
     codB, domB = codomainind(B), domainind(B)
     oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA,
@@ -271,14 +278,19 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂},
     end
     return C
 end
-function planarcontract!(C::AbstractTensorMap{S,N₁,N₂},
-                         A::AbstractTensorMap{S},
-                         (oindA, cindA)::Index2Tuple{N₃,2},
-                         B::BraidingTensor{S},
-                         (cindB, oindB)::Index2Tuple{2,2},
-                         (p1, p2)::Index2Tuple{N₁,N₂},
+function planarcontract!(C::AbstractTensorMap,
+                         A::AbstractTensorMap,
+                         (oindA, cindA)::Index2Tuple,
+                         B::BraidingTensor,
+                         (cindB, oindB)::Index2Tuple,
+                         (p1, p2)::Index2Tuple,
                          α::Number, β::Number,
-                         backend::Backend...) where {S,N₁,N₂,N₃}
+                         backend::Backend...)
+    # special case only defined for contracting 2 indices
+    length(oindB) == length(cindB) == 2 ||
+        return planarcontract!(C, A, (oindA, cindA), copy(B), (cindB, oindB), (p1, p2),
+                               α, β, backend...)
+
     codA, domA = codomainind(A), domainind(A)
     codB, domB = codomainind(B), domainind(B)
     oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA,
@@ -317,41 +329,19 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂},
     end
     return C
 end
-function planarcontract!(C::AbstractTensorMap{S,N₁,N₂},
-                         A::BraidingTensor{S},
-                         (oindA, cindA)::Index2Tuple{2,2},
-                         B::BraidingTensor{S},
-                         (cindB, oindB)::Index2Tuple{2,2},
-                         (p1, p2)::Index2Tuple{N₁,N₂},
-                         α::Number, β::Number,
-                         backend::Backend...) where {S,N₁,N₂}
-    return planarcontract!(C, copy(A), (oindA, cindA), B, (cindB, oindB), (p1, p2), α, β,
-                           backend...)
-end
 
-# Fallback cases for planarcontract!
-# TODO: implement specialised cases for contracting 0, 1, 3 and 4 indices
-function planarcontract!(C::AbstractTensorMap{S}, A::BraidingTensor{S}, pA::Index2Tuple,
-                         B::BraidingTensor{S}, pB::Index2Tuple, α::Number, β::Number,
-                         backend::Backend...) where {S}
-    return planarcontract!(C, copy(A), pA, copy(B), pB, α, β, backend...)
-end
-function planarcontract!(C::AbstractTensorMap{S}, A::BraidingTensor{S}, pA::Index2Tuple,
-                         B::AbstractTensorMap{S}, pB::Index2Tuple, α::Number, β::Number,
-                         backend::Backend...) where {S}
-    return planarcontract!(C, copy(A), pA, B, pB, α, β, backend...)
-end
-function planarcontract!(C::AbstractTensorMap{S}, A::AbstractTensorMap{S}, pA::Index2Tuple,
-                         B::BraidingTensor{S}, pB::Index2Tuple, α::Number, β::Number,
-                         backend::Backend...) where {S}
-    return planarcontract!(C, A, pA, copy(B), pB, α, β, backend...)
+# ambiguity fix:
+function planarcontract!(C::AbstractTensorMap, A::BraidingTensor, pA::Index2Tuple,
+                         B::BraidingTensor, pB::Index2Tuple, pC::Index2Tuple,
+                         α::Number, β::Number, backend::Backend...)
+    return planarcontract!(C, copy(A), pA, copy(B), pB, pC, α, β, backend...)
 end
 
-function planartrace!(C::AbstractTensorMap{S,N₁,N₂},
-                      A::BraidingTensor{S},
-                      p::Index2Tuple{N₁,N₂}, q::Index2Tuple{N₃,N₃},
+function planartrace!(C::AbstractTensorMap,
+                      A::BraidingTensor,
+                      p::Index2Tuple, q::Index2Tuple,
                       α::Number, β::Number,
-                      backend::Backend...) where {S,N₁,N₂,N₃}
+                      backend::Backend...)
     return planartrace!(C, copy(A), p, q, α, β, backend...)
 end
 
diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl
index dccb38a2..d6ebd3ea 100644
--- a/src/tensors/indexmanipulations.jl
+++ b/src/tensors/indexmanipulations.jl
@@ -1,8 +1,7 @@
 # Index manipulations
 #---------------------
 """
-    permute!(tdst::AbstractTensorMap{S,N₁,N₂}, tsrc::AbstractTensorMap{S},
-             (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}) where {S,N₁,N₂}
+    permute!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple)
         -> tdst
 
 Write into `tdst` the result of permuting the indices of `tsrc`.
@@ -10,16 +9,16 @@ The codomain and domain of `tdst` correspond to the indices in `p₁` and `p₂`
                 
 See [`permute`](@ref) for creating a new tensor and [`add_permute!`](@ref) for a more general version.
 """
-@propagate_inbounds function Base.permute!(tdst::AbstractTensorMap{S,N₁,N₂},
-                                           tsrc::AbstractTensorMap{S},
-                                           p::Index2Tuple{N₁,N₂}) where {S,N₁,N₂}
+@propagate_inbounds function Base.permute!(tdst::AbstractTensorMap,
+                                           tsrc::AbstractTensorMap,
+                                           p::Index2Tuple)
     return add_permute!(tdst, tsrc, p, true, false)
 end
 
 """
-    permute(tsrc::AbstractTensorMap{S}, (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}};
-            copy::Bool=false) where {S,N₁,N₂}
-        -> tdst::TensorMap{S,N₁,N₂}
+    permute(tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple;
+            copy::Bool=false)
+        -> tdst::TensorMap
 
 Return tensor `tdst` obtained by permuting the indices of `tsrc`.
 The codomain and domain of `tdst` correspond to the indices in `p₁` and `p₂` of `tsrc` respectively.
@@ -28,10 +27,10 @@ If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwis
 
 To permute into an existing destination, see [permute!](@ref) and [`add_permute!`](@ref)
 """
-function permute(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple{N₁,N₂};
-                 copy::Bool=false) where {S,N₁,N₂}
-    cod = ProductSpace{S,N₁}(map(n -> space(t, n), p₁))
-    dom = ProductSpace{S,N₂}(map(n -> dual(space(t, n)), p₂))
+function permute(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂};
+                 copy::Bool=false) where {N₁,N₂}
+    cod = ProductSpace{spacetype(t),N₁}(map(n -> space(t, n), p₁))
+    dom = ProductSpace{spacetype(t),N₂}(map(n -> dual(space(t, n)), p₂))
     # share data if possible
     if !copy
         if p₁ === codomainind(t) && p₂ === domainind(t)
@@ -45,14 +44,13 @@ function permute(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple{N₁,N₂};
         return permute!(similar(t, cod ← dom), t, (p₁, p₂))
     end
 end
-function permute(t::AdjointTensorMap{S}, (p₁, p₂)::Index2Tuple;
-                 copy::Bool=false) where {S}
+function permute(t::AdjointTensorMap, (p₁, p₂)::Index2Tuple; copy::Bool=false)
     p₁′ = adjointtensorindices(t, p₂)
     p₂′ = adjointtensorindices(t, p₁)
-    return adjoint(permute(adjoint(t), (p₁′, p₂′); copy=copy))
+    return adjoint(permute(adjoint(t), (p₁′, p₂′); copy))
 end
 function permute(t::AbstractTensorMap, p::IndexTuple; copy::Bool=false)
-    return permute(t, (p, ()); copy=copy)
+    return permute(t, (p, ()); copy)
 end
 
 function has_shared_permute(t::TensorMap, (p₁, p₂)::Index2Tuple)
@@ -76,8 +74,8 @@ end
 
 # Braid
 """
-    braid!(tdst::AbstractTensorMap{S,N₁,N₂}, tsrc::AbstractTensorMap{S},
-           (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}, levels::Tuple) where {S,N₁,N₂}
+    braid!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap,
+           (p₁, p₂)::Index2Tuple, levels::Tuple)
         -> tdst
 
 Write into `tdst` the result of braiding the indices of `tsrc`.
@@ -87,17 +85,17 @@ which determines whether they will braid over or under any other index with whic
 
 See [`braid`](@ref) for creating a new tensor and [`add_braid!`](@ref) for a more general version.
 """
-@propagate_inbounds function braid!(tdst::AbstractTensorMap{S,N₁,N₂},
-                                    tsrc::AbstractTensorMap{S},
-                                    (p₁, p₂)::Index2Tuple{N₁,N₂},
-                                    levels::IndexTuple) where {S,N₁,N₂}
-    return add_braid!(tdst, tsrc, (p₁, p₂), levels, true, false)
+@propagate_inbounds function braid!(tdst::AbstractTensorMap,
+                                    tsrc::AbstractTensorMap,
+                                    p::Index2Tuple,
+                                    levels::IndexTuple)
+    return add_braid!(tdst, tsrc, p, levels, true, false)
 end
 
 """
-    braid(tsrc::AbstractTensorMap{S}, (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}, levels::Tuple;
-          copy::Bool = false) where {S,N₁,N₂}
-        -> tdst::TensorMap{S,N₁,N₂}
+    braid(tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple, levels::IndexTuple;
+          copy::Bool = false)
+        -> tdst::TensorMap
 
 Return tensor `tdst` obtained by braiding the indices of `tsrc`.
 The codomain and domain of `tdst` correspond to the indices in `p₁` and `p₂` of `tsrc` respectively.
@@ -108,18 +106,18 @@ If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwis
 
 To braid into an existing destination, see [braid!](@ref) and [`add_braid!`](@ref)
 """
-function braid(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple, levels::IndexTuple;
-               copy::Bool=false) where {S}
+function braid(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple, levels::IndexTuple;
+               copy::Bool=false)
     @assert length(levels) == numind(t)
-    if BraidingStyle(sectortype(S)) isa SymmetricBraiding
+    if BraidingStyle(sectortype(t)) isa SymmetricBraiding
         return permute(t, (p₁, p₂); copy=copy)
     end
     if !copy && p₁ == codomainind(t) && p₂ == domainind(t)
         return t
     end
     # general case
-    cod = ProductSpace{S}(map(n -> space(t, n), p₁))
-    dom = ProductSpace{S}(map(n -> dual(space(t, n)), p₂))
+    cod = ProductSpace{spacetype(t)}(map(n -> space(t, n), p₁))
+    dom = ProductSpace{spacetype(t)}(map(n -> dual(space(t, n)), p₂))
     @inbounds begin
         return braid!(similar(t, cod ← dom), t, (p₁, p₂), levels)
     end
@@ -130,8 +128,8 @@ end
 _transpose_indices(t::AbstractTensorMap) = (reverse(domainind(t)), reverse(codomainind(t)))
 
 """
-    transpose!(tdst::AbstractTensorMap{S,N₁,N₂}, tsrc::AbstractTensorMap{S},
-               (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}) where {S,N₁,N₂}
+    transpose!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap,
+               (p₁, p₂)::Index2Tuple)
         -> tdst
 
 Write into `tdst` the result of transposing the indices of `tsrc`.
@@ -148,9 +146,9 @@ function LinearAlgebra.transpose!(tdst::AbstractTensorMap,
 end
 
 """
-    transpose(tsrc::AbstractTensorMap{S}, (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}};
-              copy::Bool=false) where {S,N₁,N₂}
-        -> tdst::TensorMap{S,N₁,N₂}
+    transpose(tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple;
+              copy::Bool=false)
+        -> tdst::TensorMap
 
 Return tensor `tdst` obtained by transposing the indices of `tsrc`.
 The codomain and domain of `tdst` correspond to the indices in `p₁` and `p₂` of `tsrc` respectively.
@@ -161,26 +159,26 @@ If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwis
 
 To permute into an existing destination, see [permute!](@ref) and [`add_permute!`](@ref)
 """
-function LinearAlgebra.transpose(t::AbstractTensorMap{S},
+function LinearAlgebra.transpose(t::AbstractTensorMap,
                                  (p₁, p₂)::Index2Tuple=_transpose_indices(t);
-                                 copy::Bool=false) where {S}
-    if sectortype(S) === Trivial
+                                 copy::Bool=false)
+    if sectortype(t) === Trivial
         return permute(t, (p₁, p₂); copy=copy)
     end
     if !copy && p₁ == codomainind(t) && p₂ == domainind(t)
         return t
     end
     # general case
-    cod = ProductSpace{S}(map(n -> space(t, n), p₁))
-    dom = ProductSpace{S}(map(n -> dual(space(t, n)), p₂))
+    cod = ProductSpace{spacetype(t)}(map(n -> space(t, n), p₁))
+    dom = ProductSpace{spacetype(t)}(map(n -> dual(space(t, n)), p₂))
     @inbounds begin
         return transpose!(similar(t, cod ← dom), t, (p₁, p₂))
     end
 end
 
-function LinearAlgebra.transpose(t::AdjointTensorMap{S},
+function LinearAlgebra.transpose(t::AdjointTensorMap,
                                  (p₁, p₂)::Index2Tuple=_transpose_indices(t);
-                                 copy::Bool=false) where {S}
+                                 copy::Bool=false)
     p₁′ = map(n -> adjointtensorindex(t, n), p₂)
     p₂′ = map(n -> adjointtensorindex(t, n), p₁)
     return adjoint(transpose(adjoint(t), (p₁′, p₂′); copy=copy))
@@ -268,23 +266,23 @@ twist(t::AbstractTensorMap, i::Int; inv::Bool=false) = twist!(copy(t), i; inv=in
 #-------------------------------------
 # Full implementations based on `add`
 #-------------------------------------
-@propagate_inbounds function add_permute!(tdst::AbstractTensorMap{S,N₁,N₂},
+@propagate_inbounds function add_permute!(tdst::AbstractTensorMap{E,S,N₁,N₂},
                                           tsrc::AbstractTensorMap,
                                           p::Index2Tuple{N₁,N₂},
                                           α::Number,
                                           β::Number,
-                                          backend::Backend...) where {S,N₁,N₂}
+                                          backend::Backend...) where {E,S,N₁,N₂}
     treepermuter(f₁, f₂) = permute(f₁, f₂, p[1], p[2])
     return add_transform!(tdst, tsrc, p, treepermuter, α, β, backend...)
 end
 
-@propagate_inbounds function add_braid!(tdst::AbstractTensorMap{S,N₁,N₂},
+@propagate_inbounds function add_braid!(tdst::AbstractTensorMap{E,S,N₁,N₂},
                                         tsrc::AbstractTensorMap,
                                         p::Index2Tuple{N₁,N₂},
                                         levels::IndexTuple,
                                         α::Number,
                                         β::Number,
-                                        backend::Backend...) where {S,N₁,N₂}
+                                        backend::Backend...) where {E,S,N₁,N₂}
     length(levels) == numind(tsrc) ||
         throw(ArgumentError("incorrect levels $levels for tensor map $(codomain(tsrc)) ← $(domain(tsrc))"))
 
@@ -295,23 +293,23 @@ end
     return add_transform!(tdst, tsrc, p, treebraider, α, β, backend...)
 end
 
-@propagate_inbounds function add_transpose!(tdst::AbstractTensorMap{S,N₁,N₂},
+@propagate_inbounds function add_transpose!(tdst::AbstractTensorMap{E,S,N₁,N₂},
                                             tsrc::AbstractTensorMap,
                                             p::Index2Tuple{N₁,N₂},
                                             α::Number,
                                             β::Number,
-                                            backend::Backend...) where {S,N₁,N₂}
+                                            backend::Backend...) where {E,S,N₁,N₂}
     treetransposer(f₁, f₂) = transpose(f₁, f₂, p[1], p[2])
     return add_transform!(tdst, tsrc, p, treetransposer, α, β, backend...)
 end
 
-function add_transform!(tdst::AbstractTensorMap{S,N₁,N₂},
+function add_transform!(tdst::AbstractTensorMap{E,S,N₁,N₂},
                         tsrc::AbstractTensorMap,
                         (p₁, p₂)::Index2Tuple{N₁,N₂},
                         fusiontreetransform,
                         α::Number,
                         β::Number,
-                        backend::Backend...) where {S,N₁,N₂}
+                        backend::Backend...) where {E,S,N₁,N₂}
     @boundscheck begin
         all(i -> space(tsrc, p₁[i]) == space(tdst, i), 1:N₁) ||
             throw(SpaceMismatch("source = $(codomain(tsrc))←$(domain(tsrc)),
diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl
index 4678b780..5b79d6e2 100644
--- a/src/tensors/linalg.jl
+++ b/src/tensors/linalg.jl
@@ -164,8 +164,8 @@ function Base.fill!(t::AbstractTensorMap, value::Number)
     end
     return t
 end
-function LinearAlgebra.adjoint!(tdst::AbstractTensorMap{S},
-                                tsrc::AbstractTensorMap{S}) where {S}
+function LinearAlgebra.adjoint!(tdst::AbstractTensorMap,
+                                tsrc::AbstractTensorMap)
     InnerProductStyle(tdst) === EuclideanProduct() || throw_invalid_innerproduct(:adjoint!)
     space(tdst) == adjoint(space(tsrc)) ||
         throw(SpaceMismatch("$(space(tdst)) ≠ adjoint($(space(tsrc)))"))
@@ -383,8 +383,7 @@ for f in (:sqrt, :log, :asin, :acos, :acosh, :atanh, :acoth)
 end
 
 # concatenate tensors
-function catdomain(t1::AbstractTensorMap{S,N₁,1},
-                   t2::AbstractTensorMap{S,N₁,1}) where {S,N₁}
+function catdomain(t1::T, t2::T) where {S,N₁,T<:AbstractTensorMap{<:Any,S,N₁,1}}
     codomain(t1) == codomain(t2) ||
         throw(SpaceMismatch("codomains of tensors to concatenate must match:\n" *
                             "$(codomain(t1)) ≠ $(codomain(t2))"))
@@ -401,8 +400,7 @@ function catdomain(t1::AbstractTensorMap{S,N₁,1},
     end
     return t
 end
-function catcodomain(t1::AbstractTensorMap{S,1,N₂},
-                     t2::AbstractTensorMap{S,1,N₂}) where {S,N₂}
+function catcodomain(t1::T, t2::T) where {S,N₂,T<:AbstractTensorMap{<:Any,S,1,N₂}}
     domain(t1) == domain(t2) ||
         throw(SpaceMismatch("domains of tensors to concatenate must match:\n" *
                             "$(domain(t1)) ≠ $(domain(t2))"))
@@ -422,14 +420,16 @@ end
 
 # tensor product of tensors
 """
-    ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}, ...) -> TensorMap{S}
-    otimes(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}, ...) -> TensorMap{S}
+    ⊗(t1::AbstractTensorMap, t2::AbstractTensorMap, ...) -> TensorMap
+    otimes(t1::AbstractTensorMap, t2::AbstractTensorMap, ...) -> TensorMap
 
 Compute the tensor product between two `AbstractTensorMap` instances, which results in a
 new `TensorMap` instance whose codomain is `codomain(t1) ⊗ codomain(t2)` and whose domain
 is `domain(t1) ⊗ domain(t2)`.
 """
-function ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}) where {S}
+function ⊗(t1::AbstractTensorMap, t2::AbstractTensorMap)
+    (S = spacetype(t1)) === spacetype(t2) ||
+        throw(SpaceMismatch("spacetype(t1) ≠ spacetype(t2)"))
     cod1, cod2 = codomain(t1), codomain(t2)
     dom1, dom2 = domain(t1), domain(t2)
     cod = cod1 ⊗ cod2
diff --git a/src/tensors/tensor.jl b/src/tensors/tensor.jl
index d0538949..80f9d0bb 100644
--- a/src/tensors/tensor.jl
+++ b/src/tensors/tensor.jl
@@ -3,65 +3,67 @@
 #==========================================================#
 #! format: off
 """
-    struct TensorMap{S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{S, N₁, N₂}
+    struct TensorMap{E, S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{E, S, N₁, N₂}
 
 Specific subtype of [`AbstractTensorMap`](@ref) for representing tensor maps (morphisms in
 a tensor category) whose data is stored in blocks of some subtype of `DenseMatrix`.
 """
-struct TensorMap{S<:IndexSpace, N₁, N₂, I<:Sector, A<:Union{<:DenseMatrix,SectorDict{I,<:DenseMatrix}}, F₁, F₂} <: AbstractTensorMap{S, N₁, N₂}
+struct TensorMap{E, S<:IndexSpace, N₁, N₂, I<:Sector, A<:Union{<:DenseMatrix{E},SectorDict{I,<:DenseMatrix{E}}},
+                 F₁, F₂} <: AbstractTensorMap{E, S, N₁, N₂}
     data::A
     codom::ProductSpace{S,N₁}
     dom::ProductSpace{S,N₂}
     rowr::SectorDict{I,FusionTreeDict{F₁,UnitRange{Int}}}
     colr::SectorDict{I,FusionTreeDict{F₂,UnitRange{Int}}}
-    function TensorMap{S, N₁, N₂, I, A, F₁, F₂}(data::A,
+    function TensorMap{E, S, N₁, N₂, I, A, F₁, F₂}(data::A,
                 codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂},
                 rowr::SectorDict{I,FusionTreeDict{F₁,UnitRange{Int}}},
                 colr::SectorDict{I,FusionTreeDict{F₂,UnitRange{Int}}}) where
-                    {S<:IndexSpace, N₁, N₂, I<:Sector, A<:SectorDict{I,<:DenseMatrix},
+                    {E,S<:IndexSpace, N₁, N₂, I<:Sector, A<:SectorDict{I,<:DenseMatrix{E}},
                      F₁<:FusionTree{I,N₁}, F₂<:FusionTree{I,N₂}}
-        T = scalartype(valtype(data))
-        T ⊆ field(S) || @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1)
-        return new{S,N₁,N₂,I,A,F₁,F₂}(data, codom, dom, rowr, colr)
+        E ⊆ field(S) || @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1)
+        return new{E,S,N₁,N₂,I,A,F₁,F₂}(data, codom, dom, rowr, colr)
     end
-    function TensorMap{S,N₁,N₂,Trivial,A,Nothing,Nothing}(data::A,
+    function TensorMap{E,S,N₁,N₂,Trivial,A,Nothing,Nothing}(data::A,
                                                           codom::ProductSpace{S,N₁},
                                                           dom::ProductSpace{S,N₂}) where
-             {S<:IndexSpace,N₁,N₂,A<:DenseMatrix}
-        T = scalartype(data)
-        T ⊆ field(S) ||
-            @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1)
-        return new{S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom)
+             {E,S<:IndexSpace,N₁,N₂,A<:DenseMatrix{E}}
+        E ⊆ field(S) || @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1)
+        return new{E,S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom)
     end
 end
 #! format: on
 
 """
-    Tensor{S, N, I, A, F₁, F₂} = TensorMap{S, N, 0, I, A, F₁, F₂}
+    Tensor{E, S, N, I, A, F₁, F₂} = TensorMap{E, S, N, 0, I, A, F₁, F₂}
 
 Specific subtype of [`AbstractTensor`](@ref) for representing tensors whose data is stored
 in blocks of some subtype of `DenseMatrix`.
 
-A `Tensor{S, N, I, A, F₁, F₂}` is actually a special case `TensorMap{S, N, 0, I, A, F₁, F₂}`,
+A `Tensor{E, S, N, I, A, F₁, F₂}` is actually a special case `TensorMap{E, S, N, 0, I, A, F₁, F₂}`,
 i.e. a tensor map with only a non-trivial output space.
 """
-const Tensor{S,N,I,A,F₁,F₂} = TensorMap{S,N,0,I,A,F₁,F₂}
+const Tensor{E,S,N,I,A,F₁,F₂} = TensorMap{E,S,N,0,I,A,F₁,F₂}
+
 """
-    TrivialTensorMap{S<:IndexSpace, N₁, N₂, A<:DenseMatrix} = TensorMap{S, N₁, N₂, Trivial, 
+    TrivialTensorMap{E, S, N₁, N₂, A<:DenseMatrix} = TensorMap{E, S, N₁, N₂, Trivial, 
                                                                         A, Nothing, Nothing}
 
 A special case of [`TensorMap`](@ref) for representing tensor maps with trivial symmetry,
 i.e., whose `sectortype` is `Trivial`.
 """
-const TrivialTensorMap{S,N₁,N₂,A<:DenseMatrix} = TensorMap{S,N₁,N₂,Trivial,A,
-                                                           Nothing,Nothing}
+const TrivialTensorMap{E,S,N₁,N₂,A<:DenseMatrix} = TensorMap{E,S,N₁,N₂,Trivial,A,Nothing,
+                                                             Nothing}
+
 """
-    TrivialTensor{S, N, A} = TrivialTensorMap{S, N, 0, A}
+    TrivialTensor{E, S, N, A} = TrivialTensorMap{E, S, N, 0, A}
 
 A special case of [`Tensor`](@ref) for representing tensors with trivial symmetry, i.e.,
 whose `sectortype` is `Trivial`.
 """
-const TrivialTensor{S,N,A} = TrivialTensorMap{S,N,0,A}
+const TrivialTensor{E,S,N,A} = TrivialTensorMap{E,S,N,0,A}
+
+# TODO: check if argument order should change
 """
     tensormaptype(::Type{S}, N₁::Int, N₂::Int, [::Type{T}]) where {S<:IndexSpace,T} -> ::Type{<:TensorMap}
 
@@ -72,17 +74,19 @@ function tensormaptype(::Type{S}, N₁::Int, N₂::Int, ::Type{T}) where {S,T}
     I = sectortype(S)
     if T <: DenseMatrix
         M = T
+        E = scalartype(T)
     elseif T <: Number
         M = Matrix{T}
+        E = T
     else
         throw(ArgumentError("the final argument of `tensormaptype` should either be the scalar or the storage type, i.e. a subtype of `Number` or of `DenseMatrix`"))
     end
     if I === Trivial
-        return TensorMap{S,N₁,N₂,I,M,Nothing,Nothing}
+        return TensorMap{E,S,N₁,N₂,I,M,Nothing,Nothing}
     else
         F₁ = fusiontreetype(I, N₁)
         F₂ = fusiontreetype(I, N₂)
-        return TensorMap{S,N₁,N₂,I,SectorDict{I,M},F₁,F₂}
+        return TensorMap{E,S,N₁,N₂,I,SectorDict{I,M},F₁,F₂}
     end
 end
 tensormaptype(S, N₁, N₂=0) = tensormaptype(S, N₁, N₂, Float64)
@@ -100,12 +104,12 @@ blocksectors(t::TensorMap) = keys(t.data)
 
 Return the type of the storage `A` of the tensor map.
 """
-function storagetype(::Type{<:TensorMap{<:IndexSpace,N₁,N₂,Trivial,A}}) where
-         {N₁,N₂,A<:DenseMatrix}
+function storagetype(::Type{<:TrivialTensorMap{E,S,N₁,N₂,A}}) where
+         {E,S,N₁,N₂,A}
     return A
 end
-function storagetype(::Type{<:TensorMap{<:IndexSpace,N₁,N₂,I,<:SectorDict{I,A}}}) where
-         {N₁,N₂,I<:Sector,A<:DenseMatrix}
+function storagetype(::Type{<:TensorMap{E,S,N₁,N₂,I,<:SectorDict{I,A}}}) where
+         {E,S,N₁,N₂,I<:Sector,A<:DenseMatrix}
     return A
 end
 
@@ -159,14 +163,15 @@ function TensorMap(data::AbstractDict{<:Sector,<:DenseMatrix}, codom::ProductSpa
     end
     F₁ = fusiontreetype(I, N₁)
     F₂ = fusiontreetype(I, N₂)
+    E = scalartype(valtype(data))
     if !isreal(I)
         data2 = SectorDict(c => complex(data[c]) for c in blocksectoriterator)
         A = typeof(data2)
-        return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data2, codom, dom, rowr, colr)
+        return TensorMap{E,S,N₁,N₂,I,A,F₁,F₂}(data2, codom, dom, rowr, colr)
     else
         data2 = SectorDict(c => data[c] for c in blocksectoriterator)
         A = typeof(data2)
-        return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data2, codom, dom, rowr, colr)
+        return TensorMap{E,S,N₁,N₂,I,A,F₁,F₂}(data2, codom, dom, rowr, colr)
     end
 end
 
@@ -179,7 +184,8 @@ function TensorMap(f, codom::ProductSpace{S,N₁},
         d2 = dim(dom)
         data = f((d1, d2))
         A = typeof(data)
-        return TensorMap{S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom)
+        E = scalartype(A)
+        return TensorMap{E,S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom)
     end
     blocksectoriterator = blocksectors(codom ← dom)
     rowr, rowdims = _buildblockstructure(codom, blocksectoriterator)
@@ -193,7 +199,8 @@ function TensorMap(f, codom::ProductSpace{S,N₁},
     F₁ = fusiontreetype(I, N₁)
     F₂ = fusiontreetype(I, N₂)
     A = typeof(data)
-    return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data, codom, dom, rowr, colr)
+    E = scalartype(valtype(A))
+    return TensorMap{E,S,N₁,N₂,I,A,F₁,F₂}(data, codom, dom, rowr, colr)
 end
 
 # auxiliary function
@@ -353,7 +360,8 @@ function TensorMap(data::DenseArray, codom::ProductSpace{S,N₁}, dom::ProductSp
     if sectortype(S) === Trivial
         data2 = reshape(data, (d1, d2))
         A = typeof(data2)
-        return TensorMap{S,N₁,N₂,Trivial,A,Nothing,Nothing}(data2, codom, dom)
+        E = scalartype(A)
+        return TensorMap{E,S,N₁,N₂,Trivial,A,Nothing,Nothing}(data2, codom, dom)
     else
         t = TensorMap(zeros, eltype(data), codom, dom)
         ta = convert(Array, t)
@@ -425,7 +433,7 @@ Base.similar(t::AbstractTensorMap, T::Type) = similar(t, T, space(t))
 Base.similar(t::AbstractTensorMap) = similar(t, scalartype(t), space(t))
 
 # actual implementation
-function Base.similar(t::TensorMap{S}, ::Type{T}, P::TensorMapSpace{S}) where {T,S}
+function Base.similar(t::TensorMap, ::Type{T}, P::TensorMapSpace{S}) where {T,S}
     N₁ = length(codomain(P))
     N₂ = length(domain(P))
     I = sectortype(S)
@@ -433,14 +441,14 @@ function Base.similar(t::TensorMap{S}, ::Type{T}, P::TensorMapSpace{S}) where {T
     if I === Trivial
         data = similar(t.data, T, (dim(codomain(P)), dim(domain(P))))
         A = typeof(data)
-        return TrivialTensorMap{S,N₁,N₂,A}(data, codomain(P), domain(P))
+        return TrivialTensorMap{T,S,N₁,N₂,A}(data, codomain(P), domain(P))
     end
     F₁ = fusiontreetype(I, N₁)
     F₂ = fusiontreetype(I, N₂)
     if space(t) == P
         data = SectorDict(c => similar(b, T) for (c, b) in blocks(t))
         A = typeof(data)
-        return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data, codomain(P), domain(P), t.rowr, t.colr)
+        return TensorMap{T,S,N₁,N₂,I,A,F₁,F₂}(data, codomain(P), domain(P), t.rowr, t.colr)
     end
 
     blocksectoriterator = blocksectors(P)
@@ -484,7 +492,7 @@ function Base.similar(t::TensorMap{S}, ::Type{T}, P::TensorMapSpace{S}) where {T
     data = SectorDict{I,M}(c => M(undef, (rowdims[c], coldims[c]))
                            for c in blocksectoriterator)
     A = typeof(data)
-    return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data, codomain(P), domain(P), rowr, colr)
+    return TensorMap{T,S,N₁,N₂,I,A,F₁,F₂}(data, codomain(P), domain(P), rowr, colr)
 end
 
 function Base.complex(t::AbstractTensorMap)
@@ -538,7 +546,7 @@ function block(t::TensorMap, s::Sector)
     end
 end
 
-function blocks(t::TensorMap{<:IndexSpace,N₁,N₂,Trivial}) where {N₁,N₂}
+function blocks(t::TrivialTensorMap)
     return SingletonDict(Trivial() => t.data)
 end
 blocks(t::TensorMap) = t.data
@@ -547,7 +555,7 @@ fusiontrees(t::TrivialTensorMap) = ((nothing, nothing),)
 fusiontrees(t::TensorMap) = TensorKeyIterator(t.rowr, t.colr)
 
 """
-    Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I},
+    Base.getindex(t::TensorMap
                   sectors::NTuple{N₁+N₂,I}) where {N₁,N₂,I<:Sector} 
         -> StridedViews.StridedView
     t[sectors]
@@ -561,8 +569,9 @@ respectively, then a `StridedViews.StridedView` of size
 This method is only available for the case where `FusionStyle(I) isa UniqueFusion`,
 since it assumes a  uniquely defined coupled charge.
 """
-@inline function Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I},
-                               sectors::Tuple{Vararg{I}}) where {N₁,N₂,I<:Sector}
+@inline function Base.getindex(t::TensorMap, sectors::Tuple{I,Vararg{I}}) where {I<:Sector}
+    I === sectortype(t) || throw(SectorMismatch("Not a valid sectortype for this tensor."))
+    length(sectors) == numind(t) || throw(ArgumentError("Number of sectors does not match."))
     FusionStyle(I) isa UniqueFusion ||
         throw(SectorMismatch("Indexing with sectors only possible if unique fusion"))
     s1 = TupleTools.getindices(sectors, codomainind(t))
@@ -584,9 +593,9 @@ end
 end
 
 """
-    Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I},
+    Base.getindex(t::TensorMap{E,S,N₁,N₂,I},
                   f₁::FusionTree{I,N₁},
-                  f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector}
+                  f₂::FusionTree{I,N₂}) where {E,SN₁,N₂,I<:Sector}
         -> StridedViews.StridedView
     t[f₁, f₂]
 
@@ -597,9 +606,9 @@ Return a view into the data slice of `t` corresponding to the splitting - fusion
 represents the slice of `block(t, c)` whose row indices correspond to `f₁.uncoupled` and
 column indices correspond to `f₂.uncoupled`.
 """
-@inline function Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I},
+@inline function Base.getindex(t::TensorMap{E,S,N₁,N₂,I},
                                f₁::FusionTree{I,N₁},
-                               f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector}
+                               f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂,I<:Sector}
     c = f₁.coupled
     @boundscheck begin
         c == f₂.coupled || throw(SectorMismatch())
@@ -613,10 +622,10 @@ column indices correspond to `f₂.uncoupled`.
 end
 
 """
-    Base.setindex!(t::TensorMap{<:IndexSpace,N₁,N₂,I},
+    Base.setindex!(t::TensorMap{E,S,N₁,N₂,I},
                    v,
                    f₁::FusionTree{I,N₁},
-                   f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector}
+                   f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂,I<:Sector}
     t[f₁, f₂] = v
 
 Copies `v` into the  data slice of `t` corresponding to the splitting - fusion tree pair
@@ -624,12 +633,13 @@ Copies `v` into the  data slice of `t` corresponding to the splitting - fusion t
 of size `(dims(codomain(t), f₁.uncoupled)..., dims(domain(t), f₂.uncoupled))` using
 `Base.copy!`.
 
-See also [`Base.getindex(::TensorMap{<:IndexSpace,N₁,N₂,I<:Sector}, ::FusionTree{I<:Sector,N₁}, ::FusionTree{I<:Sector,N₂})`](@ref)
+See also [`Base.getindex(::TensorMap{E,S,N₁,N₂,I<:Sector}, ::FusionTree{I<:Sector,N₁}, ::FusionTree{I<:Sector,N₂})`](@ref)
 """
-@propagate_inbounds function Base.setindex!(t::TensorMap{<:IndexSpace,N₁,N₂,I},
+@propagate_inbounds function Base.setindex!(t::TensorMap{E,S,N₁,N₂,I},
                                             v,
                                             f₁::FusionTree{I,N₁},
-                                            f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector}
+                                            f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂,
+                                                                         I<:Sector}
     return copy!(getindex(t, f₁, f₂), v)
 end
 
@@ -682,16 +692,16 @@ end
 function Base.summary(t::TensorMap)
     return print("TensorMap(", space(t), ")")
 end
-function Base.show(io::IO, t::TensorMap{S}) where {S<:IndexSpace}
+function Base.show(io::IO, t::TensorMap)
     if get(io, :compact, false)
         print(io, "TensorMap(", space(t), ")")
         return
     end
     println(io, "TensorMap(", space(t), "):")
-    if sectortype(S) == Trivial
+    if sectortype(t) == Trivial
         Base.print_array(io, t[])
         println(io)
-    elseif FusionStyle(sectortype(S)) isa UniqueFusion
+    elseif FusionStyle(sectortype(t)) isa UniqueFusion
         for (f₁, f₂) in fusiontrees(t)
             println(io, "* Data for sector ", f₁.uncoupled, " ← ", f₂.uncoupled, ":")
             Base.print_array(io, t[f₁, f₂])
@@ -708,28 +718,28 @@ end
 
 # Real and imaginary parts
 #---------------------------
-function Base.real(t::AbstractTensorMap{S}) where {S}
+function Base.real(t::AbstractTensorMap)
     # `isreal` for a `Sector` returns true iff the F and R symbols are real. This guarantees
     # that the real/imaginary part of a tensor `t` can be obtained by just taking
     # real/imaginary part of the degeneracy data.
-    if isreal(sectortype(S))
+    if isreal(sectortype(t))
         realdata = Dict(k => real(v) for (k, v) in blocks(t))
         return TensorMap(realdata, codomain(t), domain(t))
     else
-        msg = "`real` has not been implemented for `AbstractTensorMap{$(S)}`."
+        msg = "`real` has not been implemented for `$(typeof(t))`."
         throw(ArgumentError(msg))
     end
 end
 
-function Base.imag(t::AbstractTensorMap{S}) where {S}
+function Base.imag(t::AbstractTensorMap)
     # `isreal` for a `Sector` returns true iff the F and R symbols are real. This guarantees
     # that the real/imaginary part of a tensor `t` can be obtained by just taking
     # real/imaginary part of the degeneracy data.
-    if isreal(sectortype(S))
+    if isreal(sectortype(t))
         imagdata = Dict(k => imag(v) for (k, v) in blocks(t))
         return TensorMap(imagdata, codomain(t), domain(t))
     else
-        msg = "`imag` has not been implemented for `AbstractTensorMap{$(S)}`."
+        msg = "`imag` has not been implemented for `$(typeof(t))`."
         throw(ArgumentError(msg))
     end
 end
@@ -741,19 +751,20 @@ function Base.convert(::Type{TensorMap}, t::AbstractTensorMap)
     return copy!(TensorMap(undef, scalartype(t), codomain(t), domain(t)), t)
 end
 
-function Base.convert(T::Type{TensorMap{S,N₁,N₂,I,A,F₁,F₂}},
-                      t::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂,I,A,F₁,F₂}
-    if typeof(t) == T
+function Base.convert(T::Type{<:TensorMap{E,S,N₁,N₂}},
+                      t::AbstractTensorMap{<:Any,S,N₁,N₂}) where {E,S,N₁,N₂}
+    if typeof(t) === T
         return t
     else
-        data = Dict{I,storagetype(T)}(c => convert(storagetype(T), b)
+        data = Dict{sectortype(T),storagetype(T)}(c => convert(storagetype(T), b)
                                       for (c, b) in blocks(t))
         return TensorMap(data, codomain(t), domain(t))
     end
 end
 
-function Base.promote_rule(::Type{<:T1},
-                           t2::Type{<:T2}) where {S,N₁,N₂,T1<:TensorMap{S,N₁,N₂},
-                                                  T2<:TensorMap{S,N₁,N₂}}
-    return tensormaptype(S, N₁, N₂, promote_type(storagetype(T1), storagetype(T2)))
+function Base.promote_rule(::Type{<:T₁},
+                           ::Type{<:T₂}) where {S,N₁,N₂,
+                                                T₁<:TensorMap{<:Any,S,N₁,N₂},
+                                                T₂<:TensorMap{<:Any,S,N₁,N₂}}
+    return tensormaptype(S, N₁, N₂, promote_type(storagetype(T₁), storagetype(T₂)))
 end
diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl
index 9f8ed5a9..9114f767 100644
--- a/src/tensors/tensoroperations.jl
+++ b/src/tensors/tensoroperations.jl
@@ -34,9 +34,9 @@ function _canonicalize(p::IndexTuple, t::AbstractTensorMap)
 end
 
 # tensoradd!
-function TO.tensoradd!(C::AbstractTensorMap{S}, pC::Index2Tuple,
-                       A::AbstractTensorMap{S}, conjA::Symbol,
-                       α::Number, β::Number, backend::Backend...) where {S}
+function TO.tensoradd!(C::AbstractTensorMap, pC::Index2Tuple,
+                       A::AbstractTensorMap, conjA::Symbol,
+                       α::Number, β::Number, backend::Backend...)
     if conjA == :N
         A′ = A
         pC′ = _canonicalize(pC, C)
@@ -50,17 +50,17 @@ function TO.tensoradd!(C::AbstractTensorMap{S}, pC::Index2Tuple,
     return C
 end
 
-function TO.tensoradd_type(TC, ::Index2Tuple{N₁,N₂}, A::AbstractTensorMap{S},
-                           ::Symbol) where {S,N₁,N₂}
+function TO.tensoradd_type(TC, ::Index2Tuple{N₁,N₂}, A::AbstractTensorMap,
+                           ::Symbol) where {N₁,N₂}
     M = similarstoragetype(A, TC)
-    return tensormaptype(S, N₁, N₂, M)
+    return tensormaptype(spacetype(A), N₁, N₂, M)
 end
 
 function TO.tensoradd_structure(pC::Index2Tuple{N₁,N₂},
-                                A::AbstractTensorMap{S}, conjA::Symbol) where {S,N₁,N₂}
+                                A::AbstractTensorMap, conjA::Symbol) where {N₁,N₂}
     if conjA == :N
-        cod = ProductSpace{S,N₁}(space.(Ref(A), pC[1]))
-        dom = ProductSpace{S,N₂}(dual.(space.(Ref(A), pC[2])))
+        cod = ProductSpace{spacetype(A),N₁}(space.(Ref(A), pC[1]))
+        dom = ProductSpace{spacetype(A),N₂}(dual.(space.(Ref(A), pC[2])))
         return dom → cod
     else
         return TO.tensoradd_structure(adjointtensorindices(A, pC), adjoint(A), :N)
@@ -68,9 +68,9 @@ function TO.tensoradd_structure(pC::Index2Tuple{N₁,N₂},
 end
 
 # tensortrace!
-function TO.tensortrace!(C::AbstractTensorMap{S}, p::Index2Tuple,
-                         A::AbstractTensorMap{S}, q::Index2Tuple, conjA::Symbol,
-                         α::Number, β::Number, backend::Backend...) where {S}
+function TO.tensortrace!(C::AbstractTensorMap, p::Index2Tuple,
+                         A::AbstractTensorMap, q::Index2Tuple, conjA::Symbol,
+                         α::Number, β::Number, backend::Backend...)
     if conjA == :N
         A′ = A
         p′ = _canonicalize(p, C)
@@ -89,10 +89,10 @@ function TO.tensortrace!(C::AbstractTensorMap{S}, p::Index2Tuple,
 end
 
 # tensorcontract!
-function TO.tensorcontract!(C::AbstractTensorMap{S,N₁,N₂}, pAB::Index2Tuple,
-                            A::AbstractTensorMap{S}, pA::Index2Tuple, conjA::Symbol,
-                            B::AbstractTensorMap{S}, pB::Index2Tuple, conjB::Symbol,
-                            α::Number, β::Number, backend::Backend...) where {S,N₁,N₂}
+function TO.tensorcontract!(C::AbstractTensorMap, pAB::Index2Tuple,
+                            A::AbstractTensorMap, pA::Index2Tuple, conjA::Symbol,
+                            B::AbstractTensorMap, pB::Index2Tuple, conjB::Symbol,
+                            α::Number, β::Number, backend::Backend...)
     pAB′ = _canonicalize(pAB, C)
     if conjA == :N
         A′ = A
@@ -117,28 +117,30 @@ function TO.tensorcontract!(C::AbstractTensorMap{S,N₁,N₂}, pAB::Index2Tuple,
 end
 
 function TO.tensorcontract_type(TC, ::Index2Tuple{N₁,N₂},
-                                A::AbstractTensorMap{S}, pA, conjA,
-                                B::AbstractTensorMap{S}, pB, conjB) where {S,N₁,N₂}
+                                A::AbstractTensorMap, pA, conjA,
+                                B::AbstractTensorMap, pB, conjB) where {N₁,N₂}
     M = similarstoragetype(A, TC)
     M == similarstoragetype(B, TC) || throw(ArgumentError("incompatible storage types"))
-    return tensormaptype(S, N₁, N₂, M)
+    spacetype(A) == spacetype(B) || throw(SpaceMismatch("incompatible space types"))
+    return tensormaptype(spacetype(A), N₁, N₂, M)
 end
 
 function TO.tensorcontract_structure(pC::Index2Tuple{N₁,N₂},
-                                     A::AbstractTensorMap{S}, pA::Index2Tuple, conjA,
-                                     B::AbstractTensorMap{S}, pB::Index2Tuple,
-                                     conjB) where {S,N₁,N₂}
+                                     A::AbstractTensorMap, pA::Index2Tuple, conjA,
+                                     B::AbstractTensorMap, pB::Index2Tuple,
+                                     conjB) where {N₁,N₂}
+    spacetype(A) == spacetype(B) || throw(SpaceMismatch("incompatible space types"))
     spaces1 = TO.flag2op(conjA).(space.(Ref(A), pA[1]))
     spaces2 = TO.flag2op(conjB).(space.(Ref(B), pB[2]))
     spaces = (spaces1..., spaces2...)
-    cod = ProductSpace{S,N₁}(getindex.(Ref(spaces), pC[1]))
-    dom = ProductSpace{S,N₂}(dual.(getindex.(Ref(spaces), pC[2])))
+    cod = ProductSpace{spacetype(A),N₁}(getindex.(Ref(spaces), pC[1]))
+    dom = ProductSpace{spacetype(A),N₂}(dual.(getindex.(Ref(spaces), pC[2])))
     return dom → cod
 end
 
-function TO.checkcontractible(tA::AbstractTensorMap{S}, iA::Int, conjA::Symbol,
-                              tB::AbstractTensorMap{S}, iB::Int, conjB::Symbol,
-                              label) where {S}
+function TO.checkcontractible(tA::AbstractTensorMap, iA::Int, conjA::Symbol,
+                              tB::AbstractTensorMap, iB::Int, conjB::Symbol,
+                              label)
     sA = TO.tensorstructure(tA, iA, conjA)'
     sB = TO.tensorstructure(tB, iB, conjB)
     sA == sB ||
@@ -154,17 +156,27 @@ TO.tensorcost(t::AbstractTensorMap, i::Int) = dim(space(t, i))
 
 # Trace implementation
 #----------------------
-function trace_permute!(tdst::AbstractTensorMap{S,N₁,N₂},
-                        tsrc::AbstractTensorMap{S},
-                        (p₁, p₂)::Index2Tuple{N₁,N₂},
-                        (q₁, q₂)::Index2Tuple{N₃,N₃},
+function trace_permute!(tdst::AbstractTensorMap,
+                        tsrc::AbstractTensorMap,
+                        (p₁, p₂)::Index2Tuple,
+                        (q₁, q₂)::Index2Tuple,
                         α::Number,
                         β::Number,
-                        backend::Backend...) where {S,N₁,N₂,N₃}
+                        backend::Backend...)
+    # some input checks
+    (S = spacetype(tdst)) == spacetype(tsrc) ||
+        throw(SpaceMismatch("incompatible spacetypes"))
     if !(BraidingStyle(sectortype(S)) isa SymmetricBraiding)
         throw(SectorMismatch("only tensors with symmetric braiding rules can be contracted; try `@planar` instead"))
     end
+    (N₃ = length(q₁)) == length(q₂) ||
+        throw(IndexError("number of trace indices does not match"))
+
+    N₁, N₂ = length(p₁), length(p₂)
+
     @boundscheck begin
+        numout(tdst) == N₁ || throw(IndexError("number of output indices does not match"))
+        numin(tdst) == N₂ || throw(IndexError("number of input indices does not match"))
         all(i -> space(tsrc, p₁[i]) == space(tdst, i), 1:N₁) ||
             throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)),
                     tdst = $(codomain(tdst))←$(domain(tdst)), p₁ = $(p₁), p₂ = $(p₂)"))
@@ -217,15 +229,15 @@ end
 # TODO: contraction with either A or B a rank (1, 1) tensor does not require to
 # permute the fusion tree and should therefore be special cased. This will speed
 # up MPS algorithms
-function contract!(C::AbstractTensorMap{S},
-                   A::AbstractTensorMap{S},
-                   (oindA, cindA)::Index2Tuple{N₁,N₃},
-                   B::AbstractTensorMap{S},
-                   (cindB, oindB)::Index2Tuple{N₃,N₂},
+function contract!(C::AbstractTensorMap,
+                   A::AbstractTensorMap, (oindA, cindA)::Index2Tuple,
+                   B::AbstractTensorMap, (cindB, oindB)::Index2Tuple,
                    (p₁, p₂)::Index2Tuple,
-                   α::Number,
-                   β::Number,
-                   backend::Backend...) where {S,N₁,N₂,N₃}
+                   α::Number, β::Number,
+                   backend::Backend...)
+    length(cindA) == length(cindB) ||
+        throw(IndexError("number of contracted indices does not match"))
+    N₁, N₂ = length(oindA), length(oindB)
 
     # find optimal contraction scheme
     hsp = has_shared_permute
@@ -275,16 +287,17 @@ function contract!(C::AbstractTensorMap{S},
 end
 
 # TODO: also transform _contract! into new interface, and add backend support
-function _contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S},
-                    β, C::AbstractTensorMap{S},
-                    oindA::IndexTuple{N₁}, cindA::IndexTuple,
-                    oindB::IndexTuple{N₂}, cindB::IndexTuple,
-                    p₁::IndexTuple, p₂::IndexTuple) where {S,N₁,N₂}
-    if !(BraidingStyle(sectortype(S)) isa SymmetricBraiding)
+function _contract!(α, A::AbstractTensorMap, B::AbstractTensorMap,
+                    β, C::AbstractTensorMap,
+                    oindA::IndexTuple, cindA::IndexTuple,
+                    oindB::IndexTuple, cindB::IndexTuple,
+                    p₁::IndexTuple, p₂::IndexTuple)
+    if !(BraidingStyle(sectortype(C)) isa SymmetricBraiding)
         throw(SectorMismatch("only tensors with symmetric braiding rules can be contracted; try `@planar` instead"))
     end
+    N₁, N₂ = length(oindA), length(oindB)
     copyA = false
-    if BraidingStyle(sectortype(S)) isa Fermionic
+    if BraidingStyle(sectortype(A)) isa Fermionic
         for i in cindA
             if !isdual(space(A, i))
                 copyA = true
@@ -293,7 +306,7 @@ function _contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S},
     end
     A′ = permute(A, (oindA, cindA); copy=copyA)
     B′ = permute(B, (cindB, oindB))
-    if BraidingStyle(sectortype(S)) isa Fermionic
+    if BraidingStyle(sectortype(A)) isa Fermionic
         for i in domainind(A′)
             if !isdual(space(A′, i))
                 A′ = twist!(A′, i)
@@ -315,7 +328,7 @@ end
 
 # Scalar implementation
 #-----------------------
-function scalar(t::AbstractTensorMap{S}) where {S<:IndexSpace}
+function scalar(t::AbstractTensorMap)
     return dim(codomain(t)) == dim(domain(t)) == 1 ?
            first(blocks(t))[2][1, 1] : throw(DimensionMismatch())
 end
diff --git a/test/planar.jl b/test/planar.jl
index d7edf146..c9b8cc54 100644
--- a/test/planar.jl
+++ b/test/planar.jl
@@ -12,13 +12,13 @@ function force_planar(V::GradedSpace)
     return GradedSpace((c ⊠ PlanarTrivial() => dim(V, c) for c in sectors(V))..., isdual(V))
 end
 force_planar(V::ProductSpace) = mapreduce(force_planar, ⊗, V)
-function force_planar(tsrc::TensorMap{ComplexSpace})
+function force_planar(tsrc::TensorMap{<:Any,ComplexSpace})
     tdst = TensorMap(undef, scalartype(tsrc),
                      force_planar(codomain(tsrc)) ← force_planar(domain(tsrc)))
     copyto!(blocks(tdst)[PlanarTrivial()], blocks(tsrc)[Trivial()])
     return tdst
 end
-function force_planar(tsrc::TensorMap{<:GradedSpace})
+function force_planar(tsrc::TensorMap{<:Any,<:GradedSpace})
     tdst = TensorMap(undef, scalartype(tsrc),
                      force_planar(codomain(tsrc)) ← force_planar(domain(tsrc)))
     for (c, b) in blocks(tsrc)

From 8cf81940417cce04fb7f198bf872cff3bd1c2028 Mon Sep 17 00:00:00 2001
From: Jutho <Jutho@users.noreply.github.com>
Date: Tue, 9 Apr 2024 12:48:05 +0200
Subject: [PATCH 04/21] Update Changelog.md

---
 Changelog.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Changelog.md b/Changelog.md
index 12a5874a..7a9f8524 100644
--- a/Changelog.md
+++ b/Changelog.md
@@ -5,6 +5,7 @@ Features that are planned to be implemented before the release of v1.0.0, in no
 - [ ] Separate `Sectors` module
 - [ ] Make `TrivialTensorMap` and `TensorMap` be the same
 - [ ] Simplify `TensorMap` type to hide `rowr` and `colr`
+- [ ] Change block order in `rowr` / `colr` to speed up particular contractions
 - [ ] Make `AdjointTensorMap` generic
 - [ ] Rewrite planar operations in order to be AD-compatible
 - [ ] Fix rrules for fermionic tensors

From f93eb7a6048d3d050afe3a4e8f1f03540667fb0a Mon Sep 17 00:00:00 2001
From: Lukas <37111893+lkdvos@users.noreply.github.com>
Date: Sun, 21 Apr 2024 17:49:16 +0200
Subject: [PATCH 05/21] Change `copy(BraidingTensor)` to return a
 `BraidingTensor` (#2)

This changes the behaviour of copy as an instantiator for creating a TensorMap from a BraidingTensor.
The rationale is that while this does sometimes happen in Julia Base, this is always in the context of lazy wrapper types, for which it makes sense to not copy the parent and then create a new wrapper. BraidingTensor does not wrap anything, so this definition makes less sense.
---
 Changelog.md                  |  8 ++++++++
 src/tensors/braidingtensor.jl | 23 +++++++++++++----------
 src/tensors/tensor.jl         |  5 +++--
 test/braidingtensor.jl        | 34 +++++++++++++++++-----------------
 4 files changed, 41 insertions(+), 29 deletions(-)

diff --git a/Changelog.md b/Changelog.md
index 7a9f8524..be534eed 100644
--- a/Changelog.md
+++ b/Changelog.md
@@ -21,3 +21,11 @@ This adds the scalar type as a parameter to the `AbstractTensorMap` type. This i
 the contexts where different types of tensors are used (`AdjointTensor`, `BraidingTensor`,
 ...), which still have the same scalartype. Additionally, this removes many specializations
 for methods in order to reduce ambiguity errors.
+
+### `copy(BraidingTensor`
+
+This PR changes the behaviour of `copy` as an instantiator for creating a `TensorMap` from a
+`BraidingTensor`. The rationale is that while this does sometimes happen in Julia `Base`,
+this is always in the context of lazy wrapper types, for which it makes sense to not copy
+the parent and then create a new wrapper. `BraidingTensor` does not wrap anything, so this
+definition makes less sense.
diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl
index f5924bf0..98a2e08e 100644
--- a/src/tensors/braidingtensor.jl
+++ b/src/tensors/braidingtensor.jl
@@ -127,7 +127,9 @@ end
 
 Base.similar(::BraidingTensor, T::Type, P::TensorMapSpace) = TensorMap(undef, T, P)
 
-Base.copy(b::BraidingTensor) = copy!(similar(b), b)
+# efficient copy constructor
+Base.copy(b::BraidingTensor) = BraidingTensor(copy(b.V1), copy(b.V2), b.adjoint)
+
 function Base.copy!(t::TensorMap, b::BraidingTensor)
     space(t) == space(b) || throw(SectorMismatch())
     fill!(t, zero(scalartype(t)))
@@ -148,8 +150,8 @@ function Base.copy!(t::TensorMap, b::BraidingTensor)
     end
     return t
 end
-TensorMap(b::BraidingTensor) = copy(b)
-Base.convert(::Type{TensorMap}, b::BraidingTensor) = copy(b)
+TensorMap(b::BraidingTensor) = copy!(similar(b), b)
+Base.convert(::Type{TensorMap}, b::BraidingTensor) = TensorMap(b)
 
 function block(b::BraidingTensor, s::Sector)
     sectortype(b) == typeof(s) || throw(SectorMismatch())
@@ -199,7 +201,8 @@ function add_transform!(tdst::AbstractTensorMap,
                         α::Number,
                         β::Number,
                         backend::Backend...)
-    return add_transform!(tdst, copy(tsrc), (p₁, p₂), fusiontreetransform, α, β, backend...)
+    return add_transform!(tdst, TensorMap(tsrc), (p₁, p₂), fusiontreetransform, α, β,
+                          backend...)
 end
 
 # VectorInterface
@@ -213,7 +216,7 @@ end
 function TO.tensoradd!(C::AbstractTensorMap, pC::Index2Tuple,
                        A::BraidingTensor, conjA::Symbol, α::Number, β::Number,
                        backend::Backend...)
-    return TO.tensoradd!(C, pC, copy(A), conjA, α, β, backend...)
+    return TO.tensoradd!(C, pC, TensorMap(A), conjA, α, β, backend...)
 end
 
 # Planar operations
@@ -224,7 +227,7 @@ function planaradd!(C::AbstractTensorMap,
                     A::BraidingTensor, p::Index2Tuple,
                     α::Number, β::Number,
                     backend::Backend...)
-    return planaradd!(C, copy(A), p, α, β, backend...)
+    return planaradd!(C, TensorMap(A), p, α, β, backend...)
 end
 
 function planarcontract!(C::AbstractTensorMap,
@@ -237,7 +240,7 @@ function planarcontract!(C::AbstractTensorMap,
                          backend::Backend...)
     # special case only defined for contracting 2 indices
     length(oindA) == length(cindA) == 2 ||
-        return planarcontract!(C, copy(A), (oindA, cindA), B, (cindB, oindB), (p1, p2),
+        return planarcontract!(C, TensorMap(A), (oindA, cindA), B, (cindB, oindB), (p1, p2),
                                α, β, backend...)
 
     codA, domA = codomainind(A), domainind(A)
@@ -288,7 +291,7 @@ function planarcontract!(C::AbstractTensorMap,
                          backend::Backend...)
     # special case only defined for contracting 2 indices
     length(oindB) == length(cindB) == 2 ||
-        return planarcontract!(C, A, (oindA, cindA), copy(B), (cindB, oindB), (p1, p2),
+        return planarcontract!(C, A, (oindA, cindA), TensorMap(B), (cindB, oindB), (p1, p2),
                                α, β, backend...)
 
     codA, domA = codomainind(A), domainind(A)
@@ -334,7 +337,7 @@ end
 function planarcontract!(C::AbstractTensorMap, A::BraidingTensor, pA::Index2Tuple,
                          B::BraidingTensor, pB::Index2Tuple, pC::Index2Tuple,
                          α::Number, β::Number, backend::Backend...)
-    return planarcontract!(C, copy(A), pA, copy(B), pB, pC, α, β, backend...)
+    return planarcontract!(C, TensorMap(A), pA, TensorMap(B), pB, pC, α, β, backend...)
 end
 
 function planartrace!(C::AbstractTensorMap,
@@ -342,7 +345,7 @@ function planartrace!(C::AbstractTensorMap,
                       p::Index2Tuple, q::Index2Tuple,
                       α::Number, β::Number,
                       backend::Backend...)
-    return planartrace!(C, copy(A), p, q, α, β, backend...)
+    return planartrace!(C, TensorMap(A), p, q, α, β, backend...)
 end
 
 # function planarcontract!(C::AbstractTensorMap{S,N₁,N₂},
diff --git a/src/tensors/tensor.jl b/src/tensors/tensor.jl
index 80f9d0bb..82daa4bf 100644
--- a/src/tensors/tensor.jl
+++ b/src/tensors/tensor.jl
@@ -571,7 +571,8 @@ since it assumes a  uniquely defined coupled charge.
 """
 @inline function Base.getindex(t::TensorMap, sectors::Tuple{I,Vararg{I}}) where {I<:Sector}
     I === sectortype(t) || throw(SectorMismatch("Not a valid sectortype for this tensor."))
-    length(sectors) == numind(t) || throw(ArgumentError("Number of sectors does not match."))
+    length(sectors) == numind(t) ||
+        throw(ArgumentError("Number of sectors does not match."))
     FusionStyle(I) isa UniqueFusion ||
         throw(SectorMismatch("Indexing with sectors only possible if unique fusion"))
     s1 = TupleTools.getindices(sectors, codomainind(t))
@@ -757,7 +758,7 @@ function Base.convert(T::Type{<:TensorMap{E,S,N₁,N₂}},
         return t
     else
         data = Dict{sectortype(T),storagetype(T)}(c => convert(storagetype(T), b)
-                                      for (c, b) in blocks(t))
+                                                  for (c, b) in blocks(t))
         return TensorMap(data, codomain(t), domain(t))
     end
 end
diff --git a/test/braidingtensor.jl b/test/braidingtensor.jl
index b61abc61..726a28c4 100644
--- a/test/braidingtensor.jl
+++ b/test/braidingtensor.jl
@@ -13,119 +13,119 @@ for V in (V1, V2, V3)
 
     t = TensorMap(randn, V * V' * V' * V, V * V')
 
-    ττ = copy(BraidingTensor(V, V'))
+    ττ = TensorMap(BraidingTensor(V, V'))
     @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-1 -2; 1 2] * t[1 2 -3 -4; -5 -6]
     @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-1 -2; 1 2] * t[1 2 -3 -4; -5 -6]
     @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -2 -1] * t[1 2 -3 -4; -5 -6]
     @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-2 2; -1 1] * t[1 2 -3 -4; -5 -6]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V', V'))
+    ττ = TensorMap(BraidingTensor(V', V'))
     @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6]
     @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6]
     @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -3 -2] * t[-1 1 2 -4; -5 -6]
     @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-3 2; -2 1] * t[-1 1 2 -4; -5 -6]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V', V))
+    ττ = TensorMap(BraidingTensor(V', V))
     @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6]
     @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6]
     @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[-3 2; -2 1] * t[-1 1 2 -4; -5 -6]
     @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V, V'))
+    ττ = TensorMap(BraidingTensor(V, V'))
     @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-3 2; -2 1] * t[-1 1 2 -4; -5 -6]
     @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-3 2; -2 1] * t[-1 1 2 -4; -5 -6]
     @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6]
     @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[2 1; -3 -2] * t[-1 1 2 -4; -5 -6]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V, V))
+    ττ = TensorMap(BraidingTensor(V, V))
     @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[2 1; -3 -2] * t[-1 1 2 -4; -5 -6]
     @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[2 1; -3 -2] * t[-1 1 2 -4; -5 -6]
     @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6]
     @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V', V))
+    ττ = TensorMap(BraidingTensor(V', V))
     @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-3 -4; 1 2] * t[-1 -2 1 2; -5 -6]
     @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-3 -4; 1 2] * t[-1 -2 1 2; -5 -6]
     @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -4 -3] * t[-1 -2 1 2; -5 -6]
     @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-4 2; -3 1] * t[-1 -2 1 2; -5 -6]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V', V))
+    ττ = TensorMap(BraidingTensor(V', V))
     @planar2 t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ[1 2; -5 -6]
     @planar2 t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * ττ[1 2; -5 -6]
     @planar2 t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ[-6 -5; 2 1]
     @planar2 t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ'[2 -6; 1 -5]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V, V'))
+    ττ = TensorMap(BraidingTensor(V, V'))
     @planar2 t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ'[1 2; -5 -6]
     @planar2 t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * ττ'[1 2; -5 -6]
     @planar2 t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ'[-6 -5; 2 1]
     @planar2 t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ[2 -6; 1 -5]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V, V))
+    ττ = TensorMap(BraidingTensor(V, V))
     @planar2 t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * τ[-4 -6; 1 2]
     @planar2 t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * ττ[-4 -6; 1 2]
     @planar2 t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * τ[2 1; -6 -4]
     @planar2 t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * τ'[-6 2; -4 1]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V', V))
+    ττ = TensorMap(BraidingTensor(V', V))
     @planar2 t1[(); (-1, -2)] := τ[2 1; 3 4] * t[1 2 3 4; -1 -2]
     @planar2 t2[(); (-1, -2)] := ττ[2 1; 3 4] * t[1 2 3 4; -1 -2]
     @planar2 t3[(); (-1, -2)] := τ[4 3; 1 2] * t[1 2 3 4; -1 -2]
     @planar2 t4[(); (-1, -2)] := τ'[1 4; 2 3] * t[1 2 3 4; -1 -2]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V, V))
+    ττ = TensorMap(BraidingTensor(V, V))
     @planar2 t1[-1; -2] := τ[2 1; 3 4] * t[-1 1 2 3; -2 4]
     @planar2 t2[-1; -2] := ττ[2 1; 3 4] * t[-1 1 2 3; -2 4]
     @planar2 t3[-1; -2] := τ[4 3; 1 2] * t[-1 1 2 3; -2 4]
     @planar2 t4[-1; -2] := τ'[1 4; 2 3] * t[-1 1 2 3; -2 4]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V, V'))
+    ττ = TensorMap(BraidingTensor(V, V'))
     @planar2 t1[-1 -2] := τ[2 1; 3 4] * t[-1 -2 1 2; 4 3]
     @planar2 t2[-1 -2] := ττ[2 1; 3 4] * t[-1 -2 1 2; 4 3]
     @planar2 t3[-1 -2] := τ[4 3; 1 2] * t[-1 -2 1 2; 4 3]
     @planar2 t4[-1 -2] := τ'[1 4; 2 3] * t[-1 -2 1 2; 4 3]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V, V'))
+    ττ = TensorMap(BraidingTensor(V, V'))
     @planar2 t1[-1 -2; -3 -4] := τ[-1 3; 1 2] * t[1 2 3 -2; -3 -4]
     @planar2 t2[-1 -2; -3 -4] := ττ[-1 3; 1 2] * t[1 2 3 -2; -3 -4]
     @planar2 t3[-1 -2; -3 -4] := τ[2 1; 3 -1] * t[1 2 3 -2; -3 -4]
     @planar2 t4[-1 -2; -3 -4] := τ'[3 2; -1 1] * t[1 2 3 -2; -3 -4]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V', V'))
+    ττ = TensorMap(BraidingTensor(V', V'))
     @planar2 t1[-1 -2; -3 -4] := τ'[-2 3; 1 2] * t[-1 1 2 3; -3 -4]
     @planar2 t2[-1 -2; -3 -4] := ττ'[-2 3; 1 2] * t[-1 1 2 3; -3 -4]
     @planar2 t3[-1 -2; -3 -4] := τ'[2 1; 3 -2] * t[-1 1 2 3; -3 -4]
     @planar2 t4[-1 -2; -3 -4] := τ[3 2; -2 1] * t[-1 1 2 3; -3 -4]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V', V))
+    ττ = TensorMap(BraidingTensor(V', V))
     @planar2 t1[-1 -2 -3; -4] := τ[-3 3; 1 2] * t[-1 -2 1 2; -4 3]
     @planar2 t2[-1 -2 -3; -4] := ττ[-3 3; 1 2] * t[-1 -2 1 2; -4 3]
     @planar2 t3[-1 -2 -3; -4] := τ[2 1; 3 -3] * t[-1 -2 1 2; -4 3]
     @planar2 t4[-1 -2 -3; -4] := τ'[3 2; -3 1] * t[-1 -2 1 2; -4 3]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V', V))
+    ττ = TensorMap(BraidingTensor(V', V))
     @planar2 t1[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ[1 2; -4 3]
     @planar2 t2[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * ττ[1 2; -4 3]
     @planar2 t3[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ[3 -4; 2 1]
     @planar2 t4[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ'[2 3; 1 -4]
     @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
 
-    ττ = copy(BraidingTensor(V, V'))
+    ττ = TensorMap(BraidingTensor(V, V'))
     @planar2 t1[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ'[1 2; -4 3]
     @planar2 t2[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * ττ'[1 2; -4 3]
     @planar2 t3[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ'[3 -4; 2 1]

From d98087ead71e3e05c53ef9954213b9e795b6e7e2 Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 18 Mar 2024 17:12:57 +0100
Subject: [PATCH 06/21] Split Sectors module

---
 src/TensorKit.jl       |  6 ++++-
 src/sectors/groups.jl  |  1 +
 src/sectors/product.jl |  3 +--
 src/sectors/sectors.jl | 57 ++++++++++++++++++++++++++++++++++++++++++
 test/newsectors.jl     | 27 ++++++++++----------
 5 files changed, 78 insertions(+), 16 deletions(-)

diff --git a/src/TensorKit.jl b/src/TensorKit.jl
index a796c9e5..10b88f5c 100644
--- a/src/TensorKit.jl
+++ b/src/TensorKit.jl
@@ -163,11 +163,15 @@ Base.show(io::IO, ::IndexError{Nothing}) = print(io, "IndexError()")
 Base.show(io::IO, e::IndexError) = print(io, "IndexError(", e.message, ")")
 
 # typerepr
-type_repr(T::Type) = repr(T)
+# type_repr(T::Type) = repr(T)
 
 # Definitions and methods for superselection sectors (quantum numbers)
 #----------------------------------------------------------------------
 include("sectors/sectors.jl")
+using .Sectors
+import .Sectors: dim, BraidingStyle, FusionStyle, ⊠, ⊗
+import .Sectors: dual, type_repr
+import .Sectors: twist
 
 # Constructing and manipulating fusion trees and iterators thereof
 #------------------------------------------------------------------
diff --git a/src/sectors/groups.jl b/src/sectors/groups.jl
index eab121fa..da5bb7a0 100644
--- a/src/sectors/groups.jl
+++ b/src/sectors/groups.jl
@@ -17,6 +17,7 @@ type_repr(::Type{ℤ₂}) = "ℤ₂"
 type_repr(::Type{ℤ₃}) = "ℤ₃"
 type_repr(::Type{ℤ₄}) = "ℤ₄"
 type_repr(::Type{SU₂}) = "SU₂"
+type_repr(T::Type) = repr(T)
 
 const GroupTuple = Tuple{Vararg{Group}}
 
diff --git a/src/sectors/product.jl b/src/sectors/product.jl
index 3bc45526..4331504c 100644
--- a/src/sectors/product.jl
+++ b/src/sectors/product.jl
@@ -74,7 +74,7 @@ function Nsymbol(a::P, b::P, c::P) where {P<:ProductSector}
 end
 
 _firstsector(x::ProductSector) = x.sectors[1]
-_tailsector(x::ProductSector) = ProductSector(tail(x.sectors))
+_tailsector(x::ProductSector) = ProductSector(Base.tail(x.sectors))
 
 function Fsymbol(a::P, b::P, c::P, d::P, e::P, f::P) where {P<:ProductSector}
     heads = map(_firstsector, (a, b, c, d, e, f))
@@ -196,7 +196,6 @@ group representations, we have `Irrep[G₁] ⊠ Irrep[G₂] == Irrep[G₁ × G
 ⊠(s1::Sector, p2::ProductSector) = ProductSector(tuple(s1, p2.sectors...))
 ⊠(p1::ProductSector, p2::ProductSector) = ProductSector(tuple(p1.sectors..., p2.sectors...))
 
-# grow types from the left using Base.tuple_type_cons
 ⊠(I1::Type{Trivial}, I2::Type{Trivial}) = Trivial
 ⊠(I1::Type{Trivial}, I2::Type{<:ProductSector}) = I2
 ⊠(I1::Type{Trivial}, I2::Type{<:Sector}) = I2
diff --git a/src/sectors/sectors.jl b/src/sectors/sectors.jl
index bfb83312..4f59689c 100644
--- a/src/sectors/sectors.jl
+++ b/src/sectors/sectors.jl
@@ -1,6 +1,45 @@
 # Superselection sectors (quantum numbers):
 # for defining graded vector spaces and invariant subspaces of tensor products
 #==============================================================================#
+
+module Sectors
+
+# exports
+export Sector, Group, AbstractIrrep
+export Irrep
+
+export Nsymbol, Fsymbol, Rsymbol, Asymbol, Bsymbol
+export dim, sqrtdim, isqrtdim, frobeniusschur, twist, fusiontensor
+export deligneproduct
+export FusionStyle, UniqueFusion, MultipleFusion, SimpleFusion, GenericFusion,
+       MultiplicityFreeFusion
+export BraidingStyle, NoBraiding, SymmetricBraiding, Bosonic, Fermionic, Anyonic
+export SectorSet, SectorValues, findindex, vertex_ind2label, vertex_labeltype
+
+export pentagon_equation, hexagon_equation
+
+export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep
+export ProductSector
+export FermionParity, FermionNumber, FermionSpin
+export PlanarTrivial, FibonacciAnyon, IsingAnyon
+
+# unicode
+export ⊠, ⊗, ×
+export ℤ, ℤ₂, ℤ₃, ℤ₄, U₁, SU, SU₂, CU₁
+export fℤ₂, fU₁, fSU₂
+
+# imports
+using Base: SizeUnknown, HasLength, IsInfinite
+using Base: HasEltype, EltypeUnknown
+using Base.Iterators: product, filter
+using Base: tuple_type_head, tuple_type_tail
+
+using LinearAlgebra: tr
+using TensorOperations
+using HalfIntegers
+using WignerSymbols
+
+
 """
     abstract type Sector end
 
@@ -464,3 +503,21 @@ 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
+
+# utility
+_kron(A, B, C, D...) = _kron(_kron(A, B), C, D...)
+function _kron(A, B)
+    sA = size(A)
+    sB = size(B)
+    s = map(*, sA, sB)
+    C = similar(A, promote_type(eltype(A), eltype(B)), s)
+    for IA in eachindex(IndexCartesian(), A)
+        for IB in eachindex(IndexCartesian(), B)
+            I = CartesianIndex(IB.I .+ (IA.I .- 1) .* sB)
+            C[I] = A[IA] * B[IB]
+        end
+    end
+    return C
+end
+
+end
diff --git a/test/newsectors.jl b/test/newsectors.jl
index 4aa5e606..61563136 100644
--- a/test/newsectors.jl
+++ b/test/newsectors.jl
@@ -5,9 +5,10 @@ module NewSectors
 
 export NewSU2Irrep
 
-using HalfIntegers, WignerSymbols, TensorKit
+using HalfIntegers, WignerSymbols
+using TensorKit.Sectors
 
-struct NewSU2Irrep <: TensorKit.Sector
+struct NewSU2Irrep <: Sector
     j::HalfInt
     function NewSU2Irrep(j)
         j >= zero(j) || error("Not a valid SU₂ irrep")
@@ -19,26 +20,26 @@ Base.convert(::Type{NewSU2Irrep}, j::Real) = NewSU2Irrep(j)
 const _su2one = NewSU2Irrep(zero(HalfInt))
 Base.one(::Type{NewSU2Irrep}) = _su2one
 Base.conj(s::NewSU2Irrep) = s
-function TensorKit.:⊗(s1::NewSU2Irrep, s2::NewSU2Irrep)
-    return TensorKit.SectorSet{NewSU2Irrep}(abs(s1.j - s2.j):(s1.j + s2.j))
+function Sectors.:⊗(s1::NewSU2Irrep, s2::NewSU2Irrep)
+    return SectorSet{NewSU2Irrep}(abs(s1.j - s2.j):(s1.j + s2.j))
 end
 
-Base.IteratorSize(::Type{TensorKit.SectorValues{NewSU2Irrep}}) = Base.IsInfinite()
-Base.iterate(::TensorKit.SectorValues{NewSU2Irrep}, i=0) = (NewSU2Irrep(half(i)), i + 1)
+Base.IteratorSize(::Type{SectorValues{NewSU2Irrep}}) = Base.IsInfinite()
+Base.iterate(::SectorValues{NewSU2Irrep}, i=0) = (NewSU2Irrep(half(i)), i + 1)
 # Base.getindex(::SectorValues{NewSU2Irrep}, i::Int) =
 #     1 <= i ? NewSU2Irrep(half(i-1)) : throw(BoundsError(values(NewSU2Irrep), i))
 # findindex(::SectorValues{NewSU2Irrep}, s::NewSU2Irrep) = twice(s.j)+1
 
 # TensorKit.dim(s::NewSU2Irrep) = twice(s.j)+1
 #
-TensorKit.FusionStyle(::Type{NewSU2Irrep}) = GenericFusion()
-TensorKit.BraidingStyle(::Type{NewSU2Irrep}) = Bosonic()
+Sectors.FusionStyle(::Type{NewSU2Irrep}) = GenericFusion()
+Sectors.BraidingStyle(::Type{NewSU2Irrep}) = Bosonic()
 Base.isreal(::Type{NewSU2Irrep}) = true
 
-function TensorKit.Nsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
+function Sectors.Nsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
     return convert(Int, WignerSymbols.δ(sa.j, sb.j, sc.j))
 end
-function TensorKit.Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep,
+function Sectors.Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep,
                            s4::NewSU2Irrep, s5::NewSU2Irrep, s6::NewSU2Irrep)
     n1 = Nsymbol(s1, s2, s5)
     n2 = Nsymbol(s5, s3, s4)
@@ -49,13 +50,13 @@ function TensorKit.Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep,
                                                        s4.j, s3.j, s5.j, s6.j)
     return fill(f, (n1, n2, n3, n4))
 end
-function TensorKit.Rsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
+function Sectors.Rsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
     Nsymbol(sa, sb, sc) > 0 || return fill(0.0, (0, 0))
     return fill(iseven(convert(Int, sa.j + sb.j - sc.j)) ? 1.0 : -1.0, (1, 1))
 end
-TensorKit.dim(s::NewSU2Irrep) = twice(s.j) + 1
+Sectors.dim(s::NewSU2Irrep) = twice(s.j) + 1
 
-function TensorKit.fusiontensor(a::NewSU2Irrep, b::NewSU2Irrep, c::NewSU2Irrep)
+function Sectors.fusiontensor(a::NewSU2Irrep, b::NewSU2Irrep, c::NewSU2Irrep)
     C = Array{Float64}(undef, dim(a), dim(b), dim(c), 1)
     ja, jb, jc = a.j, b.j, c.j
 

From 7ebbd1218fbe163385205bd4791c561ba17ca06a Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 18 Mar 2024 18:01:43 +0100
Subject: [PATCH 07/21] convert submodule into its own module

---
 Project.toml                                  |  3 +-
 TensorKitSectors/Project.toml                 | 19 +++++
 TensorKitSectors/src/TensorKitSectors.jl      | 55 +++++++++++++
 .../src}/anyons.jl                            |  0
 TensorKitSectors/src/auxiliary.jl             | 15 ++++
 .../src}/fermions.jl                          |  0
 .../src}/groups.jl                            |  0
 .../src}/irreps.jl                            |  0
 .../src}/product.jl                           |  0
 .../src}/sectors.jl                           | 67 ----------------
 TensorKitSectors/test/newsectors.jl           | 79 +++++++++++++++++++
 TensorKitSectors/test/runtests.jl             | 53 +++++++++++++
 {test => TensorKitSectors/test}/sectors.jl    | 16 ++--
 src/TensorKit.jl                              | 10 +--
 test/newsectors.jl                            | 26 +++---
 test/runtests.jl                              |  1 -
 16 files changed, 252 insertions(+), 92 deletions(-)
 create mode 100644 TensorKitSectors/Project.toml
 create mode 100644 TensorKitSectors/src/TensorKitSectors.jl
 rename {src/sectors => TensorKitSectors/src}/anyons.jl (100%)
 create mode 100644 TensorKitSectors/src/auxiliary.jl
 rename {src/sectors => TensorKitSectors/src}/fermions.jl (100%)
 rename {src/sectors => TensorKitSectors/src}/groups.jl (100%)
 rename {src/sectors => TensorKitSectors/src}/irreps.jl (100%)
 rename {src/sectors => TensorKitSectors/src}/product.jl (100%)
 rename {src/sectors => TensorKitSectors/src}/sectors.jl (90%)
 create mode 100644 TensorKitSectors/test/newsectors.jl
 create mode 100644 TensorKitSectors/test/runtests.jl
 rename {test => TensorKitSectors/test}/sectors.jl (90%)

diff --git a/Project.toml b/Project.toml
index 5299ed3f..523731dd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
+TensorKitSectors = "7ea117e4-0799-4282-ac6e-748b4bc2d0f5"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
@@ -27,8 +28,8 @@ ChainRulesTestUtils = "1"
 Combinatorics = "1"
 FiniteDifferences = "0.12"
 HalfIntegers = "1"
-LinearAlgebra = "1"
 LRUCache = "1.0.2"
+LinearAlgebra = "1"
 PackageExtensionCompat = "1"
 Random = "1"
 Strided = "2"
diff --git a/TensorKitSectors/Project.toml b/TensorKitSectors/Project.toml
new file mode 100644
index 00000000..6ed915cd
--- /dev/null
+++ b/TensorKitSectors/Project.toml
@@ -0,0 +1,19 @@
+name = "TensorKitSectors"
+uuid = "7ea117e4-0799-4282-ac6e-748b4bc2d0f5"
+authors = ["Jutho Haegeman", "Lukas Devos"]
+version = "0.1.0"
+
+[deps]
+HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
+WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
+
+[extras]
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
+
+[targets]
+test = ["Test", "TestExtras", "TensorOperations", "Random"]
diff --git a/TensorKitSectors/src/TensorKitSectors.jl b/TensorKitSectors/src/TensorKitSectors.jl
new file mode 100644
index 00000000..8ac27194
--- /dev/null
+++ b/TensorKitSectors/src/TensorKitSectors.jl
@@ -0,0 +1,55 @@
+# Superselection sectors (quantum numbers):
+# for defining graded vector spaces and invariant subspaces of tensor products
+#==========================================================================================#
+
+module TensorKitSectors
+
+# exports
+# -------
+export Sector, Group, AbstractIrrep
+export Irrep
+
+export Nsymbol, Fsymbol, Rsymbol, Asymbol, Bsymbol
+export dim, sqrtdim, isqrtdim, frobeniusschur, twist, fusiontensor, dual
+export deligneproduct
+export FusionStyle, UniqueFusion, MultipleFusion, SimpleFusion, GenericFusion,
+       MultiplicityFreeFusion
+export BraidingStyle, NoBraiding, SymmetricBraiding, Bosonic, Fermionic, Anyonic
+export SectorSet, SectorValues, findindex, vertex_ind2label, vertex_labeltype
+
+export pentagon_equation, hexagon_equation
+
+export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep
+export ProductSector
+export FermionParity, FermionNumber, FermionSpin
+export PlanarTrivial, FibonacciAnyon, IsingAnyon
+
+# unicode exports
+# ---------------
+export ⊠, ⊗, ×
+export ℤ, ℤ₂, ℤ₃, ℤ₄, U₁, SU, SU₂, CU₁
+export fℤ₂, fU₁, fSU₂
+
+# imports
+# -------
+using Base: SizeUnknown, HasLength, IsInfinite
+using Base: HasEltype, EltypeUnknown
+using Base.Iterators: product, filter
+using Base: tuple_type_head, tuple_type_tail
+
+using LinearAlgebra: tr
+using TensorOperations
+using HalfIntegers
+using WignerSymbols
+
+# includes
+# --------
+include("auxiliary.jl")
+include("sectors.jl")
+include("groups.jl")
+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
+
+end # module TensorKitSectors
diff --git a/src/sectors/anyons.jl b/TensorKitSectors/src/anyons.jl
similarity index 100%
rename from src/sectors/anyons.jl
rename to TensorKitSectors/src/anyons.jl
diff --git a/TensorKitSectors/src/auxiliary.jl b/TensorKitSectors/src/auxiliary.jl
new file mode 100644
index 00000000..46486da8
--- /dev/null
+++ b/TensorKitSectors/src/auxiliary.jl
@@ -0,0 +1,15 @@
+# kronecker product with arbitrary number of dimensions
+_kron(A, B, C, D...) = _kron(_kron(A, B), C, D...)
+function _kron(A, B)
+    sA = size(A)
+    sB = size(B)
+    s = map(*, sA, sB)
+    C = similar(A, promote_type(eltype(A), eltype(B)), s)
+    for IA in eachindex(IndexCartesian(), A)
+        for IB in eachindex(IndexCartesian(), B)
+            I = CartesianIndex(IB.I .+ (IA.I .- 1) .* sB)
+            C[I] = A[IA] * B[IB]
+        end
+    end
+    return C
+end
diff --git a/src/sectors/fermions.jl b/TensorKitSectors/src/fermions.jl
similarity index 100%
rename from src/sectors/fermions.jl
rename to TensorKitSectors/src/fermions.jl
diff --git a/src/sectors/groups.jl b/TensorKitSectors/src/groups.jl
similarity index 100%
rename from src/sectors/groups.jl
rename to TensorKitSectors/src/groups.jl
diff --git a/src/sectors/irreps.jl b/TensorKitSectors/src/irreps.jl
similarity index 100%
rename from src/sectors/irreps.jl
rename to TensorKitSectors/src/irreps.jl
diff --git a/src/sectors/product.jl b/TensorKitSectors/src/product.jl
similarity index 100%
rename from src/sectors/product.jl
rename to TensorKitSectors/src/product.jl
diff --git a/src/sectors/sectors.jl b/TensorKitSectors/src/sectors.jl
similarity index 90%
rename from src/sectors/sectors.jl
rename to TensorKitSectors/src/sectors.jl
index 4f59689c..1cf2e3e5 100644
--- a/src/sectors/sectors.jl
+++ b/TensorKitSectors/src/sectors.jl
@@ -1,45 +1,3 @@
-# Superselection sectors (quantum numbers):
-# for defining graded vector spaces and invariant subspaces of tensor products
-#==============================================================================#
-
-module Sectors
-
-# exports
-export Sector, Group, AbstractIrrep
-export Irrep
-
-export Nsymbol, Fsymbol, Rsymbol, Asymbol, Bsymbol
-export dim, sqrtdim, isqrtdim, frobeniusschur, twist, fusiontensor
-export deligneproduct
-export FusionStyle, UniqueFusion, MultipleFusion, SimpleFusion, GenericFusion,
-       MultiplicityFreeFusion
-export BraidingStyle, NoBraiding, SymmetricBraiding, Bosonic, Fermionic, Anyonic
-export SectorSet, SectorValues, findindex, vertex_ind2label, vertex_labeltype
-
-export pentagon_equation, hexagon_equation
-
-export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep
-export ProductSector
-export FermionParity, FermionNumber, FermionSpin
-export PlanarTrivial, FibonacciAnyon, IsingAnyon
-
-# unicode
-export ⊠, ⊗, ×
-export ℤ, ℤ₂, ℤ₃, ℤ₄, U₁, SU, SU₂, CU₁
-export fℤ₂, fU₁, fSU₂
-
-# imports
-using Base: SizeUnknown, HasLength, IsInfinite
-using Base: HasEltype, EltypeUnknown
-using Base.Iterators: product, filter
-using Base: tuple_type_head, tuple_type_tail
-
-using LinearAlgebra: tr
-using TensorOperations
-using HalfIntegers
-using WignerSymbols
-
-
 """
     abstract type Sector end
 
@@ -496,28 +454,3 @@ function Base.iterate(s::SectorSet{I}, args...) where {I<:Sector}
     val, state = next
     return convert(I, s.f(val)), state
 end
-
-# possible sectors
-include("groups.jl")
-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
-
-# utility
-_kron(A, B, C, D...) = _kron(_kron(A, B), C, D...)
-function _kron(A, B)
-    sA = size(A)
-    sB = size(B)
-    s = map(*, sA, sB)
-    C = similar(A, promote_type(eltype(A), eltype(B)), s)
-    for IA in eachindex(IndexCartesian(), A)
-        for IB in eachindex(IndexCartesian(), B)
-            I = CartesianIndex(IB.I .+ (IA.I .- 1) .* sB)
-            C[I] = A[IA] * B[IB]
-        end
-    end
-    return C
-end
-
-end
diff --git a/TensorKitSectors/test/newsectors.jl b/TensorKitSectors/test/newsectors.jl
new file mode 100644
index 00000000..16cd8ef0
--- /dev/null
+++ b/TensorKitSectors/test/newsectors.jl
@@ -0,0 +1,79 @@
+# NewSU2Irrep
+# Test of a bare-bones sector implementation, which is set to be `DegenerateNonAbelian`
+# (even though it is not)
+module NewSectors
+
+export NewSU2Irrep
+
+using HalfIntegers, WignerSymbols
+using TensorKitSectors
+
+import TensorKitSectors: FusionStyle, BraidingStyle, Nsymbol, Fsymbol, Rsymbol, dim,
+                         fusiontensor, ⊗
+
+struct NewSU2Irrep <: Sector
+    j::HalfInt
+    function NewSU2Irrep(j)
+        j >= zero(j) || error("Not a valid SU₂ irrep")
+        return new(j)
+    end
+end
+Base.convert(::Type{NewSU2Irrep}, j::Real) = NewSU2Irrep(j)
+
+const _su2one = NewSU2Irrep(zero(HalfInt))
+Base.one(::Type{NewSU2Irrep}) = _su2one
+Base.conj(s::NewSU2Irrep) = s
+function ⊗(s1::NewSU2Irrep, s2::NewSU2Irrep)
+    return SectorSet{NewSU2Irrep}(abs(s1.j - s2.j):(s1.j + s2.j))
+end
+
+Base.IteratorSize(::Type{SectorValues{NewSU2Irrep}}) = Base.IsInfinite()
+Base.iterate(::SectorValues{NewSU2Irrep}, i=0) = (NewSU2Irrep(half(i)), i + 1)
+# Base.getindex(::SectorValues{NewSU2Irrep}, i::Int) =
+#     1 <= i ? NewSU2Irrep(half(i-1)) : throw(BoundsError(values(NewSU2Irrep), i))
+# findindex(::SectorValues{NewSU2Irrep}, s::NewSU2Irrep) = twice(s.j)+1
+
+# TensorKit.dim(s::NewSU2Irrep) = twice(s.j)+1
+#
+FusionStyle(::Type{NewSU2Irrep}) = GenericFusion()
+BraidingStyle(::Type{NewSU2Irrep}) = Bosonic()
+Base.isreal(::Type{NewSU2Irrep}) = true
+
+function Nsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
+    return convert(Int, WignerSymbols.δ(sa.j, sb.j, sc.j))
+end
+
+function Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep,
+                 s4::NewSU2Irrep, s5::NewSU2Irrep, s6::NewSU2Irrep)
+    n1 = Nsymbol(s1, s2, s5)
+    n2 = Nsymbol(s5, s3, s4)
+    n3 = Nsymbol(s2, s3, s6)
+    n4 = Nsymbol(s1, s6, s4)
+    f = all(==(_su2one), (s1, s2, s3, s4, s5, s6)) ? 1.0 :
+        sqrt(dim(s5) * dim(s6)) * WignerSymbols.racahW(Float64, s1.j, s2.j,
+                                                       s4.j, s3.j, s5.j, s6.j)
+    return fill(f, (n1, n2, n3, n4))
+end
+
+function Rsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
+    Nsymbol(sa, sb, sc) > 0 || return fill(0.0, (0, 0))
+    return fill(iseven(convert(Int, sa.j + sb.j - sc.j)) ? 1.0 : -1.0, (1, 1))
+end
+
+dim(s::NewSU2Irrep) = twice(s.j) + 1
+
+function fusiontensor(a::NewSU2Irrep, b::NewSU2Irrep, c::NewSU2Irrep)
+    C = Array{Float64}(undef, dim(a), dim(b), dim(c), 1)
+    ja, jb, jc = a.j, b.j, c.j
+
+    for kc in 1:dim(c), kb in 1:dim(b), ka in 1:dim(a)
+        C[ka, kb, kc, 1] = WignerSymbols.clebschgordan(ja, ja + 1 - ka, jb, jb + 1 - kb, jc,
+                                                       jc + 1 - kc)
+    end
+    return C
+end
+
+Base.hash(s::NewSU2Irrep, h::UInt) = hash(s.j, h)
+Base.isless(s1::NewSU2Irrep, s2::NewSU2Irrep) = isless(s1.j, s2.j)
+
+end
diff --git a/TensorKitSectors/test/runtests.jl b/TensorKitSectors/test/runtests.jl
new file mode 100644
index 00000000..338e2f76
--- /dev/null
+++ b/TensorKitSectors/test/runtests.jl
@@ -0,0 +1,53 @@
+using Test
+using TestExtras
+using Random
+using TensorKitSectors
+using TensorOperations
+using Base.Iterators: take, product
+using LinearAlgebra: LinearAlgebra
+
+const TKS = TensorKitSectors
+
+include("newsectors.jl")
+using .NewSectors
+
+smallset(::Type{I}) where {I<:Sector} = take(values(I), 5)
+function smallset(::Type{ProductSector{Tuple{I1,I2}}}) where {I1,I2}
+    iter = product(smallset(I1), smallset(I2))
+    s = collect(i ⊠ j for (i, j) in iter if dim(i) * dim(j) <= 6)
+    return length(s) > 6 ? rand(s, 6) : s
+end
+function smallset(::Type{ProductSector{Tuple{I1,I2,I3}}}) where {I1,I2,I3}
+    iter = product(smallset(I1), smallset(I2), smallset(I3))
+    s = collect(i ⊠ j ⊠ k for (i, j, k) in iter if dim(i) * dim(j) * dim(k) <= 6)
+    return length(s) > 6 ? rand(s, 6) : s
+end
+function randsector(::Type{I}) where {I<:Sector}
+    s = collect(smallset(I))
+    a = rand(s)
+    while a == one(a) # don't use trivial label
+        a = rand(s)
+    end
+    return a
+end
+function hasfusiontensor(I::Type{<:Sector})
+    try
+        fusiontensor(one(I), one(I), one(I))
+        return true
+    catch e
+        if e isa MethodError
+            return false
+        else
+            rethrow(e)
+        end
+    end
+end
+
+sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irrep,
+              FibonacciAnyon, IsingAnyon, FermionParity, FermionParity ⊠ FermionParity,
+              Z3Irrep ⊠ Z4Irrep, FermionParity ⊠ U1Irrep ⊠ SU2Irrep,
+              FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep,
+              NewSU2Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ NewSU2Irrep,
+              Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)
+
+include("sectors.jl")
diff --git a/test/sectors.jl b/TensorKitSectors/test/sectors.jl
similarity index 90%
rename from test/sectors.jl
rename to TensorKitSectors/test/sectors.jl
index a0702552..efc5ab37 100644
--- a/test/sectors.jl
+++ b/TensorKitSectors/test/sectors.jl
@@ -2,12 +2,12 @@ println("------------------------------------")
 println("Sectors")
 println("------------------------------------")
 @timedtestset "Sectors" verbose = true begin
-    @timedtestset "Sector properties of $(TensorKit.type_repr(I))" for I in sectorlist
-        Istr = TensorKit.type_repr(I)
+    @timedtestset "Sector properties of $(TKS.type_repr(I))" for I in sectorlist
+        Istr = TKS.type_repr(I)
         @testset "Sector $Istr: Basic properties" begin
             s = (randsector(I), randsector(I), randsector(I))
             @test eval(Meta.parse(sprint(show, I))) == I
-            @test eval(Meta.parse(TensorKit.type_repr(I))) == I
+            @test eval(Meta.parse(TKS.type_repr(I))) == I
             @test eval(Meta.parse(sprint(show, s[1]))) == s[1]
             @test @constinferred(hash(s[1])) == hash(deepcopy(s[1]))
             @test @constinferred(one(s[1])) == @constinferred(one(I))
@@ -30,21 +30,21 @@ println("------------------------------------")
                 @test !isless(s, sprev) # confirm compatibility with sort order
                 if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
                     @test_throws ArgumentError values(I)[i]
-                    @test_throws ArgumentError TensorKit.findindex(values(I), s)
+                    @test_throws ArgumentError findindex(values(I), s)
                 elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int})
                     @test s == @constinferred (values(I)[i])
-                    @test TensorKit.findindex(values(I), s) == i
+                    @test findindex(values(I), s) == i
                 end
                 sprev = s
                 i >= 10 && break
             end
             @test one(I) == first(values(I))
             if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
-                @test_throws ArgumentError TensorKit.findindex(values(I), one(I))
+                @test_throws ArgumentError findindex(values(I), one(I))
             elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int})
-                @test (@constinferred TensorKit.findindex(values(I), one(I))) == 1
+                @test (@constinferred findindex(values(I), one(I))) == 1
                 for s in smallset(I)
-                    @test (@constinferred values(I)[TensorKit.findindex(values(I), s)]) == s
+                    @test (@constinferred values(I)[findindex(values(I), s)]) == s
                 end
             end
         end
diff --git a/src/TensorKit.jl b/src/TensorKit.jl
index 10b88f5c..a7756643 100644
--- a/src/TensorKit.jl
+++ b/src/TensorKit.jl
@@ -167,11 +167,11 @@ Base.show(io::IO, e::IndexError) = print(io, "IndexError(", e.message, ")")
 
 # Definitions and methods for superselection sectors (quantum numbers)
 #----------------------------------------------------------------------
-include("sectors/sectors.jl")
-using .Sectors
-import .Sectors: dim, BraidingStyle, FusionStyle, ⊠, ⊗
-import .Sectors: dual, type_repr
-import .Sectors: twist
+
+using TensorKitSectors
+import TensorKitSectors: dim, BraidingStyle, FusionStyle, ⊠, ⊗
+import TensorKitSectors: dual, type_repr
+import TensorKitSectors: twist
 
 # Constructing and manipulating fusion trees and iterators thereof
 #------------------------------------------------------------------
diff --git a/test/newsectors.jl b/test/newsectors.jl
index 61563136..16cd8ef0 100644
--- a/test/newsectors.jl
+++ b/test/newsectors.jl
@@ -6,7 +6,10 @@ module NewSectors
 export NewSU2Irrep
 
 using HalfIntegers, WignerSymbols
-using TensorKit.Sectors
+using TensorKitSectors
+
+import TensorKitSectors: FusionStyle, BraidingStyle, Nsymbol, Fsymbol, Rsymbol, dim,
+                         fusiontensor, ⊗
 
 struct NewSU2Irrep <: Sector
     j::HalfInt
@@ -20,7 +23,7 @@ Base.convert(::Type{NewSU2Irrep}, j::Real) = NewSU2Irrep(j)
 const _su2one = NewSU2Irrep(zero(HalfInt))
 Base.one(::Type{NewSU2Irrep}) = _su2one
 Base.conj(s::NewSU2Irrep) = s
-function Sectors.:⊗(s1::NewSU2Irrep, s2::NewSU2Irrep)
+function ⊗(s1::NewSU2Irrep, s2::NewSU2Irrep)
     return SectorSet{NewSU2Irrep}(abs(s1.j - s2.j):(s1.j + s2.j))
 end
 
@@ -32,15 +35,16 @@ Base.iterate(::SectorValues{NewSU2Irrep}, i=0) = (NewSU2Irrep(half(i)), i + 1)
 
 # TensorKit.dim(s::NewSU2Irrep) = twice(s.j)+1
 #
-Sectors.FusionStyle(::Type{NewSU2Irrep}) = GenericFusion()
-Sectors.BraidingStyle(::Type{NewSU2Irrep}) = Bosonic()
+FusionStyle(::Type{NewSU2Irrep}) = GenericFusion()
+BraidingStyle(::Type{NewSU2Irrep}) = Bosonic()
 Base.isreal(::Type{NewSU2Irrep}) = true
 
-function Sectors.Nsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
+function Nsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
     return convert(Int, WignerSymbols.δ(sa.j, sb.j, sc.j))
 end
-function Sectors.Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep,
-                           s4::NewSU2Irrep, s5::NewSU2Irrep, s6::NewSU2Irrep)
+
+function Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep,
+                 s4::NewSU2Irrep, s5::NewSU2Irrep, s6::NewSU2Irrep)
     n1 = Nsymbol(s1, s2, s5)
     n2 = Nsymbol(s5, s3, s4)
     n3 = Nsymbol(s2, s3, s6)
@@ -50,13 +54,15 @@ function Sectors.Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep,
                                                        s4.j, s3.j, s5.j, s6.j)
     return fill(f, (n1, n2, n3, n4))
 end
-function Sectors.Rsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
+
+function Rsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
     Nsymbol(sa, sb, sc) > 0 || return fill(0.0, (0, 0))
     return fill(iseven(convert(Int, sa.j + sb.j - sc.j)) ? 1.0 : -1.0, (1, 1))
 end
-Sectors.dim(s::NewSU2Irrep) = twice(s.j) + 1
 
-function Sectors.fusiontensor(a::NewSU2Irrep, b::NewSU2Irrep, c::NewSU2Irrep)
+dim(s::NewSU2Irrep) = twice(s.j) + 1
+
+function fusiontensor(a::NewSU2Irrep, b::NewSU2Irrep, c::NewSU2Irrep)
     C = Array{Float64}(undef, dim(a), dim(b), dim(c), 1)
     ja, jb, jc = a.j, b.j, c.j
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 21760665..9705fc57 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -57,7 +57,6 @@ sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irre
               Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)
 
 Ti = time()
-include("sectors.jl")
 include("fusiontrees.jl")
 include("spaces.jl")
 include("tensors.jl")

From e9735b69216c43eb844641f0d36be2e6a6d6c446 Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 18 Mar 2024 18:59:31 +0100
Subject: [PATCH 08/21] Organize tests

---
 TensorKitSectors/Project.toml       |   4 +-
 TensorKitSectors/test/newsectors.jl |  19 ++-
 TensorKitSectors/test/runtests.jl   |  61 ++++----
 TensorKitSectors/test/sectors.jl    | 215 +++++++++++++---------------
 TensorKitSectors/test/testsetup.jl  |  55 +++++++
 5 files changed, 187 insertions(+), 167 deletions(-)
 create mode 100644 TensorKitSectors/test/testsetup.jl

diff --git a/TensorKitSectors/Project.toml b/TensorKitSectors/Project.toml
index 6ed915cd..7b915d17 100644
--- a/TensorKitSectors/Project.toml
+++ b/TensorKitSectors/Project.toml
@@ -11,9 +11,11 @@ WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Test", "TestExtras", "TensorOperations", "Random"]
+test = ["Test", "TestExtras", "TestItemRunner", "TensorOperations", "Random", "Reexport"]
diff --git a/TensorKitSectors/test/newsectors.jl b/TensorKitSectors/test/newsectors.jl
index 16cd8ef0..c1bbe803 100644
--- a/TensorKitSectors/test/newsectors.jl
+++ b/TensorKitSectors/test/newsectors.jl
@@ -1,11 +1,15 @@
-# NewSU2Irrep
-# Test of a bare-bones sector implementation, which is set to be `DegenerateNonAbelian`
-# (even though it is not)
+"""
+    module NewSectors
+
+Bare-bones implementation of a new sector type, `NewSU2Irrep`, which is set to be of
+`GenericFusion` type and `Bosonic` braiding style.
+"""
 module NewSectors
 
 export NewSU2Irrep
 
-using HalfIntegers, WignerSymbols
+using HalfIntegers
+using WignerSymbols
 using TensorKitSectors
 
 import TensorKitSectors: FusionStyle, BraidingStyle, Nsymbol, Fsymbol, Rsymbol, dim,
@@ -29,12 +33,7 @@ end
 
 Base.IteratorSize(::Type{SectorValues{NewSU2Irrep}}) = Base.IsInfinite()
 Base.iterate(::SectorValues{NewSU2Irrep}, i=0) = (NewSU2Irrep(half(i)), i + 1)
-# Base.getindex(::SectorValues{NewSU2Irrep}, i::Int) =
-#     1 <= i ? NewSU2Irrep(half(i-1)) : throw(BoundsError(values(NewSU2Irrep), i))
-# findindex(::SectorValues{NewSU2Irrep}, s::NewSU2Irrep) = twice(s.j)+1
 
-# TensorKit.dim(s::NewSU2Irrep) = twice(s.j)+1
-#
 FusionStyle(::Type{NewSU2Irrep}) = GenericFusion()
 BraidingStyle(::Type{NewSU2Irrep}) = Bosonic()
 Base.isreal(::Type{NewSU2Irrep}) = true
@@ -76,4 +75,4 @@ end
 Base.hash(s::NewSU2Irrep, h::UInt) = hash(s.j, h)
 Base.isless(s1::NewSU2Irrep, s2::NewSU2Irrep) = isless(s1.j, s2.j)
 
-end
+end # module NewSectors
diff --git a/TensorKitSectors/test/runtests.jl b/TensorKitSectors/test/runtests.jl
index 338e2f76..1cea3b6e 100644
--- a/TensorKitSectors/test/runtests.jl
+++ b/TensorKitSectors/test/runtests.jl
@@ -8,46 +8,33 @@ using LinearAlgebra: LinearAlgebra
 
 const TKS = TensorKitSectors
 
+include("testsetup.jl")
+using .TestSetup
 include("newsectors.jl")
 using .NewSectors
 
-smallset(::Type{I}) where {I<:Sector} = take(values(I), 5)
-function smallset(::Type{ProductSector{Tuple{I1,I2}}}) where {I1,I2}
-    iter = product(smallset(I1), smallset(I2))
-    s = collect(i ⊠ j for (i, j) in iter if dim(i) * dim(j) <= 6)
-    return length(s) > 6 ? rand(s, 6) : s
-end
-function smallset(::Type{ProductSector{Tuple{I1,I2,I3}}}) where {I1,I2,I3}
-    iter = product(smallset(I1), smallset(I2), smallset(I3))
-    s = collect(i ⊠ j ⊠ k for (i, j, k) in iter if dim(i) * dim(j) * dim(k) <= 6)
-    return length(s) > 6 ? rand(s, 6) : s
-end
-function randsector(::Type{I}) where {I<:Sector}
-    s = collect(smallset(I))
-    a = rand(s)
-    while a == one(a) # don't use trivial label
-        a = rand(s)
-    end
-    return a
-end
-function hasfusiontensor(I::Type{<:Sector})
-    try
-        fusiontensor(one(I), one(I), one(I))
-        return true
-    catch e
-        if e isa MethodError
-            return false
-        else
-            rethrow(e)
-        end
-    end
+const sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irrep,
+                    FibonacciAnyon, IsingAnyon, FermionParity,
+                    FermionParity ⊠ FermionParity,
+                    Z3Irrep ⊠ Z4Irrep, FermionParity ⊠ U1Irrep ⊠ SU2Irrep,
+                    FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep,
+                    NewSU2Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ NewSU2Irrep,
+                    Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)
+
+@testset verbose=true "$I" for I in sectorlist
+    @include("sectors.jl")
 end
 
-sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irrep,
-              FibonacciAnyon, IsingAnyon, FermionParity, FermionParity ⊠ FermionParity,
-              Z3Irrep ⊠ Z4Irrep, FermionParity ⊠ U1Irrep ⊠ SU2Irrep,
-              FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep,
-              NewSU2Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ NewSU2Irrep,
-              Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)
+@testset "Deligne product" begin
+    sectorlist′ = (Trivial, sectorlist...)
+    for I1 in sectorlist′, I2 in sectorlist′
+        a = first(smallset(I1))
+        b = first(smallset(I2))
 
-include("sectors.jl")
+        @constinferred a ⊠ b
+        @constinferred a ⊠ b ⊠ a
+        @constinferred a ⊠ b ⊠ a ⊠ b
+        @constinferred I1 ⊠ I2
+        @test typeof(a ⊠ b) == I1 ⊠ I2
+    end
+end
diff --git a/TensorKitSectors/test/sectors.jl b/TensorKitSectors/test/sectors.jl
index efc5ab37..582092c3 100644
--- a/TensorKitSectors/test/sectors.jl
+++ b/TensorKitSectors/test/sectors.jl
@@ -1,134 +1,111 @@
-println("------------------------------------")
-println("Sectors")
-println("------------------------------------")
-@timedtestset "Sectors" verbose = true begin
-    @timedtestset "Sector properties of $(TKS.type_repr(I))" for I in sectorlist
-        Istr = TKS.type_repr(I)
-        @testset "Sector $Istr: Basic properties" begin
-            s = (randsector(I), randsector(I), randsector(I))
-            @test eval(Meta.parse(sprint(show, I))) == I
-            @test eval(Meta.parse(TKS.type_repr(I))) == I
-            @test eval(Meta.parse(sprint(show, s[1]))) == s[1]
-            @test @constinferred(hash(s[1])) == hash(deepcopy(s[1]))
-            @test @constinferred(one(s[1])) == @constinferred(one(I))
-            @constinferred dual(s[1])
-            @constinferred dim(s[1])
-            @constinferred frobeniusschur(s[1])
-            @constinferred Nsymbol(s...)
-            @constinferred Rsymbol(s...)
-            @constinferred Bsymbol(s...)
-            @constinferred Fsymbol(s..., s...)
-            it = @constinferred s[1] ⊗ s[2]
-            @test eltype(it) === I
-            @test collect(it) isa Array{I}
-            @constinferred ⊗(s..., s...)
+Istr = TKS.type_repr(I)
+@testset "Sector $Istr: Basic properties" begin
+    s = (randsector(I), randsector(I), randsector(I))
+    @test eval(Meta.parse(sprint(show, I))) == I
+    @test eval(Meta.parse(TKS.type_repr(I))) == I
+    @test eval(Meta.parse(sprint(show, s[1]))) == s[1]
+    @test @constinferred(hash(s[1])) == hash(deepcopy(s[1]))
+    @test @constinferred(one(s[1])) == @constinferred(one(I))
+    @constinferred dual(s[1])
+    @constinferred dim(s[1])
+    @constinferred frobeniusschur(s[1])
+    @constinferred Nsymbol(s...)
+    @constinferred Rsymbol(s...)
+    @constinferred Bsymbol(s...)
+    @constinferred Fsymbol(s..., s...)
+    it = @constinferred s[1] ⊗ s[2]
+    @constinferred ⊗(s..., s...)
+end
+@testset "Sector $Istr: Value iterator" begin
+    @test eltype(values(I)) == I
+    sprev = one(I)
+    for (i, s) in enumerate(values(I))
+        @test !isless(s, sprev) # confirm compatibility with sort order
+        if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
+            @test_throws ArgumentError values(I)[i]
+            @test_throws ArgumentError findindex(values(I), s)
+        elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int})
+            @test s == @constinferred (values(I)[i])
+            @test findindex(values(I), s) == i
         end
-        @testset "Sector $Istr: Value iterator" begin
-            @test eltype(values(I)) == I
-            sprev = one(I)
-            for (i, s) in enumerate(values(I))
-                @test !isless(s, sprev) # confirm compatibility with sort order
-                if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
-                    @test_throws ArgumentError values(I)[i]
-                    @test_throws ArgumentError findindex(values(I), s)
-                elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int})
-                    @test s == @constinferred (values(I)[i])
-                    @test findindex(values(I), s) == i
-                end
-                sprev = s
-                i >= 10 && break
-            end
-            @test one(I) == first(values(I))
-            if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
-                @test_throws ArgumentError findindex(values(I), one(I))
-            elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int})
-                @test (@constinferred findindex(values(I), one(I))) == 1
-                for s in smallset(I)
-                    @test (@constinferred values(I)[findindex(values(I), s)]) == s
-                end
-            end
+        sprev = s
+        i >= 10 && break
+    end
+    @test one(I) == first(values(I))
+    if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
+        @test_throws ArgumentError findindex(values(I), one(I))
+    elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int})
+        @test (@constinferred findindex(values(I), one(I))) == 1
+        for s in smallset(I)
+            @test (@constinferred values(I)[findindex(values(I), s)]) == s
         end
-        if BraidingStyle(I) isa Bosonic && hasfusiontensor(I)
-            @testset "Sector $I: fusion tensor and F-move and R-move" begin
-                for a in smallset(I), b in smallset(I)
-                    for c in ⊗(a, b)
-                        X1 = permutedims(fusiontensor(a, b, c), (2, 1, 3, 4))
-                        X2 = fusiontensor(b, a, c)
-                        l = dim(a) * dim(b) * dim(c)
-                        R = LinearAlgebra.transpose(Rsymbol(a, b, c))
-                        sz = (l, convert(Int, Nsymbol(a, b, c)))
-                        @test reshape(X1, sz) ≈ reshape(X2, sz) * R
-                    end
-                end
-                for a in smallset(I), b in smallset(I), c in smallset(I)
-                    for e in ⊗(a, b), f in ⊗(b, c)
-                        for d in intersect(⊗(e, c), ⊗(a, f))
-                            X1 = fusiontensor(a, b, e)
-                            X2 = fusiontensor(e, c, d)
-                            Y1 = fusiontensor(b, c, f)
-                            Y2 = fusiontensor(a, f, d)
-                            @tensor f1[-1, -2, -3, -4] := conj(Y2[a, f, d, -4]) *
-                                                          conj(Y1[b, c, f, -3]) *
-                                                          X1[a, b, e, -1] * X2[e, c, d, -2]
-                            if FusionStyle(I) isa MultiplicityFreeFusion
-                                f2 = fill(Fsymbol(a, b, c, d, e, f) * dim(d), (1, 1, 1, 1))
-                            else
-                                f2 = Fsymbol(a, b, c, d, e, f) * dim(d)
-                            end
-                            @test isapprox(f1, f2; atol=1e-12, rtol=1e-12)
-                        end
-                    end
-                end
+    end
+end
+if BraidingStyle(I) isa Bosonic && hasfusiontensor(I)
+    @testset "Sector $I: fusion tensor and F-move and R-move" begin
+        for a in smallset(I), b in smallset(I)
+            for c in ⊗(a, b)
+                X1 = permutedims(fusiontensor(a, b, c), (2, 1, 3, 4))
+                X2 = fusiontensor(b, a, c)
+                l = dim(a) * dim(b) * dim(c)
+                R = LinearAlgebra.transpose(Rsymbol(a, b, c))
+                sz = (l, convert(Int, Nsymbol(a, b, c)))
+                @test reshape(X1, sz) ≈ reshape(X2, sz) * R
             end
         end
-        @testset "Sector $Istr: Unitarity of F-move" begin
-            for a in smallset(I), b in smallset(I), c in smallset(I)
-                for d in ⊗(a, b, c)
-                    es = collect(intersect(⊗(a, b), map(dual, ⊗(c, dual(d)))))
-                    fs = collect(intersect(⊗(b, c), map(dual, ⊗(dual(d), a))))
+        for a in smallset(I), b in smallset(I), c in smallset(I)
+            for e in ⊗(a, b), f in ⊗(b, c)
+                for d in intersect(⊗(e, c), ⊗(a, f))
+                    X1 = fusiontensor(a, b, e)
+                    X2 = fusiontensor(e, c, d)
+                    Y1 = fusiontensor(b, c, f)
+                    Y2 = fusiontensor(a, f, d)
+                    @tensor f1[-1, -2, -3, -4] := conj(Y2[a, f, d, -4]) *
+                                                  conj(Y1[b, c, f, -3]) *
+                                                  X1[a, b, e, -1] * X2[e, c, d, -2]
                     if FusionStyle(I) isa MultiplicityFreeFusion
-                        @test length(es) == length(fs)
-                        F = [Fsymbol(a, b, c, d, e, f) for e in es, f in fs]
+                        f2 = fill(Fsymbol(a, b, c, d, e, f) * dim(d), (1, 1, 1, 1))
                     else
-                        Fblocks = Vector{Any}()
-                        for e in es
-                            for f in fs
-                                Fs = Fsymbol(a, b, c, d, e, f)
-                                push!(Fblocks,
-                                      reshape(Fs,
-                                              (size(Fs, 1) * size(Fs, 2),
-                                               size(Fs, 3) * size(Fs, 4))))
-                            end
-                        end
-                        F = hvcat(length(fs), Fblocks...)
+                        f2 = Fsymbol(a, b, c, d, e, f) * dim(d)
                     end
-                    @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12)
+                    @test isapprox(f1, f2; atol=1e-12, rtol=1e-12)
                 end
             end
         end
-        @testset "Sector $Istr: Pentagon equation" begin
-            for a in smallset(I), b in smallset(I), c in smallset(I), d in smallset(I)
-                @test pentagon_equation(a, b, c, d; atol=1e-12, rtol=1e-12)
-            end
-        end
-        @testset "Sector $Istr: Hexagon equation" begin
-            for a in smallset(I), b in smallset(I), c in smallset(I)
-                @test hexagon_equation(a, b, c; atol=1e-12, rtol=1e-12)
+    end
+end
+@testset "Sector $Istr: Unitarity of F-move" begin
+    for a in smallset(I), b in smallset(I), c in smallset(I)
+        for d in ⊗(a, b, c)
+            es = collect(intersect(⊗(a, b), map(dual, ⊗(c, dual(d)))))
+            fs = collect(intersect(⊗(b, c), map(dual, ⊗(dual(d), a))))
+            if FusionStyle(I) isa MultiplicityFreeFusion
+                @test length(es) == length(fs)
+                F = [Fsymbol(a, b, c, d, e, f) for e in es, f in fs]
+            else
+                Fblocks = Vector{Any}()
+                for e in es
+                    for f in fs
+                        Fs = Fsymbol(a, b, c, d, e, f)
+                        push!(Fblocks,
+                              reshape(Fs,
+                                      (size(Fs, 1) * size(Fs, 2),
+                                       size(Fs, 3) * size(Fs, 4))))
+                    end
+                end
+                F = hvcat(length(fs), Fblocks...)
             end
+            @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12)
         end
     end
-
-    @testset "Deligne product" begin
-        sectorlist′ = (Trivial, sectorlist...)
-        for I1 in sectorlist′, I2 in sectorlist′
-            a = first(smallset(I1))
-            b = first(smallset(I2))
-
-            @constinferred a ⊠ b
-            @constinferred a ⊠ b ⊠ a
-            @constinferred a ⊠ b ⊠ a ⊠ b
-            @constinferred I1 ⊠ I2
-            @test typeof(a ⊠ b) == I1 ⊠ I2
-        end
+end
+@testset "Sector $Istr: Pentagon equation" begin
+    for a in smallset(I), b in smallset(I), c in smallset(I), d in smallset(I)
+        @test pentagon_equation(a, b, c, d; atol=1e-12, rtol=1e-12)
+    end
+end
+@testset "Sector $Istr: Hexagon equation" begin
+    for a in smallset(I), b in smallset(I), c in smallset(I)
+        @test hexagon_equation(a, b, c; atol=1e-12, rtol=1e-12)
     end
 end
diff --git a/TensorKitSectors/test/testsetup.jl b/TensorKitSectors/test/testsetup.jl
new file mode 100644
index 00000000..ea76781c
--- /dev/null
+++ b/TensorKitSectors/test/testsetup.jl
@@ -0,0 +1,55 @@
+"""
+    module TestSetup
+
+Basic utility for testing sectors.
+"""
+module TestSetup
+
+export smallset, randsector, hasfusiontensor
+export @include
+
+using TensorKitSectors, Test, Random
+using TestExtras
+using Base.Iterators: take, product
+
+smallset(::Type{I}) where {I<:Sector} = take(values(I), 5)
+function smallset(::Type{ProductSector{Tuple{I1,I2}}}) where {I1,I2}
+    iter = product(smallset(I1), smallset(I2))
+    s = collect(i ⊠ j for (i, j) in iter if dim(i) * dim(j) <= 6)
+    return length(s) > 6 ? rand(s, 6) : s
+end
+function smallset(::Type{ProductSector{Tuple{I1,I2,I3}}}) where {I1,I2,I3}
+    iter = product(smallset(I1), smallset(I2), smallset(I3))
+    s = collect(i ⊠ j ⊠ k for (i, j, k) in iter if dim(i) * dim(j) * dim(k) <= 6)
+    return length(s) > 6 ? rand(s, 6) : s
+end
+function randsector(::Type{I}) where {I<:Sector}
+    s = collect(smallset(I))
+    a = rand(s)
+    while a == one(a) # don't use trivial label
+        a = rand(s)
+    end
+    return a
+end
+function hasfusiontensor(I::Type{<:Sector})
+    try
+        fusiontensor(one(I), one(I), one(I))
+        return true
+    catch e
+        if e isa MethodError
+            return false
+        else
+            rethrow(e)
+        end
+    end
+end
+
+# include file as-is in local scope
+macro include(filename::String)
+    dir = dirname(string(__source__.file))
+    filepath = joinpath(dir, filename)
+    source = "quote; " * read(filepath, String) * "; end"
+    return esc(Meta.parse(source).args[1])
+end
+
+end # module TestSetup

From 3391897edc46146a2c5ff82814344a5ea3d7ba5e Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 18 Mar 2024 19:08:47 +0100
Subject: [PATCH 09/21] Add precompile statements

---
 TensorKitSectors/src/TensorKitSectors.jl | 12 +++++++++
 TensorKitSectors/src/precompile.jl       | 32 ++++++++++++++++++++++++
 2 files changed, 44 insertions(+)
 create mode 100644 TensorKitSectors/src/precompile.jl

diff --git a/TensorKitSectors/src/TensorKitSectors.jl b/TensorKitSectors/src/TensorKitSectors.jl
index 8ac27194..aae4d1dc 100644
--- a/TensorKitSectors/src/TensorKitSectors.jl
+++ b/TensorKitSectors/src/TensorKitSectors.jl
@@ -52,4 +52,16 @@ include("product.jl")   # direct product of different sectors
 include("fermions.jl")  # irreps with defined fermionparity and fermionic braiding
 include("anyons.jl")    # non-group sectors
 
+# precompile
+# ----------
+include("precompile.jl")
+
+function __precompile__()
+    for I in (Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep,
+              FermionParity, FermionNumber, FermionSpin, PlanarTrivial, FibonacciAnyon,
+              IsingAnyon)
+        precompile_sector(I)
+    end
+end
+
 end # module TensorKitSectors
diff --git a/TensorKitSectors/src/precompile.jl b/TensorKitSectors/src/precompile.jl
new file mode 100644
index 00000000..2c5dc41e
--- /dev/null
+++ b/TensorKitSectors/src/precompile.jl
@@ -0,0 +1,32 @@
+"""
+    precompile_sector(I::Type{<:Sector})
+
+Precompile common methods for the given sector type.
+"""
+function precompile_sector(::Type{I}) where {I<:Sector}
+    precompile(Nsymbol, (I, I, I))
+    precompile(Fsymbol, (I, I, I, I, I, I))
+    precompile(Rsymbol, (I, I, I))
+    precompile(Asymbol, (I, I, I))
+    precompile(Bsymbol, (I, I, I))
+    
+    precompile(⊗, (I,))
+    precompile(⊗, (I, I))
+    precompile(⊗, (I, I, I))
+    precompile(⊗, (I, I, I, I))
+    
+    precompile(FusionStyle, (I,))
+    precompile(BraidingStyle, (I,))
+    
+    precompile(dim, (I,))
+    precompile(sqrtdim, (I,))
+    precompile(isqrtdim, (I,))
+    precompile(dual, (I,))
+    precompile(twist, (I,))
+    precompile(frobeniusschur, (I,))
+    
+    try
+        precompile(fusiontensor, (I, I, I))
+    catch
+    end
+end

From e6a06e8c92fe4fe5694f4ae40497e8699cb2e883 Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 18 Mar 2024 20:33:42 +0100
Subject: [PATCH 10/21] =?UTF-8?q?Add=20non-unicode=20alternative=20`fusion?=
 =?UTF-8?q?product`=20for=20`=E2=8A=97`?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 TensorKitSectors/src/TensorKitSectors.jl | 2 +-
 TensorKitSectors/src/sectors.jl          | 3 +++
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/TensorKitSectors/src/TensorKitSectors.jl b/TensorKitSectors/src/TensorKitSectors.jl
index aae4d1dc..b91c4f78 100644
--- a/TensorKitSectors/src/TensorKitSectors.jl
+++ b/TensorKitSectors/src/TensorKitSectors.jl
@@ -11,7 +11,7 @@ export Irrep
 
 export Nsymbol, Fsymbol, Rsymbol, Asymbol, Bsymbol
 export dim, sqrtdim, isqrtdim, frobeniusschur, twist, fusiontensor, dual
-export deligneproduct
+export fusionproduct, deligneproduct
 export FusionStyle, UniqueFusion, MultipleFusion, SimpleFusion, GenericFusion,
        MultiplicityFreeFusion
 export BraidingStyle, NoBraiding, SymmetricBraiding, Bosonic, Fermionic, Anyonic
diff --git a/TensorKitSectors/src/sectors.jl b/TensorKitSectors/src/sectors.jl
index 1cf2e3e5..bea4f578 100644
--- a/TensorKitSectors/src/sectors.jl
+++ b/TensorKitSectors/src/sectors.jl
@@ -124,6 +124,9 @@ Return an iterable of elements of `c::I` that appear in the fusion product `a 
 Note that every element `c` should appear at most once, fusion degeneracies (if
 `FusionStyle(I) == GenericFusion()`) should be accessed via `Nsymbol(a, b, c)`.
 """
+fusionproduct(::Sector, ::Sector)
+const ⊗ = fusionproduct
+
 ⊗(::Trivial, ::Trivial) = (Trivial(),)
 ⊗(I::Sector) = (I,)
 const otimes = ⊗

From fa62472b8138622a37419407962c960f5e883097 Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 18 Mar 2024 20:41:25 +0100
Subject: [PATCH 11/21] =?UTF-8?q?Add=20`directproduct`=20for=20`=C3=97`?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 TensorKitSectors/src/TensorKitSectors.jl | 2 +-
 TensorKitSectors/src/groups.jl           | 9 +++++++++
 2 files changed, 10 insertions(+), 1 deletion(-)

diff --git a/TensorKitSectors/src/TensorKitSectors.jl b/TensorKitSectors/src/TensorKitSectors.jl
index b91c4f78..9f5b6136 100644
--- a/TensorKitSectors/src/TensorKitSectors.jl
+++ b/TensorKitSectors/src/TensorKitSectors.jl
@@ -11,7 +11,7 @@ export Irrep
 
 export Nsymbol, Fsymbol, Rsymbol, Asymbol, Bsymbol
 export dim, sqrtdim, isqrtdim, frobeniusschur, twist, fusiontensor, dual
-export fusionproduct, deligneproduct
+export fusionproduct, deligneproduct, directproduct
 export FusionStyle, UniqueFusion, MultipleFusion, SimpleFusion, GenericFusion,
        MultiplicityFreeFusion
 export BraidingStyle, NoBraiding, SymmetricBraiding, Bosonic, Fermionic, Anyonic
diff --git a/TensorKitSectors/src/groups.jl b/TensorKitSectors/src/groups.jl
index da5bb7a0..6560babd 100644
--- a/TensorKitSectors/src/groups.jl
+++ b/TensorKitSectors/src/groups.jl
@@ -23,6 +23,15 @@ const GroupTuple = Tuple{Vararg{Group}}
 
 abstract type ProductGroup{T<:GroupTuple} <: Group end
 
+@doc """
+    directproduct(G::Vararg{Type{<:Group}}) -> ProductGroup{Tuple{G...}}
+    ×(G::Vararg{Type{<:Group}}) -> ProductGroup{Tuple{G...}}
+
+Construct the direct product of a (list of) groups.
+"""
+directproduct
+const × = directproduct
+
 ×(a::Type{<:Group}, b::Type{<:Group}, c::Type{<:Group}...) = ×(×(a, b), c...)
 ×(G::Type{<:Group}) = ProductGroup{Tuple{G}}
 ×(G1::Type{ProductGroup{Tuple{}}},

From deefc45909a0bd54076a21b3c32f7c41fb993bca Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 18 Mar 2024 20:54:10 +0100
Subject: [PATCH 12/21] split "trivial.jl"

---
 TensorKitSectors/src/TensorKitSectors.jl |  1 +
 TensorKitSectors/src/sectors.jl          | 45 ------------------------
 TensorKitSectors/src/trivial.jl          | 36 +++++++++++++++++++
 3 files changed, 37 insertions(+), 45 deletions(-)
 create mode 100644 TensorKitSectors/src/trivial.jl

diff --git a/TensorKitSectors/src/TensorKitSectors.jl b/TensorKitSectors/src/TensorKitSectors.jl
index 9f5b6136..c00d9121 100644
--- a/TensorKitSectors/src/TensorKitSectors.jl
+++ b/TensorKitSectors/src/TensorKitSectors.jl
@@ -46,6 +46,7 @@ using WignerSymbols
 # --------
 include("auxiliary.jl")
 include("sectors.jl")
+include("trivial.jl")
 include("groups.jl")
 include("irreps.jl")    # irreps of symmetry groups, with bosonic braiding
 include("product.jl")   # direct product of different sectors
diff --git a/TensorKitSectors/src/sectors.jl b/TensorKitSectors/src/sectors.jl
index bea4f578..d1acbada 100644
--- a/TensorKitSectors/src/sectors.jl
+++ b/TensorKitSectors/src/sectors.jl
@@ -58,25 +58,6 @@ Base.IteratorEltype(::Type{<:SectorValues}) = HasEltype()
 Base.eltype(::Type{SectorValues{I}}) where {I<:Sector} = I
 Base.values(::Type{I}) where {I<:Sector} = SectorValues{I}()
 
-# Define a sector for ungraded vector spaces
-"""
-    Trivial
-
-Singleton type to represent the trivial sector, i.e. the trivial representation of the
-trivial group. This is equivalent to `Rep[ℤ₁]`, or the unit object of the category `Vect` of
-ordinary vector spaces.
-"""
-struct Trivial <: Sector end
-Base.show(io::IO, ::Trivial) = print(io, "Trivial()")
-
-Base.IteratorSize(::Type{SectorValues{Trivial}}) = HasLength()
-Base.length(::SectorValues{Trivial}) = 1
-Base.iterate(::SectorValues{Trivial}, i=false) = return i ? nothing : (Trivial(), true)
-function Base.getindex(::SectorValues{Trivial}, i::Int)
-    return i == 1 ? Trivial() : throw(BoundsError(values(Trivial), i))
-end
-findindex(::SectorValues{Trivial}, c::Trivial) = 1
-
 """
     one(::Sector) -> Sector
     one(::Type{<:Sector}) -> Sector
@@ -84,7 +65,6 @@ findindex(::SectorValues{Trivial}, c::Trivial) = 1
 Return the unit element within this type of sector.
 """
 Base.one(a::Sector) = one(typeof(a))
-Base.one(::Type{Trivial}) = Trivial()
 
 """
     dual(a::Sector) -> Sector
@@ -92,7 +72,6 @@ Base.one(::Type{Trivial}) = Trivial()
 Return the conjugate label `conj(a)`.
 """
 dual(a::Sector) = conj(a)
-Base.conj(::Trivial) = Trivial()
 
 """
     isreal(::Type{<:Sector}) -> Bool
@@ -109,9 +88,6 @@ function Base.isreal(I::Type{<:Sector})
         return (eltype(Fsymbol(u, u, u, u, u, u)) <: Real)
     end
 end
-Base.isreal(::Type{Trivial}) = true
-
-Base.isless(::Trivial, ::Trivial) = false
 
 # FusionStyle: the most important aspect of Sector
 #---------------------------------------------
@@ -138,7 +114,6 @@ Return an `Integer` representing the number of times `c` appears in the fusion p
 `a ⊗ b`. Could be a `Bool` if `FusionStyle(I) == UniqueFusion()` or `SimpleFusion()`.
 """
 function Nsymbol end
-Nsymbol(::Trivial, ::Trivial, ::Trivial) = true
 
 # trait to describe the fusion of superselection sectors
 abstract type FusionStyle end
@@ -166,23 +141,6 @@ There is an abstract supertype `MultipleFusion` of which both `SimpleFusion` and
 `MultiplicityFreeFusion = Union{UniqueFusion,SimpleFusion}`.
 """
 FusionStyle(a::Sector) = FusionStyle(typeof(a))
-FusionStyle(::Type{Trivial}) = UniqueFusion()
-
-# NOTE: the following inline is extremely important for performance, especially
-# in the case of UniqueFusion, because ⊗(...) is computed very often
-@inline function ⊗(a::I, b::I, c::I, rest::Vararg{I}) where {I<:Sector}
-    if FusionStyle(I) isa UniqueFusion
-        return a ⊗ first(⊗(b, c, rest...))
-    else
-        s = Set{I}()
-        for d in ⊗(b, c, rest...)
-            for e in a ⊗ d
-                push!(s, e)
-            end
-        end
-        return s
-    end
-end
 
 """
     Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I<:Sector}
@@ -204,7 +162,6 @@ it is a rank 4 array of size
 `(Nsymbol(a, b, e), Nsymbol(e, c, d), Nsymbol(b, c, f), Nsymbol(a, f, d))`.
 """
 function Fsymbol end
-Fsymbol(::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial) = 1
 
 """
     Rsymbol(a::I, b::I, c::I) where {I<:Sector}
@@ -220,7 +177,6 @@ number. Otherwise it is a square matrix with row and column size
 `Nsymbol(a,b,c) == Nsymbol(b,a,c)`.
 """
 function Rsymbol end
-Rsymbol(::Trivial, ::Trivial, ::Trivial) = 1
 
 # If a I::Sector with `fusion(I) == GenericFusion` fusion wants to have custom vertex
 # labels, a specialized method for `vertindex2label` should be added
@@ -366,7 +322,6 @@ braids are in fact equivalent to crossings (i.e. braiding twice is an identity:
 `isone(Rsymbol(b,a,c)*Rsymbol(a,b,c)) == true`) and permutations are uniquely defined.
 """
 BraidingStyle(a::Sector) = BraidingStyle(typeof(a))
-BraidingStyle(::Type{Trivial}) = Bosonic()
 
 # Pentagon and Hexagon equations
 #-------------------------------------------------------------------------------
diff --git a/TensorKitSectors/src/trivial.jl b/TensorKitSectors/src/trivial.jl
new file mode 100644
index 00000000..be60db31
--- /dev/null
+++ b/TensorKitSectors/src/trivial.jl
@@ -0,0 +1,36 @@
+"""
+    Trivial
+
+Singleton type to represent the trivial sector, i.e. the trivial representation of the
+trivial group. This is equivalent to `Rep[ℤ₁]`, or the unit object of the category `Vect` of
+ordinary vector spaces.
+"""
+struct Trivial <: Sector end
+
+Base.show(io::IO, ::Trivial) = print(io, "Trivial()")
+
+# iterators
+Base.IteratorSize(::Type{SectorValues{Trivial}}) = HasLength()
+Base.length(::SectorValues{Trivial}) = 1
+Base.iterate(::SectorValues{Trivial}, i=false) = return i ? nothing : (Trivial(), true)
+function Base.getindex(::SectorValues{Trivial}, i::Int)
+    return i == 1 ? Trivial() : throw(BoundsError(values(Trivial), i))
+end
+findindex(::SectorValues{Trivial}, c::Trivial) = 1
+
+# basic properties
+Base.one(::Type{Trivial}) = Trivial()
+Base.conj(::Trivial) = Trivial()
+
+Base.isreal(::Type{Trivial}) = true
+Base.isless(::Trivial, ::Trivial) = false
+
+# fusion rules
+⊗(::Trivial, ::Trivial) = (Trivial(),)
+Nsymbol(::Trivial, ::Trivial, ::Trivial) = true
+FusionStyle(::Type{Trivial}) = UniqueFusion()
+Fsymbol(::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial) = 1
+
+# braiding rules
+Rsymbol(::Trivial, ::Trivial, ::Trivial) = 1
+BraidingStyle(::Type{Trivial}) = Bosonic()

From e76eb51b1904de4bdaf259f899867ac9eedea09e Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Wed, 20 Mar 2024 09:07:54 +0100
Subject: [PATCH 13/21] Formatter

---
 TensorKitSectors/src/precompile.jl | 8 ++++----
 TensorKitSectors/test/runtests.jl  | 2 +-
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/TensorKitSectors/src/precompile.jl b/TensorKitSectors/src/precompile.jl
index 2c5dc41e..491595b9 100644
--- a/TensorKitSectors/src/precompile.jl
+++ b/TensorKitSectors/src/precompile.jl
@@ -9,22 +9,22 @@ function precompile_sector(::Type{I}) where {I<:Sector}
     precompile(Rsymbol, (I, I, I))
     precompile(Asymbol, (I, I, I))
     precompile(Bsymbol, (I, I, I))
-    
+
     precompile(⊗, (I,))
     precompile(⊗, (I, I))
     precompile(⊗, (I, I, I))
     precompile(⊗, (I, I, I, I))
-    
+
     precompile(FusionStyle, (I,))
     precompile(BraidingStyle, (I,))
-    
+
     precompile(dim, (I,))
     precompile(sqrtdim, (I,))
     precompile(isqrtdim, (I,))
     precompile(dual, (I,))
     precompile(twist, (I,))
     precompile(frobeniusschur, (I,))
-    
+
     try
         precompile(fusiontensor, (I, I, I))
     catch
diff --git a/TensorKitSectors/test/runtests.jl b/TensorKitSectors/test/runtests.jl
index 1cea3b6e..f32673ea 100644
--- a/TensorKitSectors/test/runtests.jl
+++ b/TensorKitSectors/test/runtests.jl
@@ -21,7 +21,7 @@ const sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewS
                     NewSU2Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ NewSU2Irrep,
                     Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)
 
-@testset verbose=true "$I" for I in sectorlist
+@testset verbose = true "$I" for I in sectorlist
     @include("sectors.jl")
 end
 

From 1ab57750b34459b86170ff929b58ac258ad22845 Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Fri, 19 Apr 2024 09:12:28 +0200
Subject: [PATCH 14/21] Change unicode alternatives to `otimes` and `times`

Incorporate changes from #aada0b2
---
 TensorKitSectors/src/TensorKitSectors.jl |  2 +-
 TensorKitSectors/src/groups.jl           |  8 ++++----
 TensorKitSectors/src/sectors.jl          | 22 ++++++++++++++++++----
 3 files changed, 23 insertions(+), 9 deletions(-)

diff --git a/TensorKitSectors/src/TensorKitSectors.jl b/TensorKitSectors/src/TensorKitSectors.jl
index c00d9121..2a5b456d 100644
--- a/TensorKitSectors/src/TensorKitSectors.jl
+++ b/TensorKitSectors/src/TensorKitSectors.jl
@@ -11,7 +11,7 @@ export Irrep
 
 export Nsymbol, Fsymbol, Rsymbol, Asymbol, Bsymbol
 export dim, sqrtdim, isqrtdim, frobeniusschur, twist, fusiontensor, dual
-export fusionproduct, deligneproduct, directproduct
+export otimes, deligneproduct, times
 export FusionStyle, UniqueFusion, MultipleFusion, SimpleFusion, GenericFusion,
        MultiplicityFreeFusion
 export BraidingStyle, NoBraiding, SymmetricBraiding, Bosonic, Fermionic, Anyonic
diff --git a/TensorKitSectors/src/groups.jl b/TensorKitSectors/src/groups.jl
index 6560babd..2d34e069 100644
--- a/TensorKitSectors/src/groups.jl
+++ b/TensorKitSectors/src/groups.jl
@@ -23,14 +23,14 @@ const GroupTuple = Tuple{Vararg{Group}}
 
 abstract type ProductGroup{T<:GroupTuple} <: Group end
 
-@doc """
-    directproduct(G::Vararg{Type{<:Group}}) -> ProductGroup{Tuple{G...}}
+"""
     ×(G::Vararg{Type{<:Group}}) -> ProductGroup{Tuple{G...}}
+    times(G::Vararg{Type{<:Group}}) -> ProductGroup{Tuple{G...}}
 
 Construct the direct product of a (list of) groups.
 """
-directproduct
-const × = directproduct
+function ×(::Vararg{Type{<:Group}}) end
+const times = ×
 
 ×(a::Type{<:Group}, b::Type{<:Group}, c::Type{<:Group}...) = ×(×(a, b), c...)
 ×(G::Type{<:Group}) = ProductGroup{Tuple{G}}
diff --git a/TensorKitSectors/src/sectors.jl b/TensorKitSectors/src/sectors.jl
index d1acbada..f279beee 100644
--- a/TensorKitSectors/src/sectors.jl
+++ b/TensorKitSectors/src/sectors.jl
@@ -100,12 +100,26 @@ Return an iterable of elements of `c::I` that appear in the fusion product `a 
 Note that every element `c` should appear at most once, fusion degeneracies (if
 `FusionStyle(I) == GenericFusion()`) should be accessed via `Nsymbol(a, b, c)`.
 """
-fusionproduct(::Sector, ::Sector)
-const ⊗ = fusionproduct
+function ⊗(::Vararg{Sector}) end
+const otimes = ⊗
 
-⊗(::Trivial, ::Trivial) = (Trivial(),)
 ⊗(I::Sector) = (I,)
-const otimes = ⊗
+
+# NOTE: the following inline is extremely important for performance, especially
+# in the case of UniqueFusion, because ⊗(...) is computed very often
+@inline function ⊗(a::I, b::I, c::I, rest::Vararg{I}) where {I<:Sector}
+    if FusionStyle(I) isa UniqueFusion
+        return a ⊗ first(⊗(b, c, rest...))
+    else
+        s = Set{I}()
+        for d in ⊗(b, c, rest...)
+            for e in a ⊗ d
+                push!(s, e)
+            end
+        end
+        return s
+    end
+end
 
 """
     Nsymbol(a::I, b::I, c::I) where {I<:Sector} -> Integer

From 8a0d082250484a37eab2109c0f39717da31c4238 Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Fri, 19 Apr 2024 09:15:49 +0200
Subject: [PATCH 15/21] Update some testset descriptions

---
 TensorKitSectors/test/runtests.jl | 2 +-
 TensorKitSectors/test/sectors.jl  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/TensorKitSectors/test/runtests.jl b/TensorKitSectors/test/runtests.jl
index f32673ea..3e1ddc7a 100644
--- a/TensorKitSectors/test/runtests.jl
+++ b/TensorKitSectors/test/runtests.jl
@@ -21,7 +21,7 @@ const sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewS
                     NewSU2Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ NewSU2Irrep,
                     Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)
 
-@testset verbose = true "$I" for I in sectorlist
+@testset "$(TensorKitSectors.type_repr(I))" for I in sectorlist
     @include("sectors.jl")
 end
 
diff --git a/TensorKitSectors/test/sectors.jl b/TensorKitSectors/test/sectors.jl
index 582092c3..53d20ee7 100644
--- a/TensorKitSectors/test/sectors.jl
+++ b/TensorKitSectors/test/sectors.jl
@@ -42,7 +42,7 @@ end
     end
 end
 if BraidingStyle(I) isa Bosonic && hasfusiontensor(I)
-    @testset "Sector $I: fusion tensor and F-move and R-move" begin
+    @testset "Sector $Istr: fusion tensor and F-move and R-move" begin
         for a in smallset(I), b in smallset(I)
             for c in ⊗(a, b)
                 X1 = permutedims(fusiontensor(a, b, c), (2, 1, 3, 4))

From b1090c322684c74c0834a6cd3e641508fe90c5ab Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 22 Apr 2024 11:53:34 +0200
Subject: [PATCH 16/21] Resolve method ambiguities

---
 TensorKitSectors/src/sectors.jl | 2 +-
 src/spaces/productspace.jl      | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/TensorKitSectors/src/sectors.jl b/TensorKitSectors/src/sectors.jl
index f279beee..f54d7364 100644
--- a/TensorKitSectors/src/sectors.jl
+++ b/TensorKitSectors/src/sectors.jl
@@ -100,7 +100,7 @@ Return an iterable of elements of `c::I` that appear in the fusion product `a 
 Note that every element `c` should appear at most once, fusion degeneracies (if
 `FusionStyle(I) == GenericFusion()`) should be accessed via `Nsymbol(a, b, c)`.
 """
-function ⊗(::Vararg{Sector}) end
+function ⊗ end
 const otimes = ⊗
 
 ⊗(I::Sector) = (I,)
diff --git a/src/spaces/productspace.jl b/src/spaces/productspace.jl
index 07b08076..e27be39b 100644
--- a/src/spaces/productspace.jl
+++ b/src/spaces/productspace.jl
@@ -208,7 +208,7 @@ Base.hash(P::ProductSpace, h::UInt) = hash(P.spaces, h)
 
 # Default construction from product of spaces
 #---------------------------------------------
-⊗(V::ElementarySpace...) = ProductSpace(V...)
+⊗(V::ElementarySpace, Vrest::ElementarySpace...) = ProductSpace(V, Vrest...)
 ⊗(P::ProductSpace) = P
 function ⊗(P1::ProductSpace{S}, P2::ProductSpace{S}) where {S<:ElementarySpace}
     return ProductSpace{S}(tuple(P1.spaces..., P2.spaces...))

From 620f180a1ee2b452943ac8af0a0065cdc5be9218 Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 22 Apr 2024 11:56:20 +0200
Subject: [PATCH 17/21] Change some formatting

---
 TensorKitSectors/src/sectors.jl | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/TensorKitSectors/src/sectors.jl b/TensorKitSectors/src/sectors.jl
index f54d7364..527c1856 100644
--- a/TensorKitSectors/src/sectors.jl
+++ b/TensorKitSectors/src/sectors.jl
@@ -131,13 +131,10 @@ function Nsymbol end
 
 # trait to describe the fusion of superselection sectors
 abstract type FusionStyle end
-struct UniqueFusion <: FusionStyle # unique fusion output when fusing two sectors
-end
+struct UniqueFusion <: FusionStyle end # unique fusion output when fusing two sectors
 abstract type MultipleFusion <: FusionStyle end
-struct SimpleFusion <: MultipleFusion # multiple fusion but multiplicity free
-end
-struct GenericFusion <: MultipleFusion # multiple fusion with multiplicities
-end
+struct SimpleFusion <: MultipleFusion end # multiple fusion but multiplicity free
+struct GenericFusion <: MultipleFusion end # multiple fusion with multiplicities
 const MultiplicityFreeFusion = Union{UniqueFusion,SimpleFusion}
 
 """

From d4964513e394f728a1aec45a560d109e0a7bcb29 Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 22 Apr 2024 11:59:25 +0200
Subject: [PATCH 18/21] Import local version of `TensorKitSectors`

---
 test/runtests.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 9705fc57..eccb7f93 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,3 +1,7 @@
+# TODO: remove this once separate folders are registered
+import Pkg;
+Pkg.develop(Pkg.PackageSpec(path=joinpath(@__DIR__, "..", "TensorKitSectors")))
+
 using Test
 using TestExtras
 using Random

From 72f552f59aee10b1fdd53e071ba51d9c9cae701e Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 22 Apr 2024 12:56:58 +0200
Subject: [PATCH 19/21] Enable TensorKitSectors github action

---
 .github/workflows/CI-Sectors.yml | 49 ++++++++++++++++++++++++++++++++
 1 file changed, 49 insertions(+)
 create mode 100644 .github/workflows/CI-Sectors.yml

diff --git a/.github/workflows/CI-Sectors.yml b/.github/workflows/CI-Sectors.yml
new file mode 100644
index 00000000..b7f8be6c
--- /dev/null
+++ b/.github/workflows/CI-Sectors.yml
@@ -0,0 +1,49 @@
+name: CI-sectors
+on:
+  push:
+    branches:
+      - 'master'
+      - 'main'
+      - 'release-'
+    tags: '*'
+  pull_request:
+  workflow_dispatch:
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  test:
+    name: TensorKitSectors - ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1' # automatically expands to the latest stable 1.x release of Julia
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v1
+      - uses: julia-actions/julia-buildpkg@latest
+        with:
+          project: TensorKitSectors
+      - uses: julia-actions/julia-runtest@latest
+        with:
+          project: TensorKitSectors
+        env:
+          JULIA_NUM_THREADS: 4
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v4
+        env:
+          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
+        with:
+          file: lcov.info
\ No newline at end of file

From f88771e06291ece0432fe41acf42954b23de6a2b Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Mon, 22 Apr 2024 13:28:18 +0200
Subject: [PATCH 20/21] Switch to local include for CI

---
 Project.toml       | 1 -
 src/TensorKit.jl   | 9 +++++----
 test/newsectors.jl | 7 ++++---
 test/runtests.jl   | 4 ----
 4 files changed, 9 insertions(+), 12 deletions(-)

diff --git a/Project.toml b/Project.toml
index 523731dd..7edb6cd5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,7 +9,6 @@ LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
-TensorKitSectors = "7ea117e4-0799-4282-ac6e-748b4bc2d0f5"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
diff --git a/src/TensorKit.jl b/src/TensorKit.jl
index a7756643..e640f065 100644
--- a/src/TensorKit.jl
+++ b/src/TensorKit.jl
@@ -168,10 +168,11 @@ Base.show(io::IO, e::IndexError) = print(io, "IndexError(", e.message, ")")
 # Definitions and methods for superselection sectors (quantum numbers)
 #----------------------------------------------------------------------
 
-using TensorKitSectors
-import TensorKitSectors: dim, BraidingStyle, FusionStyle, ⊠, ⊗
-import TensorKitSectors: dual, type_repr
-import TensorKitSectors: twist
+include(joinpath(@__DIR__, "..", "TensorKitSectors", "src", "TensorKitSectors.jl"))
+using .TensorKitSectors
+import .TensorKitSectors: dim, BraidingStyle, FusionStyle, ⊠, ⊗
+import .TensorKitSectors: dual, type_repr
+import .TensorKitSectors: twist
 
 # Constructing and manipulating fusion trees and iterators thereof
 #------------------------------------------------------------------
diff --git a/test/newsectors.jl b/test/newsectors.jl
index 16cd8ef0..ca625aa1 100644
--- a/test/newsectors.jl
+++ b/test/newsectors.jl
@@ -6,10 +6,11 @@ module NewSectors
 export NewSU2Irrep
 
 using HalfIntegers, WignerSymbols
-using TensorKitSectors
+using TensorKit
+using TensorKit: TensorKitSectors
 
-import TensorKitSectors: FusionStyle, BraidingStyle, Nsymbol, Fsymbol, Rsymbol, dim,
-                         fusiontensor, ⊗
+import TensorKit.TensorKitSectors: FusionStyle, BraidingStyle, Nsymbol, Fsymbol, Rsymbol, dim,
+                         fusiontensor, ⊗, SectorValues
 
 struct NewSU2Irrep <: Sector
     j::HalfInt
diff --git a/test/runtests.jl b/test/runtests.jl
index eccb7f93..9705fc57 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,3 @@
-# TODO: remove this once separate folders are registered
-import Pkg;
-Pkg.develop(Pkg.PackageSpec(path=joinpath(@__DIR__, "..", "TensorKitSectors")))
-
 using Test
 using TestExtras
 using Random

From 78620eb0428a18ad97877670218cf09d3b8fabdc Mon Sep 17 00:00:00 2001
From: lkdvos <lukas.devos@ugent.be>
Date: Tue, 23 Apr 2024 08:19:33 +0200
Subject: [PATCH 21/21] Remove `NewSU2Sector` tests from main package

---
 test/newsectors.jl | 80 ----------------------------------------------
 test/runtests.jl   | 13 +++-----
 2 files changed, 5 insertions(+), 88 deletions(-)
 delete mode 100644 test/newsectors.jl

diff --git a/test/newsectors.jl b/test/newsectors.jl
deleted file mode 100644
index ca625aa1..00000000
--- a/test/newsectors.jl
+++ /dev/null
@@ -1,80 +0,0 @@
-# NewSU2Irrep
-# Test of a bare-bones sector implementation, which is set to be `DegenerateNonAbelian`
-# (even though it is not)
-module NewSectors
-
-export NewSU2Irrep
-
-using HalfIntegers, WignerSymbols
-using TensorKit
-using TensorKit: TensorKitSectors
-
-import TensorKit.TensorKitSectors: FusionStyle, BraidingStyle, Nsymbol, Fsymbol, Rsymbol, dim,
-                         fusiontensor, ⊗, SectorValues
-
-struct NewSU2Irrep <: Sector
-    j::HalfInt
-    function NewSU2Irrep(j)
-        j >= zero(j) || error("Not a valid SU₂ irrep")
-        return new(j)
-    end
-end
-Base.convert(::Type{NewSU2Irrep}, j::Real) = NewSU2Irrep(j)
-
-const _su2one = NewSU2Irrep(zero(HalfInt))
-Base.one(::Type{NewSU2Irrep}) = _su2one
-Base.conj(s::NewSU2Irrep) = s
-function ⊗(s1::NewSU2Irrep, s2::NewSU2Irrep)
-    return SectorSet{NewSU2Irrep}(abs(s1.j - s2.j):(s1.j + s2.j))
-end
-
-Base.IteratorSize(::Type{SectorValues{NewSU2Irrep}}) = Base.IsInfinite()
-Base.iterate(::SectorValues{NewSU2Irrep}, i=0) = (NewSU2Irrep(half(i)), i + 1)
-# Base.getindex(::SectorValues{NewSU2Irrep}, i::Int) =
-#     1 <= i ? NewSU2Irrep(half(i-1)) : throw(BoundsError(values(NewSU2Irrep), i))
-# findindex(::SectorValues{NewSU2Irrep}, s::NewSU2Irrep) = twice(s.j)+1
-
-# TensorKit.dim(s::NewSU2Irrep) = twice(s.j)+1
-#
-FusionStyle(::Type{NewSU2Irrep}) = GenericFusion()
-BraidingStyle(::Type{NewSU2Irrep}) = Bosonic()
-Base.isreal(::Type{NewSU2Irrep}) = true
-
-function Nsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
-    return convert(Int, WignerSymbols.δ(sa.j, sb.j, sc.j))
-end
-
-function Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep,
-                 s4::NewSU2Irrep, s5::NewSU2Irrep, s6::NewSU2Irrep)
-    n1 = Nsymbol(s1, s2, s5)
-    n2 = Nsymbol(s5, s3, s4)
-    n3 = Nsymbol(s2, s3, s6)
-    n4 = Nsymbol(s1, s6, s4)
-    f = all(==(_su2one), (s1, s2, s3, s4, s5, s6)) ? 1.0 :
-        sqrt(dim(s5) * dim(s6)) * WignerSymbols.racahW(Float64, s1.j, s2.j,
-                                                       s4.j, s3.j, s5.j, s6.j)
-    return fill(f, (n1, n2, n3, n4))
-end
-
-function Rsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep)
-    Nsymbol(sa, sb, sc) > 0 || return fill(0.0, (0, 0))
-    return fill(iseven(convert(Int, sa.j + sb.j - sc.j)) ? 1.0 : -1.0, (1, 1))
-end
-
-dim(s::NewSU2Irrep) = twice(s.j) + 1
-
-function fusiontensor(a::NewSU2Irrep, b::NewSU2Irrep, c::NewSU2Irrep)
-    C = Array{Float64}(undef, dim(a), dim(b), dim(c), 1)
-    ja, jb, jc = a.j, b.j, c.j
-
-    for kc in 1:dim(c), kb in 1:dim(b), ka in 1:dim(a)
-        C[ka, kb, kc, 1] = WignerSymbols.clebschgordan(ja, ja + 1 - ka, jb, jb + 1 - kb, jc,
-                                                       jc + 1 - kc)
-    end
-    return C
-end
-
-Base.hash(s::NewSU2Irrep, h::UInt) = hash(s.j, h)
-Base.isless(s1::NewSU2Irrep, s2::NewSU2Irrep) = isless(s1.j, s2.j)
-
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index 9705fc57..34037d27 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,9 +10,6 @@ using Base.Iterators: take, product
 # const SU3Irrep = SUNIrrep{3}
 using LinearAlgebra: LinearAlgebra
 
-include("newsectors.jl")
-using .NewSectors
-
 const TK = TensorKit
 
 Random.seed!(1234)
@@ -49,11 +46,11 @@ function hasfusiontensor(I::Type{<:Sector})
     end
 end
 
-sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irrep,
-              FibonacciAnyon, IsingAnyon, FermionParity, FermionParity ⊠ FermionParity,
-              Z3Irrep ⊠ Z4Irrep, FermionParity ⊠ U1Irrep ⊠ SU2Irrep,
-              FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep,
-              NewSU2Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ NewSU2Irrep,
+sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, Z3Irrep ⊠ Z4Irrep,
+              U1Irrep, CU1Irrep, SU2Irrep,
+              FermionParity, FermionParity ⊠ FermionParity,
+              FermionParity ⊠ U1Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, # Hubbard-like
+              FibonacciAnyon, IsingAnyon,
               Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)
 
 Ti = time()