diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9dc80a2..59a6711 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,17 @@
 # News
 
+## v0.4.0 - 2025-04-22
+
+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-21
 
 - Deprecate `PauliBasis` as it is identical to `SpinBasis(1//2)^N` but without the compatibility with `CompositeBasis`.
diff --git a/Project.toml b/Project.toml
index 4b1f053..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.10"
+version = "0.4.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/QuantumInterface.jl b/src/QuantumInterface.jl
index 043add9..2a23f16 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")
diff --git a/src/abstract_types.jl b/src/abstract_types.jl
index f8667c9..96c77a0 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)
+    print(stream, "$(typeof(x).name.name)(dim=$(dimension(x.basis_l))x$(dimension(x.basis_r)))\n")
+    if multiplicable(x,x)
         print(stream, "  basis: ")
         show(stream, basis(x))
     else
diff --git a/src/bases.jl b/src/bases.jl
index 9969c1a..b031c3e 100644
--- a/src/bases.jl
+++ b/src/bases.jl
@@ -1,35 +1,93 @@
 """
-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 `dimension`, 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).
+Custom subtypes can also define composite systems by implementing
+`Base.length` and `Base.getindex`.
+
+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
 
+abstract type OperatorBasis <: 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)
 
+Return the number of subsystems of a quantum object in its tensor product
+decomposition.
+
+See also [`CompositeBasis`](@ref).
+"""
+Base.length(b::Basis) = 1
+
+"""
+    getindex(b::Basis)
+
+Get the i'th factor in the tensor product decomposition of the basis into
+subsystems.
+
+See also [`CompositeBasis`](@ref).
+"""
+Base.getindex(b::Basis, i) = i==1 ? b : throw(BoundsError(b,i))
+Base.firstindex(b::Basis) = 1
+Base.lastindex(b::Basis) = length(b)
+
+Base.iterate(b::Basis, state=1) = state > length(b) ? nothing : (b[state], state+1)
+
+"""
+    dimension(b::Basis)
+
 Total dimension of the Hilbert space.
 """
-Base.length(b::Basis) = prod(b.shape)
+dimension(b::Basis) = throw(ArgumentError("dimesion() is not defined for $(typeof(b))"))
 
 """
-    basis(a)
+    shape(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 [`CompositeBasis`](@ref).
 """
-function basis end
+shape(b::Basis) = [dimension(b[i]) for i=1:length(b)]
 
+##
+# GenericBasis, CompositeBasis, SumBasis
+##
 
 """
     GenericBasis(N)
@@ -41,33 +99,50 @@ 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
+dimension(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.
-"""
-struct CompositeBasis{S,B} <: Basis
-    shape::S
-    bases::B
+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{B<:Basis} <: Basis
+    bases::Vector{B}
+    shape::Vector{Int}
+    lengths::Vector{Int}
+    N::Int
+    D::Int
+    function CompositeBasis(bases::Vector{B}) where B<:Basis
+        # to enable this check the the lazy operators in QuantumOpticsBase need to be changed
+        #length(bases) > 1 || throw(ArgumentError("CompositeBasis must only be used for composite systems"))
+        shape_ = mapreduce(shape, vcat, bases)
+        lengths = cumsum(map(length, bases))
+        new{B}(bases, shape_, lengths, lengths[end], prod(shape_))
+    end
 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)
