From ad3b7a3d737fa081456d81fdd04477863e08f1a5 Mon Sep 17 00:00:00 2001 From: Akira Kyle Date: Thu, 29 Jun 2023 11:09:59 -0600 Subject: [PATCH] Implement KrausOperators, conversions, tensor, checks, tests --- src/QuantumOpticsBase.jl | 6 +- src/pauli.jl | 4 +- src/superoperators.jl | 253 ++++++++++++++++++++++++++++++++++++ test/test_superoperators.jl | 95 ++++++++++++-- 4 files changed, 343 insertions(+), 15 deletions(-) diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl index 0c931c72..6041e2c3 100644 --- a/src/QuantumOpticsBase.jl +++ b/src/QuantumOpticsBase.jl @@ -36,8 +36,10 @@ export Basis, GenericBasis, CompositeBasis, basis, current_time, time_shift, time_stretch, time_restrict, static_operator, #superoperators SuperOperator, DenseSuperOperator, DenseSuperOpType, - SparseSuperOperator, SparseSuperOpType, spre, spost, sprepost, liouvillian, - identitysuperoperator, + SparseSuperOperator, SparseSuperOpType, ChoiState, KrausOperators, + canonicalize, orthogonalize, make_trace_preserving, is_cptp, is_cptni, + is_completely_positive, is_trace_preserving, is_trace_nonincreasing, + spre, spost, sprepost, liouvillian, identitysuperoperator, #fock FockBasis, number, destroy, create, fockstate, coherentstate, coherentstate!, diff --git a/src/pauli.jl b/src/pauli.jl index a43c1d82..f975db4c 100644 --- a/src/pauli.jl +++ b/src/pauli.jl @@ -257,9 +257,9 @@ ChiMatrix(ptm::DensePauliTransferMatrix) = ChiMatrix(SuperOperator(ptm)) """Equality for all varieties of superoperators.""" ==(sop1::T, sop2::T) where T<:Union{DensePauliTransferMatrix, DenseSuperOpType, DenseChiMatrix} = sop1.data == sop2.data -==(sop1::Union{DensePauliTransferMatrix, DenseSuperOpType, DenseChiMatrix}, sop2::Union{DensePauliTransferMatrix, DenseSuperOpType, DenseChiMatrix}) = false +==(sop1::Union{DensePauliTransferMatrix, DenseSuperOpType, DenseChiMatrix}, sop2::Union{DensePauliTransferMatrix, DenseChiMatrix}) = false """Approximate equality for all varieties of superoperators.""" -function isapprox(sop1::T, sop2::T; kwargs...) where T<:Union{DensePauliTransferMatrix, DenseSuperOpType, DenseChiMatrix} +function isapprox(sop1::T, sop2::T; kwargs...) where T<:Union{DensePauliTransferMatrix, DenseChiMatrix} return isapprox(sop1.data, sop2.data; kwargs...) end diff --git a/src/superoperators.jl b/src/superoperators.jl index e174e09b..94cb26ef 100644 --- a/src/superoperators.jl +++ b/src/superoperators.jl @@ -1,6 +1,21 @@ +import Base: isapprox import QuantumInterface: AbstractSuperOperator import FastExpm: fastExpm +# TODO: this should belong in QuantumInterface.jl +abstract type OperatorBasis{BL<:Basis,BR<:Basis} end +abstract type SuperOperatorBasis{BL<:OperatorBasis,BR<:OperatorBasis} end + +""" + tensor(E::AbstractSuperOperator, F::AbstractSuperOperator, G::AbstractSuperOperator...) + +Tensor product ``\\mathcal{E}⊗\\mathcal{F}⊗\\mathcal{G}⊗…`` of the given super operators. +""" +tensor(a::AbstractSuperOperator, b::AbstractSuperOperator) = arithmetic_binary_error("Tensor product", a, b) +tensor(op::AbstractSuperOperator) = op +tensor(operators::AbstractSuperOperator...) = reduce(tensor, operators) + + """ SuperOperator <: AbstractSuperOperator @@ -69,6 +84,9 @@ sparse(a::SuperOperator) = SuperOperator(a.basis_l, a.basis_r, sparse(a.data)) ==(a::SuperOperator{B1,B2}, b::SuperOperator{B1,B2}) where {B1,B2} = (samebases(a,b) && a.data == b.data) ==(a::SuperOperator, b::SuperOperator) = false +isapprox(a::SuperOperator{B1,B2}, b::SuperOperator{B1,B2}; kwargs...) where {B1,B2} = + (samebases(a,b) && isapprox(a.data, b.data; kwargs...)) +isapprox(a::SuperOperator, b::SuperOperator; kwargs...) = false Base.length(a::SuperOperator) = length(a.basis_l[1])*length(a.basis_l[2])*length(a.basis_r[1])*length(a.basis_r[2]) samebases(a::SuperOperator, b::SuperOperator) = samebases(a.basis_l[1], b.basis_l[1]) && samebases(a.basis_l[2], b.basis_l[2]) && @@ -318,10 +336,14 @@ end throw(IncompatibleBases()) end +# TODO should all of PauliTransferMatrix, ChiMatrix, ChoiState, and KrausOperators subclass AbstractSuperOperator? """ ChoiState <: AbstractSuperOperator Superoperator represented as a choi state. + +The convention is chosen such that the input operators live in `(basis_l[1], basis_r[1])` while +the output operators live in `(basis_r[2], basis_r[2])`. """ mutable struct ChoiState{B1,B2,T} <: AbstractSuperOperator{B1,B2} basis_l::B1 @@ -344,6 +366,27 @@ 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)) +isapprox(a::ChoiState, b::ChoiState; kwargs...) = isapprox(SuperOperator(a), SuperOperator(b); kwargs...) + +# Container to hold each of the four bases for a Choi operator when converting it to +# an operator so that if any are CompositeBases tensor doesn't lossily collapse them +struct ChoiSubBasis{S,B<:Basis} <: Basis + shape::S + basis::B +end +ChoiSubBasis(b::Basis) = ChoiSubBasis(b.shape, b) + +# TODO: decide whether to document and export this +choi_to_operator(c::ChoiState) = Operator( + ChoiSubBasis(c.basis_l[2])⊗ChoiSubBasis(c.basis_l[1]), ChoiSubBasis(c.basis_r[2])⊗ChoiSubBasis(c.basis_r[1]), c.data) + +function tensor(a::ChoiState, b::ChoiState) + op = choi_to_operator(a) ⊗ choi_to_operator(b) + op = permutesystems(op, [1,3,2,4]) + ChoiState((a.basis_l[1] ⊗ b.basis_l[1], a.basis_l[2] ⊗ b.basis_l[2]), + (a.basis_r[1] ⊗ b.basis_r[1], a.basis_r[2] ⊗ b.basis_r[2]), op.data) +end +tensor(a::SuperOperator, b::SuperOperator) = SuperOperator(tensor(ChoiState(a), ChoiState(b))) # reshape swaps within systems due to colum major ordering # https://docs.qojulia.org/quantumobjects/operators/#tensor_order @@ -372,6 +415,216 @@ 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)...) + +""" + KrausOperators <: AbstractSuperOperator + +Superoperator represented as a list of Kraus operators. + +Note that KrausOperators can only represent linear maps taking density operators to other +(potentially unnormalized) density operators. +In contrast the `SuperOperator` or `ChoiState` representations can represent arbitrary linear maps +taking arbitrary operators defined on ``H_A \\to H_B`` to ``H_C \\to H_D``. +In otherwords, the Kraus representation is only defined for completely positive linear maps of the form +``(H_A \\to H_A) \\to (H_B \\to H_B)``. +Thus converting from `SuperOperator` or `ChoiState` to `KrausOperators` will throw an exception if the +map cannot be faithfully represented up to the specificed tolerance `tol`. + +---------------------------- +Old text: +Note unlike the SuperOperator or ChoiState types where it is possible to have +`basis_l[1] != basis_l[2]` and `basis_r[1] != basis_r[2]` +which allows representations of maps between general linear operators defined on ``H_A \\to H_B``, +a quantum channel can only act on valid density operators which live in ``H_A \\to H_A``. +Thus the Kraus representation is only defined for quantum channels which map +``(H_A \\to H_A) \\to (H_B \\to H_B)``. +""" +mutable struct KrausOperators{B1,B2,T} <: AbstractSuperOperator{B1,B2} + basis_l::B1 + basis_r::B2 + data::Vector{T} + function KrausOperators{BL,BR,T}(basis_l::BL, basis_r::BR, data::Vector{T}) where {BL,BR,T} + if (any(!samebases(basis_r, M.basis_r) for M in data) || + any(!samebases(basis_l, M.basis_l) for M in data)) + throw(DimensionMismatch("Tried to assign data with incompatible bases")) + end + + new(basis_l, basis_r, data) + end +end +KrausOperators{BL,BR}(b1::BL,b2::BR,data::Vector{T}) where {BL,BR,T} = KrausOperators{BL,BR,T}(b1,b2,data) +KrausOperators(b1::BL,b2::BR,data::Vector{T}) where {BL,BR,T} = KrausOperators{BL,BR,T}(b1,b2,data) + +dense(a::KrausOperators) = KrausOperators(a.basis_l, a.basis_r, [dense(op) for op in a.data]) +sparse(a::KrausOperators) = KrausOperators(a.basis_l, a.basis_r, [sparse(op) for op in a.data]) +dagger(a::KrausOperators) = KrausOperators(a.basis_r, a.basis_l, [dagger(op) for op in a.data]) +*(a::KrausOperators{B1,B2}, b::KrausOperators{B2,B3}) where {B1,B2,B3} = + KrausOperators(a.basis_l, b.basis_r, [A*B for A in a.data for B in b.data]) +*(a::KrausOperators, b::KrausOperators) = throw(IncompatibleBases()) +*(a::KrausOperators{BL,BR}, b::Operator{BR,BR}) where {BL,BR} = sum(op*b*dagger(op) for op in a.data) +==(a::KrausOperators, b::KrausOperators) = (SuperOperator(a) == SuperOperator(b)) +isapprox(a::KrausOperators, b::KrausOperators; kwargs...) = isapprox(SuperOperator(a), SuperOperator(b); kwargs...) +tensor(a::KrausOperators, b::KrausOperators) = + KrausOperators(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r, + [A ⊗ B for A in a.data for B in b.data]) + +""" + orthogonalize(kraus::KrausOperators; tol=√eps) + +Orthogonalize the set kraus operators by performing a qr decomposition on their vec'd operators. +Note that this is different than `canonicalize` which returns a kraus decomposition such +that the kraus operators are Hilbert–Schmidt orthorgonal. + +If the input dimension is d and output dimension is d' then the number of kraus +operators returned is guaranteed to be no greater than dd', however it may be greater +than the Kraus rank. + +`orthogonalize` should always be much faster than canonicalize as it avoids an explicit eigendecomposition +and thus also preserves sparsity if the kraus operators are sparse. +""" +function orthogonalize(kraus::KrausOperators; tol=_get_tol(kraus)) + bl, br = kraus.basis_l, kraus.basis_r + dim = length(bl)*length(br) + + A = stack(reshape(op.data, dim) for op in kraus.data; dims=1) + F = qr(A; tol=tol) + # rank(F) for some reason doesn't work but should + rank = maximum(findnz(F.R)[1]) + # Sanity checks that help illustrate what qr() returns: + # @assert (F.R ≈ (sparse(F.Q') * A[F.prow,F.pcol])[1:dim,:]) + # @assert (all(iszero(F.R[rank+1:end,:]))) + + ops = [Operator(bl, br, copy(reshape( # copy materializes reshaped view + F.R[i,invperm(F.pcol)], (length(bl), length(br))))) for i=1:rank] + return KrausOperators(bl, br, ops) +end + +""" + canonicalize(kraus::KrausOperators; tol=√eps) + +Transform the quantum channel into canonical form such that the kraus operators ``{A_k}`` +are Hilbert–Schmidt orthorgonal: + +```math +\\Tr A_i^\\dagger A_j \\sim \\delta_{i,j} +``` + +If the input dimension is d and output dimension is d' then the number of kraus +operators returned is guaranteed to be no greater than dd' and will furthermore +be equal the Kraus rank of the channel up to numerical imprecision controlled by `tol`. +""" +canonicalize(kraus::KrausOperators; tol=_get_tol(kraus)) = KrausOperators(ChoiState(kraus); tol=tol) + +# TODO: document +function make_trace_preserving(kraus; tol=_get_tol(kraus)) + m = I - sum(dagger(M)*M for M in kraus.data).data + if isa(_positive_eigen(m, tol), Number) + throw(ArgumentError("Channel must be trace nonincreasing")) + end + K = Operator(kraus.basis_l, kraus.basis_r, sqrt(Matrix(m))) + return KrausOperators(kraus.basis_l, kraus.basis_r, [kraus.data; K]) +end + +SuperOperator(kraus::KrausOperators) = + SuperOperator((kraus.basis_l, kraus.basis_l), (kraus.basis_r, kraus.basis_r), + (sum(conj(op)⊗op for op in kraus.data)).data) + +ChoiState(kraus::KrausOperators) = + ChoiState((kraus.basis_r, kraus.basis_l), (kraus.basis_r, kraus.basis_l), + (sum((M=op.data; reshape(M, (length(M), 1))*reshape(M, (1, length(M)))) + for op in kraus.data))) + +_choi_state_maps_density_ops(choi::ChoiState) = (samebases(choi.basis_l[1], choi.basis_r[1]) && + samebases(choi.basis_l[2], choi.basis_r[2])) + +# TODO: consider using https://github.com/jlapeyre/IsApprox.jl +_is_hermitian(M, tol) = ishermitian(M) || isapprox(M, M', atol=tol) +_is_identity(M, tol) = isapprox(M, I, atol=tol) + +# TODO: document +# data must be Hermitian! +function _positive_eigen(data, tol) + # LinearAlgebra's eigen returns eigenvals sorted smallest to largest for Hermitian matrices + vals, vecs = eigen(Hermitian(Matrix(data))) + vals[1] < -tol && return vals[1] + ret = [(val, vecs[:,i]) for (i, val) in enumerate(vals) if val > tol] + return ret +end + +function KrausOperators(choi::ChoiState; tol=_get_tol(choi)) + if !_choi_state_maps_density_ops(choi) + throw(DimensionMismatch("Tried to convert Choi state of something that isn't a quantum channel mapping density operators to density operators")) + end + if !_is_hermitian(choi.data, tol) + throw(ArgumentError("Tried to convert nonhermitian Choi state")) + end + bl, br = choi.basis_l[2], choi.basis_l[1] + eigs = _positive_eigen(choi.data, tol) + if isa(eigs, Number) + throw(ArgumentError("Tried to convert a non-positive semidefinite Choi state,"* + "failed for smallest eigval $(eigs), consider increasing tol=$(tol)")) + end + + ops = [Operator(bl, br, sqrt(val)*reshape(vec, length(bl), length(br))) for (val, vec) in eigs] + return KrausOperators(bl, br, ops) +end + +KrausOperators(op::SuperOperator; tol=_get_tol(op)) = KrausOperators(ChoiState(op); tol=tol) + +# TODO: document superoperator representation precident: everything of mixed type returns SuperOperator *(a::ChoiState, b::SuperOperator) = SuperOperator(a)*b *(a::SuperOperator, b::ChoiState) = a*SuperOperator(b) +*(a::KrausOperators, b::SuperOperator) = SuperOperator(a)*b +*(a::SuperOperator, b::KrausOperators) = a*SuperOperator(b) +*(a::KrausOperators, b::ChoiState) = SuperOperator(a)*SuperOperator(b) +*(a::ChoiState, b::KrausOperators) = SuperOperator(a)*SuperOperator(b) + +_get_tol(kraus::KrausOperators) = sqrt(eps(real(eltype(eltype(fieldtypes(typeof(kraus))[3]))))) +_get_tol(super::SuperOperator) = sqrt(eps(real(eltype(fieldtypes(typeof(super))[3])))) +_get_tol(super::ChoiState) = sqrt(eps(real(eltype(fieldtypes(typeof(super))[3])))) + +# TODO: document this +is_completely_positive(choi::KrausOperators; tol=_get_tol(choi)) = true + +function is_completely_positive(choi::ChoiState; tol=_get_tol(choi)) + _choi_state_maps_density_ops(choi) || return false + _is_hermitian(choi.data, tol) || return false + isa(_positive_eigen(Hermitian(choi.data), tol), Number) && return false + return true +end + +is_completely_positive(super::SuperOperator; tol=_get_tol(super)) = + is_completely_positive(ChoiState(super); tol=tol) + +# TODO: document this +is_trace_preserving(kraus::KrausOperators; tol=_get_tol(kraus)) = + _is_identity(sum(dagger(M)*M for M in kraus.data).data, tol) + +is_trace_preserving(choi::ChoiState; tol=_get_tol(choi)) = + _is_identity(ptrace(choi_to_operator(choi), 1).data, tol) + +is_trace_preserving(super::SuperOperator; tol=_get_tol(super)) = + is_trace_preserving(ChoiState(super); tol=tol) + +# TODO: document this +function is_trace_nonincreasing(kraus::KrausOperators; tol=_get_tol(kraus)) + m = I - sum(dagger(M)*M for M in kraus.data).data + _is_hermitian(m, tol) || return false + return !isa(_positive_eigen(Hermitian(m), tol), Number) +end + +function is_trace_nonincreasing(choi::ChoiState; tol=_get_tol(choi)) + m = I - ptrace(choi_to_operator(choi), 1).data + _is_hermitian(m, tol) || return false + return !isa(_positive_eigen(Hermitian(m), tol), Number) +end + +is_trace_nonincreasing(super::SuperOperator; tol=_get_tol(super)) = + is_trace_nonincreasing(ChoiState(super); tol=tol) + +# TODO: document this +is_cptp(sop; tol=_get_tol(sop)) = is_completely_positive(sop; tol=tol) && is_trace_preserving(sop; tol=tol) + +# TODO: document this +is_cptni(sop; tol=_get_tol(sop)) = is_completely_positive(sop; tol=tol) && is_trace_nonincreasing(sop; tol=tol) diff --git a/test/test_superoperators.jl b/test/test_superoperators.jl index dd968c17..0843b7f1 100644 --- a/test/test_superoperators.jl +++ b/test/test_superoperators.jl @@ -205,6 +205,7 @@ end @test tracedistance(L*ρ₀, ρ) < 1e-10 @test tracedistance(ChoiState(L)*ρ₀, ρ) < 1e-10 +# TODO: reenable these now that QuantumOptics.jl is a test dependency? # 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 @@ -248,20 +249,92 @@ 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 +# 0-2-4 binomial code encoder 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)) +enc_proj = z_l ⊗ dagger(spinup(b_logical)) + o_l ⊗ dagger(spindown(b_logical)) +dec_proj = dagger(enc_proj) +enc_sup = sprepost(enc_proj, dec_proj) +dec_sup = sprepost(dec_proj, enc_proj) +enc_kraus = KrausOperators(b_fock, b_logical, [enc_proj]) +dec_kraus = KrausOperators(b_logical, b_fock, [dec_proj]) +## testing conversions +@test dec_sup == dagger(enc_sup) +@test dec_kraus == dagger(enc_kraus) +@test ChoiState(enc_sup) == ChoiState(enc_kraus) +@test ChoiState(dec_sup) == dagger(ChoiState(enc_sup)) +@test ChoiState(dec_kraus) == dagger(ChoiState(enc_kraus)) +@test SuperOperator(ChoiState(enc_sup)) == enc_sup +@test SuperOperator(KrausOperators(enc_sup)) ≈ enc_sup +@test KrausOperators(ChoiState(enc_kraus)) ≈ enc_kraus +@test KrausOperators(SuperOperator(enc_kraus)) ≈ enc_kraus +## testing multipication +@test dec_sup*enc_sup ≈ dense(identitysuperoperator(b_logical)) +@test SuperOperator(dec_kraus*enc_kraus) ≈ dense(identitysuperoperator(b_logical)) +@test dec_sup*enc_kraus ≈ dense(identitysuperoperator(b_logical)) +@test dec_kraus*enc_sup ≈ dense(identitysuperoperator(b_logical)) +@test SuperOperator(dec_kraus*enc_kraus) ≈ dense(identitysuperoperator(b_logical)) +@test dec_sup*ChoiState(enc_sup) ≈ dense(identitysuperoperator(b_logical)) +@test ChoiState(dec_sup)*enc_sup ≈ dense(identitysuperoperator(b_logical)) +@test SuperOperator(ChoiState(dec_sup)*ChoiState(enc_sup)) ≈ dense(identitysuperoperator(b_logical)) +## testing channel checks +@test is_cptp(enc_kraus) +@test is_cptni(dec_kraus) +@test is_cptp(enc_sup) +@test is_cptni(dec_sup) +@test is_cptp(ChoiState(enc_kraus)) +@test is_cptni(ChoiState(dec_kraus)) + +# TODO: fix avg_gate_fidelity to work with all superoperator types +#@test avg_gate_fidelity(dec_sup*enc_sup) ≈ 1 +#@test avg_gate_fidelity(dec_kraus*enc_kraus) ≈ 1 +#@test avg_gate_fidelity(ChoiState(dec_sup)*ChoiState(enc_sup)) ≈ 1 + +# test qubit amplitude and phase damping channels +function ampl_damp_kraus(γ) + b = SpinBasis(1//2) + K0 = dm(spindown(b)) + sqrt(1-γ)*dm(spinup(b)) + K1 = sqrt(γ)*sigmam(b) + return KrausOperators(b,b,[K0, K1]) +end + +function phase_damp_kraus(γ) + b = SpinBasis(1//2) + K0 = dm(spindown(b)) + sqrt(1-γ)*dm(spinup(b)) + K1 = sqrt(γ)*dm(spinup(b)) + return KrausOperators(b,b,[K0, K1]) +end + +tensor_it(x, N) = tensor(fill(x, N)...) + +function test_channel(κa, κp, N, do_sparse) + tol = 1e-6 + b = SpinBasis(1//2) + La = liouvillian(identityoperator(b), [sigmam(b)]; rates=[κa]) + Lp = liouvillian(identityoperator(b), [sigmaz(b)]; rates=[κp]) + super = tensor_it(exp(La + Lp), N) + super = do_sparse ? sparse(super) : dense(super) + ka = ampl_damp_kraus(1-exp(-κa)) + kp = phase_damp_kraus(1-exp(-4*κp)) + kraus = tensor_it(ka*kp, N) + kraus = do_sparse ? sparse(kraus) : dense(kraus) + @test is_cptp(super; tol=tol) + @test is_cptp(kraus; tol=tol) + @test isapprox(SuperOperator(kraus), super; atol=tol) + @test isapprox(ChoiState(kraus), ChoiState(super); atol=tol) + @test isapprox(kraus, KrausOperators(super); atol=tol) +end + +for N=1:4 + for κa=[0, 1e-1, 1e-2, 1e-3] + for κp=[0, 1e-1, 1e-2, 1e-3] + test_channel(κa, κp, N, false) + test_channel(κa, κp, N, true) + end + end +end + end # testset