From ec6c98a2499e5d766fe5e388c4a8b091e03028e7 Mon Sep 17 00:00:00 2001
From: Yuning <neuromancer43@gmail.com>
Date: Tue, 2 Jul 2024 17:43:05 +0200
Subject: [PATCH] add dephasing matrix

---
 src/SpinShuttling.jl        | 41 ++++++++++++++++++++++++++++++++++-
 src/stochastics.jl          | 31 ++++++++++++++++++++++++--
 test/runtests.jl            | 11 +++++-----
 test/testdephasingfactor.jl | 43 +++++++++++++++++++++++++++++++++++++
 4 files changed, 118 insertions(+), 8 deletions(-)
 create mode 100644 test/testdephasingfactor.jl

diff --git a/src/SpinShuttling.jl b/src/SpinShuttling.jl
index fe4c3ed..39e3dbe 100644
--- a/src/SpinShuttling.jl
+++ b/src/SpinShuttling.jl
@@ -15,7 +15,7 @@ export ShuttlingModel, OneSpinModel, TwoSpinModel,
 OneSpinForthBackModel, TwoSpinParallelModel, RandomFunction, CompositeRandomFunction,
 OrnsteinUhlenbeckField, PinkBrownianField
 export averagefidelity, fidelity, sampling, characteristicfunction, characteristicvalue
-export covariance, covariancematrix
+export dephasingmatrix, covariance, covariancematrix
 export W
 
 """
@@ -276,6 +276,45 @@ function sampling(model::ShuttlingModel, objective::Function, M::Int; vector::Bo
     return sampling(samplingfunction, M)
 end
 
+"""
+
+# Arguments
+- `i::Int`: index of qubits, range from (1,n)
+- `p::Int`: index of spin state, range from (1,2^n)
+- `n::Int`: number of spins
+"""
+function m(i::Int,p::Int,n::Int)
+    1/2-digits(p, base=2, pad=n)[i]
+end
+
+"""
+Calculate the dephasing matrix of a given spin shuttling model.
+"""
+function dephasingmatrix(model::ShuttlingModel)::Symmetric{<:Real}
+    n=model.n
+    W=zeros(2^n,2^n)
+    for j in 1:2^n
+        W[j,j]=1
+        for k in 1:j-1
+            c=[trunc(Int,m(i,j-1,n)-m(i,k-1,n)) for i in 1:n]
+            R = CompositeRandomFunction(model.R, c)
+            W[j,k] = characteristicvalue(R)
+            W[k,j] = W[j,k]
+        end
+    end
+    return Symmetric(W)
+end
+
+function dephasingcoeffs(n::Int)::Array{Real,3}
+    M=zeros(2^n,2^n, n)
+    for j in 1:2^n
+        for k in 1:2^n
+            c=[m(i,j-1,n)-m(i,k-1,n) for i in 1:n]
+            M[j,k, :] = c
+        end
+    end
+    return M
+end
 
 """
 Sample a phase integral of the process. 
diff --git a/src/stochastics.jl b/src/stochastics.jl
index 4fcdf1a..755b8e3 100644
--- a/src/stochastics.jl
+++ b/src/stochastics.jl
@@ -14,6 +14,7 @@ struct OrnsteinUhlenbeckField <: GaussianRandomField
     σ::Real # covariance
 end
 
+
 """
 Pink-Brownian Field, the correlation function of which is
 `σ^2 * (expinti(-γ[2]abs(t₁ - t₂)) - expinti(-γ[1]abs(t₁ - t₂)))/log(γ[2]/γ[1]) * exp(-|x₁-x₂|/θ)`
@@ -56,6 +57,13 @@ end
 
 """
 Divide the covariance matrix of a direct summed random function into partitions. 
+
+# Arguments
+- `R::RandomFunction`: a direct sum of random processes R₁⊕ R₂⊕ ... ⊕ Rₙ
+- `n::Int`: number of partitions or spins
+
+# Returns
+- `Matrix{Matrix{<:Real}}`: a matrix of covariance matrices
 """
 function covariancepartition(R::RandomFunction, n::Int)::Matrix{Matrix{<:Real}}
     Λ=Matrix{Matrix{<:Real}}(undef, n, n)
@@ -69,7 +77,15 @@ function covariancepartition(R::RandomFunction, n::Int)::Matrix{Matrix{<:Real}}
     return Λ
 end
 
+"""
+
+# Arguments
+- `R::RandomFunction`: a direct sum of random processes R₁⊕ R₂⊕ ... ⊕ Rₙ
+- `n::Int`: number of partitions or spins
 