-
-tensor(b::Basis) = b
+CompositeBasis(bases::Basis...) = CompositeBasis([bases...])
+CompositeBasis(bases::Tuple) = CompositeBasis([bases...])
+
+Base.:(==)(b1::CompositeBasis, b2::CompositeBasis) = all(((i, j),) -> i == j, zip(b1.bases, b2.bases))
+Base.length(b::CompositeBasis) = b.N
+# should probably factor this out and think about "proper" compression...
+# also add a labeled basis as a "layer" on top in linear indices...
+function Base.getindex(b::CompositeBasis, i::Integer)
+    (i < 1 || i > b.N) && throw(BoundsError(b,i))
+    bases_idx = findfirst(l -> i<=l, b.lengths) 
+    inner_idx = i - (bases_idx == 1 ? 0 : b.lengths[bases_idx-1])
+    b.bases[bases_idx][inner_idx]
+end
+Base.getindex(b::CompositeBasis, indices) = [b[i] for i in indices]
+shape(b::CompositeBasis) = b.shape
+dimension(b::CompositeBasis) = b.D
 
 """
     tensor(x::Basis, y::Basis, z::Basis...)
@@ -77,43 +152,93 @@ 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...))
+tensor(b1::Basis, b2::Basis) = CompositeBasis([b1, b2])
+tensor(bases::Basis...) = reduce(tensor, bases)
+tensor(basis::Basis) = basis
+
+function tensor(b1::CompositeBasis, b2::CompositeBasis)
+    if typeof(b1.bases[end]) == typeof(b2.bases[1])
+        t = tensor(b1.bases[end], b2.bases[1])
+        if !(t isa CompositeBasis)
+            return CompositeBasis([b1.bases[1:end-1]; t;  b2.bases[2:end]])
+        end
+    end
+    return CompositeBasis([b1.bases; b2.bases])
+end
+
 function tensor(b1::CompositeBasis, b2::Basis)
-    N = length(b1.bases)
-    shape = vcat(b1.shape, length(b2))
-    bases = (b1.bases..., b2)
-    CompositeBasis(shape, bases)
+    if b1.bases[end] isa typeof(b2)
+        t = tensor(b1.bases[end], b2)
+        if !(t isa CompositeBasis)
+            return CompositeBasis([b1.bases[1:end-1]; t])
+        end
+    end
+    return CompositeBasis([b1.bases; b2])
 end
+
 function tensor(b1::Basis, b2::CompositeBasis)
-    N = length(b2.bases)
-    shape = vcat(length(b1), b2.shape)
-    bases = (b1, b2.bases...)
-    CompositeBasis(shape, bases)
+    if b2.bases[1] isa typeof(b1)
+        t = tensor(b1, b2.bases[1])
+        if !(t isa CompositeBasis)
+            return CompositeBasis([t; b2[2:end]])
+        end
+    end
+    return CompositeBasis([b1; b2.bases])
 end
-tensor(bases::Basis...) = reduce(tensor, bases)
 
 Base.:^(b::Basis, N::Integer) = tensor_pow(b, N)
 
 """
-    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([dimension(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))
+dimension(b::SumBasis) = sum(b.shape)
+
+
+
+"""
+    nsubspaces(b)
+
+Return the number of subspaces of a [`SumBasis`](@ref) in its direct sum
+decomposition.
+"""
+nsubspaces(b::SumBasis) = length(b.bases)
+
+"""
+    subspace(b, i)
+
+Return the basis for the `i`th subspace of of a [`SumBasis`](@ref).
+"""
+subspace(b::SumBasis, i) = b.bases[i]
+
+"""
+    directsum(b1::Basis, b2::Basis)
+
+Construct the [`SumBasis`](@ref) out of two sub-bases.
+"""
+directsum(b1::Basis, b2::Basis) = SumBasis([dimension(b1), dimension(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; dimension(b2)], [b1.bases; b2])
+directsum(b1::Basis, b2::SumBasis) = SumBasis([dimension(b1); b2.shape], [b1; b2.bases])
+directsum(bases::Basis...) = reduce(directsum, bases)
+directsum(basis::Basis) = basis
+
+# TODO: what to do about embed for SumBasis?
+#embed(b::SumBasis, indices, ops) = embed(b, b, indices, ops)
+
+##
+# Basis checks
+##
 
 """
 Exception that should be raised for an illegal algebraic operation.
@@ -123,33 +248,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)
@@ -157,34 +282,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
@@ -197,13 +316,13 @@ Reduced basis, state or operator on the specified subsystems.
 The `indices` argument, which can be a single integer or a vector of integers,
 specifies which subsystems are kept. At least one index must be specified.
 """
-function reduced(b::CompositeBasis, indices)
+function reduced(b::Basis, indices)
     if length(indices)==0
         throw(ArgumentError("At least one subsystem must be specified in reduced."))
     elseif length(indices)==1
-        return b.bases[indices[1]]
+        return b[indices[1]]
     else
-        return CompositeBasis(b.shape[indices], b.bases[indices])
+        return tensor(b[indices]...)
     end
 end
 
@@ -217,12 +336,14 @@ specifies which subsystems are traced out. The number of indices has to be
 smaller than the number of subsystems, i.e. it is not allowed to perform a
 full trace.
 """
-function ptrace(b::CompositeBasis, indices)
-    J = [i for i in 1:length(b.bases) if i ∉ indices]
+function ptrace(b::Basis, indices)
+    J = [i for i in 1:length(b) if i ∉ indices]
     length(J) > 0 || throw(ArgumentError("Tracing over all indices is not allowed in ptrace."))
     reduced(b, J)
 end
 
+_index_complement(b::Basis, indices) = complement(length(b), indices)
+reduced(a, indices) = ptrace(a, _index_complement(basis(a), indices))
 
 """
     permutesystems(a, perm)
@@ -232,10 +353,10 @@ Change the ordering of the subsystems of the given object.
 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)
-    (nsubsystems(b) == length(perm)) || throw(ArgumentError("Must have nsubsystems(b) == length(perm) in permutesystems"))
+function permutesystems(b::Basis, perm)
+    (length(b) == length(perm)) || throw(ArgumentError("Must have length(b) == length(perm) in permutesystems"))
     isperm(perm) || throw(ArgumentError("Must pass actual permeutation to permutesystems"))
-    CompositeBasis(b.shape[perm], b.bases[perm])
+    tensor(b.bases[perm]...)
 end
 
 
@@ -247,22 +368,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(ArgumentError("Fock cutoff and offset must be positive and cutoff must be less than offset"))
         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)
