regist

Project.toml

@@ -1,6 +1,6 @@
 name = "QuanEstimation"
 uuid = "088c8dff-a786-4a66-974c-03d3f6773f87"
-authors = ["Hauiming Yu <Huaimingyuuu@gmail.com> and contributors"]
+authors = ["Jing Liu <liujingphys@hust.edu.cn>", "Hauiming Yu <Huaimingyuuu@gmail.com>", "Mao Zhang <zhangmao2018@hust.edu.cn>"]
 version = "0.1.0"
 
 [deps]
@@ -24,7 +24,7 @@ Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
-julia = "1.6.2"
+julia = "1.7.2"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -1,6 +1,22 @@
-# QuanEstimation
+# QuanEstimation.jl
 
+<!-- needs to be modified -->
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://HuaimingYuuu.github.io/QuanEstimation.jl/stable)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://HuaimingYuuu.github.io/QuanEstimation.jl/dev)
 [![Build Status](https://github.com/HuaimingYuuu/QuanEstimation.jl/workflows/CI/badge.svg)](https://github.com/HuaimingYuuu/QuanEstimation.jl/actions)
 [![Coverage](https://codecov.io/gh/HuaimingYuuu/QuanEstimation.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/HuaimingYuuu/QuanEstimation.jl)
+<!-- needs to be modified -->
+QuanEstimation is a Python-Julia based open-source toolkit for quantum parameter estimation, which consist in the calculation of the quantum metrological tools and quantum resources, the optimization of the probe state, control and measurement in quantum metrology. Futhermore, QuanEstimation can also perform comprehensive optimization with respect to the probe state, control and measurement to generate not only optimal quantum parameter estimation schemes, but also adaptive measurement schemes.
+
+This package is a fully Julia implementation of [QuanEstimation](https://github.com/QuanEstimation/QuanEstimation).
+
+## Installation
+
+Run the command in the terminal to install QuanEstimation:  
+
+~~~
+pkg> add QuanEstimation
+~~~
+
+## License
+<a rel="license" href="http://creativecommons.org/licenses/by-nc-sa/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-nc-sa/4.0/88x31.png" /></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-nc-sa/4.0/">Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License</a>.

examples/Bayesian_CramerRao_bounds.jl

@@ -0,0 +1,58 @@
+using QuanEstimation
+using Trapz
+
+# free Hamiltonian
+function H0_func(x)
+    return 0.5*B*omega0*(sx*cos(x)+sz*sin(x))
+end
+# derivative of the free Hamiltonian on x
+function dH_func(x)
+    return [0.5*B*omega0*(-sx*sin(x)+sz*cos(x))]
+end
+# 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
+
+B, omega0 = 0.5*pi, 1.0
+sx = [0. 1.; 1. 0.0im]
+sy = [0. -im; im 0.]
+sz = [1. 0.0im; 0. -1.]
+# initial state
+rho0 = 0.5*ones(2, 2)
+# prior distribution
+x = range(-0.5*pi, stop=0.5*pi, length=100) |>Vector
+mu, eta = 0.0, 0.2
+p_tp = [p_func(x[i], mu, eta) for i in 1:length(x)]
+dp_tp = [dp_func(x[i], mu, eta) for i in 1:length(x)]
+# normalization of the distribution
+c = trapz(x, p_tp)
+p = p_tp/c
+dp = dp_tp/c
+# time length for the evolution
+tspan = range(0., stop=1., length=1000)
+# dynamics
+rho = Vector{Matrix{ComplexF64}}(undef, length(x))
+drho = Vector{Vector{Matrix{ComplexF64}}}(undef, length(x))
+for i = 1:length(x) 
+    H0_tp = H0_func(x[i])
+    dH_tp = dH_func(x[i])
+    rho_tp, drho_tp = QuanEstimation.expm(tspan, rho0, H0_tp, dH_tp)
+    rho[i], drho[i] = rho_tp[end], drho_tp[end]
+end
+
+# Classical Bayesian bounds
+f_BCRB1 = QuanEstimation.BCRB([x], p, [], rho, drho, btype=1)
+f_BCRB2 = QuanEstimation.BCRB([x], p, [], rho, drho, btype=2)
+f_BCRB3 = QuanEstimation.BCRB([x], p, dp, rho, drho, btype=3)
+f_VTB = QuanEstimation.VTB([x], p, dp, rho, drho)
+
+# Quantum Bayesian bounds
+f_BQCRB1 = QuanEstimation.BQCRB([x], p, [], rho, drho, btype=1)
+f_BQCRB2 = QuanEstimation.BQCRB([x], p, [], rho, drho, btype=2)
+f_BQCRB3 = QuanEstimation.BQCRB([x], p, dp, rho, drho, btype=3)
+f_QVTB = QuanEstimation.QVTB([x], p, dp, rho, drho)
+f_QZZB = QuanEstimation.QZZB([x], p, rho)

examples/Bayesian_estimation.jl

@@ -0,0 +1,51 @@
+using QuanEstimation
+using Random
+using StatsBase
+
+# free Hamiltonian
+function H0_func(x)
+    return 0.5*B*omega0*(sx*cos(x)+sz*sin(x))
+end
+# derivative of the free Hamiltonian on x
+function dH_func(x)
+    return [0.5*B*omega0*(-sx*sin(x)+sz*cos(x))]
+end
+
+B, omega0 = pi/2.0, 1.0
+sx = [0. 1.; 1. 0.0im]
+sy = [0. -im; im 0.]
+sz = [1. 0.0im; 0. -1.]
+# initial state
+rho0 = 0.5*ones(2, 2)
+# measurement 
+M1 = 0.5*[1.0+0.0im  1.; 1.  1.]
+M2 = 0.5*[1.0+0.0im -1.; -1.  1.]
+M = [M1, M2]
+# prior distribution
+x = range(0., stop=0.5*pi, length=100) |>Vector
+p = (1.0/(x[end]-x[1]))*ones(length(x))
+# time length for the evolution
+tspan = range(0., stop=1., length=1000)
+# dynamics
+rho = Vector{Matrix{ComplexF64}}(undef, length(x))
+for i = 1:length(x) 
+    H0_tp = H0_func(x[i])
+    dH_tp = dH_func(x[i])
+    rho_tp, drho_tp = QuanEstimation.expm(tspan, rho0, H0_tp, dH_tp)
+    rho[i] = rho_tp[end]
+end
+
+# Generation of the 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
+
+#===============Maximum a posteriori estimation===============#
+pout, xout = QuanEstimation.Bayes([x], p, rho, y; M=M, estimator="MAP",
+                                  savefile=false)
+
+#===============Maximum likelihood estimation===============#
+Lout, xout = QuanEstimation.MLE([x], rho, y, M=M; savefile=false)

examples/CMopt_NV.jl

@@ -0,0 +1,63 @@
+using QuanEstimation
+using Random
+using LinearAlgebra
+
+# initial state
+rho0 = zeros(ComplexF64, 6, 6)
+rho0[1:4:5, 1:4:5] .= 0.5
+# 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]
+# All numbers are divided by 100 in this example 
+# for better calculation accurancy
+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)])
+# derivatives of the free Hamiltonian on B1, B2 and B3
+dH = gS*S+gI*Is
+# control Hamiltonians 
+Hc = [S1, S2, S3]
+# dissipation
+decay = [[S3, 2pi/cons]]
+# generation of a set of POVM basis
+dim = size(rho0, 1)
+POVM_basis = [QuanEstimation.basis(dim, i)*QuanEstimation.basis(dim, i)' 
+              for i in 1:dim]
+# time length for the evolution
+tspan = range(0., 2., length=4000)
+# control and measurement optimization
+opt = QuanEstimation.CMopt(ctrl_bound=[-0.2,0.2], seed=1234)
+
+##==========choose comprehensive optimization algorithm==========##
+##-------------algorithm: DE---------------------##
+alg = QuanEstimation.DE(p_num=10, max_episode=1000, c=1.0, cr=0.5)
+# input the dynamics data
+dynamics = QuanEstimation.Lindblad(opt, tspan, rho0, H0, dH, Hc, 
+                                   decay=decay)   
+# objective function: CFI
+obj = QuanEstimation.CFIM_obj() 
+# run the comprehensive optimization problem
+QuanEstimation.run(opt, alg, obj, dynamics; savefile=false)
+
+##-------------algorithm: PSO---------------------##
+# alg = QuanEstimation.PSO(p_num=10, max_episode=[1000,100], c0=1.0, 
+#                          c1=2.0, c2=2.0)
+# # input the dynamics data
+# dynamics = QuanEstimation.Lindblad(opt, tspan, rho0, H0, dH, Hc, 
+#                                    decay=decay)   
+# # objective function: CFI
+# obj = QuanEstimation.CFIM_obj() 
+# # run the comprehensive optimization problem
+# QuanEstimation.run(opt, alg, obj, dynamics; savefile=false)

examples/CMopt_qubit.jl

@@ -0,0 +1,48 @@
+using QuanEstimation
+
+# initial state
+rho0 = 0.5*ones(2, 2)
+# free Hamiltonian
+omega = 1.0
+sx = [0. 1.; 1. 0.0im]
+sy = [0. -im; im 0.]
+sz = [1. 0.0im; 0. -1.]
+H0 = 0.5*omega*sz
+# derivative of the free Hamiltonian on omega
+dH = [0.5*sz]
+# control Hamiltonians 
+Hc = [sx, sy, sz]
+# dissipation
+sp = [0. 1.; 0. 0.0im]
+sm = [0. 0.; 1. 0.0im]
+decay = [[sp, 0.0], [sm, 0.1]]
+# measurement
+M1 = 0.5*[1.0+0.0im  1.; 1.  1.]
+M2 = 0.5*[1.0+0.0im -1.; -1.  1.]
+M = [M1, M2]
+# time length for the evolution
+tspan = range(0., 10., length=2500)
+# control and measurement optimization
+opt = QuanEstimation.CMopt(ctrl_bound=[-2.0,2.0], seed=1234)
+
+##==========choose comprehensive optimization algorithm==========##
+##-------------algorithm: DE---------------------##
+alg = QuanEstimation.DE(p_num=10, max_episode=1000, c=1.0, cr=0.5)
+# input the dynamics data
+dynamics = QuanEstimation.Lindblad(opt, tspan, rho0, H0, dH, Hc, 
+                                   decay=decay)   
+# objective function: CFI
+obj = QuanEstimation.CFIM_obj() 
+# run the comprehensive optimization problem
+QuanEstimation.run(opt, alg, obj, dynamics; savefile=false)
+
+##-------------algorithm: PSO---------------------##
+# alg = QuanEstimation.PSO(p_num=10, max_episode=[1000,100], c0=1.0, 
+#                          c1=2.0, c2=2.0)
+# # input the dynamics data
+# dynamics = QuanEstimation.Lindblad(opt, tspan, rho0, H0, dH, Hc, 
+#                                    decay=decay)   
+# # objective function: CFI
+# obj = QuanEstimation.CFIM_obj() 
+# # run the comprehensive optimization problem
+# QuanEstimation.run(opt, alg, obj, dynamics; savefile=false)

examples/CramerRao_bounds.jl

@@ -0,0 +1,34 @@
+using QuanEstimation
+
+# initial state
+rho0 = 0.5*ones(2, 2)
+# free Hamiltonian
+omega = 1.0
+sx = [0. 1.; 1. 0.0im]
+sy = [0. -im; im 0.]
+sz = [1. 0.0im; 0. -1.]
+H0 = 0.5*omega*sz
+# derivative of the free Hamiltonian on omega
+dH = [0.5*sz]
+# dissipation
+sp = [0. 1.; 0. 0.0im]
+sm = [0. 0.; 1. 0.0im]
+decay = [[sp, 0.0], [sm, 0.1]]
+# measurement
+M1 = 0.5*[1.0+0.0im  1.; 1.  1.]
+M2 = 0.5*[1.0+0.0im -1.; -1.  1.]
+M = [M1, M2]
+# time length for the evolution
+tspan = range(0., 50., length=2000)
+# dynamics
+rho, drho = QuanEstimation.expm(tspan, rho0, H0, dH, decay)
+# calculation of the CFI and QFI
+Im, F = Float64[], Float64[]
+for ti in 2:length(tspan)
+    # CFI
+    I_tp = QuanEstimation.CFIM(rho[ti], drho[ti], M)
+    append!(Im, I_tp)
+    # QFI
+    F_tp = QuanEstimation.QFIM(rho[ti], drho[ti])
+    append!(F, F_tp)
+end

examples/Holevo_CramerRao_bound.jl

@@ -0,0 +1,40 @@
+using QuanEstimation
+using LinearAlgebra
+
+# initial state
+psi0 = [1., 0., 0., 1.]/sqrt(2)
+rho0 = psi0*psi0'
+# free Hamiltonian
+omega1, omega2, g = 1.0, 1.0, 0.1
+sx = [0. 1.; 1. 0.0im]
+sy = [0. -im; im 0.]
+sz = [1. 0.0im; 0. -1.]
+H0 = omega1*kron(sz, I(2)) + omega2*kron(I(2), sz) + g*kron(sx, sx)
+# derivatives of the free Hamiltonian with respect to omega2 and g
+dH = [kron(I(2), sz), kron(sx, sx)]
+# dissipation
+decay = [[kron(sz, I(2)), 0.05], [kron(I(2), sz), 0.05]]
+# measurement
+m1 = [1., 0., 0., 0.]
+M1 = 0.85*m1*m1'
+M2 = 0.1*ones(4, 4)
+M = [M1, M2, I(4)-M1-M2]
+# time length for the evolution
+tspan = range(0., 10., length=1000)
+# dynamics
+rho, drho = QuanEstimation.expm(tspan, rho0, H0, dH, decay)
+# weight matrix
+W = one(zeros(2, 2))
+# calculation of the CFIM, QFIM and HCRB
+Im, F, f = [], [], Float64[]
+for ti in 2:length(tspan)
+    #CFIM
+    I_tp = QuanEstimation.CFIM(rho[ti], drho[ti], M)
+    append!(Im, [I_tp])
+    #QFIM
+    F_tp = QuanEstimation.QFIM(rho[ti], drho[ti])
+    append!(F, [F_tp])
+    #HCRB
+    f_tp = QuanEstimation.HCRB(rho[ti], drho[ti], W)
+    append!(f, f_tp)
+end