From 183dcdf270f281ec487103b763d1e9767623828b Mon Sep 17 00:00:00 2001 From: HuaiMing Date: Tue, 12 Apr 2022 09:42:44 +0800 Subject: [PATCH 1/2] wrappers --- examples/Bayesian_bounds_multiparameter.jl | 72 +++ examples/Bayesian_bounds_singleparameter.jl | 60 ++ .../Bayesian_estimation_multiparameter.jl | 50 ++ .../Bayesian_estimation_singleparameter.jl | 51 ++ examples/Comopt_Kraus_multiparameter.jl | 49 ++ examples/Comopt_Kraus_single_parameter.jl | 27 + examples/Comopt_multiparameter.jl | 89 +++ examples/Comopt_single_parameter.jl | 69 +++ examples/Copt_multi_parameter.jl | 102 ++++ examples/Copt_single_parameter.jl | 62 +++ examples/Mopt_Kraus_multiparameter.jl | 59 ++ examples/Mopt_Kraus_single_parameter.jl | 53 ++ examples/Mopt_multiparameter.jl | 71 +++ examples/Mopt_single_parameter.jl | 60 ++ examples/OBB.jl | 52 ++ examples/Sopt_Kraus_multiparameter.jl | 68 +++ examples/Sopt_Kraus_single_parameter.jl | 57 ++ examples/Sopt_multiparameter.jl | 84 +++ examples/Sopt_single_parameter.jl | 83 +++ examples/adaptive_kraus_multiparameter.jl | 57 ++ examples/adaptive_kraus_singleparameter.jl | 43 ++ examples/adaptive_multiparameter.jl | 51 ++ examples/adaptive_singleparameter.jl | 49 ++ examples/resources.jl | 40 ++ src/QuanEstimation.jl | 1 + src/algorithm/AD.jl | 9 +- src/algorithm/DE.jl | 39 +- src/algorithm/NM.jl | 3 + src/algorithm/PSO.jl | 39 +- src/algorithm/algorithm.jl | 126 ++--- src/common/common.jl | 10 +- src/dynamics/Kraus.jl | 60 -- src/dynamics/dynamics.jl | 4 +- src/dynamics/kraus/Kraus.jl | 23 + src/dynamics/kraus/KrausData.jl | 30 + src/dynamics/kraus/KrausDynamics.jl | 22 + src/dynamics/kraus/KrausWrapper.jl | 46 ++ src/dynamics/lindblad/Lindblad.jl | 31 ++ src/dynamics/lindblad/LindbladData.jl | 265 +++++++++ .../LindbladDynamics.jl} | 304 +---------- src/dynamics/lindblad/LindbladWrapper.jl | 512 ++++++++++++++++++ .../AsymptoticBound/AsymptoticBound.jl | 140 +++-- .../AsymptoticBound/AsymptoticBoundWrapper.jl | 39 ++ src/objective/AsymptoticBound/Holevo.jl | 4 +- src/optim/optim.jl | 97 ++-- src/output.jl | 4 +- src/run.jl | 14 + 47 files changed, 2765 insertions(+), 515 deletions(-) create mode 100644 examples/Bayesian_bounds_multiparameter.jl create mode 100644 examples/Bayesian_bounds_singleparameter.jl create mode 100644 examples/Bayesian_estimation_multiparameter.jl create mode 100644 examples/Bayesian_estimation_singleparameter.jl create mode 100644 examples/Comopt_Kraus_multiparameter.jl create mode 100644 examples/Comopt_Kraus_single_parameter.jl create mode 100644 examples/Comopt_multiparameter.jl create mode 100644 examples/Comopt_single_parameter.jl create mode 100644 examples/Copt_multi_parameter.jl create mode 100644 examples/Copt_single_parameter.jl create mode 100644 examples/Mopt_Kraus_multiparameter.jl create mode 100644 examples/Mopt_Kraus_single_parameter.jl create mode 100644 examples/Mopt_multiparameter.jl create mode 100644 examples/Mopt_single_parameter.jl create mode 100644 examples/OBB.jl create mode 100644 examples/Sopt_Kraus_multiparameter.jl create mode 100644 examples/Sopt_Kraus_single_parameter.jl create mode 100644 examples/Sopt_multiparameter.jl create mode 100644 examples/Sopt_single_parameter.jl create mode 100644 examples/adaptive_kraus_multiparameter.jl create mode 100644 examples/adaptive_kraus_singleparameter.jl create mode 100644 examples/adaptive_multiparameter.jl create mode 100644 examples/adaptive_singleparameter.jl create mode 100644 examples/resources.jl delete mode 100644 src/dynamics/Kraus.jl create mode 100644 src/dynamics/kraus/Kraus.jl create mode 100644 src/dynamics/kraus/KrausData.jl create mode 100644 src/dynamics/kraus/KrausDynamics.jl create mode 100644 src/dynamics/kraus/KrausWrapper.jl create mode 100644 src/dynamics/lindblad/Lindblad.jl create mode 100644 src/dynamics/lindblad/LindbladData.jl rename src/dynamics/{Lindblad.jl => lindblad/LindbladDynamics.jl} (71%) create mode 100644 src/dynamics/lindblad/LindbladWrapper.jl create mode 100644 src/objective/AsymptoticBound/AsymptoticBoundWrapper.jl 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/QuanEstimation.jl b/src/QuanEstimation.jl index ee1e6af..c19c5f4 100644 --- a/src/QuanEstimation.jl +++ b/src/QuanEstimation.jl @@ -16,6 +16,7 @@ using Interpolations const pkgpath = @__DIR__ const GLOBAL_RNG = MersenneTwister(1234) +const GLOBAL_EPS = 1e-8 include("optim/optim.jl") include("common/common.jl") include("dynamics/dynamics.jl") diff --git a/src/algorithm/AD.jl b/src/algorithm/AD.jl index 8685621..4568cb8 100644 --- a/src/algorithm/AD.jl +++ b/src/algorithm/AD.jl @@ -129,9 +129,12 @@ function update!(opt::Mopt_Rotation, alg::AbstractAD, obj, dynamics, output) dim = size(dynamics.data.ρ0)[1] M_num = length(POVM_basis) suN = suN_generator(dim) - append!(opt.Lambda, [Matrix{ComplexF64}(I,dim,dim)]) - append!(opt.Lambda, [suN[i] for i in 1:length(suN)]) - + if ismissing(opt.Lambda) + opt.Lambda = Matrix{ComplexF64}[] + append!(opt.Lambda, [Matrix{ComplexF64}(I,dim,dim)]) + append!(opt.Lambda, [suN[i] for i in 1:length(suN)]) + end + U = rotation_matrix(opt.s, opt.Lambda) M = [U*POVM_basis[i]*U' for i in 1:M_num] obj_QFIM = QFIM_Obj(obj) diff --git a/src/algorithm/DE.jl b/src/algorithm/DE.jl index c56b8c3..cca7654 100644 --- a/src/algorithm/DE.jl +++ b/src/algorithm/DE.jl @@ -1,6 +1,9 @@ #### control optimization #### function update!(opt::ControlOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg + if ismissing(ini_population) + ini_population = ([opt.ctrl,],) + end ini_population = ini_population[1] ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) @@ -76,6 +79,9 @@ end #### state optimization #### function update!(opt::StateOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg + if ismissing(ini_population) + ini_population = ([opt.ψ₀,],) + end ini_population = ini_population[1] dim = length(dynamics.data.ψ0) populations = repeat(dynamics, p_num) @@ -136,6 +142,9 @@ end #### projective measurement optimization #### function update!(opt::Mopt_Projection, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg + if ismissing(ini_population) + ini_population = ([opt.C], ) + end ini_population = ini_population[1] dim = size(dynamics.data.ρ0)[1] @@ -215,8 +224,11 @@ end #### find the optimal linear combination of a given set of POVM #### function update!(opt::Mopt_LinearComb, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg - ini_population = ini_population[1] (; B, POVM_basis, M_num) = opt + if ismissing(ini_population) + ini_population = ( [B], ) + end + ini_population = ini_population[1] basis_num = length(POVM_basis) populations = [[zeros(basis_num) for j in 1:M_num] for i in 1:p_num] @@ -295,12 +307,19 @@ end #### find the optimal rotated measurement of a given set of POVM #### function update!(opt::Mopt_Rotation, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg - ini_population = ini_population[1] (; s, POVM_basis, Lambda) = opt + if ismissing(ini_population) + ini_population = ([s,],) + end + ini_population = ini_population[1] dim = size(dynamics.data.ρ0)[1] suN = suN_generator(dim) - append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) - append!(Lambda, [suN[i] for i in 1:length(suN)]) + if ismissing(Lambda) + Lambda = Matrix{ComplexF64}[] + append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) + append!(Lambda, [suN[i] for i in 1:length(suN)]) + end + M_num = length(POVM_basis) populations = [zeros(dim^2) for i in 1:p_num] # initialization @@ -377,6 +396,9 @@ end #### state and control optimization #### function update!(opt::StateControlOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg + if ismissing(ini_population) + ini_population = ([opt.ψ₀], [opt.ctrl,]) + end psi0, ctrl0 = ini_population ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) @@ -473,6 +495,9 @@ end #### state and measurement optimization #### function update!(opt::StateMeasurementOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg + if ismissing(ini_population) + ini_population = ( [opt.ψ₀], [opt.C,]) + end psi0, measurement0 = ini_population dim = length(dynamics.data.ψ0) M_num = length(opt.C) @@ -573,6 +598,9 @@ end #### control and measurement optimization #### function update!(opt::ControlMeasurementOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg + if ismissing(ini_population) + ini_population = ([opt.ctrl,], [opt.C]) + end ctrl0, measurement0 = ini_population dim = size(dynamics.data.ρ0)[1] ctrl_length = length(dynamics.data.ctrl[1]) @@ -684,6 +712,9 @@ end #### state, control and measurement optimization #### function update!(opt::StateControlMeasurementOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg + if ismissing(ini_population) + ini_population = ([opt.ψ₀], [opt.ctrl,], [opt.C]) + end psi0, ctrl0, measurement0 = ini_population dim = length(dynamics.data.ψ0) ctrl_length = length(dynamics.data.ctrl[1]) diff --git a/src/algorithm/NM.jl b/src/algorithm/NM.jl index e423d7a..12f6b0b 100644 --- a/src/algorithm/NM.jl +++ b/src/algorithm/NM.jl @@ -1,5 +1,8 @@ function update!(opt::StateOpt, alg::NM, obj, dynamics, output) (; max_episode, state_num, ini_state, ar, ae, ac, as0, rng) = alg + if ismissing(ini_state) + ini_state = [opt.ψ₀] + end dim = length(dynamics.data.ψ0) nelder_mead = repeat(dynamics, state_num) diff --git a/src/algorithm/PSO.jl b/src/algorithm/PSO.jl index 28693e4..a32a8c1 100644 --- a/src/algorithm/PSO.jl +++ b/src/algorithm/PSO.jl @@ -1,6 +1,9 @@ #### control optimization #### function update!(opt::ControlOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg + if ismissing(ini_particle) + ini_particle = ([opt.ctrl,],) + end ini_particle = ini_particle[1] ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) @@ -100,6 +103,9 @@ end #### state optimization #### function update!(opt::StateOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg + if ismissing(ini_particle) + ini_particle = ([opt.ψ₀], ) + end ini_particle = ini_particle[1] dim = length(dynamics.data.ψ0) particles = repeat(dynamics, p_num) @@ -174,6 +180,9 @@ end #### projective measurement optimization #### function update!(opt::Mopt_Projection, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg + if ismissing(ini_particle) + ini_particle = ([opt.C], ) + end ini_particle = ini_particle[1] dim = size(dynamics.data.ρ0)[1] M_num = length(opt.C) @@ -197,7 +206,7 @@ function update!(opt::Mopt_Projection, alg::PSO, obj, dynamics, output) obj_copy = set_M(obj, M) f_ini, f_comp = objective(obj_copy, dynamics) - obj_QFIM = QFIM_Obj(obj) + obj_QFIM = QFIM_Obj(obj) f_opt, f_comp = objective(obj_QFIM, dynamics) set_f!(output, f_ini) @@ -264,9 +273,12 @@ end #### find the optimal linear combination of a given set of POVM #### function update!(opt::Mopt_LinearComb, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg - ini_particle = ini_particle[1] (; B, POVM_basis, M_num) = opt basis_num = length(POVM_basis) + if ismissing(ini_particle) + ini_particle = ( [B], ) + end + ini_particle = ini_particle[1] particles = [[zeros(basis_num) for j in 1:M_num] for i in 1:p_num] if typeof(max_episode) == Int @@ -356,13 +368,20 @@ end #### find the optimal rotated measurement of a given set of POVM #### function update!(opt::Mopt_Rotation, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg + if ismissing(ini_particle) + ini_particle = ([opt.s], ) + end ini_particle = ini_particle[1] (; s, POVM_basis, Lambda) = opt M_num = length(POVM_basis) dim = size(dynamics.data.ρ0)[1] suN = suN_generator(dim) - append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) - append!(Lambda, [suN[i] for i in 1:length(suN)]) + if ismissing(Lambda) + Lambda = Matrix{ComplexF64}[] + append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) + append!(Lambda, [suN[i] for i in 1:length(suN)]) + end + particles = repeat(s, p_num) if typeof(max_episode) == Int @@ -448,6 +467,9 @@ end #### state and control optimization #### function update!(opt::StateControlOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg + if ismissing(ini_particle) + ini_particle = ([opt.ψ₀], [opt.ctrl,]) + end psi0, ctrl0 = ini_particle dim = length(dynamics.data.ψ0) ctrl_length = length(dynamics.data.ctrl[1]) @@ -556,6 +578,9 @@ end #### state and measurement optimization #### function update!(opt::StateMeasurementOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg + if ismissing(ini_particle) + ini_particle = ([opt.ψ₀], [opt.C]) + end psi0, measurement0 = ini_particle dim = length(dynamics.data.ψ0) M_num = length(opt.C) @@ -668,6 +693,9 @@ end #### control and measurement optimization #### function update!(opt::ControlMeasurementOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg + if ismissing(ini_particle) + ini_particle = ([opt.ctrl,], [opt.C]) + end ctrl0, measurement0 = ini_particle ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) @@ -788,6 +816,9 @@ end #### state, control and measurement optimization #### function update!(opt::StateControlMeasurementOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg + if ismissing(ini_particle) + ini_particle = ([opt.ψ₀], [opt.ctrl,], [opt.C]) + end psi0, ctrl0, measurement0 = ini_particle ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) diff --git a/src/algorithm/algorithm.jl b/src/algorithm/algorithm.jl index 845bc03..7656120 100644 --- a/src/algorithm/algorithm.jl +++ b/src/algorithm/algorithm.jl @@ -1,99 +1,103 @@ abstract type AbstractAlgorithm end abstract type AbstractGRAPE <: AbstractAlgorithm end -struct GRAPE <: AbstractGRAPE - max_episode::Number - ϵ::Number +Base.@kwdef struct GRAPE <: AbstractGRAPE + max_episode::Int = 300 + ϵ::Number = 0.01 end -struct GRAPE_Adam <: AbstractGRAPE - max_episode::Number - ϵ::Number - beta1::Number - beta2::Number +Base.@kwdef struct GRAPE_Adam <: AbstractGRAPE + max_episode::Int = 300 + ϵ::Number = 0.01 + beta1::Number = 0.90 + beta2::Number = 0.99 end GRAPE(max_episode, ϵ, beta1, beta2) = GRAPE_Adam(max_episode, ϵ, beta1, beta2) +GRAPE(;max_episode=300, ϵ=0.01, beta1=0.90, beta2=0.99, Adam::Bool=true) = Adam ? GRAPE_Adam(max_episode, ϵ, beta1, beta2) : GRAPE(max_episode, ϵ) abstract type AbstractautoGRAPE <: AbstractAlgorithm end -struct autoGRAPE <: AbstractautoGRAPE - max_episode::Number - ϵ::Number +Base.@kwdef struct autoGRAPE <: AbstractautoGRAPE + max_episode::Int = 300 + ϵ::Number = 0.01 end -struct autoGRAPE_Adam <: AbstractautoGRAPE - max_episode::Number - ϵ::Number - beta1::Number - beta2::Number +Base.@kwdef struct autoGRAPE_Adam <: AbstractautoGRAPE + max_episode::Int = 300 + ϵ::Number = 0.01 + beta1::Number = 0.90 + beta2::Number = 0.99 end autoGRAPE(max_episode, ϵ, beta1, beta2) = autoGRAPE_Adam(max_episode, ϵ, beta1, beta2) +autoGRAPE(;max_episode=300, ϵ=0.01, beta1=0.90, beta2=0.99, Adam::Bool=true) = Adam ? autoGRAPE_Adam(max_episode, ϵ, beta1, beta2) : autoGRAPE(max_episode, ϵ) abstract type AbstractAD <: AbstractAlgorithm end -struct AD <: AbstractAD - max_episode::Number - ϵ::Number +Base.@kwdef struct AD <: AbstractAD + max_episode::Number = 300 + ϵ::Number = 0.01 end -struct AD_Adam <: AbstractAD - max_episode::Number - ϵ::Number - beta1::Number - beta2::Number +Base.@kwdef struct AD_Adam <: AbstractAD + max_episode::Number = 300 + ϵ::Number = 0.01 + beta1::Number = 0.90 + beta2::Number = 0.99 end -AD(max_episode, ϵ, beta1, beta2) = AD_Adam(max_episode, ϵ, beta1, beta2) - -struct PSO <: AbstractAlgorithm - max_episode::Union{T,Vector{T}} where {T<:Number} - p_num::Number - ini_particle::Tuple - c0::Number - c1::Number - c2::Number - rng::AbstractRNG +AD(max_episode, ϵ, beta1, beta2) = AD_Adam(max_episode=max_episode, ϵ=ϵ, beta1=beta1, beta2=beta2) +AD(;max_episode=300, ϵ=0.01, beta1=0.90, beta2=0.99, Adam::Bool=true) = Adam ? AD_Adam(max_episode, ϵ, beta1, beta2) : AD(max_episode, ϵ) + +Base.@kwdef struct PSO <: AbstractAlgorithm + max_episode::Union{T,Vector{T}} where {T<:Int} = [1000, 100] + p_num::Int = 10 + ini_particle::Union{Tuple, Missing} = missing + c0::Number = 1.0 + c1::Number = 2.0 + c2::Number = 2.0 + rng::AbstractRNG = GLOBAL_RNG end PSO(max_episode, p_num, ini_particle, c0, c1, c2) = PSO(max_episode, p_num, ini_particle, c0, c1, c2, GLOBAL_RNG) PSO(max_episode, p_num, ini_particle, c0, c1, c2, seed::Number) = - PSO(max_episode, p_num, ini_particle, c0, c1, c2, StableRNG(seed)) - -struct DE <: AbstractAlgorithm - max_episode::Number - p_num::Number - ini_population::Tuple - c::Number - cr::Number - rng::AbstractRNG + PSO(max_episode, p_num, ini_particle, c0, c1, c2, MersenneTwister(seed)) + +Base.@kwdef struct DE <: AbstractAlgorithm + max_episode::Number = 1000 + p_num::Number = 10 + ini_population::Union{Tuple, Missing} = missing + c::Number = 1.0 + cr::Number = 0.5 + rng::AbstractRNG = GLOBAL_RNG end -DE(max_episode, p_num, c, cr) = DE(max_episode, p_num, c, cr, GLOBAL_RNG) +DE(max_episode, p_num, ini_population, c, cr) = + DE(max_episode, p_num, ini_population, c, cr, GLOBAL_RNG) DE(max_episode, p_num, ini_population, c, cr, seed::Number) = DE(max_episode, p_num, ini_population, c, cr, MersenneTwister(seed)) -struct DDPG <: AbstractAlgorithm - max_episode::Number - layer_num::Number - layer_dim::Number - rng::AbstractRNG +Base.@kwdef struct DDPG <: AbstractAlgorithm + max_episode::Int = 500 + layer_num::Int = 3 + layer_dim::Int = 200 + rng::AbstractRNG = GLOBAL_RNG end DDPG(max_episode, layer_num, layer_dim) = DDPG(max_episode, layer_num, layer_dim, GLOBAL_RNG) DDPG(max_episode, layer_num, layer_dim, seed::Number) = - DDPG(max_episode, layer_num, layer_dim, MersenneTwister(seed)) - -struct NM <: AbstractAlgorithm - max_episode::Number - state_num::Number - ini_state::AbstractVector - ar::Number - ae::Number - ac::Number - as0::Number - rng::AbstractRNG + DDPG(max_episode, layer_num, layer_dim, StableRNG(seed)) + +Base.@kwdef struct NM <: AbstractAlgorithm + max_episode::Int = 1000 + state_num::Int = 10 + ini_state::Union{AbstractVector, Missing} = missing + ar::Number = 1.0 + ae::Number = 2.0 + ac::Number = 0.5 + as0::Number = 0.5 + rng::AbstractRNG = GLOBAL_RNG end NM(max_episode, state_num, nelder_mead, ar, ae, ac, as0) = @@ -117,4 +121,4 @@ include("DDPG.jl") include("DE.jl") include("GRAPE.jl") include("NM.jl") -include("PSO.jl") +include("PSO.jl") \ No newline at end of file diff --git a/src/common/common.jl b/src/common/common.jl index 914f615..3011bce 100644 --- a/src/common/common.jl +++ b/src/common/common.jl @@ -24,7 +24,7 @@ function Base.repeat(system, N) end function Base.repeat(system, M, N) - reshape(repeat(system, M * N), M, N) + reshape(repeat(system, M * N), M,N) end function filterZeros!(x::Matrix{T}) where {T <: Complex} @@ -277,7 +277,7 @@ function initial_ctrl!(opt, ctrl0, dynamics, p_num, rng) ctrl0 = [ctrl0[i] for i in 1:p_num] end for pj in 1:length(ctrl0) - dynamics[pj].data.ctrl = [[ctrl0[pj][i,j] for j in 1:ctrl_length] for i in 1:ctrl_num] + dynamics[pj].data.ctrl = deepcopy(ctrl0[pj]) #[[ctrl0[pj][i, j] for j in 1:ctrl_length] for i in 1:ctrl_num] end if opt.ctrl_bound[1] == -Inf || opt.ctrl_bound[2] == Inf for pj in (length(ctrl0)+1):p_num @@ -310,7 +310,7 @@ function initial_M!(measurement0, C_all, dim, p_num, M_num, rng) measurement0 = [measurement0[i] for i in 1:p_num] end for pj in 1:length(measurement0) - C_all[pj] = [[measurement0[pj][i,j] for j in 1:dim] for i in 1:M_num] + C_all[pj] = deepcopy(measurement0[pj]) #[[measurement0[pj][i,j] for j in 1:dim] for i in 1:M_num] end for pj in (length(measurement0)+1):p_num M_tp = [Vector{ComplexF64}(undef, dim) for i in 1:M_num] @@ -331,7 +331,7 @@ function initial_LinearComb!(measurement0, B_all, basis_num, M_num, p_num, rng) measurement0 = [measurement0[i] for i in 1:p_num] end for pj in 1:length(measurement0) - B_all[pj] = [[measurement0[pj][i,j] for j in 1:basis_num] for i in 1:M_num] + B_all[pj] = deepcopy(measurement0[pj]) #[[measurement0[pj][i,j] for j in 1:basis_num] for i in 1:M_num] end for pj in (length(measurement0)+1):p_num @@ -358,7 +358,7 @@ function initial_LinearComb!(measurement0, B_all, basis_num, M_num, p_num, rng) measurement0 = [measurement0[i] for i in 1:p_num] end for pj in 1:length(measurement0) - B_all[pj] = [[measurement0[pj][i,j] for j in 1:dim] for i in 1:M_num] + B_all[pj] = deepcopy(measurement0[pj]) #[[measurement0[pj][i,j] for j in 1:basis_num] for i in 1:M_num] end for pj in (length(measurement0)+1):p_num diff --git a/src/dynamics/Kraus.jl b/src/dynamics/Kraus.jl deleted file mode 100644 index 06c3de3..0000000 --- a/src/dynamics/Kraus.jl +++ /dev/null @@ -1,60 +0,0 @@ -# dynamics in Kraus rep. -struct Kraus{R} <: AbstractDynamics - data::AbstractDynamicsData - noise_type::Symbol - ctrl_type::Symbol - state_rep::Symbol -end - -mutable struct Kraus_dm <: AbstractDynamicsData - K::AbstractVector - dK::AbstractVecOrMat - ρ0::AbstractMatrix -end - -mutable struct Kraus_pure <: AbstractDynamicsData - K::AbstractVector - dK::AbstractVecOrMat - ψ0::AbstractVector -end - -# Constructor for Kraus dynamics -Kraus(K::AbstractVector, dK::AbstractVector, ρ0::AbstractMatrix) = - Kraus{dm}(Kraus_dm(K, dK, ρ0), :noiseless, :free, :dm) -Kraus(K::AbstractVector, dK::AbstractVector, ψ0::AbstractVector) = - Kraus{ket}(Kraus_pure(K, dK, ψ0), :noiseless, :free, :ket) -Kraus(K::AbstractVector, dK::AbstractMatrix, ρ0::AbstractMatrix) = - Kraus{dm}(Kraus_dm(K, [[dK[i,j,:,:][1] for j in 1:size(dK,2)] for i in 1:size(dK,1)], ρ0), :noiseless, :free, :dm) -Kraus(K::AbstractVector, dK::AbstractMatrix, ψ0::AbstractVector) = - Kraus{ket}(Kraus_pure(K,[[dK[i,j,:,:][1] for j in 1:size(dK,2)] for i in 1:size(dK,1)], ψ0), :noiseless, :free, :ket) - -function set_state(dynamics::Kraus, state::AbstractVector) - temp = deepcopy(dynamics) - temp.data.ψ0 = state - temp -end - -function set_state(dynamics::Kraus, state::AbstractMatrix) - temp = deepcopy(dynamics) - temp.data.ρ0 = state - temp -end - -#### evolution of pure states under time-independent Hamiltonian without noise and controls #### -function evolve(dynamics::Kraus{ket}) - (; K, dK, ψ0) = dynamics.data - ρ0 = ψ0 * ψ0' - ρ = [K * ρ0 * K' for K in K] |> sum - dρ = [[dK * ρ0 * K' + K * ρ0 * dK' for (K,dK) in zip(K,dK)] |> sum for dK in dK] - - ρ, dρ -end - -#### evolution of density matrix under time-independent Hamiltonian without noise and controls #### -function evolve(dynamics::Kraus{dm}) - (; K, dK, ρ0) = dynamics.data - ρ = [K * ρ0 * K' for K in K] |> sum - dρ = [[dK * ρ0 * K' + K * ρ0 * dK' for (K,dK) in zip(K,dK)] |> sum for dK in dK] - - ρ, dρ -end diff --git a/src/dynamics/dynamics.jl b/src/dynamics/dynamics.jl index 5f5e66a..91b73aa 100644 --- a/src/dynamics/dynamics.jl +++ b/src/dynamics/dynamics.jl @@ -24,5 +24,5 @@ isCtrl(::free) = false isCtrl(::controlled) = true isCtrl(dynamics::AbstractDynamics) = dynamics.ctrl_type |> eval |> isCtrl -include("Lindblad.jl") -include("Kraus.jl") +include("lindblad/Lindblad.jl") +include("kraus/Kraus.jl") diff --git a/src/dynamics/kraus/Kraus.jl b/src/dynamics/kraus/Kraus.jl new file mode 100644 index 0000000..c76514a --- /dev/null +++ b/src/dynamics/kraus/Kraus.jl @@ -0,0 +1,23 @@ +# dynamics in Kraus rep. +struct Kraus{R} <: AbstractDynamics + data::AbstractDynamicsData + noise_type::Symbol + ctrl_type::Symbol + state_rep::Symbol +end + +include("KrausData.jl") +include("KrausDynamics.jl") +include("KrausWrapper.jl") + +function set_state(dynamics::Kraus, state::AbstractVector) + temp = deepcopy(dynamics) + temp.data.ψ0 = state + temp +end + +function set_state(dynamics::Kraus, state::AbstractMatrix) + temp = deepcopy(dynamics) + temp.data.ρ0 = state + temp +end diff --git a/src/dynamics/kraus/KrausData.jl b/src/dynamics/kraus/KrausData.jl new file mode 100644 index 0000000..d0c4059 --- /dev/null +++ b/src/dynamics/kraus/KrausData.jl @@ -0,0 +1,30 @@ +abstract type KrausDynamicsData <: AbstractDynamicsData end + +mutable struct Kraus_dm <: KrausDynamicsData + K::AbstractVector + dK::AbstractVecOrMat + ρ0::AbstractMatrix +end + +mutable struct Kraus_pure <: KrausDynamicsData + K::AbstractVector + dK::AbstractVecOrMat + ψ0::AbstractVector +end + +# Constructor for Kraus dynamics +Kraus(K::AbstractVector, dK::AbstractVector, ρ0::AbstractMatrix) = + Kraus{dm}(Kraus_dm(K, dK, ρ0), :noiseless, :free, :dm) +Kraus(K::AbstractVector, dK::AbstractVector, ψ0::AbstractVector) = + Kraus{ket}(Kraus_pure(K, dK, ψ0), :noiseless, :free, :ket) +Kraus(K::AbstractVector, dK::AbstractMatrix, ρ0::AbstractMatrix) = + Kraus{dm}(Kraus_dm(K, [[dK[i,j,:,:][1] for j in 1:size(dK,2)] for i in 1:size(dK,1)], ρ0), :noiseless, :free, :dm) +Kraus(K::AbstractVector, dK::AbstractMatrix, ψ0::AbstractVector) = + Kraus{ket}(Kraus_pure(K,[[dK[i,j,:,:][1] for j in 1:size(dK,2)] for i in 1:size(dK,1)], ψ0), :noiseless, :free, :ket) + +para_type(data::KrausDynamicsData) = length(data.dK[1]) == 1 ? :single_para : :multi_para + +get_dim(k::Kraus_dm) = size(k.ρ0, 1) +get_dim(k::Kraus_pure) = size(k.ψ0, 1) + +get_para(k::KrausDynamicsData) = length(k.dK[1]) diff --git a/src/dynamics/kraus/KrausDynamics.jl b/src/dynamics/kraus/KrausDynamics.jl new file mode 100644 index 0000000..4056f95 --- /dev/null +++ b/src/dynamics/kraus/KrausDynamics.jl @@ -0,0 +1,22 @@ +#### evolution of pure states under time-independent Hamiltonian without noise and controls #### +function evolve(dynamics::Kraus{ket}) + (; K, dK, ψ0) = dynamics.data + ρ0 = ψ0 * ψ0' + K_num = length(K) + para_num = length(dK[1]) + ρ = [K[i] * ρ0 * K[i]' for i in 1:K_num] |> sum + dρ = [[dK[i][j] * ρ0 * K[i]' + K[i] * ρ0 * dK[i][j]' for i in 1:K_num] |> sum for j in 1:para_num] + + ρ, dρ +end + +#### evolution of density matrix under time-independent Hamiltonian without noise and controls #### +function evolve(dynamics::Kraus{dm}) + (; K, dK, ρ0) = dynamics.data + K_num = length(K) + para_num = length(dK[1]) + ρ = [K[i] * ρ0 * K[i]' for i in 1:K_num] |> sum + dρ = [[dK[i][j] * ρ0 * K[i]' + K[i] * ρ0 * dK[i][j]' for i in 1:K_num] |> sum for j in 1:para_num] + + ρ, dρ +end diff --git a/src/dynamics/kraus/KrausWrapper.jl b/src/dynamics/kraus/KrausWrapper.jl new file mode 100644 index 0000000..5f14516 --- /dev/null +++ b/src/dynamics/kraus/KrausWrapper.jl @@ -0,0 +1,46 @@ +function Kraus(opt::StateOpt, K, dK;rng=GLOBAL_RNG) + (;ψ₀) = opt + dim = size(K[1], 1) + if ismissing(ψ₀) + r_ini = 2*rand(rng, dim) - ones(dim) + r = r_ini ./ norm(r_ini) + ϕ = 2pi*rand(rng, dim) + ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.ψ₀ = ψ₀ + end + + K = complex.(K) + dK = [complex.(dk) for dk in dK] + ψ₀ = complex(ψ₀) + + Kraus(K, dK, ψ₀) +end + +function Kraus(opt::AbstractMopt, ρ₀::AbstractMatrix, K, dK;rng=GLOBAL_RNG, eps=GLOBAL_EPS) + dim = size(ρ₀, 1) + _ini_measurement!(opt, dim, rng; eps=eps) + + K = complex.(K) + dK = [complex.(dk) for dk in dK] + ρ₀ = complex(ρ₀) + Kraus(K, dK, ρ₀) +end + + +function Kraus(opt::CompOpt, K, dK;rng=GLOBAL_RNG, eps=GLOBAL_EPS) + (;ψ₀) = opt + dim = size(K[1], 1) + if ismissing(ψ₀) + r_ini = 2*rand(rng, dim) - ones(dim) + r = r_ini ./ norm(r_ini) + ϕ = 2pi*rand(rng, dim) + ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.ψ₀ = ψ₀ + end + + _ini_measurement!(opt, dim, rng; eps=eps) + K = complex.(K) + dK = [complex.(dk) for dk in dK] + ψ₀ = complex(ψ₀) + Kraus(K, dK, ψ₀) +end \ No newline at end of file diff --git a/src/dynamics/lindblad/Lindblad.jl b/src/dynamics/lindblad/Lindblad.jl new file mode 100644 index 0000000..361f393 --- /dev/null +++ b/src/dynamics/lindblad/Lindblad.jl @@ -0,0 +1,31 @@ + +# dynamics in Lindblad form +struct Lindblad{N,C,R,P} <: AbstractDynamics + data::AbstractDynamicsData + noise_type::Symbol + ctrl_type::Symbol + state_rep::Symbol + para_type::Symbol +end + +include("LindbladData.jl") +include("LindbladDynamics.jl") +include("LindbladWrapper.jl") + +function set_ctrl(dynamics::Lindblad, ctrl) + temp = deepcopy(dynamics) + temp.data.ctrl = ctrl + temp +end + +function set_state(dynamics::Lindblad, state::AbstractVector) + temp = deepcopy(dynamics) + temp.data.ψ0 = state + temp +end + +function set_state(dynamics::Lindblad, state::AbstractMatrix) + temp = deepcopy(dynamics) + temp.data.ρ0 = state + temp +end \ No newline at end of file diff --git a/src/dynamics/lindblad/LindbladData.jl b/src/dynamics/lindblad/LindbladData.jl new file mode 100644 index 0000000..06ca73c --- /dev/null +++ b/src/dynamics/lindblad/LindbladData.jl @@ -0,0 +1,265 @@ +abstract type LindbladDynamicsData <: AbstractDynamicsData end + + +mutable struct Lindblad_noiseless_free <: LindbladDynamicsData + H0::AbstractMatrix + dH::AbstractVector + ρ0::AbstractMatrix + tspan::AbstractVector +end + +mutable struct Lindblad_noisy_free <: LindbladDynamicsData + H0::AbstractMatrix + dH::AbstractVector + ρ0::AbstractMatrix + tspan::AbstractVector + decay_opt::AbstractVector + γ::AbstractVector +end + + +mutable struct Lindblad_noiseless_timedepend <: LindbladDynamicsData + H0::AbstractVector + dH::AbstractVector + ρ0::AbstractMatrix + tspan::AbstractVector +end + +mutable struct Lindblad_noisy_timedepend <: LindbladDynamicsData + H0::AbstractVector + dH::AbstractVector + ρ0::AbstractMatrix + tspan::AbstractVector + decay_opt::AbstractVector + γ::AbstractVector +end + +mutable struct Lindblad_noiseless_controlled <: LindbladDynamicsData + H0::AbstractVecOrMat + dH::AbstractVector + ρ0::AbstractMatrix + tspan::AbstractVector + Hc::AbstractVector + ctrl::AbstractVector +end + +mutable struct Lindblad_noisy_controlled <: LindbladDynamicsData + H0::AbstractVecOrMat + dH::AbstractVector + ρ0::AbstractMatrix + tspan::AbstractVector + decay_opt::AbstractVector + γ::AbstractVector + Hc::AbstractVector + ctrl::AbstractVector +end + + +mutable struct Lindblad_noiseless_free_pure <: LindbladDynamicsData + H0::AbstractMatrix + dH::AbstractVector + ψ0::AbstractVector + tspan::AbstractVector +end + +mutable struct Lindblad_noisy_free_pure <: LindbladDynamicsData + H0::AbstractMatrix + dH::AbstractVector + ψ0::AbstractVector + tspan::AbstractVector + decay_opt::AbstractVector + γ::AbstractVector +end + +mutable struct Lindblad_noiseless_timedepend_pure <: LindbladDynamicsData + H0::AbstractVector + dH::AbstractVector + ψ0::AbstractVector + tspan::AbstractVector +end + +mutable struct Lindblad_noisy_timedepend_pure <: LindbladDynamicsData + H0::AbstractVector + dH::AbstractVector + ψ0::AbstractVector + tspan::AbstractVector + decay_opt::AbstractVector + γ::AbstractVector +end +mutable struct Lindblad_noiseless_controlled_pure <: LindbladDynamicsData + H0::AbstractVecOrMat + dH::AbstractVector + ψ0::AbstractVector + tspan::AbstractVector + Hc::AbstractVector + ctrl::AbstractVector +end + +mutable struct Lindblad_noisy_controlled_pure <: LindbladDynamicsData + H0::AbstractVecOrMat + dH::AbstractVector + ψ0::AbstractVector + tspan::AbstractVector + decay_opt::AbstractVector + γ::AbstractVector + Hc::AbstractVector + ctrl::AbstractVector +end + +# Constructor of Lindblad dynamics +Lindblad( + H0::AbstractMatrix, + dH::AbstractVector, + ρ0::AbstractMatrix, + tspan::AbstractVector, +) = Lindblad(Lindblad_noiseless_free(H0, dH, ρ0, tspan), :noiseless, :free) + +Lindblad( + H0::AbstractMatrix, + dH::AbstractVector, + ρ0::AbstractMatrix, + tspan::AbstractVector, + decay_opt::AbstractVector, + γ::AbstractVector, +) = Lindblad(Lindblad_noisy_free(H0, dH, ρ0, tspan, decay_opt, γ), :noisy, :free) + +Lindblad( + H0::AbstractVector, + dH::AbstractVector, + ρ0::AbstractMatrix, + tspan::AbstractVector, +) = Lindblad(Lindblad_noiseless_timedepend(H0, dH, ρ0, tspan), :noiseless, :timedepend) + +Lindblad( + H0::AbstractVector, + dH::AbstractVector, + ρ0::AbstractMatrix, + tspan::AbstractVector, + decay_opt::AbstractVector, + γ::AbstractVector, +) = Lindblad( + Lindblad_noisy_timedepend(H0, dH, ρ0, tspan, decay_opt, γ), + :noisy, + :timedepend, +) + +Lindblad( + H0::AbstractVecOrMat, + dH::AbstractVector, + Hc::AbstractVector, + ctrl::AbstractVector, + ρ0::AbstractMatrix, + tspan::AbstractVector, +) = Lindblad( + Lindblad_noiseless_controlled(H0, dH, ρ0, tspan, Hc, ctrl), + :noiseless, + :controlled, +) + +Lindblad( + H0::AbstractVecOrMat, + dH::AbstractVector, + Hc::AbstractVector, + ctrl::AbstractVector, + ρ0::AbstractMatrix, + tspan::AbstractVector, + decay_opt::AbstractVector, + γ::AbstractVector, +) = Lindblad( + Lindblad_noisy_controlled(H0, dH, ρ0, tspan, decay_opt, γ, Hc, ctrl), + :noisy, + :controlled, +) + +Lindblad( + H0::AbstractMatrix, + dH::AbstractVector, + ψ0::AbstractVector, + tspan::AbstractVector, +) = Lindblad(Lindblad_noiseless_free_pure(H0, dH, ψ0, tspan), :noiseless, :free, :ket) + +Lindblad( + H0::AbstractMatrix, + dH::AbstractVector, + ψ0::AbstractVector, + tspan::AbstractVector, + decay_opt::AbstractVector, + γ::AbstractVector, +) = Lindblad(Lindblad_noisy_free_pure(H0, dH, ψ0, tspan, decay_opt, γ), :noisy, :free, :ket) + +Lindblad( + H0::AbstractVector, + dH::AbstractVector, + ψ0::AbstractVector, + tspan::AbstractVector, +) = Lindblad( + Lindblad_noiseless_timedepend(H0, dH, ψ0, tspan), + :noiseless, + :timedepend, + :ket, +) + +Lindblad( + H0::AbstractVector, + dH::AbstractVector, + ψ0::AbstractVector, + tspan::AbstractVector, + decay_opt::AbstractVector, + γ::AbstractVector, +) = Lindblad( + Lindblad_noisy_timedepend_pure(H0, dH, ψ0, tspan, decay_opt, γ), + :noisy, + :timedepend, + :ket, +) + +Lindblad( + H0::AbstractVecOrMat, + dH::AbstractVector, + Hc::AbstractVector, + ctrl::AbstractVector, + ψ0::AbstractVector, + tspan::AbstractVector, +) = Lindblad( + Lindblad_noiseless_controlled(H0, dH, ψ0, tspan, Hc, ctrl), + :noiseless, + :controlled, + :ket, +) + +Lindblad( + H0::AbstractVecOrMat, + dH::AbstractVector, + Hc::AbstractVector, + ctrl::AbstractVector, + ψ0::AbstractVector, + tspan::AbstractVector, + decay_opt::AbstractVector, + γ::AbstractVector, +) = Lindblad( + Lindblad_noisy_controlled_pure(H0, dH, ψ0, tspan, decay_opt, γ, Hc, ctrl), + :noisy, + :controlled, + :ket, +) + +Lindblad(data::LindbladDynamicsData, N, C, R) = + para_type(data) |> P -> Lindblad{((N, C, R, P) .|> eval)...}(data, N, C, R, P) +Lindblad(data::LindbladDynamicsData, N, C) = Lindblad(data, N, C, :dm) + +para_type(data::LindbladDynamicsData) = length(data.dH) == 1 ? :single_para : :multi_para + +get_dim(d::Lindblad_noiseless_free) = size(d.ρ0, 1) +get_dim(d::Lindblad_noisy_free) = size(d.ρ0, 1) +get_dim(d::Lindblad_noiseless_controlled) = size(d.ρ0, 1) +get_dim(d::Lindblad_noisy_controlled) = size(d.ρ0, 1) +get_dim(d::Lindblad_noiseless_timedepend) = size(d.ρ0, 1) +get_dim(d::Lindblad_noisy_timedepend) = size(d.ρ0, 1) +get_dim(d::Lindblad_noiseless_free_pure) = size(d.ψ0, 1) +get_dim(d::Lindblad_noisy_free_pure) = size(d.ψ0, 1) +get_dim(d::Lindblad_noiseless_controlled_pure) = size(d.ψ0, 1) +get_dim(d::Lindblad_noisy_controlled_pure) = size(d.ψ0, 1) +get_dim(d::Lindblad_noiseless_timedepend_pure) = size(d.ψ0, 1) +get_dim(d::Lindblad_noisy_timedepend_pure) = size(d.ψ0, 1) + +get_para(d::LindbladDynamicsData) = length(d.dH) \ No newline at end of file diff --git a/src/dynamics/Lindblad.jl b/src/dynamics/lindblad/LindbladDynamics.jl similarity index 71% rename from src/dynamics/Lindblad.jl rename to src/dynamics/lindblad/LindbladDynamics.jl index c4a024c..e3b058b 100644 --- a/src/dynamics/Lindblad.jl +++ b/src/dynamics/lindblad/LindbladDynamics.jl @@ -1,281 +1,3 @@ - -# dynamics in Lindblad form -struct Lindblad{N,C,R,P} <: AbstractDynamics - data::AbstractDynamicsData - noise_type::Symbol - ctrl_type::Symbol - state_rep::Symbol - para_type::Symbol -end - -Lindblad(data, N, C, R) = - para_type(data) |> P -> Lindblad{((N, C, R, P) .|> eval)...}(data, N, C, R, P) -Lindblad(data, N, C) = Lindblad(data, N, C, :dm) - -mutable struct Lindblad_noiseless_free <: AbstractDynamicsData - H0::AbstractMatrix - dH::AbstractVector - ρ0::AbstractMatrix - tspan::AbstractVector -end - -mutable struct Lindblad_noisy_free <: AbstractDynamicsData - H0::AbstractMatrix - dH::AbstractVector - ρ0::AbstractMatrix - tspan::AbstractVector - decay_opt::AbstractVector - γ::AbstractVector -end - - -mutable struct Lindblad_noiseless_timedepend <: AbstractDynamicsData - H0::AbstractVector - dH::AbstractVector - ρ0::AbstractMatrix - tspan::AbstractVector -end - -mutable struct Lindblad_noisy_timedepend <: AbstractDynamicsData - H0::AbstractVector - dH::AbstractVector - ρ0::AbstractMatrix - tspan::AbstractVector - decay_opt::AbstractVector - γ::AbstractVector -end - -mutable struct Lindblad_noiseless_controlled <: AbstractDynamicsData - H0::AbstractVecOrMat - dH::AbstractVector - ρ0::AbstractMatrix - tspan::AbstractVector - Hc::AbstractVector - ctrl::AbstractVector -end - -mutable struct Lindblad_noisy_controlled <: AbstractDynamicsData - H0::AbstractVecOrMat - dH::AbstractVector - ρ0::AbstractMatrix - tspan::AbstractVector - decay_opt::AbstractVector - γ::AbstractVector - Hc::AbstractVector - ctrl::AbstractVector -end - - -mutable struct Lindblad_noiseless_free_pure <: AbstractDynamicsData - H0::AbstractMatrix - dH::AbstractVector - ψ0::AbstractVector - tspan::AbstractVector -end - -mutable struct Lindblad_noisy_free_pure <: AbstractDynamicsData - H0::AbstractMatrix - dH::AbstractVector - ψ0::AbstractVector - tspan::AbstractVector - decay_opt::AbstractVector - γ::AbstractVector -end - -mutable struct Lindblad_noiseless_timedepend_pure <: AbstractDynamicsData - H0::AbstractVector - dH::AbstractVector - ψ0::AbstractVector - tspan::AbstractVector -end - -mutable struct Lindblad_noisy_timedepend_pure <: AbstractDynamicsData - H0::AbstractVector - dH::AbstractVector - ψ0::AbstractVector - tspan::AbstractVector - decay_opt::AbstractVector - γ::AbstractVector -end -mutable struct Lindblad_noiseless_controlled_pure <: AbstractDynamicsData - H0::AbstractVecOrMat - dH::AbstractVector - ψ0::AbstractVector - tspan::AbstractVector - Hc::AbstractVector - ctrl::AbstractVector -end - -mutable struct Lindblad_noisy_controlled_pure <: AbstractDynamicsData - H0::AbstractVecOrMat - dH::AbstractVector - ψ0::AbstractVector - tspan::AbstractVector - decay_opt::AbstractVector - γ::AbstractVector - Hc::AbstractVector - ctrl::AbstractVector -end - - -para_type(data::AbstractDynamicsData) = length(data.dH) == 1 ? :single_para : :multi_para - -# Constructor of Lindblad dynamics -Lindblad( - H0::AbstractMatrix, - dH::AbstractVector, - ρ0::AbstractMatrix, - tspan::AbstractVector, -) = Lindblad(Lindblad_noiseless_free(H0, dH, ρ0, tspan), :noiseless, :free) - -Lindblad( - H0::AbstractMatrix, - dH::AbstractVector, - ρ0::AbstractMatrix, - tspan::AbstractVector, - decay_opt::AbstractVector, - γ::AbstractVector, -) = Lindblad(Lindblad_noisy_free(H0, dH, ρ0, tspan, decay_opt, γ), :noisy, :free) - -Lindblad( - H0::AbstractVector, - dH::AbstractVector, - ρ0::AbstractMatrix, - tspan::AbstractVector, -) = Lindblad(Lindblad_noiseless_timedepend(H0, dH, ρ0, tspan), :noiseless, :timedepend) - -Lindblad( - H0::AbstractVector, - dH::AbstractVector, - ρ0::AbstractMatrix, - tspan::AbstractVector, - decay_opt::AbstractVector, - γ::AbstractVector, -) = Lindblad( - Lindblad_noisy_timedepend(H0, dH, ρ0, tspan, decay_opt, γ), - :noisy, - :timedepend, -) - -Lindblad( - H0::AbstractVecOrMat, - dH::AbstractVector, - Hc::AbstractVector, - ctrl::AbstractVector, - ρ0::AbstractMatrix, - tspan::AbstractVector, -) = Lindblad( - Lindblad_noiseless_controlled(H0, dH, ρ0, tspan, Hc, ctrl), - :noiseless, - :controlled, -) - -Lindblad( - H0::AbstractVecOrMat, - dH::AbstractVector, - Hc::AbstractVector, - ctrl::AbstractVector, - ρ0::AbstractMatrix, - tspan::AbstractVector, - decay_opt::AbstractVector, - γ::AbstractVector, -) = Lindblad( - Lindblad_noisy_controlled(H0, dH, ρ0, tspan, decay_opt, γ, Hc, ctrl), - :noisy, - :controlled, -) - -Lindblad( - H0::AbstractMatrix, - dH::AbstractVector, - ψ0::AbstractVector, - tspan::AbstractVector, -) = Lindblad(Lindblad_noiseless_free_pure(H0, dH, ψ0, tspan), :noiseless, :free, :ket) - -Lindblad( - H0::AbstractMatrix, - dH::AbstractVector, - ψ0::AbstractVector, - tspan::AbstractVector, - decay_opt::AbstractVector, - γ::AbstractVector, -) = Lindblad(Lindblad_noisy_free_pure(H0, dH, ψ0, tspan, decay_opt, γ), :noisy, :free, :ket) - -Lindblad( - H0::AbstractVector, - dH::AbstractVector, - ψ0::AbstractVector, - tspan::AbstractVector, -) = Lindblad( - Lindblad_noiseless_timedepend(H0, dH, ψ0, tspan), - :noiseless, - :timedepend, - :ket, -) - -Lindblad( - H0::AbstractVector, - dH::AbstractVector, - ψ0::AbstractVector, - tspan::AbstractVector, - decay_opt::AbstractVector, - γ::AbstractVector, -) = Lindblad( - Lindblad_noisy_timedepend_pure(H0, dH, ψ0, tspan, decay_opt, γ), - :noisy, - :timedepend, - :ket, -) - -Lindblad( - H0::AbstractVecOrMat, - dH::AbstractVector, - Hc::AbstractVector, - ctrl::AbstractVector, - ψ0::AbstractVector, - tspan::AbstractVector, -) = Lindblad( - Lindblad_noiseless_controlled(H0, dH, ψ0, tspan, Hc, ctrl), - :noiseless, - :controlled, - :ket, -) - -Lindblad( - H0::AbstractVecOrMat, - dH::AbstractVector, - Hc::AbstractVector, - ctrl::AbstractVector, - ψ0::AbstractVector, - tspan::AbstractVector, - decay_opt::AbstractVector, - γ::AbstractVector, -) = Lindblad( - Lindblad_noisy_controlled_pure(H0, dH, ψ0, tspan, decay_opt, γ, Hc, ctrl), - :noisy, - :controlled, - :ket, -) - -function set_ctrl(dynamics::Lindblad, ctrl) - temp = deepcopy(dynamics) - temp.data.ctrl = ctrl - temp -end - -function set_state(dynamics::Lindblad, state::AbstractVector) - temp = deepcopy(dynamics) - temp.data.ψ0 = state - temp -end - -function set_state(dynamics::Lindblad, state::AbstractMatrix) - temp = deepcopy(dynamics) - temp.data.ρ0 = state - temp -end - -# functions for evolve dynamics in Lindblad form function liouville_commu(H) kron(one(H), H) - kron(H |> transpose, one(H)) end @@ -498,7 +220,7 @@ function secondorder_derivative( 2 * im * Δt * dH_L[i] * ∂ρt_∂x[i] for i = 1:para_num ] + [exp_L] .* ∂2ρt_∂x end - ρt = exp(vec(H[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat, ∂2ρt_∂x |> vec2mat end @@ -534,7 +256,7 @@ function evolve(dynamics::Lindblad{noiseless,timedepend,ket}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H0[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H0[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end @@ -552,7 +274,7 @@ function evolve(dynamics::Lindblad{noiseless,free,dm}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end @@ -588,7 +310,7 @@ function evolve(dynamics::Lindblad{noiseless,timedepend,ket}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H0[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H0[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end @@ -606,7 +328,7 @@ function evolve(dynamics::Lindblad{noiseless,free,dm}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H0[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end @@ -624,7 +346,7 @@ function evolve(dynamics::Lindblad{noiseless,timedepend,dm}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H0[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H0[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end @@ -680,7 +402,7 @@ function evolve(dynamics::Lindblad{noisy,timedepend,ket}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H0[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H0[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end @@ -699,7 +421,7 @@ function evolve(dynamics::Lindblad{noisy,timedepend,dm}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H0[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H0[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end @@ -722,7 +444,7 @@ function evolve(dynamics::Lindblad{noiseless,controlled,dm}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end @@ -744,18 +466,18 @@ function evolve(dynamics::Lindblad{noisy,controlled,dm}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end #### evolution of state under time-independent Hamiltonian with noise and controls #### function evolve(dynamics::Lindblad{noisy,controlled,ket}) - (; H0, dH, ψ0, tspan, decay_opt, γ, Hc, ctrl) = dynamics.data + (; H0, dH, ψ0, tspan, decay_opt, γ, ctrl, Hc) = dynamics.data para_num = length(dH) ctrl_num = length(Hc) ctrl_interval = ((length(tspan) - 1) / length(ctrl[1])) |> Int - ctrl = [repeat(ctrl[i], 1, ctrl_interval) |> transpose |> vec for i = 1:ctrl_num] + ctrl = [repeat(dynamics.data.ctrl[i], 1, ctrl_interval) |> transpose |> vec for i = 1:ctrl_num] H = Htot(H0, Hc, ctrl) dH_L = [liouville_commu(dH[i]) for i = 1:para_num] ρt = (ψ0 * ψ0') |> vec @@ -766,7 +488,7 @@ function evolve(dynamics::Lindblad{noisy,controlled,ket}) ρt = exp_L * ρt ∂ρt_∂x = [-im * Δt * dH_L[i] * ρt for i = 1:para_num] + [exp_L] .* ∂ρt_∂x end - ρt = exp(vec(H[end])' * zero(ρt)) * ρt + # ρt = exp(vec(H[end])' * zero(ρt)) * ρt ρt |> vec2mat, ∂ρt_∂x |> vec2mat end diff --git a/src/dynamics/lindblad/LindbladWrapper.jl b/src/dynamics/lindblad/LindbladWrapper.jl new file mode 100644 index 0000000..87bc0d9 --- /dev/null +++ b/src/dynamics/lindblad/LindbladWrapper.jl @@ -0,0 +1,512 @@ +# wrapper for Lindblad dynamics with ControlOpt +function Lindblad(opt::ControlOpt, tspan, ρ₀, H0, dH, Hc=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) + (;ctrl) = opt + dim = size(ρ₀, 1) + if ismissing(dH) + dH = [zeros(ComplexF64, dim, dim)] + end + ctrl_num = length(Hc) + tnum = length(tspan) + + if ismissing(decay) + decay_opt = [zeros(ComplexF64, dim, dim)] + γ = [0.0] + else + decay_opt = [decay[1] for decay in decay] + γ = [decay[2] for decay in decay] + end + + if ismissing(Hc) + Hc = [zeros(ComplexF64, dim, dim)] + ctrl = [zeros(tnum-1)] + elseif ismissing(ctrl) + opt.ctrl = [zeros(tnum-1) for _ in 1:ctrl_num] + opt.ctrl = ctrl + else + ctrl_length = length(ctrl) + if ctrl_num < ctrl_length + throw(ArgumentError( + "There are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences: too many coefficients sequences" + )) + elseif ctrl_num < ctrl_length + throw(ArgumentError( + "Not enough coefficients sequences: there are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences. The rest of the control sequences are set to be 0." + )) + end + + ratio_num = ceil((length(tspan)-1) / length(ctrl[1])) + if length(tspan) - 1 % length(ctrl[1]) != 0 + tnum = ratio_num * length(ctrl[1]) |> Int + tspan = range(tspan[1], tspan[end], length=tnum+1) + end + end + H0 = complex(H0) + dH = complex.(dH) + ρ₀ = complex(ρ₀) + + Lindblad(H0, dH, Hc, ctrl, ρ₀, tspan, decay_opt, γ) +end + +function Lindblad(opt::StateOpt, tspan, H0, dH, Hc=missing, ctrl=missing, decay=missing;rng=GLOBAL_RNG) + (;ψ₀) = opt + dim = H0 isa AbstractVector ? size(H0[1], 1) : size(H0, 1) + if ismissing(ψ₀) + r_ini = 2*rand(rng, dim) - ones(dim) + r = r_ini ./ norm(r_ini) + ϕ = 2pi*rand(rng, dim) + ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.ψ₀ = ψ₀ + end + + if ismissing(dH) + dH = [zeros(ComplexF64, dim, dim)] + end + + if !ismissing(Hc) && !ismissing(ctrl) + ctrl_length = length(ctrl) + ctrl_num = length(Hc) + if ctrl_num < ctrl_length + throw(ArgumentError( + "Too many contrl coefficients: there are $ctrl_num control Hamiltonians but $ctrl_length control coefficients given." + )) + elseif ctrl_num < ctrl_length + throw(ArgumentError( + "Insufficient coefficients sequences: there are $ctrl_num control Hamiltonians but $ctrl_length coefficients given. The rest of the control sequences are setten to be 0." + )) + end + + if length(ctrl[1]) == 1 + hc = sum([c[1]*hc for (c,hc) in zip(ctrl,Hc)]) + + if typeof(H0) <: AbstractMatrix + H0 = complex(H0 + hc) + elseif typeof(H0) <: AbstractVector + H0 = [complex(h0+hc) for h0 in H0] + else + ## TODO wrong type of H0 + end + else + ratio_num = ceil((length(tspan)-1) / length(ctrl[1])) + + if length(tspan) - 1 % length(ctrl[1]) != 0 + tnum = ratio_num * length(ctrl[1]) |> Int + tspan = range(tspan[1], tspan[end], length=tnum+1) + if typeof(H0) <: AbstractVector + itp = interpolate((tspan,), H0, Gridded(Linear())) + H0 = itp(tspan) + end + end + + hc = [sum([c[i]*hc for (c,hc) in zip(ctrl,Hc)]) for i in 1:length(ctrl[1]) ] + + if typeof(H0) <: AbstractMatrix + H0 = [complex(H0+hc) for hc in hc] + elseif typeof(H0) <: AbstractVector + H0 = complex.(H0+hc) + else + ## TODO wrong type of H0 + end + end + end + if ismissing(decay) + γ = [0.0] + else + decay_opt = [decay[1] for decay in decay] + γ = [decay[2] for decay in decay] + end + + dH = complex.(dH) + ψ₀ = complex(ψ₀) + + if all(iszero.(γ)) # if any non-zero decay rate + return Lindblad(H0, dH, ψ₀, tspan) + else + return Lindblad(H0, dH, ψ₀, tspan, decay_opt, γ) + end +end + +function _ini_measurement!(opt::Mopt_Projection, dim::Int, rng; eps=GLOBAL_EPS) + (;C) = opt + ## initialize the Mopt target C + if ismissing(C) + M = [ComplexF64[] for _ in 1:dim] + for i in 1:dim + 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 + C = gramschmidt(M) + end + + opt.C = complex.(C) +end + +function _ini_measurement!(opt::Mopt_LinearComb, dim::Int, rng; eps=GLOBAL_EPS) + (;B, POVM_basis, M_num) = opt + if ismissing(POVM_basis) + opt.POVM_basis = SIC(dim) + else + ## TODO: accuracy ? + for P in POVM_basis + if minimum(eigvals(P)) < (-eps) + throw(ArgumentError("The given POVMs should be semidefinite!")) + end + end + if !(sum(POVM_basis) ≈ I(dim)) + throw(ArgumentError("The sum of the given POVMs should be identity matrix!")) + end + end + + if ismissing(B) + opt.B = [rand(rng, length(POVM_basis)) for _ in 1:M_num] + end +end + + +function _ini_measurement!(opt::Mopt_Rotation, dim::Int, rng; eps=GLOBAL_EPS) + (;s, POVM_basis, Lambda) = opt + if ismissing(POVM_basis) + throw(ArgumentError("The initial POVM basis should not be empty!")) + else + ## TODO: accuracy ? + for P in POVM_basis + if minimum(eigvals(P)) < -eps + throw(ArgumentError("The given POVMs should be semidefinite!")) + end + end + if !(sum(POVM_basis) ≈ I(dim)) + throw(ArgumentError("The sum of the given POVMs should be identity matrix!")) + end + end + + if ismissing(s) + opt.s = rand(rng, dim^2) + end +end + +function Lindblad(opt::AbstractMopt, tspan, ρ₀, H0, dH, Hc=missing, ctrl=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) + dim = size(ρ₀, 1) + _ini_measurement!(opt, dim, rng; eps=eps) + + if ismissing(dH) + dH = [zeros(ComplexF64, dim, dim)] + end + + if !ismissing(Hc) && !ismissing(ctrl) + ctrl_length = length(ctrl) + ctrl_num = length(Hc) + if ctrl_num < ctrl_length + throw(ArgumentError( + "Too many contrl coefficients: there are $ctrl_num control Hamiltonians but $ctrl_length control coefficients given." + )) + elseif ctrl_num < ctrl_length + throw(ArgumentError( + "Insufficient coefficients sequences: there are $ctrl_num control Hamiltonians but $ctrl_length coefficients given. The rest of the control sequences are setten to be 0." + )) + end + + if length(ctrl[1]) == 1 + hc = sum([c[1]*hc for (c,hc) in zip(ctrl,Hc)]) + + if typeof(H0) <: AbstractMatrix + H0 = complex(H0+hc) + elseif typeof(H0) <: AbstractVector + H0 = [complex(h0+hc) for h0 in H0] + else + ## TODO wrong type of H0 + end + else + ratio_num = ceil((length(tspan)-1) / length(ctrl[1])) + + if length(tspan) - 1 % length(ctrl[1]) != 0 + tnum = ratio_num * length(ctrl[1]) |> Int + tspan = range(tspan[1], tspan[end], length=tnum+1) + if typeof(H0) <: AbstractVector + itp = interpolate((tspan,), H0, Gridded(Linear())) + H0 = itp(tspan) + end + end + + hc = [sum([c[i]*hc for (c,hc) in zip(ctrl,Hc)]) for i in 1:length(ctrl[1]) ] + + if typeof(H0) <: AbstractMatrix + H0 = [complex(H0+hc) for hc in hc] + elseif typeof(H0) <: AbstractVector + H0 = complex.(H0+hc) + else + ## TODO wrong type of H0 + end + end + end + if ismissing(decay) + γ = [0.0] + else + decay_opt = [decay[1] for decay in decay] + γ = [decay[2] for decay in decay] + end + + dH = complex.(dH) + ρ₀ = complex(ρ₀) + + if all(iszero.(γ)) # if any non-zero decay rate + return Lindblad(H0, dH, ρ₀, tspan) + else + return Lindblad(H0, dH, ρ₀, tspan, decay_opt, γ) + end +end + +function _ini_measurement!(opt::CompOpt, dim::Int, rng; eps=GLOBAL_EPS) + (;C) = opt + ## initialize the Mopt target C + M = [ComplexF64[] for _ in 1:dim] + if ismissing(C) + for i in 1:dim + 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 + opt.C = gramschmidt(M) + end +end + +function Lindblad(opt::StateControlOpt, tspan, H0, dH, Hc=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) + (;ψ₀, ctrl) = opt + dim = H0 isa AbstractVector ? size(H0[1], 1) : size(H0, 1) + if ismissing(ψ₀) + r_ini = 2*rand(rng, dim) - ones(dim) + r = r_ini ./ norm(r_ini) + ϕ = 2pi*rand(rng, dim) + ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.ψ₀ = ψ₀ + end + if ismissing(dH) + dH = [zeros(ComplexF64, dim, dim)] + end + ctrl_num = length(Hc) + tnum = length(tspan) + + if ismissing(decay) + decay_opt = [zeros(ComplexF64, dim, dim)] + γ = [0.0] + else + decay_opt = [decay[1] for decay in decay] + γ = [decay[2] for decay in decay] + end + + if ismissing(Hc) + Hc = [zeros(ComplexF64, dim, dim)] + opt.ctrl = [zeros(tnum-1)] + elseif ismissing(ctrl) + ctrl = [zeros(tnum-1) for _ in 1:ctrl_num] + opt.ctrl = ctrl + else + ctrl_length = length(ctrl) + if ctrl_num < ctrl_length + throw(ArgumentError( + "There are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences: too many coefficients sequences" + )) + elseif ctrl_num < ctrl_length + throw(ArgumentError( + "Not enough coefficients sequences: there are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences. The rest of the control sequences are set to be 0." + )) + end + + ratio_num = ceil((length(tspan)-1) / length(ctrl[1])) + if length(tspan) - 1 % length(ctrl[1]) != 0 + tnum = ratio_num * length(ctrl[1]) |> Int + tspan = range(tspan[1], tspan[end], length=tnum+1) + end + end + H0 = complex(H0) + dH = complex.(dH) + ψ₀ = complex(ψ₀) + Lindblad(H0, dH, Hc, ctrl, ψ₀, tspan, decay_opt, γ) +end + +function Lindblad(opt::ControlMeasurementOpt, tspan, ρ₀, H0, dH, Hc=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) + (;ctrl) = opt + dim = size(ρ₀, 1) + _ini_measurement!(opt, dim, rng; eps=eps) + if ismissing(dH) + dH = [zeros(ComplexF64, dim, dim)] + end + ctrl_num = length(Hc) + tnum = length(tspan) + + if ismissing(decay) + decay_opt = [zeros(ComplexF64, dim, dim)] + γ = [0.0] + else + decay_opt = [decay[1] for decay in decay] + γ = [decay[2] for decay in decay] + end + + if ismissing(Hc) + Hc = [zeros(ComplexF64, dim, dim)] + opt.ctrl = [zeros(tnum-1)] + elseif ismissing(ctrl) + ctrl = [zeros(tnum-1) for _ in 1:ctrl_num] + opt.ctrl = ctrl + else + ctrl_length = length(ctrl) + if ctrl_num < ctrl_length + throw(ArgumentError( + "There are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences: too many coefficients sequences" + )) + elseif ctrl_num < ctrl_length + throw(ArgumentError( + "Not enough coefficients sequences: there are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences. The rest of the control sequences are set to be 0." + )) + end + + ratio_num = ceil((length(tspan)-1) / length(ctrl[1])) + if length(tspan) - 1 % length(ctrl[1]) != 0 + tnum = ratio_num * length(ctrl[1]) |> Int + tspan = range(tspan[1], tspan[end], length=tnum+1) + end + end + H0 = complex(H0) + dH = complex.(dH) + ρ₀ = complex(ρ₀) + + Lindblad(H0, dH, Hc, ctrl, ρ₀, tspan, decay_opt, γ) +end + +function Lindblad(opt::StateMeasurementOpt, tspan, H0, dH, Hc=missing, ctrl=missing, decay=missing;rng=GLOBAL_RNG) + (;ψ₀) = opt + dim = H0 isa AbstractVector ? size(H0[1], 1) : size(H0, 1) + _ini_measurement!(opt, dim, rng; eps=eps) + if ismissing(ψ₀) + r_ini = 2*rand(rng, dim) - ones(dim) + r = r_ini ./ norm(r_ini) + ϕ = 2pi*rand(rng, dim) + ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.ψ₀ = ψ₀ + end + + if ismissing(dH) + dH = [zeros(ComplexF64, dim, dim)] + end + + if !ismissing(Hc) && !ismissing(ctrl) + ctrl_length = length(ctrl) + ctrl_num = length(Hc) + if ctrl_num < ctrl_length + throw(ArgumentError( + "Too many contrl coefficients: there are $ctrl_num control Hamiltonians but $ctrl_length control coefficients given." + )) + elseif ctrl_num < ctrl_length + throw(ArgumentError( + "Insufficient coefficients sequences: there are $ctrl_num control Hamiltonians but $ctrl_length coefficients given. The rest of the control sequences are setten to be 0." + )) + end + + if length(ctrl[1]) == 1 + hc = sum([c[1]*hc for (c,hc) in zip(ctrl,Hc)]) + + if typeof(H0) <: AbstractMatrix + H0 = complex(H0 + hc) + elseif typeof(H0) <: AbstractVector + H0 = [complex(h0+hc) for h0 in H0] + else + ## TODO wrong type of H0 + end + else + ratio_num = ceil((length(tspan)-1) / length(ctrl[1])) + + if length(tspan) - 1 % length(ctrl[1]) != 0 + tnum = ratio_num * length(ctrl[1]) |> Int + tspan = range(tspan[1], tspan[end], length=tnum+1) + if typeof(H0) <: AbstractVector + itp = interpolate((tspan,), H0, Gridded(Linear())) + H0 = itp(tspan) + end + end + + hc = [sum([c[i]*hc for (c,hc) in zip(ctrl,Hc)]) for i in 1:length(ctrl[1]) ] + + if typeof(H0) <: AbstractMatrix + H0 = [complex(H0+hc) for hc in hc] + elseif typeof(H0) <: AbstractVector + H0 = complex.(H0+hc) + else + ## TODO wrong type of H0 + end + end + end + if ismissing(decay) + γ = [0.0] + else + decay_opt = [decay[1] for decay in decay] + γ = [decay[2] for decay in decay] + end + + dH = complex.(dH) + ψ₀ = complex(ψ₀) + + if all(iszero.(γ)) # if any non-zero decay rate + return Lindblad(H0, dH, ψ₀, tspan) + else + return Lindblad(H0, dH, ψ₀, tspan, decay_opt, γ) + end +end + +function Lindblad(opt::StateControlMeasurementOpt, tspan, H0, dH, Hc=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) + (;ctrl, ψ₀) = opt + dim = H0 isa AbstractVector ? size(H0[1], 1) : size(H0, 1) + _ini_measurement!(opt, dim, rng; eps=eps) + + if ismissing(ψ₀) + r_ini = 2*rand(rng, dim) - ones(dim) + r = r_ini ./ norm(r_ini) + ϕ = 2pi*rand(rng, dim) + ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.ψ₀ = ψ₀ + end + + if ismissing(dH) + dH = [zeros(ComplexF64, dim, dim)] + end + ctrl_num = length(Hc) + tnum = length(tspan) + + if ismissing(decay) + decay_opt = [zeros(ComplexF64, dim, dim)] + γ = [0.0] + else + decay_opt = [decay[1] for decay in decay] + γ = [decay[2] for decay in decay] + end + + if ismissing(Hc) + Hc = [zeros(ComplexF64, dim, dim)] + opt.ctrl = [zeros(tnum-1)] + elseif ismissing(ctrl) + ctrl = [zeros(tnum-1) for _ in 1:ctrl_num] + opt.ctrl = ctrl + else + ctrl_length = length(ctrl) + if ctrl_num < ctrl_length + throw(ArgumentError( + "There are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences: too many coefficients sequences" + )) + elseif ctrl_num < ctrl_length + throw(ArgumentError( + "Not enough coefficients sequences: there are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences. The rest of the control sequences are set to be 0." + )) + end + + ratio_num = ceil((length(tspan)-1) / length(ctrl[1])) + if length(tspan) - 1 % length(ctrl[1]) != 0 + tnum = ratio_num * length(ctrl[1]) |> Int + tspan = range(tspan[1], tspan[end], length=tnum+1) + end + end + H0 = complex(H0) + dH = complex.(dH) + ψ₀ = complex(ψ₀) + + Lindblad(H0, dH, Hc, ctrl, ψ₀, tspan, decay_opt, γ) +end \ No newline at end of file diff --git a/src/objective/AsymptoticBound/AsymptoticBound.jl b/src/objective/AsymptoticBound/AsymptoticBound.jl index f2f1978..040aed4 100644 --- a/src/objective/AsymptoticBound/AsymptoticBound.jl +++ b/src/objective/AsymptoticBound/AsymptoticBound.jl @@ -4,29 +4,29 @@ abstract type SLD <: AbstractLDtype end abstract type RLD <: AbstractLDtype end abstract type LLD <: AbstractLDtype end -struct QFIM_Obj{P,D} <: AbstractObj - W::AbstractMatrix - eps::Number +Base.@kwdef struct QFIM_Obj{P,D} <: AbstractObj + W::Union{AbstractMatrix, Missing} + eps::Number = GLOBAL_EPS end -struct CFIM_Obj{P} <: AbstractObj - M::AbstractVecOrMat - W::AbstractMatrix - eps::Number +Base.@kwdef struct CFIM_Obj{P} <: AbstractObj + M::Union{AbstractVecOrMat, Missing} + W::Union{AbstractMatrix, Missing} + eps::Number = GLOBAL_EPS end -struct HCRB_Obj{P} <: AbstractObj - W::AbstractMatrix - eps::Number +Base.@kwdef struct HCRB_Obj{P} <: AbstractObj + W::Union{AbstractMatrix, Missing} + eps::Number = GLOBAL_EPS end -QFIM_Obj(W, eps, syms::Symbol...) = QFIM_Obj{eval.(syms)...}(W, eps) -CFIM_Obj(M, W, eps, syms::Symbol...) = CFIM_Obj{eval.(syms)...}(M, W, eps) -HCRB_Obj(W, eps, syms::Symbol...) = HCRB_Obj{eval.(syms)...}(W, eps) +QFIM_Obj(;W=missing, eps=GLOBAL_EPS, para_type::Symbol=:single_para, LD_type::Symbol=:SLD) = QFIM_Obj{eval.([para_type, LD_type])...}(W, eps) +CFIM_Obj(;M=missing, W=missing, eps=GLOBAL_EPS, para_type::Symbol=:single_para) = CFIM_Obj{eval(para_type)}(M, W, eps) +HCRB_Obj(;W=missing, eps=GLOBAL_EPS, para_type::Symbol=:single_para) = HCRB_Obj{eval(para_type)}(W, eps) -QFIM_Obj(W, eps, str::String...) = QFIM_Obj(W, eps, Symbol.(str)...) -CFIM_Obj(M, W, eps, str::String...) = CFIM_Obj(M, W, eps, Symbol.(str)...) -HCRB_Obj(W, eps, str::String...) = HCRB_Obj(W, eps, Symbol.(str)...) +QFIM_Obj(W::AbstractMatrix, eps::Number, para_type::String, LD_type::String) = QFIM_Obj(W, eps, Symbol.([para_type, LD_type])...) +CFIM_Obj(M::AbstractVecOrMat, W::AbstractMatrix, eps::Number, para_type::String) = CFIM_Obj(M, W, eps, Symbol(para_type)) +HCRB_Obj(W::AbstractMatrix, eps::Number, para_type::String) = HCRB_Obj(W, eps, Symbol(para_type)) obj_type(::QFIM_Obj) = :QFIM obj_type(::CFIM_Obj) = :CFIM @@ -36,25 +36,24 @@ para_type(::QFIM_Obj{single_para,D}) where {D} = :single_para para_type(::QFIM_Obj{multi_para,D}) where {D} = :multi_para para_type(::CFIM_Obj{single_para}) = :single_para para_type(::CFIM_Obj{multi_para}) = :multi_para +para_type(::HCRB_Obj{single_para}) = :single_para para_type(::HCRB_Obj{multi_para}) = :multi_para LD_type(::QFIM_Obj{P,SLD}) where {P} = :SLD LD_type(::QFIM_Obj{P,RLD}) where {P} = :RLD LD_type(::QFIM_Obj{P,LLD}) where {P} = :LLD -QFIM_Obj(opt::CFIM{P}) where P = QFIM_Obj{P, SLD}(opt.W, opt.eps) -QFIM_Obj(opt::CFIM{P}, LDtype::Symbol) where P = QFIM_Obj{P, eval(LDtype)}(opt.W, opt.eps) +QFIM_Obj(opt::CFIM_Obj{P}) where P = QFIM_Obj{P, SLD}(opt.W, opt.eps) +QFIM_Obj(opt::CFIM_Obj{P}, LDtype::Symbol) where P = QFIM_Obj{P, eval(LDtype)}(opt.W, opt.eps) -function set_M(obj::CFIM_Obj{single_para}, M::AbstractMatrix) - temp = deepcopy(obj) - temp.M = M - temp -end +const obj_idx = Dict( + :QFIM => QFIM_Obj, + :CFIM => CFIM_Obj, + :HCRB => HCRB_Obj +) -function set_M(obj::CFIM_Obj{multi_para}, M::AbstractMatrix) - temp = deepcopy(obj) - temp.M = M - temp +function set_M(obj::CFIM_Obj{P}, M::AbstractVector) where P + CFIM_Obj{P}(M, obj.W, obj.eps) end function objective(obj::QFIM_Obj{single_para,SLD}, dynamics::Lindblad) @@ -64,6 +63,42 @@ function objective(obj::QFIM_Obj{single_para,SLD}, dynamics::Lindblad) return f, f end +function objective(obj::QFIM_Obj{single_para,SLD}, ρ, dρ) + (; W, eps) = obj + f = W[1] * QFIM_SLD(ρ, dρ[1]; eps = eps) + return f, f +end + +function objective(obj::QFIM_Obj{multi_para,SLD}, ρ, dρ) + (; W, eps) = obj + f = tr(W * pinv(QFIM_SLD(ρ, dρ; eps = eps))) + return f, 1.0 / f +end + +function objective(obj::QFIM_Obj{single_para,RLD}, ρ, dρ) + (; W, eps) = obj + f = W[1] * QFIM_RLD(ρ, dρ[1]; eps = eps) + return f, f +end + +function objective(obj::QFIM_Obj{multi_para,RLD}, ρ, dρ) + (; W, eps) = obj + f = tr(W * pinv(QFIM_RLD(ρ, dρ; eps = eps))) |> real + return f, 1.0 / f +end + +function objective(obj::QFIM_Obj{single_para,LLD}, ρ, dρ) + (; W, eps) = obj + f = W[1] * QFIM_LLD(ρ, dρ[1]; eps = eps) + return f, f +end + +function objective(obj::QFIM_Obj{multi_para,LLD}, ρ, dρ) + (; W, eps) = obj + f = tr(W * pinv(QFIM_LLD(ρ, dρ; eps = eps))) |> real + return f, 1.0 / f +end + function objective(obj::QFIM_Obj{multi_para,SLD}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) @@ -81,7 +116,7 @@ end function objective(obj::QFIM_Obj{multi_para,RLD}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) - f = tr(W * pinv(QFIM_RLD(ρ, dρ; eps = eps))) + f = tr(W * pinv(QFIM_RLD(ρ, dρ; eps = eps))) |> real return f, 1.0 / f end @@ -95,7 +130,7 @@ end function objective(obj::QFIM_Obj{multi_para,LLD}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) - f = tr(W * pinv(QFIM_LLD(ρ, dρ; eps = eps))) + f = tr(W * pinv(QFIM_LLD(ρ, dρ; eps = eps))) |> real return f, 1.0 / f end @@ -123,7 +158,7 @@ end function objective(obj::QFIM_Obj{multi_para,RLD}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) - f = tr(W * pinv(QFIM_RLD(ρ, dρ; eps = eps))) + f = tr(W * pinv(QFIM_RLD(ρ, dρ; eps = eps))) |> real return f, 1.0 / f end @@ -137,7 +172,19 @@ end function objective(obj::QFIM_Obj{multi_para,LLD}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) - f = tr(W * pinv(QFIM_LLD(ρ, dρ; eps = eps))) + f = tr(W * pinv(QFIM_LLD(ρ, dρ; eps = eps))) |> real + return f, 1.0 / f +end + +function objective(obj::CFIM_Obj{single_para}, ρ, dρ) + (; M, W, eps) = obj + f = W[1] * CFIM(ρ, dρ[1], M; eps = eps) + return f, f +end + +function objective(obj::CFIM_Obj{multi_para}, ρ, dρ) + (; M, W, eps) = obj + f = tr(W * pinv(CFIM(ρ, dρ, M; eps = eps))) return f, 1.0 / f end @@ -169,17 +216,23 @@ function objective(obj::CFIM_Obj{multi_para}, dynamics::Kraus) return f, 1.0 / f end +function objective(obj::HCRB_Obj{multi_para}, ρ, dρ) + (; W, eps) = obj + f = Holevo_bound_obj(ρ, dρ, W; eps = eps) + return f, 1.0 / f +end + function objective(obj::HCRB_Obj{multi_para}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) - f = Holevo_bound(ρ, dρ, W; eps = eps) + f = Holevo_bound_obj(ρ, dρ, W; eps = eps) return f, 1.0 / f end function objective(obj::HCRB_Obj{multi_para}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) - f = Holevo_bound(ρ, dρ, W; eps = eps) + f = Holevo_bound_obj(ρ, dρ, W; eps = eps) return f, 1.0 / f end @@ -254,16 +307,17 @@ function objective(opt::Mopt_Rotation, obj::CFIM_Obj{multi_para}, dynamics::Krau end ##### -function objective(::Type{Val{:expm}}, obj, dynamics) - temp = Float64[] - (; tspan, ctrl) = dynamics.data - for i = 1:length(tspan)-1 - dynamics_copy = set_ctrl(dynamics, [ctrl[1:i] for ctrl in ctrl]) - dynamics_copy.data.tspan = tspan[1:i+1] - append!(temp, objective(obj, dynamics_copy)) - end - temp -end # function objective +# function objective(::Type{Val{:expm}}, obj, dynamics) +# temp = [] +# (; tspan, ctrl) = dynamics.data +# for i = 1:length( ctrl) +# dynamics_copy = set_ctrl(dynamics, [ctrl[1:i] for ctrl in ctrl]) +# dynamics_copy.data.tspan = tspan[1:i+1] +# append!(temp, [objective(obj, dynamics_copy)]) +# end +# temp +# end # function objective include("CramerRao.jl") include("Holevo.jl") +include("AsymptoticBoundWrapper.jl") diff --git a/src/objective/AsymptoticBound/AsymptoticBoundWrapper.jl b/src/objective/AsymptoticBound/AsymptoticBoundWrapper.jl new file mode 100644 index 0000000..c16535e --- /dev/null +++ b/src/objective/AsymptoticBound/AsymptoticBoundWrapper.jl @@ -0,0 +1,39 @@ +function Objective(dynamics::AbstractDynamics, obj::QFIM_Obj{P,D}) where {P,D} + (;W, eps) = obj + if ismissing(W) + W = I(get_para(dynamics.data))|>Matrix + end + + d = LD_type(obj) |> eval + p = para_type(dynamics.data) |> eval + + return QFIM_Obj{p,d}(W,eps) +end + +function Objective(dynamics::AbstractDynamics, obj::CFIM_Obj{P}) where {P} + (;W, M, eps) = obj + if ismissing(W) + W = I(get_para(dynamics.data))|>Matrix + end + + if ismissing(M) + M = SIC(get_dim(dynamics.data)) + end + + p = para_type(dynamics.data) |> eval + + return CFIM_Obj{p}(M,W,eps) +end + +function Objective(dynamics::AbstractDynamics, obj::HCRB_Obj{P}) where {P} + (;W, eps) = obj + if ismissing(W) + W = I(get_para(dynamics.data))|>Matrix + end + + p = para_type(dynamics.data) |> eval + + return HCRB_Obj{p}(W,eps) +end + + diff --git a/src/objective/AsymptoticBound/Holevo.jl b/src/objective/AsymptoticBound/Holevo.jl index 63944d0..5dfa8d0 100644 --- a/src/objective/AsymptoticBound/Holevo.jl +++ b/src/objective/AsymptoticBound/Holevo.jl @@ -7,14 +7,14 @@ end function HCRB( ρ::Matrix{T}, ∂ρ_∂x::Vector{Matrix{T}}, - C::Matrix{Float64}, + C::Matrix{Float64}; eps = 1e-6, ) where {T<:Complex} if length(∂ρ_∂x) == 1 println( "In single parameter scenario, HCRB is equivalent to QFI. This function will return the value of QFI", ) - f = QFI(ρ, ∂ρ_∂x[1], eps) + f = QFIM_SLD(ρ, ∂ρ_∂x[1]; eps=eps) return f else Holevo_bound(ρ, ∂ρ_∂x, C; eps = eps) diff --git a/src/optim/optim.jl b/src/optim/optim.jl index 3c93225..382b20a 100644 --- a/src/optim/optim.jl +++ b/src/optim/optim.jl @@ -7,60 +7,81 @@ abstract type Rotation <: AbstractMeasurementType end abstract type Opt <: AbstractOpt end -mutable struct ControlOpt <: Opt - ctrl::AbstractVector - ctrl_bound::AbstractVector +Base.@kwdef mutable struct ControlOpt <: Opt + ctrl::Union{AbstractVector, Missing} = missing + ctrl_bound::AbstractVector = [-Inf, Inf] end -mutable struct StateOpt <: Opt - ψ0::AbstractVector +Copt = ControlOpt +ControlOpt(ctrl::Matrix{R}, ctrl_bound::AbstractVector) where R<:Number = ControlOpt([c[:] for c in eachrow(ctrl)], ctrl_bound) + +Base.@kwdef mutable struct StateOpt <: Opt + ψ₀::Union{AbstractVector, Missing} = missing end +Sopt = StateOpt + abstract type AbstractMopt <: Opt end -mutable struct Mopt_Projection <: AbstractMopt - C::AbstractVector +Base.@kwdef mutable struct Mopt_Projection <: AbstractMopt + C::Union{AbstractVector, Missing} = missing end - -mutable struct Mopt_LinearComb <: AbstractMopt - B::AbstractVector - POVM_basis::AbstractVector - M_num::Number +Base.@kwdef mutable struct Mopt_LinearComb <: AbstractMopt + B::Union{AbstractVector, Missing} = missing + POVM_basis::Union{AbstractVector, Missing} = missing + M_num::Int = 1 +end +Base.@kwdef mutable struct Mopt_Rotation <: AbstractMopt + s::Union{AbstractVector, Missing} = missing + POVM_basis::Union{AbstractVector, Missing} = missing + Lambda::Union{AbstractVector, Missing} = missing end -mutable struct Mopt_Rotation <: AbstractMopt - s::AbstractVector - POVM_basis::AbstractVector - Lambda::AbstractVector +function MeasurementOpt(;method=:Projection,kwargs...) + if method==:Projection + return Mopt_Projection(;kwargs...) + elseif method==:LinearCombination + return Mopt_LinearComb(;kwargs...) + elseif method==:Rotation + return Mopt_Rotation(;kwargs...) + end end +Mopt = MeasurementOpt abstract type CompOpt <: Opt end -mutable struct StateControlOpt <: CompOpt - ψ0::AbstractVector - ctrl::AbstractVector - ctrl_bound::AbstractVector +Base.@kwdef mutable struct StateControlOpt <: CompOpt + ψ₀::Union{AbstractVector, Missing} = missing + ctrl::Union{AbstractVector, Missing} = missing + ctrl_bound::AbstractVector = [-Inf, Inf] end -mutable struct ControlMeasurementOpt <: CompOpt - ctrl::AbstractVector - C::AbstractVector - ctrl_bound::AbstractVector -end +SCopt = StateControlOpt -mutable struct StateMeasurementOpt <: CompOpt - ψ0::AbstractVector - C::AbstractVector +Base.@kwdef mutable struct ControlMeasurementOpt <: CompOpt + ctrl::Union{AbstractVector, Missing} = missing + C::Union{AbstractVector, Missing} = missing + ctrl_bound::AbstractVector = [-Inf, Inf] +end + +CMopt = ControlMeasurementOpt + +Base.@kwdef mutable struct StateMeasurementOpt <: CompOpt + ψ₀::Union{AbstractVector, Missing} = missing + C::Union{AbstractVector, Missing} = missing end -mutable struct StateControlMeasurementOpt <: CompOpt - ctrl::AbstractVector - ψ0::AbstractVector - C::AbstractVector - ctrl_bound::AbstractVector +SMopt = StateMeasurementOpt + +Base.@kwdef mutable struct StateControlMeasurementOpt <: CompOpt + ctrl::Union{AbstractVector, Missing} = missing + ψ₀::Union{AbstractVector, Missing} = missing + C::Union{AbstractVector, Missing} = missing + ctrl_bound::AbstractVector = [-Inf, Inf] end -MeasurementOpt(M, mtype::Symbol = :Projection) = MeasurementOpt{eval(mtype)}(M) +SCMopt = StateControlMeasurementOpt + opt_target(::ControlOpt) = :Copt opt_target(::StateOpt) = :Sopt opt_target(::Mopt_Projection) = :Mopt @@ -73,14 +94,14 @@ opt_target(::StateMeasurementOpt) = :SMopt opt_target(::StateControlMeasurementOpt) = :SCMopt result(opt::ControlOpt) = [opt.ctrl] -result(opt::StateOpt) = [opt.ψ0] +result(opt::StateOpt) = [opt.ψ₀] result(opt::Mopt_Projection) = [opt.C] result(opt::Mopt_LinearComb) = [opt.B, opt.POVM_basis, opt.M_num] result(opt::Mopt_Rotation) = [opt.s] -result(opt::StateControlOpt) = [opt.ψ0, opt.ctrl] +result(opt::StateControlOpt) = [opt.ψ₀, opt.ctrl] result(opt::ControlMeasurementOpt) = [opt.ctrl, opt.C] -result(opt::StateMeasurementOpt) = [opt.ψ0, opt.C] -result(opt::StateControlMeasurementOpt) = [opt.ψ0, opt.ctrl, opt.C] +result(opt::StateMeasurementOpt) = [opt.ψ₀, opt.C] +result(opt::StateControlMeasurementOpt) = [opt.ψ₀, opt.ctrl, opt.C] #with reward result(opt, ::Type{Val{:save_reward}}) = [result(opt)..., [0.0]] diff --git a/src/output.jl b/src/output.jl index 08d4410..db22ecd 100644 --- a/src/output.jl +++ b/src/output.jl @@ -23,8 +23,8 @@ function set_io!(output::AbstractOutput, buffer...) end Output{T}(opt::AbstractOpt) where {T} = Output{T}([], [], res_file(opt), []) -Output(opt::AbstractOpt, save_type::Bool) = - save_type ? Output{savefile}(opt) : Output{no_save}(opt) +Output(opt::AbstractOpt; save::Bool=false) = + save ? Output{savefile}(opt) : Output{no_save}(opt) save_type(::Output{savefile}) = :savefile save_type(::Output{no_save}) = :no_save diff --git a/src/run.jl b/src/run.jl index 4a36dc0..84bc9c2 100644 --- a/src/run.jl +++ b/src/run.jl @@ -20,3 +20,17 @@ function run(system::QuanEstSystem) update!(optim, algorithm, obj, dynamics, output) show(obj, output) #io4 end + +function run(opt::AbstractOpt, alg::AbstractAlgorithm, obj::AbstractObj,dynamics::AbstractDynamics;savefile::Bool=false) + output = Output(opt;save=savefile) + obj = Objective(dynamics, obj) + system = QuanEstSystem(opt, alg, obj, dynamics, output) + run(system) +end + +# function run(opt::AbstractOpt, alg::AbstractAlgorithm, dynamics::AbstractDynamics; objective::Symbol, W=missing, M=missing, LD_type=:SLD, eps=GLOBAL_EPS, save::Bool=false) +# output = Output(opt;save=save) +# obj = Objective(opt, objective, W, M, LD_type, eps) +# system = QuanEstSystem(opt, alg, obj, dynamics, output) +# run(system) +# end \ No newline at end of file From 6d8e6c9544eece322dc5e1a8a93beab82a813f3b Mon Sep 17 00:00:00 2001 From: HuaiMing Date: Mon, 18 Apr 2022 18:16:04 +0800 Subject: [PATCH 2/2] revise, docs --- src/Project.toml | 28 ++ src/algorithm/AD.jl | 48 +-- src/algorithm/DE.jl | 40 +-- src/algorithm/GRAPE.jl | 28 +- src/algorithm/NM.jl | 2 +- src/algorithm/PSO.jl | 40 +-- src/algorithm/algorithm.jl | 32 +- src/common/adaptive.jl | 58 +++- src/common/common.jl | 42 +-- src/dynamics/kraus/KrausWrapper.jl | 24 +- src/dynamics/lindblad/LindbladDynamics.jl | 160 +++++++++- src/dynamics/lindblad/LindbladWrapper.jl | 110 ++++--- .../AsymptoticBound/AsymptoticBound.jl | 130 ++++---- .../AsymptoticBound/AsymptoticBoundWrapper.jl | 12 +- src/objective/AsymptoticBound/CramerRao.jl | 277 +++++++++++++++--- src/objective/AsymptoticBound/Holevo.jl | 41 ++- .../BayesianBound/BayesianCramerRao.jl | 71 ++++- src/optim/optim.jl | 39 +-- 18 files changed, 837 insertions(+), 345 deletions(-) create mode 100644 src/Project.toml 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"] diff --git a/src/algorithm/AD.jl b/src/algorithm/AD.jl index 4568cb8..793c8ac 100644 --- a/src/algorithm/AD.jl +++ b/src/algorithm/AD.jl @@ -28,18 +28,18 @@ function update!(opt::ControlOpt, alg::AbstractautoGRAPE, obj, dynamics, output) end function update_ctrl!(alg::autoGRAPE_Adam, obj, dynamics, δ) - (; ϵ, beta1, beta2) = alg + (; epsilon, beta1, beta2) = alg for ci in 1:length(δ) mt, vt = 0.0, 0.0 for ti in 1:length(δ[1]) dynamics.data.ctrl[ci][ti], mt, vt = Adam(δ[ci][ti], ti, - dynamics.data.ctrl[ci][ti], mt, vt, ϵ, beta1, beta2, obj.eps) + dynamics.data.ctrl[ci][ti], mt, vt, epsilon, beta1, beta2, obj.eps) end end end function update_ctrl!(alg::autoGRAPE, obj, dynamics, δ) - dynamics.data.ctrl += alg.ϵ*δ + dynamics.data.ctrl += alg.epsilon*δ end #### state optimization #### @@ -64,15 +64,15 @@ function update!(opt::StateOpt, alg::AbstractAD, obj, dynamics, output) end function update_state!(alg::AD_Adam, obj, dynamics, δ) - (; ϵ, beta1, beta2) = alg + (; epsilon, beta1, beta2) = alg mt, vt = 0.0, 0.0 for ti in 1:length(δ) - dynamics.data.ψ0[ti], mt, vt = Adam(δ[ti], ti, dynamics.data.ψ0[ti], mt, vt, ϵ, beta1, beta2, obj.eps) + dynamics.data.ψ0[ti], mt, vt = Adam(δ[ti], ti, dynamics.data.ψ0[ti], mt, vt, epsilon, beta1, beta2, obj.eps) end end function update_state!(alg::AD, obj, dynamics, δ) - dynamics.data.ψ0 += alg.ϵ*δ + dynamics.data.ψ0 += alg.epsilon*δ end #### find the optimal linear combination of a given set of POVM #### @@ -83,7 +83,7 @@ function update!(opt::Mopt_LinearComb, alg::AbstractAD, obj, dynamics, output) bound_LC_coeff!(opt.B) M = [sum([opt.B[i][j]*POVM_basis[j] for j in 1:basis_num]) for i in 1:M_num] - obj_QFIM = QFIM_Obj(obj) + obj_QFIM = QFIM_obj(obj) f_opt, f_comp = objective(obj_QFIM, dynamics) obj_POVM = set_M(obj, POVM_basis) f_povm, f_comp = objective(obj_POVM, dynamics) @@ -109,17 +109,17 @@ function update!(opt::Mopt_LinearComb, alg::AbstractAD, obj, dynamics, output) end function update_M!(opt::Mopt_LinearComb, alg::AD_Adam, obj, δ) - (; ϵ, beta1, beta2) = alg + (; epsilon, beta1, beta2) = alg for ci in 1:length(δ) mt, vt = 0.0, 0.0 for ti in 1:length(δ[1]) - opt.B[ci][ti], mt, vt = Adam(δ[ci][ti], ti, opt.B[ci][ti], mt, vt, ϵ, beta1, beta2, obj.eps) + opt.B[ci][ti], mt, vt = Adam(δ[ci][ti], ti, opt.B[ci][ti], mt, vt, epsilon, beta1, beta2, obj.eps) end end end function update_M!(opt::Mopt_LinearComb, alg::AD, obj, δ) - opt.B += alg.ϵ*δ + opt.B += alg.epsilon*δ end #### find the optimal rotated measurement of a given set of POVM #### @@ -129,15 +129,19 @@ function update!(opt::Mopt_Rotation, alg::AbstractAD, obj, dynamics, output) dim = size(dynamics.data.ρ0)[1] M_num = length(POVM_basis) suN = suN_generator(dim) - if ismissing(opt.Lambda) - opt.Lambda = Matrix{ComplexF64}[] - append!(opt.Lambda, [Matrix{ComplexF64}(I,dim,dim)]) - append!(opt.Lambda, [suN[i] for i in 1:length(suN)]) - end + opt.Lambda = Matrix{ComplexF64}[] + append!(opt.Lambda, [Matrix{ComplexF64}(I,dim,dim)]) + append!(opt.Lambda, [suN[i] for i in 1:length(suN)]) + + # if ismissing(Lambda) + # opt.Lambda = Matrix{ComplexF64}[] + # append!(opt.Lambda, [Matrix{ComplexF64}(I,dim,dim)]) + # append!(opt.Lambda, [suN[i] for i in 1:length(suN)]) + # end U = rotation_matrix(opt.s, opt.Lambda) M = [U*POVM_basis[i]*U' for i in 1:M_num] - obj_QFIM = QFIM_Obj(obj) + obj_QFIM = QFIM_obj(obj) f_opt, f_comp = objective(obj_QFIM, dynamics) obj_POVM = set_M(obj, POVM_basis) f_povm, f_comp = objective(obj_POVM, dynamics) @@ -164,15 +168,15 @@ function update!(opt::Mopt_Rotation, alg::AbstractAD, obj, dynamics, output) end function update_M!(opt::Mopt_Rotation, alg::AD_Adam, obj, δ) - (; ϵ, beta1, beta2) = alg + (; epsilon, beta1, beta2) = alg mt, vt = 0.0, 0.0 for ti in 1:length(δ) - opt.s[ti], mt, vt = Adam(δ[ti], ti, opt.s[ti], mt, vt, ϵ, beta1, beta2, obj.eps) + opt.s[ti], mt, vt = Adam(δ[ti], ti, opt.s[ti], mt, vt, epsilon, beta1, beta2, obj.eps) end end function update_M!(opt::Mopt_Rotation, alg::AD, obj, δ) - opt.s += alg.ϵ*δ + opt.s += alg.epsilon*δ end #### state abd control optimization #### @@ -207,16 +211,16 @@ function update!(opt::StateControlOpt, alg::AbstractAD, obj, dynamics, output) end function update_ctrl!(alg::AD_Adam, obj, dynamics, δ) - (; ϵ, beta1, beta2) = alg + (; epsilon, beta1, beta2) = alg for ci in 1:length(δ) mt, vt = 0.0, 0.0 for ti in 1:length(δ[1]) dynamics.data.ctrl[ci][ti], mt, vt = Adam(δ[ci][ti], ti, - dynamics.data.ctrl[ci][ti], mt, vt, ϵ, beta1, beta2, obj.eps) + dynamics.data.ctrl[ci][ti], mt, vt, epsilon, beta1, beta2, obj.eps) end end end function update_ctrl!(alg::AD, obj, dynamics, δ) - dynamics.data.ctrl += alg.ϵ*δ + dynamics.data.ctrl += alg.epsilon*δ end diff --git a/src/algorithm/DE.jl b/src/algorithm/DE.jl index cca7654..fc2c430 100644 --- a/src/algorithm/DE.jl +++ b/src/algorithm/DE.jl @@ -80,7 +80,7 @@ end function update!(opt::StateOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg if ismissing(ini_population) - ini_population = ([opt.ψ₀,],) + ini_population = ([opt.psi,],) end ini_population = ini_population[1] dim = length(dynamics.data.ψ0) @@ -143,12 +143,12 @@ end function update!(opt::Mopt_Projection, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg if ismissing(ini_population) - ini_population = ([opt.C], ) + ini_population = ([opt.M], ) end ini_population = ini_population[1] dim = size(dynamics.data.ρ0)[1] - M_num = length(opt.C) + M_num = length(opt.M) populations = [[zeros(ComplexF64, dim) for j in 1:M_num] for i in 1:p_num] # initialization @@ -161,7 +161,7 @@ function update!(opt::Mopt_Projection, alg::DE, obj, dynamics, output) p_out[pj], p_fit[pj] = objective(obj_copy, dynamics) end - obj_QFIM = QFIM_Obj(obj) + obj_QFIM = QFIM_obj(obj) f_opt, f_comp = objective(obj_QFIM, dynamics) M = [populations[1][i]*(populations[1][i])' for i in 1:M_num] @@ -242,7 +242,7 @@ function update!(opt::Mopt_LinearComb, alg::DE, obj, dynamics, output) p_out[pj], p_fit[pj] = objective(obj_copy, dynamics) end - obj_QFIM = QFIM_Obj(obj) + obj_QFIM = QFIM_obj(obj) f_opt, f_comp = objective(obj_QFIM, dynamics) obj_POVM = set_M(obj, POVM_basis) f_povm, f_comp = objective(obj_POVM, dynamics) @@ -314,11 +314,15 @@ function update!(opt::Mopt_Rotation, alg::DE, obj, dynamics, output) ini_population = ini_population[1] dim = size(dynamics.data.ρ0)[1] suN = suN_generator(dim) - if ismissing(Lambda) - Lambda = Matrix{ComplexF64}[] - append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) - append!(Lambda, [suN[i] for i in 1:length(suN)]) - end + Lambda = Matrix{ComplexF64}[] + append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) + append!(Lambda, [suN[i] for i in 1:length(suN)]) + + # if ismissing(Lambda) + # Lambda = Matrix{ComplexF64}[] + # append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) + # append!(Lambda, [suN[i] for i in 1:length(suN)]) + # end M_num = length(POVM_basis) populations = [zeros(dim^2) for i in 1:p_num] @@ -333,7 +337,7 @@ function update!(opt::Mopt_Rotation, alg::DE, obj, dynamics, output) p_out[pj], p_fit[pj] = objective(obj_copy, dynamics) end - obj_QFIM = QFIM_Obj(obj) + obj_QFIM = QFIM_obj(obj) f_opt, f_comp = objective(obj_QFIM, dynamics) obj_POVM = set_M(obj, POVM_basis) f_povm, f_comp = objective(obj_POVM, dynamics) @@ -397,7 +401,7 @@ end function update!(opt::StateControlOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg if ismissing(ini_population) - ini_population = ([opt.ψ₀], [opt.ctrl,]) + ini_population = ([opt.psi], [opt.ctrl,]) end psi0, ctrl0 = ini_population ctrl_length = length(dynamics.data.ctrl[1]) @@ -496,11 +500,11 @@ end function update!(opt::StateMeasurementOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg if ismissing(ini_population) - ini_population = ( [opt.ψ₀], [opt.C,]) + ini_population = ([opt.psi], [opt.M,]) end psi0, measurement0 = ini_population dim = length(dynamics.data.ψ0) - M_num = length(opt.C) + M_num = length(opt.M) populations = repeat(dynamics, p_num) # initialization @@ -599,13 +603,13 @@ end function update!(opt::ControlMeasurementOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg if ismissing(ini_population) - ini_population = ([opt.ctrl,], [opt.C]) + ini_population = ([opt.ctrl,], [opt.M]) end ctrl0, measurement0 = ini_population dim = size(dynamics.data.ρ0)[1] ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) - M_num = length(opt.C) + M_num = length(opt.M) populations = repeat(dynamics, p_num) # initialization @@ -713,13 +717,13 @@ end function update!(opt::StateControlMeasurementOpt, alg::DE, obj, dynamics, output) (; max_episode, p_num, ini_population, c, cr, rng) = alg if ismissing(ini_population) - ini_population = ([opt.ψ₀], [opt.ctrl,], [opt.C]) + ini_population = ([opt.psi], [opt.ctrl,], [opt.M]) end psi0, ctrl0, measurement0 = ini_population dim = length(dynamics.data.ψ0) ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) - M_num = length(opt.C) + M_num = length(opt.M) populations = repeat(dynamics, p_num) # initialization diff --git a/src/algorithm/GRAPE.jl b/src/algorithm/GRAPE.jl index 95d524f..5bd03af 100644 --- a/src/algorithm/GRAPE.jl +++ b/src/algorithm/GRAPE.jl @@ -1,4 +1,4 @@ -function update!(opt::ControlOpt, alg::AbstractGRAPE, obj::QFIM_Obj, dynamics, output) +function update!(opt::ControlOpt, alg::AbstractGRAPE, obj::QFIM_obj, dynamics, output) (; max_episode) = alg ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) @@ -22,7 +22,7 @@ function update!(opt::ControlOpt, alg::AbstractGRAPE, obj::QFIM_Obj, dynamics, o set_io!(output, output.f_list[end]) end -function update!(opt::ControlOpt, alg::AbstractGRAPE, obj::CFIM_Obj, dynamics, output) +function update!(opt::ControlOpt, alg::AbstractGRAPE, obj::CFIM_obj, dynamics, output) (; max_episode) = alg ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) @@ -124,7 +124,7 @@ function gradient_QFIM_analy(opt, alg::GRAPE_Adam, obj, dynamics) term1 = tr(∂xδρt_T_δV*Lx[1]) term2 = tr(∂ρt_T_δV*anti_commu) δF = ((2*term1-0.5*term2) |> real) - dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.ϵ, alg.beta1, alg.beta2, obj.eps) + dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.epsilon, alg.beta1, alg.beta2, obj.eps) end end elseif para_num == 2 @@ -151,7 +151,7 @@ function gradient_QFIM_analy(opt, alg::GRAPE_Adam, obj, dynamics) item1 = -coeff2*(F_T[2,2]*δF_all[1][1]+F_T[1,1]*δF_all[2][2]-F_T[2,1]*δF_all[1][2]-F_T[1,2]*δF_all[2][1])/coeff1^2 item2 = (obj.W[1,1]*δF_all[2][2]+obj.W[2,2]*δF_all[1][1]-obj.W[1,2]*δF_all[2][1]-obj.W[2,1]*δF_all[1][2])/coeff1 δF = -(item1+item2)*cost_function^2 - dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.ϵ, alg.beta1, alg.beta2, obj.eps) + dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.epsilon, alg.beta1, alg.beta2, obj.eps) end end else @@ -171,7 +171,7 @@ function gradient_QFIM_analy(opt, alg::GRAPE_Adam, obj, dynamics) δF = δF + obj.W[pm,pm]*(1.0/F_T[pm,pm]/F_T[pm,pm])*((2*term1-0.5*term2) |> real) end δF = δF*coeff - dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.ϵ, alg.beta1, alg.beta2, obj.eps) + dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.epsilon, alg.beta1, alg.beta2, obj.eps) end end end @@ -199,7 +199,7 @@ function gradient_QFIM_analy(opt, alg::GRAPE, obj, dynamics) term1 = tr(∂xδρt_T_δV*Lx[1]) term2 = tr(∂ρt_T_δV*anti_commu) δF = ((2*term1-0.5*term2) |> real) - dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.ϵ*δF + dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.epsilon*δF end end elseif para_num == 2 @@ -225,7 +225,7 @@ function gradient_QFIM_analy(opt, alg::GRAPE, obj, dynamics) item1 = -coeff2*(F_T[2,2]*δF_all[1][1]+F_T[1,1]*δF_all[2][2]-F_T[2,1]*δF_all[1][2]-F_T[1,2]*δF_all[2][1])/coeff1^2 item2 = (obj.W[1,1]*δF_all[2][2]+obj.W[2,2]*δF_all[1][1]-obj.W[1,2]*δF_all[2][1]-obj.W[2,1]*δF_all[1][2])/coeff1 δF = -(item1+item2)*cost_function^2 - dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.ϵ*δF + dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.epsilon*δF end end else @@ -244,7 +244,7 @@ function gradient_QFIM_analy(opt, alg::GRAPE, obj, dynamics) δF = δF + obj.W[pm,pm]*(1.0/F_T[pm,pm]/F_T[pm,pm])*((2*term1-0.5*term2) |> real) end δF = δF*coeff - dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.ϵ*δF + dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.epsilon*δF end end end @@ -281,7 +281,7 @@ function gradient_CFIM_analy_Adam(opt, alg, obj, dynamics) term1 = tr(∂xδρt_T_δV*L1_tidle) term2 = tr(∂ρt_T_δV*L2_tidle) δF = ((2*term1-term2) |> real) - dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.ϵ, alg.beta1, alg.beta2, obj.eps) + dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.epsilon, alg.beta1, alg.beta2, obj.eps) end end elseif para_num == 2 @@ -332,7 +332,7 @@ function gradient_CFIM_analy_Adam(opt, alg, obj, dynamics) item1 = -coeff2*(F_T[2,2]*δF_all[1][1]+F_T[1,1]*δF_all[2][2]-F_T[2,1]*δF_all[1][2]-F_T[1,2]*δF_all[2][1])/coeff1^2 item2 = (obj.W[1,1]*δF_all[2][2]+obj.W[2,2]*δF_all[1][1]-obj.W[1,2]*δF_all[2][1]-obj.W[2,1]*δF_all[1][2])/coeff1 δF = -(item1+item2)*cost_function^2 - dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.ϵ, alg.beta1, alg.beta2, obj.eps) + dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.epsilon, alg.beta1, alg.beta2, obj.eps) end end else @@ -367,7 +367,7 @@ function gradient_CFIM_analy_Adam(opt, alg, obj, dynamics) δF = δF + obj.W[pm,pm]*(1.0/F_T[pm,pm]/F_T[pm,pm])*((2*term1-term2) |> real) end δF = δF*coeff - dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.ϵ, alg.beta1, alg.beta2, obj.eps) + dynamics.data.ctrl[cm][tm], mt, vt = Adam(δF, tm, dynamics.data.ctrl[cm][tm], mt, vt, alg.epsilon, alg.beta1, alg.beta2, obj.eps) end end end @@ -403,7 +403,7 @@ function gradient_CFIM_analy(opt, alg, obj, dynamics) term1 = tr(∂xδρt_T_δV*L1_tidle) term2 = tr(∂ρt_T_δV*L2_tidle) δF = ((2*term1-term2) |> real) - dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.ϵ*δF + dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.epsilon*δF end end elseif para_num == 2 @@ -453,7 +453,7 @@ function gradient_CFIM_analy(opt, alg, obj, dynamics) item1 = -coeff2*(F_T[2,2]*δF_all[1][1]+F_T[1,1]*δF_all[2][2]-F_T[2,1]*δF_all[1][2]-F_T[1,2]*δF_all[2][1])/coeff1^2 item2 = (obj.W[1,1]*δF_all[2][2]+obj.W[2,2]*δF_all[1][1]-obj.W[1,2]*δF_all[2][1]-obj.W[2,1]*δF_all[1][2])/coeff1 δF = -(item1+item2)*cost_function^2 - dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.ϵ*δF + dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.epsilon*δF end end else @@ -486,7 +486,7 @@ function gradient_CFIM_analy(opt, alg, obj, dynamics) δF = δF + obj.W[pm,pm]*(1.0/F_T[pm,pm]/F_T[pm,pm])*((2*term1-term2) |> real) end δF = δF*coeff - dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.ϵ*δF + dynamics.data.ctrl[cm][tm] = dynamics.data.ctrl[cm][tm] + alg.epsilon*δF end end end diff --git a/src/algorithm/NM.jl b/src/algorithm/NM.jl index 12f6b0b..014db42 100644 --- a/src/algorithm/NM.jl +++ b/src/algorithm/NM.jl @@ -1,7 +1,7 @@ function update!(opt::StateOpt, alg::NM, obj, dynamics, output) (; max_episode, state_num, ini_state, ar, ae, ac, as0, rng) = alg if ismissing(ini_state) - ini_state = [opt.ψ₀] + ini_state = [opt.psi] end dim = length(dynamics.data.ψ0) nelder_mead = repeat(dynamics, state_num) diff --git a/src/algorithm/PSO.jl b/src/algorithm/PSO.jl index a32a8c1..e83d9ed 100644 --- a/src/algorithm/PSO.jl +++ b/src/algorithm/PSO.jl @@ -104,7 +104,7 @@ end function update!(opt::StateOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg if ismissing(ini_particle) - ini_particle = ([opt.ψ₀], ) + ini_particle = ([opt.psi], ) end ini_particle = ini_particle[1] dim = length(dynamics.data.ψ0) @@ -181,11 +181,11 @@ end function update!(opt::Mopt_Projection, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg if ismissing(ini_particle) - ini_particle = ([opt.C], ) + ini_particle = ([opt.M], ) end ini_particle = ini_particle[1] dim = size(dynamics.data.ρ0)[1] - M_num = length(opt.C) + M_num = length(opt.M) particles = [[zeros(ComplexF64, dim) for j in 1:M_num] for i in 1:p_num] if typeof(max_episode) == Int @@ -206,7 +206,7 @@ function update!(opt::Mopt_Projection, alg::PSO, obj, dynamics, output) obj_copy = set_M(obj, M) f_ini, f_comp = objective(obj_copy, dynamics) - obj_QFIM = QFIM_Obj(obj) + obj_QFIM = QFIM_obj(obj) f_opt, f_comp = objective(obj_QFIM, dynamics) set_f!(output, f_ini) @@ -295,7 +295,7 @@ function update!(opt::Mopt_LinearComb, alg::PSO, obj, dynamics, output) # initialization initial_LinearComb!(ini_particle, particles, basis_num, M_num, p_num, rng) - obj_QFIM = QFIM_Obj(obj) + obj_QFIM = QFIM_obj(obj) f_opt, f_comp = objective(obj_QFIM, dynamics) obj_POVM = set_M(obj, POVM_basis) f_povm, f_comp = objective(obj_POVM, dynamics) @@ -376,11 +376,15 @@ function update!(opt::Mopt_Rotation, alg::PSO, obj, dynamics, output) M_num = length(POVM_basis) dim = size(dynamics.data.ρ0)[1] suN = suN_generator(dim) - if ismissing(Lambda) - Lambda = Matrix{ComplexF64}[] - append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) - append!(Lambda, [suN[i] for i in 1:length(suN)]) - end + Lambda = Matrix{ComplexF64}[] + append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) + append!(Lambda, [suN[i] for i in 1:length(suN)]) + + # if ismissing(Lambda) + # Lambda = Matrix{ComplexF64}[] + # append!(Lambda, [Matrix{ComplexF64}(I,dim,dim)]) + # append!(Lambda, [suN[i] for i in 1:length(suN)]) + # end particles = repeat(s, p_num) @@ -399,7 +403,7 @@ function update!(opt::Mopt_Rotation, alg::PSO, obj, dynamics, output) particles = [zeros(dim^2) for i in 1:p_num] initial_Rotation!(ini_particle, particles, dim, p_num, rng) - obj_QFIM = QFIM_Obj(obj) + obj_QFIM = QFIM_obj(obj) f_opt, f_comp = objective(obj_QFIM, dynamics) obj_POVM = set_M(obj, POVM_basis) f_povm, f_comp = objective(obj_POVM, dynamics) @@ -468,7 +472,7 @@ end function update!(opt::StateControlOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg if ismissing(ini_particle) - ini_particle = ([opt.ψ₀], [opt.ctrl,]) + ini_particle = ([opt.psi], [opt.ctrl,]) end psi0, ctrl0 = ini_particle dim = length(dynamics.data.ψ0) @@ -579,11 +583,11 @@ end function update!(opt::StateMeasurementOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg if ismissing(ini_particle) - ini_particle = ([opt.ψ₀], [opt.C]) + ini_particle = ([opt.psi], [opt.M]) end psi0, measurement0 = ini_particle dim = length(dynamics.data.ψ0) - M_num = length(opt.C) + M_num = length(opt.M) particles = repeat(dynamics, p_num) if typeof(max_episode) == Int @@ -694,13 +698,13 @@ end function update!(opt::ControlMeasurementOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg if ismissing(ini_particle) - ini_particle = ([opt.ctrl,], [opt.C]) + ini_particle = ([opt.ctrl,], [opt.M]) end ctrl0, measurement0 = ini_particle ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) dim = size(dynamics.data.ρ0)[1] - M_num = length(opt.C) + M_num = length(opt.M) particles = repeat(dynamics, p_num) if typeof(max_episode) == Int @@ -817,13 +821,13 @@ end function update!(opt::StateControlMeasurementOpt, alg::PSO, obj, dynamics, output) (; max_episode, p_num, ini_particle, c0, c1, c2, rng) = alg if ismissing(ini_particle) - ini_particle = ([opt.ψ₀], [opt.ctrl,], [opt.C]) + ini_particle = ([opt.psi], [opt.ctrl,], [opt.M]) end psi0, ctrl0, measurement0 = ini_particle ctrl_length = length(dynamics.data.ctrl[1]) ctrl_num = length(dynamics.data.Hc) dim = length(dynamics.data.ψ0) - M_num = length(opt.C) + M_num = length(opt.M) particles = repeat(dynamics, p_num) if typeof(max_episode) == Int diff --git a/src/algorithm/algorithm.jl b/src/algorithm/algorithm.jl index 7656120..6919662 100644 --- a/src/algorithm/algorithm.jl +++ b/src/algorithm/algorithm.jl @@ -3,50 +3,50 @@ abstract type AbstractAlgorithm end abstract type AbstractGRAPE <: AbstractAlgorithm end Base.@kwdef struct GRAPE <: AbstractGRAPE max_episode::Int = 300 - ϵ::Number = 0.01 + epsilon::Number = 0.01 end Base.@kwdef struct GRAPE_Adam <: AbstractGRAPE max_episode::Int = 300 - ϵ::Number = 0.01 + epsilon::Number = 0.01 beta1::Number = 0.90 beta2::Number = 0.99 end -GRAPE(max_episode, ϵ, beta1, beta2) = GRAPE_Adam(max_episode, ϵ, beta1, beta2) -GRAPE(;max_episode=300, ϵ=0.01, beta1=0.90, beta2=0.99, Adam::Bool=true) = Adam ? GRAPE_Adam(max_episode, ϵ, beta1, beta2) : GRAPE(max_episode, ϵ) +GRAPE(max_episode, epsilon, beta1, beta2) = GRAPE_Adam(max_episode, epsilon, beta1, beta2) +GRAPE(;max_episode=300, epsilon=0.01, beta1=0.90, beta2=0.99, Adam::Bool=true) = Adam ? GRAPE_Adam(max_episode, epsilon, beta1, beta2) : GRAPE(max_episode, epsilon) abstract type AbstractautoGRAPE <: AbstractAlgorithm end Base.@kwdef struct autoGRAPE <: AbstractautoGRAPE max_episode::Int = 300 - ϵ::Number = 0.01 + epsilon::Number = 0.01 end Base.@kwdef struct autoGRAPE_Adam <: AbstractautoGRAPE max_episode::Int = 300 - ϵ::Number = 0.01 + epsilon::Number = 0.01 beta1::Number = 0.90 beta2::Number = 0.99 end -autoGRAPE(max_episode, ϵ, beta1, beta2) = autoGRAPE_Adam(max_episode, ϵ, beta1, beta2) -autoGRAPE(;max_episode=300, ϵ=0.01, beta1=0.90, beta2=0.99, Adam::Bool=true) = Adam ? autoGRAPE_Adam(max_episode, ϵ, beta1, beta2) : autoGRAPE(max_episode, ϵ) +autoGRAPE(max_episode, epsilon, beta1, beta2) = autoGRAPE_Adam(max_episode, epsilon, beta1, beta2) +autoGRAPE(;max_episode=300, epsilon=0.01, beta1=0.90, beta2=0.99, Adam::Bool=true) = Adam ? autoGRAPE_Adam(max_episode, epsilon, beta1, beta2) : autoGRAPE(max_episode, epsilon) abstract type AbstractAD <: AbstractAlgorithm end Base.@kwdef struct AD <: AbstractAD max_episode::Number = 300 - ϵ::Number = 0.01 + epsilon::Number = 0.01 end Base.@kwdef struct AD_Adam <: AbstractAD max_episode::Number = 300 - ϵ::Number = 0.01 + epsilon::Number = 0.01 beta1::Number = 0.90 beta2::Number = 0.99 end -AD(max_episode, ϵ, beta1, beta2) = AD_Adam(max_episode=max_episode, ϵ=ϵ, beta1=beta1, beta2=beta2) -AD(;max_episode=300, ϵ=0.01, beta1=0.90, beta2=0.99, Adam::Bool=true) = Adam ? AD_Adam(max_episode, ϵ, beta1, beta2) : AD(max_episode, ϵ) +AD(max_episode, epsilon, beta1, beta2) = AD_Adam(max_episode=max_episode, epsilon=epsilon, beta1=beta1, beta2=beta2) +AD(;max_episode=300, epsilon=0.01, beta1=0.90, beta2=0.99, Adam::Bool=true) = Adam ? AD_Adam(max_episode, epsilon, beta1, beta2) : AD(max_episode, epsilon) Base.@kwdef struct PSO <: AbstractAlgorithm max_episode::Union{T,Vector{T}} where {T<:Int} = [1000, 100] @@ -62,6 +62,8 @@ PSO(max_episode, p_num, ini_particle, c0, c1, c2) = PSO(max_episode, p_num, ini_particle, c0, c1, c2, GLOBAL_RNG) PSO(max_episode, p_num, ini_particle, c0, c1, c2, seed::Number) = PSO(max_episode, p_num, ini_particle, c0, c1, c2, MersenneTwister(seed)) +PSO(;max_episode::Union{T,Vector{T}} where {T<:Int}=[1000, 100], p_num::Number=10, ini_particle=missing, c0::Number=1.0, c1::Number=2.0, c2::Number=2.0, seed::Number=1234) = + PSO(max_episode, p_num, ini_particle, c0, c1, c2, MersenneTwister(seed)) Base.@kwdef struct DE <: AbstractAlgorithm max_episode::Number = 1000 @@ -76,6 +78,8 @@ DE(max_episode, p_num, ini_population, c, cr) = DE(max_episode, p_num, ini_population, c, cr, GLOBAL_RNG) DE(max_episode, p_num, ini_population, c, cr, seed::Number) = DE(max_episode, p_num, ini_population, c, cr, MersenneTwister(seed)) +DE(;max_episode::Number=1000, p_num::Number=10, ini_population=missing, c::Number=1.0, cr::Number=0.5, seed::Number=1234) = + DE(max_episode, p_num, ini_population, c, cr, MersenneTwister(seed)) Base.@kwdef struct DDPG <: AbstractAlgorithm max_episode::Int = 500 @@ -88,6 +92,8 @@ DDPG(max_episode, layer_num, layer_dim) = DDPG(max_episode, layer_num, layer_dim, GLOBAL_RNG) DDPG(max_episode, layer_num, layer_dim, seed::Number) = DDPG(max_episode, layer_num, layer_dim, StableRNG(seed)) +DDPG(;max_episode::Int=500, layer_num::Int=3, layer_dim::Int=200, seed::Number=1234) = + DDPG(max_episode, layer_num, layer_dim, StableRNG(seed)) Base.@kwdef struct NM <: AbstractAlgorithm max_episode::Int = 1000 @@ -104,6 +110,8 @@ NM(max_episode, state_num, nelder_mead, ar, ae, ac, as0) = NM(max_episode, state_num, nelder_mead, ar, ae, ac, as0, GLOBAL_RNG) NM(max_episode, state_num, nelder_mead, ar, ae, ac, as0, seed::Number) = NM(max_episode, state_num, nelder_mead, ar, ae, ac, as0, MersenneTwister(seed)) +NM(;max_episode::Int=1000, state_num::Int=10, nelder_mead=missing, ar::Number=1.0, ae::Number=2.0, ac::Number=0.5, as0::Number=0.5, seed::Number=1234) = + NM(max_episode, state_num, nelder_mead, ar, ae, ac, as0, MersenneTwister(seed)) alg_type(::AD) = :AD alg_type(::AD_Adam) = :AD diff --git a/src/common/adaptive.jl b/src/common/adaptive.jl index 98d4edd..70ac11d 100644 --- a/src/common/adaptive.jl +++ b/src/common/adaptive.jl @@ -1,6 +1,6 @@ -function adaptive(x::AbstractVector, p, rho0::AbstractMatrix, tspan, H, dH; decay::Union{Vector,Nothing}=nothing, - Hc::Union{Vector,Nothing}=nothing, ctrl::Union{Vector,Nothing}=nothing, M::Union{AbstractVector,Nothing}=nothing, - W::Union{Matrix,Nothing}=nothing, max_episode::Int=1000, eps::Float64=1e-8, save_file=false) +function adaptive(x::AbstractVector, p, rho0::AbstractMatrix, tspan, H, dH; save_file=false, max_episode::Int=1000, eps::Float64=1e-8, + Hc::Union{Vector,Nothing}=nothing, ctrl::Union{Vector,Nothing}=nothing, decay::Union{Vector,Nothing}=nothing, + M::Union{AbstractVector,Nothing}=nothing, W::Union{Matrix,Nothing}=nothing) dim = size(rho0)[1] para_num = length(x) if M == nothing @@ -193,8 +193,8 @@ function adaptive(x::AbstractVector, p, rho0::AbstractMatrix, tspan, H, dH; deca end -function adaptive(x::AbstractVector, p, rho0::AbstractMatrix, K, dK; M::Union{AbstractVector,Nothing}=nothing, - W::Union{Matrix,Nothing}=nothing, max_episode::Int=1000, eps::Float64=1e-8, save_file=false) +function adaptive(x::AbstractVector, p, rho0::AbstractMatrix, K, dK; save_file=false, max_episode::Int=1000, + eps::Float64=1e-8, M::Union{AbstractVector,Nothing}=nothing, W::Union{Matrix,Nothing}=nothing) dim = size(rho0)[1] para_num = length(x) @@ -359,6 +359,50 @@ function adaptive(x::AbstractVector, p, rho0::AbstractMatrix, K, dK; M::Union{Ab end end +mutable struct adaptMZI + x + p + rho0 +end + +function online(apt::adaptMZI; output::String = "phi") + (;x, p, rho0) = apt + N = Int(sqrt(size(rho0,1))) - 1 + a = destroy(N+1)' |> sparse + adaptMZI_online(x, p, rho0, a, output) +end + +function brgd(n) + if n == 1 + return ["0", "1"] + end + L0 = brgd(n-1) + L1 = deepcopy(L0) + reverse!(L1) + L0 = ["0"*l for l in L0] + L1 = ["1"*l for l in L1] + return deepcopy(vcat(L0,L1)) +end + +function offline(apt::adaptMZI, alg;eps = GLOBAL_EPS) + (;x,p,rho0) = apt + N = Int(sqrt(size(rho0,1))) - 1 + a = destroy(N+1)' |> sparse + comb = brgd(N)|>x->[[parse(Int, s) for s in ss] for ss in x] + if alg isa DE + (;p_num,ini_population,c,cr,rng,max_episode) = alg + if ismissing(ini_population) + ini_population = ([apt.rho0],) + end + DE_DeltaPhiOpt(x,p,rho0,a,comb,p_num,ini_population[1],c,cr,rng.seed,max_episode,eps) + elseif alg isa PSO + (;p_num,ini_particle,c0,c1,c2,rng,max_episode) = alg + if ismissing(ini_particle) + ini_particle = ([apt.rho0],) + end + PSO_DeltaPhiOpt(x,p,rho0,a,comb,p_num,ini_particle[1],c0,c1,c2,rng.seed,max_episode,eps) + end +end function adaptMZI_online(x, p, rho0, a, output) @@ -376,7 +420,7 @@ function adaptMZI_online(x, p, rho0, a, output) enter = readline() u = parse(Int64, enter) - pyx = zeros(length(x)) + pyx = zeros(length(x))|>sparse for xi in 1:length(x) a_res_tp = a_res[xi]*a_u(a, x[xi], phi, u) pyx[xi] = real(tr(rho0*a_res_tp'*a_res_tp))*(factorial(N-ei)/factorial(N)) @@ -462,7 +506,7 @@ function adaptMZI_offline(delta_phi, x, p, rho0, a, comb, eps) exp_ix = [exp(1.0im*xi) for xi in x] M_res = zeros(length(comb)) - Threads.@threads for ui in 1:length(comb) + for ui in 1:length(comb) u = comb[ui] phi = 0.0 diff --git a/src/common/common.jl b/src/common/common.jl index 3011bce..a1e451d 100644 --- a/src/common/common.jl +++ b/src/common/common.jl @@ -24,7 +24,7 @@ function Base.repeat(system, N) end function Base.repeat(system, M, N) - reshape(repeat(system, M * N), M,N) + reshape(repeat(system, M * N), M, N) end function filterZeros!(x::Matrix{T}) where {T <: Complex} @@ -152,12 +152,12 @@ function AdaptiveInput(x, func, dfunc; channel="dynamics") H = [func(xi) for xi in x_list] dH = [dfunc(xi) for xi in x_list] return H, dH - elseif channel == "kraus" + elseif channel == "Kraus" K = [func(xi) for xi in x_list] dK = [dfunc(xi) for xi in x_list] return K, dK else - throw("Supported values for channel are 'dynamics' and 'kraus'") + throw("Supported values for channel are 'dynamics' and 'Kraus'") end end @@ -176,13 +176,13 @@ function bound!(ctrl::Vector{Float64}, ctrl_bound) end end -function Adam(gt, t, para, mt, vt, ϵ, beta1, beta2, eps) +function Adam(gt, t, para, mt, vt, epsilon, beta1, beta2, eps) t = t+1 mt = beta1*mt + (1-beta1)*gt vt = beta2*vt + (1-beta2)*(gt*gt) m_cap = mt/(1-(beta1^t)) v_cap = vt/(1-(beta2^t)) - para = para+(ϵ*m_cap)/(sqrt(v_cap)+eps) + para = para+(epsilon*m_cap)/(sqrt(v_cap)+eps) return para, mt, vt end @@ -277,7 +277,7 @@ function initial_ctrl!(opt, ctrl0, dynamics, p_num, rng) ctrl0 = [ctrl0[i] for i in 1:p_num] end for pj in 1:length(ctrl0) - dynamics[pj].data.ctrl = deepcopy(ctrl0[pj]) #[[ctrl0[pj][i, j] for j in 1:ctrl_length] for i in 1:ctrl_num] + dynamics[pj].data.ctrl = ctrl0[pj] isa AbstractVector ? deepcopy(ctrl0[pj]) : [[ctrl0[pj][i, j] for j in 1:ctrl_length] for i in 1:ctrl_num] end if opt.ctrl_bound[1] == -Inf || opt.ctrl_bound[2] == Inf for pj in (length(ctrl0)+1):p_num @@ -310,7 +310,7 @@ function initial_M!(measurement0, C_all, dim, p_num, M_num, rng) measurement0 = [measurement0[i] for i in 1:p_num] end for pj in 1:length(measurement0) - C_all[pj] = deepcopy(measurement0[pj]) #[[measurement0[pj][i,j] for j in 1:dim] for i in 1:M_num] + C_all[pj] = measurement0[pj] isa AbstractVector ? deepcopy(measurement0[pj]) : [[measurement0[pj][i,j] for j in 1:dim] for i in 1:M_num] end for pj in (length(measurement0)+1):p_num M_tp = [Vector{ComplexF64}(undef, dim) for i in 1:M_num] @@ -331,7 +331,7 @@ function initial_LinearComb!(measurement0, B_all, basis_num, M_num, p_num, rng) measurement0 = [measurement0[i] for i in 1:p_num] end for pj in 1:length(measurement0) - B_all[pj] = deepcopy(measurement0[pj]) #[[measurement0[pj][i,j] for j in 1:basis_num] for i in 1:M_num] + B_all[pj] = measurement0[pj] isa AbstractVector ? deepcopy(measurement0[pj]) : [[measurement0[pj][i,j] for j in 1:dim] for i in 1:M_num] end for pj in (length(measurement0)+1):p_num @@ -353,29 +353,3 @@ function initial_Rotation!(measurement0, s_all, dim, p_num, rng) end end -function initial_LinearComb!(measurement0, B_all, basis_num, M_num, p_num, rng) - if length(measurement0) > p_num - measurement0 = [measurement0[i] for i in 1:p_num] - end - for pj in 1:length(measurement0) - B_all[pj] = deepcopy(measurement0[pj]) #[[measurement0[pj][i,j] for j in 1:basis_num] for i in 1:M_num] - end - - for pj in (length(measurement0)+1):p_num - B_all[pj] = [rand(rng, basis_num) for i in 1:M_num] - bound_LC_coeff!(B_all[pj]) - end -end - -function initial_Rotation!(measurement0, s_all, dim, p_num, rng) - if length(measurement0) > p_num - measurement0 = [measurement0[i] for i in 1:p_num] - end - for pj in 1:length(measurement0) - s_all[pj] = [measurement0[pj][i] for i in 1:dim*dim] - end - - for pj in (length(measurement0)+1):p_num - s_all[pj] = rand(rng, dim*dim) - end -end diff --git a/src/dynamics/kraus/KrausWrapper.jl b/src/dynamics/kraus/KrausWrapper.jl index 5f14516..e64bb17 100644 --- a/src/dynamics/kraus/KrausWrapper.jl +++ b/src/dynamics/kraus/KrausWrapper.jl @@ -1,19 +1,19 @@ function Kraus(opt::StateOpt, K, dK;rng=GLOBAL_RNG) - (;ψ₀) = opt + (;psi) = opt dim = size(K[1], 1) - if ismissing(ψ₀) + if ismissing(psi) r_ini = 2*rand(rng, dim) - ones(dim) r = r_ini ./ norm(r_ini) ϕ = 2pi*rand(rng, dim) - ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] - opt.ψ₀ = ψ₀ + psi = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.psi = psi end K = complex.(K) dK = [complex.(dk) for dk in dK] - ψ₀ = complex(ψ₀) + psi = complex(psi) - Kraus(K, dK, ψ₀) + Kraus(K, dK, psi) end function Kraus(opt::AbstractMopt, ρ₀::AbstractMatrix, K, dK;rng=GLOBAL_RNG, eps=GLOBAL_EPS) @@ -28,19 +28,19 @@ end function Kraus(opt::CompOpt, K, dK;rng=GLOBAL_RNG, eps=GLOBAL_EPS) - (;ψ₀) = opt + (;psi) = opt dim = size(K[1], 1) - if ismissing(ψ₀) + if ismissing(psi) r_ini = 2*rand(rng, dim) - ones(dim) r = r_ini ./ norm(r_ini) ϕ = 2pi*rand(rng, dim) - ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] - opt.ψ₀ = ψ₀ + psi = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.psi = psi end _ini_measurement!(opt, dim, rng; eps=eps) K = complex.(K) dK = [complex.(dk) for dk in dK] - ψ₀ = complex(ψ₀) - Kraus(K, dK, ψ₀) + psi = complex(psi) + Kraus(K, dK, psi) end \ No newline at end of file diff --git a/src/dynamics/lindblad/LindbladDynamics.jl b/src/dynamics/lindblad/LindbladDynamics.jl index e3b058b..e840e9e 100644 --- a/src/dynamics/lindblad/LindbladDynamics.jl +++ b/src/dynamics/lindblad/LindbladDynamics.jl @@ -51,7 +51,7 @@ function liouville_dissip_py(A::Array{T}) where {T<:Complex} end function dissipation( - Γ::Vector{Matrix{T}}, + Γ::AbstractVector, γ::Vector{R}, t::Int = 0, ) where {T<:Complex,R<:Real} @@ -59,7 +59,7 @@ function dissipation( end function dissipation( - Γ::Vector{Matrix{T}}, + Γ::AbstractVector, γ::Vector{Vector{R}}, t::Int = 0, ) where {T<:Complex,R<:Real} @@ -72,7 +72,7 @@ end function liouvillian( H::Matrix{T}, - decay_opt::Vector{Matrix{T}}, + decay_opt::AbstractVector, γ, t = 1, ) where {T<:Complex} @@ -109,14 +109,154 @@ function expL(H, dt) end function expm( + tspan::AbstractVector, + ρ0::AbstractMatrix, H0::AbstractMatrix, dH::AbstractMatrix, - Hc::AbstractVector, - ctrl::AbstractVector, + decay::Union{AbstractVector, Missing}=missing, + Hc::Union{AbstractVector, Missing}=missing, + ctrl::Union{AbstractVector, Missing}=missing, +) + dim = size(ρ0, 1) + tnum = length(tspan) + if ismissing(decay) + decay_opt = [zeros(ComplexF64, dim, dim)] + γ = [0.0] + else + decay_opt = [decay[1] for decay in decay] + γ = [decay[2] for decay in decay] + end + + if ismissing(Hc) + Hc = [zeros(ComplexF64, dim, dim)] + ctrl = [zeros(tnum-1)] + elseif ismissing(ctrl) + ctrl = [zeros(tnum-1)] + else + ctrl_num = length(Hc) + ctrl_length = length(ctrl) + if ctrl_num < ctrl_length + throw(ArgumentError( + "There are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences: too many coefficients sequences" + )) + elseif ctrl_num < ctrl_length + throw(ArgumentError( + "Not enough coefficients sequences: there are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences. The rest of the control sequences are set to be 0." + )) + end + + ratio_num = ceil((length(tspan)-1) / length(ctrl[1])) + if length(tspan) - 1 % length(ctrl[1]) != 0 + tnum = ratio_num * length(ctrl[1]) |> Int + tspan = range(tspan[1], tspan[end], length=tnum+1) + end + end + ctrl_num = length(Hc) + ctrl_interval = ((length(tspan) - 1) / length(ctrl[1])) |> Int + ctrl = [repeat(ctrl[i], 1, ctrl_interval) |> transpose |> vec for i = 1:ctrl_num] + + H = Htot(H0, Hc, ctrl) + dH_L = liouville_commu(dH) + + Δt = tspan[2] - tspan[1] + + ρt_all = [Vector{ComplexF64}(undef, (length(H0))^2) for i = 1:length(tspan)] + ∂ρt_∂x_all = [Vector{ComplexF64}(undef, (length(H0))^2) for i = 1:length(tspan)] + ρt_all[1] = ρ0 |> vec + ∂ρt_∂x_all[1] = ρt_all[1] |> zero + + decay_opt, γ = decay + for t = 2:length(tspan) + exp_L = expL(H[t-1], decay_opt, γ, Δt, t) + ρt_all[t] = exp_L * ρt_all[t-1] + ∂ρt_∂x_all[t] = -im * Δt * dH_L * ρt_all[t] + exp_L * ∂ρt_∂x_all[t-1] + end + ρt_all |> vec2mat, ∂ρt_∂x_all |> vec2mat +end + +function expm( + tspan::AbstractVector, ρ0::AbstractMatrix, + H0::AbstractMatrix, + dH::AbstractVector, + decay::Union{AbstractVector, Missing}=missing, + Hc::Union{AbstractVector, Missing}=missing, + ctrl::Union{AbstractVector, Missing}=missing +) + dim = size(ρ0, 1) + tnum = length(tspan) + if ismissing(decay) + decay_opt = [zeros(ComplexF64, dim, dim)] + γ = [0.0] + else + decay_opt = [decay[1] for decay in decay] + γ = [decay[2] for decay in decay] + end + + if ismissing(Hc) + Hc = [zeros(ComplexF64, dim, dim)] + ctrl0 = [zeros(tnum-1)] + elseif ismissing(ctrl) + ctrl0 = [zeros(tnum-1)] + else + ctrl_num = length(Hc) + ctrl_length = length(ctrl) + if ctrl_num < ctrl_length + throw(ArgumentError( + "There are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences: too many coefficients sequences" + )) + elseif ctrl_num < ctrl_length + throw(ArgumentError( + "Not enough coefficients sequences: there are $ctrl_num control Hamiltonians but $ctrl_length coefficients sequences. The rest of the control sequences are set to be 0." + )) + end + + ratio_num = ceil((length(tspan)-1) / length(ctrl[1])) + if length(tspan) - 1 % length(ctrl[1]) != 0 + tnum = ratio_num * length(ctrl[1]) |> Int + tspan = range(tspan[1], tspan[end], length=tnum+1) + end + ctrl0 = ctrl + end + para_num = length(dH) + ctrl_num = length(Hc) + ctrl_interval = ((length(tspan) - 1) / length(ctrl0[1])) |> Int + ctrl = [repeat(ctrl0[i], 1, ctrl_interval) |> transpose |> vec for i = 1:ctrl_num] + + H = Htot(H0, Hc, ctrl) + dH_L = [liouville_commu(dH[i]) for i = 1:para_num] + + Δt = tspan[2] - tspan[1] + + ρt_all = [Vector{ComplexF64}(undef, (length(H0))^2) for i = 1:length(tspan)] + ∂ρt_∂x_all = [ + [Vector{ComplexF64}(undef, (length(H0))^2) for j = 1:para_num] for + i = 1:length(tspan) + ] + ρt_all[1] = ρ0 |> vec + for pj = 1:para_num + ∂ρt_∂x_all[1][pj] = ρt_all[1] |> zero + end + + for t = 2:length(tspan) + exp_L = expL(H[t-1], decay_opt, γ, Δt, t) + ρt_all[t] = exp_L * ρt_all[t-1] + for pj = 1:para_num + ∂ρt_∂x_all[t][pj] = -im * Δt * dH_L[pj] * ρt_all[t] + exp_L * ∂ρt_∂x_all[t-1][pj] + end + end + ρt_all |> vec2mat, ∂ρt_∂x_all |> vec2mat +end + +function expm_py( tspan, + ρ0::AbstractMatrix, + H0::AbstractMatrix, + dH::AbstractMatrix, + Hc::AbstractVector, decay_opt::AbstractVector, γ, + ctrl::AbstractVector, ) where {T<:Complex,R<:Real} ctrl_num = length(Hc) @@ -141,15 +281,15 @@ function expm( ρt_all |> vec2mat, ∂ρt_∂x_all |> vec2mat end -function expm( +function expm_py( + tspan, + ρ0::AbstractMatrix, H0::AbstractMatrix, dH::AbstractVector, - Hc::AbstractVector, - ctrl::AbstractVector, - ρ0::AbstractMatrix, - tspan, decay_opt::AbstractVector, γ, + Hc::AbstractVector, + ctrl::AbstractVector, ) para_num = length(dH) diff --git a/src/dynamics/lindblad/LindbladWrapper.jl b/src/dynamics/lindblad/LindbladWrapper.jl index 87bc0d9..4a27381 100644 --- a/src/dynamics/lindblad/LindbladWrapper.jl +++ b/src/dynamics/lindblad/LindbladWrapper.jl @@ -1,5 +1,5 @@ # wrapper for Lindblad dynamics with ControlOpt -function Lindblad(opt::ControlOpt, tspan, ρ₀, H0, dH, Hc=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) +function Lindblad(opt::ControlOpt, tspan, ρ₀, H0, dH, Hc; decay=missing, rng=GLOBAL_RNG, eps=GLOBAL_EPS) (;ctrl) = opt dim = size(ρ₀, 1) if ismissing(dH) @@ -20,7 +20,7 @@ function Lindblad(opt::ControlOpt, tspan, ρ₀, H0, dH, Hc=missing, decay=missi Hc = [zeros(ComplexF64, dim, dim)] ctrl = [zeros(tnum-1)] elseif ismissing(ctrl) - opt.ctrl = [zeros(tnum-1) for _ in 1:ctrl_num] + ctrl = [zeros(tnum-1) for _ in 1:ctrl_num] opt.ctrl = ctrl else ctrl_length = length(ctrl) @@ -47,15 +47,18 @@ function Lindblad(opt::ControlOpt, tspan, ρ₀, H0, dH, Hc=missing, decay=missi Lindblad(H0, dH, Hc, ctrl, ρ₀, tspan, decay_opt, γ) end -function Lindblad(opt::StateOpt, tspan, H0, dH, Hc=missing, ctrl=missing, decay=missing;rng=GLOBAL_RNG) - (;ψ₀) = opt +Lindblad(opt::ControlOpt, tspan, ρ₀, H0, dH, Hc, decay; rng=GLOBAL_RNG, eps=GLOBAL_EPS) = + Lindblad(opt, tspan, ρ₀, H0, dH, Hc; decay=decay, rng=rng, eps=eps) + +function Lindblad(opt::StateOpt, tspan, H0, dH; Hc=missing, ctrl=missing, decay=missing, rng=GLOBAL_RNG, eps=GLOBAL_EPS) + (;psi) = opt dim = H0 isa AbstractVector ? size(H0[1], 1) : size(H0, 1) - if ismissing(ψ₀) + if ismissing(psi) r_ini = 2*rand(rng, dim) - ones(dim) r = r_ini ./ norm(r_ini) ϕ = 2pi*rand(rng, dim) - ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] - opt.ψ₀ = ψ₀ + psi = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.psi = psi end if ismissing(dH) @@ -116,19 +119,22 @@ function Lindblad(opt::StateOpt, tspan, H0, dH, Hc=missing, ctrl=missing, decay= end dH = complex.(dH) - ψ₀ = complex(ψ₀) + psi = complex(psi) if all(iszero.(γ)) # if any non-zero decay rate - return Lindblad(H0, dH, ψ₀, tspan) + return Lindblad(H0, dH, psi, tspan) else - return Lindblad(H0, dH, ψ₀, tspan, decay_opt, γ) + return Lindblad(H0, dH, psi, tspan, decay_opt, γ) end end +Lindblad(opt::StateOpt, tspan, H0, dH, Hc, ctrl, decay; rng=GLOBAL_RNG, eps=GLOBAL_EPS) = + Lindblad(opt, tspan, H0, dH; Hc=Hc, ctrl=ctrl, decay=decay, rng=rng, eps=eps) + function _ini_measurement!(opt::Mopt_Projection, dim::Int, rng; eps=GLOBAL_EPS) - (;C) = opt - ## initialize the Mopt target C - if ismissing(C) + (; M) = opt + ## initialize the Mopt target M + if ismissing(M) M = [ComplexF64[] for _ in 1:dim] for i in 1:dim r_ini = 2*rand(rng, dim) - ones(dim) @@ -136,10 +142,10 @@ function _ini_measurement!(opt::Mopt_Projection, dim::Int, rng; eps=GLOBAL_EPS) ϕ = 2pi*rand(rng, dim) M[i] = [r*exp(im*ϕ) for (r,ϕ) in zip(r,ϕ)] end - C = gramschmidt(M) + M = gramschmidt(M) end - opt.C = complex.(C) + opt.M = complex.(M) end function _ini_measurement!(opt::Mopt_LinearComb, dim::Int, rng; eps=GLOBAL_EPS) @@ -163,7 +169,6 @@ function _ini_measurement!(opt::Mopt_LinearComb, dim::Int, rng; eps=GLOBAL_EPS) end end - function _ini_measurement!(opt::Mopt_Rotation, dim::Int, rng; eps=GLOBAL_EPS) (;s, POVM_basis, Lambda) = opt if ismissing(POVM_basis) @@ -185,7 +190,7 @@ function _ini_measurement!(opt::Mopt_Rotation, dim::Int, rng; eps=GLOBAL_EPS) end end -function Lindblad(opt::AbstractMopt, tspan, ρ₀, H0, dH, Hc=missing, ctrl=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) +function Lindblad(opt::AbstractMopt, tspan, ρ₀, H0, dH; Hc=missing, ctrl=missing, decay=missing, rng=GLOBAL_RNG, eps=GLOBAL_EPS) dim = size(ρ₀, 1) _ini_measurement!(opt, dim, rng; eps=eps) @@ -256,30 +261,33 @@ function Lindblad(opt::AbstractMopt, tspan, ρ₀, H0, dH, Hc=missing, ctrl=miss end end +Lindblad(opt::AbstractMopt, tspan, ρ₀, H0, dH, Hc, ctrl, decay; rng=GLOBAL_RNG, eps=GLOBAL_EPS) = + Lindblad(opt, tspan, ρ₀, H0, dH; Hc=Hc, ctrl=ctrl, decay=decay, rng=rng, eps=eps) + function _ini_measurement!(opt::CompOpt, dim::Int, rng; eps=GLOBAL_EPS) - (;C) = opt - ## initialize the Mopt target C + (; M) = opt + ## initialize the Mopt target M M = [ComplexF64[] for _ in 1:dim] - if ismissing(C) + if ismissing(M) for i in 1:dim 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 - opt.C = gramschmidt(M) + opt.M = gramschmidt(M) end end -function Lindblad(opt::StateControlOpt, tspan, H0, dH, Hc=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) - (;ψ₀, ctrl) = opt +function Lindblad(opt::StateControlOpt, tspan, H0, dH, Hc; decay=missing, rng=GLOBAL_RNG, eps=GLOBAL_EPS) + (;psi, ctrl) = opt dim = H0 isa AbstractVector ? size(H0[1], 1) : size(H0, 1) - if ismissing(ψ₀) + if ismissing(psi) r_ini = 2*rand(rng, dim) - ones(dim) r = r_ini ./ norm(r_ini) ϕ = 2pi*rand(rng, dim) - ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] - opt.ψ₀ = ψ₀ + psi = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.psi = psi end if ismissing(dH) dH = [zeros(ComplexF64, dim, dim)] @@ -321,11 +329,14 @@ function Lindblad(opt::StateControlOpt, tspan, H0, dH, Hc=missing, decay=missing end H0 = complex(H0) dH = complex.(dH) - ψ₀ = complex(ψ₀) - Lindblad(H0, dH, Hc, ctrl, ψ₀, tspan, decay_opt, γ) + psi = complex(psi) + Lindblad(H0, dH, Hc, ctrl, psi, tspan, decay_opt, γ) end -function Lindblad(opt::ControlMeasurementOpt, tspan, ρ₀, H0, dH, Hc=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) +Lindblad(opt::StateControlOpt, tspan, H0, dH, Hc, decay; rng=GLOBAL_RNG, eps=GLOBAL_EPS) = + Lindblad(opt, tspan, H0, dH, Hc; decay=decay, rng=rng, eps=eps) + +function Lindblad(opt::ControlMeasurementOpt, tspan, ρ₀, H0, dH, Hc; decay=missing, rng=GLOBAL_RNG, eps=GLOBAL_EPS) (;ctrl) = opt dim = size(ρ₀, 1) _ini_measurement!(opt, dim, rng; eps=eps) @@ -374,16 +385,19 @@ function Lindblad(opt::ControlMeasurementOpt, tspan, ρ₀, H0, dH, Hc=missing, Lindblad(H0, dH, Hc, ctrl, ρ₀, tspan, decay_opt, γ) end -function Lindblad(opt::StateMeasurementOpt, tspan, H0, dH, Hc=missing, ctrl=missing, decay=missing;rng=GLOBAL_RNG) - (;ψ₀) = opt +Lindblad(opt::ControlMeasurementOpt, tspan, ρ₀, H0, dH, Hc, decay; rng=GLOBAL_RNG, eps=GLOBAL_EPS) = + Lindblad(opt, tspan, ρ₀, H0, dH, Hc; decay=decay, rng=rng, eps=eps) + +function Lindblad(opt::StateMeasurementOpt, tspan, H0, dH; Hc=missing, ctrl=missing, decay=missing, rng=GLOBAL_RNG) + (;psi) = opt dim = H0 isa AbstractVector ? size(H0[1], 1) : size(H0, 1) _ini_measurement!(opt, dim, rng; eps=eps) - if ismissing(ψ₀) + if ismissing(psi) r_ini = 2*rand(rng, dim) - ones(dim) r = r_ini ./ norm(r_ini) ϕ = 2pi*rand(rng, dim) - ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] - opt.ψ₀ = ψ₀ + psi = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.psi = psi end if ismissing(dH) @@ -444,26 +458,29 @@ function Lindblad(opt::StateMeasurementOpt, tspan, H0, dH, Hc=missing, ctrl=miss end dH = complex.(dH) - ψ₀ = complex(ψ₀) + psi = complex(psi) if all(iszero.(γ)) # if any non-zero decay rate - return Lindblad(H0, dH, ψ₀, tspan) + return Lindblad(H0, dH, psi, tspan) else - return Lindblad(H0, dH, ψ₀, tspan, decay_opt, γ) + return Lindblad(H0, dH, psi, tspan, decay_opt, γ) end end -function Lindblad(opt::StateControlMeasurementOpt, tspan, H0, dH, Hc=missing, decay=missing; rng=GLOBAL_RNG, eps=GLOBAL_EPS) - (;ctrl, ψ₀) = opt +Lindblad(opt::StateMeasurementOpt, tspan, H0, dH, Hc, ctrl, decay; rng=GLOBAL_RNG, eps=GLOBAL_EPS) = + Lindblad(opt, tspan, H0, dH; Hc=Hc, ctrl=ctrl, decay=decay, rng=rng, eps=eps) + +function Lindblad(opt::StateControlMeasurementOpt, tspan, H0, dH, Hc; decay=missing, rng=GLOBAL_RNG, eps=GLOBAL_EPS) + (;ctrl, psi) = opt dim = H0 isa AbstractVector ? size(H0[1], 1) : size(H0, 1) _ini_measurement!(opt, dim, rng; eps=eps) - if ismissing(ψ₀) + if ismissing(psi) r_ini = 2*rand(rng, dim) - ones(dim) r = r_ini ./ norm(r_ini) ϕ = 2pi*rand(rng, dim) - ψ₀ = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] - opt.ψ₀ = ψ₀ + psi = [r*exp(im*ϕ) for (r, ϕ) in zip(r, ϕ)] + opt.psi = psi end if ismissing(dH) @@ -506,7 +523,10 @@ function Lindblad(opt::StateControlMeasurementOpt, tspan, H0, dH, Hc=missing, de end H0 = complex(H0) dH = complex.(dH) - ψ₀ = complex(ψ₀) + psi = complex(psi) - Lindblad(H0, dH, Hc, ctrl, ψ₀, tspan, decay_opt, γ) -end \ No newline at end of file + Lindblad(H0, dH, Hc, ctrl, psi, tspan, decay_opt, γ) +end + +Lindblad(opt::StateControlMeasurementOpt, tspan, H0, dH, Hc, decay; rng=GLOBAL_RNG, eps=GLOBAL_EPS) = + Lindblad(opt, tspan, H0, dH, Hc; decay=decay, rng=rng, eps=eps) diff --git a/src/objective/AsymptoticBound/AsymptoticBound.jl b/src/objective/AsymptoticBound/AsymptoticBound.jl index 040aed4..2962885 100644 --- a/src/objective/AsymptoticBound/AsymptoticBound.jl +++ b/src/objective/AsymptoticBound/AsymptoticBound.jl @@ -4,232 +4,236 @@ abstract type SLD <: AbstractLDtype end abstract type RLD <: AbstractLDtype end abstract type LLD <: AbstractLDtype end -Base.@kwdef struct QFIM_Obj{P,D} <: AbstractObj +Base.@kwdef struct QFIM_obj{P,D} <: AbstractObj W::Union{AbstractMatrix, Missing} eps::Number = GLOBAL_EPS end -Base.@kwdef struct CFIM_Obj{P} <: AbstractObj +Base.@kwdef struct CFIM_obj{P} <: AbstractObj M::Union{AbstractVecOrMat, Missing} W::Union{AbstractMatrix, Missing} eps::Number = GLOBAL_EPS end -Base.@kwdef struct HCRB_Obj{P} <: AbstractObj +Base.@kwdef struct HCRB_obj{P} <: AbstractObj W::Union{AbstractMatrix, Missing} eps::Number = GLOBAL_EPS end -QFIM_Obj(;W=missing, eps=GLOBAL_EPS, para_type::Symbol=:single_para, LD_type::Symbol=:SLD) = QFIM_Obj{eval.([para_type, LD_type])...}(W, eps) -CFIM_Obj(;M=missing, W=missing, eps=GLOBAL_EPS, para_type::Symbol=:single_para) = CFIM_Obj{eval(para_type)}(M, W, eps) -HCRB_Obj(;W=missing, eps=GLOBAL_EPS, para_type::Symbol=:single_para) = HCRB_Obj{eval(para_type)}(W, eps) +QFIM_obj(;W=missing, eps=GLOBAL_EPS, para_type::Symbol=:single_para, LD_type::Symbol=:SLD) = QFIM_obj{eval.([para_type, LD_type])...}(W, eps) +CFIM_obj(;M=missing, W=missing, eps=GLOBAL_EPS, para_type::Symbol=:single_para) = CFIM_obj{eval(para_type)}(M, W, eps) +HCRB_obj(;W=missing, eps=GLOBAL_EPS, para_type::Symbol=:single_para) = HCRB_obj{eval(para_type)}(W, eps) -QFIM_Obj(W::AbstractMatrix, eps::Number, para_type::String, LD_type::String) = QFIM_Obj(W, eps, Symbol.([para_type, LD_type])...) -CFIM_Obj(M::AbstractVecOrMat, W::AbstractMatrix, eps::Number, para_type::String) = CFIM_Obj(M, W, eps, Symbol(para_type)) -HCRB_Obj(W::AbstractMatrix, eps::Number, para_type::String) = HCRB_Obj(W, eps, Symbol(para_type)) +QFIM_obj(W, eps, para_type::Symbol, LD_type::Symbol) = QFIM_obj{eval.([para_type, LD_type])...}(W, eps) +CFIM_obj(M, W, eps, para_type::Symbol) = CFIM_obj{eval(para_type)}(M, W, eps) +HCRB_obj(W, eps, para_type::Symbol) = HCRB_obj{eval(para_type)}(W, eps) -obj_type(::QFIM_Obj) = :QFIM -obj_type(::CFIM_Obj) = :CFIM -obj_type(::HCRB_Obj) = :HCRB +QFIM_obj(W::AbstractMatrix, eps::Number, para_type::String, LD_type::String) = QFIM_obj(W, eps, Symbol.([para_type, LD_type])...) +CFIM_obj(M::AbstractVecOrMat, W::AbstractMatrix, eps::Number, para_type::String) = CFIM_obj(M, W, eps, Symbol(para_type)) +HCRB_obj(W::AbstractMatrix, eps::Number, para_type::String) = HCRB_obj(W, eps, Symbol(para_type)) -para_type(::QFIM_Obj{single_para,D}) where {D} = :single_para -para_type(::QFIM_Obj{multi_para,D}) where {D} = :multi_para -para_type(::CFIM_Obj{single_para}) = :single_para -para_type(::CFIM_Obj{multi_para}) = :multi_para -para_type(::HCRB_Obj{single_para}) = :single_para -para_type(::HCRB_Obj{multi_para}) = :multi_para +obj_type(::QFIM_obj) = :QFIM +obj_type(::CFIM_obj) = :CFIM +obj_type(::HCRB_obj) = :HCRB -LD_type(::QFIM_Obj{P,SLD}) where {P} = :SLD -LD_type(::QFIM_Obj{P,RLD}) where {P} = :RLD -LD_type(::QFIM_Obj{P,LLD}) where {P} = :LLD +para_type(::QFIM_obj{single_para,D}) where {D} = :single_para +para_type(::QFIM_obj{multi_para,D}) where {D} = :multi_para +para_type(::CFIM_obj{single_para}) = :single_para +para_type(::CFIM_obj{multi_para}) = :multi_para +para_type(::HCRB_obj{single_para}) = :single_para +para_type(::HCRB_obj{multi_para}) = :multi_para -QFIM_Obj(opt::CFIM_Obj{P}) where P = QFIM_Obj{P, SLD}(opt.W, opt.eps) -QFIM_Obj(opt::CFIM_Obj{P}, LDtype::Symbol) where P = QFIM_Obj{P, eval(LDtype)}(opt.W, opt.eps) +LD_type(::QFIM_obj{P,SLD}) where {P} = :SLD +LD_type(::QFIM_obj{P,RLD}) where {P} = :RLD +LD_type(::QFIM_obj{P,LLD}) where {P} = :LLD + +QFIM_obj(opt::CFIM_obj{P}) where P = QFIM_obj{P, SLD}(opt.W, opt.eps) +QFIM_obj(opt::CFIM_obj{P}, LDtype::Symbol) where P = QFIM_obj{P, eval(LDtype)}(opt.W, opt.eps) const obj_idx = Dict( - :QFIM => QFIM_Obj, - :CFIM => CFIM_Obj, - :HCRB => HCRB_Obj + :QFIM => QFIM_obj, + :CFIM => CFIM_obj, + :HCRB => HCRB_obj ) -function set_M(obj::CFIM_Obj{P}, M::AbstractVector) where P - CFIM_Obj{P}(M, obj.W, obj.eps) +function set_M(obj::CFIM_obj{P}, M::AbstractVector) where P + CFIM_obj{P}(M, obj.W, obj.eps) end -function objective(obj::QFIM_Obj{single_para,SLD}, dynamics::Lindblad) +function objective(obj::QFIM_obj{single_para,SLD}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = W[1] * QFIM_SLD(ρ, dρ[1]; eps = eps) return f, f end -function objective(obj::QFIM_Obj{single_para,SLD}, ρ, dρ) +function objective(obj::QFIM_obj{single_para,SLD}, ρ, dρ) (; W, eps) = obj f = W[1] * QFIM_SLD(ρ, dρ[1]; eps = eps) return f, f end -function objective(obj::QFIM_Obj{multi_para,SLD}, ρ, dρ) +function objective(obj::QFIM_obj{multi_para,SLD}, ρ, dρ) (; W, eps) = obj f = tr(W * pinv(QFIM_SLD(ρ, dρ; eps = eps))) return f, 1.0 / f end -function objective(obj::QFIM_Obj{single_para,RLD}, ρ, dρ) +function objective(obj::QFIM_obj{single_para,RLD}, ρ, dρ) (; W, eps) = obj f = W[1] * QFIM_RLD(ρ, dρ[1]; eps = eps) return f, f end -function objective(obj::QFIM_Obj{multi_para,RLD}, ρ, dρ) +function objective(obj::QFIM_obj{multi_para,RLD}, ρ, dρ) (; W, eps) = obj f = tr(W * pinv(QFIM_RLD(ρ, dρ; eps = eps))) |> real return f, 1.0 / f end -function objective(obj::QFIM_Obj{single_para,LLD}, ρ, dρ) +function objective(obj::QFIM_obj{single_para,LLD}, ρ, dρ) (; W, eps) = obj f = W[1] * QFIM_LLD(ρ, dρ[1]; eps = eps) return f, f end -function objective(obj::QFIM_Obj{multi_para,LLD}, ρ, dρ) +function objective(obj::QFIM_obj{multi_para,LLD}, ρ, dρ) (; W, eps) = obj f = tr(W * pinv(QFIM_LLD(ρ, dρ; eps = eps))) |> real return f, 1.0 / f end -function objective(obj::QFIM_Obj{multi_para,SLD}, dynamics::Lindblad) +function objective(obj::QFIM_obj{multi_para,SLD}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = tr(W * pinv(QFIM_SLD(ρ, dρ; eps = eps))) return f, 1.0 / f end -function objective(obj::QFIM_Obj{single_para,RLD}, dynamics::Lindblad) +function objective(obj::QFIM_obj{single_para,RLD}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = W[1] * QFIM_RLD(ρ, dρ[1]; eps = eps) return f, f end -function objective(obj::QFIM_Obj{multi_para,RLD}, dynamics::Lindblad) +function objective(obj::QFIM_obj{multi_para,RLD}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = tr(W * pinv(QFIM_RLD(ρ, dρ; eps = eps))) |> real return f, 1.0 / f end -function objective(obj::QFIM_Obj{single_para,LLD}, dynamics::Lindblad) +function objective(obj::QFIM_obj{single_para,LLD}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = W[1] * QFIM_LLD(ρ, dρ[1]; eps = eps) return f, f end -function objective(obj::QFIM_Obj{multi_para,LLD}, dynamics::Lindblad) +function objective(obj::QFIM_obj{multi_para,LLD}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = tr(W * pinv(QFIM_LLD(ρ, dρ; eps = eps))) |> real return f, 1.0 / f end -function objective(obj::QFIM_Obj{single_para,SLD}, dynamics::Kraus) +function objective(obj::QFIM_obj{single_para,SLD}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = W[1] * QFIM_SLD(ρ, dρ[1]; eps = eps) return f, f end -function objective(obj::QFIM_Obj{multi_para,SLD}, dynamics::Kraus) +function objective(obj::QFIM_obj{multi_para,SLD}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = tr(W * pinv(QFIM_SLD(ρ, dρ; eps = eps))) return f, 1.0 / f end -function objective(obj::QFIM_Obj{single_para,RLD}, dynamics::Kraus) +function objective(obj::QFIM_obj{single_para,RLD}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = W[1] * QFIM_RLD(ρ, dρ[1]; eps = eps) return f, f end -function objective(obj::QFIM_Obj{multi_para,RLD}, dynamics::Kraus) +function objective(obj::QFIM_obj{multi_para,RLD}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = tr(W * pinv(QFIM_RLD(ρ, dρ; eps = eps))) |> real return f, 1.0 / f end -function objective(obj::QFIM_Obj{single_para,LLD}, dynamics::Kraus) +function objective(obj::QFIM_obj{single_para,LLD}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = W[1] * QFIM_LLD(ρ, dρ[1]; eps = eps) return f, f end -function objective(obj::QFIM_Obj{multi_para,LLD}, dynamics::Kraus) +function objective(obj::QFIM_obj{multi_para,LLD}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = tr(W * pinv(QFIM_LLD(ρ, dρ; eps = eps))) |> real return f, 1.0 / f end -function objective(obj::CFIM_Obj{single_para}, ρ, dρ) +function objective(obj::CFIM_obj{single_para}, ρ, dρ) (; M, W, eps) = obj f = W[1] * CFIM(ρ, dρ[1], M; eps = eps) return f, f end -function objective(obj::CFIM_Obj{multi_para}, ρ, dρ) +function objective(obj::CFIM_obj{multi_para}, ρ, dρ) (; M, W, eps) = obj f = tr(W * pinv(CFIM(ρ, dρ, M; eps = eps))) return f, 1.0 / f end -function objective(obj::CFIM_Obj{single_para}, dynamics::Lindblad) +function objective(obj::CFIM_obj{single_para}, dynamics::Lindblad) (; M, W, eps) = obj ρ, dρ = evolve(dynamics) f = W[1] * CFIM(ρ, dρ[1], M; eps = eps) return f, f end -function objective(obj::CFIM_Obj{multi_para}, dynamics::Lindblad) +function objective(obj::CFIM_obj{multi_para}, dynamics::Lindblad) (; M, W, eps) = obj ρ, dρ = evolve(dynamics) f = tr(W * pinv(CFIM(ρ, dρ, M; eps = eps))) return f, 1.0 / f end -function objective(obj::CFIM_Obj{single_para}, dynamics::Kraus) +function objective(obj::CFIM_obj{single_para}, dynamics::Kraus) (; M, W, eps) = obj ρ, dρ = evolve(dynamics) f = W[1] * CFIM(ρ, dρ[1], M; eps = eps) return f, f end -function objective(obj::CFIM_Obj{multi_para}, dynamics::Kraus) +function objective(obj::CFIM_obj{multi_para}, dynamics::Kraus) (; M, W, eps) = obj ρ, dρ = evolve(dynamics) f = tr(W * pinv(CFIM(ρ, dρ, M; eps = eps))) return f, 1.0 / f end -function objective(obj::HCRB_Obj{multi_para}, ρ, dρ) +function objective(obj::HCRB_obj{multi_para}, ρ, dρ) (; W, eps) = obj f = Holevo_bound_obj(ρ, dρ, W; eps = eps) return f, 1.0 / f end -function objective(obj::HCRB_Obj{multi_para}, dynamics::Lindblad) +function objective(obj::HCRB_obj{multi_para}, dynamics::Lindblad) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = Holevo_bound_obj(ρ, dρ, W; eps = eps) return f, 1.0 / f end -function objective(obj::HCRB_Obj{multi_para}, dynamics::Kraus) +function objective(obj::HCRB_obj{multi_para}, dynamics::Kraus) (; W, eps) = obj ρ, dρ = evolve(dynamics) f = Holevo_bound_obj(ρ, dρ, W; eps = eps) @@ -237,7 +241,7 @@ function objective(obj::HCRB_Obj{multi_para}, dynamics::Kraus) end #### objective function for linear combination in Mopt #### -function objective(opt::Mopt_LinearComb, obj::CFIM_Obj{single_para}, dynamics::Lindblad) +function objective(opt::Mopt_LinearComb, obj::CFIM_obj{single_para}, dynamics::Lindblad) (; W, eps) = obj M = [sum([opt.B[i][j]*opt.POVM_basis[j] for j in 1:length(opt.POVM_basis)]) for i in 1:opt.M_num] ρ, dρ = evolve(dynamics) @@ -245,7 +249,7 @@ function objective(opt::Mopt_LinearComb, obj::CFIM_Obj{single_para}, dynamics::L return f, f end -function objective(opt::Mopt_LinearComb, obj::CFIM_Obj{multi_para}, dynamics::Lindblad) +function objective(opt::Mopt_LinearComb, obj::CFIM_obj{multi_para}, dynamics::Lindblad) (; W, eps) = obj M = [sum([opt.B[i][j]*opt.POVM_basis[j] for j in 1:length(opt.POVM_basis)]) for i in 1:opt.M_num] ρ, dρ = evolve(dynamics) @@ -253,7 +257,7 @@ function objective(opt::Mopt_LinearComb, obj::CFIM_Obj{multi_para}, dynamics::Li return f, 1.0 / f end -function objective(opt::Mopt_LinearComb, obj::CFIM_Obj{single_para}, dynamics::Kraus) +function objective(opt::Mopt_LinearComb, obj::CFIM_obj{single_para}, dynamics::Kraus) (; W, eps) = obj M = [sum([opt.B[i][j]*opt.POVM_basis[j] for j in 1:length(opt.POVM_basis)]) for i in 1:opt.M_num] ρ, dρ = evolve(dynamics) @@ -261,7 +265,7 @@ function objective(opt::Mopt_LinearComb, obj::CFIM_Obj{single_para}, dynamics::K return f, f end -function objective(opt::Mopt_LinearComb, obj::CFIM_Obj{multi_para}, dynamics::Kraus) +function objective(opt::Mopt_LinearComb, obj::CFIM_obj{multi_para}, dynamics::Kraus) (; W, eps) = obj M = [sum([opt.B[i][j]*opt.POVM_basis[j] for j in 1:length(opt.POVM_basis)]) for i in 1:opt.M_num] ρ, dρ = evolve(dynamics) @@ -270,7 +274,7 @@ function objective(opt::Mopt_LinearComb, obj::CFIM_Obj{multi_para}, dynamics::Kr end #### objective function for rotation in Mopt #### -function objective(opt::Mopt_Rotation, obj::CFIM_Obj{single_para}, dynamics::Lindblad) +function objective(opt::Mopt_Rotation, obj::CFIM_obj{single_para}, dynamics::Lindblad) (; W, eps) = obj U = rotation_matrix(opt.s, opt.Lambda) M = [U*opt.POVM_basis[i]*U' for i in 1:length(opt.POVM_basis)] @@ -279,7 +283,7 @@ function objective(opt::Mopt_Rotation, obj::CFIM_Obj{single_para}, dynamics::Lin return f, f end -function objective(opt::Mopt_Rotation, obj::CFIM_Obj{multi_para}, dynamics::Lindblad) +function objective(opt::Mopt_Rotation, obj::CFIM_obj{multi_para}, dynamics::Lindblad) (; W, eps) = obj U = rotation_matrix(opt.s, opt.Lambda) M = [U*opt.POVM_basis[i]*U' for i in 1:length(opt.POVM_basis)] @@ -288,7 +292,7 @@ function objective(opt::Mopt_Rotation, obj::CFIM_Obj{multi_para}, dynamics::Lind return f, 1.0 / f end -function objective(opt::Mopt_Rotation, obj::CFIM_Obj{single_para}, dynamics::Kraus) +function objective(opt::Mopt_Rotation, obj::CFIM_obj{single_para}, dynamics::Kraus) (; W, eps) = obj U = rotation_matrix(opt.s, opt.Lambda) M = [U*opt.POVM_basis[i]*U' for i in 1:length(opt.POVM_basis)] @@ -297,7 +301,7 @@ function objective(opt::Mopt_Rotation, obj::CFIM_Obj{single_para}, dynamics::Kra return f, f end -function objective(opt::Mopt_Rotation, obj::CFIM_Obj{multi_para}, dynamics::Kraus) +function objective(opt::Mopt_Rotation, obj::CFIM_obj{multi_para}, dynamics::Kraus) (; W, eps) = obj U = rotation_matrix(opt.s, opt.Lambda) M = [U*opt.POVM_basis[i]*U' for i in 1:length(opt.POVM_basis)] diff --git a/src/objective/AsymptoticBound/AsymptoticBoundWrapper.jl b/src/objective/AsymptoticBound/AsymptoticBoundWrapper.jl index c16535e..5950a6b 100644 --- a/src/objective/AsymptoticBound/AsymptoticBoundWrapper.jl +++ b/src/objective/AsymptoticBound/AsymptoticBoundWrapper.jl @@ -1,4 +1,4 @@ -function Objective(dynamics::AbstractDynamics, obj::QFIM_Obj{P,D}) where {P,D} +function Objective(dynamics::AbstractDynamics, obj::QFIM_obj{P,D}) where {P,D} (;W, eps) = obj if ismissing(W) W = I(get_para(dynamics.data))|>Matrix @@ -7,10 +7,10 @@ function Objective(dynamics::AbstractDynamics, obj::QFIM_Obj{P,D}) where {P,D} d = LD_type(obj) |> eval p = para_type(dynamics.data) |> eval - return QFIM_Obj{p,d}(W,eps) + return QFIM_obj{p,d}(W,eps) end -function Objective(dynamics::AbstractDynamics, obj::CFIM_Obj{P}) where {P} +function Objective(dynamics::AbstractDynamics, obj::CFIM_obj{P}) where {P} (;W, M, eps) = obj if ismissing(W) W = I(get_para(dynamics.data))|>Matrix @@ -22,10 +22,10 @@ function Objective(dynamics::AbstractDynamics, obj::CFIM_Obj{P}) where {P} p = para_type(dynamics.data) |> eval - return CFIM_Obj{p}(M,W,eps) + return CFIM_obj{p}(M,W,eps) end -function Objective(dynamics::AbstractDynamics, obj::HCRB_Obj{P}) where {P} +function Objective(dynamics::AbstractDynamics, obj::HCRB_obj{P}) where {P} (;W, eps) = obj if ismissing(W) W = I(get_para(dynamics.data))|>Matrix @@ -33,7 +33,7 @@ function Objective(dynamics::AbstractDynamics, obj::HCRB_Obj{P}) where {P} p = para_type(dynamics.data) |> eval - return HCRB_Obj{p}(W,eps) + return HCRB_obj{p}(W, eps) end diff --git a/src/objective/AsymptoticBound/CramerRao.jl b/src/objective/AsymptoticBound/CramerRao.jl index d1b87e0..8c3c968 100644 --- a/src/objective/AsymptoticBound/CramerRao.jl +++ b/src/objective/AsymptoticBound/CramerRao.jl @@ -2,14 +2,36 @@ using Zygote: @adjoint const σ_x = [0.0 1.0; 1.0 0.0im] const σ_y = [0.0 -1.0im; 1.0im 0.0] const σ_z = [1.0 0.0im; 0.0 -1.0] -const eps_default = 1e-8 ############## logarrithmic derivative ############### +@doc raw""" + + SLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; rep = "original", eps = GLOBAL_EPS) Complex} + +Calculate the symmetric logarrithmic derivatives (SLDs).The SLD operator $L_a$ is defined as``\partial_{a}\rho=\frac{1}{2}(\rho L_{a}+L_{a}\rho)``, where ρ is the parameterized density matrix. + +- `ρ`: Density matrix. +- `dρ`: Derivatives of the density matrix with respect to the unknown parameters to be estimated. For example, drho[1] is the derivative vector with respect to the first parameter. +- `rep`: Representation of the SLD operator. Options can be: + - "original" (default) -- The SLD matrix will be written in terms of the same basis as the input density matrix (ρ). + - "eigen" -- The SLD matrix will be written in terms of the eigenbasis of the input ρ. +- `eps`: Machine epsilon +""" +function SLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; rep = "original", eps = GLOBAL_EPS) where {T<:Complex} + (x -> SLD(ρ, x; rep=rep, eps = eps)).(dρ) +end + +""" + + SLD(ρ::Matrix{T},dρ::Matrix{T};rep = "original",eps = GLOBAL_EPS,) where {T<:Complex} + +When applied to the case of single parameter. +""" function SLD( ρ::Matrix{T}, dρ::Matrix{T}; - eps = eps_default, rep = "original", + eps = GLOBAL_EPS, ) where {T<:Complex} dim = size(ρ)[1] @@ -20,7 +42,7 @@ function SLD( SLD_eig = zeros(T, dim, dim) for fi = 1:dim for fj = 1:dim - if abs(val[fi] + val[fj]) > eps + if abs(val[fi] + val[fj]) > eps SLD_eig[fi, fj] = 2 * (vec[:, fi]' * dρ * vec[:, fj]) / (val[fi] + val[fj]) end end @@ -31,33 +53,35 @@ function SLD( SLD = vec * (SLD_eig * vec') elseif rep == "eigen" SLD = SLD_eig + else + throw(ArgumentError("the rep should be chosen between")) end SLD end -@adjoint function SLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = eps_default) where {T<:Complex} +@adjoint function SLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = GLOBAL_EPS) where {T<:Complex} L = SLD(ρ, dρ; eps = eps) SLD_pullback = L̄ -> (Ḡ -> (-Ḡ * L - L * Ḡ, 2 * Ḡ))(SLD((ρ) |> Array, L̄ / 2)) L, SLD_pullback end -function SLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = eps_default) where {T<:Complex} +function SLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = GLOBAL_EPS) where {T<:Complex} (x -> SLD(ρ, x; eps = eps)).(dρ) end -function SLD_liouville(ρ::Matrix{T}, ∂ρ_∂x::Matrix{T}; eps = eps_default) where {T<:Complex} +function SLD_liouville(ρ::Matrix{T}, ∂ρ_∂x::Matrix{T}; eps = GLOBAL_EPS) where {T<:Complex} 2 * pinv(kron(ρ |> transpose, ρ |> one) + kron(ρ |> one, ρ), rtol = eps) * vec(∂ρ_∂x) |> vec2mat end -function SLD_liouville(ρ::Vector{T}, ∂ρ_∂x::Vector{T}; eps = eps_default) where {T<:Complex} +function SLD_liouville(ρ::Vector{T}, ∂ρ_∂x::Vector{T}; eps = GLOBAL_EPS) where {T<:Complex} SLD_liouville(ρ |> vec2mat, ∂ρ_∂x |> vec2mat; eps = eps) end function SLD_liouville( ρ::Matrix{T}, ∂ρ_∂x::Vector{Matrix{T}}; - eps = eps_default, + eps = GLOBAL_EPS, ) where {T<:Complex} (x -> SLD_liouville(ρ, x; eps = eps)).(∂ρ_∂x) @@ -68,38 +92,145 @@ function SLD_qr(ρ::Matrix{T}, ∂ρ_∂x::Matrix{T}) where {T<:Complex} vec2mat end -function RLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = eps_default) where {T<:Complex} +@doc raw""" + +RLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; rep = "original", eps = GLOBAL_EPS) where {T<:Complex} + +Calculate the right logarrithmic derivatives (RLDs). The RLD operator is defined as ``\partial_{a}\rho=\rho \mathcal{R}_a``, where ρ is the parameterized density matrix. + +- `ρ`: Density matrix. +- `dρ`: Derivatives of the density matrix with respect to the unknown parameters to be estimated. For example, drho[1] is the derivative vector with respect to the first parameter. +- `rep`: Representation of the RLD operator. Options can be: + - "original" (default) -- The RLD matrix will be written in terms of the same basis as the input density matrix (ρ). + - "eigen" -- The RLD matrix will be written in terms of the eigenbasis of the input ρ. +- `eps`: Machine epsilon + +""" +function RLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; rep = "original", eps = GLOBAL_EPS) where {T<:Complex} +(x -> RLD(ρ, x; rep=rep, eps = eps)).(dρ) +end + +""" + + RLD(ρ::Matrix{T},dρ::Matrix{T};rep = "original",eps = GLOBAL_EPS,) where {T<:Complex} + +When applied to the case of single parameter. +""" +function RLD( + ρ::Matrix{T}, + dρ::Matrix{T}; + rep = "original", + eps = GLOBAL_EPS, +) where {T<:Complex} + +dim = size(ρ)[1] + RLD = Matrix{ComplexF64}(undef, dim, dim) + + val, vec = eigen(ρ) + val = val |> real + RLD_eig = zeros(T, dim, dim) + for fi = 1:dim + for fj = 1:dim + if abs(val[fi]) > eps + RLD_eig[fi, fj] = (vec[:, fi]' * dρ * vec[:, fj]) / val[fi] + end + end + end + RLD_eig[findall(RLD_eig == Inf)] .= 0.0 + + if rep == "original" + RLD = vec * (RLD_eig * vec') + elseif rep == "eigen" + RLD = RLD_eig + end + RLD +end + +function RLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = GLOBAL_EPS) where {T<:Complex} pinv(ρ, rtol = eps) * dρ end -function RLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = eps_default) where {T<:Complex} +function RLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = GLOBAL_EPS) where {T<:Complex} (x -> RLD(ρ, x; eps = eps)).(dρ) end -function LLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = eps_default) where {T<:Complex} +@doc raw""" + +LLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; rep = "original", eps = GLOBAL_EPS) where {T<:Complex} + +Calculate the left logarrithmic derivatives (LLDs). The LLD operator is defined as ``\partial_{a}\rho=\mathcal{R}_a^{\dagger}\rho``, where ρ is the parameterized density matrix. +- `ρ`: Density matrix. +- `dρ`: Derivatives of the density matrix with respect to the unknown parameters to be estimated. For example, drho[1] is the derivative vector with respect to the first parameter. +- `rep`: Representation of the LLD operator. Options can be: + - "original" (default) -- The RLD matrix will be written in terms of the same basis as the input density matrix (ρ). + - "eigen" -- The RLD matrix will be written in terms of the eigenbasis of the input ρ. +- `eps`: Machine epsilon + +""" +function LLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; rep = "original", eps = GLOBAL_EPS) where {T<:Complex} + (x -> LLD(ρ, x; rep=rep, eps = eps)).(dρ) +end + +""" + +LLD(ρ::Matrix{T},dρ::Matrix{T};rep = "original",eps = GLOBAL_EPS,) where {T<:Complex} + +When applied to the case of single parameter. +""" +function LLD( + ρ::Matrix{T}, + dρ::Matrix{T}; + rep = "original", + eps = GLOBAL_EPS, +) where {T<:Complex} + + dim = size(ρ)[1] + LLD = Matrix{ComplexF64}(undef, dim, dim) + + val, vec = eigen(ρ) + val = val |> real + LLD_eig = zeros(T, dim, dim) + for fi = 1:dim + for fj = 1:dim + if abs(val[fj]) > eps + LLD_eig[fj, fi] = ((vec[:, fi]' * dρ * vec[:, fj]) / val[fj]) |> conj() + end + end + end + LLD_eig[findall(LLD_eig == Inf)] .= 0.0 + + if rep == "original" + LLD = vec * (LLD_eig * vec') + elseif rep == "eigen" + LLD = LLD_eig + end + LLD +end + +function LLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = GLOBAL_EPS) where {T<:Complex} (dρ * pinv(ρ, rtol = eps)) end -function LLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = eps_default) where {T<:Complex} +function LLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = GLOBAL_EPS) where {T<:Complex} (x -> LLD(ρ, x; eps = eps)).(dρ) end #========================================================# ####################### calculate QFI #################### -function QFIM_SLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = eps_default) where {T<:Complex} +function QFIM_SLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = GLOBAL_EPS) where {T<:Complex} SLD_tp = SLD(ρ, dρ; eps = eps) SLD2_tp = SLD_tp * SLD_tp F = tr(ρ * SLD2_tp) F |> real end -function QFIM_RLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = eps_default) where {T<:Complex} +function QFIM_RLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = GLOBAL_EPS) where {T<:Complex} RLD_tp = RLD(ρ, dρ; eps = eps) F = tr(ρ * RLD_tp * RLD_tp') F |> real end -function QFIM_LLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = eps_default) where {T<:Complex} +function QFIM_LLD(ρ::Matrix{T}, dρ::Matrix{T}; eps = GLOBAL_EPS) where {T<:Complex} LLD_tp = LLD(ρ, dρ; eps = eps) F = tr(ρ * LLD_tp * LLD_tp') F |> real @@ -114,7 +245,7 @@ end #==========================================================# ####################### calculate QFIM ##################### -function QFIM_SLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = eps_default) where {T<:Complex} +function QFIM_SLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = GLOBAL_EPS) where {T<:Complex} p_num = length(dρ) LD_tp = SLD(ρ, dρ; eps = eps) ( @@ -123,14 +254,14 @@ function QFIM_SLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = eps_default) wher ) .|> tr .|> real end -function QFIM_RLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = eps_default) where {T<:Complex} +function QFIM_RLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = GLOBAL_EPS) where {T<:Complex} p_num = length(dρ) LD_tp = RLD(ρ, dρ; eps = eps) LD_dag = [LD_tp[i]' for i = 1:p_num] ([ρ] .* (kron(LD_tp, reshape(LD_dag, 1, p_num)))) .|> tr .|> real end -function QFIM_LLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = eps_default) where {T<:Complex} +function QFIM_LLD(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = GLOBAL_EPS) where {T<:Complex} p_num = length(dρ) LD_tp = LLD(ρ, dρ; eps = eps) LD_dag = [LD_tp[i]' for i = 1:p_num] @@ -159,8 +290,36 @@ end #======================================================# #################### calculate CFIM #################### -function CFIM(ρ::Matrix{T}, dρ::Matrix{T}; eps = eps_default) where {T<:Complex} - M = SIC(size(ρ)[1]) +@doc raw""" + + CFIM(ρ::Matrix{T}, dρ::Vector{Matrix{T}}, M; eps = GLOBAL_EPS) where {T<:Complex} + +Calculate the classical Fisher information matrix (CFIM). + +- `ρ`: Density matrix. +- `dρ`: Derivatives of the density matrix with respect to the unknown parameters to be estimated. For example, drho[1] is the derivative vector with respect to the first parameter. +- `M`: A set of positive operator-valued measure (POVM). The default measurement is a set of rank-one symmetric informationally complete POVM (SIC-POVM). +- `eps`: Machine epsilon. +""" +function CFIM(ρ::Matrix{T}, dρ::Vector{Matrix{T}}, M; eps = GLOBAL_EPS) where {T<:Complex} + m_num = length(M) + p_num = length(dρ) + [ + real(tr(ρ * M[i])) < eps ? zeros(ComplexF64, p_num, p_num) : + (kron(tr.(dρ .* [M[i]]), reshape(tr.(dρ .* [M[i]]), 1, p_num)) / tr(ρ * M[i])) for + i = 1:m_num + ] |> + sum .|> + real +end + +""" + + CFIM(ρ::Matrix{T}, dρ::Matrix{T}, M; eps = GLOBAL_EPS) where {T<:Complex} + +When applied to the case of single parameter. Calculate the classical Fisher information(CFI). +""" +function CFIM(ρ::Matrix{T}, dρ::Matrix{T}, M; eps = GLOBAL_EPS) where {T<:Complex} m_num = length(M) F = 0.0 for i = 1:m_num @@ -174,9 +333,15 @@ function CFIM(ρ::Matrix{T}, dρ::Matrix{T}; eps = eps_default) where {T<:Comple F += cadd end real(F) -end +end + +""" -function CFIM(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = eps_default) where {T<:Complex} + CFIM(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = GLOBAL_EPS) where {T<:Complex} + +When the set of POVM is not given. Calculate the CFIM with SIC-POVM. The SIC-POVM is generated from the Weyl-Heisenberg covariant SIC-POVM fiducial state which can be downloaded from the [website](http://www.physics.umb.edu/Research/QBism/solutions.html). +""" +function CFIM(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = GLOBAL_EPS) where {T<:Complex} M = SIC(size(ρ)[1]) m_num = length(M) p_num = length(dρ) @@ -189,7 +354,14 @@ function CFIM(ρ::Matrix{T}, dρ::Vector{Matrix{T}}; eps = eps_default) where {T real end -function CFIM(ρ::Matrix{T}, dρ::Matrix{T}, M; eps = eps_default) where {T<:Complex} +""" + + CFIM(ρ::Matrix{T}, dρ::Matrix{T}; eps = GLOBAL_EPS) where {T<:Complex} + +When applied to the case of single parameter and the set of POVM is not given. Calculate the CFI with SIC-POVM. +""" +function CFIM(ρ::Matrix{T}, dρ::Matrix{T}; eps = GLOBAL_EPS) where {T<:Complex} + M = SIC(size(ρ)[1]) m_num = length(M) F = 0.0 for i = 1:m_num @@ -205,25 +377,14 @@ function CFIM(ρ::Matrix{T}, dρ::Matrix{T}, M; eps = eps_default) where {T<:Com real(F) end -function CFIM(ρ::Matrix{T}, dρ::Vector{Matrix{T}}, M; eps = eps_default) where {T<:Complex} - m_num = length(M) - p_num = length(dρ) - [ - real(tr(ρ * M[i])) < eps ? zeros(ComplexF64, p_num, p_num) : - (kron(tr.(dρ .* [M[i]]), reshape(tr.(dρ .* [M[i]]), 1, p_num)) / tr(ρ * M[i])) for - i = 1:m_num - ] |> - sum .|> - real -end - +""" -## QFI +""" function QFIM( ρ::Matrix{T}, dρ::Matrix{T}; LDtype = :SLD, - eps = eps_default, + eps = GLOBAL_EPS, ) where {T<:Complex} eval(Symbol("QFIM_" * string(LDtype)))(ρ, dρ; eps = eps) @@ -234,7 +395,7 @@ function QFIM( ρ::Matrix{T}, dρ::Vector{Matrix{T}}; LDtype = :SLD, - eps = eps_default, + eps = GLOBAL_EPS, ) where {T<:Complex} eval(Symbol("QFIM_" * string(LDtype)))(ρ, dρ; eps = eps) @@ -249,7 +410,7 @@ function QFIM( ρ::Matrix{T}, dρ::Matrix{T}; LDtype = :SLD, - eps = eps_default, + eps = GLOBAL_EPS, ) where {T<:Complex} F = QFIM(ρ, dρ; LDtype = LDtype, eps = eps) @@ -263,7 +424,7 @@ function QFIM( ρ::Matrix{T}, dρ::Vector{Matrix{T}}; LDtype = :SLD, - eps = eps_default, + eps = GLOBAL_EPS, ) where {T<:Complex} F = QFIM(ρ, dρ; LDtype = LDtype, eps = eps) @@ -272,6 +433,22 @@ function QFIM( end +""" + + QFIM_Bloch(r, dr; eps = 1e-8) + +Calculate the SLD based quantum Fisher information (QFI) or quantum +Fisher information matrix (QFIM) in Bloch representation. + +- `r`: Parameterized Bloch vector. + +- `dr`: Derivative(s) of the Bloch vector with respect to the unknown parameters to be + estimated. For example, dr[1] is the derivative vector with respect to the first + parameter. + +- `eps`: Machine epsilon. + +""" ## TODO: 👇 check type stability function QFIM_Bloch(r, dr; eps = 1e-8) para_num = length(dr) @@ -378,6 +555,26 @@ function G_Gauss(S::M, dC::VM, c::V) where {M<:AbstractMatrix,V,VM<:AbstractVect G end +""" + + QFIM_Gauss(R̄::V, dR̄::VV, D::M, dD::VM) where {V,VV,M,VM<:AbstractVecOrMat} + +Calculate the SLD based quantum Fisher information matrix (QFIM) with gaussian states. + +- `R̄` : First-order moment. + +- `dR̄`: Derivatives of the first-order moment with respect to the unknown parameters to be +estimated. For example, dR[1] is the derivative vector on the first +parameter. + +- `D`: Second-order moment. + +- `dD`: Derivatives of the second-order moment with respect to the unknown parameters to be +estimated. + +- `eps`: Machine epsilon + +""" function QFIM_Gauss(R̄::V, dR̄::VV, D::M, dD::VM) where {V,VV,M,VM<:AbstractVecOrMat} para_num = length(dR̄) quad_num = length(R̄) diff --git a/src/objective/AsymptoticBound/Holevo.jl b/src/objective/AsymptoticBound/Holevo.jl index 5dfa8d0..beba772 100644 --- a/src/objective/AsymptoticBound/Holevo.jl +++ b/src/objective/AsymptoticBound/Holevo.jl @@ -4,33 +4,48 @@ function decomposition(A) return R end +""" + + HCRB(ρ::Matrix{T},dρ::Vector{Matrix{T}},C::Matrix{Float64};eps = 1e-8,) where {T<:Complex} + +Caltulate the Holevo Cramer-Rao bound (HCRB) via the semidefinite program (SDP). + +- `ρ`: Density matrix. + +- `dρ`: Derivatives of the density matrix on the unknown parameters to be estimated. For example, drho[0] is the derivative vector on the first parameter. + +- `W`: Weight matrix. + +- `eps`: Machine epsilon. + +""" function HCRB( ρ::Matrix{T}, - ∂ρ_∂x::Vector{Matrix{T}}, + dρ::Vector{Matrix{T}}, C::Matrix{Float64}; - eps = 1e-6, + eps = 1e-8, ) where {T<:Complex} - if length(∂ρ_∂x) == 1 + if length(dρ) == 1 println( - "In single parameter scenario, HCRB is equivalent to QFI. This function will return the value of QFI", + "In the single-parameter scenario, the HCRB is equivalent to the QFI. This function will return the value of the QFI", ) - f = QFIM_SLD(ρ, ∂ρ_∂x[1]; eps=eps) + f = QFIM_SLD(ρ, dρ[1]; eps=eps) return f else - Holevo_bound(ρ, ∂ρ_∂x, C; eps = eps) + Holevo_bound(ρ, dρ, C; eps = eps) end end function Holevo_bound( ρ::Matrix{T}, - ∂ρ_∂x::Vector{Matrix{T}}, + dρ::Vector{Matrix{T}}, C::Matrix{Float64}; eps = eps_default, ) where {T<:Complex} dim = size(ρ)[1] num = dim * dim - para_num = length(∂ρ_∂x) + para_num = length(dρ) suN = suN_generator(dim) / sqrt(2) Lambda = [Matrix{ComplexF64}(I, dim, dim) / sqrt(2)] append!(Lambda, [suN[i] for i = 1:length(suN)]) @@ -38,7 +53,7 @@ function Holevo_bound( for pa = 1:para_num for ra = 2:num - vec_∂ρ[pa][ra] = (∂ρ_∂x[pa] * Lambda[ra]) |> tr |> real + vec_∂ρ[pa][ra] = (dρ[pa] * Lambda[ra]) |> tr |> real end end S = zeros(ComplexF64, num, num) @@ -66,15 +81,15 @@ function Holevo_bound( end end problem = minimize(tr(C * V), constraints) - Convex.solve!(problem, SCS.Optimizer(verbose = false)) - return evaluate(tr(C * V)), evaluate(X), evaluate(V) + Convex.solve!(problem, SCS.Optimizer, silent_solver=true) + return evaluate(tr(C * V)) end function Holevo_bound_obj( ρ::Matrix{T}, - ∂ρ_∂x::Vector{Matrix{T}}, + dρ::Vector{Matrix{T}}, C::Matrix{Float64}; eps = eps_default, ) where {T<:Complex} - return Holevo_bound(ρ, ∂ρ_∂x, C; eps = eps)[1] + return Holevo_bound(ρ, dρ, C; eps = eps)[1] end diff --git a/src/objective/BayesianBound/BayesianCramerRao.jl b/src/objective/BayesianBound/BayesianCramerRao.jl index 619736c..8c4c834 100644 --- a/src/objective/BayesianBound/BayesianCramerRao.jl +++ b/src/objective/BayesianBound/BayesianCramerRao.jl @@ -1,6 +1,55 @@ +########## Bayesian version of CFIM ########## +function BCFIM(x, p, rho, drho; M::Union{AbstractVector,Nothing}=nothing, eps=1e-8) + para_num = length(x) + if para_num == 1 + #### singleparameter scenario #### + p_num = length(p) + if M==nothing + M = SIC(size(rho[1])[1]) + end + F_tp = zeros(p_num) + for i in 1:p_num + F_tp[i] = CFIM(rho[i], drho[i], M; eps=eps) + end + F = 0.0 + arr = [p[i]*F_tp[i] for i in 1:p_num] + F = trapz(x[1], arr) + else + #### multiparameter scenario #### + if M==nothing + M = SIC(size(vec(rho)[1])[1]) + end + + xnum = length(x) + trapzm(x, integrands, slice_dim) = [trapz(tuple(x...), I) for I in [reshape(hcat(integrands...)[i,:], length.(x)...) for i in 1:slice_dim]] + Fs = [p*CFIM(rho,drho,M;eps=eps)|>vec for (p,rho,drho) in zip(p,rho,drho)] + F = trapzm(x, Fs, xnum^2) |> I->reshape(I,xnum,xnum) + end +end +########## Bayesian version of QFIM ########## +function BQFIM(x, p, rho, drho; LDtype=:SLD, eps=1e-8) + para_num = length(x) + if para_num == 1 + #### singleparameter scenario #### + p_num = length(p) + F_tp = zeros(p_num) + for i in 1:p_num + F_tp[i] = QFIM(rho[i], drho[i]; LDtype=LDtype, eps=eps) + end + F = 0.0 + arr = [p[i]*F_tp[i] for i in 1:p_num] + F = trapz(x[1], arr) + else + #### multiparameter scenario #### + xnum = length(x) + trapzm(x, integrands, slice_dim) = [trapz(tuple(x...), I) for I in [reshape(hcat(integrands...)[i,:], length.(x)...) for i in 1:slice_dim]] + Fs = [p*QFIM(rho,drho;LDtype=LDtype,eps=eps)|>vec for (p,rho,drho) in zip(p,rho,drho)] + F = trapzm(x, Fs, xnum^2) |> I->reshape(I,xnum,xnum) + end +end ########## Bayesian quantum Cramer-Rao bound ########## -function BQCRB(x, p, rho, drho; b=nothing, db=nothing, LDtype="SLD", btype=1, eps=1e-8) +function BQCRB(x, p, rho, drho; b=nothing, db=nothing, LDtype=:SLD, btype=1, eps=1e-8) para_num = length(x) if b==nothing @@ -13,7 +62,7 @@ function BQCRB(x, p, rho, drho; b=nothing, db=nothing, LDtype="SLD", btype=1, ep if para_num == 1 - #### singleparameter senario #### + #### singleparameter scenario #### p_num = length(p) if typeof(drho[1]) == Vector{Matrix{ComplexF64}} @@ -45,7 +94,7 @@ function BQCRB(x, p, rho, drho; b=nothing, db=nothing, LDtype="SLD", btype=1, ep end return F else - #### multiparameter senario #### + #### multiparameter scenario #### xnum = length(x) bs = Iterators.product(b...) @@ -80,7 +129,7 @@ function BCRB(x, p, rho, drho; M::Union{AbstractVector,Nothing}=nothing, b=nothi end if para_num == 1 - #### singleparameter senario #### + #### singleparameter scenario #### p_num = length(p) if M==nothing M = SIC(size(rho[1])[1]) @@ -114,7 +163,7 @@ function BCRB(x, p, rho, drho; M::Union{AbstractVector,Nothing}=nothing, b=nothi end return F else - #### multiparameter senario #### + #### multiparameter scenario #### if M==nothing M = SIC(size(vec(rho)[1])[1]) end @@ -140,10 +189,10 @@ function BCRB(x, p, rho, drho; M::Union{AbstractVector,Nothing}=nothing, b=nothi end end -function QVTB(x, p, dp, rho, drho; LDtype="SLD", btype=1, eps=1e-8) +function QVTB(x, p, dp, rho, drho; LDtype=:SLD, btype=1, eps=1e-8) para_num = length(x) if para_num == 1 - #### singleparameter senario #### + #### singleparameter scenario #### p_num = length(p) if typeof(drho[1]) == Vector{Matrix{ComplexF64}} drho = [drho[i][1] for i in 1:p_num] @@ -170,7 +219,7 @@ function QVTB(x, p, dp, rho, drho; LDtype="SLD", btype=1, eps=1e-8) end return I else - #### multiparameter senario #### + #### multiparameter scenario #### xnum = length(x) trapzm(x, integrands, slice_dim) = [trapz(tuple(x...), I) for I in [reshape(hcat(integrands...)[i,:], length.(x)...) for i in 1:slice_dim]] Ip(p,dp) = dp*dp'/p^2 @@ -193,7 +242,7 @@ end function VTB(x, p, dp, rho, drho; M::Union{AbstractVector,Nothing}=nothing, btype=1, eps=1e-8) para_num = length(x) if para_num == 1 - #### singleparameter senario #### + #### singleparameter scenario #### p_num = length(p) if M==nothing M = SIC(size(rho[1])[1]) @@ -222,7 +271,7 @@ function VTB(x, p, dp, rho, drho; M::Union{AbstractVector,Nothing}=nothing, btyp end return res else - #### multiparameter senario #### + #### multiparameter scenario #### if M==nothing M = SIC(size(vec(rho)[1])[1]) end @@ -267,7 +316,7 @@ function interp1(xspan, yspan, x) return y end -function OBB(x, p, dp, rho, drho, d2rho; LDtype="SLD", eps=1e-8) +function OBB(x, p, dp, rho, drho, d2rho; LDtype=:SLD, eps=1e-8) p_num = length(p) if typeof(drho[1]) == Vector{Matrix{ComplexF64}} diff --git a/src/optim/optim.jl b/src/optim/optim.jl index 382b20a..c91f135 100644 --- a/src/optim/optim.jl +++ b/src/optim/optim.jl @@ -2,7 +2,7 @@ abstract type AbstractOpt end abstract type AbstractMeasurementType end abstract type Projection <: AbstractMeasurementType end -abstract type LinearComb <: AbstractMeasurementType end +abstract type LC <: AbstractMeasurementType end abstract type Rotation <: AbstractMeasurementType end abstract type Opt <: AbstractOpt end @@ -16,7 +16,7 @@ Copt = ControlOpt ControlOpt(ctrl::Matrix{R}, ctrl_bound::AbstractVector) where R<:Number = ControlOpt([c[:] for c in eachrow(ctrl)], ctrl_bound) Base.@kwdef mutable struct StateOpt <: Opt - ψ₀::Union{AbstractVector, Missing} = missing + psi::Union{AbstractVector, Missing} = missing end Sopt = StateOpt @@ -24,7 +24,7 @@ Sopt = StateOpt abstract type AbstractMopt <: Opt end Base.@kwdef mutable struct Mopt_Projection <: AbstractMopt - C::Union{AbstractVector, Missing} = missing + M::Union{AbstractVector, Missing} = missing end Base.@kwdef mutable struct Mopt_LinearComb <: AbstractMopt B::Union{AbstractVector, Missing} = missing @@ -37,21 +37,22 @@ Base.@kwdef mutable struct Mopt_Rotation <: AbstractMopt Lambda::Union{AbstractVector, Missing} = missing end -function MeasurementOpt(;method=:Projection,kwargs...) - if method==:Projection +function MeasurementOpt(;mtype=:Projection, kwargs...) + if mtype==:Projection return Mopt_Projection(;kwargs...) - elseif method==:LinearCombination + elseif mtype==:LC return Mopt_LinearComb(;kwargs...) - elseif method==:Rotation + elseif mtype==:Rotation return Mopt_Rotation(;kwargs...) end end + Mopt = MeasurementOpt abstract type CompOpt <: Opt end Base.@kwdef mutable struct StateControlOpt <: CompOpt - ψ₀::Union{AbstractVector, Missing} = missing + psi::Union{AbstractVector, Missing} = missing ctrl::Union{AbstractVector, Missing} = missing ctrl_bound::AbstractVector = [-Inf, Inf] end @@ -60,23 +61,23 @@ SCopt = StateControlOpt Base.@kwdef mutable struct ControlMeasurementOpt <: CompOpt ctrl::Union{AbstractVector, Missing} = missing - C::Union{AbstractVector, Missing} = missing + M::Union{AbstractVector, Missing} = missing ctrl_bound::AbstractVector = [-Inf, Inf] end CMopt = ControlMeasurementOpt Base.@kwdef mutable struct StateMeasurementOpt <: CompOpt - ψ₀::Union{AbstractVector, Missing} = missing - C::Union{AbstractVector, Missing} = missing + psi::Union{AbstractVector, Missing} = missing + M::Union{AbstractVector, Missing} = missing end SMopt = StateMeasurementOpt Base.@kwdef mutable struct StateControlMeasurementOpt <: CompOpt ctrl::Union{AbstractVector, Missing} = missing - ψ₀::Union{AbstractVector, Missing} = missing - C::Union{AbstractVector, Missing} = missing + psi::Union{AbstractVector, Missing} = missing + M::Union{AbstractVector, Missing} = missing ctrl_bound::AbstractVector = [-Inf, Inf] end @@ -94,14 +95,14 @@ opt_target(::StateMeasurementOpt) = :SMopt opt_target(::StateControlMeasurementOpt) = :SCMopt result(opt::ControlOpt) = [opt.ctrl] -result(opt::StateOpt) = [opt.ψ₀] -result(opt::Mopt_Projection) = [opt.C] +result(opt::StateOpt) = [opt.psi] +result(opt::Mopt_Projection) = [opt.M] result(opt::Mopt_LinearComb) = [opt.B, opt.POVM_basis, opt.M_num] result(opt::Mopt_Rotation) = [opt.s] -result(opt::StateControlOpt) = [opt.ψ₀, opt.ctrl] -result(opt::ControlMeasurementOpt) = [opt.ctrl, opt.C] -result(opt::StateMeasurementOpt) = [opt.ψ₀, opt.C] -result(opt::StateControlMeasurementOpt) = [opt.ψ₀, opt.ctrl, opt.C] +result(opt::StateControlOpt) = [opt.psi, opt.ctrl] +result(opt::ControlMeasurementOpt) = [opt.ctrl, opt.M] +result(opt::StateMeasurementOpt) = [opt.psi, opt.M] +result(opt::StateControlMeasurementOpt) = [opt.psi, opt.ctrl, opt.M] #with reward result(opt, ::Type{Val{:save_reward}}) = [result(opt)..., [0.0]]