From 2c7cd3ef0fb7118b01883eb85246605542c96ef7 Mon Sep 17 00:00:00 2001 From: Yuning Date: Tue, 27 Aug 2024 22:14:37 +0200 Subject: [PATCH] fix analytics for FB 1/f shuttling :fire:, use adaptive integral for multiple-dimensional integral --- Manifest.toml | 13 ++++- Project.toml | 23 ++++---- src/SpinShuttling.jl | 27 ++++++--- src/analytics.jl | 11 +++- src/integration.jl | 20 +++---- test/testfidelity.jl | 132 +++++++++++++++++++------------------------ 6 files changed, 122 insertions(+), 104 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index b982324..1f6ba69 100644 --- a/Manifest.toml +++ b/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" diff --git a/Project.toml b/Project.toml index 0855b36..27929d6 100644 --- a/Project.toml +++ b/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" diff --git a/src/SpinShuttling.jl b/src/SpinShuttling.jl index 51513c1..19a1f2b 100644 --- a/src/SpinShuttling.jl +++ b/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 diff --git a/src/analytics.jl b/src/analytics.jl index af46a5b..9cd313f 100644 --- a/src/analytics.jl +++ b/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 diff --git a/src/integration.jl b/src/integration.jl index 9cf15fd..69f339f 100644 --- a/src/integration.jl +++ b/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,10 +89,10 @@ 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] @@ -100,7 +100,7 @@ function integrate(z::ArrayOrSubArray{<:Real,2}, h_x::Real, h_y::Real; method::S 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. diff --git a/test/testfidelity.jl b/test/testfidelity.jl index 9ecf591..b13c1f7 100644 --- a/test/testfidelity.jl +++ b/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)