+dimension(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
 
 
 """
@@ -270,100 +411,173 @@ 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
         N > 0 || throw(ArgumentError("N must be greater than 0"))
-        new{T}([N], N)
+        new{T}(N)
     end
 end
 
 Base.:(==)(b1::NLevelBasis, b2::NLevelBasis) = b1.N == b2.N
-
+dimension(b::NLevelBasis) = b.N
 
 """
-    SpinBasis(n)
+    SpinBasis(n, N=1)
 
-Basis for spin-n particles.
+Basis for spin-`n` particles over `N` systems.
 
-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}
+    D::T
+    N::T
+    function SpinBasis(spinnumber::Rational{T}, N=1) where T
         n = numerator(spinnumber)
         d = denominator(spinnumber)
         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)
+        D = numerator(spinnumber*2 + 1)
+        new{T}(spinnumber, D, N)
     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.:(==)(b1::SpinBasis, b2::SpinBasis) = b1.D==b2.D && b1.N == b2.N
+Base.length(b::SpinBasis) = b.N
+Base.getindex(b::SpinBasis, i) = SpinBasis(b.spinnumber, length(i))
+shape(b::SpinBasis) = fill(b.D, b.N)
+dimension(b::SpinBasis) = b.D^b.N
+function tensor(b1::SpinBasis, b2::SpinBasis)
+    if b1.spinnumber == b2.spinnumber
+        return SpinBasis(b1.spinnumber, b1.N+b2.N)
+    else
+        return CompositeBasis([b1, b2])
+    end
+end
+
+"""
+    spinnumber(b::SpinBasis)
 
+Return the spin number of the given spin basis.
 
+See [`SpinBasis`](@ref).
 """
-    SumBasis(b1, b2...)
+spinnumber(b::SpinBasis) = b.spinnumber
 
-Similar to [`CompositeBasis`](@ref) but for the [`directsum`](@ref) (⊕)
+
+##
+# Operator Bases
+##
+
+"""
+    KetBraBasis(BL,BR)
+
+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.
 """
-struct SumBasis{S,B} <: Basis
-    shape::S
-    bases::B
+struct KetBraBasis{BL<:Basis, BR<:Basis} <: OperatorBasis
+    left::BL
+    right::BR
 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...,))
+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)
+dimension(b::KetBraBasis) = dimension(b.left)*dimension(b.right)
+tensor(b1::KetBraBasis, b2::KetBraBasis) = KetBraBasis(tensor(b1.left,b2.left), tensor(b1.right, b2.right))
 
-==(b1::T, b2::T) where T<:SumBasis = equal_shape(b1.shape, b2.shape)
-==(b1::SumBasis, b2::SumBasis) = false
-length(b::SumBasis) = sum(b.shape)
+"""
+    ChoiBasis(ref_basis,out_basis)
 
+The Choi basis is used to represent superoperators in the Choi representation
+where the `ref_basis` denotes the ancillary reference system with which an input
+state will be jointly measured in order to accomplish teleportation simulation
+of the channel with the channel's output appearing in the `out_basis` system.
 """
