Merge branch 'pg/bases' of github.com:ZKSI/QI.jl into pg/bases

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.
 

src/base.jl

@@ -1,21 +1,34 @@
-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.
 - `dim`: length of the vector.
 
 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)

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
 

src/gates.jl

@@ -1,4 +1,3 @@
-
 export 𝕀, sx, sy, sz, qft, hadamard, grover
 
 sx = ComplexF64[0 1; 1 0]

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

src/ptrace.jl

@@ -1,3 +1,5 @@
+export ptrace
+
 """
 $(SIGNATURES)
 - `ρ`: quantum state.

src/ptranspose.jl

@@ -1,3 +1,4 @@
+export ptranspose
 """
 $(SIGNATURES)
 - `ρ`: quantum state.

src/reshuffle.jl

@@ -1,3 +1,5 @@
+export reshuffle
+
 # """
 #   Performs reshuffling of indices of a matrix.
 #   Given multiindexed matrix M_{(m,μ),(n,ν)} it returns

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

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

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

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)

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