diff --git a/Project.toml b/Project.toml
index 137aa050c..5f3ffb4f2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -43,7 +43,7 @@ Random = "1"
 SafeTestsets = "0.1"
 ScopedValues = "1.3.0"
 Strided = "2"
-TensorKitSectors = "=0.3.0, 0.3.2"
+TensorKitSectors = "0.3.3"
 TensorOperations = "5.1"
 Test = "1"
 TestExtras = "0.2,0.3"
diff --git a/src/TensorKit.jl b/src/TensorKit.jl
index a7342d9fb..32e7ae577 100644
--- a/src/TensorKit.jl
+++ b/src/TensorKit.jl
@@ -11,11 +11,11 @@ module TensorKit
 export Sector, AbstractIrrep, Irrep
 export FusionStyle, UniqueFusion, MultipleFusion, MultiplicityFreeFusion, SimpleFusion, GenericFusion
 export UnitStyle, SimpleUnit, GenericUnit
-export BraidingStyle, SymmetricBraiding, Bosonic, Fermionic, Anyonic, NoBraiding
+export BraidingStyle, SymmetricBraiding, Bosonic, Fermionic, Anyonic, NoBraiding, HasBraiding
 export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep
 export ProductSector
 export FermionParity, FermionNumber, FermionSpin
-export FibonacciAnyon, IsingAnyon
+export FibonacciAnyon, IsingAnyon, IsingBimodule
 
 # Export common vector space, fusion tree and tensor types
 export VectorSpace, Field, ElementarySpace # abstract vector spaces
@@ -34,7 +34,10 @@ export SpaceMismatch, SectorMismatch, IndexError # error types
 # Export general vector space methods
 export space, field, dual, dim, reduceddim, dims, fuse, flip, isdual
 export unitspace, zerospace, oplus, ominus
+export leftunitspace, rightunitspace, isunitspace
 export insertleftunit, insertrightunit, removeunit
+
+# partial order for vector spaces
 export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic
 
 # Reexport methods for sectors and properties thereof
diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl
index 4b97603ca..a7105cda6 100644
--- a/src/auxiliary/auxiliary.jl
+++ b/src/auxiliary/auxiliary.jl
@@ -80,3 +80,9 @@ function _copyto!(A::StridedView{<:Any, 1}, B::StridedView{<:Any, 2})
 
     return A
 end
+
+@static if VERSION < v"1.11" # TODO: remove once support for v1.10 is dropped
+    _allequal(f, xs) = allequal(Base.Generator(f, xs))
+else
+    _allequal(f, xs) = allequal(f, xs)
+end
diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl
index 1665cdb27..0618590a6 100644
--- a/src/fusiontrees/fusiontrees.jl
+++ b/src/fusiontrees/fusiontrees.jl
@@ -92,7 +92,7 @@ function FusionTree{I}(
         uncoupled::NTuple{N}, coupled = unit(I), isdual = ntuple(Returns(false), N)
     ) where {I <: Sector, N}
     FusionStyle(I) isa UniqueFusion ||
