diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0e99bb7..d1b3b60 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,27 @@
 # News
 
+## v0.4.0 - 2025-04-05
+
+This version implements the RFC in #40 without deprecating anything. Future versions will first add deprecation warnings before removing any existing interfaces.
+
+- Eliminate all type parameters from `StateVector`, `AbstractKet`, `AbstractBra`, `AbstractOperator`, and `AbstractSuperOperator`.
+- Implement new basis checking interface consisting of `multiplicable`, `addible`, `check_multiplicable`, `check_addible`, and `@compatiblebases`.
+- Add `basis_l` and `basis_r` to get left and right bases of operators. 
+- Implement `getindex` for `CompositeBasis` and `SumBasis`.
+- Add method interface to existing subtypes of `Bases`: `spinnumber`, `cutoff`, `offset`.
+- Add `KetBraBasis`, `ChoiRefSysBasis`, `ChoiOutSysBasis`, and `_PauliBasis`. Note that eventually `_PauliBasis` will be renamed to `PauliBasis`.
+- Change internal fields and constructors for subtypes of `Basis`.
+
+## v0.3.10 - 2025-04-06
+
+- Deprecate `PauliBasis` as it is identical to `SpinBasis(1//2)^N` but without the compatibility with `CompositeBasis`.
+
+## v0.3.9 - 2025-04-03
+
+- Deprecate `equal_bases`.
+
 ## v0.3.8 - 2025-03-08
+
 - Introduce `metrics.jl` file for metric declarations.
 
 ## v0.3.7 - 2025-01-16
diff --git a/Project.toml b/Project.toml
index b62f57d..9f743ed 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumInterface"
 uuid = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 authors = ["QuantumInterface.jl contributors"]
-version = "0.3.8"
+version = "0.4.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/QuantumInterface.jl b/src/QuantumInterface.jl
index 1b65544..73c48a1 100644
--- a/src/QuantumInterface.jl
+++ b/src/QuantumInterface.jl
@@ -88,11 +88,11 @@ function wigner end
 include("bases.jl")
 include("abstract_types.jl")
 
-include("linalg.jl")
 include("tensor.jl")
 include("embed_permute.jl")
 include("expect_variance.jl")
 include("identityoperator.jl")
+include("superoperators.jl")
 
 include("julia_base.jl")
 include("julia_linalg.jl")
@@ -103,4 +103,6 @@ include("express.jl")
 
 include("metrics.jl")
 
+include("deprecated.jl")
+
 end # module
diff --git a/src/abstract_types.jl b/src/abstract_types.jl
index f8667c9..d643084 100644
--- a/src/abstract_types.jl
+++ b/src/abstract_types.jl
@@ -1,47 +1,57 @@
 """
-Abstract base class for `Bra` and `Ket` states.
+Abstract type for all state vectors.
 
-The state vector class stores the coefficients of an abstract state
-in respect to a certain basis. These coefficients are stored in the
-`data` field and the basis is defined in the `basis`
-field.
+This type represents any abstract pure quantum state as given by an element of a
+Hilbert space with respect to a certain basis. To be compatible with methods
+defined in `QuantumInterface`, all subtypes must implement the `basis` method
+which should return a subtype of the abstract [`Basis`](@ref) type.
+
+See also [`AbstractKet`](@ref) and [`AbstractBra`](@ref).
 """
-abstract type StateVector{B,T} end
-abstract type AbstractKet{B,T} <: StateVector{B,T} end
-abstract type AbstractBra{B,T} <: StateVector{B,T} end
+abstract type StateVector end
 
 """
-Abstract base class for all operators.
+Abstract type for `Ket` states.
 
-All deriving operator classes have to define the fields
-`basis_l` and `basis_r` defining the left and right side bases.
+This subtype of [`StateVector`](@ref) is meant to represent `Ket` states which
+are related to their dual `Bra` by the conjugate transpose.
 
-For fast time evolution also at least the function
-`mul!(result::Ket,op::AbstractOperator,x::Ket,alpha,beta)` should be
-implemented. Many other generic multiplication functions can be defined in
-terms of this function and are provided automatically.
+See also [`AbstractBra`](@ref).
 """
-abstract type AbstractOperator{BL,BR} end
+abstract type AbstractKet <: StateVector end
+
+"""
+Abstract type for `Bra` states.
+
+This subtype of [`StateVector`](@ref) is meant to represent `Bra` states which
+are related to their dual `Ket` by the conjugate transpose.
 
+See also [`AbstractBra`](@ref).
 """
-Base class for all super operator classes.
+abstract type AbstractBra <: StateVector end
 
-Super operators are bijective mappings from operators given in one specific
-basis to operators, possibly given in respect to another, different basis.
-To embed super operators in an algebraic framework they are defined with a
-left hand basis `basis_l` and a right hand basis `basis_r` where each of
-them again consists of a left and right hand basis.
-```math
-A_{bl_1,bl_2} = S_{(bl_1,bl_2) ↔ (br_1,br_2)} B_{br_1,br_2}
-\\\\
-A_{br_1,br_2} = B_{bl_1,bl_2} S_{(bl_1,bl_2) ↔ (br_1,br_2)}
-```
 """
-abstract type AbstractSuperOperator{B1,B2} end
+Abstract type for all operators and super operators.
+
+This type represents any abstract mixed quantum state given by a density
+operator (or superoperator) mapping between two Hilbert spaces.  All subtypes
+must implement the [`basis_l`](@ref) and [`basis_r`](@ref) methods which return
+subtypes of [`Basis`](@ref) representing the left and right bases that the
+operator maps between. A subtype is considered compatible with multiplication by
+a subtype of [`AbstractBra`](@ref) defined in same left basis as the operator
+and a subtype of [`AbstractKet`](@ref) defined in the same right basis as the
+operator.
+
+For fast time evolution also at least the function
+`mul!(result::Ket,op::AbstractOperator,x::Ket,alpha,beta)` should be
+implemented. Many other generic multiplication functions can be defined in terms
+of this function and are provided automatically.
+"""
+abstract type AbstractOperator 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)
+    if multiplicable(x,x)
         print(stream, "  basis: ")
         show(stream, basis(x))
     else
diff --git a/src/bases.jl b/src/bases.jl
index 6e4b077..78d5629 100644
--- a/src/bases.jl
+++ b/src/bases.jl
@@ -1,35 +1,76 @@
 """
-Abstract base class for all specialized bases.
+Abstract type for all specialized bases of a Hilbert space.
 
-The Basis class is meant to specify a basis of the Hilbert space of the
-studied system. Besides basis specific information all subclasses must
-implement a shape variable which indicates the dimension of the used
-Hilbert space. For a spin-1/2 Hilbert space this would be the
-vector `[2]`. A system composed of two spins would then have a
-shape vector `[2 2]`.
+This type specifies an orthonormal basis for the Hilbert space of the given
+system. All subtypes must implement `Base.:(==)` and `Base.length`, where the
+latter should return the total dimension of the Hilbert space
 
-Composite systems can be defined with help of the [`CompositeBasis`](@ref)
-class.
+Composite systems can be defined with help of [`CompositeBasis`](@ref).
+
+All relevant properties of concrete subtypes of `Basis` defined in
+`QuantumInterface` should be accessed using their documented functions and
+should not assume anything about the internal representation of instances of
+these types (i.e. do not access the fields of the structs directly).
 """
 abstract type Basis end
 
