diff --git a/Project.toml b/Project.toml
index fcbaa25..c11392b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FusionTensors"
 uuid = "e16ca583-1f51-4df0-8e12-57d32947d33e"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.3.2"
+version = "0.4.0"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -18,16 +18,16 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
 WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
 
 [compat]
-Accessors = "0.1.39"
-BlockArrays = "1.2.0"
-BlockSparseArrays = "0.4.0"
-GradedArrays = "0.2.1"
-HalfIntegers = "1.6.0"
-LRUCache = "1.6.1"
-LinearAlgebra = "1.10.0"
-Strided = "2.2.0"
-TensorAlgebra = "0.2.4"
-TensorProducts = "0.1.3"
-TypeParameterAccessors = "0.3.0"
+Accessors = "0.1.42"
+BlockArrays = "1.6"
+BlockSparseArrays = "0.7.4"
+GradedArrays = "0.4.7"
+HalfIntegers = "1.6"
+LRUCache = "1.6"
+LinearAlgebra = "1.10"
+Strided = "2.3"
+TensorAlgebra = "0.3.8"
+TensorProducts = "0.1.7"
+TypeParameterAccessors = "0.4"
 WignerSymbols = "2.0.0"
 julia = "1.10"
diff --git a/docs/Project.toml b/docs/Project.toml
index 3ba749f..5d0dd21 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -5,5 +5,5 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
 Documenter = "1.10.0"
-FusionTensors = "0.3.0"
+FusionTensors = "0.4"
 Literate = "2.20.1"
diff --git a/src/fusion_trees/clebsch_gordan_tensors.jl b/src/fusion_trees/clebsch_gordan_tensors.jl
index b7cc21b..37575d4 100644
--- a/src/fusion_trees/clebsch_gordan_tensors.jl
+++ b/src/fusion_trees/clebsch_gordan_tensors.jl
@@ -5,17 +5,18 @@
 using HalfIntegers: half
 using WignerSymbols: clebschgordan
 
-using GradedArrays: dual
-using GradedArrays.SymmetrySectors:
+using GradedArrays:
   AbelianStyle,
   AbstractSector,
   NotAbelianStyle,
-  SymmetryStyle,
   O2,
   SU,
+  SymmetryStyle,
+  dual,
   istrivial,
   quantum_dimension,
   sector_label,
+  sectors,
   trivial,
   zero_odd
 using TensorAlgebra: contract
@@ -66,7 +67,7 @@ function clebsch_gordan_tensor(s1::O2, s2::O2, s3::O2)
   d2 = quantum_dimension(s2)
   d3 = quantum_dimension(s3)
   cgt = zeros((d1, d2, d3))
-  s3 ∉ blocklabels(s1 ⊗ s2) && return cgt
+  s3 ∉ sectors(s1 ⊗ s2) && return cgt
 
   # adapted from TensorKit
   l1 = sector_label(s1)
diff --git a/src/fusion_trees/fusiontree.jl b/src/fusion_trees/fusiontree.jl
index c5a340a..461c8e2 100644
--- a/src/fusion_trees/fusiontree.jl
+++ b/src/fusion_trees/fusiontree.jl
@@ -3,14 +3,20 @@
 # TBD
 # compatibility with TensorKit conventions?
 
-using GradedArrays: GradedArrays, AbstractGradedUnitRange, flip, isdual, sector_type
-using GradedArrays.SymmetrySectors:
+using GradedArrays:
+  GradedArrays,
   ×,
+  AbstractGradedUnitRange,
   AbstractSector,
   SectorProduct,
-  SymmetrySectors,
   arguments,
+  flip,
+  flip_dual,
+  isdual,
   nsymbol,
+  sector_multiplicities,
+  sector_type,
+  sectors,
   to_gradedrange,
   trivial
 using TensorAlgebra: flatten_tuples
@@ -41,7 +47,7 @@ using TensorProducts: ⊗
 #
 #
 # The interface uses AbstractGradedArrays as input for interface simplicity
-# however only blocklabels are used and blocklengths are never read.
+# however only sectors are used and blocklengths are never read.
 
 struct SectorFusionTree{S,N,M}
   leaves::NTuple{N,S}  # TBD rename outer_sectors or leave_sectors?
@@ -98,8 +104,7 @@ function GradedArrays.flip(f::SectorFusionTree)
   )
 end
 
-# SymmetrySectors interface
-function SymmetrySectors.:×(f1::SectorFusionTree, f2::SectorFusionTree)
+function GradedArrays.:×(f1::SectorFusionTree, f2::SectorFusionTree)
   @assert arrows(f1) == arrows(f2)
   product_leaves = .×(leaves(f1), leaves(f2))
   product_root_sector = root_sector(f1) × root_sector(f2)
@@ -120,7 +125,7 @@ function SymmetrySectors.:×(f1::SectorFusionTree, f2::SectorFusionTree)
   )
 end
 
-function SymmetrySectors.arguments(f::SectorFusionTree{<:SectorProduct})
+function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct})
   transposed_indices = outer_multiplicity_split.(
     Base.tail(leaves(f)),
     branch_sectors(f),
@@ -144,11 +149,11 @@ function SymmetrySectors.arguments(f::SectorFusionTree{<:SectorProduct})
   )
 end
 
-function SymmetrySectors.arguments(f::SectorFusionTree{<:SectorProduct,0})
+function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct,0})
   return map(arg -> SectorFusionTree((), (), arg, (), ()), arguments(root_sector(f)))
 end
 
-function SymmetrySectors.arguments(f::SectorFusionTree{<:SectorProduct,1})
+function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct,1})
   arguments_root = arguments(root_sector(f))
   arguments_leave = arguments(only(leaves(f)))
   # use map(keys) to stay agnostic with respect to SectorProduct implementation
@@ -164,7 +169,7 @@ fusiontree_eltype(::Type{<:AbstractSector}) = Float64
 function build_trees(legs::Vararg{AbstractGradedUnitRange})
   # construct all authorized trees for each outer block in legs
   tree_arrows = isdual.(legs)
