Merge 807aa09 into b46911d

src/QuantumOptics.jl

@@ -35,7 +35,8 @@ export bases, Basis, GenericBasis, CompositeBasis, basis,
         steadystate,
         timecorrelations,
         semiclassical,
-        stochastic
+        stochastic, PauliBasis, PauliTransferMatrix, DensePauliTransferMatrix,
+        ChiMatrix, DenseChiMatrix
 
 
 include("sortedindices.jl")
@@ -86,6 +87,7 @@ module stochastic
     using .stochastic_definitions
 end
 include("printing.jl")
+include("pauli.jl")
 
 using .bases
 using .states
@@ -109,6 +111,7 @@ using .metrics
 using .spectralanalysis
 using .timecorrelations
 using .printing
+using .pauli
 
 
 end # module

src/pauli.jl

@@ -0,0 +1,231 @@
+module pauli
+
+export PauliBasis, PauliTransferMatrix, DensePauliTransferMatrix,
+        ChiMatrix, DenseChiMatrix
+
+import Base: ==
+
+using ..bases, ..spin, ..superoperators
+using ..operators: identityoperator, AbstractOperator
+using ..superoperators: SuperOperator
+using ..operators_dense: DenseOperator
+using ..spin: sigmax, sigmay, sigmaz
+using SparseArrays: sparse
+using LinearAlgebra: tr
+
+"""
+    PauliBasis(num_qubits::Int)
+
+Basis for an N-qubit space where `num_qubits` specifies the number of qubits.
+The dimension of the basis is 2²ᴺ.
+"""
+mutable struct PauliBasis{B<:Tuple{Vararg{Basis}}} <: Basis
+    shape::Vector{Int}
+    bases::B
+    function PauliBasis(num_qubits::Int)
+        type = Tuple{(SpinBasis{1//2} for _ in 1:num_qubits)...}
+        shape = [2 for _ in 1:num_qubits]
+        bases = Tuple(SpinBasis(1//2) for _ in 1:num_qubits)
+        return new{type}(shape, bases)
+    end
+end
+==(pb1::PauliBasis, pb2::PauliBasis) = length(pb1.bases) == length(pb2.bases)
+
+"""
+    Base class for Pauli transfer matrix classes.
+"""
+abstract type PauliTransferMatrix{B1<:Tuple{PauliBasis, PauliBasis}, B2<:Tuple{PauliBasis, PauliBasis}} end
+
+
+"""
+    DensePauliTransferMatrix(B1, B2, data)
+
+DensePauliTransferMatrix stored as a dense matrix.
+"""
+mutable struct DensePauliTransferMatrix{B1<:Tuple{PauliBasis, PauliBasis},
+                                        B2<:Tuple{PauliBasis, PauliBasis},
+                                        T<:Matrix{Float64}} <: PauliTransferMatrix{B1, B2}
+    basis_l::B1
+    basis_r::B2
+    data::T
+    function DensePauliTransferMatrix(basis_l::BL, basis_r::BR, data::T) where {BL<:Tuple{PauliBasis, PauliBasis},
+                                                                                BR<:Tuple{PauliBasis, PauliBasis},
+                                                                                T<:Matrix{Float64}}
+        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
+           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
+            throw(DimensionMismatch())
+        end
+        new{BL, BR, T}(basis_l, basis_r, data)
+    end
+end
+
+PauliTransferMatrix(ptm::DensePauliTransferMatrix{B, B, Matrix{Float64}}) where B <: Tuple{PauliBasis, PauliBasis} = ptm
+
+"""
+    Base class for χ (process) matrix classes.
+"""
+abstract type ChiMatrix{B1<:Tuple{PauliBasis, PauliBasis}, B2<:Tuple{PauliBasis, PauliBasis}} end
+
+"""
+    DenseChiMatrix(b, b, data)
+
+DenseChiMatrix stored as a dense matrix.
+"""
+mutable struct DenseChiMatrix{B1<:Tuple{PauliBasis, PauliBasis},
+                              B2<:Tuple{PauliBasis, PauliBasis},
+                              T<:Matrix{ComplexF64}} <: PauliTransferMatrix{B1, B2}
+    basis_l::B1
+    basis_r::B2
+    data::T
+    function DenseChiMatrix(basis_l::BL, basis_r::BR, data::T) where {BL<:Tuple{PauliBasis, PauliBasis},
+                                                                                BR<:Tuple{PauliBasis, PauliBasis},
+                                                                                T<:Matrix{ComplexF64}}
+        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
+           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
+            throw(DimensionMismatch())
+        end
+        new{BL, BR, T}(basis_l, basis_r, data)
+    end
+end
+
+ChiMatrix(chi_matrix::DenseChiMatrix{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis} = chi_matrix
+
+# TODO MAKE A GENERATOR FUNCTION
+"""
+    pauli_operators(num_qubits::Int)
+
+Generate a list of N-qubit Pauli operators.
+"""
+function pauli_operators(num_qubits::Int)
+    pauli_funcs = (identityoperator, sigmax, sigmay, sigmaz)
+    po = []
+    for paulis in Iterators.product((pauli_funcs for _ in 1:num_qubits)...)
+        basis_vector = reduce(⊗, f(SpinBasis(1//2)) for f in paulis)
+        push!(po, basis_vector)
+    end
+    return po
+end
+
+"""
+    pauli_basis_vectors(num_qubits::Int)
+
+Generate a matrix of basis vectors in the Pauli representation given a number
+of qubits.
+"""
+function pauli_basis_vectors(num_qubits::Int)
+    po = pauli_operators(num_qubits)
+    sop_dim = 4 ^ num_qubits
+    return mapreduce(x -> sparse(reshape(x.data, sop_dim)), (x, y) -> [x y], po)
+end
+
+"""
+    PauliTransferMatrix(sop::DenseSuperOperator)
+
+Convert a superoperator to its representation as a Pauli transfer matrix.
+"""
+function PauliTransferMatrix(sop::DenseSuperOperator{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    num_qubits = length(sop.basis_l[1].bases)
+    pbv = pauli_basis_vectors(num_qubits)
+    sop_dim = 4 ^ num_qubits
+    data = Matrix{Float64}(undef, (sop_dim, sop_dim))
+    data .= real.(pbv' * sop.data * pbv / √sop_dim)
+    return DensePauliTransferMatrix(sop.basis_l, sop.basis_r, data)
+end
+
+SuperOperator(unitary::DenseOperator{B, B, Matrix{ComplexF64}}) where B <: PauliBasis = spre(unitary) * spost(unitary')
+SuperOperator(sop::DenseSuperOperator{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis} = sop
+
+"""
+    SuperOperator(ptm::DensePauliTransferMatrix)
+
+Convert a Pauli transfer matrix to its representation as a superoperator.
+"""
+function SuperOperator(ptm::DensePauliTransferMatrix{B, B, Matrix{Float64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    num_qubits = length(ptm.basis_l[1].bases)
+    pbv = pauli_basis_vectors(num_qubits)
+    sop_dim = 4 ^ num_qubits
+    data = Matrix{ComplexF64}(undef, (sop_dim, sop_dim))
+    data .= pbv * ptm.data * pbv' / √sop_dim
+    return DenseSuperOperator(ptm.basis_l, ptm.basis_r, data)
+end
+
+"""
+    PauliTransferMatrix(unitary::DenseOperator)
+
+Convert an operator, presumably a unitary operator, to its representation as a
+Pauli transfer matrix.
+"""
+PauliTransferMatrix(unitary::DenseOperator{B, B, Matrix{ComplexF64}}) where B <: PauliBasis = PauliTransferMatrix(SuperOperator(unitary))
+
+"""
+    ChiMatrix(unitary::DenseOperator)
+
+Convert an operator, presumably a unitary operator, to its representation as a χ matrix.
+"""
+function ChiMatrix(unitary::DenseOperator{B, B, Matrix{ComplexF64}}) where B <: PauliBasis
+    num_qubits = length(unitary.basis_l.bases)
+    pbv = pauli_basis_vectors(num_qubits)
+    aj = pbv' * reshape(unitary.data, 4 ^ num_qubits)
+    return DenseChiMatrix((unitary.basis_l, unitary.basis_l), (unitary.basis_r, unitary.basis_r), aj * aj' / (2 ^ num_qubits))
+end
+
+"""
+    ChiMatrix(sop::DenseSuperOperator)
+
+Convert a superoperator to its representation as a Chi matrix.
+"""
+function ChiMatrix(sop::DenseSuperOperator{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    num_qubits = length(sop.basis_l)
+    sop_dim = 4 ^ num_qubits
+    po = pauli_operators(num_qubits)
+    data = Matrix{ComplexF64}(undef, (sop_dim, sop_dim))
+    for (idx, jdx) in Iterators.product(1:sop_dim, 1:sop_dim)
+        data[idx, jdx] = tr((spre(po[idx]) * spost(po[jdx])).data' * sop.data) / √sop_dim
+    end
+    return DenseChiMatrix(sop.basis_l, sop.basis_r, data)
+end
+
+"""
+    PauliTransferMatrix(chi_matrix::DenseChiMatrix)
+
+Convert a χ matrix to its representation as a Pauli transfer matrix.
+"""
+function PauliTransferMatrix(chi_matrix::DenseChiMatrix{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    num_qubits = length(chi_matrix.basis_l)
+    sop_dim = 4 ^ num_qubits
+    po = pauli_operators(num_qubits)
+    data = Matrix{Float64}(undef, (sop_dim, sop_dim))
+    for (idx, jdx) in Iterators.product(1:sop_dim, 1:sop_dim)
+        data[idx, jdx] = tr(mapreduce(x -> po[idx] * po[x[1]] * po[jdx] * po[x[2]] * chi_matrix.data[x[1], x[2]],
+                                      +,
+                                      Iterators.product(1:16, 1:16)).data) / sop_dim |> real
+    end
+    return DensePauliTransferMatrix(chi_matrix.basis_l, chi_matrix.basis_r, data)
+end
+
+"""
+    SuperOperator(chi_matrix::DenseChiMatrix)
+
+Convert a χ matrix to its representation as a superoperator.
+"""
+function SuperOperator(chi_matrix::DenseChiMatrix{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    return SuperOperator(PauliTransferMatrix(chi_matrix))
+end
+
+"""
+    ChiMatrix(ptm::DensePauliTransferMatrix)
+
+Convert a Pauli transfer matrix to its representation as a χ matrix.
+"""
+function ChiMatrix(ptm::DensePauliTransferMatrix{B, B, Matrix{Float64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    return ChiMatrix(SuperOperator(ptm))
+end
+
+"""Equality for all varieties of superoperators."""
+==(sop1::Union{DensePauliTransferMatrix, DenseSuperOperator, DenseChiMatrix},
+   sop2::Union{DensePauliTransferMatrix, DenseSuperOperator, DenseChiMatrix}) = ((typeof(sop1) == typeof(sop2)) &
+                                                                                 (sop1.basis_l == sop2.basis_l) &
+                                                                                 (sop1.basis_r == sop2.basis_r) &
+                                                                                  isapprox(sop1.data, sop2.data))
+
+end # end module

src/spin.jl

@@ -34,6 +34,18 @@ SpinBasis(spinnumber::Int) = SpinBasis(convert(Rational{Int}, spinnumber))
 
 ==(b1::SpinBasis, b2::SpinBasis) = b1.spinnumber==b2.spinnumber
 
+"""
+    identityoperator(b::SpinBasis)
+
+Pauli ``I`` operator for the given Spin basis.
+"""
+function identityoperator(b::SpinBasis)
+    N = length(b)
+    diag = ComplexF64[complex(1) for m=b.spinnumber:-1:-b.spinnumber]
+    data = spdiagm(0 => diag)
+    SparseOperator(b, data)
+end
+
 """
     sigmax(b::SpinBasis)
 

test/runtests.jl

@@ -46,6 +46,7 @@ names = [
     "test_stochastic_schroedinger.jl",
     "test_stochastic_master.jl",
     "test_stochastic_semiclassical.jl",
+    "test_pauli.jl",
 
     "test_printing.jl"
 ]

test/test_pauli.jl

@@ -0,0 +1,72 @@
+using LinearAlgebra
+using Test
+
+using QuantumOptics
+
+@testset "pauli" begin
+
+@test_throws MethodError PauliBasis(1.4)
+
+# Test conversion of unitary matrices to superoperators.
+q2 = PauliBasis(2)
+q3 = PauliBasis(3)
+CZ = DenseOperator(q2, q2, diagm(0 => [1,1,1,-1]))
+CZ_sop = SuperOperator(CZ)
+
+# Test conversion of unitary matrices to superoperators.
+@test diag(CZ_sop.data) ==  ComplexF64[1,1,1,-1,1,1,1,-1,1,1,1,-1,-1,-1,-1,1]
+@test CZ_sop.basis_l == CZ_sop.basis_r == (q2, q2)
+
+# Test conversion of superoperator to Pauli transfer matrix.
+CZ_ptm = PauliTransferMatrix(CZ_sop)
+
+# Test DensePauliTransferMatrix constructor.
+@test_throws DimensionMismatch DensePauliTransferMatrix((q2, q2), (q3, q3), CZ_ptm.data)
+@test DensePauliTransferMatrix((q2, q2), (q2, q2), CZ_ptm.data) == CZ_ptm
+
+@test all(isapprox.(CZ_ptm.data[[1,30,47,52,72,91,117,140,166,185,205,210,227,256]], 1))
+@test all(isapprox.(CZ_ptm.data[[106,151]], -1))
+
+@test CZ_ptm == PauliTransferMatrix(ChiMatrix(CZ))
+
+# Test construction of non-symmetric unitary.
+CNOT = DenseOperator(q2, q2, diagm(0 => [1,1,0,0], 1 => [0,0,1], -1 => [0,0,1]))
+CNOT_sop = SuperOperator(CNOT)
+CNOT_chi = ChiMatrix(CNOT)
+CNOT_ptm = PauliTransferMatrix(CNOT)
+
+@test CNOT_sop.basis_l == CNOT_sop.basis_r == (q2, q2)
+@test CNOT_chi.basis_l == CNOT_chi.basis_r == (q2, q2)
+@test CNOT_ptm.basis_l == CNOT_ptm.basis_r == (q2, q2)
+
+@test all(isapprox.(imag.(CNOT_sop.data), 0))
+@test all(isapprox.(imag.(CNOT_chi.data), 0))
+@test all(isapprox.(imag.(CNOT_ptm.data), 0))
+
+@test all(isapprox.(CNOT_sop.data[[1,18,36,51,69,86,104,119,141,158,176,191,201,218,236,251]], 1))
+@test all(isapprox.(CNOT_chi.data[[1,2,13,17,18,29,193,194,205,222]], 1))
+@test all(isapprox.(CNOT_chi.data[[14,30,206,209,210,221]], -1))
+@test all(isapprox.(CNOT_ptm.data[[1,18,47,64,70,85,108,138,153,183,205,222,227,244,]], 1))
+@test all(isapprox.(CNOT_ptm.data[[123,168]], -1))
+
+# Test DenseChiMatrix constructor.
+@test_throws DimensionMismatch DenseChiMatrix((q2, q2), (q3, q3), CNOT_chi.data)
+@test DenseChiMatrix((q2, q2), (q2, q2), CNOT_chi.data) == CNOT_chi
+
+# Test equality and conversion among all three bases.
+cphase = Complex{Float64}[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 exp(1im*.6)]
+
+CPHASE = DenseOperator(q2, cphase)
+
+CPHASE_sop = SuperOperator(CPHASE)
+CPHASE_chi = ChiMatrix(CPHASE)
+CPHASE_ptm = PauliTransferMatrix(CPHASE)
+
+@test ChiMatrix(CPHASE_sop) == CPHASE_chi
+@test ChiMatrix(CPHASE_ptm) == CPHASE_chi
+@test SuperOperator(CPHASE_chi) == CPHASE_sop
+@test SuperOperator(CPHASE_ptm) == CPHASE_sop
+@test PauliTransferMatrix(CPHASE_sop) == CPHASE_ptm
+@test PauliTransferMatrix(CPHASE_chi) == CPHASE_ptm
+
+end # testset