diff --git a/src/SpinShuttling.jl b/src/SpinShuttling.jl index 19a1f2b..cbfc378 100644 --- a/src/SpinShuttling.jl +++ b/src/SpinShuttling.jl @@ -200,14 +200,14 @@ end """ Calculate the dephasing matrix of a given spin shuttling model. """ -function dephasingmatrix(model::ShuttlingModel)::Matrix{<:Real} +function dephasingmatrix(model::ShuttlingModel; method::Symbol=:simpson)::Matrix{<: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] - W[j, k] = dephasingfactor(model, c) + W[j, k] = dephasingfactor(model, c, method=method) W[k, j] = W[j, k] end end @@ -220,24 +220,24 @@ Calculate the dephasingfactor according to a special combinator of the noise seq # Arguments - `model::ShuttlingModel`: The spin shuttling model - `c::Vector{Int}`: The combinator of the noise sequence, which should have the same length as the number of spins. - +- `method::Symbol`: The method to calculate the characteristic value, `:trapezoid`, `:simpson`, `:adaptive` """ -# 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 +function dephasingfactor(model::ShuttlingModel, c::Vector{Int}; method::Symbol=:simpson)::Real # 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] + @assert method in [:trapezoid, :simpson, :adaptive] + if method==:adaptive + 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) + else + R = CompositeRandomFunction(model.R, c) + return characteristicvalue(R, method=method) end - return exp.(-sum((c * c').*M)/2) end function dephasingcoeffs(n::Int)::Array{Real,3} @@ -259,9 +259,9 @@ of the covariance matrix. # Arguments - `model::ShuttlingModel`: The spin shuttling model """ -function statefidelity(model::ShuttlingModel)::Real +function statefidelity(model::ShuttlingModel; method::Symbol=:simpson)::Real Ψ= model.Ψ - w=dephasingmatrix(model) + w=dephasingmatrix(model, method=method) ρt=w.*(Ψ*Ψ') return Ψ'*ρt*Ψ end diff --git a/src/stochastics.jl b/src/stochastics.jl index 2046abe..72a5266 100644 --- a/src/stochastics.jl +++ b/src/stochastics.jl @@ -202,12 +202,12 @@ Compute the characteristic functional of the process from the numerical quadrature of the covariance matrix. Using Simpson's rule by default. """ -function characteristicfunction(R::RandomFunction)::Tuple{Vector{<:Real},Vector{<:Number}} +function characteristicfunction(R::RandomFunction; method::Symbol=:simpson)::Tuple{Vector{<:Real},Vector{<:Number}} # need further optimization dt = R.P[2][1] - R.P[1][1] N = size(R.Σ, 1) @assert N % 2 == 1 - χ(j::Int) = exp.(1im * integrate(view(R.μ, 1:j), dt)) * exp.(-integrate(view(R.Σ, 1:j, 1:j), dt, dt) / 2) + χ(j::Int) = exp.(1im * integrate(view(R.μ, 1:j), dt, method=method)) * exp.(-integrate(view(R.Σ, 1:j, 1:j), dt, dt,method=method) / 2) t = [p[1] for p in R.P[2:2:N-1]] f = [χ(j) for j in 3:2:N] # only for simpson's rule return (t, f) @@ -218,7 +218,9 @@ Compute the final phase of the characteristic functional of the process from the numerical quadrature of the covariance matrix. Using Simpson's rule by default. """ -function characteristicvalue(R::RandomFunction)::Number +function characteristicvalue(R::RandomFunction; method::Symbol=:simpson)::Number dt = R.P[2][1] - R.P[1][1] - return exp.(1im * integrate(R.μ, dt)) * exp.(-integrate((@view R.Σ[:, :]), dt, dt) / 2) + f1=exp.(1im * integrate(R.μ, dt, method=method)) + f2=exp.(-integrate((@view R.Σ[:, :]), dt, dt, method=method) / 2) + return f1*f2 end diff --git a/test/testfidelity.jl b/test/testfidelity.jl index b13c1f7..68bcd2a 100644 --- a/test/testfidelity.jl +++ b/test/testfidelity.jl @@ -111,7 +111,7 @@ end if visualize display(heatmap(collect(model.R.Σ), title="cross covariance matrix, two spin EPR")) end - f1=statefidelity(model) + f1=statefidelity(model, method=:adaptive) f2, f2_err=sampling(model, statefidelity, M) f3=1/2*(1+W(T0, T1, L,B)) @test isapprox(f1, f3,rtol=3e-2) diff --git a/test/testintegration.jl b/test/testintegration.jl index b6dd407..05b31b3 100644 --- a/test/testintegration.jl +++ b/test/testintegration.jl @@ -36,7 +36,7 @@ using QuadGK: quadgk h_y=(y_range[2]-y_range[1])/(y_range[3]-1); @test isapprox(integrate(g, (-0.0,2.0,21),(-2,2.0,11)), 64/3; rtol=1e-8); @test z isa Matrix{<:Real} && h_x isa Real && h_y isa Real - @test isapprox(integrate(z, h_x, h_y), 64/3; rtol=1e-8); + @test isapprox(integrate(z, h_x, h_y), 64/3; rtol=1e-5); end ## @testset "cone" begin