diff --git a/Project.toml b/Project.toml
index a542085..640d3b2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,17 +5,23 @@ version = "0.1.0"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
-JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+BlockTensorKit = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f"
+TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 
 [compat]
 Aqua = "0.8.9"
 Artifacts = "1.10, 1"
-JSON3 = "1.14.1"
+BlockTensorKit = "0.1.10"
+DelimitedFiles = "1.9"
 SafeTestsets = "0.1"
-TensorKitSectors = "0.1.2"
+TensorKit = "0.14.9"
+TensorKitSectors = "0.1.7"
 Test = "1.10"
 TestExtras = "0.3"
+TupleTools = "1.1"
 julia = "1.10"
 
 [extras]
diff --git a/src/MultiTensorKit.jl b/src/MultiTensorKit.jl
index 34d66bc..9f59e4b 100644
--- a/src/MultiTensorKit.jl
+++ b/src/MultiTensorKit.jl
@@ -1,11 +1,21 @@
 module MultiTensorKit
 
 export BimoduleSector, A4Object
+export leftoneunit, rightoneunit
 
-using JSON3
+using DelimitedFiles
 using Artifacts
 using TensorKitSectors
 
+using TupleTools
+using TupleTools: insertafter
+
+using BlockTensorKit
+import BlockTensorKit: SumSpace
+
+using TensorKit
+import TensorKit: hasblock, dim
+
 include("bimodulesector.jl")
 
 end
diff --git a/src/bimodulesector.jl b/src/bimodulesector.jl
index fd52d31..7e6197d 100644
--- a/src/bimodulesector.jl
+++ b/src/bimodulesector.jl
@@ -2,6 +2,12 @@ struct BimoduleSector{Name} <: Sector
     i::Int
     j::Int
     label::Int
+    function BimoduleSector{:A4}(i::Int, j::Int, label::Int)
+        i <= size(BimoduleSector{:A4}) && j <= size(BimoduleSector{:A4}) ||
+            throw(DomainError("object outside the matrix A4"))
+        return label <= _numlabels(BimoduleSector{:A4}, i, j) ? new{:A4}(i, j, label) :
+               throw(DomainError("label outside category A4($i, $j)"))
+    end
 end
 BimoduleSector{Name}(data::NTuple{3,Int}) where {Name} = BimoduleSector{Name}(data...)
 const A4Object = BimoduleSector{:A4}
@@ -16,14 +22,16 @@ function Base.convert(::Type{BimoduleSector{Name}}, d::NTuple{3,Int}) where {Nam
     return BimoduleSector{Name}(d...)
 end
 
-Base.IteratorSize(::Type{SectorValues{<:BimoduleSector}}) = Base.SizeUnknown()
+Base.size(::Type{A4Object}) = 7
+
+Base.IteratorSize(::Type{<:SectorValues{<:BimoduleSector}}) = Base.SizeUnknown()
 
 # TODO: generalize?
-# TODO: numlabels?
 function Base.iterate(iter::SectorValues{A4Object}, (I, label)=(1, 1))
-    I > 12 * 12 && return nothing
-    i, j = CartesianIndices((12, 12))[I].I
-    maxlabel = numlabels(A4Object, i, j)
+    s = size(A4Object)
+    I > s * s && return nothing
+    i, j = CartesianIndices((s, s))[I].I
+    maxlabel = _numlabels(A4Object, i, j)
     return if label > maxlabel
         iterate(iter, (I + 1, 1))
     else
@@ -31,8 +39,14 @@ function Base.iterate(iter::SectorValues{A4Object}, (I, label)=(1, 1))
     end
 end
 
+function Base.length(::SectorValues{A4Object})
+    s = size(A4Object)
+    return sum(_numlabels(A4Object, i, j) for i in 1:s, j in 1:s)
+end
+
 TensorKitSectors.FusionStyle(::Type{A4Object}) = GenericFusion()
 TensorKitSectors.BraidingStyle(::Type{A4Object}) = NoBraiding()
+TensorKitSectors.sectorscalartype(::Type{A4Object}) = ComplexF64
 
 function TensorKitSectors.:⊗(a::A4Object, b::A4Object)
     @assert a.j == b.i
@@ -42,30 +56,35 @@ function TensorKitSectors.:⊗(a::A4Object, b::A4Object)
                     if (a_l == a.label && b_l == b.label)]
 end
 
-function _numlabels(::Type{A4Object}, i, j)
-    return Ncache = _get_Ncache(A4Object)
+function _numlabels(::Type{T}, i, j) where {T<:BimoduleSector}
+    return length(_get_dual_cache(T)[2][i, j])
 end
 