+"""
+    basis(a)
+
+Return the basis of a quantum object.
+
+If it's ambiguous, e.g. if an operator has a different left and right basis, an
+[`IncompatibleBases`](@ref) error is thrown.
+
+See [`StateVector`](@ref) and [`AbstractOperator`](@ref)
+"""
+function basis end
+
+"""
+    basis_l(a)
+
+Return the left basis of an operator.
+"""
+function basis_l end
+
+"""
+    basis_r(a)
+
+Return the right basis of an operator.
+"""
+function basis_r end
+
+
 """
     length(b::Basis)
 
 Total dimension of the Hilbert space.
 """
-Base.length(b::Basis) = prod(b.shape)
+Base.length(b::Basis) = throw(ArgumentError("Base.length() is not defined for $(typeof(b))"))
 
 """
-    basis(a)
+    size(b::Basis)
 
-Return the basis of an object.
+A vector containing the local dimensions of each Hilbert space in its tensor
+product decomposition into subsystems.
 
-If it's ambiguous, e.g. if an operator has a different left and right basis,
-an [`IncompatibleBases`](@ref) error is thrown.
+See also [`nsubsystems`](@ref) and [`CompositeBasis`](@ref).
 """
-function basis end
+Base.size(b::Basis) = [length(b)]
+
+"""
+    getindex(b::Basis)
 
+Get the i'th factor in the tensor product decomposition of the basis into
+subsystems.
+
+See also [`nsubsystems`](@ref) and [`CompositeBasis`](@ref).
+"""
+Base.getindex(b::Basis, i) = i==1 ? b : throw(BoundsError("attempted to access a nonexistent subsystem basis"))
+
+##
+# GenericBasis, CompositeBasis, SumBasis
+##
 
 """
     GenericBasis(N)
@@ -41,33 +82,32 @@ bases of state vectors and operators are correct for algebraic operations.
 The preferred way is to specify special bases for different systems.
 """
 struct GenericBasis{S} <: Basis
-    shape::S
+    dim::S
 end
-GenericBasis(N::Integer) = GenericBasis([N])
-
-Base.:(==)(b1::GenericBasis, b2::GenericBasis) = equal_shape(b1.shape, b2.shape)
-
+Base.:(==)(b1::GenericBasis, b2::GenericBasis) = b1.dim == b2.dim
+Base.length(b::GenericBasis) = b.dim
 
 """
     CompositeBasis(b1, b2...)
 
 Basis for composite Hilbert spaces.
 
-Stores the subbases in a vector and creates the shape vector directly
-from the shape vectors of these subbases. Instead of creating a CompositeBasis
-directly `tensor(b1, b2...)` or `b1 ⊗ b2 ⊗ …` can be used.
+Stores the subbases in a vector and creates the shape vector directly from the
+dimensions of these subbases. Instead of creating a CompositeBasis directly,
+`tensor(b1, b2...)` or `b1 ⊗ b2 ⊗ …` should be used.
 """
-struct CompositeBasis{S,B} <: Basis
-    shape::S
-    bases::B
+struct CompositeBasis{B<:Basis,S<:Integer} <: Basis
+    shape::Vector{S}
+    bases::Vector{B}
 end
-CompositeBasis(bases) = CompositeBasis([length(b) for b ∈ bases], bases)
-CompositeBasis(bases::Basis...) = CompositeBasis((bases...,))
-CompositeBasis(bases::Vector) = CompositeBasis((bases...,))
-
-Base.:(==)(b1::T, b2::T) where T<:CompositeBasis = equal_shape(b1.shape, b2.shape)
+CompositeBasis(bases) = CompositeBasis([length(b) for b in bases], bases)
+CompositeBasis(bases::Basis...) = CompositeBasis([bases...])
+CompositeBasis(bases::Tuple) = CompositeBasis([bases...])
 
-tensor(b::Basis) = b
+Base.:(==)(b1::CompositeBasis, b2::CompositeBasis) = all(((i, j),) -> i == j, zip(b1.bases, b2.bases))
+Base.length(b::CompositeBasis) = prod(b.shape)
+Base.size(b::CompositeBasis) = b.shape
+Base.getindex(b::CompositeBasis, i) = b.bases[i]
 
 """
     tensor(x::Basis, y::Basis, z::Basis...)
@@ -77,21 +117,12 @@ Create a [`CompositeBasis`](@ref) from the given bases.
 Any given CompositeBasis is expanded so that the resulting CompositeBasis never
 contains another CompositeBasis.
 """
-tensor(b1::Basis, b2::Basis) = CompositeBasis([length(b1); length(b2)], (b1, b2))
-tensor(b1::CompositeBasis, b2::CompositeBasis) = CompositeBasis([b1.shape; b2.shape], (b1.bases..., b2.bases...))
-function tensor(b1::CompositeBasis, b2::Basis)
-    N = length(b1.bases)
-    shape = vcat(b1.shape, length(b2))
-    bases = (b1.bases..., b2)
-    CompositeBasis(shape, bases)
-end
-function tensor(b1::Basis, b2::CompositeBasis)
-    N = length(b2.bases)
-    shape = vcat(length(b1), b2.shape)
-    bases = (b1, b2.bases...)
-    CompositeBasis(shape, bases)
-end
+tensor(b1::Basis, b2::Basis) = CompositeBasis([length(b1), length(b2)], [b1, b2])
+tensor(b1::CompositeBasis, b2::CompositeBasis) = CompositeBasis([b1.shape; b2.shape], [b1.bases; b2.bases])
+tensor(b1::CompositeBasis, b2::Basis) = CompositeBasis([b1.shape; length(b2)], [b1.bases; b2])
+tensor(b1::Basis, b2::CompositeBasis) = CompositeBasis([length(b1); b2.shape], [b1; b2.bases])
 tensor(bases::Basis...) = reduce(tensor, bases)
+tensor(basis::Basis) = basis
 
 function Base.:^(b::Basis, N::Integer)
     if N < 1
@@ -101,41 +132,39 @@ function Base.:^(b::Basis, N::Integer)
 end
 
 """
-    equal_shape(a, b)
+    SumBasis(b1, b2...)
 
-Check if two shape vectors are the same.
+Similar to [`CompositeBasis`](@ref) but for the [`directsum`](@ref) (⊕)
 """
-function equal_shape(a, b)
-    if a === b
-        return true
-    end
-    if length(a) != length(b)
-        return false
-    end
-    for i=1:length(a)
-        if a[i]!=b[i]
-            return false
-        end
-    end
-    return true
+struct SumBasis{S<:Integer,B<:Basis} <: Basis
+    shape::Vector{S}
+    bases::Vector{B}
 end
+SumBasis(bases) = SumBasis([length(b) for b in bases], bases)
+SumBasis(bases::Basis...) = SumBasis([bases...])
+SumBasis(bases::Tuple) = SumBasis([bases...])
+
+Base.:(==)(b1::SumBasis, b2::SumBasis) = all(((i, j),) -> i == j, zip(b1.bases, b2.bases))
+Base.length(b::SumBasis) = sum(b.shape)
+Base.getindex(b::SumBasis, i) = getindex(b.bases, i)
 
 """
-    equal_bases(a, b)
+    directsum(b1::Basis, b2::Basis)
 
-Check if two subbases vectors are identical.
+Construct the [`SumBasis`](@ref) out of two sub-bases.
 """
