join files together

src/Ket.jl

@@ -15,8 +15,7 @@ include("games.jl")
 include("measurements.jl")
 include("entropy.jl")
 include("norms.jl")
-include("partial_tra.jl")
-include("permute_sys.jl")
+include("multilinear.jl")
 include("supermaps.jl")
 include("entanglement.jl")
 

src/partial_tra.jl → src/multilinear.jl

@@ -184,3 +184,55 @@ Takes the partial transpose of matrix `X` with subsystem dimensions `dims` on th
 """
 partial_transpose(X::AbstractMatrix, transp::Integer, dims::Vector{<:Integer}) = partial_transpose(X, [transp], dims)
 export partial_transpose
+
+"""
+    permute_systems(X::AbstractVector, perm::Vector, dims::Vector)
+
+Permutes the order of the subsystems of vector `X` with subsystem dimensions `dims` according to the permutation `perm`.
+"""
+function permute_systems(X::AbstractVector{T}, perm::Vector{<:Integer}, dims::Vector{<:Integer}) where {T}
+    perm == 1:length(perm) && return X
+
+    dX = length(X)
+    dimsY = dims[perm]
+
+    Y = Vector{T}(undef, dX)
+
+    ti = Vector{Int64}(undef, length(dims))
+
+    for i in 1:dX
+        _tidx!(ti, i, dims)
+        permute!(ti, perm)
+        Yi = _idx(ti, dimsY)
+        Y[Yi] = X[i]
+    end
+    return Y
+end
+
+"""
+    permute_systems(X::AbstractMatrix, perm::Vector, dims::Vector)
+
+Permutes the order of the subsystems of the square matrix `X`, which is composed by square subsystems of dimensions `dims`, according to the permutation `perm`.
+"""
+function permute_systems(X::AbstractMatrix, perm::Vector{<:Integer}, dims::Vector{<:Integer})
+    perm == 1:length(perm) && return X
+
+    idxperm = permute_systems(axes(X, 1), perm, dims)
+    return X[idxperm, idxperm]
+end
+
+"""
+    permute_systems(X::AbstractMatrix, perm::Vector, dims::Matrix)
+
+Permutes the order of the subsystems of the matrix `X`, which is composed by subsystems of dimensions `dims`, according to the permutation `perm`.
+`dims` should be a n x 2 matrix where `dims[i, 1]` is the number of rows of subsystem i, and `dims[i,2]` is its number of columns. 
+"""
+function permute_systems(X::AbstractMatrix, perm::Vector{<:Integer}, dims::Matrix{<:Integer})
+    perm == 1:length(perm) && return X
+
+    rowperm = permute_systems(axes(X, 1), perm, dims[:, 1])
+    colperm = permute_systems(axes(X, 2), perm, dims[:, 2])
+
+    return X[rowperm, colperm]
+end
+export permute_systems

test/permute_sys.jl → test/multilinear.jl

@@ -1,3 +1,70 @@
+@testset "Multilinear algebra" begin
+@testset "Partial trace      " begin
+    d1, d2, d3 = 2, 3, 4
+    for R in [Float64, Double64, Float128, BigFloat]
+        for T in [R, Complex{R}]
+            a = randn(T, d1, d1)
+            b = randn(T, d2, d2)
+            c = randn(T, d3, d3)
+            ab = kron(a, b)
+            ac = kron(a, c)
+            bc = kron(b, c)
+            abc = kron(ab, c)
+            @test partial_trace(ab, [1, 2], [d1, d2])[1] ≈ tr(ab)
+            @test partial_trace(ab, 2, [d1, d2]) ≈ a * tr(b)
+            @test partial_trace(ab, 1, [d1, d2]) ≈ b * tr(a)
+            @test partial_trace(ab, Int64[], [d1, d2]) ≈ ab
+            @test partial_trace(abc, [1, 2, 3], [d1, d2, d3])[1] ≈ tr(abc)
+            @test partial_trace(abc, [2, 3], [d1, d2, d3]) ≈ a * tr(b) * tr(c)
+            @test partial_trace(abc, [1, 3], [d1, d2, d3]) ≈ b * tr(a) * tr(c)
+            @test partial_trace(abc, [1, 2], [d1, d2, d3]) ≈ c * tr(a) * tr(b)
+            @test partial_trace(abc, 3, [d1, d2, d3]) ≈ ab * tr(c)
+            @test partial_trace(abc, 2, [d1, d2, d3]) ≈ ac * tr(b)
+            @test partial_trace(abc, 1, [d1, d2, d3]) ≈ bc * tr(a)
+            @test partial_trace(abc, Int64[], [d1, d2, d3]) ≈ abc
+        end
+    end
+    for wrapper in [Symmetric, Hermitian]
+        M = wrapper(randn(ComplexF64, (d1 * d2 * d3, d1 * d2 * d3)))
+        x = Matrix(M)
+        @test partial_trace(M, 2, [d1, d2, d3]) ≈ partial_trace(x, 2, [d1, d2, d3])
+        @test partial_trace(M, [1, 3], [d1, d2, d3]) ≈ partial_trace(x, [1, 3], [d1, d2, d3])
+    end
+end
+
+@testset "Partial transpose  " begin
+    d1, d2, d3 = 2, 3, 4
+    for R in [Float64, Double64, Float128, BigFloat]
+        for T in [R, Complex{R}]
+            a = randn(T, d1, d1)
+            b = randn(T, d2, d2)
+            c = randn(T, d3, d3)
+            ab = kron(a, b)
+            ac = kron(a, c)
+            bc = kron(b, c)
+            abc = kron(ab, c)
+            @test partial_transpose(ab, [1, 2], [d1, d2]) ≈ transpose(ab)
+            @test partial_transpose(ab, 2, [d1, d2]) ≈ kron(a, transpose(b))
+            @test partial_transpose(ab, 1, [d1, d2]) ≈ kron(transpose(a), b)
+            @test partial_transpose(ab, Int64[], [d1, d2]) ≈ ab
+            @test partial_transpose(abc, [1, 2, 3], [d1, d2, d3]) ≈ transpose(abc)
+            @test partial_transpose(abc, [2, 3], [d1, d2, d3]) ≈ kron(a, transpose(b), transpose(c))
+            @test partial_transpose(abc, [1, 3], [d1, d2, d3]) ≈ kron(transpose(a), b, transpose(c))
+            @test partial_transpose(abc, [1, 2], [d1, d2, d3]) ≈ kron(transpose(a), transpose(b), c)
+            @test partial_transpose(abc, 3, [d1, d2, d3]) ≈ kron(ab, transpose(c))
+            @test partial_transpose(abc, 2, [d1, d2, d3]) ≈ kron(a, transpose(b), c)
+            @test partial_transpose(abc, 1, [d1, d2, d3]) ≈ kron(transpose(a), bc)
+            @test partial_transpose(abc, Int64[], [d1, d2, d3]) ≈ abc
+        end
+    end
+    for wrapper in [Symmetric, Hermitian]
+        M = wrapper(randn(ComplexF64, (d1 * d2 * d3, d1 * d2 * d3)))
+        x = Matrix(M)
+        @test partial_transpose(M, 2, [d1, d2, d3]) ≈ partial_transpose(x, 2, [d1, d2, d3])
+        @test partial_transpose(M, [1, 3], [d1, d2, d3]) ≈ partial_transpose(x, [1, 3], [d1, d2, d3])
+    end
+end
+
 @testset "Permute systems    " begin
     @testset "Vectors" begin
         d1, d2, d3 = 2, 3, 4
@@ -74,3 +141,6 @@
         end
     end
 end
+end
+
+#TODO add test with JuMP variables

test/runtests.jl

@@ -11,8 +11,7 @@ include("entanglement.jl")
 include("measurements.jl")
 include("nonlocal.jl")
 include("norms.jl")
-include("partial_tra.jl")
-include("permute_sys.jl")
+include("multilinear.jl")
 include("random.jl")
 include("states.jl")
 include("supermaps.jl")