wrappers

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

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

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)

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)

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)

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)