-        error("fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`")
+        throw(ArgumentError("fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`"))
     return FusionTree{I}(
         map(s -> convert(I, s), uncoupled), convert(I, coupled), isdual,
         _abelianinner(map(s -> convert(I, s), (uncoupled..., dual(coupled))))
diff --git a/src/fusiontrees/iterator.jl b/src/fusiontrees/iterator.jl
index 50f0781f6..341929420 100644
--- a/src/fusiontrees/iterator.jl
+++ b/src/fusiontrees/iterator.jl
@@ -17,6 +17,9 @@ function fusiontrees(uncoupled::Tuple{Vararg{I}}, coupled::I) where {I <: Sector
     return fusiontrees(uncoupled, coupled, isdual)
 end
 function fusiontrees(uncoupled::Tuple{I, Vararg{I}}) where {I <: Sector}
+    UnitStyle(I) isa GenericUnit && throw(
+        ArgumentError("Must specify coupled sector when UnitStyle is GenericUnit.")
+    )
     coupled = unit(I)
     isdual = ntuple(n -> false, length(uncoupled))
     return fusiontrees(uncoupled, coupled, isdual)
diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index 7df490b30..d3be67a5f 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -90,10 +90,8 @@ field(::Type{<:GradedSpace}) = ℂ
 InnerProductStyle(::Type{<:GradedSpace}) = EuclideanInnerProduct()
 
 function dim(V::GradedSpace)
-    return reduce(
-        +, dim(V, c) * dim(c) for c in sectors(V);
-        init = zero(dim(unit(sectortype(V))))
-    )
+    init = 0 * dim(first(allunits(sectortype(V))))
+    return sum(c -> dim(c) * dim(V, c), sectors(V); init = init)
 end
 function dim(V::GradedSpace{I, <:AbstractDict}, c::I) where {I <: Sector}
     return get(V.dims, isdual(V) ? dual(c) : c, 0)
@@ -123,8 +121,11 @@ function flip(V::GradedSpace{I}) where {I <: Sector}
     end
 end
 
-unitspace(S::Type{<:GradedSpace{I}}) where {I <: Sector} = S(unit(I) => 1)
-zerospace(S::Type{<:GradedSpace{I}}) where {I <: Sector} = S(unit(I) => 0)
+function unitspace(S::Type{<:GradedSpace{I}}) where {I <: Sector}
+    return S(unit => 1 for unit in allunits(I))
+end
+zerospace(S::Type{<:GradedSpace}) = S()
+
 # TODO: the following methods can probably be implemented more efficiently for
 # `FiniteGradedSpace`, but we don't expect them to be used often in hot loops, so
 # these generic definitions (which are still quite efficient) are good for now.
diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl
index b4126401a..d032a1586 100644
--- a/src/spaces/homspace.jl
+++ b/src/spaces/homspace.jl
@@ -93,11 +93,16 @@ function blocksectors(W::HomSpace)
     N₁ = length(codom)
     N₂ = length(dom)
     I = sectortype(W)
-    # TODO: is sort! still necessary now that blocksectors of ProductSpace is sorted?
-    if N₂ <= N₁
-        return sort!(filter!(c -> hasblock(codom, c), blocksectors(dom)))
+    if N₁ == N₂ == 0
+        return allunits(I)
+    elseif N₁ == 0
+        return filter!(isunit, collect(blocksectors(dom))) # module space cannot end in empty space
+    elseif N₂ == 0
+        return filter!(isunit, collect(blocksectors(codom)))
+    elseif N₂ <= N₁
+        return filter!(c -> hasblock(codom, c), collect(blocksectors(dom)))
     else
-        return sort!(filter!(c -> hasblock(dom, c), blocksectors(codom)))
+        return filter!(c -> hasblock(dom, c), collect(blocksectors(codom)))
     end
 end
 
@@ -227,7 +232,7 @@ function TensorOperations.tensorcontract(
 end
 
 """
-    insertleftunit(W::HomSpace, i=numind(W) + 1; conj=false, dual=false)
+    insertleftunit(W::HomSpace, i = numind(W) + 1; conj = false, dual = false)
 
 Insert a trivial vector space, isomorphic to the underlying field, at position `i`,
 which can be specified as an `Int` or as `Val(i)` for improved type stability.
@@ -237,7 +242,8 @@ See also [`insertrightunit`](@ref insertrightunit(::HomSpace, ::Val{i}) where {i
 [`removeunit`](@ref removeunit(::HomSpace, ::Val{i}) where {i}).
 """
 function insertleftunit(
-        W::HomSpace, ::Val{i} = Val(numind(W) + 1); conj::Bool = false, dual::Bool = false
+        W::HomSpace, ::Val{i} = Val(numind(W) + 1);
+        conj::Bool = false, dual::Bool = false
     ) where {i}
     if i ≤ numout(W)
         return insertleftunit(codomain(W), Val(i); conj, dual) ← domain(W)
@@ -247,7 +253,7 @@ function insertleftunit(
 end
 
 """
-    insertrightunit(W::HomSpace, i=numind(W); conj=false, dual=false)
+    insertrightunit(W::HomSpace, i = numind(W); conj = false, dual = false)
 
 Insert a trivial vector space, isomorphic to the underlying field, after position `i`,
 which can be specified as an `Int` or as `Val(i)` for improved type stability.
@@ -257,7 +263,8 @@ See also [`insertleftunit`](@ref insertleftunit(::HomSpace, ::Val{i}) where {i})
 [`removeunit`](@ref removeunit(::HomSpace, ::Val{i}) where {i}).
 """
 function insertrightunit(
-        W::HomSpace, ::Val{i} = Val(numind(W)); conj::Bool = false, dual::Bool = false
+        W::HomSpace, ::Val{i} = Val(numind(W));
+        conj::Bool = false, dual::Bool = false
     ) where {i}
     if i ≤ numout(W)
         return insertrightunit(codomain(W), Val(i); conj, dual) ← domain(W)
diff --git a/src/spaces/productspace.jl b/src/spaces/productspace.jl
index a0e149764..4b75ce737 100644
--- a/src/spaces/productspace.jl
+++ b/src/spaces/productspace.jl
@@ -1,6 +1,6 @@
 """
-    struct ProductSpace{S<:ElementarySpace, N} <: CompositeSpace{S}
-    ProductSpace(spaces::NTuple{N, S}) where {S<:ElementarySpace, N}
+    struct ProductSpace{S <: ElementarySpace, N} <: CompositeSpace{S}
+    ProductSpace(spaces::NTuple{N, S}) where {S <: ElementarySpace, N}
 
 A `ProductSpace` is a tensor product space of `N` vector spaces of type `S <: ElementarySpace`.
 Only tensor products between [`ElementarySpace`](@ref) objects of the same type are allowed.
@@ -87,7 +87,7 @@ end
 
 # more specific methods
 """
-    sectors(P::ProductSpace{S, N}) where {S<:ElementarySpace}
+    sectors(P::ProductSpace{S, N}) where {S <: ElementarySpace}
 
 Return an iterator over all possible combinations of sectors (represented as an
 `NTuple{N, sectortype(S)}`) that can appear within the tensor product space `P`.
@@ -151,7 +151,7 @@ function blocksectors(P::ProductSpace{S, N}) where {S, N}
     end
     bs = Vector{I}()
     if N == 0
-        push!(bs, unit(I))
+        append!(bs, allunits(I))
     elseif N == 1
         for s in sectors(P)
             push!(bs, first(s))
@@ -196,7 +196,7 @@ hasblock(P::ProductSpace, c::Sector) = !isempty(fusiontrees(P, c))
     blockdim(P::ProductSpace, c::Sector)
 
 Return the total dimension of a coupled sector `c` in the product space, by summing over
-all `dim(P, s)` for all tuples of sectors `s::NTuple{N, <:Sector}` that can fuse to  `c`,
+all `dim(P, s)` for all tuples of sectors `s::NTuple{N, Sector}` that can fuse to  `c`,
 counted with the correct multiplicity (i.e. number of ways in which `s` can fuse to `c`).
 
 See also [`hasblock`](@ref) and [`blocksectors`](@ref).
@@ -228,8 +228,8 @@ end
 
 # unit element with respect to the monoidal structure of taking tensor products
 """
-    one(::S) where {S<:ElementarySpace} -> ProductSpace{S, 0}
-    one(::ProductSpace{S}) where {S<:ElementarySpace} -> ProductSpace{S, 0}
+    one(::S) where {S <: ElementarySpace} -> ProductSpace{S, 0}
+    one(::ProductSpace{S}) where {S <: ElementarySpace} -> ProductSpace{S, 0}
 
 Return a tensor product of zero spaces of type `S`, i.e. this is the unit object under the
 tensor product operation, such that `V ⊗ one(V) == V`.
@@ -248,7 +248,7 @@ fuse(P::ProductSpace{S, 0}) where {S <: ElementarySpace} = unitspace(S)
 fuse(P::ProductSpace{S}) where {S <: ElementarySpace} = fuse(P.spaces...)
 
 """
-    insertleftunit(P::ProductSpace, i::Int=length(P) + 1; conj=false, dual=false)
+    insertleftunit(P::ProductSpace, i::Int = length(P) + 1; conj = false, dual = false)
 
 Insert a trivial vector space, isomorphic to the underlying field, at position `i`,
 which can be specified as an `Int` or as `Val(i)` for improved type stability.
@@ -260,7 +260,18 @@ function insertleftunit(
         P::ProductSpace, ::Val{i} = Val(length(P) + 1);
         conj::Bool = false, dual::Bool = false
     ) where {i}
-    u = unitspace(spacetype(P))
+    N = length(P)
+    I = sectortype(P)
+    if UnitStyle(I) isa SimpleUnit
+        u = unitspace(spacetype(P))
+    else
+        N > 0 || throw(ArgumentError("cannot insert a sensible unit space in the empty product space"))
+        if i == N + 1
+            u = rightunitspace(P[N])
+        else
+            u = leftunitspace(P[i])
+        end
+    end
     if dual
         u = TensorKit.dual(u)
     end
@@ -271,7 +282,7 @@ function insertleftunit(
 end
 
 """
-    insertrightunit(P::ProductSpace, i=lenght(P); conj=false, dual=false)
+    insertrightunit(P::ProductSpace, i = length(P); conj = false, dual = false)
 
 Insert a trivial vector space, isomorphic to the underlying field, after position `i`,
 which can be specified as an `Int` or as `Val(i)` for improved type stability.
@@ -283,7 +294,14 @@ function insertrightunit(
         P::ProductSpace, ::Val{i} = Val(length(P));
         conj::Bool = false, dual::Bool = false
     ) where {i}
-    u = unitspace(spacetype(P))
+    N = length(P)
+    I = sectortype(P)
+    if UnitStyle(I) isa SimpleUnit
+        u = unitspace(spacetype(P))
+    else
+        N > 0 || throw(ArgumentError("cannot insert a sensible unit space in the empty product space"))
+        u = rightunitspace(P[i])
+    end
     if dual
         u = TensorKit.dual(u)
     end
@@ -305,7 +323,8 @@ and [`insertrightunit`](@ref insertrightunit(::ProductSpace, ::Val{i}) where {i}
 """
 function removeunit(P::ProductSpace, ::Val{i}) where {i}
     1 ≤ i ≤ length(P) || _boundserror(P, i)
-    isisomorphic(P[i], unitspace(P[i])) || _nontrivialspaceerror(P, i)
+    I = sectortype(P)
+    isunitspace(P[i]) || _nontrivialspaceerror(P, i)
     return ProductSpace{spacetype(P)}(TupleTools.deleteat(P.spaces, i))
 end
 
diff --git a/src/spaces/vectorspaces.jl b/src/spaces/vectorspaces.jl
index c1d8448fb..5ef171c31 100644
--- a/src/spaces/vectorspaces.jl
+++ b/src/spaces/vectorspaces.jl
@@ -41,8 +41,8 @@ generally, objects in linear monoidal categories.
 abstract type VectorSpace end
 
 """
-    field(a) -> Type{𝔽<:Field}
-    field(::Type{T}) -> Type{𝔽<:Field}
+    field(a) -> Type{𝔽 <: Field}
+    field(::Type{T}) -> Type{𝔽 <: Field}
 
 Return the type of field over which object `a` (e.g. a vector space or a tensor) is defined.
 This also works in type domain.
@@ -130,10 +130,13 @@ Always returns `false` for spaces where `V == conj(V)`, i.e. vector spaces over
 """ isconj(::ElementarySpace)
 
 """
-    unitspace(V::S) where {S<:ElementarySpace} -> S
+    unitspace(V::S) where {S <: ElementarySpace} -> S
 
 Return the corresponding vector space of type `S` that represents the trivial
 one-dimensional space, i.e. the space that is isomorphic to the corresponding field.
+For vector spaces where `I = sectortype(S)` has a semi-simple unit structure
+(`UnitStyle(I) == GenericUnit()`), this returns a multi-dimensional space corresponding to all unit sectors:
+`dim(unitspace(V), s) == 1` for all `s in allunits(I)`. 
 
 !!! note
     `unitspace(V)`is different from `one(V)`. The latter returns the empty product space
@@ -144,7 +147,7 @@ Base.oneunit(V::ElementarySpace) = unitspace(V)
 Base.oneunit(::Type{V}) where {V <: ElementarySpace} = unitspace(V)
 
 """
-    zerospace(V::S) where {S<:ElementarySpace} -> S
+    zerospace(V::S) where {S <: ElementarySpace} -> S
 
 Return the corresponding vector space of type `S` that represents the zero-dimensional or empty space.
 This is the zero element of the direct sum of vector spaces.
@@ -154,6 +157,68 @@ zerospace(V::ElementarySpace) = zerospace(typeof(V))
 Base.zero(V::ElementarySpace) = zerospace(V)
 Base.zero(::Type{V}) where {V <: ElementarySpace} = zerospace(V)
 
+"""
+    leftunitspace(V::S) where {S <: ElementarySpace} -> S
+
+Return the corresponding vector space of type `S` that represents the trivial
+one-dimensional space, i.e. the space that is isomorphic to the corresponding field. For vector spaces 
+of type `GradedSpace{I}`, this one-dimensional space contains the unique left unit of the objects in `Sector` `I` present
+in the vector space.
+"""
+function leftunitspace(V::ElementarySpace)
+    I = sectortype(V)
+    if UnitStyle(I) isa SimpleUnit
+        return unitspace(typeof(V))
+    else
+        !isempty(sectors(V)) || throw(ArgumentError("Cannot determine the left unit of an empty space"))
+        _allequal(leftunit, sectors(V)) ||
+            throw(ArgumentError("sectors of $V do not have the same left unit"))
+
+        sector = leftunit(first(sectors(V)))
+        return spacetype(V)(sector => 1)
+    end
+end
+
+"""
+    rightunitspace(V::S) where {S <: ElementarySpace} -> S
+
+Return the corresponding vector space of type `ElementarySpace` that represents the trivial
+one-dimensional space, i.e. the space that is isomorphic to the corresponding field. For vector spaces 
+of type `GradedSpace{I}`, this corresponds to the right unit of the objects in `Sector` `I` present
+in the vector space.
+"""
+function rightunitspace(V::ElementarySpace)
+    I = sectortype(V)
+    if UnitStyle(I) isa SimpleUnit
+        return unitspace(typeof(V))
+    else
+        !isempty(sectors(V)) || throw(ArgumentError("Cannot determine the right unit of an empty space"))
+        _allequal(rightunit, sectors(V)) ||
+            throw(ArgumentError("sectors of $V do not have the same right unit"))
+
+        sector = rightunit(first(sectors(V)))
+        return spacetype(V)(sector => 1)
+    end
+end
+
+"""
+    isunitspace(V::S) where {S <: ElementarySpace} -> Bool
+
+Return whether the elementary space `V` is a unit space, i.e. is isomorphic to the
+trivial one-dimensional space. For vector spaces of type `GradedSpace{I}` where `Sector` `I` has a
+semi-simple unit structure, this returns `true` if `V` is isomorphic to either the left, right or
+semi-simple unit space.
+"""
+function isunitspace(V::ElementarySpace)
+    I = sectortype(V)
+    return if isa(UnitStyle(I), SimpleUnit)
+        isisomorphic(V, unitspace(V))
+    else
+        (dim(V) == 0 || !all(isunit, sectors(V))) && return false
+        return true
+    end
+end
+
 """
     ⊕(V₁::S, V₂::S, V₃::S...) where {S<:ElementarySpace} -> S
     oplus(V₁::S, V₂::S, V₃::S...) where {S<:ElementarySpace} -> S
diff --git a/src/tensors/diagonal.jl b/src/tensors/diagonal.jl
index 904e40cba..356f5dbf4 100644
--- a/src/tensors/diagonal.jl
+++ b/src/tensors/diagonal.jl
@@ -306,7 +306,7 @@ function LinearAlgebra.pinv(d::DiagonalTensorMap; kwargs...)
     if iszero(atol)
         rtol = get(kwargs, :rtol, zero(real(T)))
     else
-        rtol = sqrt(eps(real(float(oneunit(T))))) * length(d.data)
+        rtol = sqrt(eps(real(float(one(T))))) * length(d.data)
     end
     pdata = let tol = max(atol, rtol * maximum(abs, d.data))
         map(x -> abs(x) < tol ? zero(x) : pinv(x), d.data)
diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl
index 3cfe08a61..0c918e85e 100644
--- a/src/tensors/indexmanipulations.jl
+++ b/src/tensors/indexmanipulations.jl
@@ -267,7 +267,7 @@ function has_shared_twist(t, inds)
     if BraidingStyle(I) == NoBraiding()
         for i in inds
             cs = sectors(space(t, i))
-            all(isone, cs) || throw(SectorMismatch(lazy"Cannot twist sectors $cs"))
+            all(isunit, cs) || throw(SectorMismatch(lazy"Cannot twist sectors $cs"))
         end
         return true
     elseif BraidingStyle(I) == Bosonic()
diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl
index d09609edc..747a44320 100644
--- a/src/tensors/linalg.jl
+++ b/src/tensors/linalg.jl
@@ -296,7 +296,7 @@ function LinearAlgebra.rank(
         t::AbstractTensorMap;
         atol::Real = 0, rtol::Real = atol > 0 ? 0 : _default_rtol(t)
     )
-    r = dim(one(sectortype(t))) * 0
+    r = 0 * dim(first(allunits(sectortype(t))))
     dim(t) == 0 && return r
     S = LinearAlgebra.svdvals(t)
     tol = max(atol, rtol * maximum(first, values(S)))
diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl
index d22fa45d2..4a1d627c6 100644
--- a/src/tensors/tensoroperations.jl
+++ b/src/tensors/tensoroperations.jl
@@ -424,8 +424,9 @@ end
 
 # Scalar implementation
 #-----------------------
-function scalar(t::AbstractTensorMap)
-    # TODO: should scalar only work if N₁ == N₂ == 0?
-    return dim(codomain(t)) == dim(domain(t)) == 1 ?
-        first(blocks(t))[2][1, 1] : throw(DimensionMismatch())
+function scalar(t::AbstractTensorMap{T, S, 0, 0}) where {T, S}
+    Bs = collect(blocks(t))
+    inds = findall(!iszero ∘ last, Bs)
+    isempty(inds) && return zero(scalartype(t))
+    return only(last(Bs[only(inds)]))
 end
diff --git a/test/setup.jl b/test/setup.jl
index 0af9ce633..719f3ad36 100644
--- a/test/setup.jl
+++ b/test/setup.jl
@@ -1,8 +1,9 @@
 module TestSetup
 
 export smallset, randsector, hasfusiontensor, force_planar
+export random_fusion
 export sectorlist
-export Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VSU₂U₁, Vfib
+export Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VSU₂U₁, Vfib, VIB_diag, VIB_M
 
 using Random
 using TensorKit
@@ -25,14 +26,15 @@ end
 function randsector(::Type{I}) where {I <: Sector}
     s = collect(smallset(I))
     a = rand(s)
-    while a == one(a) # don't use trivial label
+    while isunit(a) # don't use trivial label
         a = rand(s)
     end
     return a
 end
 function hasfusiontensor(I::Type{<:Sector})
+    isa(UnitStyle(I), GenericUnit) && return false
     try
-        TensorKit.fusiontensor(one(I), one(I), one(I))
+        TensorKit.fusiontensor(unit(I), unit(I), unit(I))
         return true
     catch e
         if e isa MethodError
@@ -74,6 +76,18 @@ function force_planar(tsrc::TensorMap{<:Any, <:GradedSpace})
     return tdst
 end
 
+function random_fusion(I::Type{<:Sector}, ::Val{N}) where {N} # for fusion tree tests
+    N == 1 && return (randsector(I),)
+    tail = random_fusion(I, Val(N - 1))
+    s = randsector(I)
+    counter = 0
+    while isempty(⊗(s, first(tail))) && counter < 20
+        counter += 1
+        s = (counter < 20) ? randsector(I) : leftunit(first(tail))
+    end
+    return (s, tail...)
+end
+
 sectorlist = (
     Z2Irrep, Z3Irrep, Z4Irrep, Z3Irrep ⊠ Z4Irrep,
     U1Irrep, CU1Irrep, SU2Irrep,
@@ -81,6 +95,7 @@ sectorlist = (
     FermionParity ⊠ U1Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, # Hubbard-like
     FibonacciAnyon, IsingAnyon,
     Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon,
+    IsingBimodule, IsingBimodule ⊠ SU2Irrep, IsingBimodule ⊠ IsingBimodule,
 )
 
 # spaces
@@ -159,4 +174,25 @@ Vfib = (
     Vect[FibonacciAnyon](:I => 2, :τ => 2),
 )
 
+C0, C1 = IsingBimodule(1, 1, 0), IsingBimodule(1, 1, 1)
+D0, D1 = IsingBimodule(2, 2, 0), IsingBimodule(2, 2, 1)
+M, Mop = IsingBimodule(1, 2, 0), IsingBimodule(2, 1, 0)
+VIB_diag = (
+    Vect[IsingBimodule](C0 => 1, C1 => 2),
+    Vect[IsingBimodule](C0 => 2, C1 => 1),
+    Vect[IsingBimodule](C0 => 3, C1 => 1),
+    Vect[IsingBimodule](C0 => 2, C1 => 3),
+    Vect[IsingBimodule](C0 => 3, C1 => 2),
+)
+
+# not a random ordering! designed to make V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 work (tensors)
+# while V1 ⊗ V2 ← V4 isn't empty (factorizations)
+VIB_M = (
+    Vect[IsingBimodule](C0 => 1, C1 => 2),
+    Vect[IsingBimodule](M => 3),
+    Vect[IsingBimodule](C0 => 2, C1 => 3),
+    Vect[IsingBimodule](M => 4),
+    Vect[IsingBimodule](D0 => 3, D1 => 4),
+)
+
 end
diff --git a/test/symmetries/fusiontrees.jl b/test/symmetries/fusiontrees.jl
index c001f4e26..e5dbf8132 100644
--- a/test/symmetries/fusiontrees.jl
+++ b/test/symmetries/fusiontrees.jl
@@ -13,7 +13,7 @@ using .TestSetup
 @timedtestset "Fusion trees for $(TensorKit.type_repr(I))" verbose = true for I in sectorlist
     Istr = TensorKit.type_repr(I)
     N = 5
-    out = ntuple(n -> randsector(I), N)
+    out = random_fusion(I, Val(N))
     isdual = ntuple(n -> rand(Bool), N)
     in = rand(collect(⊗(out...)))
     numtrees = length(fusiontrees(out, in, isdual))
@@ -33,54 +33,66 @@ using .TestSetup
         @test eval(Meta.parse(sprint(show, f; context = (:module => @__MODULE__)))) == f
     end
     @testset "Fusion tree $Istr: constructor properties" begin
-        u = unit(I)
-        @constinferred FusionTree((), u, (), (), ())
-        @constinferred FusionTree((u,), u, (false,), (), ())
-        @constinferred FusionTree((u, u), u, (false, false), (), (1,))
-        @constinferred FusionTree((u, u, u), u, (false, false, false), (u,), (1, 1))
-        @constinferred FusionTree(
-            (u, u, u, u), u, (false, false, false, false), (u, u), (1, 1, 1)
-        )
-        @test_throws MethodError FusionTree((u, u, u), u, (false, false), (u,), (1, 1))
-        @test_throws MethodError FusionTree(
-            (u, u, u), u, (false, false, false), (u, u), (1, 1)
-        )
-        @test_throws MethodError FusionTree(
-            (u, u, u), u, (false, false, false), (u,), (1, 1, 1)
-        )
-        @test_throws MethodError FusionTree((u, u, u), u, (false, false, false), (), (1,))
-
-        f = FusionTree((u, u, u), u, (false, false, false), (u,), (1, 1))
-        @test sectortype(f) == I
-        @test length(f) == 3
-        @test FusionStyle(f) == FusionStyle(I)
-        @test BraidingStyle(f) == BraidingStyle(I)
-
-        if FusionStyle(I) isa UniqueFusion
-            @constinferred FusionTree((), u, ())
-            @constinferred FusionTree((u,), u, (false,))
-            @constinferred FusionTree((u, u), u, (false, false))
-            @constinferred FusionTree((u, u, u), u)
-            @constinferred FusionTree((u, u, u, u))
-            @test_throws MethodError FusionTree((u, u), u, (false, false, false))
-        else
-            errstr = "fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`"
-            @test_throws errstr FusionTree((), u, ())
-            @test_throws errstr FusionTree((u,), u, (false,))
-            @test_throws errstr FusionTree((u, u), u, (false, false))
-            @test_throws errstr FusionTree((u, u, u), u)
-            @test_throws errstr FusionTree((u, u, u, u))
+        for u in allunits(I)
+            @constinferred FusionTree((), u, (), (), ())
+            @constinferred FusionTree((u,), u, (false,), (), ())
+            @constinferred FusionTree((u, u), u, (false, false), (), (1,))
+            @constinferred FusionTree((u, u, u), u, (false, false, false), (u,), (1, 1))
+            @constinferred FusionTree(
+                (u, u, u, u), u, (false, false, false, false), (u, u), (1, 1, 1)
+            )
+            @test_throws MethodError FusionTree((u, u, u), u, (false, false), (u,), (1, 1))
+            @test_throws MethodError FusionTree(
+                (u, u, u), u, (false, false, false), (u, u), (1, 1)
+            )
+            @test_throws MethodError FusionTree(
+                (u, u, u), u, (false, false, false), (u,), (1, 1, 1)
+            )
+            @test_throws MethodError FusionTree((u, u, u), u, (false, false, false), (), (1,))
+
+            f = FusionTree((u, u, u), u, (false, false, false), (u,), (1, 1))
+            @test sectortype(f) == I
+            @test length(f) == 3
+            @test FusionStyle(f) == FusionStyle(I)
+            @test BraidingStyle(f) == BraidingStyle(I)
+
+            if FusionStyle(I) isa UniqueFusion
+                @constinferred FusionTree((), u, ())
+                @constinferred FusionTree((u,), u, (false,))
+                @constinferred FusionTree((u, u), u, (false, false))
+                @constinferred FusionTree((u, u, u), u)
+                if UnitStyle(I) isa SimpleUnit
+                    @constinferred FusionTree((u, u, u, u))
+                else
+                    @test_throws ArgumentError FusionTree((u, u, u, u))
+                end
+                @test_throws MethodError FusionTree((u, u), u, (false, false, false))
+            else
+                @test_throws ArgumentError FusionTree((), u, ())
+                @test_throws ArgumentError FusionTree((u,), u, (false,))
+                @test_throws ArgumentError FusionTree((u, u), u, (false, false))
+                @test_throws ArgumentError FusionTree((u, u, u), u)
+                if I <: ProductSector && UnitStyle(I) isa GenericUnit
+                    @test_throws DomainError FusionTree((u, u, u, u))
+                else
+                    @test_throws ArgumentError FusionTree((u, u, u, u))
+                end
+            end
         end
     end
     @testset "Fusion tree $Istr: insertat" begin
         N = 4
-        out2 = ntuple(n -> randsector(I), N)
+        out2 = random_fusion(I, Val(N))
         in2 = rand(collect(⊗(out2...)))
         isdual2 = ntuple(n -> rand(Bool), N)
         f2 = rand(collect(fusiontrees(out2, in2, isdual2)))
         for i in 1:N
-            out1 = ntuple(n -> randsector(I), N)
-            out1 = Base.setindex(out1, in2, i)
+            out1 = random_fusion(I, Val(N)) # guaranteed good fusion
+            out1 = Base.setindex(out1, in2, i) # can lead to poor fusion
+            while isempty(⊗(out1...)) # TODO: better way to do this?
+                out1 = random_fusion(I, Val(N))
+                out1 = Base.setindex(out1, in2, i)
+            end
             in1 = rand(collect(⊗(out1...)))
             isdual1 = ntuple(n -> rand(Bool), N)
             isdual1 = Base.setindex(isdual1, false, i)
@@ -93,25 +105,27 @@ using .TestSetup
             @test length(TK.insertat(f1b, 1, f1a)) == 1
             @test first(TK.insertat(f1b, 1, f1a)) == (f1 => 1)
 
-            levels = ntuple(identity, N)
-            function _reinsert_partial_tree(t, f)
-                (t′, c′) = first(TK.insertat(t, 1, f))
-                @test c′ == one(c′)
-                return t′
-            end
-            braid_i_to_1 = braid(f1, levels, (i, (1:(i - 1))..., ((i + 1):N)...))
-            trees2 = Dict(_reinsert_partial_tree(t, f2) => c for (t, c) in braid_i_to_1)
-            trees3 = empty(trees2)
-            p = (((N + 1):(N + i - 1))..., (1:N)..., ((N + i):(2N - 1))...)
-            levels = ((i:(N + i - 1))..., (1:(i - 1))..., ((i + N):(2N - 1))...)
-            for (t, coeff) in trees2
-                for (t′, coeff′) in braid(t, levels, p)
-                    trees3[t′] = get(trees3, t′, zero(coeff′)) + coeff * coeff′
+            if UnitStyle(I) isa SimpleUnit
+                levels = ntuple(identity, N)
+                function _reinsert_partial_tree(t, f)
+                    (t′, c′) = first(TK.insertat(t, 1, f))
+                    @test c′ == one(c′)
+                    return t′
+                end
+                braid_i_to_1 = braid(f1, levels, (i, (1:(i - 1))..., ((i + 1):N)...))
+                trees2 = Dict(_reinsert_partial_tree(t, f2) => c for (t, c) in braid_i_to_1)
+                trees3 = empty(trees2)
+                p = (((N + 1):(N + i - 1))..., (1:N)..., ((N + i):(2N - 1))...)
+                levels = ((i:(N + i - 1))..., (1:(i - 1))..., ((i + N):(2N - 1))...)
+                for (t, coeff) in trees2
+                    for (t′, coeff′) in braid(t, levels, p)
+                        trees3[t′] = get(trees3, t′, zero(coeff′)) + coeff * coeff′
+                    end
+                end
+                for (t, coeff) in trees3
+                    coeff′ = get(trees, t, zero(coeff))
+                    @test isapprox(coeff′, coeff; atol = 1.0e-12, rtol = 1.0e-12)
                 end
-            end
-            for (t, coeff) in trees3
-                coeff′ = get(trees, t, zero(coeff))
-                @test isapprox(coeff′, coeff; atol = 1.0e-12, rtol = 1.0e-12)
             end
 
             if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
@@ -223,13 +237,13 @@ using .TestSetup
             end
         end
     end
-    @testset "Fusion tree $Istr: elementary artin braid" begin
+    (BraidingStyle(I) isa HasBraiding) && @testset "Fusion tree $Istr: elementary artin braid" begin
         N = length(out)
         isdual = ntuple(n -> rand(Bool), N)
         for in in ⊗(out...)
             for i in 1:(N - 1)
                 for f in fusiontrees(out, in, isdual)
-                    d1 = @constinferred TK.artin_braid(f, i)
+                    d1 = @constinferred TK.artin_braid(f, i) # braid between random objects
                     @test norm(values(d1)) ≈ 1
                     d2 = empty(d1)
                     for (f1, coeff1) in d1
@@ -279,7 +293,7 @@ using .TestSetup
             end
         end
     end
-    @testset "Fusion tree $Istr: braiding and permuting" begin
+    (BraidingStyle(I) isa HasBraiding) && @testset "Fusion tree $Istr: braiding and permuting" begin
         f = rand(collect(it))
         p = tuple(randperm(N)...)
         ip = invperm(p)
@@ -314,11 +328,20 @@ using .TestSetup
 
     @testset "Fusion tree $Istr: merging" begin
         N = 3
-        out1 = ntuple(n -> randsector(I), N)
+        out1 = random_fusion(I, Val(N))
+        out2 = random_fusion(I, Val(N))
         in1 = rand(collect(⊗(out1...)))
-        f1 = rand(collect(fusiontrees(out1, in1)))
-        out2 = ntuple(n -> randsector(I), N)
         in2 = rand(collect(⊗(out2...)))
+        tp = ⊗(in1, in2) # messy solution but it works
+        while isempty(tp)
+            out1 = random_fusion(I, Val(N))
+            out2 = random_fusion(I, Val(N))
+            in1 = rand(collect(⊗(out1...)))
+            in2 = rand(collect(⊗(out2...)))
+            tp = ⊗(in1, in2)
+        end
+
+        f1 = rand(collect(fusiontrees(out1, in1)))
         f2 = rand(collect(fusiontrees(out2, in2)))
 
         @constinferred TK.merge(f1, f2, first(in1 ⊗ in2), 1)
@@ -332,50 +355,52 @@ using .TestSetup
                 for (f, coeff) in TK.merge(f1, f2, c, μ)
         )
 
-        for c in in1 ⊗ in2
-            R = Rsymbol(in1, in2, c)
-            for μ in 1:Nsymbol(in1, in2, c)
-                trees1 = TK.merge(f1, f2, c, μ)
-
-                # test merge and braid interplay
-                trees2 = Dict{keytype(trees1), complex(valtype(trees1))}()
-                trees3 = Dict{keytype(trees1), complex(valtype(trees1))}()
-                for ν in 1:Nsymbol(in2, in1, c)
-                    for (t, coeff) in TK.merge(f2, f1, c, ν)
-                        trees2[t] = get(trees2, t, zero(valtype(trees2))) + coeff * R[μ, ν]
+        if BraidingStyle(I) isa HasBraiding
+            for c in in1 ⊗ in2
+                R = Rsymbol(in1, in2, c)
+                for μ in 1:Nsymbol(in1, in2, c)
+                    trees1 = TK.merge(f1, f2, c, μ)
+
+                    # test merge and braid interplay
+                    trees2 = Dict{keytype(trees1), complex(valtype(trees1))}()
+                    trees3 = Dict{keytype(trees1), complex(valtype(trees1))}()
+                    for ν in 1:Nsymbol(in2, in1, c)
+                        for (t, coeff) in TK.merge(f2, f1, c, ν)
+                            trees2[t] = get(trees2, t, zero(valtype(trees2))) + coeff * R[μ, ν]
+                        end
                     end
-                end
-                perm = ((N .+ (1:N))..., (1:N)...)
-                levels = ntuple(identity, 2 * N)
-                for (t, coeff) in trees1
-                    for (t′, coeff′) in braid(t, levels, perm)
-                        trees3[t′] = get(trees3, t′, zero(valtype(trees3))) + coeff * coeff′
+                    perm = ((N .+ (1:N))..., (1:N)...)
+                    levels = ntuple(identity, 2 * N)
+                    for (t, coeff) in trees1
+                        for (t′, coeff′) in braid(t, levels, perm)
+                            trees3[t′] = get(trees3, t′, zero(valtype(trees3))) + coeff * coeff′
+                        end
+                    end
+                    for (t, coeff) in trees3
+                        coeff′ = get(trees2, t, zero(coeff))
+                        @test isapprox(coeff, coeff′; atol = 1.0e-12, rtol = 1.0e-12)
                     end
-                end
-                for (t, coeff) in trees3
-                    coeff′ = get(trees2, t, zero(coeff))
-                    @test isapprox(coeff, coeff′; atol = 1.0e-12, rtol = 1.0e-12)
-                end
 
-                # test via conversion
-                if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
-                    Af1 = convert(Array, f1)
-                    Af2 = convert(Array, f2)
-                    Af0 = convert(
-                        Array,
-                        FusionTree((f1.coupled, f2.coupled), c, (false, false), (), (μ,))
-                    )
-                    _Af = TensorOperations.tensorcontract(
-                        1:(N + 2), Af1, [1:N; -1], Af0, [-1; N + 1; N + 2]
-                    )
-                    Af = TensorOperations.tensorcontract(
-                        1:(2N + 1), Af2, [N .+ (1:N); -1], _Af, [1:N; -1; 2N + 1]
-                    )
-                    Af′ = zero(Af)
-                    for (f, coeff) in trees1
-                        Af′ .+= coeff .* convert(Array, f)
+                    # test via conversion
+                    if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
+                        Af1 = convert(Array, f1)
+                        Af2 = convert(Array, f2)
+                        Af0 = convert(
+                            Array,
+                            FusionTree((f1.coupled, f2.coupled), c, (false, false), (), (μ,))
+                        )
+                        _Af = TensorOperations.tensorcontract(
+                            1:(N + 2), Af1, [1:N; -1], Af0, [-1; N + 1; N + 2]
+                        )
+                        Af = TensorOperations.tensorcontract(
+                            1:(2N + 1), Af2, [N .+ (1:N); -1], _Af, [1:N; -1; 2N + 1]
+                        )
+                        Af′ = zero(Af)
+                        for (f, coeff) in trees1
+                            Af′ .+= coeff .* convert(Array, f)
+                        end
+                        @test Af ≈ Af′
                     end
-                    @test Af ≈ Af′
                 end
             end
         end
@@ -386,15 +411,30 @@ using .TestSetup
     else
         N = 4
     end
-    out = ntuple(n -> randsector(I), N)
-    numtrees = count(n -> true, fusiontrees((out..., map(dual, out)...)))
-    while !(0 < numtrees < 100)
-        out = ntuple(n -> randsector(I), N)
+    if UnitStyle(I) isa SimpleUnit
+        out = random_fusion(I, Val(N))
         numtrees = count(n -> true, fusiontrees((out..., map(dual, out)...)))
+        while !(0 < numtrees < 100)
+            out = random_fusion(I, Val(N))
+            numtrees = count(n -> true, fusiontrees((out..., map(dual, out)...)))
+        end
+        incoming = rand(collect(⊗(out...)))
+        f1 = rand(collect(fusiontrees(out, incoming, ntuple(n -> rand(Bool), N))))
+        f2 = rand(collect(fusiontrees(out[randperm(N)], incoming, ntuple(n -> rand(Bool), N))))
+    else
+        out = random_fusion(I, Val(N))
+        out2 = random_fusion(I, Val(N))
+        tp = ⊗(out...)
+        tp2 = ⊗(out2...)
+        while isempty(intersect(tp, tp2)) # guarantee fusion to same coloring
+            out2 = random_fusion(I, Val(N))
+            tp2 = ⊗(out2...)
+        end
+        @test_throws ArgumentError fusiontrees((out..., map(dual, out)...))
+        incoming = rand(collect(intersect(tp, tp2)))
+        f1 = rand(collect(fusiontrees(out, incoming, ntuple(n -> rand(Bool), N))))
+        f2 = rand(collect(fusiontrees(out2, incoming, ntuple(n -> rand(Bool), N)))) # no permuting
     end
-    incoming = rand(collect(⊗(out...)))
-    f1 = rand(collect(fusiontrees(out, incoming, ntuple(n -> rand(Bool), N))))
-    f2 = rand(collect(fusiontrees(out[randperm(N)], incoming, ntuple(n -> rand(Bool), N))))
 
     @testset "Double fusion tree $Istr: repartitioning" begin
         for n in 0:(2 * N)
@@ -473,32 +513,11 @@ using .TestSetup
                 end
 
                 if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
-                    Af1 = convert(Array, f1)
-                    Af2 = convert(Array, f2)
-                    sz1 = size(Af1)
-                    sz2 = size(Af2)
-                    d1 = prod(sz1[1:(end - 1)])
-                    d2 = prod(sz2[1:(end - 1)])
-                    dc = sz1[end]
-                    A = reshape(
-                        reshape(Af1, (d1, dc)) * reshape(Af2, (d2, dc))',
-                        (sz1[1:(end - 1)]..., sz2[1:(end - 1)]...)
-                    )
+                    A = convert(Array, (f1, f2))
                     Ap = permutedims(A, (p1..., p2...))
                     A2 = zero(Ap)
                     for ((f1′, f2′), coeff) in d
-                        Af1′ = convert(Array, f1′)
-                        Af2′ = convert(Array, f2′)
-                        sz1′ = size(Af1′)
-                        sz2′ = size(Af2′)
-                        d1′ = prod(sz1′[1:(end - 1)])
-                        d2′ = prod(sz2′[1:(end - 1)])
-                        dc′ = sz1′[end]
-                        A2 += coeff * reshape(
-                            reshape(Af1′, (d1′, dc′)) *
-                                reshape(Af2′, (d2′, dc′))',
-                            (sz1′[1:(end - 1)]..., sz2′[1:(end - 1)]...)
-                        )
+                        A2 .+= coeff .* convert(Array, (f1′, f2′))
                     end
                     @test Ap ≈ A2
                 end
@@ -577,7 +596,7 @@ using .TestSetup
     @testset "Double fusion tree $Istr: planar trace" begin
         d1 = transpose(f1, f1, (N + 1, 1:N..., ((2N):-1:(N + 3))...), (N + 2,))
         f1front, = TK.split(f1, N - 1)
-        T = TensorKitSectors._Fscalartype(I)
+        T = sectorscalartype(I)
         d2 = Dict{typeof((f1front, f1front)), T}()
         for ((f1′, f2′), coeff′) in d1
             for ((f1′′, f2′′), coeff′′) in
diff --git a/test/symmetries/spaces.jl b/test/symmetries/spaces.jl
index ede4ae52f..75e9fd0b2 100644
--- a/test/symmetries/spaces.jl
+++ b/test/symmetries/spaces.jl
@@ -87,7 +87,9 @@ end
     @test @constinferred(axes(V)) == Base.OneTo(d)
     @test ℝ^d == ℝ[](d) == CartesianSpace(d) == typeof(V)(d)
     W = @constinferred ℝ^1
+    @test @constinferred(isunitspace(W))
     @test @constinferred(unitspace(V)) == W == unitspace(typeof(V))
+    @test @constinferred(leftunitspace(V)) == W == @constinferred(rightunitspace(V))
     @test @constinferred(zerospace(V)) == ℝ^0 == zerospace(typeof(V))
     @test @constinferred(⊕(V, zerospace(V))) == V
     @test @constinferred(⊕(V, V)) == ℝ^(2d)
@@ -133,7 +135,9 @@ end
     @test @constinferred(axes(V)) == Base.OneTo(d)
     @test ℂ^d == Vect[Trivial](d) == Vect[](Trivial() => d) == ℂ[](d) == typeof(V)(d)
     W = @constinferred ℂ^1
+    @test @constinferred(isunitspace(W))
     @test @constinferred(unitspace(V)) == W == unitspace(typeof(V))
+    @test @constinferred(leftunitspace(V)) == W == @constinferred(rightunitspace(V))
     @test @constinferred(zerospace(V)) == ℂ^0 == zerospace(typeof(V))
     @test @constinferred(⊕(V, zerospace(V))) == V
     @test @constinferred(⊕(V, V)) == ℂ^(2d)
@@ -187,7 +191,7 @@ end
 
 @timedtestset "ElementarySpace: $(type_repr(Vect[I]))" for I in sectorlist
     if Base.IteratorSize(values(I)) === Base.IsInfinite()
-        set = unique(vcat(unit(I), [randsector(I) for k in 1:10]))
+        set = unique(vcat(allunits(I)..., [randsector(I) for k in 1:10]))
         gen = (c => 2 for c in set)
     else
         gen = (values(I)[k] => (k + 1) for k in 1:length(values(I)))
@@ -220,13 +224,23 @@ end
     @test eval_show(typeof(V)) == typeof(V)
     # space with no sectors
     @test dim(@constinferred(zerospace(V))) == 0
-    # space with a single sector
-    W = @constinferred GradedSpace(unit(I) => 1)
-    @test W == GradedSpace(unit(I) => 1, randsector(I) => 0)
+    # space with unit(s), always test as if multifusion
+    W = @constinferred GradedSpace(unit => 1 for unit in allunits(I))
+    dict = Dict(unit => 1 for unit in allunits(I))
+    @test W == GradedSpace(dict)
+    @test W == GradedSpace(push!(dict, randsector(I) => 0))
+    @test @constinferred(zerospace(V)) == GradedSpace(unit => 0 for unit in allunits(I))
+    randunit = rand(collect(allunits(I)))
+    @test_throws ArgumentError("Sector $(randunit) appears multiple times") GradedSpace(randunit => 1, randunit => 3)
+
+    @test isunitspace(W)
     @test @constinferred(unitspace(V)) == W == unitspace(typeof(V))
-    @test @constinferred(zerospace(V)) == GradedSpace(unit(I) => 0)
-    # randsector never returns trivial sector, so this cannot error
-    @test_throws ArgumentError GradedSpace(unit(I) => 1, randsector(I) => 0, unit(I) => 3)
+    if UnitStyle(I) isa SimpleUnit
+        @test @constinferred(leftunitspace(V)) == W == @constinferred(rightunitspace(V))
+    else
+        @test_throws ArgumentError leftunitspace(V)
+        @test_throws ArgumentError rightunitspace(V)
+    end
     @test eval_show(W) == W
     @test isa(V, VectorSpace)
     @test isa(V, ElementarySpace)
@@ -265,7 +279,9 @@ end
     @test V == @constinferred infimum(V, ⊕(V, V))
     @test V ≺ ⊕(V, V)
     @test !(V ≻ ⊕(V, V))
-    @test infimum(V, GradedSpace(unit(I) => 3)) == GradedSpace(unit(I) => 2)
+
+    u = first(allunits(I))
+    @test infimum(V, GradedSpace(u => 3)) == GradedSpace(u => 2)
     @test_throws SpaceMismatch (⊕(V, V'))
 end
 
diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl
index 2c1638a30..b25a9dbbd 100644
--- a/test/tensors/factorizations.jl
+++ b/test/tensors/factorizations.jl
@@ -9,17 +9,17 @@ spacelist = try
     if ENV["CI"] == "true"
         println("Detected running on CI")
         if Sys.iswindows()
-            (Vtr, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂)
+            (Vtr, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VIB_diag)
         elseif Sys.isapple()
-            (Vtr, Vℤ₃, VfU₁, VfSU₂)
+            (Vtr, Vℤ₃, VfU₁, VfSU₂, VIB_M)
         else
-            (Vtr, VU₁, VCU₁, VSU₂, VfSU₂)
+            (Vtr, VU₁, VCU₁, VSU₂, VfSU₂, VIB_diag, VIB_M)
         end
     else
-        (Vtr, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂)
+        (Vtr, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VIB_diag, VIB_M)
     end
 catch
-    (Vtr, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂)
+    (Vtr, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VIB_diag, VIB_M)
 end
 
 eltypes = (Float32, ComplexF64)
@@ -33,11 +33,13 @@ for V in spacelist
     @timedtestset "Factorizations with symmetry: $Istr" verbose = true begin
         V1, V2, V3, V4, V5 = V
         W = V1 ⊗ V2
+        @assert !isempty(blocksectors(W))
+        @assert !isempty(intersect(blocksectors(V4), blocksectors(W)))
 
         @testset "QR decomposition" begin
             for T in eltypes,
                     t in (
-                        rand(T, W, W), rand(T, W, W)', rand(T, W, V1), rand(T, V1, W)',
+                        rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)',
                         DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 
@@ -64,7 +66,7 @@ for V in spacelist
 
             # empty tensor
             for T in eltypes
-                t = rand(T, V1 ⊗ V2, zero(V1))
+                t = rand(T, V1 ⊗ V2, zerospace(V1))
 
                 Q, R = @constinferred qr_full(t)
                 @test Q * R ≈ t
@@ -90,7 +92,7 @@ for V in spacelist
         @testset "LQ decomposition" begin
             for T in eltypes,
                     t in (
-                        rand(T, W, W), rand(T, W, W)', rand(T, W, V1), rand(T, V1, W)',
+                        rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)',
                         DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 
@@ -113,7 +115,7 @@ for V in spacelist
 
             for T in eltypes
                 # empty tensor
-                t = rand(T, zero(V1), V1 ⊗ V2)
+                t = rand(T, zerospace(V1), V1 ⊗ V2)
 
                 L, Q = @constinferred lq_full(t)
                 @test L * Q ≈ t
@@ -139,7 +141,7 @@ for V in spacelist
         @testset "Polar decomposition" begin
             for T in eltypes,
                     t in (
-                        rand(T, W, W), rand(T, W, W)', rand(T, W, V1), rand(T, V1, W)',
+                        rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)',
                         DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 
@@ -155,7 +157,7 @@ for V in spacelist
             end
 
             for T in eltypes,
-                    t in (rand(T, W, W), rand(T, W, W)', rand(T, V1, W), rand(T, W, V1)')
+                    t in (rand(T, W, W), rand(T, W, W)', rand(T, V4, W), rand(T, W, V4)')
 
                 @assert codomain(t) ≾ domain(t)
                 p, wᴴ = @constinferred right_polar(t)
@@ -173,8 +175,8 @@ for V in spacelist
             for T in eltypes,
                     t in (
                         rand(T, W, W), rand(T, W, W)',
-                        rand(T, W, V1), rand(T, V1, W),
-                        rand(T, W, V1)', rand(T, V1, W)',
+                        rand(T, W, V4), rand(T, V4, W),
+                        rand(T, W, V4)', rand(T, V4, W)',
                         DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 
@@ -220,7 +222,7 @@ for V in spacelist
             end
 
             # empty tensor
-            for T in eltypes, t in (rand(T, W, zero(V1)), rand(T, zero(V1), W))
+            for T in eltypes, t in (rand(T, W, zerospace(V1)), rand(T, zerospace(V1), W))
                 U, S, Vᴴ = @constinferred svd_full(t)
                 @test U * S * Vᴴ ≈ t
                 @test isunitary(U)
@@ -236,8 +238,8 @@ for V in spacelist
             for T in eltypes,
                     t in (
                         randn(T, W, W), randn(T, W, W)',
-                        randn(T, W, V1), randn(T, V1, W),
-                        randn(T, W, V1)', randn(T, V1, W)',
+                        randn(T, W, V4), randn(T, V4, W),
+                        randn(T, W, V4)', randn(T, V4, W)',
                         DiagonalTensorMap(randn(T, reduceddim(V1)), V1),
                     )
 
@@ -249,13 +251,15 @@ for V in spacelist
                 @test isisometric(U)
                 @test isisometric(Vᴴ; side = :right)
 
-                trunc = truncrank(dim(domain(S)) ÷ 2)
+                # dimension of S is a float for IsingBimodule
+                nvals = round(Int, dim(domain(S)) / 2)
+                trunc = truncrank(nvals)
                 U1, S1, Vᴴ1, ϵ1 = @constinferred svd_trunc(t; trunc)
                 @test t * Vᴴ1' ≈ U1 * S1
                 @test isisometric(U1)
                 @test isisometric(Vᴴ1; side = :right)
                 @test norm(t - U1 * S1 * Vᴴ1) ≈ ϵ1 atol = eps(real(T))^(4 / 5)
-                @test dim(domain(S1)) <= trunc.howmany
+                @test dim(domain(S1)) <= nvals
 
                 λ = minimum(minimum, values(LinearAlgebra.diag(S1)))
                 trunc = trunctol(; atol = λ - 10eps(λ))
@@ -286,14 +290,14 @@ for V in spacelist
                 @test norm(t - U4 * S4 * Vᴴ4) ≈ ϵ4 atol = eps(real(T))^(4 / 5)
                 @test ϵ4 ≤ ϵ2
 
-                trunc = truncrank(dim(domain(S)) ÷ 2) & trunctol(; atol = λ - 10eps(λ))
+                trunc = truncrank(nvals) & trunctol(; atol = λ - 10eps(λ))
                 U5, S5, Vᴴ5, ϵ5 = @constinferred svd_trunc(t; trunc)
                 @test t * Vᴴ5' ≈ U5 * S5
                 @test isisometric(U5)
                 @test isisometric(Vᴴ5; side = :right)
                 @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5)
                 @test minimum(minimum, values(LinearAlgebra.diag(S5))) >= λ
-                @test dim(domain(S5)) ≤ dim(domain(S)) ÷ 2
+                @test dim(domain(S5)) ≤ nvals
             end
         end
 
@@ -317,9 +321,10 @@ for V in spacelist
                 @test @constinferred isposdef(vdv)
                 t isa DiagonalTensorMap || @test !isposdef(t) # unlikely for non-hermitian map
 
-                d, v = @constinferred eig_trunc(t; trunc = truncrank(dim(domain(t)) ÷ 2))
+                nvals = round(Int, dim(domain(t)) / 2)
+                d, v = @constinferred eig_trunc(t; trunc = truncrank(nvals))
                 @test t * v ≈ v * d
-                @test dim(domain(d)) ≤ dim(domain(t)) ÷ 2
+                @test dim(domain(d)) ≤ nvals
 
                 t2 = (t + t')
                 D, V = eigen(t2)
@@ -348,9 +353,9 @@ for V in spacelist
                 @test isposdef(t - λ * one(t) + 0.1 * one(t))
                 @test !isposdef(t - λ * one(t) - 0.1 * one(t))
 
-                d, v = @constinferred eigh_trunc(t; trunc = truncrank(dim(domain(t)) ÷ 2))
+                d, v = @constinferred eigh_trunc(t; trunc = truncrank(nvals))
                 @test t * v ≈ v * d
-                @test dim(domain(d)) ≤ dim(domain(t)) ÷ 2
+                @test dim(domain(d)) ≤ nvals
             end
         end
 
@@ -358,26 +363,28 @@ for V in spacelist
             for T in eltypes,
                     t in (
                         rand(T, W, W), rand(T, W, W)',
-                        rand(T, W, V1), rand(T, V1, W),
-                        rand(T, W, V1)', rand(T, V1, W)',
+                        rand(T, W, V4), rand(T, V4, W),
+                        rand(T, W, V4)', rand(T, V4, W)',
                         DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 
                 d1, d2 = dim(codomain(t)), dim(domain(t))
-                @test rank(t) == min(d1, d2)
+                r = rank(t)
+                @test r == min(d1, d2)
+                @test typeof(r) == typeof(d1)
                 M = left_null(t)
-                @test @constinferred(rank(M)) + rank(t) == d1
+                @test @constinferred(rank(M)) + r ≈ d1
                 Mᴴ = right_null(t)
-                @test rank(Mᴴ) + rank(t) == d2
+                @test rank(Mᴴ) + r ≈ d2
             end
             for T in eltypes
                 u = unitary(T, V1 ⊗ V2, V1 ⊗ V2)
                 @test @constinferred(cond(u)) ≈ one(real(T))
                 @test @constinferred(rank(u)) == dim(V1 ⊗ V2)
 
-                t = rand(T, zero(V1), W)
+                t = rand(T, zerospace(V1), W)
                 @test rank(t) == 0
-                t2 = rand(T, zero(V1) * zero(V2), zero(V1) * zero(V2))
+                t2 = rand(T, zerospace(V1) * zerospace(V2), zerospace(V1) * zerospace(V2))
                 @test rank(t2) == 0
                 @test cond(t2) == 0.0
             end
@@ -425,7 +432,7 @@ for V in spacelist
             for T in eltypes,
                     t in (
                         randn(T, W, W), randn(T, W, W)',
-                        randn(T, W, V1), randn(T, V1, W)',
+                        randn(T, W, V4), randn(T, V4, W)',
                     )
                 t2 = project_isometric(t)
                 @test isisometric(t2)
diff --git a/test/tensors/tensors.jl b/test/tensors/tensors.jl
index 34934eb34..43481e7d1 100644
--- a/test/tensors/tensors.jl
+++ b/test/tensors/tensors.jl
@@ -11,22 +11,23 @@ spacelist = try
     if get(ENV, "CI", "false") == "true"
         println("Detected running on CI")
         if Sys.iswindows()
-            (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂)
+            (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VIB_diag)
         elseif Sys.isapple()
-            (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VfU₁, VfSU₂, VSU₂U₁) #, VSU₃)
+            (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VfU₁, VfSU₂, VSU₂U₁, VIB_M) #, VSU₃)
         else
-            (Vtr, Vℤ₂, Vfℤ₂, VU₁, VCU₁, VSU₂, VfSU₂, VSU₂U₁) #, VSU₃)
+            (Vtr, Vℤ₂, Vfℤ₂, VU₁, VCU₁, VSU₂, VfSU₂, VSU₂U₁, VIB_diag, VIB_M) #, VSU₃)
         end
     else
-        (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VSU₂U₁) #, VSU₃)
+        (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VSU₂U₁, VIB_diag, VIB_M) #, VSU₃)
     end
 catch
-    (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VSU₂U₁) #, VSU₃)
+    (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VSU₂U₁, VIB_diag, VIB_M) #, VSU₃)
 end
 
 for V in spacelist
     I = sectortype(first(V))
     Istr = type_repr(I)
+    symmetricbraiding = BraidingStyle(I) isa SymmetricBraiding
     println("---------------------------------------")
     println("Tensors with symmetry: $Istr")
     println("---------------------------------------")
@@ -45,18 +46,20 @@ for V in spacelist
                 @test typeof(t) == TensorMap{T, spacetype(t), 5, 0, Vector{T}}
                 # blocks
                 bs = @constinferred blocks(t)
-                (c, b1), state = @constinferred Nothing iterate(bs)
-                @test c == first(blocksectors(W))
-                next = @constinferred Nothing iterate(bs, state)
-                b2 = @constinferred block(t, first(blocksectors(t)))
-                @test b1 == b2
-                @test eltype(bs) === Pair{typeof(c), typeof(b1)}
-                @test typeof(b1) === TensorKit.blocktype(t)
-                @test typeof(c) === sectortype(t)
+                if !isempty(blocksectors(t)) # multifusion space ending on module gives empty data
+                    (c, b1), state = @constinferred Nothing iterate(bs)
+                    @test c == first(blocksectors(W))
+                    next = @constinferred Nothing iterate(bs, state)
+                    b2 = @constinferred block(t, first(blocksectors(t)))
+                    @test b1 == b2
+                    @test eltype(bs) === Pair{typeof(c), typeof(b1)}
+                    @test typeof(b1) === TensorKit.blocktype(t)
+                    @test typeof(c) === sectortype(t)
+                end
             end
         end
         @timedtestset "Tensor Dict conversion" begin
-            W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
+            W = V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5
             for T in (Int, Float32, ComplexF64)
                 t = @constinferred rand(T, W)
                 d = convert(Dict, t)
@@ -96,7 +99,7 @@ for V in spacelist
             end
         end
         @timedtestset "Basic linear algebra" begin
-            W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
+            W = V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5
             for T in (Float32, ComplexF64)
                 t = @constinferred rand(T, W)
                 @test scalartype(t) == T
@@ -134,41 +137,47 @@ for V in spacelist
                 @test dot(t2, t) ≈ conj(dot(t2', t'))
                 @test dot(t2, t) ≈ dot(t', t2')
 
-                i1 = @constinferred(isomorphism(T, V1 ⊗ V2, V2 ⊗ V1))
-                i2 = @constinferred(isomorphism(Vector{T}, V2 ⊗ V1, V1 ⊗ V2))
-                @test i1 * i2 == @constinferred(id(T, V1 ⊗ V2))
-                @test i2 * i1 == @constinferred(id(Vector{T}, V2 ⊗ V1))
+                if UnitStyle(I) isa SimpleUnit || !isempty(blocksectors(V2 ⊗ V1))
+                    i1 = @constinferred(isomorphism(T, V1 ⊗ V2, V2 ⊗ V1)) # can't reverse fusion here when modules are involved
+                    i2 = @constinferred(isomorphism(Vector{T}, V2 ⊗ V1, V1 ⊗ V2))
+                    @test i1 * i2 == @constinferred(id(T, V1 ⊗ V2))
+                    @test i2 * i1 == @constinferred(id(Vector{T}, V2 ⊗ V1))
+                end
 
-                w = @constinferred(isometry(T, V1 ⊗ (oneunit(V1) ⊕ oneunit(V1)), V1))
+                w = @constinferred isometry(T, V1 ⊗ (rightunitspace(V1) ⊕ rightunitspace(V1)), V1)
                 @test dim(w) == 2 * dim(V1 ← V1)
                 @test w' * w == id(Vector{T}, V1)
                 @test w * w' == (w * w')^2
             end
         end
         @timedtestset "Trivial space insertion and removal" begin
-            W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
+            W = V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5
             for T in (Float32, ComplexF64)
                 t = @constinferred rand(T, W)
                 t2 = @constinferred insertleftunit(t)
                 @test t2 == @constinferred insertrightunit(t)
-                @test numind(t2) == numind(t) + 1
                 @test space(t2) == insertleftunit(space(t))
-                @test scalartype(t2) === T
-                @test t.data === t2.data
                 @test @constinferred(removeunit(t2, $(numind(t2)))) == t
                 t3 = @constinferred insertleftunit(t; copy = true)
                 @test t3 == @constinferred insertrightunit(t; copy = true)
+                @test @constinferred(removeunit(t3, $(numind(t3)))) == t
+
+                @test numind(t2) == numind(t) + 1
+                @test scalartype(t2) === T
+                @test t.data === t2.data
+
                 @test t.data !== t3.data
                 for (c, b) in blocks(t)
                     @test b == block(t3, c)
                 end
-                @test @constinferred(removeunit(t3, $(numind(t3)))) == t
+
                 t4 = @constinferred insertrightunit(t, 3; dual = true)
-                @test numin(t4) == numin(t) && numout(t4) == numout(t) + 1
+                @test numin(t4) == numin(t) + 1 && numout(t4) == numout(t)
                 for (c, b) in blocks(t)
                     @test b == block(t4, c)
                 end
                 @test @constinferred(removeunit(t4, 4)) == t
+
                 t5 = @constinferred insertleftunit(t, 4; dual = true)
                 @test numin(t5) == numin(t) + 1 && numout(t5) == numout(t)
                 for (c, b) in blocks(t)
@@ -223,7 +232,7 @@ for V in spacelist
             @test Base.promote_typeof(t, tc) == typeof(tc)
             @test Base.promote_typeof(tc, t) == typeof(tc + t)
         end
-        @timedtestset "Permutations: test via inner product invariance" begin
+        symmetricbraiding && @timedtestset "Permutations: test via inner product invariance" begin
             W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5
             t = rand(ComplexF64, W)
             t′ = randn!(similar(t))
@@ -237,7 +246,7 @@ for V in spacelist
                     @test dot(t2′, t2) ≈ dot(t′, t) ≈ dot(transpose(t2′), transpose(t2))
                 end
 
-                t3 = VERSION < v"1.7" ? repartition(t, k) : @constinferred repartition(t, $k)
+                t3 = @constinferred repartition(t, $k)
                 @test norm(t3) ≈ norm(t)
                 t3′ = @constinferred repartition!(similar(t3), t′)
                 @test norm(t3′) ≈ norm(t′)
@@ -269,28 +278,46 @@ for V in spacelist
             end
         end
         @timedtestset "Full trace: test self-consistency" begin
-            t = rand(ComplexF64, V1 ⊗ V2' ⊗ V2 ⊗ V1')
-            t2 = permute(t, ((1, 2), (4, 3)))
-            s = @constinferred tr(t2)
-            @test conj(s) ≈ tr(t2')
-            if !isdual(V1)
-                t2 = twist!(t2, 1)
-            end
-            if isdual(V2)
-                t2 = twist!(t2, 2)
+            if symmetricbraiding
+                t = rand(ComplexF64, V1 ⊗ V2' ⊗ V2 ⊗ V1')
+                t2 = permute(t, ((1, 2), (4, 3)))
+                s = @constinferred tr(t2)
+                @test conj(s) ≈ tr(t2')
+                if !isdual(V1)
+                    t2 = twist!(t2, 1)
+                end
+                if isdual(V2)
+                    t2 = twist!(t2, 2)
+                end
+                ss = tr(t2)
+                @tensor s2 = t[a, b, b, a]
+                @tensor t3[a, b] := t[a, c, c, b]
+                @tensor s3 = t3[a, a]
+                @test ss ≈ s2
+                @test ss ≈ s3
             end
-            ss = tr(t2)
-            @tensor s2 = t[a, b, b, a]
-            @tensor t3[a, b] := t[a, c, c, b]
-            @tensor s3 = t3[a, a]
+            t = rand(ComplexF64, V1 ⊗ V2 ← V1 ⊗ V2) # avoid permutes
+            ss = @constinferred tr(t)
+            @test conj(ss) ≈ tr(t')
+            @planar s2 = t[a b; a b]
+            @planar t3[a; b] := t[a c; b c]
+            @planar s3 = t3[a; a]
+
             @test ss ≈ s2
             @test ss ≈ s3
         end
         @timedtestset "Partial trace: test self-consistency" begin
-            t = rand(ComplexF64, V1 ⊗ V2' ⊗ V3 ⊗ V2 ⊗ V1' ⊗ V3')
-            @tensor t2[a, b] := t[c, d, b, d, c, a]
-            @tensor t4[a, b, c, d] := t[d, e, b, e, c, a]
-            @tensor t5[a, b] := t4[a, b, c, c]
+            if symmetricbraiding
+                t = rand(ComplexF64, V1 ⊗ V2 ⊗ V3 ← V1 ⊗ V2 ⊗ V3)
+                @tensor t2[a; b] := t[c d b; c d a]
+                @tensor t4[a b; c d] := t[e d c; e b a]
+                @tensor t5[a; b] := t4[a c; b c]
+                @test t2 ≈ t5
+            end
+            t = rand(ComplexF64, V3 ⊗ V4 ⊗ V5 ← V3 ⊗ V4 ⊗ V5) # compatible with module fusion
+            @planar t2[a; b] := t[c a d; c b d]
+            @planar t4[a b; c d] := t[e a b; e c d]
+            @planar t5[a; b] := t4[a c; b c]
             @test t2 ≈ t5
         end
         if BraidingStyle(I) isa Bosonic && hasfusiontensor(I)
@@ -301,7 +328,8 @@ for V in spacelist
                 @test t3 ≈ convert(Array, t2)
             end
         end
-        @timedtestset "Trace and contraction" begin
+        #TODO: find version that works for all multifusion cases
+        symmetricbraiding && @timedtestset "Trace and contraction" begin
             t1 = rand(ComplexF64, V1 ⊗ V2 ⊗ V3)
             t2 = rand(ComplexF64, V2' ⊗ V4 ⊗ V1')
             t3 = t1 ⊗ t2
@@ -326,14 +354,14 @@ for V in spacelist
                 @test HrA12array ≈ convert(Array, HrA12)
             end
         end
-        @timedtestset "Index flipping: test flipping inverse" begin
+        (BraidingStyle(I) isa HasBraiding) && @timedtestset "Index flipping: test flipping inverse" begin
             t = rand(ComplexF64, V1 ⊗ V1' ← V1' ⊗ V1)
             for i in 1:4
                 @test t ≈ flip(flip(t, i), i; inv = true)
                 @test t ≈ flip(flip(t, i; inv = true), i)
             end
         end
-        @timedtestset "Index flipping: test via explicit flip" begin
+        symmetricbraiding && @timedtestset "Index flipping: test via explicit flip" begin
             t = rand(ComplexF64, V1 ⊗ V1' ← V1' ⊗ V1)
             F1 = unitary(flip(V1), V1)
 
@@ -346,7 +374,7 @@ for V in spacelist
             @tensor tf[a, b; c, d] := conj(F1[d, d']) * t[a, b; c, d']
             @test twist!(flip(t, 4), 4) ≈ tf
         end
-        @timedtestset "Index flipping: test via contraction" begin
+        symmetricbraiding && @timedtestset "Index flipping: test via contraction" begin
             t1 = rand(ComplexF64, V1 ⊗ V2 ⊗ V3 ← V4)
             t2 = rand(ComplexF64, V2' ⊗ V5 ← V4' ⊗ V1)
             @tensor ta[a, b] := t1[x, y, a, z] * t2[y, b, z, x]
@@ -360,7 +388,7 @@ for V in spacelist
         end
         @timedtestset "Multiplication of isometries: test properties" begin
             W2 = V4 ⊗ V5
-            W1 = W2 ⊗ (oneunit(V1) ⊕ oneunit(V1))
+            W1 = W2 ⊗ (unitspace(V1) ⊕ unitspace(V1))
             for T in (Float64, ComplexF64)
                 t1 = randisometry(T, W1, W2)
                 t2 = randisometry(T, W2 ← W2)
@@ -426,7 +454,7 @@ for V in spacelist
             end
         end
         @timedtestset "diag/diagm" begin
-            W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
+            W = V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5
             t = randn(ComplexF64, W)
             d = LinearAlgebra.diag(t)
             D = LinearAlgebra.diagm(codomain(t), domain(t), d)
@@ -505,8 +533,13 @@ for V in spacelist
         end
         @timedtestset "Tensor product: test via norm preservation" begin
             for T in (Float32, ComplexF64)
-                t1 = rand(T, V2 ⊗ V3 ⊗ V1, V1 ⊗ V2)
-                t2 = rand(T, V2 ⊗ V1 ⊗ V3, V1 ⊗ V1)
+                if UnitStyle(I) isa SimpleUnit || !isempty(blocksectors(V2 ⊗ V1))
+                    t1 = rand(T, V2 ⊗ V3 ⊗ V1, V1 ⊗ V2)
+                    t2 = rand(T, V2 ⊗ V1 ⊗ V3, V1 ⊗ V1)
+                else
+                    t1 = rand(T, V3 ⊗ V4 ⊗ V5, V1 ⊗ V2)
+                    t2 = rand(T, V5' ⊗ V4' ⊗ V3', V2' ⊗ V1')
+                end
                 t = @constinferred (t1 ⊗ t2)
                 @test norm(t) ≈ norm(t1) * norm(t2)
             end
@@ -528,7 +561,7 @@ for V in spacelist
                 end
             end
         end
-        @timedtestset "Tensor product: test via tensor contraction" begin
+        symmetricbraiding && @timedtestset "Tensor product: test via tensor contraction" begin
             for T in (Float32, ComplexF64)
                 t1 = rand(T, V2 ⊗ V3 ⊗ V1)
                 t2 = rand(T, V2 ⊗ V1 ⊗ V3)
@@ -537,10 +570,15 @@ for V in spacelist
                 @test t ≈ t′
             end
         end
-        @timedtestset "Tensor absorpsion" begin
+        @timedtestset "Tensor absorption" begin
             # absorbing small into large
-            t1 = zeros(V1 ⊕ V1, V2 ⊗ V3)
-            t2 = rand(V1, V2 ⊗ V3)
+            if UnitStyle(I) isa SimpleUnit || !isempty(blocksectors(V2 ⊗ V3))
+                t1 = zeros(V1 ⊕ V1, V2 ⊗ V3)
+                t2 = rand(V1, V2 ⊗ V3)
+            else
+                t1 = zeros(V1 ⊕ V2, V3 ⊗ V4 ⊗ V5)
+                t2 = rand(V1, V3 ⊗ V4 ⊗ V5)
+            end
             t3 = @constinferred absorb(t1, t2)
             @test norm(t3) ≈ norm(t2)
             @test norm(t1) == 0
@@ -549,8 +587,13 @@ for V in spacelist
             @test t3 ≈ t4
 
             # absorbing large into small
-            t1 = rand(V1 ⊕ V1, V2 ⊗ V3)
-            t2 = zeros(V1, V2 ⊗ V3)
+            if UnitStyle(I) isa SimpleUnit || !isempty(blocksectors(V2 ⊗ V3))
+                t1 = rand(V1 ⊕ V1, V2 ⊗ V3)
+                t2 = zeros(V1, V2 ⊗ V3)
+            else
+                t1 = rand(V1 ⊕ V2, V3 ⊗ V4 ⊗ V5)
+                t2 = zeros(V1, V3 ⊗ V4 ⊗ V5)
+            end
             t3 = @constinferred absorb(t2, t1)
             @test norm(t3) < norm(t1)
             @test norm(t2) == 0