diff --git a/README.md b/README.md
index e10d612..026c060 100644
--- a/README.md
+++ b/README.md
@@ -4,7 +4,7 @@
 [![DOI](https://zenodo.org/badge/23916883.svg)](https://zenodo.org/badge/latestdoi/23916883)
 # QuantumInformation
 
-A Julia package for numerical computation in quantum information theory.
+A Julia package for numerical computation in quantum information theory. [Published in PLoS ONE](https://doi.org/10.1371/journal.pone.0209358).
 
 Numerical investigations are prevalent in quantum information theory. Numerical experiments can be used to find counter examples for theorems, to test hypotheses or to gain insight about quantum objects and operations.
 
diff --git a/src/base.jl b/src/base.jl
index 5e76fc8..ba112e1 100644
--- a/src/base.jl
+++ b/src/base.jl
@@ -1,11 +1,19 @@
-function ket(::Type{T}, val::Int, dim::Int) where T<:Number
+export ket, bra, ketbra, proj, res, unres, max_mixed, max_entangled,
+    werner_state, permutesystems
+
+function ket(::Type{T}, val::Int, dim::Int) where T<:AbstractVector{<:Number}
     dim > 0 ? () : throw(ArgumentError("Vector dimension has to be nonnegative"))
     1 <= val <= dim ? () : throw(ArgumentError("Label have to be smaller than vector dimension"))
-    ψ = zeros(T, dim)
-    ψ[val] = one(T)
+    ψ = zero(T(undef, dim))
+    ψ[val] = one(eltype(T))
     ψ
 end
 
+function ket(::Type{T}, val::Int, dim::Int) where T<:Number
+    @warn "This method is deprecated and will be removed. Use calls like `ket(Matrix{ComplexF64}, 1, 2)`."
+    ket(Vector{T}, val, dim)
+end
+
 """
 $(SIGNATURES)
 - `val`: non-zero entry - label.
@@ -13,9 +21,14 @@ $(SIGNATURES)
 
 Return complex column vector \$|val\\rangle\$ of unit norm describing quantum state.
 """
-ket(val::Int, dim::Int) = ket(ComplexF64, val, dim)
+ket(val::Int, dim::Int) = ket(Vector{ComplexF64}, val, dim)
 
-bra(::Type{T}, val::Int, dim::Int) where T<:Number = ket(T, val, dim)'
+function bra(::Type{T}, val::Int, dim::Int) where T<:Number 
+    @warn "This method is deprecated and will be removed. Use calls like `bra(Matrix{ComplexF64}, 1, 2)`."
+    ket(Vector{T}, val, dim)'
+end
+
+bra(::Type{T}, val::Int, dim::Int) where T<:AbstractVector{<:Number} = ket(T, val, dim)'
 
 """
 $(SIGNATURES)
@@ -24,16 +37,21 @@ $(SIGNATURES)
 
 Return Hermitian conjugate \$\\langle val| = |val\\rangle^\\dagger\$ of the ket with the same label.
 """
-bra(val::Int, dim::Int) = bra(ComplexF64, val, dim)
+bra(val::Int, dim::Int) = bra(Vector{ComplexF64}, val, dim)
 
-function ketbra(::Type{T}, valk::Int, valb::Int, dim::Int) where T<:Number
+function ketbra(::Type{T}, valk::Int, valb::Int, dim::Int) where T<:AbstractMatrix{<:Number}
     dim > 0 ? () : throw(ArgumentError("Vector dimension has to be nonnegative"))
     1 <= valk <= dim && 1 <= valb <= dim ? () : throw(ArgumentError("Ket and bra labels have to be smaller than operator dimension"))
-    ρ = zeros(T, dim, dim)
-    ρ[valk,valb] = one(T)
+    ρ = zero(T(undef, dim, dim))
+    ρ[valk,valb] = one(eltype(T))
     ρ
 end
 
+function ketbra(::Type{T}, valk::Int, valb::Int, dim::Int) where T<:Number
+    @warn "This method is deprecated and will be removed. Use calls like `ketbra(Matrix{ComplexF64}, 1, 1, 2)`."
+    ketbra(Matrix{T}, valk, valb, dim)
+end
+
 """
 $(SIGNATURES)
 - `valk`: non-zero entry - label.
@@ -42,7 +60,7 @@ $(SIGNATURES)
 
 # Return outer product \$|valk\\rangle\\langle vakb|\$ of states \$|valk\\rangle\$ and \$|valb\\rangle\$.
 """
-ketbra(valk::Int, valb::Int, dim::Int) = ketbra(ComplexF64, valk, valb, dim)
+ketbra(valk::Int, valb::Int, dim::Int) = ketbra(Matrix{ComplexF64}, valk, valb, dim)
 
 """
 $(SIGNATURES)
diff --git a/src/channels.jl b/src/channels.jl
index c0a973a..2923713 100644
--- a/src/channels.jl
+++ b/src/channels.jl
@@ -61,7 +61,7 @@ struct SuperOperator{T<:AbstractMatrix{<:Number}} <: AbstractQuantumOperation{T}
             throw(ArgumentError("Superoperator matrix has invalid dimensions"))
         end
         odim, idim = sr, sc
-        new{T1}(T1(m), idim, odim)
+        new{T1}(convert(T1, m), idim, odim)
     end
 end
 
@@ -103,7 +103,7 @@ struct DynamicalMatrix{T<:AbstractMatrix{<:Number}} <: AbstractQuantumOperation{
         if r!=c || r!=idim*odim
             throw(ArgumentError("DynamicalMatrix matrix has invalid dimensions"))
         end
-        new(T1(m), idim, odim)
+        new(convert(T1, m), idim, odim)
     end
 end
 
@@ -139,13 +139,13 @@ struct UnitaryChannel{T<:AbstractMatrix{<:Number}} <: AbstractQuantumOperation{T
     function UnitaryChannel{T1}(m::T2) where {T1<:AbstractMatrix{<:Number}, T2<:AbstractMatrix{<:Number}}
         odim, idim = size(m)
         idim == odim ? () : throw(ArgumentError("UnitaryChannel matrix has to be square"))
-        new{T1}(T1(m), idim, odim)
+        new{T1}(convert(T1,m), idim, odim)
     end
 end
 
 function UnitaryChannel{T1}(m::T2, idim::Int, odim::Int) where {T1<:AbstractMatrix{<:Number}, T2<:AbstractMatrix{<:Number}}
     (odim, idim) == size(m) ? () : throw(ArgumentError("Matrix size and operator dimensions mismatch"))
-    UnitaryChannel{T1}(T1(m))
+    UnitaryChannel{T1}(convert(T1,m))
 end
 
 """
@@ -201,7 +201,7 @@ struct PostSelectionMeasurement{T<:AbstractMatrix{<:Number}} <: AbstractQuantumO
     odim::Int
     function PostSelectionMeasurement{T1}(m::T2) where {T1<:AbstractMatrix{<:Number}, T2<:AbstractMatrix{<:Number}}
         odim, idim = size(m)
-        new{T1}(T1(m), idim, odim)
+        new{T1}(convert(T1,m), idim, odim)
     end
 end
 
@@ -311,7 +311,7 @@ Transforms list of Kraus operators into super-operator matrix.
 """
 function Base.convert(::Type{SuperOperator{T1}}, Φ::KrausOperators{T2}) where {T1<:AbstractMatrix{<:Number}, T2<:AbstractMatrix{<:Number}}
     m = sum(k⊗(conj.(k)) for k in Φ.matrices)
-    SuperOperator{T1}(T1(m), Φ.idim, Φ.odim)
+    SuperOperator{T1}(convert(T1,m), Φ.idim, Φ.odim)
 end
 
 """
@@ -325,7 +325,7 @@ function Base.convert(::Type{Stinespring{T1}}, Φ::KrausOperators{T2}) where {T1
     ko = orthogonalize(Φ)
     # TODO: improvement: transform to stacking
     m = sum(k ⊗ ket(i, ko.idim*ko.odim) for (i, k) in enumerate(ko.matrices))
-    Stinespring{T1}(T1(m), ko.idim, ko.odim)
+    Stinespring{T1}(convert(T1,m), ko.idim, ko.odim)
 end
 
 """
@@ -337,7 +337,7 @@ Transforms list of Kraus operators into dynamical matrix.
 """
 function Base.convert(::Type{DynamicalMatrix{T1}}, Φ::KrausOperators{T2}) where {T1<:AbstractMatrix{<:Number}, T2<:AbstractMatrix{<:Number}}
     m = sum(res(k) * res(k)' for k in Φ.matrices)
-    DynamicalMatrix{T1}(T1(m), Φ.idim, Φ.odim)
+    DynamicalMatrix{T1}(convert(T1,m), Φ.idim, Φ.odim)
 end
 
 """
@@ -360,7 +360,7 @@ Transforms super-operator matrix into dynamical matrix.
 """
 function Base.convert(::Type{DynamicalMatrix{T1}}, Φ::SuperOperator{T2}) where {T1<:AbstractMatrix{<:Number}, T2<:AbstractMatrix{<:Number}}
     m = reshuffle(Φ.matrix, [Φ.odim Φ.odim; Φ.idim Φ.idim])
-    DynamicalMatrix{T1}(T1(m), Φ.idim, Φ.odim)
+    DynamicalMatrix{T1}(convert(T1,m), Φ.idim, Φ.odim)
 end
 
 """
@@ -421,11 +421,11 @@ Transforms dynamical matrix into super-operator matrix.
 """
 function Base.convert(::Type{SuperOperator{T1}}, Φ::DynamicalMatrix{T2}) where {T1<:AbstractMatrix{<:Number}, T2<:AbstractMatrix{<:Number}}
     m = reshuffle(Φ.matrix, [Φ.odim Φ.idim; Φ.odim Φ.idim])
-    SuperOperator{T1}(T1(m), Φ.idim, Φ.odim)
+    SuperOperator{T1}(convert(T1,m), Φ.idim, Φ.odim)
 end
 
 function Base.convert(::Type{KrausOperators{T1}}, Φ::UnitaryChannel{T2}) where {T1<:AbstractMatrix{<:Number}, T2<:AbstractMatrix{<:Number}}
-    KrausOperators{T1}(T1[T1(Φ.matrix)], Φ.idim, Φ.odim)
+    KrausOperators{T1}(T1[convert(T1,Φ.matrix)], Φ.idim, Φ.odim)
 end
 
 function Base.convert(::Type{KrausOperators{T1}}, Φ::IdentityChannel{T2}) where {T1<:AbstractMatrix{N1}, T2<:AbstractMatrix{N2}} where {N1<:Number, N2<:Number}
@@ -442,14 +442,14 @@ function Base.convert(::Type{KrausOperators{T1}}, Φ::POVMMeasurement{T2}) where
     for (i, p) in enumerate(Φ.matrices)
         sqrtp = sqrt(p)
         k = ket(i, Φ.odim)*sum(bra(j, Φ.idim)*sqrtp for j in 1:Φ.idim)
-        push!(v, T1(k))
+        push!(v, convert(T1,k))
     end
     KrausOperators{T1}(v, Φ.idim, Φ.odim)
 end
 
 function Base.convert(::Type{KrausOperators{T1}}, Φ::PostSelectionMeasurement{T2}) where {T1<:AbstractMatrix{<:Number}, T2<:AbstractMatrix{<:Number}}
     m = Φ.matrix
-    v = T1[T1(m)]
+    v = T1[convert(T1,m)]
     KrausOperators{T1}(v, Φ.idim, Φ.odim)
 end
 
diff --git a/src/gates.jl b/src/gates.jl
index d79a81e..6ff4331 100644
--- a/src/gates.jl
+++ b/src/gates.jl
@@ -1,4 +1,3 @@
-
 export 𝕀, sx, sy, sz, qft, hadamard, grover
 
 sx = ComplexF64[0 1; 1 0]
diff --git a/src/matrixbases.jl b/src/matrixbases.jl
index e5676cb..fafa99c 100644
--- a/src/matrixbases.jl
+++ b/src/matrixbases.jl
@@ -1,11 +1,21 @@
+import .Base: length, iterate
+export AbstractMatrixBasisIterator, HermitianBasisIterator, AbstractBasis,
+    AbstractMatrixBasis, HermitianBasis, hermitianbasis, represent, combine
 abstract type AbstractMatrixBasisIterator{T<:AbstractMatrix} end
 struct HermitianBasisIterator{T} <: AbstractMatrixBasisIterator{T} 
     dim::Int
 end
-import .Base: length, iterate
-export hermitianbasis, matrixtocoeffs, coeffstomatrix
 
-hermitianbasis(dim::Int) = hermitianbasis(Matrix{ComplexF64}, dim)
+abstract type AbstractBasis end
+abstract type AbstractMatrixBasis{T} <: AbstractBasis where T<:AbstractMatrix{<:Number} end 
+
+struct HermitianBasis{T} <: AbstractMatrixBasis{T} 
+    iterator::HermitianBasisIterator{T}
+
+    function HermitianBasis{T}(dim::Integer) where T<:AbstractMatrix{<:Number}
+        new(HermitianBasisIterator{T}(dim)) 
+    end
+end
 
 """
 $(SIGNATURES)
@@ -15,6 +25,9 @@ Returns elementary hermitian matrices of dimension `dim` x `dim`.
 """
 hermitianbasis(T::Type{<:AbstractMatrix{<:Number}}, dim::Int) = HermitianBasisIterator{T}(dim)
 
+hermitianbasis(dim::Int) = hermitianbasis(Matrix{ComplexF64}, dim)
+
+
 function iterate(itr::HermitianBasisIterator{T}, state=(1,1)) where T<:AbstractMatrix{<:Number}
     dim = itr.dim
     (a, b) = state
@@ -32,10 +45,10 @@ end
 
 length(itr::HermitianBasisIterator) = itr.dim^2
 
-function matrixtocoeffs(basis::T, m::Matrix{<:Number}) where T<:AbstractMatrixBasisIterator
-    tr.([m] .* basis) 
+function represent(basis::T, m::Matrix{<:Number}) where T<:AbstractMatrixBasis
+    tr.([m] .* basis.iterator) 
 end
 
-function coeffstomatrix(basis::T, v::Vector{<:Number}) where T<:AbstractMatrixBasisIterator
-    sum(basis .* v)
+function combine(basis::T, v::Vector{<:Number}) where T<:AbstractMatrixBasis
+    sum(basis.iterator .* v)
 end
\ No newline at end of file
diff --git a/src/ptrace.jl b/src/ptrace.jl
index 04001f7..ef43f42 100644
--- a/src/ptrace.jl
+++ b/src/ptrace.jl
@@ -1,3 +1,5 @@
+export ptrace
+
 """
 $(SIGNATURES)
 - `ρ`: quantum state.
diff --git a/src/ptranspose.jl b/src/ptranspose.jl
index 3d257ca..ad64b17 100644
--- a/src/ptranspose.jl
+++ b/src/ptranspose.jl
@@ -1,3 +1,4 @@
+export ptranspose
 """
 $(SIGNATURES)
 - `ρ`: quantum state.
diff --git a/src/reshuffle.jl b/src/reshuffle.jl
index 1ed3cd3..4360ec1 100644
--- a/src/reshuffle.jl
+++ b/src/reshuffle.jl
@@ -1,3 +1,5 @@
+export reshuffle
+
 # """
 #   Performs reshuffling of indices of a matrix.
 #   Given multiindexed matrix M_{(m,μ),(n,ν)} it returns
diff --git a/test/base.jl b/test/base.jl
index ccf3a81..206268a 100644
--- a/test/base.jl
+++ b/test/base.jl
@@ -5,8 +5,8 @@
     ψ = ComplexF64[1, 0, 0, 0]
     @test norm(ϕ - ψ) ≈ 0.
 
-    @test typeof(ket(Float64, 1, 4)) == Vector{Float64}
-    @test typeof(ket(ComplexF64, 1, 4)) == Vector{ComplexF64}
+    @test typeof(ket(Vector{Float64}, 1, 4)) == Vector{Float64}
+    @test typeof(ket(Vector{ComplexF64}, 1, 4)) == Vector{ComplexF64}
 
     @test_throws ArgumentError ket(4, 3)
 end
@@ -18,8 +18,8 @@ end
 
     @test_throws ArgumentError bra(4,3)
 
-    @test typeof(bra(Float64, 1, 4)) == LinearAlgebra.Adjoint{Float64,Array{Float64,1}}
-    @test typeof(bra(ComplexF64, 1, 4)) == LinearAlgebra.Adjoint{Complex{Float64},Array{Complex{Float64},1}}
+    @test typeof(bra(Vector{Float64}, 1, 4)) == LinearAlgebra.Adjoint{Float64,Array{Float64,1}}
+    @test typeof(bra(Vector{ComplexF64}, 1, 4)) == LinearAlgebra.Adjoint{Complex{Float64},Array{Complex{Float64},1}}
 end
 
 @testset "ketbra" begin
@@ -30,8 +30,8 @@ end
 
     @test_throws ArgumentError ketbra(4, 4, 3)
 
-    @test typeof(ketbra(Float64, 1, 1, 4)) == Matrix{Float64}
-    @test typeof(ketbra(ComplexF64, 1, 1, 4)) == Matrix{ComplexF64}
+    @test typeof(ketbra(Matrix{Float64}, 1, 1, 4)) == Matrix{Float64}
+    @test typeof(ketbra(Matrix{ComplexF64}, 1, 1, 4)) == Matrix{ComplexF64}
 end
 
 @testset "proj" begin
diff --git a/test/channels.jl b/test/channels.jl
index 8ed3a16..de565c9 100644
--- a/test/channels.jl
+++ b/test/channels.jl
@@ -133,6 +133,9 @@ end
 @testset "UnitaryChannel" begin
     @test_throws ArgumentError UnitaryChannel(ones(4, 5))
     @test_throws ArgumentError UnitaryChannel(ones(4, 4), 4, 5)
+
+    c = UnitaryChannel(Diagonal(ComplexF64[1 -1.0im]))
+    @test c isa UnitaryChannel{<:Diagonal}
 end
 
 @testset "POVMMeasurement" begin
diff --git a/test/matrixbases.jl b/test/matrixbases.jl
index a6b0bd5..fbb4444 100644
--- a/test/matrixbases.jl
+++ b/test/matrixbases.jl
@@ -1,8 +1,21 @@
-@testset "base_matrices" begin
+@testset "MatrixBases" begin
+
+@testset "HermitianBasisIterator" begin
     d = 4
-    m = collect(Matrix{ComplexF64}, base_matrices(4))
-    for i=1:d, j=1:d
-        v = tr(m[i]' * m[j])
-        i == j ? @test(v == 1.) : @test(v == 0.)
-    end
+    m = collect(HermitianBasisIterator{Matrix{ComplexF64}}(d))
+    @test [tr(m[i]' * m[j]) for i=1:d, j=1:d] ≈ Matrix{Float64}(I, d, d)
 end
+
+@testset "represent, combine" begin
+    d = 4
+    A = reshape(collect(1:16), 4, 4)
+    vA = represent(HermitianBasis{Matrix{ComplexF64}}(d), A)
+    Ap = combine(HermitianBasis{Matrix{ComplexF64}}(d), vA)
+    @test A ≈ Ap
+    B = A*A'
+    vB = represent(HermitianBasis{Matrix{ComplexF64}}(d), B)
+    Bp = combine(HermitianBasis{Matrix{ComplexF64}}(d), vB)
+    @test B ≈ Bp
+end
+
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index a59ac36..8183d73 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,8 +6,8 @@ using SparseArrays
 using Test
 
 my_tests = ["utils.jl", "base.jl", "ptrace.jl", "ptranspose.jl", "reshuffle.jl",
-            "channels.jl", "functionals.jl", "gates.jl", "permute_systems.jl",
-            "randomqobjects.jl", "convex.jl"]
+            "channels.jl", "functionals.jl", "gates.jl", "matrixbases.jl",
+            "permute_systems.jl", "randomqobjects.jl", "convex.jl"]
 
 for my_test in my_tests
     include(my_test)
diff --git a/test/utils.jl b/test/utils.jl
index 16d5b87..3b192cf 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -13,12 +13,13 @@ end
 end
 
 @testset "renormalize" begin
-    v = randn(10)
+    rng = MersenneTwister(1234);
+    v = randn(rng, 10)
     renormalize!(v)
 
     @test norm(v) ≈ 1 atol=1e-13
 
-    A = randn(10, 10)
+    A = randn(rng, 10, 10)
     renormalize!(A)
 
     @test tr(A) ≈ 1 atol=1e-13