diff --git a/examples/Bayesian_bounds_multiparameter.jl b/examples/Bayesian_bounds_multiparameter.jl new file mode 100644 index 0000000..44930e1 --- /dev/null +++ b/examples/Bayesian_bounds_multiparameter.jl @@ -0,0 +1,72 @@ +using Trapz +include("../src/QuanEstimation.jl") + +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +# free Hamiltonian +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +function H0_func(x) + return 0.5*x[2]*(sx*cos(x[1])+sz*sin(x[1])) +end +function dH_func(x) + return [0.5*x[2]*(-sx*sin(x[1])+sz*cos(x[1])), 0.5*(sx*cos(x[1])+sz*sin(x[1]))] +end +# dynamics +tspan = range(0.0, stop=1.0, length=1000) +# dissipation +decay_opt = [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])] +gamma = [0.0] +# prior distribution +function p_func(x, y, mu_x, mu_y, sigmax, sigmay, r) + term1 = ((x-mu_x)/sigmax)^2-2*r*(((x-mu_x)/sigmax))*(((y-mu_y)/sigmay))+(((y-mu_y)/sigmay))^2 + term2 = exp(-term1/2.0/(1-r^2)) + term3 = 2*pi*sigmax*sigmay*sqrt(1-r^2) + return term2/term3 +end +function dp_func(x, y, mu_x, mu_y, sigmax, sigmay, r) + term1 = -(2*((x-mu_x)/sigmax^2)-2*r*((y-mu_y)/sigmay)/sigmax)/2.0/(1-r^2) + term2 = -(2*((y-mu_y)/sigmay^2)-2*r*((x-mu_x)/sigmax)/sigmay)/2.0/(1-r^2) + p = p_func(x, y, mu_x, mu_y, sigmax, sigmay, r) + return [term1*p, term2*p] +end +x = [range(-pi/2.0, stop=pi/2.0, length=100), range(pi/2-0.1, stop=pi/2+0.1, length=10)].|>Vector +sigmax, sigmay = 0.5, 1.0 +mu_x, mu_y = 0.0, 0.0 +r = 0.5 +para_num = length(x) + +p_tp = Matrix{Float64}(undef, length.(x)...) +dp_tp = Matrix{Vector{Float64}}(undef, length.(x)...) +rho = Matrix{Matrix{ComplexF64}}(undef, length.(x)...) +drho = Matrix{Vector{Matrix{ComplexF64}}}(undef, length.(x)...) +for i = 1:length(x[1]) + for j = 1:length(x[2]) + x_tp = [x[1][i], x[2][j]] + H0_tp = H0_func(x_tp) + dH_tp = dH_func(x_tp) + rho_tp, drho_tp = QuanEstimation.expm(H0_tp, dH_tp, [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])], [zeros(length(tspan)-1)], rho0, tspan, decay_opt, gamma) + rho[i,j] = rho_tp[end] + drho[i,j] = drho_tp[end] + p_tp[i,j] = p_func(x[1][i], x[2][j], mu_x, mu_y, sigmax, sigmay, r) + dp_tp[i,j] = dp_func(x[1][i], x[2][j], mu_x, mu_y, sigmax, sigmay, r) + end +end +c = trapz(tuple(x...), p_tp) +p = p_tp/c +dp = dp_tp/c + +f_BCRB1 = QuanEstimation.BCRB(x, p, rho, drho, M=nothing, btype=1) +f_BCRB2 = QuanEstimation.BCRB(x, p, rho, drho, M=nothing, btype=2) +f_VTB1 = QuanEstimation.VTB(x, p, dp, rho, drho, M=nothing, btype=1) +f_VTB2 = QuanEstimation.VTB(x, p, dp, rho, drho, M=nothing, btype=2) + +f_BQCRB1 = QuanEstimation.BQCRB(x, p, rho, drho, btype=1) +f_BQCRB2 = QuanEstimation.BQCRB(x, p, rho, drho, btype=2) +f_QVTB1 = QuanEstimation.QVTB(x, p, dp, rho, drho, btype=1) +f_QVTB2 = QuanEstimation.QVTB(x, p, dp, rho, drho, btype=2) + +for f in [f_BCRB1,f_BCRB2,f_VTB1,f_VTB2,f_BQCRB1,f_BQCRB2,f_QVTB1,f_QVTB2] + println(f) +end \ No newline at end of file diff --git a/examples/Bayesian_bounds_singleparameter.jl b/examples/Bayesian_bounds_singleparameter.jl new file mode 100644 index 0000000..9fb1a0c --- /dev/null +++ b/examples/Bayesian_bounds_singleparameter.jl @@ -0,0 +1,60 @@ +using Trapz +include("../src/QuanEstimation.jl") + +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +# free Hamiltonian +B = pi/2.0 +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +function H0_func(x) + return 0.5*B*(sx*cos(x)+sz*sin(x)) +end +function dH_func(x) + return [0.5*B*(-sx*sin(x)+sz*cos(x))] +end +# dynamics +decay_opt = [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])] +gamma = [0.0] +tspan = range(0.0, stop=1.0, length=1000) +# prior distribution +function p_func(x, mu, eta) + return exp(-(x-mu)^2/(2*eta^2))/(eta*sqrt(2*pi)) +end +function dp_func(x, mu, eta) + return -(x-mu)*exp(-(x-mu)^2/(2*eta^2))/(eta^3*sqrt(2*pi)) +end +x = [range(-pi/2.0, stop=pi/2.0, length=100)].|>Vector +mu = 0.0 +eta = 0.2 +p_tp = [p_func(x[1][i], mu, eta) for i in 1:length(x[1])] +dp_tp = [dp_func(x[1][i], mu, eta) for i in 1:length(x[1])] +c = trapz(x[1], p_tp) +p = p_tp/c +dp = dp_tp/c + +rho = Vector{Matrix{ComplexF64}}(undef, length(x[1])) +drho = Vector{Vector{Matrix{ComplexF64}}}(undef, length(x[1])) +for i = 1:length(x[1]) + H0_tp = H0_func(x[1][i]) + dH_tp = dH_func(x[1][i]) + rho_tp, drho_tp = QuanEstimation.expm(H0_tp, dH_tp, [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])], [zeros(length(tspan)-1)], rho0, tspan, decay_opt, gamma) + rho[i] = rho_tp[end] + drho[i] = drho_tp[end] +end + +f_BCRB1 = QuanEstimation.BCRB(x, p, rho, drho, M=nothing, btype=1) +f_BCRB2 = QuanEstimation.BCRB(x, p, rho, drho, M=nothing, btype=2) +f_VTB1 = QuanEstimation.VTB(x, p, dp, rho, drho, M=nothing, btype=1) +f_VTB2 = QuanEstimation.VTB(x, p, dp, rho, drho, M=nothing, btype=2) + +f_BQCRB1 = QuanEstimation.BQCRB(x, p, rho, drho, btype=1) +f_BQCRB2 = QuanEstimation.BQCRB(x, p, rho, drho, btype=2) +f_QVTB1 = QuanEstimation.QVTB(x, p, dp, rho, drho, btype=1) +f_QVTB2 = QuanEstimation.QVTB(x, p, dp, rho, drho, btype=2) +f_QZZB = QuanEstimation.QZZB(x, p, rho) + +for f in [f_BCRB1,f_BCRB2,f_VTB1,f_VTB2,f_BQCRB1,f_BQCRB2,f_QVTB1,f_QVTB2,f_QZZB] + println(f) +end \ No newline at end of file diff --git a/examples/Bayesian_estimation_multiparameter.jl b/examples/Bayesian_estimation_multiparameter.jl new file mode 100644 index 0000000..baac975 --- /dev/null +++ b/examples/Bayesian_estimation_multiparameter.jl @@ -0,0 +1,50 @@ +using Random +using StatsBase +include("../src/QuanEstimation.jl") + +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +# free Hamiltonian +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +function H0_func(x) + return 0.5*x[2]*(sx*cos(x[1])+sz*sin(x[1])) +end +function dH_func(x) + return [0.5*x[2]*(-sx*sin(x[1])+sz*cos(x[1])), 0.5*(sx*cos(x[1])+sz*sin(x[1]))] +end +# measurement +M1 = 0.5*[1.0+0.0im 1.0; 1.0 1.0] +M2 = 0.5*[1.0+0.0im -1.0; -1.0 1.0] +M = [M1, M2] +# dynamics +tspan = range(0.0, stop=1.0, length=1000) +# dissipation +decay_opt = [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])] +gamma = [0.0] +# prior distribution +x = [range(0.0, stop=pi/2.0, length=100), range(pi/2-0.1, stop=pi/2+0.1, length=10)].|>Vector +p = (1.0/(x[1][end]-x[1][1]))*(1.0/(x[2][end]-x[2][1]))*ones((length(x[1]), length(x[2]))) + +rho = Matrix{Matrix{ComplexF64}}(undef, length.(x)...) +for i = 1:length(x[1]) + for j = 1:length(x[2]) + x_tp = [x[1][i], x[2][j]] + H0_tp = H0_func(x_tp) + dH_tp = dH_func(x_tp) + rho_tp, drho_tp = QuanEstimation.expm(H0_tp, dH_tp, [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])], [zeros(length(tspan)-1)],rho0, tspan, decay_opt, gamma) + rho[i,j] = rho_tp[end] + end +end + +Random.seed!(1234) +y = [0 for i in 1:500] +res_rand = sample(1:length(y), 125, replace=false) +for i in 1:length(res_rand) + y[res_rand[i]] = 1 +end +pout, xout = QuanEstimation.Bayes(x, p, rho, y; M=M, savefile=false) +println(xout) +Lout, xout = QuanEstimation.MLE(x, rho, y; M=M, savefile=false) +println(xout) diff --git a/examples/Bayesian_estimation_singleparameter.jl b/examples/Bayesian_estimation_singleparameter.jl new file mode 100644 index 0000000..b73b913 --- /dev/null +++ b/examples/Bayesian_estimation_singleparameter.jl @@ -0,0 +1,51 @@ +using Random +using StatsBase +include("../src/QuanEstimation.jl") + +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +# Hamiltonian +B = pi/2.0 +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +function H0_func(x) + return 0.5*B*(sx*cos(x)+sz*sin(x)) +end +function dH_func(x) + return [0.5*B*(-sx*sin(x)+sz*cos(x))] +end +# measurement +M1 = 0.5*[1.0+0.0im 1.0; 1.0 1.0] +M2 = 0.5*[1.0+0.0im -1.0; -1.0 1.0] +M = [M1, M2] +# dynamics +decay_opt = [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])] +gamma = [0.0] +tspan = range(0.0, stop=1.0, length=1000) + +#### prior distribution #### +x = [range(0.0, stop=pi/2.0, length=100)].|>Vector +p = (1.0/(x[1][end]-x[1][1]))*ones(length(x[1])) +dp = zeros(length(x[1])) + +rho = Vector{Matrix{ComplexF64}}(undef, length(x[1])) +for i = 1:length(x[1]) + H0_tp = H0_func(x[1][i]) + dH_tp = dH_func(x[1][i]) + rho_tp, drho_tp = QuanEstimation.expm(H0_tp, dH_tp, [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])], [zeros(length(tspan)-1)], rho0, tspan, decay_opt, gamma) + rho[i] = rho_tp[end] +end + +# generate experimental results +Random.seed!(1234) +y = [0 for i in 1:500] +res_rand = sample(1:length(y), 125, replace=false) +for i in 1:length(res_rand) + y[res_rand[i]] = 1 +end + +pout, xout = QuanEstimation.Bayes(x, p, rho, y; M=M, savefile=false) +println(xout) +Lout, xout = QuanEstimation.MLE(x, rho, y, M=M; savefile=false) +println(xout) diff --git a/examples/Comopt_Kraus_multiparameter.jl b/examples/Comopt_Kraus_multiparameter.jl new file mode 100644 index 0000000..099230b --- /dev/null +++ b/examples/Comopt_Kraus_multiparameter.jl @@ -0,0 +1,49 @@ +using Random +using StableRNGs +include("../src/QuanEstimation.jl") + +# Kraus operators for the generalized qubit amplitude damping +n, p = 0.1, 0.1 +ψ0 = [1., 0] +ψ1 = [0., 1] +K0 =sqrt(1-n)*(ψ0*ψ0' + sqrt(1-p)*ψ1*ψ1') +K1 = sqrt(p-p*n)*ψ0*ψ1' +K2 = sqrt(n)*(sqrt(1-p)*ψ0*ψ0' + ψ1*ψ1') +K3 = sqrt(p*n)*ψ1*ψ0' +K = [K0, K1, K2, K3] + +dK0_n = -0.5*(ψ0*ψ0'+sqrt(1-p)*ψ1*ψ1')/sqrt(1-n) +dK1_n = -0.5*p*ψ0*ψ1'/sqrt(p-p*n) +dK2_n = 0.5*(sqrt(1-p)*ψ0*ψ0'+ψ1*ψ1')/sqrt(n) +dK3_n = 0.5*p*ψ1*ψ0'/sqrt(p*n) +dK0_p = -0.5*sqrt(1-n)*ψ1*ψ1'/sqrt(1-p) +dK1_p = 0.5*(1-n)*ψ0*ψ1'/sqrt(p-p*n) +dK2_p = -0.5*sqrt(n)*ψ0*ψ0'/sqrt(1-p) +dK3_p = -0.5*sqrt(n)*ψ0*ψ0'/sqrt(1-p) +dK3_p = 0.5*n*ψ1*ψ0'/sqrt(p*n) +dK = [[dK0_n, dK0_p], [dK1_n, dK1_p], [dK2_n, dK2_p], [dK3_n, dK3_p]] + +# measurement +dim = size(K[1],1) +POVM_basis = QuanEstimation.SIC(dim) +M_num = dim + +# state measurement optimization +opt = QuanEstimation.SMopt() + +# initial state +ρ₀ = 0.5*ones(2,2) + +# dynamics +dynamics = QuanEstimation.Kraus(opt, K, dK) + +# set objective +obj = QuanEstimation.CFIM_Obj() + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) diff --git a/examples/Comopt_Kraus_single_parameter.jl b/examples/Comopt_Kraus_single_parameter.jl new file mode 100644 index 0000000..960018e --- /dev/null +++ b/examples/Comopt_Kraus_single_parameter.jl @@ -0,0 +1,27 @@ +using Random +using StableRNGs +include("../src/QuanEstimation.jl") + +# Kraus operators for the amplitude damping channel +γ = 0.1 +K1 = [1. 0; 0 sqrt(1-γ)] +K2 = [0. sqrt(γ); 0 0] +K = [K1, K2] + +dK1 = [1. 0; 0 -sqrt(1-γ)/2] +dK2 = [0. sqrt(γ); 0 0] +dK = [[dK1], [dK2]] + +opt = QuanEstimation.SMopt() +dynamics = QuanEstimation.Kraus(opt, K, dK) + +# set objective +obj = QuanEstimation.CFIM_Obj() + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) diff --git a/examples/Comopt_multiparameter.jl b/examples/Comopt_multiparameter.jl new file mode 100644 index 0000000..2a6929a --- /dev/null +++ b/examples/Comopt_multiparameter.jl @@ -0,0 +1,89 @@ +using Random +using LinearAlgebra +include("../src/QuanEstimation.jl") + +# initial state +ρ₀ = zeros(ComplexF64, 6, 6) +ρ₀[1:4:5, 1:4:5].=0.5 +dim = size(ρ₀, 1) + +# Hamiltonian +sx = [0. 1; 1 0] +sy = [0. -im; im 0] +sz = [1. 0; 0 -1] +s1 = [0. 1 0; 1 0 1; 0 1 0] / sqrt(2) +s2 = [0. -im 0; im 0 -im; 0 im 0] / sqrt(2) +s3 = [1. 0 0; 0 0 0; 0 0 -1] +Is = I1, I2, I3 =[ kron(I(3), sx), kron(I(3), sy), kron(I(3), sz)] +S = S1, S2, S3 = [kron(s1, I(2)), kron(s2, I(2)), kron(s3, I(2))] +B = B1, B2, B3 = [5.0e-4, 5.0e-4, 5.0e-4] +cons = 100 +D = (2pi * 2.87 * 1000) / cons +gS = (2pi * 28.03 * 1000) / cons +gI = (2pi * 4.32) / cons +A1 = (2pi * 3.65) / cons +A2 = (2pi * 3.03) / cons +H0 = sum([ + D*kron(s3^2, I(2)), + sum(gS * B .* S), + sum(gI * B .* Is), + A1 * (kron(s1, sx) + kron(s2, sy)), + A2 * kron(s3, sz) +]) +dH = gS * S + gI * Is +Hc = [S1, S2, S3] + +# dissipation +decay = [[S3, 2pi/cons]] + +# measurement +M = [QuanEstimation.basis(dim, i)*QuanEstimation.basis(dim, i)' for i in 1:dim] + +#dynamics +tspan = range(0, 2, length=4000) + +# initial control coefficients +cnum = 10 +ctrl0 = [zeros(cnum) for _ in 1:length(Hc)] + +# state control optimization +opt = QuanEstimation.SCopt() +dynamics = QuanEstimation.Lindblad(opt, tspan, H0, dH, Hc, decay) + +# control measurement optimization +opt = QuanEstimation.CMopt() +dynamics = QuanEstimation.Lindblad(opt, tspan, ρ₀, H0, dH, Hc, decay) + +# state measurement optimization +opt = QuanEstimation.SMopt() +dynamics = QuanEstimation.Lindblad(opt, tspan, H0, dH, Hc, ctrl0, decay) + +# # state control measurement optimization +# opt = QuanEstimation.SCMopt() +# dynamics = QuanEstimation.Lindblad(opt, tspan, H0, dH, Hc, decay) + +# control algorithm: AD +if opt isa QuanEstimation.SCopt + alg = QuanEstimation.AD(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) + obj = QuanEstimation.QFIM_Obj() + obj = QuanEstimation.CFIM_Obj(M=M) + # obj = QuanEstimation.HCRB_Obj() + + QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) +end + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) \ No newline at end of file diff --git a/examples/Comopt_single_parameter.jl b/examples/Comopt_single_parameter.jl new file mode 100644 index 0000000..fb90a30 --- /dev/null +++ b/examples/Comopt_single_parameter.jl @@ -0,0 +1,69 @@ +using Random +using StableRNGs +include("../src/QuanEstimation.jl") + +# initial state +ρ₀ = 0.5 * ones(2,2) +# Hamiltonian +ω₀ = 1.0 +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +H0 = 0.5 * ω₀ * sz +dH = [0.5 * sz] +Hc = [sx, sy, sz] +# measurement +M1 = 0.5*[1.0+0.0im 1.0; 1.0 1.0] +M2 = 0.5*[1.0+0.0im -1.0; -1.0 1.0] +M = [M1, M2] +# dissipation +sp = [0 1; 0 0.0im] +sm = [0 0; 1 0.0im] +decay = [[sp, 0.0], [sm, 0.1]] +# dynamics +tspan = range(0.0, 10.0, length=2500) +# initial control coefficients +cnum = length(tspan) - 1 +ctrl0 = [zeros(cnum) for _ in 1:length(Hc)] + +# state control optimization +opt = QuanEstimation.SCopt() +dynamics = QuanEstimation.Lindblad(opt, tspan, H0, dH, Hc, decay) + +# control measurement optimization +opt = QuanEstimation.CMopt() +dynamics = QuanEstimation.Lindblad(opt, tspan, ρ₀, H0, dH, Hc, decay) + +# state measurement optimization +opt = QuanEstimation.SMopt() +dynamics = QuanEstimation.Lindblad(opt, tspan, H0, dH, Hc, ctrl0, decay) + +# state control measurement optimization +opt = QuanEstimation.SCMopt() +dynamics = QuanEstimation.Lindblad(opt, tspan, H0, dH, Hc, decay) + +# control algorithm: AD +if opt isa QuanEstimation.SCopt + alg = QuanEstimation.AD(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) + obj = QuanEstimation.QFIM_Obj() + obj = QuanEstimation.CFIM_Obj(M=M) + # obj = QuanEstimation.HCRB_Obj() + + QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) +end + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) \ No newline at end of file diff --git a/examples/Copt_multi_parameter.jl b/examples/Copt_multi_parameter.jl new file mode 100644 index 0000000..115428d --- /dev/null +++ b/examples/Copt_multi_parameter.jl @@ -0,0 +1,102 @@ +using Random +using LinearAlgebra +include("../src/QuanEstimation.jl") + +# initial state +ρ₀ = zeros(ComplexF64, 6, 6) +ρ₀[1:4:5, 1:4:5].=0.5 +dim = size(ρ₀, 1) + +# Hamiltonian +sx = [0. 1; 1 0] +sy = [0. -im; im 0] +sz = [1. 0; 0 -1] +s1 = [0. 1 0; 1 0 1; 0 1 0] / sqrt(2) +s2 = [0. -im 0; im 0 -im; 0 im 0] / sqrt(2) +s3 = [1. 0 0; 0 0 0; 0 0 -1] +Is = I1, I2, I3 =[ kron(I(3), sx), kron(I(3), sy), kron(I(3), sz)] +S = S1, S2, S3 = [kron(s1, I(2)), kron(s2, I(2)), kron(s3, I(2))] +B = B1, B2, B3 = [5.0e-4, 5.0e-4, 5.0e-4] +cons = 100 +D = (2pi * 2.87 * 1000) / cons +gS = (2pi * 28.03 * 1000) / cons +gI = (2pi * 4.32) / cons +A1 = (2pi * 3.65) / cons +A2 = (2pi * 3.03) / cons +H0 = sum([ + D*kron(s3^2, I(2)), + sum(gS * B .* S), + sum(gI * B .* Is), + A1 * (kron(s1, sx) + kron(s2, sy)), + A2 * kron(s3, sz) +]) +dH = gS * S + gI * Is +Hc = [S1, S2, S3] + +# dissipation +decay = [[S3, 2pi/cons]] + +# measurement +M = [QuanEstimation.basis(dim, i)*QuanEstimation.basis(dim, i)' for i in 1:dim] + +#dynamics +tspan = range(0, 2, length=4000) + +# initial contrl coefficients +rng = MersenneTwister(1234) +cnum = 10 +Hc_num = length(Hc) +ini_1 = [zeros(cnum) for _ in 1:Hc_num] +ini_2 = 0.2 .* [ones(cnum) for _ in 1:Hc_num] +ini_3 = -0.2 .* [ones(cnum) for _ in 1:Hc_num] +ini_4 = [[range(-0.2, 0.2, length=cnum)...] for _ in 1:Hc_num] +ini_5 = [[range(-0.2, 0., length=cnum)...] for _ in 1:Hc_num] +ini_6 = [[range(0., 0.2, length=cnum)...] for _ in 1:Hc_num] +ini_7 = [-0.2*ones(cnum) + 0.01*rand(rng, cnum) for _ in 1:Hc_num] +ini_8 = [-0.2*ones(cnum) + 0.01*rand(rng, cnum) for _ in 1:Hc_num] +ini_9 = [-0.2*ones(cnum) + 0.05*rand(rng, cnum) for _ in 1:Hc_num] +ini_10 = [-0.2*ones(cnum) + 0.05*rand(rng, cnum) for _ in 1:Hc_num] +ctrl0 = [Symbol("ini_", i)|>eval for i in 1:10] + +opt = QuanEstimation.Copt(ctrl=ini_1, ctrl_bound=[-0.2, 0.2]) +dynamics = QuanEstimation.Lindblad(opt, tspan, ρ₀, H0, dH, decay, Hc) + +# # control algorithm: GRAPE +# alg = QuanEstimation.GRAPE(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +# obj = QuanEstimation.QFIM_Obj() +# obj = QuanEstimation.CFIM_Obj(M=M) +# # obj = QuanEstimation.HCRB_Obj() + +# QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# # control algorithm: auto-GRAPE +# alg = QuanEstimation.autoGRAPE(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +# obj = QuanEstimation.QFIM_Obj() +# obj = QuanEstimation.CFIM_Obj(M=M) +# # obj = QuanEstimation.HCRB_Obj() + +# QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, ini_particle=(ctrl0, ), c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, ini_population=(ctrl0, ), c=1.0, cr=0.5, rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DDPG +alg = QuanEstimation.DDPG(max_episode=100, layer_num=4, layer_dim=250, rng=StableRNG(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) \ No newline at end of file diff --git a/examples/Copt_single_parameter.jl b/examples/Copt_single_parameter.jl new file mode 100644 index 0000000..92738de --- /dev/null +++ b/examples/Copt_single_parameter.jl @@ -0,0 +1,62 @@ +using Random +using StableRNGs +include("../src/QuanEstimation.jl") + +# initial state +ρ₀ = 0.5 * ones(2,2) +# Hamiltonian +ω₀ = 1.0 +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +H0 = 0.5 * ω₀ * sz +dH = [0.5 * sz] +Hc = [sx, sy, sz] +# measurement +M1 = 0.5*[1.0+0.0im 1.0; 1.0 1.0] +M2 = 0.5*[1.0+0.0im -1.0; -1.0 1.0] +M = [M1, M2] +# dissipation +sp = [0 1; 0 0.0im] +sm = [0 0; 1 0.0im] +decay = [[sp, 0.0], [sm, 0.1]] +# dynamics +tspan = range(0.0, 10.0, length=2500) +# initial control coefficients +cnum = length(tspan) - 1 +ctrl0 = [zeros(cnum) for _ in 1:length(Hc)] + +opt = QuanEstimation.Copt() +dynamics = QuanEstimation.Lindblad(opt, tspan, ρ₀, H0, dH, decay, Hc) + +# control algorithm: GRAPE +alg = QuanEstimation.GRAPE(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +# control algorithm: auto-GRAPE +alg = QuanEstimation.autoGRAPE(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +# control algorithm: DDPG +alg = QuanEstimation.DDPG(max_episode=100, layer_num=4, layer_dim=250, rng=StableRNG(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run!(opt, alg, obj, dynamics;savefile=false) \ No newline at end of file diff --git a/examples/Mopt_Kraus_multiparameter.jl b/examples/Mopt_Kraus_multiparameter.jl new file mode 100644 index 0000000..7e3aee7 --- /dev/null +++ b/examples/Mopt_Kraus_multiparameter.jl @@ -0,0 +1,59 @@ +using Random +using StableRNGs +include("../src/QuanEstimation.jl") + +# Kraus operators for the generalized qubit amplitude damping +n, p = 0.1, 0.1 +ψ0 = [1., 0] +ψ1 = [0., 1] +K0 =sqrt(1-n)*(ψ0*ψ0' + sqrt(1-p)*ψ1*ψ1') +K1 = sqrt(p-p*n)*ψ0*ψ1' +K2 = sqrt(n)*(sqrt(1-p)*ψ0*ψ0' + ψ1*ψ1') +K3 = sqrt(p*n)*ψ1*ψ0' +K = [K0, K1, K2, K3] + +dK0_n = -0.5*(ψ0*ψ0'+sqrt(1-p)*ψ1*ψ1')/sqrt(1-n) +dK1_n = -0.5*p*ψ0*ψ1'/sqrt(p-p*n) +dK2_n = 0.5*(sqrt(1-p)*ψ0*ψ0'+ψ1*ψ1')/sqrt(n) +dK3_n = 0.5*p*ψ1*ψ0'/sqrt(p*n) +dK0_p = -0.5*sqrt(1-n)*ψ1*ψ1'/sqrt(1-p) +dK1_p = 0.5*(1-n)*ψ0*ψ1'/sqrt(p-p*n) +dK2_p = -0.5*sqrt(n)*ψ0*ψ0'/sqrt(1-p) +dK3_p = -0.5*sqrt(n)*ψ0*ψ0'/sqrt(1-p) +dK3_p = 0.5*n*ψ1*ψ0'/sqrt(p*n) +dK = [[dK0_n, dK0_p], [dK1_n, dK1_p], [dK2_n, dK2_p], [dK3_n, dK3_p]] + +# measurement +dim = size(K[1],1) +POVM_basis = QuanEstimation.SIC(dim) +M_num = dim + +# Measurement optimization -- Projection +opt = QuanEstimation.Mopt() + +# Measurement optimization -- Linear combination +opt = QuanEstimation.Mopt(method=:LinearCombination, POVM_basis=POVM_basis, M_num=M_num) + +# Measurement optimization -- Rotation +opt = QuanEstimation.Mopt(method=:Rotation, POVM_basis=POVM_basis) + +# initial state +ρ₀ = 0.5*ones(2,2) + +# dynamics +dynamics = QuanEstimation.Kraus(opt, ρ₀, K, dK) + +# set objective +obj = QuanEstimation.CFIM_Obj() + +# control algorithm: AD +alg = QuanEstimation.AD(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) diff --git a/examples/Mopt_Kraus_single_parameter.jl b/examples/Mopt_Kraus_single_parameter.jl new file mode 100644 index 0000000..72e096a --- /dev/null +++ b/examples/Mopt_Kraus_single_parameter.jl @@ -0,0 +1,53 @@ +using Random +using StableRNGs +include("../src/QuanEstimation.jl") + +# Kraus operators for the amplitude damping channel +γ = 0.1 +K1 = [1. 0; 0 sqrt(1-γ)] +K2 = [0. sqrt(γ); 0 0] +K = [K1, K2] + +dK1 = [1. 0; 0 -sqrt(1-γ)/2] +dK2 = [0. sqrt(γ); 0 0] +dK = [[dK1], [dK2]] + +# measurement +dim = size(K[1],1) +POVM_basis = QuanEstimation.SIC(dim) +M_num = dim + +# Measurement optimization -- Projection +opt = QuanEstimation.Mopt() + +# Measurement optimization -- Linear combination +opt = QuanEstimation.Mopt(method=:LinearCombination, POVM_basis=POVM_basis, M_num=M_num) + +# Measurement optimization -- Rotation +opt = QuanEstimation.Mopt(method=:Rotation, POVM_basis=POVM_basis) + +# initial state +ρ₀ = 0.5*ones(2,2) + +# dynamics +dynamics = QuanEstimation.Kraus(opt, ρ₀, K, dK) + +# set objective +obj = QuanEstimation.CFIM_Obj() + +# control algorithm: AD +alg = QuanEstimation.AD(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: NM +alg = QuanEstimation.NM(max_episode=100, state_num=10, ar=1.0, ae=2.0, ac=0.5, as0=0.5,rng=MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + diff --git a/examples/Mopt_multiparameter.jl b/examples/Mopt_multiparameter.jl new file mode 100644 index 0000000..47679a9 --- /dev/null +++ b/examples/Mopt_multiparameter.jl @@ -0,0 +1,71 @@ +using Random +using LinearAlgebra +include("../src/QuanEstimation.jl") + +# initial state +ρ₀ = zeros(ComplexF64, 6, 6) +ρ₀[1:4:5, 1:4:5] .= 0.5 +dim = size(ρ₀, 1) + +# Hamiltonian +sx = [0. 1; 1 0] +sy = [0. -im; im 0] +sz = [1. 0; 0 -1] +s1 = [0. 1 0; 1 0 1; 0 1 0] / sqrt(2) +s2 = [0. -im 0; im 0 -im; 0 im 0] / sqrt(2) +s3 = [1. 0 0; 0 0 0; 0 0 -1] +Is = I1, I2, I3 =[ kron(I(3), sx), kron(I(3), sy), kron(I(3), sz)] +S = S1, S2, S3 = [kron(s1, I(2)), kron(s2, I(2)), kron(s3, I(2))] +B = B1, B2, B3 = [5.0e-4, 5.0e-4, 5.0e-4] +cons = 100 +D = (2pi * 2.87 * 1000) / cons +gS = (2pi * 28.03 * 1000) / cons +gI = (2pi * 4.32) / cons +A1 = (2pi * 3.65) / cons +A2 = (2pi * 3.03) / cons +H0 = sum([ + D*kron(s3^2, I(2)), + sum(gS * B .* S), + sum(gI * B .* Is), + A1 * (kron(s1, sx) + kron(s2, sy)), + A2 * kron(s3, sz) +]) +dH = gS * S + gI * Is +Hc = [S1, S2, S3] + +# dissipation +decay = [[S3, 2pi/cons]] + +# measurement +POVM_basis = [QuanEstimation.basis(dim, i)*QuanEstimation.basis(dim, i)' for i in 1:dim] +M_num = length(POVM_basis) + +#dynamics +tspan = range(0.0, 10.0, length=2500) + +# Measurement optimization -- Projection +opt = QuanEstimation.Mopt() + +# Measurement optimization -- Linear combination +opt = QuanEstimation.Mopt(method=:LinearCombination, POVM_basis=POVM_basis, M_num=M_num) + +# Measurement optimization -- Rotation +opt = QuanEstimation.Mopt(method=:Rotation, POVM_basis=POVM_basis) + +# set objective +obj = QuanEstimation.CFIM_Obj() + +# set dynamics +dynamics = QuanEstimation.Lindblad(opt, tspan ,ρ₀, H0, dH, decay) + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) +# +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# measurement optimization algorithm: AD +alg = QuanEstimation.AD(Adam=true, max_episode=100, ϵ=0.01, beta1=0.90, beta2=0.99) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) \ No newline at end of file diff --git a/examples/Mopt_single_parameter.jl b/examples/Mopt_single_parameter.jl new file mode 100644 index 0000000..7209efa --- /dev/null +++ b/examples/Mopt_single_parameter.jl @@ -0,0 +1,60 @@ +using Random +using StableRNGs +using LinearAlgebra +include("../src/QuanEstimation.jl") + +# initial state +ρ₀ = 0.5 * ones(2,2) +dim = size(ρ₀,1) +# Hamiltonian +ω₀ = 1.0 +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +H0 = 0.5 * ω₀ * sz +dH = [0.5 * sz] +# measurement +M_num = dim +rng = MersenneTwister(1234) +M= [ComplexF64[] for _ in 1:M_num] +for i in 1:M_num + r_ini = 2*rand(rng, dim) - ones(dim) + r = r_ini / norm(r_ini) + ϕ = 2pi*rand(rng, dim) + M[i] = [r*exp(im*ϕ) for (r,ϕ) in zip(r,ϕ)] +end +Measurement = QuanEstimation.gramschmidt(M) +POVM_basis = [m*m' for m in Measurement] +# dissipation +sp = [0 1; 0 0.0im] +sm = [0 0; 1 0.0im] +decay = [[sp, 0.0], [sm, 0.1]] +# dynamics +tspan = range(0.0, 10.0, length=2500) + +# Measurement optimization -- Projection +opt = QuanEstimation.Mopt(C=Measurement) + +# Measurement optimization -- Linear combination +opt = QuanEstimation.Mopt(method=:LinearCombination, POVM_basis=POVM_basis, M_num=M_num) + +# Measurement optimization -- Rotation +opt = QuanEstimation.Mopt(method=:Rotation, POVM_basis=POVM_basis) + +# set objective +obj = QuanEstimation.CFIM_Obj() + +# set dynamics +dynamics = QuanEstimation.Lindblad(opt, tspan ,ρ₀, H0, dH, decay) + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# measurement optimization algorithm: AD +alg = QuanEstimation.AD(Adam=true, max_episode=100, ϵ=0.01, beta1=0.90, beta2=0.99) +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) \ No newline at end of file diff --git a/examples/OBB.jl b/examples/OBB.jl new file mode 100644 index 0000000..7ab4210 --- /dev/null +++ b/examples/OBB.jl @@ -0,0 +1,52 @@ +using Trapz +include("../src/QuanEstimation.jl") + +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +# free Hamiltonian +B = pi/2.0 +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +function H0_func(x) + return 0.5*B*(sx*cos(x)+sz*sin(x)) +end +function dH_func(x) + return [0.5*B*(-sx*sin(x)+sz*cos(x))] +end +function d2H_func(x) + return [0.5*B*(-sx*cos(x)-sz*sin(x))] +end +# dynamics +decay_opt = [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])] +gamma = [0.0] +tspan = range(0.0, stop=1.0, length=1000) +# prior distribution +function p_func(x, mu, eta) + return exp(-(x-mu)^2/(2*eta^2))/(eta*sqrt(2*pi)) +end +function dp_func(x, mu, eta) + return -(x-mu)*exp(-(x-mu)^2/(2*eta^2))/(eta^3*sqrt(2*pi)) +end +x = [range(-pi/2.0, stop=pi/2.0, length=100)].|>Vector +mu, eta = 0.0, 0.5 +p_tp = [p_func(x[1][i], mu, eta) for i in 1:length(x[1])] +dp_tp = [dp_func(x[1][i], mu, eta) for i in 1:length(x[1])] +c = trapz(x[1], p_tp) +p, dp = p_tp/c, dp_tp/c + +rho = Vector{Matrix{ComplexF64}}(undef, length(x[1])) +drho = Vector{Vector{Matrix{ComplexF64}}}(undef, length(x[1])) +d2rho = Vector{Vector{Matrix{ComplexF64}}}(undef, length(x[1])) +for i = 1:length(x[1]) + H0_tp = H0_func(x[1][i]) + dH_tp = dH_func(x[1][i]) + d2H_tp = d2H_func(x[1][i]) + rho_tp, drho_tp, d2rho_tp = QuanEstimation.secondorder_derivative(H0_tp, dH_tp, d2H_tp, rho0, decay_opt, gamma, + [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])], [zeros(length(tspan)-1)],tspan) + rho[i] = rho_tp + drho[i] = drho_tp + d2rho[i] = d2rho_tp +end +f_OBB = QuanEstimation.OBB(x, p, dp, rho, drho, d2rho) +println(f_OBB) \ No newline at end of file diff --git a/examples/Sopt_Kraus_multiparameter.jl b/examples/Sopt_Kraus_multiparameter.jl new file mode 100644 index 0000000..d52ae5c --- /dev/null +++ b/examples/Sopt_Kraus_multiparameter.jl @@ -0,0 +1,68 @@ +using Random +using StableRNGs +include("../src/QuanEstimation.jl") + +# Kraus operators for the generalized qubit amplitude damping +n, p = 0.1, 0.1 +ψ0 = [1., 0] +ψ1 = [0., 1] +K0 =sqrt(1-n)*(ψ0*ψ0' + sqrt(1-p)*ψ1*ψ1') +K1 = sqrt(p-p*n)*ψ0*ψ1' +K2 = sqrt(n)*(sqrt(1-p)*ψ0*ψ0' + ψ1*ψ1') +K3 = sqrt(p*n)*ψ1*ψ0' +K = [K0, K1, K2, K3] + +dK0_n = -0.5*(ψ0*ψ0'+sqrt(1-p)*ψ1*ψ1')/sqrt(1-n) +dK1_n = -0.5*p*ψ0*ψ1'/sqrt(p-p*n) +dK2_n = 0.5*(sqrt(1-p)*ψ0*ψ0'+ψ1*ψ1')/sqrt(n) +dK3_n = 0.5*p*ψ1*ψ0'/sqrt(p*n) +dK0_p = -0.5*sqrt(1-n)*ψ1*ψ1'/sqrt(1-p) +dK1_p = 0.5*(1-n)*ψ0*ψ1'/sqrt(p-p*n) +dK2_p = -0.5*sqrt(n)*ψ0*ψ0'/sqrt(1-p) +dK3_p = -0.5*sqrt(n)*ψ0*ψ0'/sqrt(1-p) +dK3_p = 0.5*n*ψ1*ψ0'/sqrt(p*n) +dK = [[dK0_n, dK0_p], [dK1_n, dK1_p], [dK2_n, dK2_p], [dK3_n, dK3_p]] + +opt = QuanEstimation.Sopt() +dynamics = QuanEstimation.Kraus(opt, K, dK) + +# control algorithm: AD +alg = QuanEstimation.AD(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DDPG +alg = QuanEstimation.DDPG(max_episode=100, layer_num=4, layer_dim=250, rng=StableRNG(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: NM +alg = QuanEstimation.NM(max_episode=100, state_num=10, ar=1.0, ae=2.0, ac=0.5, as0=0.5,rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + diff --git a/examples/Sopt_Kraus_single_parameter.jl b/examples/Sopt_Kraus_single_parameter.jl new file mode 100644 index 0000000..0b8fd13 --- /dev/null +++ b/examples/Sopt_Kraus_single_parameter.jl @@ -0,0 +1,57 @@ +using Random +using StableRNGs +include("../src/QuanEstimation.jl") + +# Kraus operators for the amplitude damping channel +γ = 0.1 +K1 = [1. 0; 0 sqrt(1-γ)] +K2 = [0. sqrt(γ); 0 0] +K = [K1, K2] + +dK1 = [1. 0; 0 -sqrt(1-γ)/2] +dK2 = [0. sqrt(γ); 0 0] +dK = [[dK1], [dK2]] + +opt = QuanEstimation.Sopt() +dynamics = QuanEstimation.Kraus(opt, K, dK) + +# control algorithm: AD +alg = QuanEstimation.AD(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DDPG +alg = QuanEstimation.DDPG(max_episode=100, layer_num=4, layer_dim=250, rng=StableRNG(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: NM +alg = QuanEstimation.NM(max_episode=100, state_num=10, ar=1.0, ae=2.0, ac=0.5, as0=0.5,rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + diff --git a/examples/Sopt_multiparameter.jl b/examples/Sopt_multiparameter.jl new file mode 100644 index 0000000..3353b50 --- /dev/null +++ b/examples/Sopt_multiparameter.jl @@ -0,0 +1,84 @@ +using Random +using StableRNGs +using LinearAlgebra +using SparseArrays +include("../src/QuanEstimation.jl") + +N = 4 +# initial state +j, θ, ϕ = N÷2, 0.5pi, 0.5pi +Jp = Matrix(spdiagm(1=>[sqrt(j*(j+1)-m*(m+1)) for m in j:-1:-j][2:end])) +Jm = Jp' +ψ₀ = exp(0.5*θ*exp(im*ϕ)*Jm - 0.5*θ*exp(-im*ϕ)*Jp)*QuanEstimation.basis(Int(2*j+1), 1) +dim = length(ψ₀) +# free Hamiltonian +Λ = 1.0 +g = 0.5 +h = 0.1 + +Jx = 0.5*(Jp + Jm) +Jy = -0.5im*(Jp - Jm) +Jz = spdiagm(j:-1:-j) +H0 = -Λ*(Jx*Jx + g*Jy*Jy) / N + g * Jy^2 / N - h*Jz +dH = [-Λ*Jy*Jy/N, -Jz] +# measurement +M_num = N +rng = MersenneTwister(1234) +M_tp= [ComplexF64[] for _ in 1:M_num] +for i in 1:M_num + r_ini = 2*rand(rng, dim) - ones(dim) + r = r_ini / norm(r_ini) + φ = 2pi*rand(rng, dim) + M_tp[i] = [r*exp(im*φ) for (r,φ) in zip(r,φ)] +end +M_basis = QuanEstimation.gramschmidt(M_tp) +M = [m_basis*m_basis' for m_basis in M_basis] +# dissipation +decay = [[Jz, 0.1]] +# dynamics +tspan = range(0.0, 10.0, length=2500) + +W = [1/3 0.0;0.0 2/3] + +opt = QuanEstimation.Sopt(ψ₀=ψ₀) +dynamics = QuanEstimation.Lindblad(opt, tspan, H0, dH, decay) + +# control algorithm: AD +alg = QuanEstimation.AD(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +obj = QuanEstimation.QFIM_Obj(W=W) +obj = QuanEstimation.CFIM_Obj(M=M, W=W) +# obj = QuanEstimation.HCRB_Obj(W=W) + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DDPG +alg = QuanEstimation.DDPG(max_episode=100, layer_num=4, layer_dim=250, rng=StableRNG(1234)) +obj = QuanEstimation.QFIM_Obj(W=W) +obj = QuanEstimation.CFIM_Obj(M=M, W=W) +# obj = QuanEstimation.HCRB_Obj(W=W) + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj(W=W) +obj = QuanEstimation.CFIM_Obj(M=M, W=W) +# obj = QuanEstimation.HCRB_Obj(W=W) + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj(W=W) +obj = QuanEstimation.CFIM_Obj(M=M, W=W) +# obj = QuanEstimation.HCRB_Obj(W=W) + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: NM +alg = QuanEstimation.NM(max_episode=100, state_num=10, ar=1.0, ae=2.0, ac=0.5, as0=0.5,rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) diff --git a/examples/Sopt_single_parameter.jl b/examples/Sopt_single_parameter.jl new file mode 100644 index 0000000..33ada5d --- /dev/null +++ b/examples/Sopt_single_parameter.jl @@ -0,0 +1,83 @@ +using Random +using StableRNGs +using LinearAlgebra +using SparseArrays +include("../src/QuanEstimation.jl") + +N = 8 +# initial state +j, θ, ϕ = N÷2, 0.5pi, 0.5pi +Jp = Matrix(spdiagm(1=>[sqrt(j*(j+1)-m*(m+1)) for m in j:-1:-j][2:end])) +Jm = Jp' +ψ₀ = exp(0.5*θ*exp(im*ϕ)*Jm - 0.5*θ*exp(-im*ϕ)*Jp)*QuanEstimation.basis(Int(2*j+1), 1) +dim = length(ψ₀) +# free Hamiltonian +Λ = 1.0 +g = 0.5 +h = 0.1 + +Jx = 0.5*(Jp + Jm) +Jy = -0.5im*(Jp - Jm) +Jz = spdiagm(j:-1:-j) +H0 = -Λ*(Jx*Jx + g*Jy*Jy) / N - h*Jz +dH = [-Λ*Jy*Jy/N] +# measurement +M_num = N +rng = MersenneTwister(1234) +M_tp= [ComplexF64[] for _ in 1:M_num] +for i in 1:M_num + r_ini = 2*rand(rng, dim) - ones(dim) + r = r_ini / norm(r_ini) + φ = 2pi*rand(rng, dim) + M_tp[i] = [r*exp(im*φ) for (r,φ) in zip(r,φ)] +end +M_basis = QuanEstimation.gramschmidt(M_tp) +M = [m_basis*m_basis' for m_basis in M_basis] +# dissipation +decay = [[Jz, 0.1]] +# dynamics +tspan = range(0.0, 10.0, length=2500) + +opt = QuanEstimation.Sopt(ψ₀=ψ₀) +dynamics = QuanEstimation.Lindblad(opt, tspan, H0, dH, decay) + +# control algorithm: AD +alg = QuanEstimation.AD(Adam=true, max_episode=50, ϵ=0.01, beta1=0.90, beta2=0.99) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DDPG +alg = QuanEstimation.DDPG(max_episode=100, layer_num=4, layer_dim=250, rng=StableRNG(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: PSO +alg = QuanEstimation.PSO(max_episode=[100, 10], p_num=10, c0=1.0, c1=2.0, c2=2.0, rng= MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: DE +alg = QuanEstimation.DE(max_episode=100, p_num=10, c=1.0, cr=0.5, rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj(M=M) +# obj = QuanEstimation.HCRB_Obj() + + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) + +# control algorithm: NM +alg = QuanEstimation.NM(max_episode=100, state_num=10, ar=1.0, ae=2.0, ac=0.5, as0=0.5,rng=MersenneTwister(1234)) +obj = QuanEstimation.QFIM_Obj() +obj = QuanEstimation.CFIM_Obj() +# obj = QuanEstimation.HCRB_Obj() + +QuanEstimation.run(opt, alg, obj, dynamics;savefile=false) diff --git a/examples/adaptive_kraus_multiparameter.jl b/examples/adaptive_kraus_multiparameter.jl new file mode 100644 index 0000000..f56a274 --- /dev/null +++ b/examples/adaptive_kraus_multiparameter.jl @@ -0,0 +1,57 @@ +using Random +using StatsBase +include("../src/QuanEstimation.jl") + +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +psi0 = [1, 0] +psi1 = [0, 1] +function K_func(x) + n, p = x[1], x[2] + K0 = sqrt(1-n)*(psi0*psi0'+sqrt(1-p)*psi1*psi1') + K1 = sqrt(p-p*n)*psi0*psi1' + K2 = sqrt(n)*(sqrt(1-p)*psi0*psi0'+psi1*psi1') + K3 = sqrt(p*n)*psi1*psi0' + return [K0, K1, K2, K3] +end +function dK_func(x) + n, p = x[1], x[2] + dK0_n = -0.5*(psi0*psi0'+sqrt(1-p)*psi1*psi1')/sqrt(1-n) + dK1_n = -0.5*p*psi0*psi1'/sqrt(p-p*n) + dK2_n = 0.5*(sqrt(1-p)*psi0*psi0'+psi1*psi1')/sqrt(n) + dK3_n = 0.5*p*psi1*psi0'/sqrt(p*n) + dK0_p = -0.5*sqrt(1-n)*psi1*psi1'/sqrt(1-p) + dK1_p = 0.5*(1-n)*psi0*psi1'/sqrt(p-p*n) + dK2_p = -0.5*sqrt(n)*psi0*psi0'/sqrt(1-p) + dK3_p = -0.5*sqrt(n)*psi0*psi0'/sqrt(1-p) + dK3_p = 0.5*n*psi1*psi0'/sqrt(p*n) + return [[dK0_n, dK0_p], [dK1_n, dK1_p], [dK2_n, dK2_p], [dK3_n, dK3_p]] +end +# measurement +M1 = 0.5*[1.0+0.0im 1.0; 1.0 1.0] +M2 = 0.5*[1.0+0.0im -1.0; -1.0 1.0] +M = [M1, M2] +# prior distribution +x = [range(0.1, stop=0.9, length=100), range(0.1, stop=0.9, length=10)].|>Vector +p = (1.0/(x[1][end]-x[1][1]))*(1.0/(x[2][end]-x[2][1]))*ones((length(x[1]), length(x[2]))) + +rho = Matrix{Matrix{ComplexF64}}(undef, length.(x)...) +for i = 1:length(x[1]) + for j = 1:length(x[2]) + x_tp = [x[1][i], x[2][j]] + K_tp = K_func(x_tp) + rho[i,j] = [K * rho0 * K' for K in K_tp] |> sum + end +end +# Bayesian estimation +Random.seed!(1234) +y = [0 for i in 1:500] +res_rand = sample(1:length(y), 125, replace=false) +for i in 1:length(res_rand) + y[res_rand[i]] = 1 +end +pout, xout = QuanEstimation.Bayes(x, p, rho, y, M=M, savefile=false) +# adaptive +p = pout +K, dK = QuanEstimation.AdaptiveInput(x, K_func, dK_func; channel="kraus") +QuanEstimation.adaptive(x, p, rho0, K, dK, M=M, max_episode=10) diff --git a/examples/adaptive_kraus_singleparameter.jl b/examples/adaptive_kraus_singleparameter.jl new file mode 100644 index 0000000..d193043 --- /dev/null +++ b/examples/adaptive_kraus_singleparameter.jl @@ -0,0 +1,43 @@ +using Random +using StatsBase +include("../src/QuanEstimation.jl") + +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +# free Hamiltonian +function K_func(x) + K1 = [1.0 0.0; 0.0 sqrt(1 - x[1])] + K2 = [0.0 sqrt(x[1]); 0.0 0.0] + return [K1, K2] +end +function dK_func(x) + dK1 = [1.0 0.0; 0.0 -0.5 / sqrt(1 - x[1])] + dK2 = [0.0 0.5 / sqrt(x[1]); 0.0 0.0] + return [[dK1],[dK2]] +end +# measurement +M1 = 0.5*[1.0+0.0im 1.0; 1.0 1.0] +M2 = 0.5*[1.0+0.0im -1.0; -1.0 1.0] +M = [M1, M2] +#### prior distribution #### +x = [range(0.1, stop=0.9, length=100)].|>Vector +p = (1.0/(x[1][end]-x[1][1]))*ones(length(x[1])) + +rho = Vector{Matrix{ComplexF64}}(undef, length(x[1])) +for i = 1:length(x[1]) + K_tp = K_func([x[1][i]]) + rho[i] = [K * rho0 * K' for K in K_tp] |> sum +end + +# Bayesian estimation +Random.seed!(1234) +y = [0 for i in 1:500] +res_rand = sample(1:length(y), 125, replace=false) +for i in 1:length(res_rand) + y[res_rand[i]] = 1 +end +pout, xout = QuanEstimation.Bayes(x, p, rho, y, M=M, savefile=false) + +p = pout +K, dK = QuanEstimation.AdaptiveInput(x, K_func, dK_func; channel="kraus") +QuanEstimation.adaptive(x, p, rho0, K, dK; M=M, max_episode=10) diff --git a/examples/adaptive_multiparameter.jl b/examples/adaptive_multiparameter.jl new file mode 100644 index 0000000..c1115cc --- /dev/null +++ b/examples/adaptive_multiparameter.jl @@ -0,0 +1,51 @@ +using Random +using StatsBase +include("../src/QuanEstimation.jl") + +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +# free Hamiltonian +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +function H0_func(x) + return 0.5*x[2]*(sx*cos(x[1])+sz*sin(x[1])) +end +function dH_func(x) + return [0.5*x[2]*(-sx*sin(x[1])+sz*cos(x[1])), 0.5*(sx*cos(x[1])+sz*sin(x[1]))] +end +# measurement +M1 = 0.5*[1.0+0.0im 1.0; 1.0 1.0] +M2 = 0.5*[1.0+0.0im -1.0; -1.0 1.0] +M = [M1, M2] +# dynamics +tspan = range(0.0, stop=1.0, length=1000) +# dissipation +decay_opt = [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])] +gamma = [0.0] +# prior distribution +x = [range(0.0, stop=pi/2.0, length=100), range(pi/2-0.1, stop=pi/2+0.1, length=10)].|>Vector +p = (1.0/(x[1][end]-x[1][1]))*(1.0/(x[2][end]-x[2][1]))*ones((length(x[1]), length(x[2]))) + +rho = Matrix{Matrix{ComplexF64}}(undef, length.(x)...) +for i = 1:length(x[1]) + for j = 1:length(x[2]) + x_tp = [x[1][i], x[2][j]] + H0_tp = H0_func(x_tp) + dH_tp = dH_func(x_tp) + rho_tp, drho_tp = QuanEstimation.expm(H0_tp, dH_tp, [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])], [zeros(length(tspan)-1)],rho0, tspan,decay_opt, gamma) + rho[i,j] = rho_tp[end] + end +end +# Bayesian estimation +Random.seed!(1234) +y = [0 for i in 1:500] +res_rand = sample(1:length(y), 125, replace=false) +for i in 1:length(res_rand) + y[res_rand[i]] = 1 +end +pout, xout = QuanEstimation.Bayes(x, p, rho, y, M=M, savefile=false) +# adaptive +p = pout +H, dH = QuanEstimation.AdaptiveInput(x, H0_func, dH_func; channel="dynamics") +QuanEstimation.adaptive(x, p, rho0, tspan, H, dH, M=M, max_episode=10) diff --git a/examples/adaptive_singleparameter.jl b/examples/adaptive_singleparameter.jl new file mode 100644 index 0000000..0a2a26d --- /dev/null +++ b/examples/adaptive_singleparameter.jl @@ -0,0 +1,49 @@ +using Random +using StatsBase +include("../src/QuanEstimation.jl") + +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +# free Hamiltonian +B = pi/2.0 +sx = [0.0 1.0; 1.0 0.0im] +sy = [0.0 -1.0im; 1.0im 0.0] +sz = [1.0 0.0im; 0.0 -1.0] +function H0_func(x) + return 0.5*B*(sx*cos(x[1])+sz*sin(x[1])) +end +function dH_func(x) + return [0.5*B*(-sx*sin(x[1])+sz*cos(x[1]))] +end +# measurement +M1 = 0.5*[1.0+0.0im 1.0; 1.0 1.0] +M2 = 0.5*[1.0+0.0im -1.0; -1.0 1.0] +M = [M1, M2] +# dynamics +decay_opt = [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])] +gamma = [0.0] +tspan = range(0.0, stop=1.0, length=1000) |>Vector +#### prior distribution #### +x = [range(-0.25*pi+0.1, stop=3.0*pi/4.0-0.1, length=100)].|>Vector +p = (1.0/(x[1][end]-x[1][1]))*ones(length(x[1])) + +rho = Vector{Matrix{ComplexF64}}(undef, length(x[1])) +for i = 1:length(x[1]) + H0_tp = H0_func(x[1][i]) + dH_tp = dH_func(x[1][i]) + rho_tp, drho_tp = QuanEstimation.expm(H0_tp, dH_tp, [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])], [zeros(length(tspan)-1)], rho0, tspan, decay_opt, gamma) + rho[i] = rho_tp[end] +end + +# Bayesian estimation +Random.seed!(1234) +y = [0 for i in 1:500] +res_rand = sample(1:length(y), 125, replace=false) +for i in 1:length(res_rand) + y[res_rand[i]] = 1 +end +pout, xout = QuanEstimation.Bayes(x, p, rho, y, M=M, savefile=false) + +p = pout +H, dH = QuanEstimation.AdaptiveInput(x, H0_func, dH_func; channel="dynamics") +QuanEstimation.adaptive(x, p, rho0, tspan, H, dH; M=M, max_episode=10) diff --git a/examples/resources.jl b/examples/resources.jl new file mode 100644 index 0000000..3599d7f --- /dev/null +++ b/examples/resources.jl @@ -0,0 +1,40 @@ +using LinearAlgebra +include("../src/QuanEstimation.jl") + +#### spin squeezing #### + +rho_CSS = [0.25 -0.35355339im -0.25; 0.35355339im 0.5 -0.35355339im; -0.25 0.35355339im 0.25] + +xi = QuanEstimation.SpinSqueezing(rho_CSS; basis="Dicke", output="KU") +println(xi) + +#### Target time #### +# initial state +rho0 = 0.5*[1.0 1.0+0.0im; 1.0 1.0] +# free Hamiltonian +omega0 = 1.0 +sz = [1.0 0.0im; 0.0 -1.0] +H0 = 0.5 * omega0 * sz +dH = [0.5 * sz] +# measurement +M1 = 0.5*[1.0+0.0im 1.0; 1.0 1.0] +M2 = 0.5*[1.0+0.0im -1.0; -1.0 1.0] +M = [M1, M2] +# dissipation +sp = [0.0im 1.0; 0.0 0.0] +sm = [0.0im 0.0; 1.0 0.0] +decay_opt = [sp, sm] +gamma = [0.0, 0.1] +# dynamics +tspan = range(0.0, stop=50.0, length=2000) |>Vector +rho, drho = QuanEstimation.expm(H0, dH, [zeros(ComplexF64,size(rho0)[1],size(rho0)[1])], [zeros(length(tspan)-1)], rho0, tspan, decay_opt, gamma) +drho = [drho[i][1] for i in 1:2000] +QFI = [] +for ti in 2:2000 + QFI_tp = QuanEstimation.QFIM_SLD(rho[ti], drho[ti]) + append!(QFI, QFI_tp) +end +println(QFI[243]) +println(tspan[243]) +t = QuanEstimation.TargetTime(20.0, tspan, QuanEstimation.QFIM_SLD, rho, drho, eps=1e-8) +println(t) diff --git a/src/Project.toml b/src/Project.toml new file mode 100644 index 0000000..3567aca --- /dev/null +++ b/src/Project.toml @@ -0,0 +1,28 @@ +name = "QuanEstimation" +uuid = "088c8dff-a786-4a66-974c-03d3f6773f87" +authors = ["Hauiming Yu and contributors"] +version = "0.1.0" + +[deps] +DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" +BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" +Convex = "f65535da-76fb-5f13-bab9-19810c17039a" +Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +ReinforcementLearning = "158674fc-8238-5cab-b5ba-03dfc80d1318" +SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" + +[compat] +julia = "1.7.2" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"]