Merge 1dceb50 into a9df934

docs/source/api.rst

@@ -26,6 +26,8 @@ Composite bases
 
 .. jl:autofunction:: bases.jl ptrace
 
+.. jl:autofunction:: bases.jl permutebases
+
 
 Subspace bases
 ^^^^^^^^^^^^^^
@@ -122,6 +124,8 @@ DenseOperators
 
 .. jl:autofunction:: operators.jl dense_identityoperator
 
+.. jl:autofunction:: operators.jl permutesystems
+
 
 .. _section-api-sparseoperators:
 
@@ -178,6 +182,8 @@ Metrics
 
 .. jl:autofunction:: metrics.jl tracedistance
 
+.. jl:autofunction:: metrics.jl tracedistance_general
+
 
 
 Systems

src/QuantumOptics.jl

@@ -3,7 +3,7 @@ __precompile__()
 module QuantumOptics
 
 export bases, Basis, GenericBasis, CompositeBasis,
-        tensor, ⊗,
+        tensor, ⊗, permutesystems,
         states, StateVector, Bra, Ket, basis_bra, basis_ket,
                 dagger, normalize, normalize!,
         operators, Operator, DenseOperator, projector,
@@ -22,7 +22,7 @@ export bases, Basis, GenericBasis, CompositeBasis,
                 positionoperator, momentumoperator, laplace_x, laplace_p,
         nlevel, NLevelBasis, transition, nlevelstate,
         nparticlebasis, BosonicNParticleBasis, FermionicNParticleBasis, nparticleoperator_1, nparticleoperator_2,
-        metrics, tracedistance,
+        metrics, tracedistance, tracedistance_general,
         spectralanalysis, operatorspectrum, operatorspectrum_hermitian,
                 eigenstates, eigenstates_hermitian, groundstate,
         timeevolution_simple,

src/bases.jl

@@ -3,7 +3,7 @@ module bases
 import Base.==
 
 export Basis, GenericBasis, CompositeBasis,
-       tensor, ⊗, dualbasis, ptrace,
+       tensor, ⊗, dualbasis, ptrace, permutesystems,
        equal_shape, equal_bases, multiplicable,
        IncompatibleBases,
        check_equal, check_multiplicable
@@ -187,4 +187,27 @@ function ptrace(b::CompositeBasis, indices::Vector{Int})
 end
 
 
+function permutesystems(bases::Vector{Basis}, perm::Vector{Int})
+    @assert length(bases) == length(perm)
+    @assert issubset(Set(1:length(bases)), Set(perm))
+    bases[perm]
+end
+
+"""
+Change the ordering of the subbases in a CompositeBasis.
+
+For a permutation vector [2,1,3] and a given composite basis [b1, b2, b3]
+this function results in [b2,b1,b3].
+
+Arguments
+---------
+
+basis
+    A composite basis.
+perm
+    Vector defining the new ordering of the sub-bases.
+"""
+bases.permutesystems(basis::CompositeBasis, perm::Vector{Int}) = CompositeBasis(permutesystems(basis.bases, perm)...)
+
+
 end # module

src/metrics.jl

@@ -2,19 +2,47 @@ module metrics
 
 using ..operators, ..operators_sparse
 
-export tracedistance
+export tracedistance, tracedistance_general
 
+"""
+Trace distance between two density operators.
 
+It uses the identity
+
+.. math:
+
+    T(\\rho, \\sigma) = \\frac{1}{2} \\sum_i |\\lambda_i|
+
+where :math:`\\lambda_i` are the eigenvalues of the matrix
+:math:`\\rho - \\sigma`. This works only if :math:`rho` and :math:`sigma`
+are density operators. For trace distances between general operators use
+:jl:func:`tracedistance_general`.
+"""
 function tracedistance(rho::DenseOperator, sigma::DenseOperator)
-    delta = (rho - sigma).data
-    @assert size(delta, 1) == size(delta, 2)
-    for i=1:size(delta,1)
-        delta[i,i] = real(delta[i,i])
+    delta = (rho - sigma)
+    @assert length(delta.basis_l) == length(delta.basis_r)
+    data = delta.data
+    for i=1:size(data,1)
+        data[i,i] = real(data[i,i])
     end
-    s = eigvals(Hermitian(delta))
+    s = eigvals(Hermitian(data))
     return 0.5*sum(abs(s))
 end
 