+# Returns
+- `Vector{Vector{<:Real}}`: a vector of mean vectors
+"""
 function meanpartition(R::RandomFunction, n::Int)::Vector{Vector{<:Real}}
     N=length(R.P)÷n
     return [R.μ[(i-1)*N+1: i*N] for i in 1:n]
@@ -79,6 +95,13 @@ end
 Create a new random function composed by a linear combination of random processes.
 The input random function represents the direct sum of these processes. 
 The output random function is a tensor contraction from the input.
+
+# Arguments
+- `R::RandomFunction`: a direct sum of random processes R₁⊕ R₂⊕ ... ⊕ Rₙ
+- `c::Vector{Int}`: a vector of coefficients
+
+# Returns
+- `RandomFunction`: a new random function composed by a linear combination of random processes
 """
 function CompositeRandomFunction(R::RandomFunction, c::Vector{Int})::RandomFunction
     n=length(c)
@@ -89,19 +112,23 @@ function CompositeRandomFunction(R::RandomFunction, c::Vector{Int})::RandomFunct
     return RandomFunction(μ, t, Σ, cholesky(Σ))
 end
 
+
 function CompositeRandomFunction(P::Vector{<:Point}, process::GaussianRandomField, c::Vector{Int})::RandomFunction
     return CompositeRandomFunction(RandomFunction(P, process), c)
 end
 
 """
 Generate a random time series from a Gaussian random field.
+
+`R()` generates a random time series from a Gaussian random field `R`
+`R(randseq)` generates a random time series from a Gaussian random field `R` with a given random sequence `randseq`.
 """
 function (R::RandomFunction)(randseq::Vector{<:Real})
     return R.μ .+ R.C.L*randseq
 end
 
 (R::RandomFunction)()=R(randn(size(R.Σ, 1)))
-
+ 
 
 """
 Covariance function of Gaussian random field.
@@ -157,7 +184,7 @@ Auto-Covariance matrix of a Gaussian random process.
 - `Symmetric{Real}`: auto-covariance matrix
 
 """
-function covariancematrix(P::Vector{<:Point}, process::GaussianRandomField)::Symmetric
+function covariancematrix(P::Vector{<:Point}, process::GaussianRandomField)::Symmetric{<:Real}
     N = length(P)
     A = Matrix{Real}(undef, N, N)
     @threads for i in 1:N
diff --git a/test/runtests.jl b/test/runtests.jl
index 1c287c8..45ff3b5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,8 +3,9 @@ using SpinShuttling
 using UnicodePlots
 using ProgressMeter
 
-include("testintegration.jl")
-include("testfidelity.jl")
-include("teststochastics.jl")
-include("testspectrum.jl")
-include("testthreads.jl")
\ No newline at end of file
+# include("testintegration.jl")
+# include("testfidelity.jl")
+# include("teststochastics.jl")
+# include("testspectrum.jl")
+# include("testthreads.jl")
+include("testdephasingfactor.jl")
\ No newline at end of file
diff --git a/test/testdephasingfactor.jl b/test/testdephasingfactor.jl
new file mode 100644
index 0000000..c250fe2
--- /dev/null
+++ b/test/testdephasingfactor.jl
@@ -0,0 +1,43 @@
+## 
+@testset begin "test dephasing matrix"
+    T=400; L=10; σ = sqrt(2) / 20; M = 5000; N=601; κₜ=1/20;κₓ=1/0.1;
+
+    B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
+    model=OneSpinModel(T,L,N,B)
+    # test customize println
+    println(model)
+
+    f=averagefidelity(model)
+    w=dephasingmatrix(model)
+
+    rho=model.Ψ*model.Ψ'
+    @test rho[1,1]+rho[2,2]≈ 1
+    @test w==w'
+
+    println(w)
+    f_c=(model.Ψ'*(w.*rho)*model.Ψ)
+
+    @test f≈ f_c
+end
+
+## 
+@testset begin "test two spin dephasing matrix"
+    L=10; σ =sqrt(2)/20; M=5000; N=501; T1=200; T0=25*0.05; κₜ=1/20; κₓ=1/0.1;
+    B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
+    model=TwoSpinModel(T0, T1, L, N, B)
+    println(model)
+    f=averagefidelity(model)
+    w=dephasingmatrix(model)
+
+    rho=model.Ψ*model.Ψ'
+    @test sum([rho[i,i] for i in 1:4 ]) ≈ 1
+    @test w==w'
+
+    println(w)
+    f_c=(model.Ψ'*(w.*rho)*model.Ψ)
+
+    println(f)
+    println(f_c)
+    @test f ≈ f_c
+end
+