Pl/channelbases

src/matrixbases.jl

@@ -1,6 +1,7 @@
 import .Base: length, iterate
 export AbstractMatrixBasisIterator, HermitianBasisIterator, AbstractBasis,
-    AbstractMatrixBasis, HermitianBasis, hermitianbasis, represent, combine
+       ChannelBasisIterator, AbstractChannelBasis, ChannelBasis, AbstractMatrixBasis,
+       HermitianBasis, hermitianbasis, represent, combine, channelbasis
 abstract type AbstractMatrixBasisIterator{T<:AbstractMatrix} end
 struct HermitianBasisIterator{T} <: AbstractMatrixBasisIterator{T} 
     dim::Int
@@ -17,17 +18,17 @@ struct HermitianBasis{T} <: AbstractMatrixBasis{T}
     end
 end
 
+
+
 """
 $(SIGNATURES)
 - `dim`: dimensions of the matrix.
-
 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
@@ -55,6 +56,72 @@ function represent(basis::Type{T}, m::Matrix{<:Number}) where T<:AbstractMatrixB
     represent(basis{typeof(m)}(d), m)
 end
 
-function combine(basis::T, v::Vector{<:Number}) where T<:AbstractMatrixBasis
+function combine(basis::AbstractMatrixBasis{T}, v::Vector{<:Number}) where T<:AbstractMatrix{<:Number}
     sum(basis.iterator .* v)
+end
+
+"""
+$(SIGNATURES)
+
+"""
+abstract type AbstractMatrixBasisIterator{T<:AbstractMatrix} end
+struct ChannelBasisIterator{T} <: AbstractMatrixBasisIterator{T} 
+    idim::Int 
+    odim::Int
+    hitr::HermitianBasisIterator{T}
+    function ChannelBasisIterator{T}(idim::Int, odim::Int) where T<:AbstractMatrix{<:Number}
+        new(idim, odim, HermitianBasisIterator{T}(idim))
+    end
+end
+
+abstract type AbstractChannelBasis{T} <: AbstractMatrixBasis{T} end 
+struct ChannelBasis{T} <: AbstractMatrixBasis{T} 
+    iterator::ChannelBasisIterator{T}
+    function ChannelBasis{T}(idim::Integer, odim::Integer) where T<:AbstractMatrix{<:Number}
+        new(ChannelBasisIterator{T}(idim, odim)) 
+    end
+end
+channelbasis(T::Type{<:AbstractMatrix{<:Number}}, idim::Int, odim::Int=idim) = ChannelBasis{T}(idim, odim)
+
+channelbasis(idim::Int, odim::Int=idim) = channelbasis(Matrix{ComplexF64}, idim, odim)
+
+function iterate(itr::ChannelBasisIterator{T}, state=(1,1,1,1)) where T<:AbstractMatrix{<:Number}
+    idim = itr.idim
+    odim = itr.odim
+    hitr = itr.hitr
+    (a, c, b, d) = state
+    (a == odim && c == odim && d == 2) && return nothing
+
+    Tn = eltype(T)
+    if a > c
+        x = (ketbra(T, a, c, odim) ⊗ ketbra(T, b, d, idim) + ketbra(T, c, a, odim) ⊗ ketbra(T, d, b, idim)) / sqrt(Tn(2))
+    elseif a < c  
+        x = (im * ketbra(T, a, c, odim) ⊗ ketbra(T, b, d, idim) - im * ketbra(T, c, a, odim) ⊗ ketbra(T, d, b, idim)) / sqrt(Tn(2))
+    elseif a < odim 
+        H = iterate(hitr, (b, d))[1]
+        x = (diagm(0 => vcat(ones(Tn, a), Tn[-a], zeros(Tn, odim - a-1))) ⊗ H) / sqrt(Tn(a + a^2))
+    else 
+        x = Matrix{Tn}(I, idim * odim, idim * odim) / sqrt(Tn(idim * odim))
+    end
+    if d < idim
+        newstate = (a, c, b, d+1)
+    elseif d == idim && b < idim
+        newstate = (a, c, b+1, 1)
+    elseif d == idim && b == idim && c < odim
+        newstate = (a, c+1, 1, 1)
+    else  
+        newstate = (a+1, 1, 1, 1)
+    end 
+    return x, newstate
+end 
+length(itr::ChannelBasisIterator) = itr.idim^2 * itr.odim^2 - itr.idim^2 + 1
+
+function represent(basis::ChannelBasis{T1}, Φ::AbstractQuantumOperation{T2}) where T1<:AbstractMatrix{<:Number} where T2<:AbstractMatrix{<:Number}
+    J = convert(DynamicalMatrix{T2}, Φ)
+    represent(basis, J.matrix) 
+end
+
+function combine(basis::ChannelBasis{T}, v::Vector{<:Number}) where T<:AbstractMatrix
+    m = sum(basis.iterator .* v)
+    DynamicalMatrix{T}(m, basis.iterator.idim, basis.iterator.odim)
 end

test/matrixbases.jl

@@ -6,6 +6,14 @@
     @test [tr(m[i]' * m[j]) for i=1:d, j=1:d] ≈ Matrix{Float64}(I, d, d)
 end
 
+@testset "ChannelBasisIterator" begin
+    d1 = 2
+    d2 = 2
+    d = d1^2 * d2^2 - d1^2 + 1
+    m = collect(ChannelBasisIterator{Matrix{ComplexF64}}(d1,d2))
+    @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), d, d) + reshape(collect(1:16), d, d)'
@@ -27,9 +35,41 @@ end
     @test length(vC) == prod(size(C))
 end
 
+@testset "represent, combine" begin
+    d1 = 2
+    d2 = 4
+    A = reshape(collect(1:16), d1 * d2, d1) * reshape(collect(1:16), d1 * d2, d1)'
+    B = Matrix{Float64}(I, d2, d2) ⊗ (ptrace(A, [d2, d1], 1))^(-1/2)
+    A = B * A * B'
+    vA = represent(ChannelBasis{Matrix{ComplexF64}}(d1, d2), A)
+    Ap = combine(ChannelBasis{Matrix{ComplexF64}}(d1, d2), vA)
+    @test A ≈ Ap.matrix
+
+    A = reshape(collect(1:64), d1 * d2, d1 * d2) * reshape(collect(1:64), d1 * d2, d1 * d2)' + Matrix{Float64}(I, d1 * d2, d1 * d2)
+    B = Matrix{Float64}(I, d1, d1) ⊗ (ptrace(A, [d1, d2], 1))^(-1/2)
+    B = B * A * B'
+    vB = represent(ChannelBasis{Matrix{ComplexF64}}(d2, d1), B)
+    Bp = combine(ChannelBasis{Matrix{ComplexF64}}(d2, d1), vB)
+    @test B ≈ Bp.matrix
+
+    #vB = represent(ChannelBasis{Matrix{ComplexF32}}(d1,d2), B)
+    #@test eltype(vB) == Float32
+
+    # C = Float16[1 2; 3 4]
+    # C += C'
+    # vC = represent(ChannelBasis, C)
+    # @test eltype(vC) == eltype(C)
+    # @test length(vC) == prod(size(C))
+end
+
 @testset "hermitainbasis" begin
     @test hermitianbasis(Matrix{Float32}, 2) == HermitianBasisIterator{Matrix{Float32}}(2)
     @test hermitianbasis(2) == HermitianBasisIterator{Matrix{ComplexF64}}(2)
 end
 
+# @testset "channelbasis" begin
+#    @test channelbasis(Matrix{Float32}, 2,2) == ChannelBasisIterator{Matrix{Float32}}(2,2)
+#    @test channelbasis(2,2) == ChannelBasisIterator{Matrix{ComplexF64}}(2,2)
+# end
+
 end

test/runtests.jl

@@ -8,7 +8,6 @@ using Test
 my_tests = ["utils.jl", "base.jl", "ptrace.jl", "ptranspose.jl", "reshuffle.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)
 end