fix analytics for FB 1/f shuttling 🔥, use adaptive integral for multi…

…ple-dimensional integral

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.2"
 manifest_format = "2.0"
-project_hash = "2e69f859a57c2a1283c3102869192c6f0c03e3c9"
+project_hash = "343c9e0872ed7e6e255309cc978ec15f5c1cf0bc"
 
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
@@ -42,6 +42,11 @@ git-tree-sha1 = "362a287c3aa50601b0bc359053d5c2468f0e7ce0"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
 version = "0.12.11"
 
+[[deps.Combinatorics]]
+git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
+uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+version = "1.0.2"
+
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
 git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248"
@@ -102,6 +107,12 @@ git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 version = "0.8.5"
 
+[[deps.HCubature]]
+deps = ["Combinatorics", "DataStructures", "LinearAlgebra", "QuadGK", "StaticArrays"]
+git-tree-sha1 = "10f37537bbd83e52c63abf6393f209dbd641fedc"
+uuid = "19dc6840-f33b-545b-b366-655c7e3ffd49"
+version = "1.6.0"
+
 [[deps.InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

Project.toml

@@ -1,11 +1,12 @@
-name = "SpinShuttling"
-uuid = "eff49e45-ae6d-44ae-9170-2124a87971c9"
-authors = ["Yuning Zhang"]
-version = "0.2.8"
-
-[deps]
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
+name = "SpinShuttling"
+uuid = "eff49e45-ae6d-44ae-9170-2124a87971c9"
+authors = ["Yuning Zhang"]
+version = "0.2.8"
+
+[deps]
+HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"

src/SpinShuttling.jl

@@ -3,7 +3,7 @@ module SpinShuttling
 using LinearAlgebra
 using Statistics
 using SpecialFunctions
-using QuadGK
+using QuadGK, HCubature
 using UnicodePlots: lineplot, lineplot!
 using Base.Threads
 
@@ -222,10 +222,22 @@ Calculate the dephasingfactor according to a special combinator of the noise seq
 - `c::Vector{Int}`: The combinator of the noise sequence, which should have the same length as the number of spins.
 
 """
+# function dephasingfactor(model::ShuttlingModel, c::Vector{Int})::Real
+#     # @assert length(c) == model.n
+#     R = CompositeRandomFunction(model.R, c)
+#     return characteristicvalue(R)
+# end
+
 function dephasingfactor(model::ShuttlingModel, c::Vector{Int})::Real
-    # @assert length(c) == model.n
-    R = CompositeRandomFunction(model.R, c)
-    return characteristicvalue(R)
+    # rewrite the function to use the hcubature methods
+    M = zeros(model.n, model.n)
+    @threads for i in 1:model.n
+        for j in 1:model.n
+            K(t::AbstractVector)= covariance((t[1], model.X[i](t[1])...), (t[2], model.X[j](t[2])...), model.B)
+            M[i, j] = hcubature(K, [0.0, 0.0], [model.T, model.T], rtol=1e-5)[1]
+        end
+    end
+    return exp.(-sum((c * c').*M)/2)
 end
 
 function dephasingcoeffs(n::Int)::Array{Real,3}
@@ -355,13 +367,14 @@ function W(T::Real, L::Real, B::OrnsteinUhlenbeckField; path=:straight)::Real
 end
 
 function W(T::Real, L::Real, B::PinkLorentzianField; path=:straight)::Real
-    β = T .* B.γ
     γ = L * B.κ
     if path == :straight
+        β = T .* B.γ
         return exp(-B.σ^2 * T^2 * F3(β, γ))
     elseif path == :forthback
-        β /= 2
-        return exp(-B.σ^2 * T^2 * (2*F3(β, γ)+F4(β, γ)))
+        T/=2
+        β = T .* B.γ 
+        return exp(-B.σ^2 * T^2*(2*F3(β, γ)+F4(β, γ)))
     end
 end
 

src/analytics.jl

@@ -75,8 +75,15 @@ end
 Ancillary function for the dephasing of the Pink-Lorentzian noise.
 """
 function F4(β::Tuple{Real,Real}, γ::Real)::Real
-    F(β::Real) = -((1 + exp(-2β)*(1 - 2β) + 2*exp(-β)*(-1 + β))/(4*β^2))+expinti(-2β)-expinti(-β)/2 
-    F(β::Real, γ::Real)::Real = 1 / γ^2 * ((exp(-2β)-1)*γ/β + log((β+γ)/β) +(2α-1)*expinti(-2β)+2*exp(-γ)*expinti(-β)-expinti(-β-γ)-exp(-2γ)*expinti(γ-β)+exp(-2γ)*expinti(2*(γ-β)))
+    F(β::Real) =  -(
+        exp(-2 * β) * (exp(β) - 1) * (2 * β + exp(β) - 1)
+    ) / (2 * β^2) + 
+    2 * expinti(-2 * β) - 
+    expinti(-β)
+
+    F(β::Real, γ::Real)::Real = (
+        (exp(-2 * β) * γ / β) -
+        (γ / β) + log(β + γ) - log(β) +(2 * γ - 1) * expinti(-2 * β) +2 * exp(-γ) * expinti(-β) -expinti(-β - γ) -exp(-2 * γ) * expinti(γ - β) +exp(-2 * γ) * expinti(2 * γ - 2 * β)) / γ^2
     if γ == 0 # pure 1/f noise
         return (F(β[2]) - F(β[1])) / log(β[2] / β[1])
     else

src/integration.jl

@@ -7,10 +7,10 @@ increment `h`
 # Arguments
 - `y::Vector{<:Real}`: f(x).
 - `h::Real`: the step of integral.
-- `method::Symbol=:simpson`: the method of integration 
+- `method::Symbol=:trapezoid`: the method of integration 
 
 """
-function integrate(y::ArrayOrSubArray{<:Real,1}, h::Real; method::Symbol=:simpson)::Real
+function integrate(y::ArrayOrSubArray{<:Real,1}, h::Real; method::Symbol=:trapezoid)::Real
     n = length(y) - 1
     if n%2==1
         println("warning: `y` length (number of intervals) must be odd, switching to first order integration.")
@@ -37,7 +37,7 @@ end
 - `a::Real`: The lower bound of the integration range.
 - `b::Real`: The upper bound of the integration range.
 - `n::Int`: The number of points at which to evaluate `f(x)`. For Simpson's rule, `n` should be an odd number.
-- `method::Symbol=:simpson`: The numerical integration method to use. Options are `:simpson` for Simpson's rule and `:trapezoid` for the Trapezoidal rule.
+- `method::Symbol=:trapezoid`: The numerical integration method to use. Options are `:simpson` for Simpson's rule and `:trapezoid` for the Trapezoidal rule.
 
 # Returns
 - `Real`: The estimated value of the integral of `f(x)` over [a, b].
@@ -48,13 +48,13 @@ result = integrate(sin, 0, π, 101, method=:simpson)
 ```
 """
 function integrate(f::Function, a::Real, b::Real, n::Int;
-    method::Symbol=:simpson)::Real
+    method::Symbol=:trapezoid)::Real
     h = (b - a) / (n - 1)
     y = f.(range(a, b, n))
     return integrate(y, h, method=method)
 end
 
-function integrate(f::Function, x_range::Tuple{Real,Real,Int}; method::Symbol=:simpson)::Real
+function integrate(f::Function, x_range::Tuple{Real,Real,Int}; method::Symbol=:trapezoid)::Real
     return integrate(f, x_range..., method=method)
 end
 
@@ -65,7 +65,7 @@ end
 - `f::Function`: The function to be integrated. It should accept two real numbers (x, y) and return a real number.
 - `x_range::Tuple{Real, Real, Int}`: A tuple representing the range and number of points for the x-dimension: (x_min, x_max, num_points_x).
 - `y_range::Tuple{Real, Real, Int}`: A tuple representing the range and number of points for the y-dimension: (y_min, y_max, num_points_y).
-- `method::Symbol=:simpson`: The numerical integration method to use. Options are `:simpson` for Simpson's rule and `:trapezoid` for the Trapezoidal rule.
+- `method::Symbol=:trapezoid`: The numerical integration method to use. Options are `:simpson` for Simpson's rule and `:trapezoid` for the Trapezoidal rule.
 
 # Returns
 - `Real`: The estimated value of the integral of `f(x, y)` over the defined 2D area.
@@ -76,7 +76,7 @@ result = integrate((x, y) -> x * y, (0, 1, 50), (0, 1, 50), method=:simpson)
 ```
 """
 function integrate(f::Function, x_range::Tuple{Real,Real,Int}, y_range::Tuple{Real,Real,Int};
-    method::Symbol=:simpson)::Real
+    method::Symbol=:trapezoid)::Real
     g = y -> integrate(x -> f(x, y), x_range, method=method)
     return integrate(g, y_range, method=method)
 end
@@ -89,18 +89,18 @@ increments `h_x` and `h_y`.
 - `z::Matrix{<:Real}`: f(x, y).
 - `h_x::Real`: the step size of the integral in the x direction.
 - `h_y::Real`: the step size of the integral in the y direction.
-- `method::Symbol=:simpson`: the method of integration.
+- `method::Symbol=:trapezoid`: the method of integration.
 ...
 """
-function integrate(z::ArrayOrSubArray{<:Real,2}, h_x::Real, h_y::Real; method::Symbol=:simpson)::Real
+function integrate(z::ArrayOrSubArray{<:Real,2}, h_x::Real, h_y::Real; method::Symbol=:trapezoid)::Real
     nrows, ncols = size(z)
     # Integrate along x direction for each y
     integral_x_direction = [integrate((@view z[:, j]), h_x, method=method) for j = 1:ncols]
     # Integrate the result in y direction
     return integrate(integral_x_direction, h_y, method=method)
 end
 
-integrate(z::ArrayOrSubArray{<:Real,2}, h; method::Symbol=:simpson) = integrate(z, h, h; method=method)
+integrate(z::ArrayOrSubArray{<:Real,2}, h; method::Symbol=:trapezoid) = integrate(z, h, h; method=method)
 
 """
 Special methods for the double integral on symmetric matrix with singularity on diagonal entries.

test/testfidelity.jl

@@ -1,57 +1,65 @@
-##
+# ##
 visualize=true
 
-##
+#
 @testset begin "test single spin shuttling fidelity"
-    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)
-
-    f1=statefidelity(model)
-    f2, f2_err=sampling(model, statefidelity, M)
-    f3=1/2*(1+W(T,L,B))
-    @test isapprox(f1, f3,rtol=3e-2)
-    @test isapprox(f2, f3, rtol=3e-2) 
-    println("NI:", f1)
-    println("MC:", f2)
-    println("TH:", f3)
+    T=100; L=10; σ = sqrt(2) / 20; M = 5000; N=601; κₜ=1/20;κₓ=1/0.1;
+    γ = (1e-9, 1e3) # MHz
+    B1 = OrnsteinUhlenbeckField(0, [κₜ, κₓ], σ)
+    B2 = PinkLorentzianField(0, κₓ, σ, γ)
+    for B in (B1,B2)
+        model=OneSpinModel(T,L,N,B)
+        # test customize println
+        println(model)
+
+        f1=statefidelity(model)
+        f2, f2_err=sampling(model, statefidelity, M)
+        f3=1/2*(1+W(T,L,B))
+        @test isapprox(f1, f3,rtol=3e-2)
+        @test isapprox(f2, f3, rtol=3e-2) 
+        println("NI:", f1)
+        println("MC:", f2)
+        println("TH:", f3)
+    end
 end
 
 #
 @testset begin "test single spin forth-back shuttling fidelity"
-    T=200; L=10; σ = sqrt(2) / 20; M = 5000; N=501; κₜ=1/20;κₓ=10; 
+    T=20; σ = sqrt(2) / 20; M = 5000; N=801; κₜ=1/20;κₓ=10;
+    L=0.1; #α=1 
+    γ = (1e-9, 1e3) # MHz 
     # exponential should be smaller than 100
-    B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
+    B1 = OrnsteinUhlenbeckField(0, [κₜ, κₓ], σ)
+    B2 = PinkLorentzianField(0, κₓ, σ, γ)
 
-    if visualize
-        t=range(1e-2*T,T, 10)
-        f_mc=[sampling(OneSpinForthBackModel(T,L,N,B), statefidelity, M)[1] for T in t]
-        f_ni=[statefidelity(OneSpinForthBackModel(T,L,N,B)) for T in t]
-        f_th=[(1+W(T,L,B,path=:forthback))/2 for T in t]
-        fig=lineplot(t, f_mc,
-            xlabel="t", ylabel="F", name="monte-carlo sampling",
-            # ribbon=@. sqrt(f_mc_err/M)
-            )
-            lineplot!(fig, t, f_ni, name="numerical integration")
-            lineplot!(fig, t, f_th, name="theoretical fidelity")
-        display(fig)
-    end
+    for B in (B1,B2)
+        if visualize
+            t=range(1e-2*T,T, 10)
+            f_mc=[sampling(OneSpinForthBackModel(T,L,N,B), statefidelity, M)[1] for T in t]
+            f_ni=[statefidelity(OneSpinForthBackModel(T,L,N,B)) for T in t]
+            f_th=[(1+W(T,L,B,path=:forthback))/2 for T in t]
+            fig=lineplot(t, f_mc,
+                xlabel="t", ylabel="F", name="monte-carlo sampling",
+                # ribbon=@. sqrt(f_mc_err/M)
+                )
+                lineplot!(fig, t, f_ni, name="numerical integration")
+                lineplot!(fig, t, f_th, name="theoretical fidelity")
+            display(fig)
+        end
 
-    model=OneSpinForthBackModel(T,L,N,B)
-    # test customize println
-    println(model)
+        model=OneSpinForthBackModel(T,L,N,B)
+        # test customize println
+        println(model)
 
-    f1=statefidelity(model)
-    f2, f2_err=sampling(model, statefidelity, M)
-    f3=1/2*(1+W(T, L, B, path=:forthback))
-    @test isapprox(f1, f3,rtol=3e-2)
-    @test isapprox(f2, f3, rtol=3e-2) 
-    println("NI:", f1)
-    println("MC:", f2)
-    println("TH:", f3)
+        f1=statefidelity(model)
+        f2, f2_err=sampling(model, statefidelity, M)
+        f3=1/2*(1+W(T, L, B, path=:forthback))
+        @test isapprox(f1, f3,rtol=3e-2)
+        @test isapprox(f2, f3, rtol=3e-2) 
+        println("NI:", f1)
+        println("MC:", f2)
+        println("TH:", f3)
+    end
 end
 
 
@@ -114,42 +122,20 @@ end
 end
 
 ##
-@testset begin "test two spin parallel shuttling fidelity"
-    L=10; σ =sqrt(2)/20; M=5000; N=501; T=200; κₜ=1/20; κₓ=1/0.1;
-    D=0.3;
-    B=OrnsteinUhlenbeckField(0,[κₜ,κₓ,κₓ],σ)
-    model=TwoSpinParallelModel(T, D, L, N, B)
-    if visualize
-        display(heatmap(collect(model.R.Σ), title="cross covariance matrix, two spin EPR"))
-    end
-    f1=statefidelity(model)
-    f2, f2_err=sampling(model, statefidelity, M)
-    w=exp(-σ^2 / (8 *κₜ*κₓ*κₓ) / κₜ^2 *(1-exp(-κₓ*D)) * SpinShuttling.P1(κₜ*T, κₓ*L))
-    f3=1/2*(1+w)
-    @test isapprox(f1, f3,rtol=3e-2)
-    @test isapprox(f2, f3, rtol=3e-2) 
-    println("NI:", f1)
-    println("MC:", f2)
-    println("TH:", f3)
-end
-
-# @testset "test two spin parallel shuttling fidelity with 1/f noise" begin
-#     σ = sqrt(2)/20; M = 5000; N=501; L=10; γ=(1e-9,1e3); # MHz
-#     # 0.01 ~ 100 μs
-#     # v = 0.1 ~ 1000 m/s
-#     v=1; T=L/v; κₓ=10;
-#     # T=10 
-#     # 1/T=0.1 N/T=20 
+# @testset begin "test two spin parallel shuttling fidelity"
+#     L=10; σ =sqrt(2)/20; M=5000; N=501; T=200; κₜ=1/20; κₓ=1/0.1;
 #     D=0.3;
-#     B=PinkLorentzianField(0,[κₓ,κₓ],σ, γ)
+#     B=OrnsteinUhlenbeckField(0,[κₜ,κₓ,κₓ],σ)
 #     model=TwoSpinParallelModel(T, D, L, N, B)
-
+#     if visualize
+#         display(heatmap(collect(model.R.Σ), title="cross covariance matrix, two spin EPR"))
+#     end
 #     f1=statefidelity(model)
 #     f2, f2_err=sampling(model, statefidelity, M)
-#     w=exp(-B.σ^2 * T^2 * 2(1-exp(-κₓ*D)) * SpinShuttling.F3(γ.*T, κₓ*L))
+#     w=exp(-σ^2 / (8 *κₜ*κₓ*κₓ) / κₜ^2 *(1-exp(-κₓ*D)) * SpinShuttling.P1(κₜ*T, κₓ*L))
 #     f3=1/2*(1+w)
 #     @test isapprox(f1, f3,rtol=3e-2)
-#     @test isapprox(f2, f3, rtol=3e-2)
+#     @test isapprox(f2, f3, rtol=3e-2) 
 #     println("NI:", f1)
 #     println("MC:", f2)
 #     println("TH:", f3)