+# User-friendly functions
+# -------------------
+#TODO: add functions to identify categories
+
 # Data from files
 # ---------------
 const artifact_path = joinpath(artifact"fusiondata", "MultiTensorKit.jl-data-v0.1.5")
 
 function extract_Nsymbol(::Type{A4Object})
-    filename = joinpath(artifact_path, "A4", "Nsymbol.json")
-    isfile(filename) || throw(LoadError(filename, 0, "Nsymbol file not found for $Name"))
-    json_string = read(filename, String)
-    Narray = copy(JSON3.read(json_string))
-    return map(reshape(Narray, 12, 12, 12)) do x
-        y = Dict{NTuple{3,Int},Int}()
-        for (k, v) in x
-            a, b, c = parse.(Int, split(string(k)[2:(end - 1)], ", "))
-            y[(a, b, c)] = v
-        end
-        return y
+    filename = joinpath(artifact_path, "A4", "Nsymbol.txt")
+    isfile(filename) || throw(LoadError(filename, 0, "Nsymbol file not found for A4"))
+    Narray = readdlm(filename) # matrix with 7 columns
+
+    data_dict = Dict{NTuple{3,Int},Dict{NTuple{3,Int},Int}}()
+    for row in eachrow(Narray)
+        i, j, k, a, b, c, N = Int.(@view(row[1:size(A4Object)]))
+        colordict = get!(data_dict, (i, j, k), Dict{NTuple{3,Int},Int}())
+        push!(colordict, (a, b, c) => N)
     end
+
+    return data_dict
 end
 
-const Ncache = IdDict{Type{<:BimoduleSector},Array{Dict{NTuple{3,Int},Int},3}}()
+const Ncache = IdDict{Type{<:BimoduleSector},
+                      Dict{NTuple{3,Int},Dict{NTuple{3,Int},Int}}}()
 
 function _get_Ncache(::Type{T}) where {T<:BimoduleSector}
     global Ncache
@@ -82,8 +101,7 @@ function TensorKitSectors.Nsymbol(a::I, b::I, c::I) where {I<:A4Object}
     return get(_get_Ncache(I)[i, j, k], (a.label, b.label, c.label), 0)
 end
 
-# TODO: can we define dual for modules?
-const Dualcache = IdDict{Type{<:BimoduleSector},Vector{Tuple{Int,Vector{Int}}}}()
+const Dualcache = IdDict{Type{<:BimoduleSector},Tuple{Vector{Int64},Matrix{Vector{Int64}}}}()
 
 function _get_dual_cache(::Type{T}) where {T<:BimoduleSector}
     global Dualcache
@@ -94,75 +112,121 @@ end
 
 function extract_dual(::Type{A4Object})
     N = _get_Ncache(A4Object)
-    ncats = size(N, 1)
-    return map(1:ncats) do i
+    ncats = size(A4Object)
+    Is = zeros(Int, ncats)
+
+    map(1:ncats) do i
         Niii = N[i, i, i]
-        nobj = maximum(first, keys(Niii))
-
-        # find identity object:
-        # I x I -> I, a x I -> a, I x a -> a
-        I = findfirst(1:nobj) do u
-            get(Niii, (u, u, u), 0) == 1 || return false
-            for j in 1:nobj
-                get(Niii, (j, u, j), 0) == 1 || return false
-                get(Niii, (u, j, j), 0) == 1 || return false
+        nobji = maximum(first, keys(N[i, i, i]))
+        # want to return a leftone and rightone for each entry in multifusion cat
+        # leftone/rightone needs to at least be the unit object within a fusion cat
+        Is[i] = findfirst(1:nobji) do a
+            get(Niii, (a, a, a), 0) == 1 || return false # I x I -> I
+            for othera in 1:nobji
+                get(Niii, (othera, a, othera), 0) == 1 || return false # a x I -> a
+                get(Niii, (a, othera, othera), 0) == 1 || return false # I x a -> a
+            end
+
+            # check leftone
+            map(1:ncats) do j
+                nobjj = maximum(first, keys(N[j, j, j]))
+                for b in 1:nobjj
+                    get(N[i, j, j], (a, b, b), 0) == 1 || return false # I = leftone(b)
+                end
+            end
+
+            # check rightone
+            map(1:ncats) do k
+                nobjk = maximum(first, keys(N[k, k, k]))
+                for c in 1:nobjk
+                    get(N[k, i, k], (c, a, c), 0) == 1 || return false # I = rightone(c)
+                end
             end
             return true
         end