-  return mapreduce(vcat, Iterators.product(blocklabels.(legs)...)) do it
+  return mapreduce(vcat, Iterators.product(sectors.(flip_dual.(legs))...)) do it
     return build_trees(it, tree_arrows)
   end
 end
@@ -240,7 +245,9 @@ function fuse_next_sector(
   parent_tree::SectorFusionTree, branch_sector::AbstractSector, level_arrow::Bool
 )
   new_space = to_gradedrange(root_sector(parent_tree) ⊗ branch_sector)
-  return mapreduce(vcat, zip(blocklabels(new_space), blocklengths(new_space))) do (la, n)
+  return mapreduce(
+    vcat, zip(sectors(new_space), sector_multiplicities(new_space))
+  ) do (la, n)
     return [
       append_tree_leave(parent_tree, branch_sector, level_arrow, la, outer_mult) for
       outer_mult in 1:n
diff --git a/src/fusiontensor/array_cast.jl b/src/fusiontensor/array_cast.jl
index edfa619..cbc5891 100644
--- a/src/fusiontensor/array_cast.jl
+++ b/src/fusiontensor/array_cast.jl
@@ -3,8 +3,7 @@
 using BlockArrays: AbstractBlockArray, BlockedArray, blockedrange, blocklengths, findblock
 
 using BlockSparseArrays: BlockSparseArrays, BlockSparseArray
-using GradedArrays: AbstractGradedUnitRange, blocklabels
-using GradedArrays.SymmetrySectors: block_dimensions, quantum_dimension
+using GradedArrays: AbstractGradedUnitRange, quantum_dimension
 using TensorAlgebra: contract
 
 # =================================  High level interface  =================================
@@ -20,11 +19,13 @@ end
 function to_fusiontensor(
   array::AbstractArray,
   codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}},
-  domain_legs::Tuple{Vararg{AbstractGradedUnitRange}},
+  domain_legs::Tuple{Vararg{AbstractGradedUnitRange}};
+  atol::Real=0,
+  rtol::Real=rtoldefault(array),
 )
-  bounds = block_dimensions.((codomain_legs..., domain_legs...))
+  bounds = blocklengths.((codomain_legs..., domain_legs...))
   blockarray = BlockedArray(array, bounds...)
-  return to_fusiontensor(blockarray, codomain_legs, domain_legs)
+  return to_fusiontensor(blockarray, codomain_legs, domain_legs; atol, rtol)
 end
 
 rtoldefault(a::AbstractArray) = rtoldefault(eltype(a))
