Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Create ChoiState type and conversions to/from SuperOperator #115

Merged
merged 4 commits into from
Nov 17, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "QuantumOpticsBase"
uuid = "4f57444f-1401-5e15-980d-4471b28d5678"
version = "0.5.3"
version = "0.5.4"

[deps]
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
Expand Down
63 changes: 61 additions & 2 deletions src/superoperators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,9 @@ mutable struct SuperOperator{B1,B2,T} <: AbstractSuperOperator{B1,B2}
basis_r::B2
data::T
function SuperOperator{BL,BR,T}(basis_l::BL, basis_r::BR, data::T) where {BL,BR,T}
if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
length(basis_r[1])*length(basis_r[2]) != size(data, 2)
if (length(basis_l) != 2 || length(basis_r) != 2 ||
length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
length(basis_r[1])*length(basis_r[2]) != size(data, 2))
throw(DimensionMismatch("Tried to assign data of size $(size(data)) to Hilbert spaces of sizes $(length.(basis_l)), $(length.(basis_r))"))
end
new(basis_l, basis_r, data)
Expand Down Expand Up @@ -316,3 +317,61 @@ end
}
throw(IncompatibleBases())
end

"""
ChoiState <: AbstractSuperOperator

Superoperator represented as a choi state.
"""
mutable struct ChoiState{B1,B2,T} <: AbstractSuperOperator{B1,B2}
basis_l::B1
basis_r::B2
data::T
function ChoiState{BL,BR,T}(basis_l::BL, basis_r::BR, data::T) where {BL,BR,T}
if (length(basis_l) != 2 || length(basis_r) != 2 ||
length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
length(basis_r[1])*length(basis_r[2]) != size(data, 2))
throw(DimensionMismatch("Tried to assign data of size $(size(data)) to Hilbert spaces of sizes $(length.(basis_l)), $(length.(basis_r))"))
end
new(basis_l, basis_r, data)
end
end
ChoiState(b1::BL, b2::BR, data::T) where {BL,BR,T} = ChoiState{BL,BR,T}(b1, b2, data)

dense(a::ChoiState) = ChoiState(a.basis_l, a.basis_r, Matrix(a.data))
sparse(a::ChoiState) = ChoiState(a.basis_l, a.basis_r, sparse(a.data))
dagger(a::ChoiState) = ChoiState(dagger(SuperOperator(a)))
*(a::ChoiState, b::ChoiState) = ChoiState(SuperOperator(a)*SuperOperator(b))
*(a::ChoiState, b::Operator) = SuperOperator(a)*b
==(a::ChoiState, b::ChoiState) = (SuperOperator(a) == SuperOperator(b))

# reshape swaps within systems due to colum major ordering
# https://docs.qojulia.org/quantumobjects/operators/#tensor_order
function _super_choi((l1, l2), (r1, r2), data)
data = reshape(data, map(length, (l2, l1, r2, r1)))
(l1, l2), (r1, r2) = (r2, l2), (r1, l1)
data = permutedims(data, (1, 3, 2, 4))
data = reshape(data, map(length, (l1⊗l2, r1⊗r2)))
return (l1, l2), (r1, r2), data
end

function _super_choi((r2, l2), (r1, l1), data::SparseMatrixCSC)
data = _permutedims(data, map(length, (l2, r2, l1, r1)), (1, 3, 2, 4))
data = reshape(data, map(length, (l1⊗l2, r1⊗r2)))
# sparse(data) is necessary since reshape of a sparse array returns a
# ReshapedSparseArray which is not a subtype of AbstractArray and so
# _permutedims fails to acces the ".m" field
# https://github.com/qojulia/QuantumOpticsBase.jl/pull/83
# https://github.com/JuliaSparse/SparseArrays.jl/issues/24
# permutedims in SparseArrays.jl only implements perm (2,1) and so
# _permutedims should probably be upstreamed
# https://github.com/JuliaLang/julia/issues/26534
return (l1, l2), (r1, r2), sparse(data)
end

ChoiState(op::SuperOperator) = ChoiState(_super_choi(op.basis_l, op.basis_r, op.data)...)
SuperOperator(op::ChoiState) = SuperOperator(_super_choi(op.basis_l, op.basis_r, op.data)...)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

