diff --git a/Project.toml b/Project.toml
index 127bb445..ce3720bb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "QuantumOpticsBase"
 uuid = "4f57444f-1401-5e15-980d-4471b28d5678"
-version = "0.4.1"
+version = "0.4.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -19,7 +19,7 @@ Adapt = "1, 2, 3.3"
 FFTW = "1.2"
 FillArrays = "0.13, 1"
 LRUCache = "1"
-QuantumInterface = "0.1.0"
+QuantumInterface = "0.2.0"
 Strided = "1, 2"
 UnsafeArrays = "1"
 julia = "1.3"
diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl
index 6e06a898..a75f8cf7 100644
--- a/src/QuantumOpticsBase.jl
+++ b/src/QuantumOpticsBase.jl
@@ -3,8 +3,11 @@ module QuantumOpticsBase
 using SparseArrays, LinearAlgebra, LRUCache, Strided, UnsafeArrays, FillArrays
 import LinearAlgebra: mul!, rmul!
 
-import QuantumInterface: dagger, directsum, ⊕, dm, embed, expect, permutesystems,
-        projector, ptrace, reduced, tensor, ⊗
+import QuantumInterface: dagger, directsum, ⊕, dm, embed, expect, identityoperator,
+        permutesystems, projector, ptrace, reduced, tensor, ⊗, variance
+
+# index helpers
+import QuantumInterface: complement, remove, shiftremove, reducedindices!, check_indices, check_sortedindices, check_embed_indices
 
 export Basis, GenericBasis, CompositeBasis, basis,
         tensor, ⊗, permutesystems, @samebases,
@@ -55,7 +58,6 @@ export Basis, GenericBasis, CompositeBasis, basis,
         SumBasis, directsum, ⊕, LazyDirectSum, getblock, setblock!,
         qfunc, wigner, coherentspinstate, qfuncsu2, wignersu2
 
-include("sortedindices.jl")
 include("polynomials.jl")
 include("bases.jl")
 include("states.jl")
diff --git a/src/operators.jl b/src/operators.jl
index e3c54a5d..8f55932e 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -13,94 +13,7 @@ abstract type DataOperator{BL,BR} <: AbstractOperator{BL,BR} end
 
 
 # Common error messages