+    end
+
+    allduals = Matrix{Vector{Int}}(undef, ncats, ncats) # ncats square matrix of vectors
+    for i in 1:ncats
+        nobji = maximum(first, keys(N[i, i, i]))
+        for j in 1:ncats
+            allduals[i, j] = Int[]
 
-        # find duals
-        # a x abar -> I
-        duals = map(1:nobj) do j
-            return findfirst(1:nobj) do k
-                return get(Niii, (j, k, I), 0) == 1
+            nobjj = maximum(first, keys(N[j, j, j]))
+            # the nested vectors contain the duals of the objects in 𝒞_ij, which are in C_ji 
+            Niji = N[i, j, i] # 𝒞_ij x 𝒞_ji -> C_ii
+            Njij = N[j, i, j] # 𝒞_ji x 𝒞_ij -> C_jj
+            for i_ob in 1:nobji, j_ob in 1:nobjj
+                get(Niji, (i_ob, j_ob, Is[i]), 0) == 1 || continue # leftone(c_ij) ∈ c_ij x c_ji
+                get(Njij, (j_ob, i_ob, Is[j]), 0) == 1 || continue # rightone(c_ij) ∈ c_ji x c_ij
+                push!(allduals[i, j], j_ob)
             end
         end
-
-        return I, duals
     end
+    return Is, allduals
 end
 
 function Base.one(a::BimoduleSector)
-    a.i == a.j || error("don't know how to define one for modules")
-    return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[a.i][1])
+    a.i == a.j || error("unit object for module categories is ill-defined")
+    return typeof(a)(a.i, a.i, _get_dual_cache(typeof(a))[1][a.i])
+end
+
+function TensorKitSectors.leftone(a::BimoduleSector)
+    return typeof(a)(a.i, a.i, _get_dual_cache(typeof(a))[1][a.i])
+end
+
+function TensorKitSectors.rightone(a::BimoduleSector)
+    return typeof(a)(a.j, a.j, _get_dual_cache(typeof(a))[1][a.j])
 end
 
 function Base.conj(a::BimoduleSector)
-    a.i == a.j || error("don't know how to define dual for modules")
-    return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[a.i][2][a.label])
+    return typeof(a)(a.j, a.i, _get_dual_cache(typeof(a))[2][a.i, a.j][a.label])
 end
 
 function extract_Fsymbol(::Type{A4Object})
-    return mapreduce((x, y) -> cat(x, y; dims=1), 1:12) do i
-        filename = joinpath(artifact_path, "A4", "Fsymbol_$i.json")
-        @assert isfile(filename) "cannot find $filename"
-        json_string = read(filename, String)
-        Farray_part = copy(JSON3.read(json_string))
-        return map(enumerate(Farray_part[Symbol(i)])) do (I, x)
-            j, k, l = Tuple(CartesianIndices((12, 12, 12))[I])
-            y = Dict{NTuple{6,Int},Array{ComplexF64,4}}()
-            for (key, v) in x
-                a, b, c, d, e, f = parse.(Int, split(string(key)[2:(end - 1)], ", "))
-                a_ob, b_ob, c_ob, d_ob, e_ob, f_ob = A4Object.(((i, j, a), (j, k, b),
-                                                                (k, l, c), (i, l, d),
-                                                                (i, k, e), (j, l, f)))
-                result = Array{ComplexF64,4}(undef,
-                                             (Nsymbol(a_ob, b_ob, e_ob),
-                                              Nsymbol(e_ob, c_ob, d_ob),
-                                              Nsymbol(b_ob, c_ob, f_ob),
-                                              Nsymbol(a_ob, f_ob, d_ob)))
-                map!(result, reshape(v, size(result))) do cmplxdict
-                    return complex(cmplxdict[:re], cmplxdict[:im])
-                end
-
-                y[(a, b, c, d, e, f)] = result
+    result = Dict{NTuple{4,Int},Dict{NTuple{6,Int},Array{ComplexF64,4}}}()
+    filename = joinpath(artifact_path, "A4", "Fsymbol.txt")
+    @assert isfile(filename) "cannot find $filename"
+    Farray = readdlm(filename)
+    for ((i, j, k, l), colordict) in convert_Fs(Farray)
+        result[(i, j, k, l)] = Dict{NTuple{6,Int},Array{ComplexF64,4}}()
+        for ((a, b, c, d, e, f), Fvals) in colordict
+            a_ob, b_ob, c_ob, d_ob, e_ob, f_ob = A4Object.(((i, j, a), (j, k, b),
+                                                            (k, l, c), (i, l, d),
+                                                            (i, k, e), (j, l, f)))
+            result[(i, j, k, l)][(a, b, c, d, e, f)] = zeros(ComplexF64,
+                                                             Nsymbol(a_ob, b_ob, e_ob),
+                                                             Nsymbol(e_ob, c_ob, d_ob),
+                                                             Nsymbol(b_ob, c_ob, f_ob),
+                                                             Nsymbol(a_ob, f_ob, d_ob))
+            for (I, v) in Fvals
+                result[(i, j, k, l)][(a, b, c, d, e, f)][I] = v
             end
         end
     end