does it make sense to also add a Ket(::ChoiState) constructor, so that we can actually extract it as a state?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think it makes sense since in general the ChoiState is mixed iff the corresponding channel is non-unitary, so it is akin to asking for a Ket(::Operator) method which has multiple ways in which it could make sense, e.g. closest pure state under a given operator norms.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, noted. Then maybe a dm(::ChoiState) would make sense, but no need to block the PR on this.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have something like that added as part of #177. The tricky part is ensuring the input and ancilla bases make sense and can be appropriately manipulated!


*(a::ChoiState, b::SuperOperator) = SuperOperator(a)*b
*(a::SuperOperator, b::ChoiState) = a*SuperOperator(b)

36 changes: 36 additions & 0 deletions test/test_superoperators.jl
Original file line number Diff line number Diff line change
@@ -1,13 +1,15 @@
using Test
using QuantumOpticsBase
using SparseArrays, LinearAlgebra
import QuantumOpticsBase: ChoiState # Remove when ChoiState is publicly exported

@testset "superoperators" begin

# Test creation
b = GenericBasis(3)
@test_throws DimensionMismatch DenseSuperOperator((b, b), (b, b), zeros(ComplexF64, 3, 3))
@test_throws DimensionMismatch SparseSuperOperator((b, b), (b, b), spzeros(ComplexF64, 3, 3))
@test_throws DimensionMismatch ChoiState((b, b), (b, b), zeros(ComplexF64, 3, 3))

# Test copy, sparse and dense
b1 = GenericBasis(2)
Expand Down Expand Up @@ -153,29 +155,45 @@ J = [Ja, Jc]
ρ₀ = dm(Ψ₀)

@test identitysuperoperator(spinbasis)*sx == sx
@test ChoiState(identitysuperoperator(spinbasis))*sx == sx
@test identitysuperoperator(sparse(spre(sx)))*sx == sparse(sx)
@test ChoiState(identitysuperoperator(sparse(spre(sx))))*sx == sparse(sx)
@test sparse(ChoiState(identitysuperoperator(spre(sx))))*sx == sparse(sx)
@test identitysuperoperator(dense(spre(sx)))*sx == dense(sx)
@test ChoiState(identitysuperoperator(dense(spre(sx))))*sx == dense(sx)
@test dense(ChoiState(identitysuperoperator(spre(sx))))*sx == dense(sx)

op1 = DenseOperator(spinbasis, [1.2+0.3im 0.7+1.2im;0.3+0.1im 0.8+3.2im])
op2 = DenseOperator(spinbasis, [0.2+0.1im 0.1+2.3im; 0.8+4.0im 0.3+1.4im])
@test tracedistance(spre(op1)*op2, op1*op2) < 1e-12
@test tracedistance(ChoiState(spre(op1))*op2, op1*op2) < 1e-12
@test tracedistance(spost(op1)*op2, op2*op1) < 1e-12
@test tracedistance(ChoiState(spost(op1))*op2, op2*op1) < 1e-12

@test spre(sparse(op1))*op2 == op1*op2
@test ChoiState(spre(sparse(op1)))*op2 == op1*op2
@test spost(sparse(op1))*op2 == op2*op1
@test ChoiState(spost(sparse(op1)))*op2 == op2*op1
@test spre(sparse(dagger(op1)))*op2 == dagger(op1)*op2
@test ChoiState(spre(sparse(dagger(op1))))*op2 == dagger(op1)*op2
@test spre(dense(dagger(op1)))*op2 ≈ dagger(op1)*op2
@test ChoiState(spre(dense(dagger(op1))))*op2 ≈ dagger(op1)*op2
@test sprepost(sparse(op1), op1)*op2 ≈ op1*op2*op1
@test ChoiState(sprepost(sparse(op1), op1))*op2 ≈ op1*op2*op1

@test spre(sparse(op1))*sparse(op2) == sparse(op1*op2)
@test ChoiState(spre(sparse(op1)))*sparse(op2) == sparse(op1*op2)
@test spost(sparse(op1))*sparse(op2) == sparse(op2*op1)
@test ChoiState(spost(sparse(op1)))*sparse(op2) == sparse(op2*op1)
@test sprepost(sparse(op1), sparse(op1))*sparse(op2) ≈ sparse(op1*op2*op1)
@test ChoiState(sprepost(sparse(op1), sparse(op1)))*sparse(op2) ≈ sparse(op1*op2*op1)