@@ -72,7 +73,7 @@ function to_fusiontensor_no_checknorm(
 end
 
 function to_array(ft::FusionTensor)
-  bounds = block_dimensions.((codomain_axes(ft)..., domain_axes(ft)...))
+  bounds = blocklengths.((codomain_axes(ft)..., domain_axes(ft)...))
   bsa = BlockSparseArray{eltype(ft)}(undef, blockedrange.(bounds))
   for (f1, f2) in keys(trees_block_mapping(ft))
     b = findblock(ft, f1, f2)
diff --git a/src/fusiontensor/base_interface.jl b/src/fusiontensor/base_interface.jl
index 3343d8b..e4ba44f 100644
--- a/src/fusiontensor/base_interface.jl
+++ b/src/fusiontensor/base_interface.jl
@@ -17,8 +17,6 @@ function Base.:*(left::FusionTensor, right::FusionTensor)
   return FusionTensor(new_data_matrix, codomain_axes(left), domain_axes(right))
 end
 
-Base.:+(ft::FusionTensor) = ft
-
 # tensor addition is a block data_matrix add.
 function Base.:+(left::FusionTensor, right::FusionTensor)
   checkaxes(axes(left), axes(right))
@@ -89,37 +87,34 @@ function Base.similar(ft::FusionTensor, T::Type)
 
   # some fusion trees have Float64 eltype: need compatible type
   @assert promote_type(T, fusiontree_eltype(sector_type(ft))) === T
-  mat = similar(data_matrix(ft), T)
-  initialize_allowed_sectors!(mat)
+  mat = initialize_data_matrix(T, codomain_axis(ft), domain_axis(ft))
   return FusionTensor(mat, axes(ft), trees_block_mapping(ft))
 end
 
-# trigger explicit error in TensorAlgebra.contract
-# TBD impose some convention? Remove?
 function Base.similar(
-  ft::FusionTensor,
-  T::Type,
-  new_axes::Tuple{AbstractGradedUnitRange,Vararg{AbstractGradedUnitRange}},
+  ::FusionTensor, ::Type, t::Tuple{AbstractGradedUnitRange,Vararg{AbstractGradedUnitRange}}
 )
-  throw(DimensionMismatch("Need bituple of axes"))
+  throw(MethodError(similar, (typeof(t),)))
 end
-function Base.similar(ft::FusionTensor, T::Type, new_axes::Tuple{})
-  throw(DimensionMismatch("Need bituple of axes"))
+function Base.similar(::FusionTensor, ::Type, ::Tuple{})
+  throw(MethodError(similar, (Tuple{},)))
 end
 
-function Base.similar(ft::FusionTensor, T::Type, new_axes::Tuple{<:Tuple,<:Tuple})
+function Base.similar(
+  ft::FusionTensor, ::Type{T}, new_axes::Tuple{<:Tuple,<:Tuple}
+) where {T}
   return similar(ft, T, tuplemortar(new_axes))
 end
-function Base.similar(::FusionTensor, T::Type, new_axes::BlockedTuple{2})
+function Base.similar(::FusionTensor, ::Type{T}, new_axes::BlockedTuple{2}) where {T}
   return FusionTensor(T, new_axes)
 end
 
 Base.show(io::IO, ft::FusionTensor) = print(io, "$(ndims(ft))-dim FusionTensor")
 
 function Base.show(io::IO, ::MIME"text/plain", ft::FusionTensor)
-  println(io, "$(ndims(ft))-dim FusionTensor with axes:")
+  print(io, "$(ndims(ft))-dim FusionTensor with axes:")
   for ax in axes(ft)
-    println(io, ax)
+    print(io, "\n", ax)
   end
   return nothing
 end
diff --git a/src/fusiontensor/fusiontensor.jl b/src/fusiontensor/fusiontensor.jl
index 64d937d..4d0b5a0 100644
--- a/src/fusiontensor/fusiontensor.jl
+++ b/src/fusiontensor/fusiontensor.jl
@@ -6,59 +6,80 @@ using BlockSparseArrays:
   AbstractBlockSparseMatrix, BlockSparseArray, eachblockstoredindex, to_block_indices
 using GradedArrays:
   AbstractGradedUnitRange,
-  blocklabels,
-  blockmergesort,
+  SectorProduct,
+  TrivialSector,
   dual,
   flip,
+  flip_dual,
   gradedrange,
   isdual,
-  map_blocklabels,
+  map_sectors,
+  sector_multiplicity,
   sector_type,
+  sectormergesort,
+  sectors,
   space_isequal
-using GradedArrays.SymmetrySectors: SectorProduct, TrivialSector
 using TensorAlgebra: BlockedTuple, tuplemortar
 using TensorProducts: tensor_product
+using TypeParameterAccessors: type_parameters
 
-struct FusionTensor{T,N,Axes,Mat,Mapping} <: AbstractArray{T,N}
+struct FusionTensor{T,N,Axes,Mat<:AbstractMatrix{T},Mapping} <: AbstractArray{T,N}
   data_matrix::Mat
   axes::Axes
   trees_block_mapping::Mapping
 
   # inner constructor to impose constraints on types
-  function FusionTensor(
-    mat::AbstractMatrix,
-    legs::BlockedTuple{2,<:Any,<:Tuple{Vararg{AbstractGradedUnitRange}}},
-    trees_block_mapping::Dict,
-  )
-    S = sector_type(axes(mat, 1))
-    @assert sector_type(axes(mat, 2)) === S
+  function FusionTensor{T,N,Axes,Mat,Mapping}(
+    mat, legs, trees_block_mapping
+  ) where {T,N,Axes,Mat,Mapping}
+    S = isempty(legs) ? TrivialSector : sector_type(first(legs))
     @assert keytype(trees_block_mapping) <:
       Tuple{<:SectorFusionTree{S},<:SectorFusionTree{S}}
     @assert all(sector_type.(Tuple(legs)) .=== S)
-    return new{
-      eltype(mat),length(legs),typeof(legs),typeof(mat),typeof(trees_block_mapping)
-    }(
-      mat, legs, trees_block_mapping
-    )
+    return new{T,N,Axes,Mat,Mapping}(mat, legs, trees_block_mapping)
   end
 end
 
+function FusionTensor(
+  mat::AbstractMatrix,
+  legs::BlockedTuple{2,<:Any,<:Tuple{Vararg{AbstractGradedUnitRange}}},
+  trees_block_mapping::Dict{<:Tuple{<:SectorFusionTree,<:SectorFusionTree}},
+)
+  return FusionTensor{
+    eltype(mat),length(legs),typeof(legs),typeof(mat),typeof(trees_block_mapping)
+  }(
+    mat, legs, trees_block_mapping
+  )
+end
+
 # getters
 data_matrix(ft::FusionTensor) = ft.data_matrix
 trees_block_mapping(ft::FusionTensor) = ft.trees_block_mapping
 
 # misc access
-codomain_axes(ft::FusionTensor) = first(blocks(axes(ft)))
-domain_axes(ft::FusionTensor) = last(blocks(axes(ft)))
+codomain_axes(ft::FusionTensor) = axes(ft)[Block(1)]
+domain_axes(ft::FusionTensor) = axes(ft)[Block(2)]
 ndims_codomain(ft::FusionTensor) = length(codomain_axes(ft))
 ndims_domain(ft::FusionTensor) = length(domain_axes(ft))
 
-matrix_size(ft::FusionTensor) = quantum_dimension.(axes(data_matrix(ft)))
-matrix_row_axis(ft::FusionTensor) = first(axes(data_matrix(ft)))
-matrix_column_axis(ft::FusionTensor) = last(axes(data_matrix(ft)))
+dummy_axis(ft::FusionTensor) = dummy_axis(sector_type(ft))
+dummy_axis(::Type{S}) where {S<:AbstractSector} = gradedrange([trivial(S) => 1])
+
+function codomain_axis(ft::FusionTensor)
+  if ndims_codomain(ft) == 0
+    return dummy_axis(ft)
+  end
+  return ⊗(codomain_axes(ft)...)
+end
+function domain_axis(ft::FusionTensor)
+  if ndims_domain(ft) == 0
+    return dummy_axis(ft)
+  end
+  return dual(⊗(dual.(domain_axes(ft))...))
+end
 function charge_block_size(ft::FusionTensor, f1::SectorFusionTree, f2::SectorFusionTree)
   b = Tuple(findblock(ft, f1, f2))
-  return ntuple(i -> Int(length(axes(ft)[i][b[i]])), ndims(ft))
+  return ntuple(i -> sector_multiplicity(axes(ft, i)[b[i]]), ndims(ft))
 end
 
 # GradedArrays interface
@@ -77,13 +98,13 @@ function BlockArrays.findblock(ft::FusionTensor, f1::SectorFusionTree, f2::Secto
   return Block(b1..., b2...)
 end
 # TBD move to GradedArrays? rename findfirst_sector?
-function find_sector_block(s::AbstractSector, l::AbstractGradedUnitRange)
-  return findfirst(==(s), blocklabels(l))
+function find_sector_block(s::AbstractSector, g::AbstractGradedUnitRange)
+  return findfirst(==(s), sectors(flip_dual(g)))
 end
 
 function sanitize_axes(raw_legs::Tuple{Vararg{AbstractGradedUnitRange}})
-  legs = promote_sectors(typeof(first(raw_legs)), raw_legs)
-  @assert all(check_unique_blocklabels.(legs))
+  legs = promote_sectors(raw_legs)
+  @assert all(check_unique_sectors.(legs))
   return legs
 end
 sanitize_axes(legs::BlockedTuple{2,(0, 0)}) = TrivialSector, legs
@@ -92,29 +113,22 @@ function sanitize_axes(raw_legs::BlockedTuple{2})
   return sector_type(first(flat_legs)), BlockedTuple(flat_legs, Val(blocklengths(raw_legs)))
 end
 
-function check_unique_blocklabels(g::AbstractGradedUnitRange)
-  return length(unique(blocklabels(g))) == blocklength(g)
+function check_unique_sectors(g::AbstractGradedUnitRange)
+  return length(unique(sectors(g))) == blocklength(g)
 end
 
-function promote_sectors(
-  ::Type{<:AbstractGradedUnitRange{LA}}, legs::Tuple{Vararg{AbstractGradedUnitRange{LA}}}
-) where {LA}  # nothing to do
-  return legs
-end
-
-function promote_sectors(
-  ::Type{<:AbstractGradedUnitRange}, legs::Tuple{Vararg{AbstractGradedUnitRange}}
-)
+promote_sectors(legs::NTuple{<:Any,<:AbstractGradedUnitRange}) = legs # nothing to do
+function promote_sectors(legs)
   T = promote_sector_type(legs)
   # fuse with trivial to insert all missing arguments inside each GradedAxis
-  # avoid depending on SymmetrySectors internals
+  # avoid depending on GradedArrays internals
   s0 = trivial(T)
-  return map_blocklabels.(s -> only(blocklabels(tensor_product(s0, s))), legs)
+  return map_sectors.(s -> only(sectors(to_gradedrange(tensor_product(s0, s)))), legs)
 end
 
 function promote_sector_type(legs)
   # fuse trivial sectors to produce unified type
-  # avoid depending on SymmetrySectors internals
+  # avoid depending on GradedArrays internals
   return sector_type(tensor_product(trivial.(legs)...))
 end
 
@@ -130,8 +144,6 @@ end
 function FusionTensor(mat::AbstractMatrix, legs::BlockedTuple{2})
   # init with empty data_matrix to construct trees_block_mapping
   ft = FusionTensor(eltype(mat), legs)
-  @assert space_isequal(matrix_row_axis(ft), axes(mat, 1))
-  @assert space_isequal(matrix_column_axis(ft), axes(mat, 2))
   for b in eachblockstoredindex(mat)
     @assert b in eachblockstoredindex(data_matrix(ft))  # check matrix block is allowed
     data_matrix(ft)[b] = mat[b]
@@ -150,9 +162,8 @@ end
 # empty matrix
 function FusionTensor(elt::Type, raw_legs::BlockedTuple{2})
   S, legs = sanitize_axes(raw_legs)
-
-  row_axis, codomain_trees_to_ranges = fuse_axes(S, first(blocks(legs)))
-  col_axis, domain_trees_to_ranges = flip_domain(fuse_axes(S, dual.(last(blocks(legs))))...)
+  row_axis, codomain_trees_to_ranges = fuse_axes(S, legs[Block(1)])
+  col_axis, domain_trees_to_ranges = flip_domain(fuse_axes(S, dual.(legs[Block(2)]))...)
 
   mat = initialize_data_matrix(elt, row_axis, col_axis)
   tree_to_block_mapping = intersect_codomain_domain(
@@ -181,17 +192,17 @@ function fusion_trees_external_multiplicities(
 end
 
 function block_fusion_trees_external_multiplicities(it::Tuple{Vararg{AbstractUnitRange}})
-  block_sectors = only.(blocklabels.(it))
-  block_mult = prod(length.(it))
+  block_sectors = only.(sectors.(flip_dual.(it)))
+  block_mult = prod(sector_multiplicity.(it))
   return build_trees(block_sectors, isdual.(it)) .=> block_mult
 end
 
 function compute_inner_ranges(fusion_trees_mult)
-  fused_leg = blockmergesort(
+  fused_leg = sectormergesort(
     gradedrange(root_sector.(first.(fusion_trees_mult)) .=> last.(fusion_trees_mult))
   )
-  range_mapping = Dict{fieldtype(eltype(fusion_trees_mult), 1),typeof(Block(1)[1:1])}()
-  fused_sectors = blocklabels(fused_leg)
+  range_mapping = Dict{type_parameters(eltype(fusion_trees_mult), 1),typeof(Block(1)[1:1])}()
+  fused_sectors = sectors(fused_leg)
   shifts = ones(Int, blocklength(fused_leg))
   for (f, m) in fusion_trees_mult
     s = root_sector(f)
@@ -225,24 +236,25 @@ end
 
 function initialize_data_matrix(
   elt::Type{<:Number},
-  mat_row_axis::AbstractGradedUnitRange,
-  mat_col_axis::AbstractGradedUnitRange,
+  codomain_axis::AbstractGradedUnitRange,
+  domain_axis::AbstractGradedUnitRange,
 )
+  @assert sector_type(codomain_axis) == sector_type(domain_axis)
   # non-abelian fusion trees have float eltype: need compatible type
-  promoted = promote_type(elt, fusiontree_eltype(sector_type(mat_row_axis)))
-  mat = BlockSparseArray{promoted}(undef, mat_row_axis, mat_col_axis)
-  initialize_allowed_sectors!(mat)
-  return mat
-end
-
-function initialize_allowed_sectors!(mat::AbstractMatrix)
-  row_sectors = blocklabels(axes(mat, 1))
-  col_sectors = blocklabels(dual(axes(mat, 2)))
+  promoted = promote_type(elt, fusiontree_eltype(sector_type(domain_axis)))
+  mat = BlockSparseArray{promoted}(
+    undef,
+    blockedrange(sector_multiplicities(codomain_axis)),
+    blockedrange(sector_multiplicities(domain_axis)),
+  )
+  row_sectors = sectors(codomain_axis)
+  col_sectors = sectors(domain_axis)
   row_block_indices = findall(in(col_sectors), row_sectors)
   col_block_indices = findall(in(row_sectors), col_sectors)
   for rc in zip(row_block_indices, col_block_indices)
     mat[Block(rc)] = mat[Block(rc)]
   end
+  return mat
 end
 
 checkaxes_dual(axes1, axes2) = checkaxes(axes1, dual.(axes2))
diff --git a/src/fusiontensor/linear_algebra_interface.jl b/src/fusiontensor/linear_algebra_interface.jl
index 9b0a57e..4a631e5 100644
--- a/src/fusiontensor/linear_algebra_interface.jl
+++ b/src/fusiontensor/linear_algebra_interface.jl
@@ -5,8 +5,7 @@ using LinearAlgebra: LinearAlgebra, mul!, norm, tr
 using BlockArrays: Block, blocks
 
 using BlockSparseArrays: eachblockstoredindex
-using GradedArrays: blocklabels
-using GradedArrays.SymmetrySectors: quantum_dimension
+using GradedArrays: quantum_dimension, sectors
 
 # allow to contract with different eltype and let BlockSparseArray ensure compatibility
 # impose matching type and number of axes at compile time
@@ -37,7 +36,7 @@ end
 
 function LinearAlgebra.norm(ft::FusionTensor)
   m = data_matrix(ft)
-  row_sectors = blocklabels(matrix_row_axis(ft))
+  row_sectors = sectors(codomain_axis(ft))
   n2 = sum(eachblockstoredindex(m); init=zero(real(eltype(ft)))) do b
     return quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * norm(m[b])^2
   end
@@ -46,7 +45,7 @@ end
 
 function LinearAlgebra.tr(ft::FusionTensor)
   m = data_matrix(ft)
-  row_sectors = blocklabels(matrix_row_axis(ft))
+  row_sectors = sectors(codomain_axis(ft))
   return sum(eachblockstoredindex(m); init=zero(eltype(ft))) do b
     return quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * tr(m[b])
   end
diff --git a/src/permutedims/permutedims.jl b/src/permutedims/permutedims.jl
index da97e6a..252eb5c 100644
--- a/src/permutedims/permutedims.jl
+++ b/src/permutedims/permutedims.jl
@@ -7,12 +7,11 @@ using TensorAlgebra: BlockedPermutation, permmortar, blockpermute
 
 function naive_permutedims(ft, biperm::BlockedPermutation{2})
   @assert ndims(ft) == length(biperm)
-  new_codomain_legs, new_domain_legs = blockpermute(axes(ft), biperm)
 
   # naive permute: cast to dense, permutedims, cast to FusionTensor
   arr = Array(ft)
   permuted_arr = permutedims(arr, Tuple(biperm))
-  permuted = to_fusiontensor(permuted_arr, new_codomain_legs, new_domain_legs)
+  permuted = to_fusiontensor(permuted_arr, blocks(axes(ft)[biperm])...)
   return permuted
 end
 
@@ -30,7 +29,7 @@ function fusiontensor_permutedims(
 end
 
 function fusiontensor_permutedims(ft, biperm::BlockedPermutation{2})
-  @assert ndims(ft) == length(biperm)
+  ndims(ft) == length(biperm) || throw(ArgumentError("Invalid permutation length"))
 
   # early return for identity operation. Do not copy. Also handle tricky 0-dim case.
   if ndims_codomain(ft) == first(blocklengths(biperm))  # compile time
@@ -39,8 +38,7 @@ function fusiontensor_permutedims(ft, biperm::BlockedPermutation{2})
     end
   end
 
-  new_codomain_legs, new_domain_legs = blockpermute(axes(ft), biperm)
-  new_ft = FusionTensor(eltype(ft), new_codomain_legs, new_domain_legs)
+  new_ft = FusionTensor(eltype(ft), axes(ft)[biperm])
   fusiontensor_permutedims!(new_ft, ft, Tuple(biperm))
   return new_ft
 end
diff --git a/src/permutedims/unitaries.jl b/src/permutedims/unitaries.jl
index ad09d41..6fbc6c7 100644
--- a/src/permutedims/unitaries.jl
+++ b/src/permutedims/unitaries.jl
@@ -3,7 +3,7 @@
 using BlockArrays: Block, findblock
 using LRUCache: LRU
 
-using GradedArrays.SymmetrySectors: quantum_dimension
+using GradedArrays: quantum_dimension
 
 const unitary_cache = LRU{Any,Float64}(; maxsize=10000)  # TBD size
 
diff --git a/test/Project.toml b/test/Project.toml
index f1f3e8d..222a336 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -13,13 +13,13 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Aqua = "0.8.11"
-BlockArrays = "1.5.0"
-BlockSparseArrays = "0.4.0"
-FusionTensors = "0.3.0"
-GradedArrays = "0.2"
+BlockArrays = "1.6"
+BlockSparseArrays = "0.7"
+FusionTensors = "0.4"
+GradedArrays = "0.4"
 LinearAlgebra = "1.10.0"
 SafeTestsets = "0.1.0"
 Suppressor = "0.2.8"
-TensorAlgebra = "0.2.7"
-TensorProducts = "0.1.3"
+TensorAlgebra = "0.3"
+TensorProducts = "0.1"
 Test = "1.10.0"
diff --git a/test/setup.jl b/test/setup.jl
index db8cf00..0607a81 100644
--- a/test/setup.jl
+++ b/test/setup.jl
@@ -1,5 +1,5 @@
 using LinearAlgebra: Adjoint
-
+using BlockArrays: blocklengths
 using BlockSparseArrays: BlockSparseMatrix, eachblockstoredindex
 using FusionTensors:
   FusionTensor,
@@ -8,23 +8,11 @@ using FusionTensors:
   domain_axes,
   checkaxes,
   checkaxes_dual,
-  matrix_column_axis,
-  matrix_row_axis,
+  domain_axis,
+  codomain_axis,
   ndims_codomain,
   ndims_domain
-using GradedArrays: blocklabels, dual, space_isequal
-
-function check_data_matrix_axes(
-  mat::BlockSparseMatrix, codomain_legs::Tuple, domain_legs::Tuple
-)
-  ft0 = FusionTensor(Float64, codomain_legs, domain_legs)
-  @assert space_isequal(matrix_row_axis(ft0), axes(mat, 1))
-  @assert space_isequal(matrix_column_axis(ft0), axes(mat, 2))
-end
-
-function check_data_matrix_axes(mat::Adjoint, codomain_legs::Tuple, domain_legs::Tuple)
-  return check_data_matrix_axes(adjoint(mat), dual.(domain_legs), dual.(codomain_legs))
-end
+using GradedArrays: dual, sectors, sector_multiplicities, space_isequal
 
 function check_sanity(ft::FusionTensor)
   nca = ndims_domain(ft)
@@ -42,15 +30,14 @@ function check_sanity(ft::FusionTensor)
 
   m = data_matrix(ft)
   @assert ndims(m) == 2 "invalid data_matrix ndims"
-  row_axis = matrix_row_axis(ft)
-  column_axis = matrix_column_axis(ft)
-  @assert row_axis === axes(m, 1) "invalid row_axis"
-  @assert column_axis === axes(m, 2) "invalid column_axis"
-  check_data_matrix_axes(m, codomain_axes(ft), domain_axes(ft))
+  row_axis = codomain_axis(ft)
+  column_axis = domain_axis(ft)
+  @assert sector_multiplicities(row_axis) == blocklengths(axes(m, 1)) "invalid row_axis"
+  @assert sector_multiplicities(column_axis) == blocklengths(axes(m, 2)) "invalid column_axis"
 
   for b in eachblockstoredindex(m)
     ir, ic = Int.(Tuple(b))
-    @assert blocklabels(row_axis)[ir] == blocklabels(dual(column_axis))[ic] "forbidden block"
+    @assert sectors(row_axis)[ir] == sectors(dual(column_axis))[ic] "forbidden block"
   end
   return nothing
 end
diff --git a/test/test_array_cast.jl b/test/test_array_cast.jl
index 1d7f87f..049b0d7 100644
--- a/test/test_array_cast.jl
+++ b/test/test_array_cast.jl
@@ -4,8 +4,7 @@ using Test: @test, @test_throws, @testset
 using BlockArrays: Block, BlockedArray, blocksize
 
 using FusionTensors: FusionTensor, data_matrix, to_fusiontensor
-using GradedArrays: dual, gradedrange
-using GradedArrays.SymmetrySectors: O2, SectorProduct, SU2, TrivialSector, U1
+using GradedArrays: O2, SectorProduct, SU2, TrivialSector, U1, dual, gradedrange
 using TensorProducts: tensor_product
 
 include("setup.jl")
diff --git a/test/test_basics.jl b/test/test_basics.jl
index dbfc91e..e150d97 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -9,13 +9,21 @@ using FusionTensors:
   FusionTensor,
   checkaxes,
   checkaxes_dual,
-  matrix_column_axis,
-  matrix_row_axis,
-  matrix_size,
+  codomain_axis,
+  domain_axis,
   ndims_domain,
   ndims_codomain
-using GradedArrays: blockmergesort, dual, flip, gradedrange, sector_type, space_isequal
-using GradedArrays.SymmetrySectors: U1, SU2, SectorProduct, TrivialSector, Z
+using GradedArrays:
+  U1,
+  SU2,
+  SectorProduct,
+  TrivialSector,
+  Z,
+  dual,
+  flip,
+  gradedrange,
+  sector_type,
+  space_isequal
 using TensorAlgebra: tuplemortar
 using TensorProducts: tensor_product
 
@@ -28,8 +36,8 @@ include("setup.jl")
   # check dual convention when initializing data_matrix
   ft0 = FusionTensor(Float64, (g1,), (g2,))
   @test ft0 isa FusionTensor
-  @test space_isequal(matrix_row_axis(ft0), g1)
-  @test space_isequal(matrix_column_axis(ft0), g2)
+  @test space_isequal(codomain_axis(ft0), g1)
+  @test space_isequal(domain_axis(ft0), g2)
 
   m = BlockSparseArray{Float64}(undef, g1, g2)
   ft1 = FusionTensor(m, (g1,), (g2,))
@@ -43,9 +51,9 @@ include("setup.jl")
   @test checkaxes(domain_axes(ft1), (g2,))
   @test ndims_codomain(ft1) == 1
   @test ndims_domain(ft1) == 1
-  @test matrix_size(ft1) == (6, 5)
-  @test space_isequal(matrix_row_axis(ft1), g1)
-  @test space_isequal(matrix_column_axis(ft1), g2)
+  @test size(data_matrix(ft1)) == (6, 5)
+  @test space_isequal(codomain_axis(ft1), g1)
+  @test space_isequal(domain_axis(ft1), g2)
   @test isnothing(check_sanity(ft0))
   @test isnothing(check_sanity(ft1))
   @test sector_type(ft1) === U1{Int}
@@ -117,9 +125,9 @@ end
   @test axes(ft) == tuplemortar(((g1, g2), (g3, g4)))
   @test ndims_codomain(ft) == 2
   @test ndims_domain(ft) == 2
-  @test matrix_size(ft) == (30, 12)
-  @test space_isequal(matrix_row_axis(ft), gr)
-  @test space_isequal(matrix_column_axis(ft), gc)
+  @test size(data_matrix(ft)) == (30, 12)
+  @test space_isequal(codomain_axis(ft), gr)
+  @test space_isequal(domain_axis(ft), gc)
   @test isnothing(check_sanity(ft))
 
   @test ndims(ft) == 4
@@ -179,46 +187,46 @@ end
   ft4 = ft3 + ft3
   @test codomain_axes(ft4) === codomain_axes(ft3)
   @test domain_axes(ft4) === domain_axes(ft3)
-  @test space_isequal(matrix_row_axis(ft4), matrix_row_axis(ft3))
-  @test space_isequal(matrix_column_axis(ft4), matrix_column_axis(ft3))
+  @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
+  @test space_isequal(domain_axis(ft4), domain_axis(ft3))
   @test isnothing(check_sanity(ft4))
 
   ft4 = ft3 - ft3
   @test codomain_axes(ft4) === codomain_axes(ft3)
   @test domain_axes(ft4) === domain_axes(ft3)
-  @test space_isequal(matrix_row_axis(ft4), matrix_row_axis(ft3))
-  @test space_isequal(matrix_column_axis(ft4), matrix_column_axis(ft3))
+  @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
+  @test space_isequal(domain_axis(ft4), domain_axis(ft3))
   @test isnothing(check_sanity(ft4))
 
   ft4 = 2 * ft3
   @test codomain_axes(ft4) === codomain_axes(ft3)
   @test domain_axes(ft4) === domain_axes(ft3)
-  @test space_isequal(matrix_row_axis(ft4), matrix_row_axis(ft3))
-  @test space_isequal(matrix_column_axis(ft4), matrix_column_axis(ft3))
+  @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
+  @test space_isequal(domain_axis(ft4), domain_axis(ft3))
   @test isnothing(check_sanity(ft4))
   @test eltype(ft4) == Float64
 
   ft4 = 2.0 * ft3
   @test codomain_axes(ft4) === codomain_axes(ft3)
   @test domain_axes(ft4) === domain_axes(ft3)
-  @test space_isequal(matrix_row_axis(ft4), matrix_row_axis(ft3))
-  @test space_isequal(matrix_column_axis(ft4), matrix_column_axis(ft3))
+  @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
+  @test space_isequal(domain_axis(ft4), domain_axis(ft3))
   @test isnothing(check_sanity(ft4))
   @test eltype(ft4) == Float64
 
   ft4 = ft3 / 2.0
   @test codomain_axes(ft4) === codomain_axes(ft3)
   @test domain_axes(ft4) === domain_axes(ft3)
-  @test space_isequal(matrix_row_axis(ft4), matrix_row_axis(ft3))
-  @test space_isequal(matrix_column_axis(ft4), matrix_column_axis(ft3))
+  @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
+  @test space_isequal(domain_axis(ft4), domain_axis(ft3))
   @test isnothing(check_sanity(ft4))
   @test eltype(ft4) == Float64
 
   ft5 = 2.0im * ft3
   @test codomain_axes(ft5) === codomain_axes(ft3)
   @test domain_axes(ft5) === domain_axes(ft3)
-  @test space_isequal(matrix_row_axis(ft5), matrix_row_axis(ft3))
-  @test space_isequal(matrix_column_axis(ft5), matrix_column_axis(ft3))
+  @test space_isequal(codomain_axis(ft5), codomain_axis(ft3))
+  @test space_isequal(domain_axis(ft5), domain_axis(ft3))
   @test isnothing(check_sanity(ft4))
   @test eltype(ft5) == ComplexF64
 
@@ -230,8 +238,8 @@ end
   @test isnothing(check_sanity(ft6))
   @test codomain_axes(ft6) === codomain_axes(ft5)
   @test domain_axes(ft6) === domain_axes(ft5)
-  @test space_isequal(matrix_row_axis(ft6), matrix_row_axis(ft5))
-  @test space_isequal(matrix_column_axis(ft6), matrix_column_axis(ft5))
+  @test space_isequal(codomain_axis(ft6), codomain_axis(ft5))
+  @test space_isequal(domain_axis(ft6), domain_axis(ft5))
   @test eltype(ft6) == ComplexF64
 
   ad = adjoint(ft3)
@@ -259,11 +267,11 @@ end
   ft = FusionTensor(Float64, (g1,), (dual(g2), dual(g3)))
   @test sector_type(ft) === S
   gr = gradedrange([SectorProduct(U1(1), SU2(0), Z{2}(0)) => 1])
-  @test space_isequal(matrix_row_axis(ft), gr)
+  @test space_isequal(codomain_axis(ft), gr)
   gc = gradedrange([
     SectorProduct(U1(2), SU2(0), Z{2}(1)) => 1, SectorProduct(U1(2), SU2(1), Z{2}(1)) => 1
   ])
-  @test space_isequal(matrix_column_axis(ft), dual(gc))
+  @test space_isequal(domain_axis(ft), dual(gc))
 
   gA = gradedrange([SectorProduct(; A=U1(1)) => 1])
   gB = gradedrange([SectorProduct(; B=SU2(1//2)) => 1])
@@ -273,6 +281,6 @@ end
 
   ft = FusionTensor(Float64, (gA, gB), (dual(gA), dual(gB), gC))
   @test sector_type(ft) === S
-  @test space_isequal(matrix_row_axis(ft), gABC)
-  @test space_isequal(matrix_column_axis(ft), dual(gABC))
+  @test space_isequal(codomain_axis(ft), gABC)
+  @test space_isequal(domain_axis(ft), dual(gABC))
 end
diff --git a/test/test_contraction.jl b/test/test_contraction.jl
index da60b66..25a0041 100644
--- a/test/test_contraction.jl
+++ b/test/test_contraction.jl
@@ -3,9 +3,8 @@ using Test: @test, @testset, @test_broken
 
 using BlockSparseArrays: BlockSparseArray
 using FusionTensors: FusionTensor, domain_axes, codomain_axes
-using GradedArrays: dual, gradedrange
-using GradedArrays.SymmetrySectors: U1
-using TensorAlgebra: contract
+using GradedArrays: U1, dual, gradedrange
+using TensorAlgebra: contract, tuplemortar
 
 include("setup.jl")
 
@@ -49,18 +48,22 @@ end
   ft3 = FusionTensor(Float64, dual.((g3, g4)), dual.((g1, g2)))
 
   ft4, legs = contract(ft1, (1, 2, 3, 4), ft2, (3, 4, 5))
-  @test legs == (1, 2, 5)
+  @test legs == tuplemortar(((1, 2), (5,)))
   @test isnothing(check_sanity(ft4))
   @test domain_axes(ft4) === domain_axes(ft2)
   @test codomain_axes(ft4) === codomain_axes(ft1)
 
   ft5 = contract((1, 2, 5), ft1, (1, 2, 3, 4), ft2, (3, 4, 5))
   @test isnothing(check_sanity(ft5))
+  @test ft4 ≈ ft5
 
-  # biperm is not allowed
-  @test_broken contract(((1, 2), (5,)), ft1, (1, 2, 3, 4), ft2, (3, 4, 5))
+  ft6 = contract(tuplemortar(((1, 2), (5,))), ft1, (1, 2, 3, 4), ft2, (3, 4, 5))
+  @test isnothing(check_sanity(ft6))
+  @test ft4 ≈ ft6
 
   @test permutedims(ft1, (), (1, 2, 3, 4)) * permutedims(ft3, (3, 4, 1, 2), ()) isa
     FusionTensor{Float64,0}
-  @test_broken contract(ft1, (1, 2, 3, 4), ft3, (3, 4, 1, 2))
+  ft7, legs = contract(ft1, (1, 2, 3, 4), ft3, (3, 4, 1, 2))
+  @test legs == tuplemortar(((), ()))
+  @test ft7 isa FusionTensor{Float64,0}
 end
diff --git a/test/test_fusion_trees.jl b/test/test_fusion_trees.jl
index 4f7de0c..565997b 100644
--- a/test/test_fusion_trees.jl
+++ b/test/test_fusion_trees.jl
@@ -12,8 +12,8 @@ using FusionTensors:
   leaves,
   outer_multiplicity_indices,
   root_sector
-using GradedArrays: blocklabels, sector_type
-using GradedArrays.SymmetrySectors: ×, SectorProduct, SU, SU2, TrivialSector, arguments
+using GradedArrays:
+  ×, SectorProduct, SU, SU2, TrivialSector, arguments, dual, flip, gradedrange, sector_type
 
 @testset "Trivial fusion trees" begin
   q = TrivialSector()
@@ -95,6 +95,41 @@ end
 end
 
 @testset "SU(3) SectorFusionTree" begin
+  # convention: irreps are already dualed if needed, arrows do not affect them. They only
+  # affect the basis on which the tree projects for self-dual irreps.
+  f3 = SU{3}((1, 0))
+  c3 = dual(f3)
+
+  trees = build_trees((f3, f3), (false, false))
+  @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))]
+
+  trees = build_trees((f3, f3), (true, false))
+  @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))]
+
+  trees = build_trees((f3, f3), (false, true))
+  @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))]
+
+  trees = build_trees((f3, f3), (true, true))
+  @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))]
+
+  trees = build_trees((f3, c3), (false, false))
+  @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))]
+
+  trees = build_trees((c3, c3), (false, false))
+  @test root_sector.(trees) == [SU{3}((1, 0)), SU{3}((2, 2))]
+
+  # test GradedUnitRange interface
+  g = gradedrange([SU{3}((1, 0)) => 1])
+  trees = build_trees(g, g)
+  @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))]
+  trees = build_trees(g, flip(g))
+  @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))]
+  trees = build_trees(g, flip(dual(g)))
+  @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))]
+  trees = build_trees(g, dual(g))
+  @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))]
+
+  # test outer outer_multiplicity > 1
   a8 = SU{3}((2, 1))
   trees = build_trees((a8, a8), (false, false))
   @test length(trees) == 6