+    return result
+end
+
+function convert_Fs(Farray_part::Matrix{Float64}) # Farray_part is a matrix with 16 columns
+    data_dict = Dict{NTuple{4,Int},
+                     Dict{NTuple{6,Int},Vector{Pair{CartesianIndex{4},ComplexF64}}}}()
+    # want to make a Dict with keys (i,j,k,l) and vals 
+    # a Dict with keys (a,b,c,d,e,f) and vals 
+    # a pair of (mu, nu, rho, sigma) and the F value
+    for row in eachrow(Farray_part)
+        i, j, k, l, a, b, c, d, e, f, mu, nu, rho, sigma = Int.(@view(row[1:14]))
+        v = complex(row[15], row[16])
+        colordict = get!(data_dict, (i, j, k, l),
+                         Dict{NTuple{6,Int},Vector{Pair{CartesianIndex{4},ComplexF64}}}())
+        Fdict = get!(colordict, (a, b, c, d, e, f),
+                     Vector{Pair{CartesianIndex{4},ComplexF64}}())
+        push!(Fdict, CartesianIndex(mu, nu, rho, sigma) => v)
+    end
+    return data_dict
 end
 
 const Fcache = IdDict{Type{<:BimoduleSector},
-                      Array{Dict{NTuple{6,Int},Array{ComplexF64,4}},4}}()
+                      Dict{NTuple{4,Int64},Dict{NTuple{6,Int64},Array{ComplexF64,4}}}}()
 
 function _get_Fcache(::Type{T}) where {T<:BimoduleSector}
     global Fcache
@@ -173,16 +237,148 @@ end
 
 function TensorKitSectors.Fsymbol(a::I, b::I, c::I, d::I, e::I,
                                   f::I) where {I<:A4Object}
-    # TODO: should this error or return 0?
-    (a.j == b.i && b.j == c.i && a.i == d.i && c.j == d.j &&
-     a.i == e.i && b.j == e.j && f.i == a.j && f.j == c.j) ||
-        throw(ArgumentError("invalid fusion channel"))
+    # required to keep track of multiplicities where F-move is partially unallowed
+    # also deals with invalid fusion channels
+    Nabe = Nsymbol(a, b, e)
+    Necd = Nsymbol(e, c, d)
+    Nbcf = Nsymbol(b, c, f)
+    Nafd = Nsymbol(a, f, d)
+
+    zero_array = zeros(sectorscalartype(I), Nabe, Necd, Nbcf, Nafd)
+    Nabe > 0 && Necd > 0 && Nbcf > 0 && Nafd > 0 ||
+        return zero_array
 
     i, j, k, l = a.i, a.j, b.j, c.j