+
+"""
+Trace distance between two operators.
+
+.. math:
+
+    T(\\rho, \\sigma) = \\frac{1}{2}
+            Tr\\{\\ \\sqrt{(\\rho-\\sigma)^\\dagger(\\rho-\\sigma)}\\}
+"""
+function tracedistance_general(rho::DenseOperator, sigma::DenseOperator)
+    delta = (rho - sigma)
+    return 0.5*trace(sqrtm((dagger(delta)*delta).data))
+end
+
 function tracedistance{T<:Operator}(rho::T, sigma::T)
     throw(ArgumentError("tracedistance not implemented for $(T). Use dense operators instead."))
 end

src/operators.jl

@@ -283,5 +283,24 @@ Partial trace of the given state vector over the specified indices.
 bases.ptrace(a::Ket, indices) = bases.ptrace(tensor(a, dagger(a)), indices)
 bases.ptrace(a::Bra, indices) = bases.ptrace(tensor(dagger(a), a), indices)
 
+"""
+Change the ordering of the subsystems of the given operator.
+
+Arguments
+---------
+rho
+    A dense operator represented in a composite basis.
+perm
+    Vector defining the new ordering of the subsystems.
+"""
+function bases.permutesystems(rho::DenseOperator, perm::Vector{Int})
+    @assert length(rho.basis_l.bases) == length(rho.basis_r.bases) == length(perm)
+    @assert issubset(Set(1:length(rho.basis_l.bases)), Set(perm))
+    data = reshape(rho.data, [reverse(rho.basis_l.shape); reverse(rho.basis_r.shape)]...)
+    dataperm = length(perm) - reverse(perm) + 1
+    data = permutedims(data, [dataperm; dataperm + length(perm)])
+    data = reshape(data, length(rho.basis_l), length(rho.basis_r))
+    DenseOperator(permutesystems(rho.basis_l, perm), permutesystems(rho.basis_r, perm), data)
+end
 
 end # module

src/operators_lazy.jl

@@ -38,8 +38,8 @@ type LazyTensor <: LazyOperator
     operators::Dict{Int,Operator}
 
     function LazyTensor(basis_l::CompositeBasis, basis_r::CompositeBasis, operators::Dict{Int,Operator}, factor::Number=1.)
-        N = length(basis_l)
-        @assert N==length(basis_r)
+        N = length(basis_l.bases)
+        @assert N==length(basis_r.bases)
         @assert maximum(keys(operators))<=N
         @assert 1<=minimum(keys(operators))
         for (i,op) = operators
@@ -333,4 +333,14 @@ function operators.gemv!(alpha, a::Bra, b::LazyProduct, beta, result::Bra)
     operators.gemv!(alpha, tmp1, b.operators[end], beta, result)
 end
 
+function operators.permutesystems(op::LazyTensor, perm::Vector{Int})
+    b_l = permutesystems(op.basis_l, perm)
+    b_r = permutesystems(op.basis_r, perm)
+    operators = Dict{Int,Operator}(findfirst(perm, i)=>op_i for (i, op_i) in op.operators)
+    LazyTensor(b_l, b_r, operators, op.factor)
+end
+
+operators.permutesystems(op::LazySum, perm::Vector{Int}) = LazySum(op.factors, Operator[permutesystems(op_i, perm) for op_i in op.operators])
+operators.permutesystems(op::LazyProduct, perm::Vector{Int}) = LazyProduct(op.factor, Operator[permutesystems(op_i, perm) for op_i in op.operators])
+
 end

src/operators_sparse.jl

@@ -98,4 +98,24 @@ function diagonaloperator{T <: Number}(b::Basis, diag::Vector{T})
   SparseOperator(b, spdiagm(complex(float(diag))))
 end
 
+
+"""
+Change the ordering of the subsystems of the given operator.
+
+Arguments
+---------
+rho
+    A sparse operator represented in a composite basis.
+perm
+    Vector defining the new ordering of the subsystems.
+"""
+function operators.permutesystems(rho::SparseOperator, perm::Vector{Int})
+    @assert length(rho.basis_l.bases) == length(rho.basis_r.bases) == length(perm)
+    @assert issubset(Set(1:length(rho.basis_l.bases)), Set(perm))
+    shape = [reverse(rho.basis_l.shape); reverse(rho.basis_r.shape)]
+    dataperm = length(perm) - reverse(perm) + 1
+    data = sparsematrix.permutedims(rho.data, shape, [dataperm; dataperm + length(perm)])
+    SparseOperator(permutesystems(rho.basis_l, perm), permutesystems(rho.basis_r, perm), data)
+end
+
 end # module

src/sparsematrix.jl