-function equal_bases(a, b)
-    if a===b
-        return true
-    end
-    for i=1:length(a)
-        if a[i]!=b[i]
-            return false
-        end
-    end
-    return true
-end
+directsum(b1::Basis, b2::Basis) = SumBasis([length(b1), length(b2)], [b1, b2])
+directsum(b1::SumBasis, b2::SumBasis) = SumBasis([b1.shape, b2.shape], [b1.bases; b2.bases])
+directsum(b1::SumBasis, b2::Basis) = SumBasis([b1.shape; length(b2)], [b1.bases; b2])
+directsum(b1::Basis, b2::SumBasis) = SumBasis([length(b1); b2.shape], [b1; b2.bases])
+directsum(bases::Basis...) = reduce(directsum, bases)
+directsum(basis::Basis) = basis
+
+embed(b::SumBasis, indices, ops) = embed(b, b, indices, ops)
+
+##
+# Basis checks
+##
 
 """
 Exception that should be raised for an illegal algebraic operation.
@@ -145,33 +174,33 @@ mutable struct IncompatibleBases <: Exception end
 const BASES_CHECK = Ref(true)
 
 """
-    @samebases
+    @compatiblebases
 
-Macro to skip checks for same bases. Useful for `*`, `expect` and similar
+Macro to skip checks for compatible bases. Useful for `*`, `expect` and similar
 functions.
 """
-macro samebases(ex)
+macro compatiblebases(ex)
     return quote
-        BASES_CHECK.x = false
+        BASES_CHECK[] = false
         local val = $(esc(ex))
-        BASES_CHECK.x = true
+        BASES_CHECK[] = true
         val
     end
 end
 
 """
-    samebases(a, b)
+    samebases(b1::Basis, b2::Basis)
 
-Test if two objects have the same bases.
+Test if two bases are the same. Equivalant to `==`. See
+[`check_samebases`](@ref).
 """
 samebases(b1::Basis, b2::Basis) = b1==b2
-samebases(b1::Tuple{Basis, Basis}, b2::Tuple{Basis, Basis}) = b1==b2 # for checking superoperators
 
 """
     check_samebases(a, b)
 
-Throw an [`IncompatibleBases`](@ref) error if the objects don't have
-the same bases.
+Throw an [`IncompatibleBases`](@ref) error if the bases are not the same. See
+[`samebases`](@ref).
 """
 function check_samebases(b1, b2)
     if BASES_CHECK[] && !samebases(b1, b2)
@@ -179,34 +208,28 @@ function check_samebases(b1, b2)
     end
 end
 
-
 """
-    multiplicable(a, b)
+    check_addible(a, b)
 
-Check if two objects are multiplicable.
+Throw an [`IncompatibleBases`](@ref) error if the objects are not addible as
+determined by `addible(a, b)`.  Disabled by use of [`@compatiblebases`](@ref)
+anywhere further up in the call stack.
 """
-multiplicable(b1::Basis, b2::Basis) = b1==b2
-
-function multiplicable(b1::CompositeBasis, b2::CompositeBasis)
-    if !equal_shape(b1.shape,b2.shape)
-        return false
-    end
-    for i=1:length(b1.shape)
-        if !multiplicable(b1.bases[i], b2.bases[i])
-            return false
-        end
+function check_addible(a, b)
+    if BASES_CHECK[] && !addible(a, b)
+        throw(IncompatibleBases())
     end
-    return true
 end
 
 """
     check_multiplicable(a, b)
 
-Throw an [`IncompatibleBases`](@ref) error if the objects are
-not multiplicable.
+Throw an [`IncompatibleBases`](@ref) error if the objects are not multiplicable
+as determined by `multiplicable(a, b)`.  Disabled by use of
+[`@compatiblebases`](@ref) anywhere further up in the call stack.
 """
-function check_multiplicable(b1, b2)
-    if BASES_CHECK[] && !multiplicable(b1, b2)
+function check_multiplicable(a, b)
+    if BASES_CHECK[] && !multiplicable(a, b)
         throw(IncompatibleBases())
     end
 end
@@ -245,6 +268,8 @@ function ptrace(b::CompositeBasis, indices)
     reduced(b, J)
 end
 
+_index_complement(b::CompositeBasis, indices) = complement(length(b.bases), indices)
+reduced(a, indices) = ptrace(a, _index_complement(basis(a), indices))
 
 """
     permutesystems(a, perm)
@@ -255,8 +280,8 @@ For a permutation vector `[2,1,3]` and a given object with basis `[b1, b2, b3]`
 this function results in `[b2, b1, b3]`.
 """
 function permutesystems(b::CompositeBasis, perm)
-    @assert length(b.bases) == length(perm)
-    @assert isperm(perm)
+    (nsubsystems(b) == length(perm)) || throw(ArgumentError("Must have nsubsystems(b) == length(perm) in permutesystems"))
+    isperm(perm) || throw(ArgumentError("Must pass actual permeutation to permutesystems"))
     CompositeBasis(b.shape[perm], b.bases[perm])
 end
 
@@ -269,22 +294,42 @@ end
     FockBasis(N,offset=0)
 
 Basis for a Fock space where `N` specifies a cutoff, i.e. what the highest
-included fock state is. Similarly, the `offset` defines the lowest included
-fock state (default is 0). Note that the dimension of this basis is `N+1-offset`.
+included fock state is. Similarly, the `offset` defines the lowest included fock
+state (default is 0). Note that the dimension of this basis is `N+1-offset`.
+The [`cutoff`](@ref) and [`offset`](@ref) functions can be used to obtain the
+respective properties of a given `FockBasis`.
 """
-struct FockBasis{T} <: Basis
-    shape::Vector{T}
+struct FockBasis{T<:Integer} <: Basis
     N::T
     offset::T
     function FockBasis(N::T,offset::T=0) where T
         if N < 0 || offset < 0 || N <= offset
             throw(DimensionMismatch())
         end
-        new{T}([N-offset+1], N, offset)
+        new{T}(N, offset)
     end
 end
 
 Base.:(==)(b1::FockBasis, b2::FockBasis) = (b1.N==b2.N && b1.offset==b2.offset)
+Base.length(b::FockBasis) = b.N - b.offset + 1
+
+"""
+    cutoff(b::FockBasis)
+
+Return the fock cutoff of the given fock basis.
+
+See [`FockBasis`](@ref).
+"""
+cutoff(b::FockBasis) = b.N
+
+"""
+    offset(b::FockBasis)
+
+Return the offset of the given fock basis.
+
+See [`FockBasis`](@ref).
+"""
+offset(b::FockBasis) = b.offset
 
 
 """
@@ -292,121 +337,127 @@ Base.:(==)(b1::FockBasis, b2::FockBasis) = (b1.N==b2.N && b1.offset==b2.offset)
 
 Basis for a system consisting of N states.
 """
-struct NLevelBasis{T} <: Basis
-    shape::Vector{T}
+struct NLevelBasis{T<:Integer} <: Basis
     N::T