-arithmetic_unary_error(funcname, x::AbstractOperator) = throw(ArgumentError("$funcname is not defined for this type of operator: $(typeof(x)).\nTry to convert to another operator type first with e.g. dense() or sparse()."))
-arithmetic_binary_error(funcname, a::AbstractOperator, b::AbstractOperator) = throw(ArgumentError("$funcname is not defined for this combination of types of operators: $(typeof(a)), $(typeof(b)).\nTry to convert to a common operator type first with e.g. dense() or sparse()."))
-addnumbererror() = throw(ArgumentError("Can't add or subtract a number and an operator. You probably want 'op + identityoperator(op)*x'."))
-
-length(a::AbstractOperator) = length(a.basis_l)::Int*length(a.basis_r)::Int
-basis(a::AbstractOperator) = (check_samebases(a); a.basis_l)
-
-# Ensure scalar broadcasting
-Base.broadcastable(x::AbstractOperator) = Ref(x)
-
-# Arithmetic operations
-+(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Addition", a, b)
-+(a::Number, b::AbstractOperator) = addnumbererror()
-+(a::AbstractOperator, b::Number) = addnumbererror()
-+(a::AbstractOperator) = a
-
--(a::AbstractOperator) = arithmetic_unary_error("Negation", a)
--(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Subtraction", a, b)
--(a::Number, b::AbstractOperator) = addnumbererror()
--(a::AbstractOperator, b::Number) = addnumbererror()
-
-*(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Multiplication", a, b)
-^(a::AbstractOperator, b::Integer) = Base.power_by_squaring(a, b)
-
-
-dagger(a::AbstractOperator) = arithmetic_unary_error("Hermitian conjugate", a)
-Base.adjoint(a::AbstractOperator) = dagger(a)
-
-conj(a::AbstractOperator) = arithmetic_unary_error("Complex conjugate", a)
-conj!(a::AbstractOperator) = conj(a::AbstractOperator)
-
-# dense(a::AbstractOperator) = arithmetic_unary_error("Conversion to dense", a)
-
-transpose(a::AbstractOperator) = arithmetic_unary_error("Transpose", a)
-
-"""
-    ishermitian(op::AbstractOperator)
-
-Check if an operator is Hermitian.
-"""
-ishermitian(op::AbstractOperator) = arithmetic_unary_error(ishermitian, op)
-
-
-"""
-    tensor(x::AbstractOperator, y::AbstractOperator, z::AbstractOperator...)
-
-Tensor product ``\\hat{x}⊗\\hat{y}⊗\\hat{z}⊗…`` of the given operators.
-"""
-tensor(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Tensor product", a, b)
-tensor(op::AbstractOperator) = op
-tensor(operators::AbstractOperator...) = reduce(tensor, operators)
-
-
-"""
-    embed(basis1[, basis2], indices::Vector, operators::Vector)
-
-Tensor product of operators where missing indices are filled up with identity operators.
-"""
-function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
-               indices, operators::Vector{T}) where T<:AbstractOperator
-
-    @assert check_embed_indices(indices)
-
-    N = length(basis_l.bases)
-    @assert length(basis_r.bases) == N
-    @assert length(indices) == length(operators)
-
-    # Embed all single-subspace operators.
-    idxop_sb = [x for x in zip(indices, operators) if x[1] isa Integer]
-    indices_sb = [x[1] for x in idxop_sb]
-    ops_sb = [x[2] for x in idxop_sb]
-
-    for (idxsb, opsb) in zip(indices_sb, ops_sb)
-        (opsb.basis_l == basis_l.bases[idxsb]) || throw(IncompatibleBases())
-        (opsb.basis_r == basis_r.bases[idxsb]) || throw(IncompatibleBases())
-    end
-
-    S = length(operators) > 0 ? mapreduce(eltype, promote_type, operators) : Any
-    embed_op = tensor([i ∈ indices_sb ? ops_sb[indexin(i, indices_sb)[1]] : identityoperator(T, S, basis_l.bases[i], basis_r.bases[i]) for i=1:N]...)
-
-    # Embed all joint-subspace operators.
-    idxop_comp = [x for x in zip(indices, operators) if x[1] isa Array]
-    for (idxs, op) in idxop_comp
-        embed_op *= embed(basis_l, basis_r, idxs, op)
-    end
-
-    return embed_op
-end
+using QuantumInterface: arithmetic_binary_error, arithmetic_unary_error, addnumbererror
 
 
 """
@@ -156,12 +69,6 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
 
     return unpermuted_op
 end
-# The dictionary implementation works for non-DataOperators
-embed(basis_l::CompositeBasis, basis_r::CompositeBasis, indices, op::T) where T<:AbstractOperator = embed(basis_l, basis_r, Dict(indices=>op))
-
-embed(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Integer, op::AbstractOperator) = embed(basis_l, basis_r, index, [op])
-embed(basis::CompositeBasis, indices, operators::Vector{T}) where {T<:AbstractOperator} = embed(basis, basis, indices, operators)
-embed(basis::CompositeBasis, indices, op::AbstractOperator) = embed(basis, basis, indices, op)
 
 function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
                 index::Integer, op::T) where T<:DataOperator
@@ -195,70 +102,6 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
     return Operator(basis_l, basis_r, data)
 end
 
-"""
-    embed(basis1[, basis2], operators::Dict)
-
-`operators` is a dictionary `Dict{Vector{Int}, AbstractOperator}`. The integer vector
-specifies in which subsystems the corresponding operator is defined.
-"""
-function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
-               operators::Dict{<:Vector{<:Integer}, T}) where T<:AbstractOperator
-    @assert length(basis_l.bases) == length(basis_r.bases)
-    N = length(basis_l.bases)::Int # type assertion to help type inference
-    if length(operators) == 0
-        return identityoperator(T, basis_l, basis_r)
-    end
-    indices, operator_list = zip(operators...)
-    operator_list = [operator_list...;]
-    S = mapreduce(eltype, promote_type, operator_list)
-    indices_flat = [indices...;]::Vector{Int} # type assertion to help type inference
-    start_indices_flat = [i[1] for i in indices]
-    complement_indices_flat = Int[i for i=1:N if i ∉ indices_flat]
-    operators_flat = AbstractOperator[]
-    if all(([minimum(I):maximum(I);]==I)::Bool for I in indices) # type assertion to help type inference
-        for i in 1:N
-            if i in complement_indices_flat
-                push!(operators_flat, identityoperator(T, S, basis_l.bases[i], basis_r.bases[i]))
-            elseif i in start_indices_flat
-                push!(operators_flat, operator_list[indexin(i, start_indices_flat)[1]])
-            end
-        end
-        return tensor(operators_flat...)
-    else
-        complement_operators = [identityoperator(T, S, basis_l.bases[i], basis_r.bases[i]) for i in complement_indices_flat]
-        op = tensor([operator_list; complement_operators]...)
-        perm = sortperm([indices_flat; complement_indices_flat])
-        return permutesystems(op, perm)
-    end
-end
-embed(basis_l::CompositeBasis, basis_r::CompositeBasis, operators::Dict{<:Integer, T}; kwargs...) where {T<:AbstractOperator} = embed(basis_l, basis_r, Dict([i]=>op_i for (i, op_i) in operators); kwargs...)
-embed(basis::CompositeBasis, operators::Dict{<:Integer, T}; kwargs...) where {T<:AbstractOperator} = embed(basis, basis, operators; kwargs...)
-embed(basis::CompositeBasis, operators::Dict{<:Vector{<:Integer}, T}; kwargs...) where {T<:AbstractOperator} = embed(basis, basis, operators; kwargs...)
-
-
-"""
-    tr(x::AbstractOperator)
-
-Trace of the given operator.
-"""
-tr(x::AbstractOperator) = arithmetic_unary_error("Trace", x)
-
-ptrace(a::AbstractOperator, index) = arithmetic_unary_error("Partial trace", a)
-
-"""
-    normalize(op)
-
-Return the normalized operator so that its `tr(op)` is one.
-"""
-normalize(op::AbstractOperator) = op/tr(op)
-
-"""
-    normalize!(op)
-
-In-place normalization of the given operator so that its `tr(x)` is one.
-"""
-normalize!(op::AbstractOperator) = throw(ArgumentError("normalize! is not defined for this type of operator: $(typeof(op)).\n You may have to fall back to the non-inplace version 'normalize()'."))
-
 """
     expect(op, state)
 
@@ -267,30 +110,15 @@ Expectation value of the given operator `op` for the specified `state`.
 `state` can either be a (density) operator or a ket.
 """
 expect(op::AbstractOperator{B,B}, state::Ket{B}) where B = dot(state.data, (op * state).data)
-expect(op::AbstractOperator{B1,B2}, state::AbstractOperator{B2,B2}) where {B1,B2} = tr(op*state)
 
-"""
-    expect(index, op, state)
-
-If an `index` is given, it assumes that `op` is defined in the subsystem specified by this number.
-"""
-function expect(indices, op::AbstractOperator{B1,B2}, state::AbstractOperator{B3,B3}) where {B1,B2,B3<:CompositeBasis}
-    N = length(state.basis_l.shape)
-    indices_ = complement(N, indices)
-    expect(op, ptrace(state, indices_))
-end
 function expect(indices, op::AbstractOperator{B,B}, state::Ket{B2}) where {B,B2<:CompositeBasis}
     N = length(state.basis.shape)
     indices_ = complement(N, indices)
     expect(op, ptrace(state, indices_))
 end
 
-expect(index::Integer, op::AbstractOperator{B1,B2}, state::AbstractOperator{B3,B3}) where {B1,B2,B3<:CompositeBasis} = expect([index], op, state)
 expect(index::Integer, op::AbstractOperator{B,B}, state::Ket{B2}) where {B,B2<:CompositeBasis} = expect([index], op, state)
 
-expect(op::AbstractOperator, states::Vector) = [expect(op, state) for state=states]
-expect(indices, op::AbstractOperator, states::Vector) = [expect(indices, op, state) for state=states]
-
 """
     variance(op, state)
 
@@ -302,60 +130,14 @@ function variance(op::AbstractOperator{B,B}, state::Ket{B}) where B
     x = op*state
     state.data'*(op*x).data - (state.data'*x.data)^2
 end
-function variance(op::AbstractOperator{B,B}, state::AbstractOperator{B,B}) where B
-    expect(op*op, state) - expect(op, state)^2
-end
-
-"""
-    variance(index, op, state)
 
-If an `index` is given, it assumes that `op` is defined in the subsystem specified by this number
-"""
-function variance(indices, op::AbstractOperator{B,B}, state::AbstractOperator{BC,BC}) where {B,BC<:CompositeBasis}
-    N = length(state.basis_l.shape)
-    indices_ = complement(N, indices)
-    variance(op, ptrace(state, indices_))
-end
 function variance(indices, op::AbstractOperator{B,B}, state::Ket{BC}) where {B,BC<:CompositeBasis}
     N = length(state.basis.shape)
     indices_ = complement(N, indices)
     variance(op, ptrace(state, indices_))
 end
 
-variance(index::Integer, op::AbstractOperator{B,B}, state::AbstractOperator{BC,BC}) where {B,BC<:CompositeBasis} = variance([index], op, state)
 variance(index::Integer, op::AbstractOperator{B,B}, state::Ket{BC}) where {B,BC<:CompositeBasis} = variance([index], op, state)
-variance(op::AbstractOperator, states::Vector) = [variance(op, state) for state=states]
-variance(indices, op::AbstractOperator, states::Vector) = [variance(indices, op, state) for state=states]
-
-
-"""
-    exp(op::AbstractOperator)
-
-Operator exponential.
-"""
-exp(op::AbstractOperator) = throw(ArgumentError("exp() is not defined for this type of operator: $(typeof(op)).\nTry to convert to dense operator first with dense()."))
-
-permutesystems(a::AbstractOperator, perm) = arithmetic_unary_error("Permutations of subsystems", a)
-
-"""
-    identityoperator(a::Basis[, b::Basis])
-    identityoperator(::Type{<:AbstractOperator}, a::Basis[, b::Basis])
-    identityoperator(::Type{<:Number}, a::Basis[, b::Basis])
-    identityoperator(::Type{<:AbstractOperator}, ::Type{<:Number}, a::Basis[, b::Basis])
-
-Return an identityoperator in the given bases. One can optionally specify the container
-type which has to a subtype of [`AbstractOperator`](@ref) as well as the number type
-to be used in the identity matrix.
-"""
-identityoperator(::Type{T}, ::Type{S}, b1::Basis, b2::Basis) where {T<:AbstractOperator,S} = throw(ArgumentError("Identity operator not defined for operator type $T."))
-identityoperator(::Type{T}, ::Type{S}, b::Basis) where {T<:AbstractOperator,S} = identityoperator(T,S,b,b)
-identityoperator(::Type{T}, bases::Basis...) where T<:AbstractOperator = identityoperator(T,eltype(T),bases...)
-identityoperator(op::T) where {T<:AbstractOperator} = identityoperator(T, op.basis_l, op.basis_r)
-
-# Catch case where eltype cannot be inferred from type; this is a bit hacky
-identityoperator(::Type{T}, ::Type{Any}, b1::Basis, b2::Basis) where T<:AbstractOperator = identityoperator(T, ComplexF64, b1, b2)
-
-one(x::Union{<:Basis,<:AbstractOperator}) = identityoperator(x)
 
 # Helper functions to check validity of arguments
 function check_ptrace_arguments(a::AbstractOperator, indices)
@@ -383,17 +165,5 @@ function check_ptrace_arguments(a::StateVector, indices)
     check_indices(length(basis(a).shape), indices)
 end
 
-samebases(a::AbstractOperator) = samebases(a.basis_l, a.basis_r)::Bool
-samebases(a::AbstractOperator, b::AbstractOperator) = samebases(a.basis_l, b.basis_l)::Bool && samebases(a.basis_r, b.basis_r)::Bool
-check_samebases(a::AbstractOperator) = check_samebases(a.basis_l, a.basis_r)
-
 multiplicable(a::AbstractOperator, b::Ket) = multiplicable(a.basis_r, b.basis)
 multiplicable(a::Bra, b::AbstractOperator) = multiplicable(a.basis, b.basis_l)
-multiplicable(a::AbstractOperator, b::AbstractOperator) = multiplicable(a.basis_r, b.basis_l)
-
-Base.size(op::AbstractOperator) = (length(op.basis_l),length(op.basis_r))
-function Base.size(op::AbstractOperator, i::Int)
-    i < 1 && throw(ErrorException("dimension index is < 1"))
-    i > 2 && return 1
-    i==1 ? length(op.basis_l) : length(op.basis_r)
-end
diff --git a/src/operators_dense.jl b/src/operators_dense.jl
index 4184fec3..d1e9a23b 100644
--- a/src/operators_dense.jl
+++ b/src/operators_dense.jl
@@ -204,9 +204,6 @@ function ptrace(psi::Bra, indices)
     return Operator(b_, b_, result)
 end
 
-_index_complement(b::CompositeBasis, indices) = complement(length(b.bases), indices)
-reduced(a, indices) = ptrace(a, _index_complement(basis(a), indices))
-
 normalize!(op::Operator) = (rmul!(op.data, 1.0/tr(op)); op)
 
 function expect(op::DataOperator{B,B}, state::Ket{B}) where B
diff --git a/src/operators_sparse.jl b/src/operators_sparse.jl
index 168488f0..7c12a240 100644
--- a/src/operators_sparse.jl
+++ b/src/operators_sparse.jl
@@ -24,12 +24,6 @@ SparseOperator(::Type{T},b::Basis) where T = SparseOperator(b,b,spzeros(T,length
 SparseOperator(b1::Basis, b2::Basis) = SparseOperator(ComplexF64, b1, b2)
 SparseOperator(b::Basis) = SparseOperator(ComplexF64, b, b)
 
-"""
-    sparse(op::AbstractOperator)
-
-Convert an arbitrary operator into a [`SparseOperator`](@ref).
-"""
-sparse(a::AbstractOperator) = throw(ArgumentError("Direct conversion from $(typeof(a)) not implemented. Use sparse(full(op)) instead."))
 sparse(a::DataOperator) = Operator(a.basis_l, a.basis_r, sparse(a.data))
 
 function ptrace(op::SparseOpPureType, indices)
@@ -56,7 +50,7 @@ function permutesystems(rho::SparseOpPureType{B1,B2}, perm) where {B1<:Composite
     @assert length(rho.basis_l.bases) == length(rho.basis_r.bases) == length(perm)
     @assert isperm(perm)
     shape = [rho.basis_l.shape; rho.basis_r.shape]
-    data = permutedims(rho.data, shape, [perm; perm .+ length(perm)])
+    data = _permutedims(rho.data, shape, [perm; perm .+ length(perm)])
     SparseOperator(permutesystems(rho.basis_l, perm), permutesystems(rho.basis_r, perm), data)
 end
 
@@ -74,10 +68,7 @@ identityoperator(::Type{T}, ::Type{S}, b1::Basis, b2::Basis) where {T<:DataOpera
 identityoperator(::Type{T}, ::Type{S}, b::Basis) where {T<:DataOperator,S<:Number} =
     Operator(b, b, Eye{S}(length(b)))
 
-identityoperator(::Type{T}, b1::Basis, b2::Basis) where T<:Number = identityoperator(DataOperator, T, b1, b2)
-identityoperator(::Type{T}, b::Basis) where T<:Number = identityoperator(DataOperator, T, b)
-identityoperator(b1::Basis, b2::Basis) = identityoperator(ComplexF64, b1, b2)
-identityoperator(b::Basis) = identityoperator(ComplexF64, b)
+identityoperator(::Type{T}, b1::Basis, b2::Basis) where T<:Number = identityoperator(DataOperator, T, b1, b2) # XXX This is purposeful type piracy over QuantumInterface, that hardcodes the use of QuantumOpticsBase.DataOperator in identityoperator. Also necessary for backward compatibility.
 
 """
     diagonaloperator(b::Basis)
diff --git a/src/particle.jl b/src/particle.jl
index 31114a95..e0f5ff81 100644
--- a/src/particle.jl
+++ b/src/particle.jl
@@ -1,6 +1,8 @@
-import Base: position
+import Base: position # TODO is it a good idea to abuse the Base namespace for a different use?
 using FFTW
 
+abstract type ParticleBasis <: Basis end
+
 """
     PositionBasis(xmin, xmax, Npoints)
     PositionBasis(b::MomentumBasis)
@@ -16,7 +18,7 @@ of ``x_{min}`` and ``x_{max}`` are due to the periodic boundary conditions
 more or less arbitrary and are chosen to be
 ``-\\pi/dp`` and ``\\pi/dp`` with ``dp=(p_{max}-p_{min})/N``.
 """
-struct PositionBasis{X1,X2,T,F} <: Basis
+struct PositionBasis{X1,X2,T,F} <: ParticleBasis
     shape::Vector{T}
     xmin::F
     xmax::F
@@ -45,7 +47,7 @@ of ``p_{min}`` and ``p_{max}`` are due to the periodic boundary conditions
 more or less arbitrary and are chosen to be
 ``-\\pi/dx`` and ``\\pi/dx`` with ``dx=(x_{max}-x_{min})/N``.
 """
-struct MomentumBasis{P1,P2,T,F} <: Basis
+struct MomentumBasis{P1,P2,T,F} <: ParticleBasis
     shape::Vector{T}
     pmin::F
     pmax::F
@@ -168,7 +170,7 @@ function position(::Type{T}, b::PositionBasis) where T
     d = convert.(T, samplepoints(b))
     return SparseOperator(b, spdiagm(0 => d))
 end
-position(b::Basis) = position(ComplexF64, b)
+position(b::ParticleBasis) = position(ComplexF64, b)
 
 """
     position([T=ComplexF64,] b:MomentumBasis)
diff --git a/src/printing.jl b/src/printing.jl
index a8e7ed6e..95cadc1b 100644
--- a/src/printing.jl
+++ b/src/printing.jl
@@ -19,47 +19,6 @@ function set_printing(; standard_order::Bool=_std_order, rounding_tol::Real=_rou
 end
 set_printing(standard_order=false, rounding_tol=1e-17)
 
-function show(stream::IO, x::GenericBasis)
-    if length(x.shape) == 1
-        write(stream, "Basis(dim=$(x.shape[1]))")
-    else
-        s = replace(string(x.shape), " " => "")
-        write(stream, "Basis(shape=$s)")
-    end
-end
-
-function show(stream::IO, x::CompositeBasis)
-    write(stream, "[")
-    for i in 1:length(x.bases)
-        show(stream, x.bases[i])
-        if i != length(x.bases)
-            write(stream, " ⊗ ")
-        end
-    end
-    write(stream, "]")
-end
-
-function show(stream::IO, x::SpinBasis)
-    d = denominator(x.spinnumber)
-    n = numerator(x.spinnumber)
-    if d == 1
-        write(stream, "Spin($n)")
-    else
-        write(stream, "Spin($n/$d)")
-    end
-end
-
-function show(stream::IO, x::FockBasis)
-    if iszero(x.offset)
-        write(stream, "Fock(cutoff=$(x.N))")
-    else
-        write(stream, "Fock(cutoff=$(x.N), offset=$(x.offset))")
-    end
-end
-
-function show(stream::IO, x::NLevelBasis)
-    write(stream, "NLevel(N=$(x.N))")
-end
 
 function show(stream::IO, x::PositionBasis)
     write(stream, "Position(xmin=$(x.xmin), xmax=$(x.xmax), N=$(x.N))")
@@ -77,17 +36,6 @@ function show(stream::IO, x::ManyBodyBasis)
     write(stream, "ManyBody(onebodybasis=$(x.onebodybasis), states:$(length(x.occupations)))")
 end
 
-function show(stream::IO, x::SumBasis)
-    write(stream, "[")
-    for i in 1:length(x.bases)
-        show(stream, x.bases[i])
-        if i != length(x.bases)
-            write(stream, " ⊕ ")
-        end
-    end
-    write(stream, "]")
-end
-
 summary(io::IO, x::Ket) = print(io, "Ket(dim=$(length(x.basis)))\n  basis: $(x.basis)\n")
 function show(stream::IO, x::Ket)
     summary(stream, x)
@@ -108,21 +56,6 @@ function show(stream::IO, x::Bra)
     end
 end
 
-function summary(stream::IO, x::AbstractOperator)
-    print(stream, "$(typeof(x).name.name)(dim=$(length(x.basis_l))x$(length(x.basis_r)))\n")
-    if samebases(x)
-        print(stream, "  basis: ")
-        show(stream, basis(x))
-    else
-        print(stream, "  basis left:  ")
-        show(stream, x.basis_l)
-        print(stream, "\n  basis right: ")
-        show(stream, x.basis_r)
-    end
-end
-
-show(stream::IO, x::AbstractOperator) = summary(stream, x)
-
 function show(stream::IO, x::DataOperator)
     summary(stream, x)
     print(stream, "\n")
diff --git a/src/sortedindices.jl b/src/sortedindices.jl
deleted file mode 100644
index fa71fe00..00000000
--- a/src/sortedindices.jl
+++ /dev/null
@@ -1,117 +0,0 @@
-"""
-6, [1, 4] => [2, 3, 5, 6]
-"""
-function complement(N, indices)
-    L = length(indices)
-    x = Vector{Int}(undef, N - L)
-    i_ = 1 # Position in the x vector
-    j = 1 # Position in indices vector
-    for i=1:N
-        if j > L || indices[j]!=i
-            x[i_] = i
-            i_ += 1
-        else
-            j += 1
-        end
-    end
-    x
-end
-
-
-"""
-[1, 4, 5], [2, 4, 7] => [1, 5]
-"""
-function remove(ind1, ind2)
-    x = Int[]
-    for i in ind1
-        if i ∉ ind2
-            push!(x, i)
-        end
-    end
-    x
-end
-
-"""
-[1, 4, 5], [2, 4, 7] => [1, 3]
-"""
-function shiftremove(ind1, ind2)
-    x = Int[]
-    for i in ind1
-        if i ∉ ind2
-            counter = 0
-            for i2 in ind2
-                if i2 < i
-                    counter += 1
-                else
-                    break
-                end
-            end
-            push!(x, i-counter)
-        end
-    end
-    x
-end
-
-function reducedindices(I_, I)
-    N = length(I_)
-    x = Vector{Int}(undef, N)
-    for n in 1:N
-        x[n] = findfirst(isequal(I_[n]), I)
-    end
-    x
-end
-
-function reducedindices!(I_, I)
-    for n in 1:length(I_)
-        I_[n] = findfirst(isequal(I_[n]), I)
-    end
-end
-
-"""
-Check if all indices are unique and smaller than or equal to imax.
-"""
-function check_indices(imax, indices)
-    N = length(indices)
-    for n=1:N
-        i = indices[n]
-        @assert 0 < i <= imax
-        for m in n+1:N
-            @assert i != indices[m]
-        end
-    end
-end
-
-"""
-Check if the indices are sorted, unique and smaller than or equal to imax.
-"""
-function check_sortedindices(imax, indices)
-    N = length(indices)
-    if N == 0
-        return nothing
-    end
-    i_ = indices[1]
-    @assert 0 < i_ <= imax
-    for i in indices[2:end]
-        @assert 0 < i <= imax
-        @assert i > i_
-    end
-end
-
-"""
-    check_embed_indices(indices::Array)
-
-Determine whether a collection of indices, written as a list of (integers or lists of integers) is unique.
-This assures that the embedded operators are in non-overlapping subspaces.
-"""
-function check_embed_indices(indices)
-    # short circuit return when `indices` is empty.
-    length(indices) == 0 && return true
-
-    err_str = "Variable `indices` comes in an unexpected form. Expecting `Array{Union{Int, Array{Int, 1}}, 1}`"
-    @assert all(x isa Array || x isa Int for x in indices) err_str
-
-    # flatten the indices and check for uniqueness
-    # use a custom flatten because it's ≈ 4x  faster than Base.Iterators.flatten
-    flatten(arr) = mapreduce(x -> x isa Vector ? flatten(x) : x, append!, arr, init=[])
-    allunique(flatten(indices))
-end
diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
index 19981228..45e08c97 100644
--- a/src/sparsematrix.jl
+++ b/src/sparsematrix.jl
@@ -239,31 +239,10 @@ function gemv!(alpha, v::AbstractVector, M::SparseMatrixCSC, beta, result::Abstr
     end
 end
 
-function sub2sub(shape1, shape2, I)
-    linearindex = LinearIndices(shape1)[I.I...]
-    CartesianIndices(shape2)[linearindex]
-end
 
-function ptrace(x, shape_nd, indices)
-    shape_nd = (shape_nd...,)
-    N = div(length(shape_nd), 2)
-    shape_2d = (x.m, x.n)
-    shape_nd_after = ([i ∈ indices || i-N ∈ indices ? 1 : shape_nd[i] for i=1:2*N]...,)
-    shape_2d_after = (prod(shape_nd_after[1:N]), prod(shape_nd_after[N+1:end]))
-    I_nd_after_max = CartesianIndex(shape_nd_after...)
-    y = spzeros(eltype(x), shape_2d_after...)
-    for I in eachindex(x)
-        I_nd = sub2sub(shape_2d, shape_nd, I)
-        if I_nd.I[indices] != I_nd.I[indices .+ N]
-            continue
-        end
-        I_after = sub2sub(shape_nd_after, shape_2d_after, min(I_nd, I_nd_after_max))
-        y[I_after] += x[I]
-    end
-    y
-end
 
-function permutedims(x, shape, perm)
+
+function _permutedims(x::AbstractSparseMatrix, shape, perm) # TODO upstream as `permutedims` to SparseArrays to avoid piracy -- used in a single location in operators_sparse
     shape = (shape...,)
     shape_perm = ([shape[i] for i in perm]...,)
     y = spzeros(eltype(x), x.m, x.n)
diff --git a/src/spinors.jl b/src/spinors.jl
index b78e2e51..f1ab83ae 100644
--- a/src/spinors.jl
+++ b/src/spinors.jl
@@ -1,45 +1,4 @@
-"""
-    SumBasis(b1, b2...)
-
-Similar to [`CompositeBasis`](@ref) but for the [`directsum`](@ref) (⊕)
-"""
-struct SumBasis{S,B} <: Basis
-    shape::S
-    bases::B
-end
-SumBasis(bases) = SumBasis(Int[length(b) for b in bases], bases)
-SumBasis(shape, bases::Vector) = (tmp = (bases...,); SumBasis(shape, tmp))
-SumBasis(bases::Vector) = SumBasis((bases...,))
-SumBasis(bases::Basis...) = SumBasis((bases...,))
-
-==(b1::T, b2::T) where T<:SumBasis = equal_shape(b1.shape, b2.shape)
-==(b1::SumBasis, b2::SumBasis) = false
-length(b::SumBasis) = sum(b.shape)
-
-"""
-    directsum(b1::Basis, b2::Basis)
-
-Construct the [`SumBasis`](@ref) out of two sub-bases.
-"""
-directsum(b1::Basis, b2::Basis) = SumBasis(Int[length(b1); length(b2)], Basis[b1, b2])
-directsum(b::Basis) = b
-directsum(b::Basis...) = reduce(directsum, b)
-function directsum(b1::SumBasis, b2::Basis)
-    shape = [b1.shape;length(b2)]
-    bases = [b1.bases...;b2]
-    return SumBasis(shape, (bases...,))
-end
-function directsum(b1::Basis, b2::SumBasis)
-    shape = [length(b1);b2.shape]
-    bases = [b1;b2.bases...]
-    return SumBasis(shape, (bases...,))
-end
-function directsum(b1::SumBasis, b2::SumBasis)
-    shape = [b1.shape;b2.shape]
-    bases = [b1.bases...;b2.bases...]
-    return SumBasis(shape, (bases...,))
-end
-directsum() = GenericBasis(0)
+using QuantumInterface: SumBasis
 
 """
     directsum(x::Ket, y::Ket)
@@ -50,7 +9,6 @@ basis given by the corresponding [`SumBasis`](@ref).
 **NOTE**: The resulting state is not normalized!
 """
 directsum(x::Ket, y::Ket) = Ket(directsum(x.basis, y.basis), [x.data; y.data])
-directsum(x::StateVector...) = reduce(directsum, x)
 
 """
     getblock(x::Ket{<:SumBasis}, i)
@@ -96,7 +54,6 @@ function directsum(a::SparseOpType, b::SparseOpType)
     data[size(a,1)+1:end, size(a,2)+1:end] = b.data
     return Operator(directsum(a.basis_l, b.basis_l), directsum(a.basis_r, b.basis_r), data)
 end
-directsum(a::AbstractOperator...) = reduce(directsum, a)
 
 """
     setblock!(op::DataOperator{<:SumBasis,<:SumBasis}, val::DataOperator, i, j)
@@ -192,8 +149,6 @@ function embed(basis_l::SumBasis, basis_r::SumBasis,
     return embedded_op
 end
 
-embed(b::SumBasis, indices, ops) = embed(b, b, indices, ops)
-
 
 """
     LazyDirectSum <: AbstractOperator
diff --git a/src/states.jl b/src/states.jl
index 0fb3129e..2bd198d1 100644
--- a/src/states.jl
+++ b/src/states.jl
@@ -51,10 +51,6 @@ Ket{B}(b::B) where B = Ket{B}(ComplexF64, b)
 Bra(b::Basis) = Bra(ComplexF64, b)
 Ket(b::Basis) = Ket(ComplexF64, b)
 
-copy(a::T) where {T<:StateVector} = T(a.basis, copy(a.data))
-length(a::StateVector) = length(a.basis)::Int
-basis(a::StateVector) = a.basis
-
 ==(x::Ket{B}, y::Ket{B}) where {B} = (samebases(x, y) && x.data==y.data)
 ==(x::Bra{B}, y::Bra{B}) where {B} = (samebases(x, y) && x.data==y.data)
 ==(x::Ket, y::Ket) = false
@@ -76,13 +72,10 @@ Base.isapprox(x::Bra, y::Bra; kwargs...) = false
 -(a::Ket, b::Ket) = throw(IncompatibleBases())
 -(a::Bra, b::Bra) = throw(IncompatibleBases())
 
--(a::T) where {T<:StateVector} = T(a.basis, -a.data)
-
 *(a::Bra{B}, b::Ket{B}) where {B} = transpose(a.data)*b.data
 *(a::Bra, b::Ket) = throw(IncompatibleBases())
 *(a::Number, b::Ket) = Ket(b.basis, a*b.data)
 *(a::Number, b::Bra) = Bra(b.basis, a*b.data)
-*(a::StateVector, b::Number) = b*a
 
 /(a::Ket, b::Number) = Ket(a.basis, a.data ./ b)
 /(a::Bra, b::Number) = Bra(a.basis, a.data ./ b)
@@ -104,32 +97,8 @@ Tensor product ``|x⟩⊗|y⟩⊗|z⟩⊗…`` of the given states.
 """
 tensor(a::Ket, b::Ket) = Ket(tensor(a.basis, b.basis), kron(b.data, a.data))
 tensor(a::Bra, b::Bra) = Bra(tensor(a.basis, b.basis), kron(b.data, a.data))
-tensor(state::StateVector) = state
 tensor(states::Ket...) = reduce(tensor, states)
 tensor(states::Bra...) = reduce(tensor, states)
-tensor(states::Vector{T}) where T<:StateVector = reduce(tensor, states)
-
-# Normalization functions
-"""
-    norm(x::StateVector)
-
-Norm of the given bra or ket state.
-"""
-norm(x::StateVector) = norm(x.data)
-
-"""
-    normalize(x::StateVector)
-
-Return the normalized state so that `norm(x)` is one.
-"""
-normalize(x::StateVector) = x/norm(x)
-
-"""
-    normalize!(x::StateVector)
-
-In-place normalization of the given bra or ket so that `norm(x)` is one.
-"""
-normalize!(x::StateVector) = (normalize!(x.data); x)
 
 function permutesystems(state::T, perm) where T<:Ket
     @assert length(state.basis.bases) == length(perm)
@@ -201,16 +170,6 @@ end
 samebases(a::Ket{B}, b::Ket{B}) where {B} = samebases(a.basis, b.basis)::Bool
 samebases(a::Bra{B}, b::Bra{B}) where {B} = samebases(a.basis, b.basis)::Bool
 
-# Array-like functions
-Base.size(x::StateVector) = size(x.data)
-@inline Base.axes(x::StateVector) = axes(x.data)
-Base.ndims(x::StateVector) = 1
-Base.ndims(::Type{<:StateVector}) = 1
-Base.eltype(x::StateVector) = eltype(x.data)
-
-# Broadcasting
-Base.broadcastable(x::StateVector) = x
-
 # Custom broadcasting style
 abstract type StateVectorStyle{B} <: Broadcast.BroadcastStyle end
 struct KetStyle{B} <: StateVectorStyle{B} end
diff --git a/test/test_aqua.jl b/test/test_aqua.jl
index 1e1985f7..d6b49a57 100644
--- a/test/test_aqua.jl
+++ b/test/test_aqua.jl
@@ -5,7 +5,10 @@ using FillArrays
 
 @testset "aqua" begin
     Aqua.test_all(QuantumOpticsBase;
-                  piracy=false,  # TODO: Due to Base methods in QuantumOpticsBase, for types defined in QuantumInterface
-                  ambiguities=(exclude=[FillArrays.reshape],)  # Due to https://github.com/JuliaArrays/FillArrays.jl/issues/105#issuecomment-1518406018
+                  ambiguities=(exclude=[FillArrays.reshape],),  # Due to https://github.com/JuliaArrays/FillArrays.jl/issues/105#issuecomment-1518406018
+                  piracy=(broken=true,)
                   )
+    # manual piracy check to exclude identityoperator
+    pirates = [pirate for pirate in Aqua.Piracy.hunt(QuantumOpticsBase) if pirate.name != :identityoperator]
+    @test isempty(pirates)
 end # testset
diff --git a/test/test_sortedindices.jl b/test/test_sortedindices.jl
index 6c46b420..fba61bf2 100644
--- a/test/test_sortedindices.jl
+++ b/test/test_sortedindices.jl
@@ -1,10 +1,9 @@
 using Test
-using QuantumOpticsBase
-
+using QuantumInterface
 
 @testset "sortedindices" begin
 
-s = QuantumOpticsBase
+s = QuantumInterface
 
 @test s.complement(6, [1, 4]) == [2, 3, 5, 6]