-    directsum(b1::Basis, b2::Basis)
+struct ChoiBasis{BL<:Basis, BR<:Basis} <: OperatorBasis
+    ref::BL
+    out::BR
+end
+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)
+dimension(b::ChoiBasis) = dimension(b.ref)*dimension(b.out)
+tensor(b1::ChoiBasis, b2::ChoiBasis) = ChoiBasis(tensor(b1.ref,b2.ref), tensor(b1.out, b2.out))
 
-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...,))
+    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} <: OperatorBasis
+    N::T
 end
-function directsum(b1::Basis, b2::SumBasis)
-    shape = [length(b1);b2.shape]
-    bases = [b1;b2.bases...]
-    return SumBasis(shape, (bases...,))
+basis_l(b::PauliBasis) = SpinBasis(1//2)^b.N
+basis_r(b::PauliBasis) = SpinBasis(1//2)^b.N
+Base.:(==)(b1::PauliBasis, b2::PauliBasis) = b1.N == b2.N
+dimension(b::PauliBasis) = 4^b.N
+tensor(b1::PauliBasis, b2::PauliBasis) = PauliBasis(b1.N+b2.N)
+
+"""
+    ChiBasis(N)
+
+The basis for a Chi process matrix, which is just the Choi state in the Pauli
+operator basis. However we do not use the `ChoiBasis`, partly to have easier
+dispatch on types, and partly because there's no sensible way to distingish
+between the "reference" and "output" systems as that information is lost in the
+computational to Pauli basis transformation (i.e. two indices into one).
+
+TODO explain better why dimension base is 2, see sec III.E.
+"""
+struct ChiBasis{T<:Integer} <: OperatorBasis
+    Nl::T
+    Nr::T
 end
-function directsum(b1::SumBasis, b2::SumBasis)
-    shape = [b1.shape;b2.shape]
-    bases = [b1.bases...;b2.bases...]
-    return SumBasis(shape, (bases...,))
+basis_l(b::ChiBasis) = SpinBasis(1//2)^b.Nl
+basis_r(b::ChiBasis) = SpinBasis(1//2)^b.Nr
+Base.:(==)(b1::ChiBasis, b2::ChiBasis) = (b1.Nl == b2.Nl && b1.Nr == b2.Nr)
+dimension(b::ChiBasis) = 2^(b.Nl+b.Nr)
+tensor(b1::ChiBasis, b2::ChiBasis) = ChiBasis(b1.Nl+b2.Nl, b1.Nr+b2.Nr)
+
+"""
+    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} <: OperatorBasis
+    shape::Vector{T}
 end
+HWPauliBasis(N::Integer) = HWPauliBasis([N])
+basis_l(b::HWPauliBasis) = tensor(NLevelBasis.(b.shape)...)
+basis_r(b::HWPauliBasis) = tensor(NLevelBasis.(b.shape)...)
+Base.:(==)(b1::HWPauliBasis, b2::HWPauliBasis) = b1.shape == b2.shape
+dimension(b::HWPauliBasis) = prod([n^2 for n in b.shape])
+tensor(b1::HWPauliBasis, b2::HWPauliBasis) = HWPauliBasis([b1.shape; b2.shape])
 
-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)
@@ -385,6 +599,9 @@ function show(stream::IO, x::SpinBasis)
     else
         write(stream, "Spin($n/$d)")
     end
+    if x.N > 1
+        write(stream, "^$(x.N)")
+    end
 end
 
 function show(stream::IO, x::FockBasis)
@@ -409,3 +626,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
index ad26bea..70da1af 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -11,16 +11,80 @@ function equal_bases(a, b)
     return true
 end
 
-struct PauliBasis{S,B} <: Basis
-    shape::S
-    bases::B
-    function PauliBasis(num_qubits::T) where {T<:Integer}
-        Base.depwarn("`PauliBasis` will be removed in future versions as it does not represent the Pauli operator basis. SpinBasis(1//2)^N or NLevelBasis(2)^N should likely be used instead.", :PauliBasis)
-        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)
+# 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
 
-Base.:(==)(pb1::PauliBasis, pb2::PauliBasis) = length(pb1.bases) == length(pb2.bases)
+# 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 9945def..96abeef 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::Basis, br::Basis,
                operators::Dict{<:Vector{<:Integer}, T}) where T<:AbstractOperator
-    (nsubsystems(basis_l) == nsubsystems(basis_r)) || throw(ArgumentError("Must have nsubsystems(bl) == nsubsystems(br) in embed"))
-    N = length(basis_l.bases)::Int # type assertion to help type inference
+    (length(bl) == length(br)) || throw(ArgumentError("Must have length(bl) == length(br) in embed"))
+    N = length(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,29 +22,29 @@ 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)
     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...)
+embed(bl::Basis, br::Basis, operators::Dict{<:Integer, T}; kwargs...) where {T<:AbstractOperator} = embed(bl, br, Dict([i]=>op_i for (i, op_i) in operators); kwargs...)
+embed(b::Basis, operators::Dict{<:Integer, T}; kwargs...) where {T<:AbstractOperator} = embed(b, b, operators; kwargs...)
+embed(b::Basis, operators::Dict{<:Vector{<:Integer}, T}; kwargs...) where {T<:AbstractOperator} = embed(b, b, operators; kwargs...)
 
 # 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(bl::Basis, br::Basis, indices, op::T) where T<:AbstractOperator = embed(bl, br, 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)
+embed(bl::Basis, br::Basis, index::Integer, op::AbstractOperator) = embed(bl, br, index, [op])
+embed(b::Basis, indices, operators::Vector{T}) where {T<:AbstractOperator} = embed(b, b, indices, operators)
+embed(b::Basis, indices, op::AbstractOperator) = embed(b, b, indices, op)
 
 
 """
@@ -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::Basis, br::Basis,
                indices, operators::Vector{T}) where T<:AbstractOperator
 
     check_embed_indices(indices) || throw(ArgumentError("Must have unique indices in embed"))
-    (nsubsystems(basis_l) == nsubsystems(basis_r)) || throw(ArgumentError("Must have nsubsystems(bl) == nsubsystems(br) in embed"))
+    (length(bl) == length(br)) || throw(ArgumentError("Must have length(bl) == length(br) in embed"))
     (length(indices) == length(operators)) || throw(ArgumentError("Must have length(indices) == length(operators) in embed"))
 
-    N = nsubsystems(basis_l)
+    N = length(bl)
 
     # Embed all single-subspace operators.
     idxop_sb = [x for x in zip(indices, operators) if x[1] isa Integer]
@@ -67,26 +67,20 @@ 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()) # FIXME issue #12
-        (opsb.basis_r == basis_r.bases[idxsb]) || throw(IncompatibleBases()) # FIXME issue #12
+        (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
 end
 
 permutesystems(a::AbstractOperator, perm) = arithmetic_unary_error("Permutations of subsystems", a)
-
-nsubsystems(s::AbstractKet) = 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
diff --git a/src/expect_variance.jl b/src/expect_variance.jl
index 86f5e64..77f5aa2 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(length(basis(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(length(basis(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 47b0056..6ff5096 100644
--- a/src/julia_base.jl
+++ b/src/julia_base.jl
@@ -8,11 +8,13 @@ addnumbererror() = throw(ArgumentError("Can't add or subtract a number and an op
 # States
 ##
 
--(a::T) where {T<:StateVector} = T(a.basis, -a.data) # FIXME issue #12
+==(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)) # FIXME issue #12
-length(a::StateVector) = length(a.basis)::Int # FIXME issue #12
-basis(a::StateVector) = a.basis # FIXME issue #12
+copy(a::T) where {T<:StateVector} = T(basis(a), copy(a.data)) # FIXME issue #12
+length(a::StateVector) = dimension(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
@@ -26,15 +28,18 @@ Base.eltype(x::StateVector) = eltype(x.data) # FIXME issue #12
 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  # FIXME issue #12
-basis(a::AbstractOperator) = (check_samebases(a); a.basis_l) # FIXME issue #12
-basis(a::AbstractSuperOperator) = (check_samebases(a); a.basis_l[1]) # FIXME issue #12
+length(a::AbstractOperator) = dimension(basis_l(a))::Int*dimension(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)) # FIXME issue #12
+Base.size(op::AbstractOperator) = (dimension(basis_l(op)),dimension(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) # FIXME issue #12
+    i==1 ? dimension(basis_l(op)) : dimension(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 2cbbcc2..1926883 100644
--- a/src/julia_linalg.jl
+++ b/src/julia_linalg.jl
@@ -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 71833e3..0000000
--- a/src/linalg.jl
+++ /dev/null
@@ -1,10 +0,0 @@
-samebases(a::AbstractOperator) = samebases(a.basis_l, a.basis_r)::Bool # FIXME issue #12
-samebases(a::AbstractOperator, b::AbstractOperator) = samebases(a.basis_l, b.basis_l)::Bool && samebases(a.basis_r, b.basis_r)::Bool # FIXME issue #12
-check_samebases(a::Union{AbstractOperator, AbstractSuperOperator}) = check_samebases(a.basis_l, a.basis_r) # FIXME issue #12
-multiplicable(a::AbstractOperator, b::AbstractOperator) = multiplicable(a.basis_r, b.basis_l) # FIXME issue #12
-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/superoperators.jl b/src/superoperators.jl
new file mode 100644
index 0000000..7f69ad0
--- /dev/null
+++ b/src/superoperators.jl
@@ -0,0 +1,163 @@
+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::AbstractKet) = 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))."))
+kraus(op::Vector{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.
+"""
+function spre(op::AbstractOperator)
+    multiplicable(op, op) || throw(ArgumentError("It's not clear what spre of a non-square operator should be. See issue #113"))
+    sprepost(op, identityoperator(op))
+end
+
+"""
+    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.
+"""
+function spost(op::AbstractOperator)
+    multiplicable(op, op) || throw(ArgumentError("It's not clear what spost of a non-square operator should be. See issue #113"))
+    sprepost(identityoperator(op), op)
+end
+
+"""
+    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/test_bases.jl b/test/test_bases.jl
index c480d41..c5545d9 100644
--- a/test/test_bases.jl
+++ b/test/test_bases.jl
@@ -1,19 +1,19 @@
 using Test
-using QuantumInterface: tensor, ⊗, ptrace, reduced, permutesystems, multiplicable
+using QuantumInterface: dimension, shape, tensor, ⊗, ptrace, reduced, permutesystems, multiplicable
 using QuantumInterface: GenericBasis, CompositeBasis, NLevelBasis, FockBasis
 
 @testset "basis" begin
 
-shape1 = [5]
-shape2 = [2, 3]
-shape3 = [6]
+d1 = 5
+d2 = 2
+d3 = 6
 
-b1 = GenericBasis(shape1)
-b2 = GenericBasis(shape2)
-b3 = GenericBasis(shape3)
+b1 = GenericBasis(d1)
+b2 = GenericBasis(d2)
+b3 = GenericBasis(d3)
 
-@test b1.shape == shape1
-@test b2.shape == shape2
+@test dimension(b1) == d1
+@test dimension(b2) == d2
 @test b1 != b2
 @test b1 != FockBasis(2)
 @test b1 == b1
@@ -22,20 +22,20 @@ b3 = GenericBasis(shape3)
 comp_b1 = tensor(b1, b2)
 comp_uni = b1 ⊗ b2
 comp_b2 = tensor(b1, b1, b2)
-@test comp_b1.shape == [prod(shape1), prod(shape2)]
-@test comp_uni.shape == [prod(shape1), prod(shape2)]
-@test comp_b2.shape == [prod(shape1), prod(shape1), prod(shape2)]
+@test shape(comp_b1) == [d1, d2]
+@test shape(comp_uni) == [d1, d2]
+@test shape(comp_b2) == [d1, d1, d2]
 
 @test b1^3 == CompositeBasis(b1, b1, b1)
 @test (b1⊗b2)^2 == CompositeBasis(b1, b2, b1, b2)
 @test_throws DomainError b1^(0)
 
 comp_b1_b2 = tensor(comp_b1, comp_b2)
-@test comp_b1_b2.shape == [prod(shape1), prod(shape2), prod(shape1), prod(shape1), prod(shape2)]
+@test shape(comp_b1_b2) == [d1, d2, d1, d1, d2]
 @test comp_b1_b2 == CompositeBasis(b1, b2, b1, b1, b2)
 
 @test_throws ArgumentError tensor()
-@test comp_b2.shape == tensor(b1, comp_b1).shape
+@test shape(comp_b2) == shape(tensor(b1, comp_b1))
 @test comp_b2 == tensor(b1, comp_b1)
 
 @test_throws ArgumentError ptrace(comp_b1, [1, 2])
@@ -50,6 +50,6 @@ comp2 = tensor(b2, b1, b3)
 @test permutesystems(comp1, [2,1,3]) == comp2
 
 @test [b1, b2] != [b1, b3]
-@test !multiplicable(comp1, b1 ⊗ b2 ⊗ NLevelBasis(prod(b3.shape)))
+@test comp1 != b1 ⊗ b2 ⊗ NLevelBasis(prod(shape(b3)))
 
 end # testset