-    return get(_get_Fcache(I)[i, j, k, l],
-               (a.label, b.label, c.label, d.label, e.label, f.label)) do
-        return zeros(sectorscalartype(A4Object),
-                     (Nsymbol(a, b, e), Nsymbol(e, c, d), Nsymbol(b, c, f),
-                      Nsymbol(a, f, d)))
+    colordict = _get_Fcache(I)[i, j, k, l]
+    return get!(colordict, (a.label, b.label, c.label, d.label, e.label, f.label), zero_array)
+end
+
+# interface with TensorKit where necessary
+#-----------------------------------------
+
+function TensorKit.blocksectors(W::TensorMapSpace{S,N₁,N₂}) where
+         {S<:Union{Vect[A4Object],
+                   SumSpace{Vect[A4Object]}},N₁,N₂}
+    codom = codomain(W)
+    dom = domain(W)
+    if N₁ == 0 && N₂ == 0 # 0x0-dimensional TensorMap is just a scalar, return all units
+        # this is a problem in full contractions where the coloring outside is 𝒞
+        return NTuple{size(A4Object),A4Object}(one(A4Object(i, i, 1))
+                                               for i in 1:size(A4Object)) # have to return all units b/c no info on W in this case
+    elseif N₁ == 0
+        @assert N₂ != 0 "one of Type A4Object doesn't exist"
+        return filter!(isone, collect(blocksectors(dom)))
+    elseif N₂ == 0
+        @assert N₁ != 0 "one of Type A4Object doesn't exist"
+        return filter!(isone, collect(blocksectors(codom)))
+    elseif N₂ <= N₁ # keep intersection
+        return filter!(c -> hasblock(codom, c), collect(blocksectors(dom)))
+    else
+        return filter!(c -> hasblock(dom, c), collect(blocksectors(codom)))
+    end
+end
+
+#TODO: is this needed?
+# function TensorKit.scalar(t::AbstractTensorMap{T,S,0,0}) where {T,
+#                                                                 S<:GradedSpace{<:BimoduleSector}}
+#     _vector = findall(!iszero, t.data) # should have 0 or 1 elements, since only one of the blocks could be non-zero
+#     if isempty(_vector)
+#         return zero(scalartype(t))
+#     end
+#     unit = one(A4Object(only(_vector), only(_vector), 1))
+#     return only(block(t, unit))
+# end
+
+# TODO: definition for zero of GradedSpace?
+
+function dim(V::GradedSpace{<:BimoduleSector})
+    T = Base.promote_op(*, Int, real(sectorscalartype(sectortype(V))))
+    return reduce(+, dim(V, c) * dim(c) for c in sectors(V); init=zero(T))
+end
+
+# limited oneunit 
+function Base.oneunit(S::GradedSpace{<:BimoduleSector})
+    allequal(a.i for a in sectors(S)) && allequal(a.j for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S are not all equal"))
+    first(sectors(S)).i == first(sectors(S)).j ||
+        throw(ArgumentError("sectors of $S are non-diagonal"))
+    sector = one(first(sectors(S)))
+    return spacetype(S)(sector => 1)
+end
+
+function Base.oneunit(S::SumSpace{<:GradedSpace{<:BimoduleSector}})
+    @assert !isempty(S) "Cannot determine type of empty space"
+    return SumSpace(oneunit(first(S.spaces)))
+end
+
+# oneunit for spaces whose elements all belong to the same sector
+function rightoneunit(S::GradedSpace{<:BimoduleSector})
+    allequal(a.j for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S do not have the same rightone"))
+
+    allequal(a.i for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S are not all equal"))
+
+    sector = rightone(first(sectors(S)))
+    return spacetype(S)(sector => 1)
+end
+
+function rightoneunit(S::SumSpace{<:GradedSpace{<:BimoduleSector}})
+    @assert !isempty(S) "Cannot determine type of empty space"
+    return SumSpace(rightoneunit(first(S.spaces)))
+end
+
+function leftoneunit(S::GradedSpace{<:BimoduleSector})
+    allequal(a.i for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S do not have the same leftone"))
+
+    allequal(a.j for a in sectors(S)) ||
+        throw(ArgumentError("sectors of $S are not all equal"))
+
+    sector = leftone(first(sectors(S)))
+    return spacetype(S)(sector => 1)
+end
+
+function leftoneunit(S::SumSpace{<:GradedSpace{<:BimoduleSector}})
+    @assert !isempty(S) "Cannot determine type of empty space"
+    return SumSpace(leftoneunit(first(S.spaces)))
+end
+
+function TensorKit.insertrightunit(P::ProductSpace{V,N}, ::Val{i};
+                                   conj::Bool=false,
+                                   dual::Bool=false) where {i,V<:GradedSpace{I},N} where {I<:BimoduleSector}
+    i > N && error("cannot insert a sensible right unit onto $P at index $(i+1)")
+    # possible change to rightone of correct space for N = 0
+    u = N > 0 ? rightoneunit(P[i]) : error("no unit object in $P")
+    if dual
+        u = TensorKit.dual(u)
     end
+    if conj
+        u = TensorKit.conj(u)
+    end
+    return ProductSpace(TupleTools.insertafter(P.spaces, i, (u,)))
+end
+
+# possible TODO: overwrite defaults at level of HomSpace and TensorMap?
+function TensorKit.insertleftunit(P::ProductSpace{V,N}, ::Val{i}; # want no defaults?
+                                  conj::Bool=false,
+                                  dual::Bool=false) where {i,V<:GradedSpace{I},N} where {I<:BimoduleSector}
+    i > N && error("cannot insert a sensible left unit onto $P at index $i") # do we want this to error in the diagonal case?
+    u = N > 0 ? leftoneunit(P[i]) : error("no unit object in $P")
+    if dual
+        u = TensorKit.dual(u)
+    end
+    if conj
+        u = TensorKit.conj(u)
+    end
+    return ProductSpace(TupleTools.insertafter(P.spaces, i - 1, (u,)))
+end
+
+function TensorKit.scalar(t::AbstractTensorMap{T,S,0,0}) where {T,
+                                                                S<:GradedSpace{<:BimoduleSector}}
+    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/test_A4.jl b/test/test_A4.jl
index 696ca5b..6e54328 100644
--- a/test/test_A4.jl
+++ b/test/test_A4.jl
@@ -1,17 +1,18 @@
 using MultiTensorKit
-using TensorKitSectors
+using TensorKitSectors, TensorKit
 using Test, TestExtras
 
 I = A4Object
+Istr = TensorKitSectors.type_repr(I)
+r = size(I)
 
-@testset "Basic type properties" begin
-    Istr = TensorKitSectors.type_repr(I)
+@testset "$Istr Basic type properties" verbose = true begin
     @test eval(Meta.parse(sprint(show, I))) == I
     @test eval(Meta.parse(TensorKitSectors.type_repr(I))) == I
 end
 
-@testset "Fusion Category $i" for i in 1:12
-    objects = A4Object.(i, i, MultiTensorKit._get_dual_cache(I)[i][2])
+@testset "$Istr Fusion Category $i" for i in 1:r
+    objects = I.(i, i, MultiTensorKit._get_dual_cache(I)[2][i, i])
 
     @testset "Basic properties" begin
         s = rand(objects, 3)
@@ -45,10 +46,262 @@ end
             end
         end
     end
+end
+
+for i in 1:r, j in 1:r
+    i != j || continue # skip if fusion category
+    @testset "$Istr right module category $i, $j" begin
+        mod_objects = I.(i, j, MultiTensorKit._get_dual_cache(I)[2][i, j])
+        fusion_objects = I.(j, j, MultiTensorKit._get_dual_cache(I)[2][j, j])
+
+        @testset "Unitarity of module F-move $i, $j" begin # written as MxD, but with i<->j checks MopxC
+            for A in mod_objects, α in fusion_objects, β in fusion_objects
+                for B in ⊗(A, α, β)
+                    Cs = collect(intersect(⊗(A, α), map(dual, ⊗(β, dual(B)))))
+                    γs = collect(intersect(⊗(α, β), map(dual, ⊗(dual(B), A))))
+                    Fblocks = Vector{Any}()
+                    for C in Cs
+                        for γ in γs
+                            Fs = Fsymbol(A, α, β, B, C, γ)
+                            push!(Fblocks,
+                                  reshape(Fs,
+                                          (size(Fs, 1) * size(Fs, 2),
+                                           size(Fs, 3) * size(Fs, 4))))
+                        end
+                    end
+                    F = hvcat(length(γs), Fblocks...)
+                    @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12)
+                end
+            end
+        end
+    end
+end
+
+for i in 1:7, j in 1:7 # C x M (or D x Mop with i<->j)
+    i != j || continue # skip if fusion category
+    @testset "$Istr left module category $i, $j unitarity check" begin
+        mod_objects = I.(i, j, MultiTensorKit._get_dual_cache(I)[2][i, j])
+        fusion_objects = I.(i, i, MultiTensorKit._get_dual_cache(I)[2][i, i])
+
+        @testset "Unitarity of left module F-move" begin
+            for a in fusion_objects, b in fusion_objects, A in mod_objects  # written for M as left C-module category
+                for B in ⊗(a, b, A)
+                    cs = collect(intersect(⊗(a, b), map(dual, ⊗(A, dual(B))))) # equivalent of es
+                    Cs = collect(intersect(⊗(b, A), map(dual, ⊗(dual(B), a)))) # equivalent of fs
+                    Fblocks = Vector{Any}()
+                    for c in cs
+                        for C in Cs
+                            Fs = Fsymbol(a, b, A, B, c, C)
+                            push!(Fblocks,
+                                  reshape(Fs,
+                                          (size(Fs, 1) * size(Fs, 2),
+                                           size(Fs, 3) * size(Fs, 4))))
+                        end
+                    end
+                    F = hvcat(length(Cs), Fblocks...)
+                    @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12)
+                end
+            end
+        end
+    end
+end
+
+for i in 1:7, j in 1:7 # bimodule check unitarity (C x M x D or D x Mop x C)
+    i != j || continue # skip if fusion category
+    @testset "$Istr bimodule category $i, $j" begin
+        C_objects = I.(i, i, MultiTensorKit._get_dual_cache(I)[2][i, i])
+        mod_objects = I.(i, j, MultiTensorKit._get_dual_cache(I)[2][i, j])
+        D_objects = I.(j, j, MultiTensorKit._get_dual_cache(I)[2][j, j])
+
+        @testset "Unitarity of bimodule F-move" begin # written as CMD
+            for a in C_objects, A in mod_objects, α in D_objects
+                for B in ⊗(a, A, α)
+                    Cs = collect(intersect(⊗(a, A), map(dual, ⊗(α, dual(B))))) # equivalent of es
+                    Ds = collect(intersect(⊗(A, α), map(dual, ⊗(dual(B), a)))) # equivalent of fs
+                    Fblocks = Vector{Any}()
+                    for C in Cs
+                        for D in Ds
+                            Fs = Fsymbol(a, A, α, B, C, D)
+                            push!(Fblocks,
+                                  reshape(Fs,
+                                          (size(Fs, 1) * size(Fs, 2),
+                                           size(Fs, 3) * size(Fs, 4))))
+                        end
+                    end
+                    F = hvcat(length(Ds), Fblocks...)
+                    @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12)
+                end
+            end
+        end
+    end
+end
+
+for i in 1:7, j in 1:7 # M x Mop x M -> M (or Mop x M x Mop -> Mop)
+    i != j || continue # skip if not fusion category
+    @testset "$Istr Module category $i,$j and opposite $j,$i" begin
+        M_objects = I.(i, j, MultiTensorKit._get_dual_cache(I)[2][i, j])
+        Mop_objects = I.(j, i, MultiTensorKit._get_dual_cache(I)[2][j, i])
+        C_objects = I.(i, i, MultiTensorKit._get_dual_cache(I)[2][i, i])
+        D_objects = I.(j, j, MultiTensorKit._get_dual_cache(I)[2][j, j])
+
+        @testset "Unitarity of mixed module F-move I" begin # written for C x M x Mop -> C but also holds for D x Mop x M -> D
+            for α in C_objects, A in M_objects, Aop in Mop_objects
+                for β in ⊗(α, A, Aop)
+                    Cs = collect(intersect(⊗(α, A), map(dual, ⊗(Aop, dual(β))))) # equivalent of es
+                    γs = collect(intersect(⊗(A, Aop), map(dual, ⊗(dual(β), α)))) # equivalent of fs
+                    Fblocks = Vector{Any}()
+                    for C in Cs
+                        for γ in γs
+                            Fs = Fsymbol(α, A, Aop, β, C, γ)
+                            push!(Fblocks,
+                                  reshape(Fs,
+                                          (size(Fs, 1) * size(Fs, 2),
+                                           size(Fs, 3) * size(Fs, 4))))
+                        end
+                    end
+                    F = hvcat(length(γs), Fblocks...)
+                    @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12)
+                end
+            end
+        end
+
+        @testset "Unitarity of mixed module F-move II" begin # written for M x Mop x C -> C but also holds for Mop x M x D -> D
+            for A in M_objects, Aop in Mop_objects, α in C_objects
+                for β in ⊗(A, Aop, α)
+                    γs = collect(intersect(⊗(A, Aop), map(dual, ⊗(α, dual(β))))) # equivalent of es
+                    Bops = collect(intersect(⊗(Aop, α), map(dual, ⊗(dual(β), A)))) # equivalent of fs
+                    Fblocks = Vector{Any}()
+                    for γ in γs
+                        for Bop in Bops
+                            Fs = Fsymbol(A, Aop, α, β, γ, Bop)
+                            push!(Fblocks,
+                                  reshape(Fs,
+                                          (size(Fs, 1) * size(Fs, 2),
+                                           size(Fs, 3) * size(Fs, 4))))
+                        end
+                    end
+                    F = hvcat(length(Bops), Fblocks...)
+                    @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12)
+                end
+            end
+        end
+
+        @testset "Unitarity of mixed module F-move III" begin # written for Mop x C x M -> D, but also holds for M x D x Mop -> C
+            for Aop in Mop_objects, α in C_objects, A in M_objects
+                for a in ⊗(Aop, α, A)
+                    Bops = collect(intersect(⊗(Aop, α), map(dual, ⊗(A, dual(a))))) # equivalent of es
+                    Bs = collect(intersect(⊗(α, A), map(dual, ⊗(dual(a), Aop)))) # equivalent of fs
+                    Fblocks = Vector{Any}()
+                    for Bop in Bops
+                        for B in Bs
+                            Fs = Fsymbol(Aop, α, A, a, Bop, B)
+                            push!(Fblocks,
+                                  reshape(Fs,
+                                          (size(Fs, 1) * size(Fs, 2),
+                                           size(Fs, 3) * size(Fs, 4))))
+                        end
+                    end
+                    F = hvcat(length(Bs), Fblocks...)
+                    @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12)
+                end
+            end
+        end
+
+        @testset "Unitarity of pure module F-move" begin
+            for A in M_objects, Aop in Mop_objects, B in M_objects # written for M x Mop x M -> M but also holds for Mop x M x Mop -> Mop
+                for C in ⊗(A, Aop, B)
+                    cs = collect(intersect(⊗(A, Aop), map(dual, ⊗(B, dual(C))))) # equivalent of es
+                    γs = collect(intersect(⊗(Aop, B), map(dual, ⊗(dual(C), A)))) # equivalent of fs
+                    Fblocks = Vector{Any}()
+                    for c in cs
+                        for γ in γs
+                            Fs = Fsymbol(A, Aop, B, C, c, γ)
+                            push!(Fblocks,
+                                  reshape(Fs,
+                                          (size(Fs, 1) * size(Fs, 2),
+                                           size(Fs, 3) * size(Fs, 4))))
+                        end
+                    end
+                    F = hvcat(length(γs), Fblocks...)
+                    @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12)
+                end
+            end
+        end
+    end
+end
+
+@testset "Triangle equation" begin
+    objects = collect(values(I))
+    for a in objects, b in objects
+        a.j == b.i || continue # skip if not compatible
+        @test triangle_equation(a, b; atol=1e-12, rtol=1e-12)
+    end
+end
+
+@testset "$Istr Pentagon equation" begin
+    objects = collect(values(I))
+    for a in objects
+        for b in objects
+            a.j == b.i || continue # skip if not compatible
+            for c in objects
+                b.j == c.i || continue # skip if not compatible
+                for d in objects
+                    c.j == d.i || continue # skip if not compatible
+                    @test pentagon_equation(a, b, c, d; atol=1e-12, rtol=1e-12)
+                end
+            end
+        end
+    end
+end
+
+@testset "$Istr ($i, $j) units and duals" for i in 1:r, j in 1:r
+    Cij_obs = I.(i, j, MultiTensorKit._get_dual_cache(I)[2][i, j])
+
+    s = rand(Cij_obs)
+    @test eval(Meta.parse(sprint(show, s))) == s
+    @test @constinferred(hash(s)) == hash(deepcopy(s))
+    @test i == j ? isone(@constinferred(one(s))) :
+          (isone(@constinferred(leftone(s))) && isone(@constinferred(rightone(s))))
+    @constinferred dual(s)
+    @test dual(s) == I.(j, i, MultiTensorKit._get_dual_cache(I)[2][i, j][s.label])
+    @test dual(dual(s)) == s
+end
+
+@testset "$Istr ($i, $j) left and right units" for i in 1:r, j in 1:r
+    Cij_obs = I.(i, j, MultiTensorKit._get_dual_cache(I)[2][i, j])
+
+    s = rand(Cij_obs, 1)[1]
+    sp = Vect[I](s => 1)
+    W = sp ← sp
+    for T in (Float32, ComplexF64)
+        t = @constinferred rand(T, W)
+
+        for a in 1:2
+            tl = @constinferred insertleftunit(t, Val(a))
+            @test numind(tl) == numind(t) + 1
+            @test space(tl) == insertleftunit(space(t), a)
+            @test scalartype(tl) === T
+            @test t.data === tl.data
+            @test @constinferred(removeunit(tl, $(a))) == t
+
+            tr = @constinferred insertrightunit(t, Val(a))
+            @test numind(tr) == numind(t) + 1
+            @test space(tr) == insertrightunit(space(t), a)
+            @test scalartype(tr) === T
+            @test t.data === tr.data
+            @test @constinferred(removeunit(tr, $(a + 1))) == t
+        end
+
+        @test_throws ErrorException insertleftunit(t) # default should error here
+        @test insertrightunit(t) isa TensorMap
+        @test_throws ErrorException insertleftunit(t, numind(t) + 1) # same as default
+        @test_throws ErrorException insertrightunit(t, numind(t) + 1) # not same as default
 
-    @testset "Pentagon equation" begin
-        for a in objects, b in objects, c in objects, d in objects
-            @test pentagon_equation(a, b, c, d; atol=1e-12, rtol=1e-12)
+        t2 = @constinferred insertrightunit(t; copy=true)
+        @test t.data !== t2.data
+        for (c, b) in blocks(t)
+            @test b == block(t2, c)
         end
+        @test @constinferred(removeunit(t2, $(numind(t2)))) == t
     end
 end
diff --git a/test/test_aqua.jl b/test/test_aqua.jl
index 14e5b7c..f57c8ba 100644
--- a/test/test_aqua.jl
+++ b/test/test_aqua.jl
@@ -3,5 +3,5 @@ using Aqua: Aqua
 using Test: @testset
 
 @testset "Code quality (Aqua.jl)" begin
-  Aqua.test_all(MultiTensorKit)
+    Aqua.test_all(MultiTensorKit)
 end