@@ -59,5 +59,19 @@ function gemv!(alpha::Complex128, v::Vector{Complex128}, M::SparseMatrixCSC{Comp
     nothing
 end
 
+function permutedims(x, shape, perm)
+    shape = (shape...)
+    shape_perm = ([shape[i] for i in perm]...)
+    y = spzeros(Complex128, x.m, x.n)
+    for I in eachindex(x)
+        linear_index = sub2ind((x.m, x.n), I.I...)
+        cartesian_index = ind2sub(shape, linear_index)
+        cartesian_index_perm = [cartesian_index[i] for i=perm]
+        linear_index_perm = sub2ind(shape_perm, cartesian_index_perm...)
+        J = ind2sub((x.m, x.n), linear_index_perm)
+        y[J...] = x[I.I...]
+    end
+    y
+end
 
 end # module

src/states.jl

@@ -102,4 +102,25 @@ basis_ket(b::Basis, indices::Array{Int}) = Ket(b, basis_vector(b.shape, indices)
 basis_ket(b::Basis, index::Int) = basis_ket(b, [index])
 
 
+"""
+Change the ordering of the subsystems of the given state.
+
+Arguments
+---------
+state
+    A state represented in a composite basis.
+perm
+    Vector defining the new ordering of the subsystems.
+"""
+function bases.permutesystems{T<:StateVector}(state::T, perm::Vector{Int})
+    @assert length(state.basis.bases) == length(perm)
+    @assert issubset(Set(1:length(state.basis.bases)), Set(perm))
+    data = reshape(state.data, reverse(state.basis.shape)...)
+    dataperm = length(perm) - reverse(perm) + 1
+    data = permutedims(data, dataperm)
+    data = reshape(data, length(data))
+    T(permutesystems(state.basis, perm), data)
+end
+
+
 end # module

test/test_bases.jl

@@ -3,9 +3,11 @@ using QuantumOptics
 
 shape1 = [5]
 shape2 = [2, 3]
+shape3 = [6]
 
 b1 = GenericBasis(shape1)
 b2 = GenericBasis(shape2)
+b3 = GenericBasis(shape3)
 
 @test b1.shape == shape1
 @test b2.shape == shape2
@@ -27,3 +29,7 @@ comp_b1_b2 = tensor(comp_b1, comp_b2)
 @test ptrace(comp_b2, [1]) == ptrace(comp_b2, [2]) == comp_b1
 @test ptrace(comp_b2, [1, 2]) == ptrace(comp_b1, [1])
 @test ptrace(comp_b2, [2, 3]) == ptrace(comp_b1, [2])
+
+comp1 = tensor(b1, b2, b3)
+comp2 = tensor(b2, b1, b3)
+@test permutesystems(comp1, [2,1,3]) == comp2

test/test_lazyoperators.jl

@@ -129,3 +129,47 @@ operators.gemv!(Complex(1.), op, psi_ket, Complex(0.), result_ket)
 @test_approx_eq_eps 0. norm(full(op)*psi_ket - result_ket) 1e-12
 operators.gemv!(Complex(1.), psi_bra, op, Complex(0.), result_bra)
 @test_approx_eq_eps 0. norm(psi_bra*full(op) - result_bra) 1e-12
+
+
+# Test permutating systems
+b1a = NLevelBasis(2)
+b1b = SpinBasis(3//2)
+b2a = SpinBasis(1//2)
+b2b = FockBasis(7)
+b3a = FockBasis(2)
+b3b = NLevelBasis(4)
+
+b_l = b1a⊗b2a⊗b3a
+b_r = b1b⊗b2b⊗b3b
+
+srand(0)
+rho1 = DenseOperator(b1a, b1b, rand(Complex128, length(b1a), length(b1b)))
+rho2 = DenseOperator(b2a, b2b, rand(Complex128, length(b2a), length(b2b)))
+rho3 = DenseOperator(b3a, b3b, rand(Complex128, length(b3a), length(b3b)))
+
+rho123 = LazyTensor(b_l, b_r, [1,2,3], [rho1,rho2,rho3])
+rho213_ = permutesystems(rho123, [2, 1, 3])
+
+@test_approx_eq_eps 0. tracedistance_general(full(rho213_), rho2⊗rho1⊗rho3) 1e-5
+
+
+rho1 = DenseOperator(b_l, b_r, rand(Complex128, length(b_l), length(b_r)))
+rho2 = DenseOperator(b_r, b_l, rand(Complex128, length(b_r), length(b_l)))
+rho3 = DenseOperator(b_l, b_r, rand(Complex128, length(b_l), length(b_r)))
+
+rho123 = LazyProduct(rho1, rho2, rho3)
+rho213_ = permutesystems(rho123, [2, 1, 3])
+rho213 = permutesystems(rho1*rho2*rho3, [2, 1, 3])
+
+@test_approx_eq_eps 0. tracedistance_general(full(rho213_), rho213) 1e-5
+
+
+rho1 = DenseOperator(b_l, b_r, rand(Complex128, length(b_l), length(b_r)))
+rho2 = DenseOperator(b_l, b_r, rand(Complex128, length(b_l), length(b_r)))
+rho3 = DenseOperator(b_l, b_r, rand(Complex128, length(b_l), length(b_r)))
+
+rho123 = LazySum(rho1, rho2, rho3)
+rho213_ = permutesystems(rho123, [2, 1, 3])
+rho213 = permutesystems(rho1+rho2+rho3, [2, 1, 3])
+
+@test_approx_eq_eps 0. tracedistance_general(full(rho213_), rho213) 1e-5

test/test_metrics.jl

@@ -1,14 +1,23 @@
 using Base.Test
 using QuantumOptics
 
-b = SpinBasis(1//2)
+b1 = SpinBasis(1//2)
+b2 = FockBasis(6)
 
-psi1 = spinup(b)
-psi2 = spindown(b)
+psi1 = spinup(b1) ⊗ coherentstate(b2, 0.1)
+psi2 = spindown(b1) ⊗ fockstate(b2, 2)
 
 rho = tensor(psi1, dagger(psi1))
 sigma = tensor(psi2, dagger(psi2))
 
 @test tracedistance(rho, sigma) == 1.
 @test tracedistance(rho, rho) == 0.
 @test tracedistance(sigma, sigma) == 0.
+
+@test_approx_eq_eps tracedistance_general(rho, sigma) 1. 1e-6
+@test tracedistance_general(rho, rho) == 0.
+@test tracedistance_general(sigma, sigma) == 0.
+
+rho = spinup(b1) ⊗ dagger(coherentstate(b2, 0.1))
+@test_throws AssertionError tracedistance(rho, rho)
+@test tracedistance_general(rho, rho) == 0.

test/test_operators.jl

@@ -93,3 +93,18 @@ operators.gemv!(complex(1.0), at, xket, complex(0.), result_ket)
 result_bra = deepcopy(xbra)
 operators.gemv!(complex(1.0), xbra, at, complex(0.), result_bra)
 @test_approx_eq 0. norm(result_bra-xbra*at)
+
+# Test permutating systems
+b1a = NLevelBasis(2)
+b1b = SpinBasis(3//2)
+b2a = SpinBasis(1//2)
+b2b = FockBasis(7)
+b3a = FockBasis(2)
+b3b = NLevelBasis(4)
+
+srand(0)
+rho1 = DenseOperator(b1a, b1b, rand(Complex128, length(b1a), length(b1b)))
+rho2 = DenseOperator(b2a, b2b, rand(Complex128, length(b2a), length(b2b)))
+rho3 = DenseOperator(b3a, b3b, rand(Complex128, length(b3a), length(b3b)))
+
+@test_approx_eq_eps 0. tracedistance_general(permutesystems(rho1⊗rho2⊗rho3, [2, 1, 3]), rho2⊗rho1⊗rho3) 1e-5

test/test_sparseoperators.jl

@@ -37,8 +37,22 @@ a = A()
 
 @test_throws ArgumentError sparse(a)
 
-
 @test diagonaloperator(b, [1, 1, 1, 1]) == I
 @test diagonaloperator(b, [1., 1., 1., 1.]) == I
 @test diagonaloperator(b, [1im, 1im, 1im, 1im]) == 1im*I
 @test diagonaloperator(b, [0:3;]) == number(b)
+
+# Test permutating systems
+b1a = NLevelBasis(2)
+b1b = SpinBasis(3//2)
+b2a = SpinBasis(1//2)
+b2b = FockBasis(7)
+b3a = FockBasis(2)
+b3b = NLevelBasis(4)
+
+srand(0)
+rho1 = SparseOperator(b1a, b1b, sparse(rand(Complex128, length(b1a), length(b1b))))
+rho2 = SparseOperator(b2a, b2b, sparse(rand(Complex128, length(b2a), length(b2b))))
+rho3 = SparseOperator(b3a, b3b, sparse(rand(Complex128, length(b3a), length(b3b))))
+
+@test_approx_eq_eps 0. tracedistance_general(full(permutesystems(rho1⊗rho2⊗rho3, [2, 1, 3])), full(rho2⊗rho1⊗rho3)) 1e-5