-    function NLevelBasis(N::T) where T<:Integer
+    function NLevelBasis(N::T) where T
         if N < 1
             throw(DimensionMismatch())
         end
-        new{T}([N], N)
+        new{T}(N)
     end
 end
 
 Base.:(==)(b1::NLevelBasis, b2::NLevelBasis) = b1.N == b2.N
-
-
-"""
-    PauliBasis(num_qubits::Int)
-
-Basis for an N-qubit space where `num_qubits` specifies the number of qubits.
-The dimension of the basis is 2²ᴺ.
-"""
-struct PauliBasis{S,B} <: Basis
-    shape::S
-    bases::B
-    function PauliBasis(num_qubits::T) where {T<:Integer}
-        shape = [2 for _ in 1:num_qubits]
-        bases = Tuple(SpinBasis(1//2) for _ in 1:num_qubits)
-        return new{typeof(shape),typeof(bases)}(shape, bases)
-    end
-end
-
-Base.:(==)(pb1::PauliBasis, pb2::PauliBasis) = length(pb1.bases) == length(pb2.bases)
-
+Base.length(b::NLevelBasis) = b.N
 
 """
     SpinBasis(n)
 
 Basis for spin-n particles.
 
-The basis can be created for arbitrary spinnumbers by using a rational number,
-e.g. `SpinBasis(3//2)`. The Pauli operators are defined for all possible
-spin numbers.
+The basis can be created for arbitrary spin numbers by using a rational number,
+e.g. `SpinBasis(3//2)`. The Pauli operators are defined for all possible spin
+numbers. The [`spinnumber`](@ref) function can be used to get the spin number
+for a `SpinBasis`.
 """
-struct SpinBasis{S,T} <: Basis
-    shape::Vector{T}
+struct SpinBasis{T<:Integer} <: Basis
     spinnumber::Rational{T}
-    function SpinBasis{S}(spinnumber::Rational{T}) where {S,T<:Integer}
+    function SpinBasis(spinnumber::Rational{T}) where T
         n = numerator(spinnumber)
         d = denominator(spinnumber)
-        @assert d==2 || d==1
-        @assert n >= 0
+        d==2 || d==1 || throw(ArgumentError("Can only construct integer or half-integer spin basis"))
+        n >= 0 || throw(ArgumentError("Can only construct positive spin basis"))
         N = numerator(spinnumber*2 + 1)
-        new{spinnumber,T}([N], spinnumber)
+        new{T}(spinnumber)
     end
 end
-SpinBasis(spinnumber::Rational) = SpinBasis{spinnumber}(spinnumber)
 SpinBasis(spinnumber) = SpinBasis(convert(Rational{Int}, spinnumber))
 
 Base.:(==)(b1::SpinBasis, b2::SpinBasis) = b1.spinnumber==b2.spinnumber
-
+Base.length(b::SpinBasis) = numerator(b.spinnumber*2 + 1)
 
 """
-    SumBasis(b1, b2...)
+    spinnumber(b::SpinBasis)
 
-Similar to [`CompositeBasis`](@ref) but for the [`directsum`](@ref) (⊕)
+Return the spin number of the given spin basis.
+
+See [`SpinBasis`](@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...,))
+spinnumber(b::SpinBasis) = b.spinnumber
+
 
-==(b1::T, b2::T) where T<:SumBasis = equal_shape(b1.shape, b2.shape)
-==(b1::SumBasis, b2::SumBasis) = false
-length(b::SumBasis) = sum(b.shape)
+##
+# Operator Bases
+##
 
 """
-    directsum(b1::Basis, b2::Basis)
+    KetBraBasis(BL,BR)
 
-Construct the [`SumBasis`](@ref) out of two sub-bases.
+The "Ket-Bra" operator basis is the standard representation for the left and
+right bases of superoperators. This basis is formed by "vec'ing" the
+outer-product "Ket-Bra" basis for an operator with a left Bra basis and right
+Ket basis which practically means flipping the Bra to a Ket. The operator itself
+is then represented as a "Super-Bra" in this basis and corresponds to
+column-stacking its matrix.
 """
-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...,))
+struct KetBraBasis{BL<:Basis, BR<:Basis} <: Basis
+    left::BL
+    right::BR
 end
-function directsum(b1::Basis, b2::SumBasis)
-    shape = [length(b1);b2.shape]
-    bases = [b1;b2.bases...]
-    return SumBasis(shape, (bases...,))
+KetBraBasis(b::Basis) = KetBraBasis(b,b)
+basis_l(b::KetBraBasis) = b.left
+basis_r(b::KetBraBasis) = b.right
+Base.:(==)(b1::KetBraBasis, b2::KetBraBasis) = (b1.left == b2.left && b1.right == b2.right)
+Base.length(b::KetBraBasis) = length(b.left)*length(b.right)
+
+struct ChoiBasis{BL<:Basis, BR<:Basis} <: Basis
+    ref::BL
+    out::BR
 end
-function directsum(b1::SumBasis, b2::SumBasis)
-    shape = [b1.shape;b2.shape]
-    bases = [b1.bases...;b2.bases...]
-    return SumBasis(shape, (bases...,))
+basis_l(b::ChoiBasis) = b.ref
+basis_r(b::ChoiBasis) = b.out
+Base.:(==)(b1::ChoiBasis, b2::ChoiBasis) = (b1.ref == b2.ref && b1.out == b2.out)
+Base.length(b::ChoiBasis) = length(b.ref)*length(b.out)
+
+"""
+    PauliBasis(N)
+
+The standard Pauli operator basis for an `N` qubit space. This consists of
+tensor products of the Pauli matrices I, X, Y, Z, in that order for each qubit.
+The dimension of the basis is 2²ᴺ.
+"""
+struct PauliBasis{T<:Integer} <: Basis
+    N::T
 end
+Base.:(==)(b1::PauliBasis, b2::PauliBasis) = b1.N == b2.N
+Base.length(b::PauliBasis) = 4^b.N
+
+"""
+    HWPauliBasis(N)
+
+The Hesienberg-Weyl Pauli operator basis consisting of the
+N represents the underlying Hilbert
+space dimension, not the operator basis dimension. For N>2, this representes the
+operator basis formed by the generalized Pauli matrices, also called the clock
+and shift matrices. The ordering is the usual one: when the index is written in
+base-N and thus has only two digits, the least significant bit gives powers of Z
+(the clock matrix), and most significant bit gives powers of X (the shfit matrix).
+"""
+struct HWPauliBasis{T<:Integer} <: Basis
+    shape::Vector{T}
+end
+HWPauliBasis(N::Integer) = HWPauliBasis([N])
+Base.:(==)(b1::HWPauliBasis, b2::HWPauliBasis) = b1.shape == b2.shape
+Base.length(b::HWPauliBasis) = prod(b.shape)
+Base.getindex(b::HWPauliBasis, i) = HWPauliBasis([b.shape[i]])
 
-embed(b::SumBasis, indices, ops) = embed(b, b, indices, ops)
 
 ##
 # show methods
 ##
 
 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
+    write(stream, "Basis(dim=$(x.dim))")
 end
 
 function show(stream::IO, x::CompositeBasis)
@@ -452,3 +503,15 @@ function show(stream::IO, x::SumBasis)
     end
     write(stream, "]")
 end
+
+function show(stream::IO, x::KetBraBasis)
+    write(stream, "KetBra(left=$(x.left), right=$(x.right))")
+end
+
+function show(stream::IO, x::PauliBasis)
+    write(stream, "Pauli(N=$(x.N)")
+end
+
+function show(stream::IO, x::HWPauliBasis)
+    write(stream, "Pauli($(x.shape)")
+end
diff --git a/src/deprecated.jl b/src/deprecated.jl
new file mode 100644
index 0000000..70da1af
--- /dev/null
+++ b/src/deprecated.jl
@@ -0,0 +1,90 @@
+function equal_bases(a, b)
+    Base.depwarn("`==` should be preferred over `equal_bases`!", :equal_bases)
+    if a===b
+        return true
+    end
+    for i=1:length(a)
+        if a[i]!=b[i]
+            return false
+        end
+    end
+    return true
+end
+
+# TODO: figure out how to deprecate abstract type
+abstract type AbstractSuperOperator end
+
+function basis(a::AbstractSuperOperator)
+    Base.depwarn("`AbstractSuperOperator` will be removed in a future version", :equal_shape)
+    (check_samebases(a); a.basis_l[1]) # FIXME issue #12
+end
+
+function check_samebases(a::AbstractSuperOperator)
+    Base.depwarn("`AbstractSuperOperator` will be removed in a future version", :equal_shape)
+    check_samebases(a.basis_l[1], a.basis_r[1]) # FIXME issue #12
+    check_samebases(a.basis_l[2], a.basis_r[2]) # FIXME issue #12
+end
+
+function equal_shape(a, b)
+    Base.depwarn("`equal_shape` will be removed in a future version", :equal_shape)
+    if a === b
+        return true
+    end
+    if length(a) != length(b)
+        return false
+    end
+    for i=1:length(a)
+        if a[i]!=b[i]
+            return false
+        end
+    end
+    return true
+end
+
+# TODO: figure out how to deprecate a macro
+macro samebases(ex)
+    return quote
+        BASES_CHECK.x = false
+        local val = $(esc(ex))
+        BASES_CHECK.x = true
+        val
+    end
+end
+
+function samebases(b1::Tuple{Basis, Basis}, b2::Tuple{Basis, Basis})
+    Base.depwarn("`samebases(b1:Basis, b2:Basis)`, `addible` or `multiplicable` should be used instead of `samebases(b1::Tuple{Basis, Basis}, b2::Tuple{Basis, Basis})`. ", :check_samebases)
+    b1==b2 # for checking superoperators
+end
+
+function samebases(a::AbstractOperator)
+    Base.depwarn("`multiplicable(a,a)` should be used instead of `samebases(a::AbstractOperator)`. ", :check_samebases)
+    samebases(a.basis_l, a.basis_r)::Bool # FIXME issue #12
+end
+
+function samebases(a::AbstractOperator, b::AbstractOperator)
+    Base.depwarn("``addible(a,b)` should be used instead of `samebases(a::AbstractOperator, b::AbstractOperator)`. ", :check_samebases)
+    samebases(a.basis_l, b.basis_l)::Bool && samebases(a.basis_r, b.basis_r)::Bool # FIXME issue #12
+end
+
+function check_samebases(a::Union{AbstractOperator, AbstractSuperOperator})
+    Base.depwarn("`check_multiplicable(a,a)` should be used instead of `check_samebases(a::Union{AbstractOperator, AbstractSuperOperator})`. ", :check_samebases)
+    check_samebases(a.basis_l, a.basis_r) # FIXME issue #12
+end
+
+function multiplicable(b1::Basis, b2::Basis)
+    Base.depwarn("`==` between bases should be used instead of `multiplicable`.", :check_samebases)
+    b1==b2
+end
+
+function multiplicable(b1::CompositeBasis, b2::CompositeBasis)
+    Base.depwarn("`==` between bases should be used instead of `multiplicable`.", :check_samebases)
+    if !equal_shape(b1.shape,b2.shape)
+        return false
+    end
+    for i=1:length(b1.shape)
+        if !multiplicable(b1.bases[i], b2.bases[i])
+            return false
+        end
+    end
+    return true
+end
diff --git a/src/embed_permute.jl b/src/embed_permute.jl
index c2cc4ca..2a6be55 100644
--- a/src/embed_permute.jl
+++ b/src/embed_permute.jl
@@ -5,12 +5,12 @@
 `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,
+function embed(bl::CompositeBasis, br::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
+    (nsubsystems(bl) == nsubsystems(br)) || throw(ArgumentError("Must have nsubsystems(bl) == nsubsystems(br) in embed"))
+    N = nsubsystems(bl)::Int # type assertion to help type inference
     if length(operators) == 0
-        return identityoperator(T, basis_l, basis_r)
+        return identityoperator(T, bl, br)
     end
     indices, operator_list = zip(operators...)
     operator_list = [operator_list...;]
@@ -22,14 +22,14 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
     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]))
+                push!(operators_flat, identityoperator(T, S, bl[i], br[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]
+        complement_operators = [identityoperator(T, S, bl[i], br[i]) for i in complement_indices_flat]
         op = tensor([operator_list; complement_operators]...)
         perm = sortperm([indices_flat; complement_indices_flat])
         return permutesystems(op, perm)
@@ -52,14 +52,14 @@ embed(basis::CompositeBasis, indices, op::AbstractOperator) = embed(basis, basis
 
 Tensor product of operators where missing indices are filled up with identity operators.
 """
-function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
+function embed(bl::CompositeBasis, br::CompositeBasis,
                indices, operators::Vector{T}) where T<:AbstractOperator
 
-    @assert check_embed_indices(indices)
+    check_embed_indices(indices) || throw(ArgumentError("Must have unique indices in embed"))
+    (nsubsystems(bl) == nsubsystems(br)) || throw(ArgumentError("Must have nsubsystems(bl) == nsubsystems(br) in embed"))
+    (length(indices) == length(operators)) || throw(ArgumentError("Must have length(indices) == length(operators) in embed"))
 
-    N = length(basis_l.bases)
-    @assert length(basis_r.bases) == N
-    @assert length(indices) == length(operators)
+    N = nsubsystems(bl)
 
     # Embed all single-subspace operators.
     idxop_sb = [x for x in zip(indices, operators) if x[1] isa Integer]
@@ -67,17 +67,17 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
     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())
+        (basis_l(opsb) == bl[idxsb]) || throw(IncompatibleBases())
+        (basis_r(opsb) == br[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_op = tensor([i ∈ indices_sb ? ops_sb[indexin(i, indices_sb)[1]] : identityoperator(T, S, bl[i], br[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)
+        embed_op *= embed(bl, br, idxs, op)
     end
 
     return embed_op
@@ -85,8 +85,26 @@ end
 
 permutesystems(a::AbstractOperator, perm) = arithmetic_unary_error("Permutations of subsystems", a)
 
-nsubsystems(s::AbstractKet) = nsubsystems(basis(s))
+"""
+    nsubsystems(a)
+
+Return the number of subsystems of a quantum object in its tensor product
+decomposition.
+
+See also [`CompositeBasis`](@ref).
+"""
+nsubsystems(s::StateVector) = nsubsystems(basis(s))
 nsubsystems(s::AbstractOperator) = nsubsystems(basis(s))
 nsubsystems(b::CompositeBasis) = length(b.bases)
 nsubsystems(b::Basis) = 1
 nsubsystems(::Nothing) = 1 # TODO Exists because of QuantumSavory; Consider removing this and reworking the functions that depend on it. E.g., a reason to have it when performing a project_traceout measurement on a state that contains only one subsystem
+
+"""
+    nsubspaces(a)
+
+Return the number of subspaces of a quantum object in its direct sum
+decomposition.
+
+See also [`SumBasis`](@ref).
+"""
+nsubspaces(b::SumBasis) = length(b.bases)
diff --git a/src/expect_variance.jl b/src/expect_variance.jl
index 86f5e64..678fb5f 100644
--- a/src/expect_variance.jl
+++ b/src/expect_variance.jl
@@ -3,33 +3,35 @@
 
 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
+expect(indices, op::AbstractOperator, state::AbstractOperator) =
+    expect(op, ptrace(state, complement(nsubsystems(state), indices)))
+
+expect(index::Integer, op::AbstractOperator, state::AbstractOperator) = expect([index], op, state)
 
-expect(index::Integer, op::AbstractOperator{B1,B2}, state::AbstractOperator{B3,B3}) where {B1,B2,B3<: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]
 
-expect(op::AbstractOperator{B1,B2}, state::AbstractOperator{B2,B2}) where {B1,B2} = tr(op*state)
+expect(op::AbstractOperator, state::AbstractOperator) =
+    (check_multiplicable(state, state); check_multiplicable(op,state); tr(op*state))
 
 """
     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
+variance(indices, op::AbstractOperator, state::AbstractOperator) =
+    variance(op, ptrace(state, complement(nsubsystems(state), indices)))
+
+variance(index::Integer, op::AbstractOperator, state::AbstractOperator) = variance([index], op, state)
 
-variance(index::Integer, op::AbstractOperator{B,B}, state::AbstractOperator{BC,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]
 
-function variance(op::AbstractOperator{B,B}, state::AbstractOperator{B,B}) where B
-    expect(op*op, state) - expect(op, state)^2
+function variance(op::AbstractOperator, state::AbstractOperator)
+    check_multiplicable(op,op)
+    check_multiplicable(state,state)
+    check_multiplicable(op,state)
+    @compatiblebases expect(op*op, state) - expect(op, state)^2
 end
diff --git a/src/julia_base.jl b/src/julia_base.jl
index 9a0532d..1c0b1c4 100644
--- a/src/julia_base.jl
+++ b/src/julia_base.jl
@@ -8,33 +8,38 @@ addnumbererror() = throw(ArgumentError("Can't add or subtract a number and an op
 # States
 ##
 
--(a::T) where {T<:StateVector} = T(a.basis, -a.data)
+==(a::AbstractKet, b::AbstractBra) = false
+==(a::AbstractBra, b::AbstractKet) = false
+-(a::T) where {T<:StateVector} = T(basis(a), -a.data) # FIXME issue #12
 *(a::StateVector, b::Number) = b*a
-copy(a::T) where {T<:StateVector} = T(a.basis, copy(a.data))
-length(a::StateVector) = length(a.basis)::Int
-basis(a::StateVector) = a.basis
+copy(a::T) where {T<:StateVector} = T(basis(a), copy(a.data)) # FIXME issue #12
+length(a::StateVector) = length(basis(a))::Int
+basis(a::StateVector) = throw(ArgumentError("basis() is not defined for this type of state vector: $(typeof(a))."))
 directsum(x::StateVector...) = reduce(directsum, x)
 
 # Array-like functions
-Base.size(x::StateVector) = size(x.data)
-@inline Base.axes(x::StateVector) = axes(x.data)
+Base.size(x::StateVector) = size(x.data) # FIXME issue #12
+@inline Base.axes(x::StateVector) = axes(x.data) # FIXME issue #12
 Base.ndims(x::StateVector) = 1
 Base.ndims(::Type{<:StateVector}) = 1
-Base.eltype(x::StateVector) = eltype(x.data)
+Base.eltype(x::StateVector) = eltype(x.data) # FIXME issue #12
 
 # Broadcasting
 Base.broadcastable(x::StateVector) = x
 
 Base.adjoint(a::StateVector) = dagger(a)
+dagger(a::StateVector) = arithmetic_unary_error("Hermitian conjugate", a)
 
 
 ##
 # Operators
 ##
 
-length(a::AbstractOperator) = length(a.basis_l)::Int*length(a.basis_r)::Int
-basis(a::AbstractOperator) = (check_samebases(a); a.basis_l)
-basis(a::AbstractSuperOperator) = (check_samebases(a); a.basis_l[1])
+length(a::AbstractOperator) = length(basis_l(a))::Int*length(basis_r(a))::Int
+basis(a::AbstractOperator) = (check_samebases(basis_l(a), basis_r(a)); basis_l(a))
+basis_l(a::AbstractOperator) = throw(ArgumentError("basis_l() is not defined for this type of operator: $(typeof(a))."))
+basis_r(a::AbstractOperator) = throw(ArgumentError("basis_r() is not defined for this type of operator: $(typeof(a))."))
+directsum(a::AbstractOperator...) = reduce(directsum, a)
 
 # Ensure scalar broadcasting
 Base.broadcastable(x::AbstractOperator) = Ref(x)
@@ -53,6 +58,37 @@ Base.broadcastable(x::AbstractOperator) = Ref(x)
 *(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Multiplication", a, b)
 ^(a::AbstractOperator, b::Integer) = Base.power_by_squaring(a, b)
 
+"""
+    addible(a, b)
+
+Check if any two subtypes of `StateVector` or `AbstractOperator`
+ can be added together.
+
+Spcefically this checks whether the left basis of a is equal
+to the left basis of b and whether the right basis of a is equal
+to the right basis of b.
+"""
+addible(a::Union{<:StateVector,<:AbstractOperator},
+        b::Union{<:StateVector,<:AbstractOperator}) = false
+addible(a::AbstractBra, b::AbstractBra) = (basis(a) == basis(b))
+addible(a::AbstractKet, b::AbstractKet) = (basis(a) == basis(b))
+addible(a::AbstractOperator, b::AbstractOperator) =
+    (basis_l(a) == basis_l(b)) && (basis_r(a) == basis_r(b))
+
+
+"""
+    multiplicable(a, b)
+
+Check if any two subtypes of `StateVector` or `AbstractOperator`,
+can be multiplied in the given order.
+"""
+multiplicible(a::Union{<:StateVector,<:AbstractOperator},
+              b::Union{<:StateVector,<:AbstractOperator}) = false
+multiplicable(a::AbstractBra, b::AbstractKet) = (basis(a) == basis(b))
+multiplicable(a::AbstractOperator, b::AbstractKet) = (basis_r(a) == basis(b))
+multiplicable(a::AbstractBra, b::AbstractOperator) = (basis(a) == basis_l(b))
+multiplicable(a::AbstractOperator, b::AbstractOperator) = (basis_r(a) == basis_l(b))
+
 """
     exp(op::AbstractOperator)
 
@@ -60,14 +96,17 @@ 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()."))
 
-Base.size(op::AbstractOperator) = (length(op.basis_l),length(op.basis_r))
+Base.size(op::AbstractOperator) = (length(basis_l(op)),length(basis_r(op)))
 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)
+    i==1 ? length(basis_l(op)) : length(basis_r(op))
 end
 
 Base.adjoint(a::AbstractOperator) = dagger(a)
+dagger(a::AbstractOperator) = arithmetic_unary_error("Hermitian conjugate", a)
 
 conj(a::AbstractOperator) = arithmetic_unary_error("Complex conjugate", a)
 conj!(a::AbstractOperator) = conj(a::AbstractOperator)
+
+ptrace(a::AbstractOperator, index) = arithmetic_unary_error("Partial trace", a)
diff --git a/src/julia_linalg.jl b/src/julia_linalg.jl
index d2f4d3d..1926883 100644
--- a/src/julia_linalg.jl
+++ b/src/julia_linalg.jl
@@ -17,7 +17,7 @@ tr(x::AbstractOperator) = arithmetic_unary_error("Trace", x)
 
 Norm of the given bra or ket state.
 """
-norm(x::StateVector) = norm(x.data)
+norm(x::StateVector) = norm(x.data) # FIXME issue #12
 
 """
     normalize(x::StateVector)
@@ -31,7 +31,7 @@ normalize(x::StateVector) = x/norm(x)
 
 In-place normalization of the given bra or ket so that `norm(x)` is one.
 """
-normalize!(x::StateVector) = (normalize!(x.data); x)
+normalize!(x::StateVector) = (normalize!(x.data); x) # FIXME issue #12
 
 """
     normalize(op)
@@ -46,3 +46,10 @@ normalize(op::AbstractOperator) = op/tr(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()'."))
+
+"""
+    transpose(op)
+
+Transpose of the given operator.
+"""
+transpose(a::AbstractOperator) = arithmetic_unary_error("Transpose", a)
diff --git a/src/linalg.jl b/src/linalg.jl
deleted file mode 100644
index 8bb47cd..0000000
--- a/src/linalg.jl
+++ /dev/null
@@ -1,10 +0,0 @@
-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::Union{AbstractOperator, AbstractSuperOperator}) = check_samebases(a.basis_l, a.basis_r)
-multiplicable(a::AbstractOperator, b::AbstractOperator) = multiplicable(a.basis_r, b.basis_l)
-dagger(a::AbstractOperator) = arithmetic_unary_error("Hermitian conjugate", a)
-transpose(a::AbstractOperator) = arithmetic_unary_error("Transpose", a)
-directsum(a::AbstractOperator...) = reduce(directsum, a)
-ptrace(a::AbstractOperator, index) = arithmetic_unary_error("Partial trace", a)
-_index_complement(b::CompositeBasis, indices) = complement(length(b.bases), indices)
-reduced(a, indices) = ptrace(a, _index_complement(basis(a), indices))
diff --git a/src/sortedindices.jl b/src/sortedindices.jl
index fa71fe0..36b89b5 100644
--- a/src/sortedindices.jl
+++ b/src/sortedindices.jl
@@ -74,9 +74,9 @@ function check_indices(imax, indices)
     N = length(indices)
     for n=1:N
         i = indices[n]
-        @assert 0 < i <= imax
+        (0 < i <= imax) || throw(ArgumentError("indices exceed allowable range"))
         for m in n+1:N
-            @assert i != indices[m]
+            (i != indices[m]) || throw(ArgumentError("indices not unique"))
         end
     end
 end
@@ -90,10 +90,10 @@ function check_sortedindices(imax, indices)
         return nothing
     end
     i_ = indices[1]
-    @assert 0 < i_ <= imax
+    (0 < i_ <= imax) || throw(ArgumentError("indices exceed allowable range"))
     for i in indices[2:end]
-        @assert 0 < i <= imax
-        @assert i > i_
+        (0 < i <= imax) || throw(ArgumentError("indices exceed allowable range"))
+        (i > i_) || throw(ArgumentError("indices not sorted"))
     end
 end
 
@@ -107,8 +107,7 @@ 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
+    all(x isa Array || x isa Int for x in indices) || throw(ArgumentError("Variable `indices` comes in an unexpected form. Expecting `Array{Union{Int, Array{Int, 1}}, 1}`"))
 
     # flatten the indices and check for uniqueness
     # use a custom flatten because it's ≈ 4x  faster than Base.Iterators.flatten
diff --git a/src/superoperators.jl b/src/superoperators.jl
new file mode 100644
index 0000000..0a98f12
--- /dev/null
+++ b/src/superoperators.jl
@@ -0,0 +1,155 @@
+import LinearAlgebra: vec
+
+"""
+    vec(op)
+
+Create a vectorized operator compatible with application with superoperators.
+"""
+vec(op::AbstractOperator) = throw(ArgumentError("vec() is not defined for this type of operator: $(typeof(op))."))
+
+"""
+    unvec(op)
+
+Converta a vectorized operator to a normal operator.
+"""
+unvec(op::AbstractOperator) = throw(ArgumentError("unvec() is not defined for this type of operator: $(typeof(op))."))
+
+"""
+    super(op)
+
+Converts to superoperator representation
+"""
+super(op::AbstractOperator) = throw(ArgumentError("super() is not defined for this type of operator: $(typeof(op))."))
+
+"""
+    choi(op)
+
+Converts to choi state representation
+"""
+choi(op::AbstractOperator) = throw(ArgumentError("choi() is not defined for this type of operator: $(typeof(op))."))
+
+"""
+    kraus(op)
+
+Converts to kraus operator sum representation
+"""
+kraus(op::AbstractOperator) = throw(ArgumentError("kraus() is not defined for this type of operator: $(typeof(op))."))
+
+"""
+    stinespring(op)
+
+Converts to (Pauli) chi process matrix representation
+"""
+stinespring(op::AbstractOperator) = throw(ArgumentError("stinespring() is not defined for this type of operator: $(typeof(op))."))
+
+
+"""
+    pauli(op)
+
+Converts to pauli vector and transfer matrix representation
+"""
+pauli(op::AbstractOperator) = throw(ArgumentError("pauli() is not defined for this type of operator: $(typeof(op))."))
+
+"""
+    chi(op)
+
+Converts to (Pauli) chi process matrix representation
+"""
+chi(op::AbstractOperator) = throw(ArgumentError("chi() is not defined for this type of operator: $(typeof(op))."))
+
+"""
+    spre(op)
+
+Create a super-operator equivalent for right side operator multiplication.
+
+For operators ``A``, ``B`` the relation
+
+```math
+    \\mathrm{spre}(A) B = A B
+```
+
+holds. `op` can be a dense or a sparse operator.
+"""
+spre(op::AbstractOperator) = throw(ArgumentError("spre() is not defined for this type of operator: $(typeof(op))."))
+
+"""
+    spost(op)
+
+Create a super-operator equivalent for left side operator multiplication.
+
+For operators ``A``, ``B`` the relation
+
+```math
+    \\mathrm{spost}(A) B = B A
+```
+
+holds. `op` can be a dense or a sparse operator.
+"""
+spost(op::AbstractOperator) = throw(ArgumentError("spost() is not defined for this type of operator: $(typeof(op))."))
+
+"""
+    sprepost(op)
+
+Create a super-operator equivalent for left and right side operator multiplication.
+
+For operators ``A``, ``B``, ``C`` the relation
+
+```math
+    \\mathrm{sprepost}(A, B) C = A C B
+```
+
+holds. `A` ond `B` can be dense or a sparse operators.
+"""
+sprepost(A::AbstractOperator, B::AbstractOperator) = throw(ArgumentError("sprepost() is not defined for these types of operator: $(typeof(A)) and  $(typeof(B))."))
+
+function _check_input(H::AbstractOperator, J::Vector, Jdagger::Vector, rates)
+    for j=J
+        @assert isa(j, AbstractOperator)
+    end
+    for j=Jdagger
+        @assert isa(j, AbstractOperator)
+    end
+    @assert length(J)==length(Jdagger)
+    if isa(rates, Matrix{<:Number})
+        @assert size(rates, 1) == size(rates, 2) == length(J)
+    elseif isa(rates, Vector{<:Number})
+        @assert length(rates) == length(J)
+    end
+end
+
+"""
+    liouvillian(H, J; rates, Jdagger)
+
+Create a super-operator equivalent to the master equation so that ``\\dot ρ = S ρ``.
+
+The super-operator ``S`` is defined by
+
+```math
+S ρ = -\\frac{i}{ħ} [H, ρ] + \\sum_i J_i ρ J_i^† - \\frac{1}{2} J_i^† J_i ρ - \\frac{1}{2} ρ J_i^† J_i
+```
+
+# Arguments
+* `H`: Hamiltonian.
+* `J`: Vector containing the jump operators.
+* `rates`: Vector or matrix specifying the coefficients for the jump operators.
+* `Jdagger`: Vector containing the hermitian conjugates of the jump operators. If they
+             are not given they are calculated automatically.
+"""
+function liouvillian(H, J; rates=ones(length(J)), Jdagger=dagger.(J))
+    _check_input(H, J, Jdagger, rates)
+    L = spre(-1im*H) + spost(1im*H)
+    if isa(rates, AbstractMatrix)
+        for i=1:length(J), j=1:length(J)
+            jdagger_j = rates[i,j]/2*Jdagger[j]*J[i]
+            L -= spre(jdagger_j) + spost(jdagger_j)
+            L += spre(rates[i,j]*J[i]) * spost(Jdagger[j])
+        end
+    elseif isa(rates, AbstractVector)
+        for i=1:length(J)
+            jdagger_j = rates[i]/2*Jdagger[i]*J[i]
+            L -= spre(jdagger_j) + spost(jdagger_j)
+            L += spre(rates[i]*J[i]) * spost(Jdagger[i])
+        end
+    end
+    return L
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0bccf25..826fe33 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -26,6 +26,7 @@ end
 println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREADS = $(Sys.CPU_THREADS)`...")
 
 @doset "sortedindices"
+@doset "bases"
 #VERSION >= v"1.9" && @doset "doctests"
 get(ENV,"JET_TEST","")=="true" && @doset "jet"
 VERSION >= v"1.9" && @doset "aqua"
diff --git a/test/test_bases.jl b/test/test_bases.jl
new file mode 100644
index 0000000..ecdeac3
--- /dev/null
+++ b/test/test_bases.jl
@@ -0,0 +1,55 @@
+using Test
+using QuantumInterface: tensor, ⊗, ptrace, reduced, permutesystems, multiplicable
+using QuantumInterface: GenericBasis, CompositeBasis, NLevelBasis, FockBasis
+
+@testset "basis" begin
+
+d1 = 5
+d2 = 2
+d3 = 6
+
+b1 = GenericBasis(d1)
+b2 = GenericBasis(d2)
+b3 = GenericBasis(d3)
+
+@test length(b1) == d1
+@test length(b2) == d2
+@test b1 != b2
+@test b1 != FockBasis(2)
+@test b1 == b1
+
+@test tensor(b1) == b1
+comp_b1 = tensor(b1, b2)
+comp_uni = b1 ⊗ b2
+comp_b2 = tensor(b1, b1, b2)
+@test size(comp_b1) == [d1, d2]
+@test size(comp_uni) == [d1, d2]
+@test size(comp_b2) == [d1, d1, d2]
+
+@test b1^3 == CompositeBasis(b1, b1, b1)
+@test (b1⊗b2)^2 == CompositeBasis(b1, b2, b1, b2)
+@test_throws ArgumentError b1^(0)
+
+comp_b1_b2 = tensor(comp_b1, comp_b2)
+@test size(comp_b1_b2) == [d1, d2, d1, d1, d2]
+@test comp_b1_b2 == CompositeBasis(b1, b2, b1, b1, b2)
+
+@test_throws ArgumentError tensor()
+@test size(comp_b2) == size(tensor(b1, comp_b1))
+@test comp_b2 == tensor(b1, comp_b1)
+
+@test_throws ArgumentError ptrace(comp_b1, [1, 2])
+@test ptrace(comp_b2, [1]) == ptrace(comp_b2, [2]) == comp_b1 == ptrace(comp_b2, 1)
+@test ptrace(comp_b2, [1, 2]) == ptrace(comp_b1, [1])
+@test ptrace(comp_b2, [2, 3]) == ptrace(comp_b1, [2])
+@test ptrace(comp_b2, [2, 3]) == reduced(comp_b2, [1])
+@test_throws ArgumentError reduced(comp_b1, [])
+
+comp1 = tensor(b1, b2, b3)
+comp2 = tensor(b2, b1, b3)
+@test permutesystems(comp1, [2,1,3]) == comp2
+
+@test [b1, b2] != [b1, b3]
+@test comp1 != b1 ⊗ b2 ⊗ NLevelBasis(prod(size(b3)))
+
+end # testset
diff --git a/test/test_sortedindices.jl b/test/test_sortedindices.jl
index ecdf308..d44b9be 100644
--- a/test/test_sortedindices.jl
+++ b/test/test_sortedindices.jl
@@ -18,15 +18,15 @@ x = [3, 5]
 s.reducedindices!(x, [2, 3, 5, 6])
 @test x == [2, 3]
 
-@test_throws AssertionError s.check_indices(5, [1, 6])
-@test_throws AssertionError s.check_indices(5, [0, 2])
+@test_throws ArgumentError s.check_indices(5, [1, 6])
+@test_throws ArgumentError s.check_indices(5, [0, 2])
 @test s.check_indices(5, Int[]) == nothing
 @test s.check_indices(5, [1, 3]) == nothing
 @test s.check_indices(5, [3, 1]) == nothing
 
-@test_throws AssertionError s.check_sortedindices(5, [1, 6])
-@test_throws AssertionError s.check_sortedindices(5, [3, 1])
-@test_throws AssertionError s.check_sortedindices(5, [0, 2])
+@test_throws ArgumentError s.check_sortedindices(5, [1, 6])
+@test_throws ArgumentError s.check_sortedindices(5, [3, 1])
+@test_throws ArgumentError s.check_sortedindices(5, [0, 2])
 @test s.check_sortedindices(5, Int[]) == nothing
 @test s.check_sortedindices(5, [1, 3]) == nothing