v0.7.5 update (#101)

* ✨ add `split_pauli_expression` 
* 🐛 also fix a bug of lamb shift test
* ➖ remove TensorOperations as dependency
* 🔖  bump version v0.7.5

Project.toml

@@ -1,7 +1,7 @@
 name = "OpenQuantumBase"
 uuid = "5dd19120-8766-11e9-1ef9-27094038a7db"
 authors = ["neversakura <[email protected]>"]
-version = "0.7.4"
+version = "0.7.5"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -14,7 +14,6 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
 [compat]
 Distributions = "0.22, 0.23, 0.24, 0.25"
@@ -25,7 +24,6 @@ QuadGK = "2.4"
 SpecialFunctions = "0.8, 0.9, 0.10, 1.0, 2.0"
 StaticArrays = "0.11, 0.12, 1.0"
 StatsBase = "0.30, 0.31, 0.32, 0.33"
-TensorOperations = "3.0"
 julia = "1.6"
 
 [extras]

src/OpenQuantumBase.jl

@@ -3,12 +3,11 @@ module OpenQuantumBase
 using DocStringExtensions
 
 import LinearAlgebra: mul!, lmul!, axpy!, ishermitian, Hermitian, eigen, eigen!,
-    tr, diag, Diagonal, norm, I, Bidiagonal
+    tr, diag, Diagonal, norm, I, Bidiagonal, qr, diagm
 import StaticArrays: SMatrix, MMatrix
 import SparseArrays: sparse, issparse, spzeros, SparseMatrixCSC
 import LinearAlgebra.BLAS: her!, gemm!
 import QuadGK: quadgk!, quadgk
-import TensorOperations: tensortrace
 
 """
 $(TYPEDEF)
@@ -116,6 +115,8 @@ export q_translate, standard_driver, local_field_term, two_local_term,
 export check_positivity, check_unitary, check_density_matrix, partial_trace,
     matrix_decompose, low_level_matrix, fidelity, inst_population, gibbs_state, purity, check_pure_state, find_degenerate
 
+export haar_unitary
+
 export construct_interpolations, gradient, log_uniform
 
 export hamiltonian_from_function, evaluate, issparse, get_cache, update_cache!,

src/base_util.jl

@@ -103,4 +103,29 @@ function bloch_to_state(θ::Real, ϕ::Real)
         throw(ArgumentError("ϕ is out of range 0≤ϕ≤2π."))
     end
     cos(θ / 2) * PauliVec[3][1] + exp(1.0im * ϕ) * sin(θ / 2) * PauliVec[3][2]
-end
+end
+
+"""
+$(SIGNATURES)
+
+Split a Pauli expression (a string like "-0.1X1X2+Z1") into a list of substrings
+that represent each Pauli clause. Each substring contains three parts: the sign;
+the prefix; and the Pauli string.
+
+# Example
+```julia-repl
+julia> split_pauli_expression("-0.1X1X2+Z2")
+2-element Vector{Any}:
+SubString{String}["-", "0.1", "X1X2"]
+SubString{String}["+", "", "Z2"]
+```
+"""
+function split_pauli_expression(p_str)
+    p_str = filter((x)->!isspace(x), p_str)
+    re_str = r"(\+|-|^)([0-9im.]*|^)([XYZ0-9]+)"
+    res = []
+    for m in eachmatch(re_str, p_str)
+        push!(res, [m[1], m[2], m[3]])
+    end
+    res
+end

src/math_util.jl

@@ -199,7 +199,9 @@ end
 """
 $(SIGNATURES)
 
-Calculate the partial trace of the density matrix `ρ`. Assume `ρ` is in the tensor space of ``ℋ₁⊗ℋ₂⊗…``, `sys_dim` is an array of the corresponding sub-system dimensions. `dim_2_keep` is an array of the indices whose corresponding subsystems are not traced out.
+Calculate the partial trace of the density matrix `ρ`. Assume `ρ` is in the tensor space of ``ℋ₁⊗ℋ₂⊗…``,
+`sys_dim` is an array of the corresponding sub-system dimensions. `dim_2_keep` is an array of the indices
+whose corresponding subsystems are not traced out.
 
 # Examples
 ```julia-repl
@@ -214,14 +216,24 @@ function partial_trace(ρ::Matrix, sys_dim::AbstractVector{Int}, dim_2_keep::Abs
     size(ρ, 1) != size(ρ, 2) && throw(ArgumentError("ρ is not a square matrix."))
     prod(sys_dim) != size(ρ, 1) && throw(ArgumentError("System dimensions do not multiply to density matrix dimension."))
 
-    N = length(sys_dim)
-    sys_dim = reverse(sys_dim)
-    ρ = reshape(ρ, sys_dim..., sys_dim...)
-    dim_2_keep = N .+ 1 .- dim_2_keep
-    dim_2_trace = [i for i in 1:N if !(i in dim_2_keep)]
-    traction_idx = collect(1:2*N)
-    traction_idx[dim_2_trace .+ N] .= traction_idx[dim_2_trace]
-    tensortrace(ρ, traction_idx)
+	N = length(sys_dim)
+	iden = CartesianIndex(ones(Int, 2*N)...)
+	sys_dim = reverse(sys_dim)
+	ρ = reshape(ρ, sys_dim..., sys_dim...)
+	dim_2_keep = N .+ 1 .- dim_2_keep
+	dim_2_trace = [i for i in 1:N if !(i in dim_2_keep)]
+	re_dim = copy(sys_dim)
+	re_dim[dim_2_trace] .= 1
+	tr_dim = copy(sys_dim)
+	tr_dim[dim_2_keep] .= 1
+	res = zeros(ComplexF64, re_dim..., re_dim...)
+	for I in CartesianIndices(size(res))
+		for k in CartesianIndices((tr_dim...,))
+			delta = CartesianIndex(k, k)
+			res[I] += ρ[I + delta - iden]
+		end
+	end
+	reshape(res, sys_dim[dim_2_keep]..., sys_dim[dim_2_keep]...)
 end
 
 """
@@ -316,3 +328,16 @@ Check whether the input `ρ` is a pure state: ``tr(ρ²)≈1``. This function wi
 function check_pure_state(ρ; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps())
     check_density_matrix(ρ) && isapprox(purity(ρ), 1, atol=atol, rtol=rtol)
 end
+
+"""
+$(SIGNATURES)
+
+Generate a Haar random matrix of dimension `(dim, dim)`.
+"""
+function haar_unitary(dim)
+   M = randn(dim,dim)+im*randn(dim, dim)
+   q,r = qr(M)
+   L = diag(r)
+   L=L./abs.(L)
+   q*diagm(0=>L)
+end

test/coupling_bath_interaction/bath.jl

@@ -77,5 +77,5 @@ lambfun_5 = OpenQuantumBase.build_lambshift([0.0, 0.01], true, cbath, Dict(:orde
 @test isapprox(lambfun_3[1,2](0.0), 0.03551, atol=1e-4)
 @test lambfun_2[1,1](0) == 0
 @test lambfun_3[2,2](0) == 0
-@test lambfun_4[1,2](0) ≈ 0.08681 atol=1e-4
-@test lambfun_5[2, 1](0.01) ≈ 0.12268 atol=1e-4
+@test lambfun_4[1,2](0) ≈ 0.03551 atol=1e-4
+@test lambfun_5[2,1](0.01) ≈ 0.05019 atol=1e-4

test/utilities/base_util.jl

@@ -15,4 +15,17 @@ using OpenQuantumBase, Test
 
 @test creation_operator(3) ≈ [0 0 0; 1 0 0; 0 sqrt(2) 0]
 @test annihilation_operator(3) ≈ [0 1 0; 0 0 sqrt(2); 0 0 0]
-@test number_operator(3) ≈ [0 0 0; 0 1 0; 0 0 2]
+@test number_operator(3) ≈ [0 0 0; 0 1 0; 0 0 2]
+
+pauli_exp = "-0.1X1X2 + Z2"
+res = OpenQuantumBase.split_pauli_expression(pauli_exp)
+@test res[1] == ["-", "0.1", "X1X2"]
+@test res[2] == ["+", "", "Z2"]
+pauli_exp = "Y2X1-2Z1"
+res = OpenQuantumBase.split_pauli_expression(pauli_exp)
+@test res[1] == ["", "", "Y2X1"]
+@test res[2] == ["-", "2", "Z1"]
+pauli_exp = "X1X2 + 1.0imZ1"
+res = OpenQuantumBase.split_pauli_expression(pauli_exp)
+@test res[2] == ["+", "1.0im", "Z1"]
+

test/utilities/math_util.jl

@@ -1,4 +1,4 @@
-using OpenQuantumBase, Test
+using OpenQuantumBase, LinearAlgebra, Test
 
 # === matrix decomposition
 v = 1.0 * σx + 2.0 * σy + 3.0 * σz
@@ -90,4 +90,7 @@ gibbs = gibbs_state(σz, 12)
     0.0+0.0im -0.9+0.0im 0.0+0.0im 0.0+0.0im
     0.0+0.0im 0.0+0.0im -1.1+0.0im 0.0+0.0im
     0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im]
-@test_logs (:warn, "Subspace dimension bigger than total dimension.") low_level_matrix(σz ⊗ σz + 0.1σz ⊗ σi, 5)
+@test_logs (:warn, "Subspace dimension bigger than total dimension.") low_level_matrix(σz ⊗ σz + 0.1σz ⊗ σi, 5)
+
+u = haar_unitary(2)
+@test u'*u ≈ I