diff --git a/test/test_linear_algebra.jl b/test/test_linear_algebra.jl
index 94218f5..9408e4c 100644
--- a/test/test_linear_algebra.jl
+++ b/test/test_linear_algebra.jl
@@ -5,8 +5,7 @@ using BlockArrays: BlockArrays
 
 using BlockSparseArrays: BlockSparseArrays
 using FusionTensors: FusionTensor, to_fusiontensor
-using GradedArrays: dual, gradedrange
-using GradedArrays.SymmetrySectors: U1, SU2, TrivialSector
+using GradedArrays: SU2, TrivialSector, U1, dual, gradedrange
 
 include("setup.jl")
 
diff --git a/test/test_permutedims.jl b/test/test_permutedims.jl
index 226523e..1d7c9a3 100644
--- a/test/test_permutedims.jl
+++ b/test/test_permutedims.jl
@@ -1,17 +1,16 @@
-using Test: @test, @testset, @test_broken
+using Test: @test, @testset, @test_broken, @test_throws
 
 using FusionTensors:
   FusionTensor,
   data_matrix,
   checkaxes,
-  matrix_column_axis,
-  matrix_row_axis,
+  codomain_axis,
+  domain_axis,
   naive_permutedims,
   ndims_domain,
   ndims_codomain,
   to_fusiontensor
-using GradedArrays: dual, gradedrange, space_isequal
-using GradedArrays.SymmetrySectors: O2, U1, SectorProduct, SU2
+using GradedArrays: O2, U1, SectorProduct, SU2, dual, gradedrange, space_isequal
 using TensorAlgebra: permmortar, tuplemortar
 
 include("setup.jl")
@@ -46,9 +45,12 @@ include("setup.jl")
 
       ft4 = permutedims(ft3, (2, 3), (4, 1))
       @test checkaxes(axes(ft1), axes(ft4))
-      @test space_isequal(matrix_column_axis(ft1), matrix_column_axis(ft4))
-      @test space_isequal(matrix_row_axis(ft1), matrix_row_axis(ft4))
+      @test space_isequal(codomain_axis(ft1), codomain_axis(ft4))
+      @test space_isequal(domain_axis(ft1), domain_axis(ft4))
       @test ft4 ≈ ft1
+
+      @test_throws MethodError permutedims(ft1, (2, 3, 4, 1))
+      @test_throws ArgumentError permutedims(ft1, (2, 3), (5, 4, 1))
     end
   end