@test sprepost(op1, op2) ≈ spre(op1)*spost(op2)
b1 = FockBasis(1)
b2 = FockBasis(5)
op = fockstate(b1, 0) ⊗ dagger(fockstate(b2, 0))
@test sprepost(dagger(op), op)*dm(fockstate(b1, 0)) == dm(fockstate(b2, 0))
@test ChoiState(sprepost(dagger(op), op))*dm(fockstate(b1, 0)) == dm(fockstate(b2, 0))
@test_throws ArgumentError spre(op)
@test_throws ArgumentError spost(op)

Expand All @@ -185,12 +203,14 @@ for j=J
ρ .+= j*ρ₀*dagger(j) - 0.5*dagger(j)*j*ρ₀ - 0.5*ρ₀*dagger(j)*j
end
@test tracedistance(L*ρ₀, ρ) < 1e-10
@test tracedistance(ChoiState(L)*ρ₀, ρ) < 1e-10

# tout, ρt = timeevolution.master([0.,1.], ρ₀, H, J; reltol=1e-7)
# @test tracedistance(exp(dense(L))*ρ₀, ρt[end]) < 1e-6
# @test tracedistance(exp(sparse(L))*ρ₀, ρt[end]) < 1e-6

@test dense(spre(op1)) == spre(op1)
@test dense(ChoiState(spre(op1))) == ChoiState(spre(op1))

@test L/2.0 == 0.5*L == L*0.5
@test -L == SparseSuperOperator(L.basis_l, L.basis_r, -L.data)
Expand Down Expand Up @@ -228,4 +248,20 @@ N = exp(log(2) * sparse(L)) # 50% loss channel
@test (0.5 - real(tr(ρ^2))) < 1e-10 # one photon state becomes maximally mixed
@test tracedistance(ρ, normalize(dm(fockstate(b, 0)) + dm(fockstate(b, 1)))) < 1e-10

# Testing 0-2-4 binomial code encoder
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is a really neat test!

b_logical = SpinBasis(1//2)
b_fock = FockBasis(5)
z_l = normalize(fockstate(b_fock, 0) + fockstate(b_fock, 4))
o_l = fockstate(b_fock, 2)
encoder_kraus = z_l ⊗ dagger(spinup(b_logical)) + o_l ⊗ dagger(spindown(b_logical))
encoder_sup = sprepost(encoder_kraus, dagger(encoder_kraus))
decoder_sup = sprepost(dagger(encoder_kraus), encoder_kraus)
@test SuperOperator(ChoiState(encoder_sup)).data == encoder_sup.data
@test decoder_sup == dagger(encoder_sup)
@test ChoiState(decoder_sup) == dagger(ChoiState(encoder_sup))
@test decoder_sup*encoder_sup ≈ dense(identitysuperoperator(b_logical))
@test decoder_sup*ChoiState(encoder_sup) ≈ dense(identitysuperoperator(b_logical))
@test ChoiState(decoder_sup)*encoder_sup ≈ dense(identitysuperoperator(b_logical))
@test SuperOperator(ChoiState(decoder_sup)*ChoiState(encoder_sup)) ≈ dense(identitysuperoperator(b_logical))

end # testset
6 changes: 5 additions & 1 deletion test/test_time_dependent_operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,11 @@ using LinearAlgebra, Random
o_t_tup = TimeDependentSum(Tuple, o_t)
@test QOB.static_operator(o_t_tup(t)).factors == [1.0im, t*3.0 + 0.0im]
@test all(QOB.static_operator(o_t_tup(t)).operators .== (a, n))
@test (@allocated set_time!(o_t_tup, t)) == 0
if VERSION.minor == 11 # issue #178 https://github.com/qojulia/QuantumOpticsBase.jl/issues/178
@test_broken (@allocated set_time!(o_t_tup, t)) == 0
else
@test (@allocated set_time!(o_t_tup, t)) == 0
end

o_t2 = TimeDependentSum(f1=>a, f2=>n)
@test o_t(t) == o_t2(t)
Expand Down
Loading