Add super_to_choi and choi_to_super

src/QuantumOpticsBase.jl

@@ -35,8 +35,9 @@ 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, is_trace_preserving, is_valid_channel,
+ spre, spost, sprepost, liouvillian, identitysuperoperator,
  #fock
  FockBasis, number, destroy, create,
  fockstate, coherentstate, coherentstate!,

src/superoperators.jl

@@ -11,8 +11,9 @@
  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)
@@ -316,3 +317,126 @@
  }
  throw(IncompatibleBases())
 end
+
+# TODO should all of PauliTransferMatrix, ChiMatrix, ChoiState, and KrausOperators subclass AbstractSuperOperator?
+
+"""
+ KrausOperators(B1, B2, data)
+
+Superoperator represented as a list of Kraus operators.
+Note unlike the SuperOperator or ChoiState types where
+its 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::T
+ function KrausOperators{BL,BR,T}(basis_l::BL, basis_r::BR, data::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::T) where {BL,BR,T} = KrausOperators{BL,BR,T}(b1,b2,data)
+KrausOperators(b1::BL,b2::BR,data::T) where {BL,BR,T} = KrausOperators{BL,BR,T}(b1,b2,data)
+KrausOperators(b,data) = KrausOperators(b,b,data)
+
+function is_trace_preserving(kraus::KrausOperators; tol=1e-9)
+ m = I(length(kraus.basis_r)) - sum(dagger(M)*M for M in kraus.data).data
+ m[abs.(m) .< tol] .= 0
+ return iszero(m)
+end
+
+function is_valid_channel(kraus::KrausOperators; tol=1e-9)
+ m = I(length(kraus.basis_r)) - sum(dagger(M)*M for M in kraus.data).data
+ eigs = eigvals(Matrix(m))
+ eigs[@. abs(eigs) < tol || eigs > 0] .= 0
+ return iszero(eigs)
+end
+
+"""
+ ChoiState(B1, B2, data)
+
+Superoperator represented as a choi state stored as a sparse or dense matrix.
+"""
+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; tol=1e-9) 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
+ if any(abs.(data - data') .> tol)
+ @warn "Trying to construct ChoiState from non-hermitian data"
+ end
+ new(basis_l, basis_r, Hermitian(data))
+ end
+end
+ChoiState{BL,BR}(b1::BL,b2::BR,data::T; tol=1e-9) where {BL,BR,T} = ChoiState{BL,BR,T}(b1,b2,data;tol=tol)
+ChoiState(b1::BL,b2::BR,data::T; tol=1e-9) where {BL,BR,T} = ChoiState{BL,BR,T}(b1,b2,data; tol=tol)
+ChoiState(b,data; tol=tol) = ChoiState(b,b,data; tol=tol)
+
+# TODO: document why we have super_to_choi return non-trace one density matrices.
+# https://forest-benchmarking.readthedocs.io/en/latest/superoperator_representations.html
+# Note the similarity to permutesystems in operators_dense.jl
+
+# 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::Matrix)
+ 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 be upstreamed
+ # https://github.com/JuliaLang/julia/issues/26534
+ return (l1, l2), (r1, r2), sparse(data)
+end
+
+ChoiState(op::SuperOperator; tol=1e-9) = ChoiState(_super_choi(op.basis_l, op.basis_r, op.data)...; tol=tol)
+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)
+
+SuperOperator(op::ChoiState) = SuperOperator(_super_choi(op.basis_l, op.basis_r, op.data)...)
+ChoiState(kraus::KrausOperators) = ChoiState(SuperOperator(kraus))
+
+function KrausOperators(choi::ChoiState; tol=1e-9)
+ if (!samebases(choi.basis_l[1], choi.basis_l[2]) ||
+ !samebases(choi.basis_r[1], choi.basis_r[2]))
+ throw(DimensionMismatch("Tried to convert choi state of something that isn't a quantum channel mapping density operators to density operators"))
+ end
+ bl, br = choi.basis_l[1], choi.basis_r[1]
+ #ishermitian(choi.data) || @warn "ChoiState is not hermitian"
+ # TODO: figure out how to do this with sparse matrices using e.g. Arpack.jl or ArnoldiMethod.jl
+ vals, vecs = eigen(Hermitian(Matrix(choi.data)))
+ for val in vals
+ (abs(val) > tol && val < 0) && @warn "eigval $(val) < 0 but abs(eigval) > tol=$(tol)"
+ end
+ ops = [Operator(bl, br, sqrt(val)*reshape(vecs[:,i], length(bl), length(br)))
+ for (i, val) in enumerate(vals) if abs(val) > tol && val > 0]
+ return KrausOperators(bl, br, ops)
+end
+
+KrausOperators(op::SuperOperator; tol=1e-9) = KrausOperators(ChoiState(op; tol=tol); tol=tol)