Implement "permutesystems" to change the ordering of subsystems.

* Implemented for dense, sparse and lazy operators.
* docs and tests.

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:
 

src/QuantumOptics.jl

@@ -3,12 +3,12 @@ __precompile__()
 module QuantumOptics
 
 export bases, Basis, GenericBasis, CompositeBasis,
-        tensor, ⊗,
+        tensor, ⊗, permutebases,
         states, StateVector, Bra, Ket, basis_bra, basis_ket,
                 dagger, normalize, normalize!,
         operators, Operator, DenseOperator, projector,
                 expect, identityoperator, dense_identityoperator,
-                ptrace, embed,
+                ptrace, embed, permutesystems,
         operators_lazy, LazyOperator, LazyTensor, LazySum, LazyProduct,
         operators_sparse, SparseOperator,
                 sparse_identityoperator, diagonaloperator,

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/operators.jl

@@ -11,7 +11,7 @@ using ..bases, ..states
 export Operator, DenseOperator,
        tensor, dagger, projector,
        expect, identityoperator, dense_identityoperator,
-       embed
+       embed, permutesystems
 
 
 """
@@ -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_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

test/test_states.jl

@@ -54,3 +54,21 @@ c = normalize(Bra(CompositeBasis(basis, basis), [1im, 2im, 0, 0]))
 @test_approx_eq 0. norm(tensor(a, b) - c)
 @test_approx_eq_eps 0. tracedistance(operators.ptrace(c, 1), tensor(dagger(b), b)) 1e-15
 @test_approx_eq_eps 0. tracedistance(operators.ptrace(c, 2), tensor(dagger(a), a)) 1e-15
+
+
+# Test permutating systems
+b1 = NLevelBasis(2)
+b2 = SpinBasis(1//2)
+b3 = FockBasis(2)
+
+srand(0)
+psi1 = normalize(Ket(b1, rand(Complex128, length(b1))))
+psi2 = normalize(Ket(b2, rand(Complex128, length(b2))))
+psi3 = normalize(Ket(b3, rand(Complex128, length(b3))))
+
+psi123 = psi1 ⊗ psi2 ⊗ psi3
+psi213 = psi2 ⊗ psi1 ⊗ psi3
+
+c = dagger(psi213)*permutesystems(psi123, [2,1,3])
+
+@test_approx_eq_eps 1. c 1e-5