From 0101b588c5ff0a1aeffc4c0cb3217182d98ab4a0 Mon Sep 17 00:00:00 2001 From: Hauiming Yu Date: Thu, 12 Aug 2021 19:24:37 +0800 Subject: [PATCH 1/2] first commit --- .github/workflows | 1 + .gitignore | 6 + LICENSE | 21 +++ Project.toml | 22 +++ README.md | 7 +- docs/Manifest.toml | 97 +++++++++++ docs/Project.toml | 3 + docs/make.jl | 23 +++ docs/src/index.md | 12 ++ src/AsymptoticBound/CramerRao.jl | 158 +++++++++++++++++ src/Common/common.jl | 104 ++++++++++++ src/Common/utils.jl | 77 +++++++++ src/Control/GRAPE.jl | 280 +++++++++++++++++++++++++++++++ src/Control/PSO.jl | 118 +++++++++++++ src/Dynamics/Adam.jl | 20 +++ src/Dynamics/dynamcs.jl | 279 ++++++++++++++++++++++++++++++ src/QuanEstimation.jl | 25 +++ test/runtests.jl | 6 + 18 files changed, 1258 insertions(+), 1 deletion(-) create mode 160000 .github/workflows create mode 100644 .gitignore create mode 100644 LICENSE create mode 100644 Project.toml create mode 100644 docs/Manifest.toml create mode 100644 docs/Project.toml create mode 100644 docs/make.jl create mode 100644 docs/src/index.md create mode 100644 src/AsymptoticBound/CramerRao.jl create mode 100644 src/Common/common.jl create mode 100644 src/Common/utils.jl create mode 100644 src/Control/GRAPE.jl create mode 100644 src/Control/PSO.jl create mode 100644 src/Dynamics/Adam.jl create mode 100644 src/Dynamics/dynamcs.jl create mode 100644 src/QuanEstimation.jl create mode 100644 test/runtests.jl diff --git a/.github/workflows b/.github/workflows new file mode 160000 index 0000000..9208e86 --- /dev/null +++ b/.github/workflows @@ -0,0 +1 @@ +Subproject commit 9208e8663ef7707ad74ebf4ea44b1f4e045b25c9 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..eadbc49 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +*.jl.*.cov +*.jl.cov +*.jl.mem +/Manifest.toml +/deps/build.log +/docs/build/ diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..92d9e90 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021 Hauiming Yu and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..f9bddce --- /dev/null +++ b/Project.toml @@ -0,0 +1,22 @@ +name = "QuanEstimation" +uuid = "088c8dff-a786-4a66-974c-03d3f6773f87" +authors = ["Hauiming Yu and contributors"] +version = "0.1.0" + +[deps] +DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[compat] +julia = "1.6.2" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/README.md b/README.md index dc9c6a9..6af71d0 100644 --- a/README.md +++ b/README.md @@ -1 +1,6 @@ -# QuanEstimation.jl \ No newline at end of file +# QuanEstimation + +[![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) diff --git a/docs/Manifest.toml b/docs/Manifest.toml new file mode 100644 index 0000000..7c574be --- /dev/null +++ b/docs/Manifest.toml @@ -0,0 +1,97 @@ +# This file is machine-generated - editing it directly is not advised + +[[ANSIColoredPrinters]] +git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" +uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9" +version = "0.0.1" + +[[Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.8.5" + +[[Documenter]] +deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] +git-tree-sha1 = "350dced36c11f794c6c4da5dc6493ec894e50c16" +uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +version = "0.27.5" + +[[IOCapture]] +deps = ["Logging", "Random"] +git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a" +uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" +version = "0.2.2" + +[[InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[JSON]] +deps = ["Dates", "Mmap", "Parsers", "Unicode"] +git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.21.2" + +[[LibGit2]] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" + +[[Parsers]] +deps = ["Dates"] +git-tree-sha1 = "477bf42b4d1496b454c10cce46645bb5b8a0cf2c" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "2.0.2" + +[[Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[QuanEstimation]] +path = ".." +uuid = "088c8dff-a786-4a66-974c-03d3f6773f87" +version = "0.1.0" + +[[REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[Random]] +deps = ["Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" + +[[Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..5f54614 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,3 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +QuanEstimation = "088c8dff-a786-4a66-974c-03d3f6773f87" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..d356f6d --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,23 @@ +using QuanEstimation +using Documenter + +DocMeta.setdocmeta!(QuanEstimation, :DocTestSetup, :(using QuanEstimation); recursive=true) + +makedocs(; + modules=[QuanEstimation], + authors="Hauiming Yu and contributors", + repo="https://github.com/HuaimingYuuu/QuanEstimation.jl/blob/{commit}{path}#{line}", + sitename="QuanEstimation.jl", + format=Documenter.HTML(; + prettyurls=get(ENV, "CI", "false") == "true", + canonical="https://HuaimingYuuu.github.io/QuanEstimation.jl", + assets=String[], + ), + pages=[ + "Home" => "index.md", + ], +) + +deploydocs(; + repo="github.com/HuaimingYuuu/QuanEstimation.jl", +) diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..7d82778 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,12 @@ +```@meta +CurrentModule = QuanEstimation +``` + +# QuanEstimation + +```@index +``` + +```@autodocs +Modules = [QuanEstimation] +``` diff --git a/src/AsymptoticBound/CramerRao.jl b/src/AsymptoticBound/CramerRao.jl new file mode 100644 index 0000000..64ec15e --- /dev/null +++ b/src/AsymptoticBound/CramerRao.jl @@ -0,0 +1,158 @@ +function CFI(ρ, dρ, M) + m_num = length(M) + p = zero(ComplexF64) + dp = zero(ComplexF64) + F = 0. + for i in 1:m_num + mp = M[i] + p += tr(ρ * mp) + dp = tr(dρ * mp) + cadd = 0. + if p != 0 + cadd = (dp^2) / p + end + F += cadd + end + real(F) +end +function CFI(M::Vector{Matrix{T}}, H::Vector{Matrix{T}}, ∂H_∂x::Matrix{T}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, times) where {T <: Complex,R <: Real} + dim = size(H[1])[1] + Δt = times[2] - times[1] + ρt = evolute(H[1], Liouville_operator, γ, times, times[1]) * (ρ_initial |> vec) + ∂ρt_∂x = -im * Δt * liouville_commu(∂H_∂x) * ρt + for t in 2:length(times) + expL = evolute(H[t], Liouville_operator, γ, times, times[t]) + ρt= expL * ρt + ∂ρt_∂x= -im * Δt * liouville_commu(∂H_∂x) * ρt + expL * ∂ρt_∂x + end + CFI(ρt|> vec2mat, ∂ρt_∂x|> vec2mat, M) +end + +function CFIM(ρ, dρ, M) + m_num = length(M) + cfim = [tr.(kron(dρ', dρ).*M[i]) / tr(ρ*M[i]) for i in 1:m_num] |> sum +end +function CFIM(M::Vector{Matrix{T}}, H::Vector{Matrix{T}}, ∂H_∂x::Matrix{T}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, times) where {T <: Complex,R <: Real} + dim = size(H[1])[1] + Δt = times[2] - times[1] + para_num = length(∂H_∂x) + ρt = evolute(H[1], Liouville_operator, γ, times, times[1]) * ρ_initial[:] + ∂ρt_∂x = [-im * Δt * liouville_commu(∂H_∂x[i]) * ρt for i in 1:para_num] + for t in 2:length(times) + expL = evolute(H[t], Liouville_operator, γ, times, times[t]) + ρt= expL * ρt + ∂ρt_∂x= [-im * Δt * liouville_commu(∂H_∂x[i]) * ρt for i in 1:para_num] + [expL] .* ∂ρt_∂x + end + CFIM(ρt|> vec2mat, ∂ρt_∂x|> vec2mat, M) +end + + +function SLD(ρ::Matrix{T}, ∂ρ_∂x::Matrix{T}) where {T <: Complex} + 2 * pinv(kron(ρ |> transpose, ρ |> one) + kron(ρ |> one, ρ)) * vec(∂ρ_∂x) |> vec2mat +end +function SLD(ρ::Vector{T},∂ρ_∂x::Vector{T}) where {T <: Complex} + SLD(ρ |> vec2mat, ∂ρ_∂x |> vec2mat) +end +function SLD(ρ::Matrix{T}, ∂ρ_∂x::Vector{Matrix{T}}) where {T <: Complex} + (x->SLD(ρ, x)).(∂ρ_∂x) +end +function SLD_qr(ρ::Matrix{T}, ∂ρ_∂x::Matrix{T}) where {T <: Complex} + 2 * (qr(kron(ρ |> transpose, ρ |> one) + kron(ρ |> one, ρ), Val(true)) \ vec(∂ρ_∂x)) |> vec2mat +end +function SLD_eig(ρ::Array{T}, dρ::Array{T})::Array{T} where {T <: Complex} + dim = size(ρ)[1] + if typeof(dρ) == Array{T,2} + purity = tr(ρ * ρ) + SLD_res = zeros(T, dim, dim) + if abs(1 - purity) < 1e-8 + SLD_res = 2 * dρ + else + val, vec_mat = eigen(ρ) + for fi in 1:dim + for fj in 1:dim + coeff = 2 / (val[fi] + val[fj]) + SLD_res[fi, fj] = coeff * (vec_mat[:,fi]' * (dρ * vec_mat[:, fj])) + end + end + SLD_res[findall(SLD_res == Inf)] .= 0. + SLD_res = vec_mat * (SLD_res * vec_mat') + end + else + # multi-parameter scenario + purity = tr(ρ * ρ) + if abs(1 - purity) < 1e-8 + SLD_res = [2 * dρ[i] for i in 1:length(dρ)] + else + # SLD_res = [zeros(T,dim,dim) for i in 1:length(dρ)] + dim = ndims(ρ) + val, vec_mat = eigens(ρ) + for para_i in 1:length(dρ) + SLD_tp = zeros(T, dim, dim) + for fi in 1:dim + for fj in 1:dim + coeff = 2. / (val[fi] + val[fj]) + SLD_tp[fi][fj] = coeff * (vec[fi]' * (dρ[para_i] * vec[fj])) + end + end + SLD_tp[findall(SLD_rp == Inf)] .= 0. + SLD_res[para_i] = vec_mat * (SLD_tp * vec_mat') + end + end + end + SLD_res +end +function RLD(ρ::Matrix{T}, dρ::Matrix{T}) where {T <: Complex} + dρ * pinv(ρ) +end +function QFI_RLD(ρ, dρ) + RLD_tp = RLD(ρ, dρ) + F = tr(ρ * RLD_tp * RLD_tp') + F |> real +end +function QFI(ρ, dρ) + SLD_tp = SLD(ρ, dρ) + SLD2_tp = SLD_tp * SLD_tp + F = tr(ρ * SLD2_tp) + F |> real +end +function QFIM(ρ, dρ) + SLD_tp = SLD(ρ, dρ) + [0.5*ρ] .* (kron(SLD_tp, SLD_tp') + kron(SLD_tp', SLD_tp)).|> tr .|>real +end +function QFI(H::Vector{Matrix{T}}, ∂H_∂x::Matrix{T}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, times) where {T <: Complex,R <: Real} + dim = size(H[1])[1] + Δt = times[2] - times[1] + ρt = evolute(H[1], Liouville_operator, γ, times, times[1]) * (ρ_initial |> vec) + ∂ρt_∂x = -im * Δt * liouville_commu(∂H_∂x) * ρt + for t in 2:length(times) + expL = evolute(H[t], Liouville_operator, γ, times, times[t]) + ρt= expL * ρt + ∂ρt_∂x= -im * Δt * liouville_commu(∂H_∂x) * ρt + expL * ∂ρt_∂x + end + QFI(ρt|> vec2mat, ∂ρt_∂x|> vec2mat) +end +function QFIM(H::Vector{Matrix{T}}, ∂H_∂x::Vector{Matrix{T}}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, times) where {T <: Complex,R <: Real} + dim = size(H[1])[1] + Δt = times[2] - times[1] + para_num = length(∂H_∂x) + ρt = evolute(H[1], Liouville_operator, γ, times, times[1]) * ρ_initial[:] + ∂ρt_∂x = [-im * Δt * liouville_commu(∂H_∂x[i]) * ρt for i in 1:para_num] + for t in 2:length(times) + expL = evolute(H[t], Liouville_operator, γ, times, times[t]) + ρt= expL * ρt + ∂ρt_∂x= [-im * Δt * liouville_commu(∂H_∂x[i]) * ρt for i in 1:para_num] + [expL] .* ∂ρt_∂x + end + QFIM(ρt|> vec2mat, ∂ρt_∂x|> vec2mat) +end +function CFI(M, system) + CFI(M,Htot(system.freeHamiltonian, system.control_Hamiltonian, system.control_coefficients), system.Hamiltonian_derivative[1], system.ρ_initial, system.Liouville_operator, system.γ, system.times) +end +function CFIM(M, system) + CFIM(M,Htot(system.freeHamiltonian, system.control_Hamiltonian, system.control_coefficients), system.Hamiltonian_derivative, system.ρ_initial, system.Liouville_operator, system.γ, system.times) +end +function QFI(system) + QFI(Htot(system.freeHamiltonian, system.control_Hamiltonian, system.control_coefficients), system.Hamiltonian_derivative[1], system.ρ_initial, system.Liouville_operator, system.γ, system.times) +end +function QFIM(system) + QFIM(Htot(system.freeHamiltonian, system.control_Hamiltonian, system.control_coefficients), system.Hamiltonian_derivative, system.ρ_initial, system.Liouville_operator, system.γ, system.times) +end \ No newline at end of file diff --git a/src/Common/common.jl b/src/Common/common.jl new file mode 100644 index 0000000..c2391cb --- /dev/null +++ b/src/Common/common.jl @@ -0,0 +1,104 @@ +function RunMixed( grape, particle_num= 10, c0=0.5, c1=0.5, c2=0.5, prec_rough = 0.1, sd=1234, episode=400) + println("Combined strategies") + println("searching initial controls with particle swarm optimization ") + tnum = length(grape.times) + ctrl_num = length(grape.control_Hamiltonian) + particles, velocity = repeat(grape, particle_num), [[10*ones(tnum) for i in 1:ctrl_num] for j in 1:particle_num] + pbest = [[zeros(tnum) for i in 1:ctrl_num] for j in 1:particle_num] + gbest = [zeros(tnum) for i in 1:ctrl_num] + velocity_best = [zeros(Float64, tnum) for i in 1:ctrl_num] + p_fit = zeros(particle_num) + qfi_ini = QFI(grape) + println("initial QFI is $(qfi_ini)") + for ei in 1:episode + fit_pre = 0.0 + fit = 0.0 + for pi in 1:particle_num + propagate!(particles[pi]) + f_now = QFI(particles[pi]) + + if f_now > p_fit[pi] + p_fit[pi] = f_now + for di in 1:ctrl_num + for ni in 1:tnum + pbest[pi][di][ni] = particles[pi].control_coefficients[di][ni] + end + end + end + end + for pj in 1:particle_num + if p_fit[pj] > fit + fit = p_fit[pj] + for dj in 1:ctrl_num + for nj in 1:tnum + gbest[dj][nj] = particles[pj].control_coefficients[dj][nj] + velocity_best[dj][nj] = velocity[pj][dj][nj] + end + end + end + end + Random.seed!(sd) + for pk in 1:particle_num + for dk in 1:ctrl_num + for ck in 1:tnum + velocity[pk][dk][ck] = c0*velocity[pk][dk][ck] + c1*rand()*(pbest[pk][dk][ck] - particles[pk].control_coefficients[dk][ck]) + + c2*rand()*(gbest[dk][ck] - particles[pk].control_coefficients[dk][ck]) + particles[pk].control_coefficients[dk][ck] = particles[pk].control_coefficients[dk][ck] + velocity[pk][dk][ck] + end + end + end + fit_pre = fit + if abs(fit-fit_pre) < prec_rough + grape.control_coefficients = gbest + println("PSO strategy finished, switching to GRAPE") + Run(grape) + return nothing + end + print("current QFI is $fit ($ei epochs) \r") + end +end + +function RunODE(grape) + println("quantum parameter estimation") + if length(grape.Hamiltonian_derivative) == 1 + println("single parameter estimation scenario") + qfi_ini = QFI(grape) + qfi_list = [qfi_ini] + println("initial QFI is $(qfi_ini)") + gradient_QFIM_ODE!(grape) + while true + qfi_now = QFI(grape) + gradient_QFIM_ODE!(grape) + if 0 < (qfi_now - qfi_ini) < 1e-4 + println("\n Iteration over, data saved.") + println("Final QFI is ", qfi_now) + save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "qfi", qfi_list) + break + else + qfi_ini = qfi_now + append!(qfi_list,qfi_now) + print("current QFI is ", qfi_now, " ($(qfi_list|>length) epochs) \r") + end + end + else + println("multiple parameters estimation scenario") + f_ini =1/(grape |> QFIM |> inv |> tr) + f_list = [f_ini] + println("initial 1/tr(F^-1) is $(f_ini)") + gradient_QFIM!(grape) + while true + f_now = 1/(grape |> QFIM |> inv |> tr) + gradient_QFIM!(grape) + if 0< f_now - f_ini < 1e-4 + println("\n Iteration over, data saved.") + println("Final 1/tr(F^-1) is ", f_now) + save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "f", f_list) + break + else + f_ini = f_now + append!(f_list,f_now) + print("current 1/tr(F^-1) is ", f_now, " ($(f_list|>length) epochs) \r") + end + end + end +end diff --git a/src/Common/utils.jl b/src/Common/utils.jl new file mode 100644 index 0000000..4b7630c --- /dev/null +++ b/src/Common/utils.jl @@ -0,0 +1,77 @@ +using SparseArrays +""" +Pauli matrces +""" +sigmax() = [.0im 1.;1. 0.] +sigmay() = [0. -1.0im;1.0im 0.] +sigmaz() = [1.0 .0im;0. -1.] +sigmap() = [.0im 1.;0. 0.] +sigmam() = [.0im 0.;1. 0.] +sigmax(i, N) = kron(I(2^(i - 1)), sigmax(), I(2^(N - i))) +sigmay(i, N) = kron(I(2^(i - 1)), sigmay(), I(2^(N - i))) +sigmaz(i, N) = kron(I(2^(i - 1)), sigmaz(), I(2^(N - i))) +sigmap(i, N) = kron(I(2^(i - 1)), sigmap(), I(2^(N - i))) +sigmam(i, N) = kron(I(2^(i - 1)), sigmam(), I(2^(N - i))) + +""" + destroy() + +""" +function destroy(M) + spdiagm(M, M, 1 => map(x -> x |> sqrt, 1:(M - 1))) +end + +""" + vec2mat() +transform a vector to its matrix form +""" +function vec2mat(x::Vector{T}) where {T <: Number} + reshape(x, x |> length |> sqrt |> Int, :) +end +function vec2mat(x) + vec2mat.(x) +end +function vec2mat(x::Matrix) + throw(ErrorException("vec2mating a matrix of size $(size(x))")) +end + +""" + repeat() + +""" +function Base.repeat(system, N) + [system for i in 1:N] +end +function Base.repeat(system, M, N) + reshape(repeat(system, M * N), M, N) +end + +""" + filterZeros!() + +""" +function filterZeros!(x::Matrix{T}) where {T <: Complex} + x[abs.(x) .< eps()] .= zero(T) + x +end +function filterZeros!(x) + filterZeros!.(x) +end + +""" + t2Num() + +""" +function t2Num(t0, dt, t) + Int(round((t - t0) / dt)) + 1 +end + +""" + basis() + +""" +function basis(dim, si, ::T)::Array{T} where {T <: Complex} + result = zeros(T, dim) + result[si] = 1.0 + result +end \ No newline at end of file diff --git a/src/Control/GRAPE.jl b/src/Control/GRAPE.jl new file mode 100644 index 0000000..6f5ef50 --- /dev/null +++ b/src/Control/GRAPE.jl @@ -0,0 +1,280 @@ +abstract type ControlSystem end +mutable struct GrapeControl{T <: Complex,M <: Real} <: ControlSystem + freeHamiltonian::Matrix{T} + Hamiltonian_derivative::Vector{Matrix{T}} + ρ_initial::Matrix{T} + times::StepRangeLen{M,Base.TwicePrecision{M},Base.TwicePrecision{M}} + Liouville_operator::Vector{Matrix{T}} + γ::Vector{M} + control_Hamiltonian::Vector{Matrix{T}} + control_coefficients::Vector{Vector{M}} + ϵ::M + ρ::Vector{Matrix{T}} + ∂ρ_∂x::Vector{Vector{Matrix{T}}} + GrapeControl(freeHamiltonian::Matrix{T}, Hamiltonian_derivative::Vector{Matrix{T}}, ρ_initial::Matrix{T}, + times::StepRangeLen{M,Base.TwicePrecision{M},Base.TwicePrecision{M}}, + Liouville_operator::Vector{Matrix{T}},γ::Vector{M}, control_Hamiltonian::Vector{Matrix{T}}, + control_coefficients::Vector{Vector{M}}, ϵ=0.01, ρ=Vector{Matrix{T}}(undef, 1), + ∂ρ_∂x=Vector{Vector{Matrix{T}}}(undef, 1),∂ρ_∂V=Vector{Vector{Matrix{T}}}(undef, 1)) where {T <: Complex,M <: Real} = + new{T,M}(freeHamiltonian, Hamiltonian_derivative, ρ_initial, times, Liouville_operator, γ, control_Hamiltonian, + control_coefficients, ϵ, ρ, ∂ρ_∂x) +end + +function gradient_CFI!(grape::GrapeControl{T}, M) where {T <: Complex} + δI = gradient(x->CFI(M, Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients)[1].|>real + grape.control_coefficients += grape.ϵ*δI +end +function gradient_CFIM!(grape::GrapeControl{T}, M) where {T <: Complex} + δI = gradient(x->CFIM(M, Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients).|>real + grape.control_coefficients += grape.ϵ*δI +end +function gradient_CFI_ADAM!(grape::GrapeControl{T}, M) where {T <: Complex} + δI = gradient(x->CFI(M, Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients)[1].|>real + Adam!(grape, δI) +end +function gradient_CFIM_ADAM!(grape::GrapeControl{T}, M, mt, vt) where {T <: Complex} + δI = gradient(x->CFIM(M, Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients).|>real + Adam!(grape, δI) +end + +function gradient_QFI!(grape::GrapeControl{T}) where {T <: Complex} + δF = gradient(x->QFI(Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients)[1].|>real + grape.control_coefficients += grape.ϵ*δF +end +function gradient_QFIM!(grape::GrapeControl{T}) where {T <: Complex} + δF = gradient(x->1/(QFIM(Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative, grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times) |> pinv |> tr |>real), grape.control_coefficients).|>real |>sum + grape.control_coefficients += grape.ϵ*δF +end +function gradient_QFI_ADAM!(grape::GrapeControl{T}) where {T <: Complex} + δF = gradient(x->QFI(Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients)[1].|>real + Adam!(grape, δF) +end +function gradient_QFIM_ADAM!(grape::GrapeControl{T}) where {T <: Complex} + δF = gradient(x->1/(QFIM(Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative, grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times) |> pinv |> tr |>real), grape.control_coefficients).|>real |>sum + Adam!(grape, δF) +end +function gradient_QFI_analitical_ADAM!(grape::GrapeControl{T}) where {T <: Complex} + δF = propagate_analitical!(grape) + Adam!(grape, δF) +end +function gradient_QFIM_ODE(grape::GrapeControl) + H = Htot(grape.freeHamiltonian, grape.control_Hamiltonian, grape.control_coefficients) + Δt = grape.times[2] - grape.times[1] + t_num = length(grape.times) + para_num = length(grape.Hamiltonian_derivative) + ctrl_num = length(grape.control_Hamiltonian) + tspan(j) = (grape.times[1], grape.times[j]) + tspan() = (grape.times[1], grape.times[end]) + u0 = grape.ρ_initial |> vec + evo(p, t) = evolute(p[t2Num(tspan()[1], Δt, t)], grape.Liouville_operator, grape.γ, grape.times, t2Num(tspan()[1], Δt, t)) + f(u, p, t) = evo(p, t) * u + prob = DiscreteProblem(f, u0, tspan(), H,dt=Δt) + ρt = solve(prob).u + ∂ρt_∂x = Vector{Vector{Vector{eltype(u0)}}}(undef, 1) + for para in 1:para_num + devo(p, t) = -1.0im * Δt * liouville_commu(grape.Hamiltonian_derivative[para]) * evo(p, t) + du0 = devo(H, tspan()[1]) * u0 + g(du, p, t) = evo(p, t) * du + devo(p, t) * ρt[t2Num(tspan()[1], Δt, t)] + dprob = DiscreteProblem(g, du0, tspan(), H,dt=Δt) + ∂ρt_∂x[para] = solve(dprob).u + end + δρt_δV = Matrix{Vector{Vector{eltype(u0)}}}(undef,ctrl_num,length(grape.times)) + for ctrl in 1:ctrl_num + for j in 1:t_num + devo(p, t) = -1.0im * Δt * liouville_commu(grape.control_Hamiltonian[ctrl]) * evo(p, t) + du0 = devo(H, tspan()[1]) * u0 + g(du, p, t) = evo(p, t) * du + devo(p, t) * ρt[t2Num(tspan()[1], Δt, t)] + dprob = DiscreteProblem(g, du0, tspan(j), H,dt=Δt) + δρt_δV[ctrl,j] = solve(dprob).u + end + end + ∂xδρt_δV = Array{Vector{eltype(u0)}, 3}(undef,para_num, ctrl_num,length(grape.times)) + for para in 1:para_num + for ctrl in 1:ctrl_num + dxevo = -1.0im * Δt * liouville_commu(grape.Hamiltonian_derivative[para]) + dkevo = -1.0im * Δt * liouville_commu(grape.control_Hamiltonian[ctrl]) + for j in 1:t_num + g(du, p, t) = dxevo * dkevo * evo(p, t) * ρt[t2Num(tspan()[1], Δt, t)] + + dxevo * evo(p, t) * δρt_δV[ctrl, j][t2Num(tspan()[1], Δt, t)] + + dkevo * evo(p, t) * ∂ρt_∂x[para][t2Num(tspan()[1], Δt, t)] + + evo(p, t) * du + du0 = dxevo * dkevo * evo(H,tspan()[1]) * ρt[t2Num(tspan()[1], Δt, tspan()[1])] + dprob = DiscreteProblem(g, du0, tspan(j), H, dt=Δt) + ∂xδρt_δV[para, ctrl, j] = solve(dprob).u[end] + end + end + end + δF = grape.control_coefficients .|> zero + for para in 1:para_num + SLD_tp = SLD(ρt[end], ∂ρt_∂x[para][end]) + for ctrl in 1:ctrl_num + for j in 1:t_num + δF[ctrl][j] -= 2 * tr((∂xδρt_δV[para,ctrl,j]|> vec2mat) * SLD_tp) - + tr((δρt_δV[ctrl, j][end] |> vec2mat) * SLD_tp^2) |> real + end + end + end + δF +end +function gradient_QFIM_ODE!(grape::GrapeControl{T}) where {T <: Complex} + grape.control_coefficients += grape.ϵ * gradient_QFIM_ODE(grape) +end + +function Run(M, grape::GrapeControl{T}) where {T<: Complex} + println("classical parameter estimation") + if length(grape.Hamiltonian_derivative) == 1 + println("single parameter estimation scenario") + cfi_ini = CFI(M, grape) + cfi_list = [cfi_ini] + println("initial CFI is $(cfi_ini)") + gradient_CFI!(M, grape) + while true + cfi_now = CFI(M, grape) + gradient_CFI!(M, grape) + if 0 < (cfi_now - cfi_ini) < 1e-4 + println("\n Iteration over, data saved.") + println("Final CFI is ", cfi_now) + save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "cfi", cfi_list) + break + else + cfi_ini = cfi_now + append!(cfi_list,cfi_now) + print("current CFI is ", cfi_now, " ($(cfi_list|>length) epochs) \r") + end + end + else + println("multiple parameters estimation scenario") + f_ini =1/(grape |> CFIM |> inv |> tr) + f_list = [f_ini] + println("initial 1/tr(F^-1) is $(f_ini)") + gradient_CFIM!(M, grape) + while true + f_now = 1/(grape |> CFIM |> inv |> tr) + gradient_CFIM!(M, grape) + if 0< f_now - f_ini < 1e-4 + println("\n Iteration over, data saved.") + println("Final 1/tr(I^-1) is ", f_now) + save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "f", f_list) + break + else + f_ini = f_now + append!(f_list,f_now) + print("current 1/tr(I^-1) is ", f_now, " ($(f_list|>length) epochs) \r") + end + end + end +end +function RunADAM(grape) + println("AutoGrape strategies") + println("quantum parameter estimation") + if length(grape.Hamiltonian_derivative) == 1 + println("single parameter estimation scenario") + qfi_ini = QFI(grape) + qfi_list = [qfi_ini] + println("initial QFI is $(qfi_ini)") + gradient_QFI_ADAM!(grape) + while true + qfi_now = QFI(grape) + gradient_QFI_ADAM!(grape) + if 0 < (qfi_now - qfi_ini) < 1e-4 + println("\n Iteration over, data saved.") + println("Final QFI is ", qfi_now) + save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "qfi", qfi_list) + break + else + qfi_ini = qfi_now + append!(qfi_list,qfi_now) + print("current QFI is ", qfi_now, " ($(qfi_list|>length) epochs) \r") + end + end + else + println("multiple parameters estimation scenario") + f_ini =1/(grape |> QFIM |> inv |> tr) + f_list = [f_ini] + println("initial 1/tr(F^-1) is $(f_ini)") + gradient_QFIM_ADAM!(grape) + while true + f_now = 1/(grape |> QFIM |> inv |> tr) + gradient_QFIM_ADAM!(grape) + if 0< f_now - f_ini < 1e-4 + println("\n Iteration over, data saved.") + println("Final 1/tr(F^-1) is ", f_now) + save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "f", f_list) + break + else + f_ini = f_now + append!(f_list,f_now) + print("current 1/tr(F^-1) is ", f_now, " ($(f_list|>length) epochs) \r") + end + end + end +end +function RunAnaliticalADAM(grape) + println("Analitical strategies") + println("quantum parameter estimation") + if length(grape.Hamiltonian_derivative) == 1 + println("single parameter estimation scenario") + qfi_ini = QFI(grape) + qfi_list = [qfi_ini] + println("initial QFI is $(qfi_ini)") + gradient_QFI_analitical_ADAM!(grape) + while true + qfi_now = QFI(grape) + gradient_QFI_analitical_ADAM!(grape) + if 0 < (qfi_now - qfi_ini) < 1e-4 + println("\n Iteration over, data saved.") + println("Final QFI is ", qfi_now) + save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "qfi", qfi_list) + break + else + qfi_ini = qfi_now + append!(qfi_list,qfi_now) + print("current QFI is ", qfi_now, " ($(qfi_list|>length) epochs) \r") + end + end + end +end +function RunADAM(M, grape::GrapeControl{T}) where {T<: Complex} + println("classical parameter estimation") + if length(grape.Hamiltonian_derivative) == 1 + println("single parameter estimation scenario") + cfi_ini = CFI(M, grape) + cfi_list = [cfi_ini] + println("initial CFI is $(cfi_ini)") + gradient_CFI_ADAM!(M, grape) + while true + cfi_now = CFI(M, grape) + gradient_CFI_ADAM!(M, grape) + if 0 < (cfi_now - cfi_ini) < 1e-4 + println("\n Iteration over, data saved.") + println("Final CFI is ", cfi_now) + save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "cfi", cfi_list) + break + else + cfi_ini = cfi_now + append!(cfi_list,cfi_now) + print("current CFI is ", cfi_now, " ($(cfi_list|>length) epochs) \r") + end + end + else + println("multiple parameters estimation scenario") + f_ini =1/(grape |> CFIM |> inv |> tr) + f_list = [f_ini] + println("initial 1/tr(F^-1) is $(f_ini)") + gradient_CFIM!(M, grape) + while true + f_now = 1/(grape |> CFIM |> inv |> tr) + gradient_CFIM!(M, grape) + if 0< f_now - f_ini < 1e-4 + println("\n Iteration over, data saved.") + println("Final 1/tr(I^-1) is ", f_now) + save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "f", f_list) + break + else + f_ini = f_now + append!(f_list,f_now) + print("current 1/tr(I^-1) is ", f_now, " ($(f_list|>length) epochs) \r") + end + end + end +end \ No newline at end of file diff --git a/src/Control/PSO.jl b/src/Control/PSO.jl new file mode 100644 index 0000000..3da14ff --- /dev/null +++ b/src/Control/PSO.jl @@ -0,0 +1,118 @@ +function RunPSO(grape::GrapeControl{T}; particle_num= 10, c0=0.5, c1=0.5, c2=0.5,v0=0.1, sd=1234, episode=400) where {T<: Complex} + println("PSO strategies") + println("searching optimal controls with particle swarm optimization ") + tnum = length(grape.times) + ctrl_num = length(grape.control_Hamiltonian) + particles = repeat(grape, particle_num) + velocity = v0.*rand(ctrl_num, tnum, particle_num)|>SharedArray + pbest = zeros(ctrl_num, tnum, particle_num) + gbest = zeros(ctrl_num, tnum) |> SharedArray + velocity_best = zeros(ctrl_num,tnum) + p_fit = zeros(particle_num) + qfi_ini = QFI(grape) + println("initial QFI is $(qfi_ini)") + fit_pre = 0.0 + fit = 0.0 + for ei in 1:episode + @inbounds @threads for pi in 1:particle_num + propagate!(particles[pi]) + f_now = QFI(particles[pi]) + + if f_now > p_fit[pi] + p_fit[pi] = f_now + for di in 1:ctrl_num + for ni in 1:tnum + pbest[di,ni,pi] = particles[pi].control_coefficients[di][ni] + end + end + end + end + @inbounds @threads for pj in 1:particle_num + if p_fit[pj] > fit + fit = p_fit[pj] + for dj in 1:ctrl_num + @inbounds @threads for nj in 1:tnum + gbest[dj, nj] = particles[pj].control_coefficients[dj][nj] + velocity_best[dj, nj] = velocity[dj, nj, pj] + end + end + end + end + Random.seed!(sd) + @inbounds @threads for pk in 1:particle_num + for dk in 1:ctrl_num + @inbounds @threads for ck in 1:tnum + velocity[dk, ck, pk] = c0*velocity[dk, ck, pk] + c1*rand()*(pbest[dk, ck, pk] - particles[pk].control_coefficients[dk][ck]) + + c2*rand()*(gbest[dk, ck] - particles[pk].control_coefficients[dk][ck]) + particles[pk].control_coefficients[dk][ck] = particles[pk].control_coefficients[dk][ck] + velocity[dk, ck, pk] + end + end + end + # println(particles[1].control_coefficients[1][1]) + if ei == episode#abs(fit-fit_pre) < 1e-5 + grape.control_coefficients = gbest + println("Final QFI is ", fit) + return nothing + end + fit_pre = fit + print("current QFI is $fit ($ei epochs) \r") + end +end +function RunPSOAdam(grape::GrapeControl{T}; particle_num= 10, c0=0.5, c1=0.5, c2=0.5,v0=0.01, sd=1234, episode=400) where {T<: Complex} + println("PSO strategies") + println("searching optimal controls with particle swarm optimization ") + tnum = length(grape.times) + ctrl_num = length(grape.control_Hamiltonian) + particles, velocity = repeat(grape, particle_num), [[v0*ones(tnum) for i in 1:ctrl_num] for j in 1:particle_num] + pbest = [[zeros(tnum) for i in 1:ctrl_num] for j in 1:particle_num] + gbest = [zeros(tnum) for i in 1:ctrl_num] + velocity_best = [zeros(Float64, tnum) for i in 1:ctrl_num] + p_fit = zeros(particle_num) + qfi_ini = QFI(grape) + println("initial QFI is $(qfi_ini)") + fit_pre = 0.0 + fit = 0.0 + for ei in 1:episode + for pi in 1:particle_num + propagate!(particles[pi]) + f_now = QFI(particles[pi]) + + if f_now > p_fit[pi] + p_fit[pi] = f_now + for di in 1:ctrl_num + for ni in 1:tnum + pbest[pi][di][ni] = particles[pi].control_coefficients[di][ni] + end + end + end + end + for pj in 1:particle_num + if p_fit[pj] > fit + fit = p_fit[pj] + for dj in 1:ctrl_num + for nj in 1:tnum + gbest[dj][nj] = particles[pj].control_coefficients[dj][nj] + velocity_best[dj][nj] = velocity[pj][dj][nj] + end + end + end + end + Random.seed!(sd) + for pk in 1:particle_num + for dk in 1:ctrl_num + for ck in 1:tnum + velocity[pk][dk][ck] = c0*velocity[pk][dk][ck] + c1*(rand()-0.5)*(pbest[pk][dk][ck] - particles[pk].control_coefficients[dk][ck]) + + c2*(rand()-0.5)*(gbest[dk][ck] - particles[pk].control_coefficients[dk][ck]) + end + end + Adam!(particles[pk], velocity[pk]) + end + if abs(fit-fit_pre) < 1e-5 + grape.control_coefficients = gbest + println("Final QFI is ", fit) + return nothing + end + fit_pre = fit + print("current QFI is $fit ($ei epochs) \r") + end +end \ No newline at end of file diff --git a/src/Dynamics/Adam.jl b/src/Dynamics/Adam.jl new file mode 100644 index 0000000..a396344 --- /dev/null +++ b/src/Dynamics/Adam.jl @@ -0,0 +1,20 @@ +""" + Adam() + +""" +function Adam(gt, t, para, m_t, v_t, alpha=0.01, beta1=0.90, beta2=0.99, epsilon=1e-8) + t = t+1 + m_t = beta1*m_t + (1-beta1)*gt + v_t = beta2*v_t + (1-beta2)*(gt*gt) + m_cap = m_t/(1-(beta1^t)) + v_cap = v_t/(1-(beta2^t)) + para = para+(alpha*m_cap)/(sqrt(v_cap)+epsilon) + return para, m_t, v_t +end +function Adam!(grape, δ, mt=0.0, vt=0.0) + for ctrl in 1:length(δ) + for ti in 1:length(grape.times) + grape.control_coefficients[ctrl][ti], mt, vt = Adam(δ[ctrl][ti], ti, grape.control_coefficients[ctrl][ti], mt, vt) + end + end +end \ No newline at end of file diff --git a/src/Dynamics/dynamcs.jl b/src/Dynamics/dynamcs.jl new file mode 100644 index 0000000..15dcd6b --- /dev/null +++ b/src/Dynamics/dynamcs.jl @@ -0,0 +1,279 @@ + +""" + liouville_commu() + +""" +function liouville_commu(H) + kron(one(H), H) - kron(H |> transpose, one(H)) +end + +""" + liouville_commu_dissip() + +""" +function liouville_dissip(Γ) + kron(Γ |> conj, Γ) - 0.5 * kron((Γ |> transpose) * Γ |> conj, Γ |> one) - 0.5 * kron(Γ |> one, Γ' * Γ) +end + +""" + dissipation() + +""" +function dissipation(Γ::Vector{Matrix{T}}, γ::Vector{R}, t::Real) where {T <: Complex,R <: Real} + [γ[i] * liouville_dissip(Γ[i]) for i in 1:length(Γ)] |> sum +end +function dissipation(Γ::Vector{Matrix{T}}, γ::Vector{Vector{R}}, t::Real) where {T <: Complex,R <: Real} + [γ[i][t] * liouville_dissip(Γ[i]) for i in 1:length(Γ)] |> sum +end + +""" + free_evolution() + +""" +function free_evolution(H0) + -1.0im * liouville_commu(H0) +end + +""" + liouvillian() + +""" +function liouvillian(H::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, t::Real) where {T <: Complex} + freepart = liouville_commu(H) + dissp = norm(γ) +1 ≈ 1 ? freepart|>zero : dissipation(Liouville_operator, γ, t) + -1.0im * freepart + dissp +end + +""" + Htot() +""" +function Htot(H0::Matrix{T}, control_Hamiltonian::Vector{Matrix{T}}, control_coefficients) where {T <: Complex} + Htot = [H0] .+ ([control_coefficients[i] .* [control_Hamiltonian[i]] for i in 1:length(control_coefficients)] |> sum ) +end + +""" + evolute() + +""" +function evolute(H, Liouville_operator, γ, times, t) + tj = Int(round((t - times[1]) / (times[2] - times[1]))) + 1 + dt = times[2] - times[1] + Ld = dt * liouvillian(H, Liouville_operator, γ, tj) + exp(Ld) +end + +""" + evolute_ODE!() + +""" +function evolute_ODE!(grape::GrapeControl) + H(p) = Htot(grape.freeHamiltonian, grape.control_Hamiltonian, p) + dt = grape.times[2] - grape.times[1] + tspan = (grape.times[1], grape.times[end]) + u0 = grape.ρ_initial + Γ = grape.Liouville_operator + f(u, p, t) = -im * (H(p)[t2Num(tspan[1], dt, t)] * u + u * H(p)[t2Num(tspan[1], dt, t)]) + + ([grape.γ[i] * (Γ[i] * u * Γ[i]' - (Γ[i]' * Γ[i] * u + u * Γ[i]' * Γ[i] )) for i in 1:length(Γ)] |> sum) + prob = ODEProblem(f, u0, tspan, grape.control_coefficients, saveat=dt) + sol = solve(prob) + sol.u +end + +""" + propagate() + +""" +function propagate(H0::Matrix{T}, ∂H_∂x::Vector{Matrix{T}}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, + γ, control_Hamiltonian::Vector{Matrix{T}}, control_coefficients::Vector{Vector{R}}, times) where {T <: Complex,R <: Real} + dim = size(H0)[1] + para_num = length(∂H_∂x) + H = Htot(H0, control_Hamiltonian, control_coefficients) + ρt = [Vector{ComplexF64}(undef, dim^2) for i in 1:length(times)] + ∂ρt_∂x = [[Vector{ComplexF64}(undef, dim^2) for i in 1:length(times)] for para in 1:para_num] + Δt = times[2] - times[1] + ρt[1] = evolute(H[1], Liouville_operator, γ, times, times[1]) * (ρ_initial |> vec) + for para in 1:para_num + ∂ρt_∂x[para][1] = -im * Δt * liouville_commu(∂H_∂x[para]) * ρt[1] + end + for t in 2:length(times) + expL = evolute(H[t], Liouville_operator, γ, times, times[t]) + ρt[t] = expL * ρt[t - 1] + for para in para_num + ∂ρt_∂x[para][t] = -im * Δt * liouville_commu(∂H_∂x[para]) * ρt[t] + expL * ∂ρt_∂x[para][t - 1] + end + end + ρt .|> vec2mat, ∂ρt_∂x .|> vec2mat +end + +""" + propagate!() + +""" +function propagate!(grape::GrapeControl) + grape.ρ, grape.∂ρ_∂x = propagate(grape.freeHamiltonian, grape.Hamiltonian_derivative, grape.ρ_initial, + grape.Liouville_operator, grape.γ, grape.control_Hamiltonian, + grape.control_coefficients, grape.times ) +end + +""" + propagate_analytical() +""" +function propagate_analytical(H0::Matrix{T}, ∂H_∂x::Vector{Matrix{T}}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, + γ, control_Hamiltonian::Vector{Matrix{T}}, control_coefficients::Vector{Vector{R}}, times) where {T <: Complex,R <: Real} + + dim = size(H0)[1] + tnum = length(times) + para_num = length(∂H_∂x) + ctrl_num = length(control_Hamiltonian) + Δt = times[2] - times[1] + + H = Htot(H0, control_Hamiltonian, control_coefficients) + + ρt = [Vector{ComplexF64}(undef, dim^2) for i in 1:tnum] + ∂ρt_∂x = [[Vector{ComplexF64}(undef, dim^2) for para in 1:para_num] for i in 1:tnum] + δρt_δV = [[] for ctrl in 1:ctrl_num] + ∂xδρt_δV = [[[] for ctrl in 1:ctrl_num] for i in 1:para_num] + ∂H_L = [Matrix{ComplexF64}(undef, dim^2,dim^2) for i in 1:para_num] + Hc_L = [Matrix{ComplexF64}(undef, dim^2,dim^2) for i in 1:ctrl_num] + + ρt[1] = evolute(H[1], Liouville_operator, γ, times, times[1]) * (ρ_initial |> vec) + for pi in 1:para_num + ∂ρt_∂x[1][pi] = -im * Δt * ∂H_L[pi] * (ρ_initial |> vec) + ∂H_L[pi] = liouville_commu(∂H_∂x[pi]) + for ci in 1:ctrl_num + append!(δρt_δV[ci], [-im*Δt*Hc_L[ci]*ρt[1]]) + append!(∂xδρt_δV[pi][ci], [-im*Δt*Hc_L[ci]*∂ρt_∂x[1][pi]]) + end + end + + for cj in 1:ctrl_num + Hc_L[cj] = liouville_commu(control_Hamiltonian[cj]) + end + + for ti in 2:tnum + expL = evolute(H[ti], Liouville_operator, γ, times, times[ti]) + ρt[ti] = expL * ρt[ti-1] + for pk in 1:para_num + ∂ρt_∂x[ti][pk] = -im * Δt * ∂H_L[pk] * ρt[ti] + expL * ∂ρt_∂x[ti-1][pk] + for ck in 1:ctrl_num + for tk in 1:ti + δρt_δV_first = popfirst!(δρt_δV[ck]) + ∂xδρt_δV_first = popfirst!(∂xδρt_δV[pk][ck]) + δρt_δV_tp = expL * δρt_δV_first + ∂xδρt_δV_tp = -im * Δt * ∂H_L[pk] * expL * δρt_δV_first + expL * ∂xδρt_δV_first + append!(δρt_δV[ck], [δρt_δV_tp]) + append!(∂xδρt_δV[pk][ck], [∂xδρt_δV_tp]) + end + δρt_δV_last = -im * Δt * Hc_L[ck] * ρt[ti] + ∂xδρt_δV_last = -im * Δt * Hc_L[ck] * ∂ρt_∂x[ti][pk] + append!(δρt_δV[ck], [δρt_δV_last]) + append!(∂xδρt_δV[pk][ck], [∂xδρt_δV_last]) + end + end + end + + ρt_T = ρt[end] |> vec2mat + ∂ρt_T = [(∂ρt_∂x[end][para] |> vec2mat) for para in 1:para_num] + F_T = QFIM(ρt_T, ∂ρt_T) + + if para_num == 1 + Lx = SLD_eig(ρt_T, ∂ρt_T[1]) + anti_commu = 2*Lx[1]*Lx[1] + δF = [[0.0 for i in 1:tnum] for ctrl in 1:ctrl_num] + for tm in 1:tnum + for cm in 1:ctrl_num + ∂ρt_T_δV = δρt_δV[cm][tm] |> vec2mat + ∂xδρt_T_δV = ∂xδρt_δV[1][cm][tm] |> vec2mat + term1 = tr(∂xδρt_T_δV*Lx[1]) + term2 = tr(∂ρt_T_δV*anti_commu) + δF[cm][tm] = ((2*term1-0.5*term2) |> real) + end + end + + elseif para_num == 2 + F_det = F_T[1,1] * F_T[2,2] - F_T[1,2] * F_T[2,1] + δF = [[0.0 for i in 1:tnum] for ctrl in 1:ctrl_num] + for tm in 1:tnum + for cm in 1:ctrl_num + for pm in 1:para_num + ∂ρt_T_δV = δρt_δV[cm][tm] |> vec2mat + ∂xδρt_T_δV = ∂xδρt_δV[pm][cm][tm] |> vec2mat + term1 = tr(∂xδρt_T_δV * Lx[pm]) + anti_commu = 2 * Lx[pm] * Lx[pm] + term2 = tr(∂ρt_T_δV * anti_commu) + δF[cm][tm] = δF[cm][tm] - ((2*term1-0.5*term2) |> real) + end + δF[cm][tm] = δF[cm][tm] / F_det + end + end + + else + δF = [[0.0 for i in 1:tnum] for ctrl in 1:ctrl_num] + for tm in 1:tnum + for cm in 1:ctrl_num + for pm in 1:para_num + ∂ρt_T_δV = δρt_δV[cm][tm] |> vec2mat + ∂xδρt_T_δV = ∂xδρt_δV[pm][cm][tm] |> vec2mat + term1 = tr(∂xδρt_T_δV * Lx[pm]) + anti_commu = 2 * Lx[pm] * Lx[pm] + term2 = tr(∂ρt_T_δV * anti_commu) + δF[cm][tm] = δF[cm][tm] - ((2*term1-0.5*term2) |> real) / (F_T[pm][pm] * F_T[pm][pm]) + end + end + end + end + ρt |> vec2mat |> filterZeros!, ∂ρt_∂x |> vec2mat |> filterZeros!, δF +end + +""" + propagate_analitical!() +""" +function propagate_analitical!(grape::GrapeControl) + grape.ρ, grape.∂ρ_∂x, δF = propagate_analitical(grape.freeHamiltonian, grape.Hamiltonian_derivative, grape.ρ_initial, + grape.Liouville_operator, grape.γ, grape.control_Hamiltonian, + grape.control_coefficients, grape.times ) + δF +end + +""" + propagate_ODEAD!() +""" +function propagate_ODEAD!(grape::GrapeControl) + H(p) = Htot(grape.freeHamiltonian, grape.control_Hamiltonian, p) + dt = grape.times[2] - grape.times[1] + tspan = (grape.times[1], grape.times[end]) + u0 = grape.ρ_initial + Γ = grape.Liouville_operator + f(u, p, t) = -im * (H(p)[t2Num(tspan[1], dt, t)] * u + u * H(p)[t2Num(tspan[1], dt, t)]) + + ([grape.γ[i] * (Γ[i] * u * Γ[i]' - (Γ[i]' * Γ[i] * u + u * Γ[i]' * Γ[i] )) for i in 1:length(Γ)] |> sum) + p = grape.control_coefficients + prob = ODEProblem(f, u0, tspan, p, saveat=dt) + u = solve(prob).u + du = Zygote.jacobian(solve(remake(prob, u0=u, p), sensealg=QuadratureAdjoint())) + u, du +end + +""" + propagate_L_ODE!() +""" +function propagate_L_ODE!(grape::GrapeControl) + H = Htot(grape.freeHamiltonian, grape.control_Hamiltonian, grape.control_coefficients) + Δt = grape.times[2] - grape.times[1] + tspan = (grape.times[1], grape.times[end]) + u0 = grape.ρ_initial |> vec + evo(p, t) = evolute(p[t2Num(tspan[1], Δt, t)], grape.Liouville_operator, grape.γ, grape.times, t2Num(tspan[1], Δt, t)) + f(u, p, t) = evo(p, t) * u + prob = DiscreteProblem(f, u0, tspan, H,dt=Δt) + ρt = solve(prob).u + ∂ρt_∂x = Vector{Vector{Vector{eltype(u0)}}}(undef, 1) + for para in 1:length(grape.Hamiltonian_derivative) + devo(p, t) = -1.0im * Δt * liouville_commu(grape.Hamiltonian_derivative[para]) * evo(p, t) + du0 = devo(H, tspan[1]) * u0 + g(du, p, t) = evo(p, t) * du + devo(p, t) * ρt[t2Num(tspan[1], Δt, t)] + dprob = DiscreteProblem(g, du0, tspan, H,dt=Δt) + ∂ρt_∂x[para] = solve(dprob).u + end + + grape.ρ, grape.∂ρ_∂x = ρt |> vec2mat, ∂ρt_∂x |> vec2mat +end + diff --git a/src/QuanEstimation.jl b/src/QuanEstimation.jl new file mode 100644 index 0000000..50650fe --- /dev/null +++ b/src/QuanEstimation.jl @@ -0,0 +1,25 @@ +module QuanEstimation + +using LinearAlgebra +using Zygote +using DifferentialEquations +using JLD +using Random +using SharedArrays +using Base.Threads +using SparseArrays + +include("AsymptoticBound/CramerRao.jl") +include("Common/common.jl") +include("Common/utils.jl") +include("Control/GRAPE.jl") +include("Control/PSO.jl") +include("Dynamics/Adam.jl") +include("Dynamics/dynamcs.jl") +# include("QuanResources/") + +export sigmax, sigmay, sigmaz, sigmam +export GrapeControl, evolute, propagate!, QFI, CFI, gradient_CFI!,gradient_QFIM! +export Run, RunODE, RunMixed, RunADAM, RunPSO + +end diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..d251e80 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,6 @@ +using QuanEstimation +using Test + +@testset "QuanEstimation.jl" begin + # Write your tests here. +end From d044cae23dcdf3b7cb0d7a9574a16fdc5c06eaa0 Mon Sep 17 00:00:00 2001 From: HuaiMing Date: Tue, 8 Mar 2022 20:30:48 +0800 Subject: [PATCH 2/2] restruct --- .gitignore | 2 + src/AsymptoticBound/CramerRao.jl | 158 ----------------- src/Common/common.jl | 104 ------------ src/Common/utils.jl | 77 --------- src/Control/GRAPE.jl | 280 ------------------------------- src/Control/PSO.jl | 118 ------------- src/Dynamics/Adam.jl | 20 --- src/Dynamics/dynamcs.jl | 279 ------------------------------ src/QuanEstimation.jl | 24 +-- src/opt_target.jl | 20 +++ src/struct.jl | 12 ++ 11 files changed, 35 insertions(+), 1059 deletions(-) delete mode 100644 src/AsymptoticBound/CramerRao.jl delete mode 100644 src/Common/common.jl delete mode 100644 src/Common/utils.jl delete mode 100644 src/Control/GRAPE.jl delete mode 100644 src/Control/PSO.jl delete mode 100644 src/Dynamics/Adam.jl delete mode 100644 src/Dynamics/dynamcs.jl create mode 100644 src/opt_target.jl create mode 100644 src/struct.jl diff --git a/.gitignore b/.gitignore index eadbc49..85c3f68 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,5 @@ /Manifest.toml /deps/build.log /docs/build/ + +/.history diff --git a/src/AsymptoticBound/CramerRao.jl b/src/AsymptoticBound/CramerRao.jl deleted file mode 100644 index 64ec15e..0000000 --- a/src/AsymptoticBound/CramerRao.jl +++ /dev/null @@ -1,158 +0,0 @@ -function CFI(ρ, dρ, M) - m_num = length(M) - p = zero(ComplexF64) - dp = zero(ComplexF64) - F = 0. - for i in 1:m_num - mp = M[i] - p += tr(ρ * mp) - dp = tr(dρ * mp) - cadd = 0. - if p != 0 - cadd = (dp^2) / p - end - F += cadd - end - real(F) -end -function CFI(M::Vector{Matrix{T}}, H::Vector{Matrix{T}}, ∂H_∂x::Matrix{T}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, times) where {T <: Complex,R <: Real} - dim = size(H[1])[1] - Δt = times[2] - times[1] - ρt = evolute(H[1], Liouville_operator, γ, times, times[1]) * (ρ_initial |> vec) - ∂ρt_∂x = -im * Δt * liouville_commu(∂H_∂x) * ρt - for t in 2:length(times) - expL = evolute(H[t], Liouville_operator, γ, times, times[t]) - ρt= expL * ρt - ∂ρt_∂x= -im * Δt * liouville_commu(∂H_∂x) * ρt + expL * ∂ρt_∂x - end - CFI(ρt|> vec2mat, ∂ρt_∂x|> vec2mat, M) -end - -function CFIM(ρ, dρ, M) - m_num = length(M) - cfim = [tr.(kron(dρ', dρ).*M[i]) / tr(ρ*M[i]) for i in 1:m_num] |> sum -end -function CFIM(M::Vector{Matrix{T}}, H::Vector{Matrix{T}}, ∂H_∂x::Matrix{T}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, times) where {T <: Complex,R <: Real} - dim = size(H[1])[1] - Δt = times[2] - times[1] - para_num = length(∂H_∂x) - ρt = evolute(H[1], Liouville_operator, γ, times, times[1]) * ρ_initial[:] - ∂ρt_∂x = [-im * Δt * liouville_commu(∂H_∂x[i]) * ρt for i in 1:para_num] - for t in 2:length(times) - expL = evolute(H[t], Liouville_operator, γ, times, times[t]) - ρt= expL * ρt - ∂ρt_∂x= [-im * Δt * liouville_commu(∂H_∂x[i]) * ρt for i in 1:para_num] + [expL] .* ∂ρt_∂x - end - CFIM(ρt|> vec2mat, ∂ρt_∂x|> vec2mat, M) -end - - -function SLD(ρ::Matrix{T}, ∂ρ_∂x::Matrix{T}) where {T <: Complex} - 2 * pinv(kron(ρ |> transpose, ρ |> one) + kron(ρ |> one, ρ)) * vec(∂ρ_∂x) |> vec2mat -end -function SLD(ρ::Vector{T},∂ρ_∂x::Vector{T}) where {T <: Complex} - SLD(ρ |> vec2mat, ∂ρ_∂x |> vec2mat) -end -function SLD(ρ::Matrix{T}, ∂ρ_∂x::Vector{Matrix{T}}) where {T <: Complex} - (x->SLD(ρ, x)).(∂ρ_∂x) -end -function SLD_qr(ρ::Matrix{T}, ∂ρ_∂x::Matrix{T}) where {T <: Complex} - 2 * (qr(kron(ρ |> transpose, ρ |> one) + kron(ρ |> one, ρ), Val(true)) \ vec(∂ρ_∂x)) |> vec2mat -end -function SLD_eig(ρ::Array{T}, dρ::Array{T})::Array{T} where {T <: Complex} - dim = size(ρ)[1] - if typeof(dρ) == Array{T,2} - purity = tr(ρ * ρ) - SLD_res = zeros(T, dim, dim) - if abs(1 - purity) < 1e-8 - SLD_res = 2 * dρ - else - val, vec_mat = eigen(ρ) - for fi in 1:dim - for fj in 1:dim - coeff = 2 / (val[fi] + val[fj]) - SLD_res[fi, fj] = coeff * (vec_mat[:,fi]' * (dρ * vec_mat[:, fj])) - end - end - SLD_res[findall(SLD_res == Inf)] .= 0. - SLD_res = vec_mat * (SLD_res * vec_mat') - end - else - # multi-parameter scenario - purity = tr(ρ * ρ) - if abs(1 - purity) < 1e-8 - SLD_res = [2 * dρ[i] for i in 1:length(dρ)] - else - # SLD_res = [zeros(T,dim,dim) for i in 1:length(dρ)] - dim = ndims(ρ) - val, vec_mat = eigens(ρ) - for para_i in 1:length(dρ) - SLD_tp = zeros(T, dim, dim) - for fi in 1:dim - for fj in 1:dim - coeff = 2. / (val[fi] + val[fj]) - SLD_tp[fi][fj] = coeff * (vec[fi]' * (dρ[para_i] * vec[fj])) - end - end - SLD_tp[findall(SLD_rp == Inf)] .= 0. - SLD_res[para_i] = vec_mat * (SLD_tp * vec_mat') - end - end - end - SLD_res -end -function RLD(ρ::Matrix{T}, dρ::Matrix{T}) where {T <: Complex} - dρ * pinv(ρ) -end -function QFI_RLD(ρ, dρ) - RLD_tp = RLD(ρ, dρ) - F = tr(ρ * RLD_tp * RLD_tp') - F |> real -end -function QFI(ρ, dρ) - SLD_tp = SLD(ρ, dρ) - SLD2_tp = SLD_tp * SLD_tp - F = tr(ρ * SLD2_tp) - F |> real -end -function QFIM(ρ, dρ) - SLD_tp = SLD(ρ, dρ) - [0.5*ρ] .* (kron(SLD_tp, SLD_tp') + kron(SLD_tp', SLD_tp)).|> tr .|>real -end -function QFI(H::Vector{Matrix{T}}, ∂H_∂x::Matrix{T}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, times) where {T <: Complex,R <: Real} - dim = size(H[1])[1] - Δt = times[2] - times[1] - ρt = evolute(H[1], Liouville_operator, γ, times, times[1]) * (ρ_initial |> vec) - ∂ρt_∂x = -im * Δt * liouville_commu(∂H_∂x) * ρt - for t in 2:length(times) - expL = evolute(H[t], Liouville_operator, γ, times, times[t]) - ρt= expL * ρt - ∂ρt_∂x= -im * Δt * liouville_commu(∂H_∂x) * ρt + expL * ∂ρt_∂x - end - QFI(ρt|> vec2mat, ∂ρt_∂x|> vec2mat) -end -function QFIM(H::Vector{Matrix{T}}, ∂H_∂x::Vector{Matrix{T}}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, times) where {T <: Complex,R <: Real} - dim = size(H[1])[1] - Δt = times[2] - times[1] - para_num = length(∂H_∂x) - ρt = evolute(H[1], Liouville_operator, γ, times, times[1]) * ρ_initial[:] - ∂ρt_∂x = [-im * Δt * liouville_commu(∂H_∂x[i]) * ρt for i in 1:para_num] - for t in 2:length(times) - expL = evolute(H[t], Liouville_operator, γ, times, times[t]) - ρt= expL * ρt - ∂ρt_∂x= [-im * Δt * liouville_commu(∂H_∂x[i]) * ρt for i in 1:para_num] + [expL] .* ∂ρt_∂x - end - QFIM(ρt|> vec2mat, ∂ρt_∂x|> vec2mat) -end -function CFI(M, system) - CFI(M,Htot(system.freeHamiltonian, system.control_Hamiltonian, system.control_coefficients), system.Hamiltonian_derivative[1], system.ρ_initial, system.Liouville_operator, system.γ, system.times) -end -function CFIM(M, system) - CFIM(M,Htot(system.freeHamiltonian, system.control_Hamiltonian, system.control_coefficients), system.Hamiltonian_derivative, system.ρ_initial, system.Liouville_operator, system.γ, system.times) -end -function QFI(system) - QFI(Htot(system.freeHamiltonian, system.control_Hamiltonian, system.control_coefficients), system.Hamiltonian_derivative[1], system.ρ_initial, system.Liouville_operator, system.γ, system.times) -end -function QFIM(system) - QFIM(Htot(system.freeHamiltonian, system.control_Hamiltonian, system.control_coefficients), system.Hamiltonian_derivative, system.ρ_initial, system.Liouville_operator, system.γ, system.times) -end \ No newline at end of file diff --git a/src/Common/common.jl b/src/Common/common.jl deleted file mode 100644 index c2391cb..0000000 --- a/src/Common/common.jl +++ /dev/null @@ -1,104 +0,0 @@ -function RunMixed( grape, particle_num= 10, c0=0.5, c1=0.5, c2=0.5, prec_rough = 0.1, sd=1234, episode=400) - println("Combined strategies") - println("searching initial controls with particle swarm optimization ") - tnum = length(grape.times) - ctrl_num = length(grape.control_Hamiltonian) - particles, velocity = repeat(grape, particle_num), [[10*ones(tnum) for i in 1:ctrl_num] for j in 1:particle_num] - pbest = [[zeros(tnum) for i in 1:ctrl_num] for j in 1:particle_num] - gbest = [zeros(tnum) for i in 1:ctrl_num] - velocity_best = [zeros(Float64, tnum) for i in 1:ctrl_num] - p_fit = zeros(particle_num) - qfi_ini = QFI(grape) - println("initial QFI is $(qfi_ini)") - for ei in 1:episode - fit_pre = 0.0 - fit = 0.0 - for pi in 1:particle_num - propagate!(particles[pi]) - f_now = QFI(particles[pi]) - - if f_now > p_fit[pi] - p_fit[pi] = f_now - for di in 1:ctrl_num - for ni in 1:tnum - pbest[pi][di][ni] = particles[pi].control_coefficients[di][ni] - end - end - end - end - for pj in 1:particle_num - if p_fit[pj] > fit - fit = p_fit[pj] - for dj in 1:ctrl_num - for nj in 1:tnum - gbest[dj][nj] = particles[pj].control_coefficients[dj][nj] - velocity_best[dj][nj] = velocity[pj][dj][nj] - end - end - end - end - Random.seed!(sd) - for pk in 1:particle_num - for dk in 1:ctrl_num - for ck in 1:tnum - velocity[pk][dk][ck] = c0*velocity[pk][dk][ck] + c1*rand()*(pbest[pk][dk][ck] - particles[pk].control_coefficients[dk][ck]) - + c2*rand()*(gbest[dk][ck] - particles[pk].control_coefficients[dk][ck]) - particles[pk].control_coefficients[dk][ck] = particles[pk].control_coefficients[dk][ck] + velocity[pk][dk][ck] - end - end - end - fit_pre = fit - if abs(fit-fit_pre) < prec_rough - grape.control_coefficients = gbest - println("PSO strategy finished, switching to GRAPE") - Run(grape) - return nothing - end - print("current QFI is $fit ($ei epochs) \r") - end -end - -function RunODE(grape) - println("quantum parameter estimation") - if length(grape.Hamiltonian_derivative) == 1 - println("single parameter estimation scenario") - qfi_ini = QFI(grape) - qfi_list = [qfi_ini] - println("initial QFI is $(qfi_ini)") - gradient_QFIM_ODE!(grape) - while true - qfi_now = QFI(grape) - gradient_QFIM_ODE!(grape) - if 0 < (qfi_now - qfi_ini) < 1e-4 - println("\n Iteration over, data saved.") - println("Final QFI is ", qfi_now) - save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "qfi", qfi_list) - break - else - qfi_ini = qfi_now - append!(qfi_list,qfi_now) - print("current QFI is ", qfi_now, " ($(qfi_list|>length) epochs) \r") - end - end - else - println("multiple parameters estimation scenario") - f_ini =1/(grape |> QFIM |> inv |> tr) - f_list = [f_ini] - println("initial 1/tr(F^-1) is $(f_ini)") - gradient_QFIM!(grape) - while true - f_now = 1/(grape |> QFIM |> inv |> tr) - gradient_QFIM!(grape) - if 0< f_now - f_ini < 1e-4 - println("\n Iteration over, data saved.") - println("Final 1/tr(F^-1) is ", f_now) - save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "f", f_list) - break - else - f_ini = f_now - append!(f_list,f_now) - print("current 1/tr(F^-1) is ", f_now, " ($(f_list|>length) epochs) \r") - end - end - end -end diff --git a/src/Common/utils.jl b/src/Common/utils.jl deleted file mode 100644 index 4b7630c..0000000 --- a/src/Common/utils.jl +++ /dev/null @@ -1,77 +0,0 @@ -using SparseArrays -""" -Pauli matrces -""" -sigmax() = [.0im 1.;1. 0.] -sigmay() = [0. -1.0im;1.0im 0.] -sigmaz() = [1.0 .0im;0. -1.] -sigmap() = [.0im 1.;0. 0.] -sigmam() = [.0im 0.;1. 0.] -sigmax(i, N) = kron(I(2^(i - 1)), sigmax(), I(2^(N - i))) -sigmay(i, N) = kron(I(2^(i - 1)), sigmay(), I(2^(N - i))) -sigmaz(i, N) = kron(I(2^(i - 1)), sigmaz(), I(2^(N - i))) -sigmap(i, N) = kron(I(2^(i - 1)), sigmap(), I(2^(N - i))) -sigmam(i, N) = kron(I(2^(i - 1)), sigmam(), I(2^(N - i))) - -""" - destroy() - -""" -function destroy(M) - spdiagm(M, M, 1 => map(x -> x |> sqrt, 1:(M - 1))) -end - -""" - vec2mat() -transform a vector to its matrix form -""" -function vec2mat(x::Vector{T}) where {T <: Number} - reshape(x, x |> length |> sqrt |> Int, :) -end -function vec2mat(x) - vec2mat.(x) -end -function vec2mat(x::Matrix) - throw(ErrorException("vec2mating a matrix of size $(size(x))")) -end - -""" - repeat() - -""" -function Base.repeat(system, N) - [system for i in 1:N] -end -function Base.repeat(system, M, N) - reshape(repeat(system, M * N), M, N) -end - -""" - filterZeros!() - -""" -function filterZeros!(x::Matrix{T}) where {T <: Complex} - x[abs.(x) .< eps()] .= zero(T) - x -end -function filterZeros!(x) - filterZeros!.(x) -end - -""" - t2Num() - -""" -function t2Num(t0, dt, t) - Int(round((t - t0) / dt)) + 1 -end - -""" - basis() - -""" -function basis(dim, si, ::T)::Array{T} where {T <: Complex} - result = zeros(T, dim) - result[si] = 1.0 - result -end \ No newline at end of file diff --git a/src/Control/GRAPE.jl b/src/Control/GRAPE.jl deleted file mode 100644 index 6f5ef50..0000000 --- a/src/Control/GRAPE.jl +++ /dev/null @@ -1,280 +0,0 @@ -abstract type ControlSystem end -mutable struct GrapeControl{T <: Complex,M <: Real} <: ControlSystem - freeHamiltonian::Matrix{T} - Hamiltonian_derivative::Vector{Matrix{T}} - ρ_initial::Matrix{T} - times::StepRangeLen{M,Base.TwicePrecision{M},Base.TwicePrecision{M}} - Liouville_operator::Vector{Matrix{T}} - γ::Vector{M} - control_Hamiltonian::Vector{Matrix{T}} - control_coefficients::Vector{Vector{M}} - ϵ::M - ρ::Vector{Matrix{T}} - ∂ρ_∂x::Vector{Vector{Matrix{T}}} - GrapeControl(freeHamiltonian::Matrix{T}, Hamiltonian_derivative::Vector{Matrix{T}}, ρ_initial::Matrix{T}, - times::StepRangeLen{M,Base.TwicePrecision{M},Base.TwicePrecision{M}}, - Liouville_operator::Vector{Matrix{T}},γ::Vector{M}, control_Hamiltonian::Vector{Matrix{T}}, - control_coefficients::Vector{Vector{M}}, ϵ=0.01, ρ=Vector{Matrix{T}}(undef, 1), - ∂ρ_∂x=Vector{Vector{Matrix{T}}}(undef, 1),∂ρ_∂V=Vector{Vector{Matrix{T}}}(undef, 1)) where {T <: Complex,M <: Real} = - new{T,M}(freeHamiltonian, Hamiltonian_derivative, ρ_initial, times, Liouville_operator, γ, control_Hamiltonian, - control_coefficients, ϵ, ρ, ∂ρ_∂x) -end - -function gradient_CFI!(grape::GrapeControl{T}, M) where {T <: Complex} - δI = gradient(x->CFI(M, Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients)[1].|>real - grape.control_coefficients += grape.ϵ*δI -end -function gradient_CFIM!(grape::GrapeControl{T}, M) where {T <: Complex} - δI = gradient(x->CFIM(M, Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients).|>real - grape.control_coefficients += grape.ϵ*δI -end -function gradient_CFI_ADAM!(grape::GrapeControl{T}, M) where {T <: Complex} - δI = gradient(x->CFI(M, Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients)[1].|>real - Adam!(grape, δI) -end -function gradient_CFIM_ADAM!(grape::GrapeControl{T}, M, mt, vt) where {T <: Complex} - δI = gradient(x->CFIM(M, Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients).|>real - Adam!(grape, δI) -end - -function gradient_QFI!(grape::GrapeControl{T}) where {T <: Complex} - δF = gradient(x->QFI(Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients)[1].|>real - grape.control_coefficients += grape.ϵ*δF -end -function gradient_QFIM!(grape::GrapeControl{T}) where {T <: Complex} - δF = gradient(x->1/(QFIM(Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative, grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times) |> pinv |> tr |>real), grape.control_coefficients).|>real |>sum - grape.control_coefficients += grape.ϵ*δF -end -function gradient_QFI_ADAM!(grape::GrapeControl{T}) where {T <: Complex} - δF = gradient(x->QFI(Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative[1], grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times), grape.control_coefficients)[1].|>real - Adam!(grape, δF) -end -function gradient_QFIM_ADAM!(grape::GrapeControl{T}) where {T <: Complex} - δF = gradient(x->1/(QFIM(Htot(grape.freeHamiltonian, grape.control_Hamiltonian, x), grape.Hamiltonian_derivative, grape.ρ_initial, grape.Liouville_operator, grape.γ, grape.times) |> pinv |> tr |>real), grape.control_coefficients).|>real |>sum - Adam!(grape, δF) -end -function gradient_QFI_analitical_ADAM!(grape::GrapeControl{T}) where {T <: Complex} - δF = propagate_analitical!(grape) - Adam!(grape, δF) -end -function gradient_QFIM_ODE(grape::GrapeControl) - H = Htot(grape.freeHamiltonian, grape.control_Hamiltonian, grape.control_coefficients) - Δt = grape.times[2] - grape.times[1] - t_num = length(grape.times) - para_num = length(grape.Hamiltonian_derivative) - ctrl_num = length(grape.control_Hamiltonian) - tspan(j) = (grape.times[1], grape.times[j]) - tspan() = (grape.times[1], grape.times[end]) - u0 = grape.ρ_initial |> vec - evo(p, t) = evolute(p[t2Num(tspan()[1], Δt, t)], grape.Liouville_operator, grape.γ, grape.times, t2Num(tspan()[1], Δt, t)) - f(u, p, t) = evo(p, t) * u - prob = DiscreteProblem(f, u0, tspan(), H,dt=Δt) - ρt = solve(prob).u - ∂ρt_∂x = Vector{Vector{Vector{eltype(u0)}}}(undef, 1) - for para in 1:para_num - devo(p, t) = -1.0im * Δt * liouville_commu(grape.Hamiltonian_derivative[para]) * evo(p, t) - du0 = devo(H, tspan()[1]) * u0 - g(du, p, t) = evo(p, t) * du + devo(p, t) * ρt[t2Num(tspan()[1], Δt, t)] - dprob = DiscreteProblem(g, du0, tspan(), H,dt=Δt) - ∂ρt_∂x[para] = solve(dprob).u - end - δρt_δV = Matrix{Vector{Vector{eltype(u0)}}}(undef,ctrl_num,length(grape.times)) - for ctrl in 1:ctrl_num - for j in 1:t_num - devo(p, t) = -1.0im * Δt * liouville_commu(grape.control_Hamiltonian[ctrl]) * evo(p, t) - du0 = devo(H, tspan()[1]) * u0 - g(du, p, t) = evo(p, t) * du + devo(p, t) * ρt[t2Num(tspan()[1], Δt, t)] - dprob = DiscreteProblem(g, du0, tspan(j), H,dt=Δt) - δρt_δV[ctrl,j] = solve(dprob).u - end - end - ∂xδρt_δV = Array{Vector{eltype(u0)}, 3}(undef,para_num, ctrl_num,length(grape.times)) - for para in 1:para_num - for ctrl in 1:ctrl_num - dxevo = -1.0im * Δt * liouville_commu(grape.Hamiltonian_derivative[para]) - dkevo = -1.0im * Δt * liouville_commu(grape.control_Hamiltonian[ctrl]) - for j in 1:t_num - g(du, p, t) = dxevo * dkevo * evo(p, t) * ρt[t2Num(tspan()[1], Δt, t)] + - dxevo * evo(p, t) * δρt_δV[ctrl, j][t2Num(tspan()[1], Δt, t)] + - dkevo * evo(p, t) * ∂ρt_∂x[para][t2Num(tspan()[1], Δt, t)] + - evo(p, t) * du - du0 = dxevo * dkevo * evo(H,tspan()[1]) * ρt[t2Num(tspan()[1], Δt, tspan()[1])] - dprob = DiscreteProblem(g, du0, tspan(j), H, dt=Δt) - ∂xδρt_δV[para, ctrl, j] = solve(dprob).u[end] - end - end - end - δF = grape.control_coefficients .|> zero - for para in 1:para_num - SLD_tp = SLD(ρt[end], ∂ρt_∂x[para][end]) - for ctrl in 1:ctrl_num - for j in 1:t_num - δF[ctrl][j] -= 2 * tr((∂xδρt_δV[para,ctrl,j]|> vec2mat) * SLD_tp) - - tr((δρt_δV[ctrl, j][end] |> vec2mat) * SLD_tp^2) |> real - end - end - end - δF -end -function gradient_QFIM_ODE!(grape::GrapeControl{T}) where {T <: Complex} - grape.control_coefficients += grape.ϵ * gradient_QFIM_ODE(grape) -end - -function Run(M, grape::GrapeControl{T}) where {T<: Complex} - println("classical parameter estimation") - if length(grape.Hamiltonian_derivative) == 1 - println("single parameter estimation scenario") - cfi_ini = CFI(M, grape) - cfi_list = [cfi_ini] - println("initial CFI is $(cfi_ini)") - gradient_CFI!(M, grape) - while true - cfi_now = CFI(M, grape) - gradient_CFI!(M, grape) - if 0 < (cfi_now - cfi_ini) < 1e-4 - println("\n Iteration over, data saved.") - println("Final CFI is ", cfi_now) - save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "cfi", cfi_list) - break - else - cfi_ini = cfi_now - append!(cfi_list,cfi_now) - print("current CFI is ", cfi_now, " ($(cfi_list|>length) epochs) \r") - end - end - else - println("multiple parameters estimation scenario") - f_ini =1/(grape |> CFIM |> inv |> tr) - f_list = [f_ini] - println("initial 1/tr(F^-1) is $(f_ini)") - gradient_CFIM!(M, grape) - while true - f_now = 1/(grape |> CFIM |> inv |> tr) - gradient_CFIM!(M, grape) - if 0< f_now - f_ini < 1e-4 - println("\n Iteration over, data saved.") - println("Final 1/tr(I^-1) is ", f_now) - save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "f", f_list) - break - else - f_ini = f_now - append!(f_list,f_now) - print("current 1/tr(I^-1) is ", f_now, " ($(f_list|>length) epochs) \r") - end - end - end -end -function RunADAM(grape) - println("AutoGrape strategies") - println("quantum parameter estimation") - if length(grape.Hamiltonian_derivative) == 1 - println("single parameter estimation scenario") - qfi_ini = QFI(grape) - qfi_list = [qfi_ini] - println("initial QFI is $(qfi_ini)") - gradient_QFI_ADAM!(grape) - while true - qfi_now = QFI(grape) - gradient_QFI_ADAM!(grape) - if 0 < (qfi_now - qfi_ini) < 1e-4 - println("\n Iteration over, data saved.") - println("Final QFI is ", qfi_now) - save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "qfi", qfi_list) - break - else - qfi_ini = qfi_now - append!(qfi_list,qfi_now) - print("current QFI is ", qfi_now, " ($(qfi_list|>length) epochs) \r") - end - end - else - println("multiple parameters estimation scenario") - f_ini =1/(grape |> QFIM |> inv |> tr) - f_list = [f_ini] - println("initial 1/tr(F^-1) is $(f_ini)") - gradient_QFIM_ADAM!(grape) - while true - f_now = 1/(grape |> QFIM |> inv |> tr) - gradient_QFIM_ADAM!(grape) - if 0< f_now - f_ini < 1e-4 - println("\n Iteration over, data saved.") - println("Final 1/tr(F^-1) is ", f_now) - save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "f", f_list) - break - else - f_ini = f_now - append!(f_list,f_now) - print("current 1/tr(F^-1) is ", f_now, " ($(f_list|>length) epochs) \r") - end - end - end -end -function RunAnaliticalADAM(grape) - println("Analitical strategies") - println("quantum parameter estimation") - if length(grape.Hamiltonian_derivative) == 1 - println("single parameter estimation scenario") - qfi_ini = QFI(grape) - qfi_list = [qfi_ini] - println("initial QFI is $(qfi_ini)") - gradient_QFI_analitical_ADAM!(grape) - while true - qfi_now = QFI(grape) - gradient_QFI_analitical_ADAM!(grape) - if 0 < (qfi_now - qfi_ini) < 1e-4 - println("\n Iteration over, data saved.") - println("Final QFI is ", qfi_now) - save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "qfi", qfi_list) - break - else - qfi_ini = qfi_now - append!(qfi_list,qfi_now) - print("current QFI is ", qfi_now, " ($(qfi_list|>length) epochs) \r") - end - end - end -end -function RunADAM(M, grape::GrapeControl{T}) where {T<: Complex} - println("classical parameter estimation") - if length(grape.Hamiltonian_derivative) == 1 - println("single parameter estimation scenario") - cfi_ini = CFI(M, grape) - cfi_list = [cfi_ini] - println("initial CFI is $(cfi_ini)") - gradient_CFI_ADAM!(M, grape) - while true - cfi_now = CFI(M, grape) - gradient_CFI_ADAM!(M, grape) - if 0 < (cfi_now - cfi_ini) < 1e-4 - println("\n Iteration over, data saved.") - println("Final CFI is ", cfi_now) - save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "cfi", cfi_list) - break - else - cfi_ini = cfi_now - append!(cfi_list,cfi_now) - print("current CFI is ", cfi_now, " ($(cfi_list|>length) epochs) \r") - end - end - else - println("multiple parameters estimation scenario") - f_ini =1/(grape |> CFIM |> inv |> tr) - f_list = [f_ini] - println("initial 1/tr(F^-1) is $(f_ini)") - gradient_CFIM!(M, grape) - while true - f_now = 1/(grape |> CFIM |> inv |> tr) - gradient_CFIM!(M, grape) - if 0< f_now - f_ini < 1e-4 - println("\n Iteration over, data saved.") - println("Final 1/tr(I^-1) is ", f_now) - save("controls.jld", "controls", grape.control_coefficients, "time_span", grape.times, "f", f_list) - break - else - f_ini = f_now - append!(f_list,f_now) - print("current 1/tr(I^-1) is ", f_now, " ($(f_list|>length) epochs) \r") - end - end - end -end \ No newline at end of file diff --git a/src/Control/PSO.jl b/src/Control/PSO.jl deleted file mode 100644 index 3da14ff..0000000 --- a/src/Control/PSO.jl +++ /dev/null @@ -1,118 +0,0 @@ -function RunPSO(grape::GrapeControl{T}; particle_num= 10, c0=0.5, c1=0.5, c2=0.5,v0=0.1, sd=1234, episode=400) where {T<: Complex} - println("PSO strategies") - println("searching optimal controls with particle swarm optimization ") - tnum = length(grape.times) - ctrl_num = length(grape.control_Hamiltonian) - particles = repeat(grape, particle_num) - velocity = v0.*rand(ctrl_num, tnum, particle_num)|>SharedArray - pbest = zeros(ctrl_num, tnum, particle_num) - gbest = zeros(ctrl_num, tnum) |> SharedArray - velocity_best = zeros(ctrl_num,tnum) - p_fit = zeros(particle_num) - qfi_ini = QFI(grape) - println("initial QFI is $(qfi_ini)") - fit_pre = 0.0 - fit = 0.0 - for ei in 1:episode - @inbounds @threads for pi in 1:particle_num - propagate!(particles[pi]) - f_now = QFI(particles[pi]) - - if f_now > p_fit[pi] - p_fit[pi] = f_now - for di in 1:ctrl_num - for ni in 1:tnum - pbest[di,ni,pi] = particles[pi].control_coefficients[di][ni] - end - end - end - end - @inbounds @threads for pj in 1:particle_num - if p_fit[pj] > fit - fit = p_fit[pj] - for dj in 1:ctrl_num - @inbounds @threads for nj in 1:tnum - gbest[dj, nj] = particles[pj].control_coefficients[dj][nj] - velocity_best[dj, nj] = velocity[dj, nj, pj] - end - end - end - end - Random.seed!(sd) - @inbounds @threads for pk in 1:particle_num - for dk in 1:ctrl_num - @inbounds @threads for ck in 1:tnum - velocity[dk, ck, pk] = c0*velocity[dk, ck, pk] + c1*rand()*(pbest[dk, ck, pk] - particles[pk].control_coefficients[dk][ck]) - + c2*rand()*(gbest[dk, ck] - particles[pk].control_coefficients[dk][ck]) - particles[pk].control_coefficients[dk][ck] = particles[pk].control_coefficients[dk][ck] + velocity[dk, ck, pk] - end - end - end - # println(particles[1].control_coefficients[1][1]) - if ei == episode#abs(fit-fit_pre) < 1e-5 - grape.control_coefficients = gbest - println("Final QFI is ", fit) - return nothing - end - fit_pre = fit - print("current QFI is $fit ($ei epochs) \r") - end -end -function RunPSOAdam(grape::GrapeControl{T}; particle_num= 10, c0=0.5, c1=0.5, c2=0.5,v0=0.01, sd=1234, episode=400) where {T<: Complex} - println("PSO strategies") - println("searching optimal controls with particle swarm optimization ") - tnum = length(grape.times) - ctrl_num = length(grape.control_Hamiltonian) - particles, velocity = repeat(grape, particle_num), [[v0*ones(tnum) for i in 1:ctrl_num] for j in 1:particle_num] - pbest = [[zeros(tnum) for i in 1:ctrl_num] for j in 1:particle_num] - gbest = [zeros(tnum) for i in 1:ctrl_num] - velocity_best = [zeros(Float64, tnum) for i in 1:ctrl_num] - p_fit = zeros(particle_num) - qfi_ini = QFI(grape) - println("initial QFI is $(qfi_ini)") - fit_pre = 0.0 - fit = 0.0 - for ei in 1:episode - for pi in 1:particle_num - propagate!(particles[pi]) - f_now = QFI(particles[pi]) - - if f_now > p_fit[pi] - p_fit[pi] = f_now - for di in 1:ctrl_num - for ni in 1:tnum - pbest[pi][di][ni] = particles[pi].control_coefficients[di][ni] - end - end - end - end - for pj in 1:particle_num - if p_fit[pj] > fit - fit = p_fit[pj] - for dj in 1:ctrl_num - for nj in 1:tnum - gbest[dj][nj] = particles[pj].control_coefficients[dj][nj] - velocity_best[dj][nj] = velocity[pj][dj][nj] - end - end - end - end - Random.seed!(sd) - for pk in 1:particle_num - for dk in 1:ctrl_num - for ck in 1:tnum - velocity[pk][dk][ck] = c0*velocity[pk][dk][ck] + c1*(rand()-0.5)*(pbest[pk][dk][ck] - particles[pk].control_coefficients[dk][ck]) - + c2*(rand()-0.5)*(gbest[dk][ck] - particles[pk].control_coefficients[dk][ck]) - end - end - Adam!(particles[pk], velocity[pk]) - end - if abs(fit-fit_pre) < 1e-5 - grape.control_coefficients = gbest - println("Final QFI is ", fit) - return nothing - end - fit_pre = fit - print("current QFI is $fit ($ei epochs) \r") - end -end \ No newline at end of file diff --git a/src/Dynamics/Adam.jl b/src/Dynamics/Adam.jl deleted file mode 100644 index a396344..0000000 --- a/src/Dynamics/Adam.jl +++ /dev/null @@ -1,20 +0,0 @@ -""" - Adam() - -""" -function Adam(gt, t, para, m_t, v_t, alpha=0.01, beta1=0.90, beta2=0.99, epsilon=1e-8) - t = t+1 - m_t = beta1*m_t + (1-beta1)*gt - v_t = beta2*v_t + (1-beta2)*(gt*gt) - m_cap = m_t/(1-(beta1^t)) - v_cap = v_t/(1-(beta2^t)) - para = para+(alpha*m_cap)/(sqrt(v_cap)+epsilon) - return para, m_t, v_t -end -function Adam!(grape, δ, mt=0.0, vt=0.0) - for ctrl in 1:length(δ) - for ti in 1:length(grape.times) - grape.control_coefficients[ctrl][ti], mt, vt = Adam(δ[ctrl][ti], ti, grape.control_coefficients[ctrl][ti], mt, vt) - end - end -end \ No newline at end of file diff --git a/src/Dynamics/dynamcs.jl b/src/Dynamics/dynamcs.jl deleted file mode 100644 index 15dcd6b..0000000 --- a/src/Dynamics/dynamcs.jl +++ /dev/null @@ -1,279 +0,0 @@ - -""" - liouville_commu() - -""" -function liouville_commu(H) - kron(one(H), H) - kron(H |> transpose, one(H)) -end - -""" - liouville_commu_dissip() - -""" -function liouville_dissip(Γ) - kron(Γ |> conj, Γ) - 0.5 * kron((Γ |> transpose) * Γ |> conj, Γ |> one) - 0.5 * kron(Γ |> one, Γ' * Γ) -end - -""" - dissipation() - -""" -function dissipation(Γ::Vector{Matrix{T}}, γ::Vector{R}, t::Real) where {T <: Complex,R <: Real} - [γ[i] * liouville_dissip(Γ[i]) for i in 1:length(Γ)] |> sum -end -function dissipation(Γ::Vector{Matrix{T}}, γ::Vector{Vector{R}}, t::Real) where {T <: Complex,R <: Real} - [γ[i][t] * liouville_dissip(Γ[i]) for i in 1:length(Γ)] |> sum -end - -""" - free_evolution() - -""" -function free_evolution(H0) - -1.0im * liouville_commu(H0) -end - -""" - liouvillian() - -""" -function liouvillian(H::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, γ, t::Real) where {T <: Complex} - freepart = liouville_commu(H) - dissp = norm(γ) +1 ≈ 1 ? freepart|>zero : dissipation(Liouville_operator, γ, t) - -1.0im * freepart + dissp -end - -""" - Htot() -""" -function Htot(H0::Matrix{T}, control_Hamiltonian::Vector{Matrix{T}}, control_coefficients) where {T <: Complex} - Htot = [H0] .+ ([control_coefficients[i] .* [control_Hamiltonian[i]] for i in 1:length(control_coefficients)] |> sum ) -end - -""" - evolute() - -""" -function evolute(H, Liouville_operator, γ, times, t) - tj = Int(round((t - times[1]) / (times[2] - times[1]))) + 1 - dt = times[2] - times[1] - Ld = dt * liouvillian(H, Liouville_operator, γ, tj) - exp(Ld) -end - -""" - evolute_ODE!() - -""" -function evolute_ODE!(grape::GrapeControl) - H(p) = Htot(grape.freeHamiltonian, grape.control_Hamiltonian, p) - dt = grape.times[2] - grape.times[1] - tspan = (grape.times[1], grape.times[end]) - u0 = grape.ρ_initial - Γ = grape.Liouville_operator - f(u, p, t) = -im * (H(p)[t2Num(tspan[1], dt, t)] * u + u * H(p)[t2Num(tspan[1], dt, t)]) + - ([grape.γ[i] * (Γ[i] * u * Γ[i]' - (Γ[i]' * Γ[i] * u + u * Γ[i]' * Γ[i] )) for i in 1:length(Γ)] |> sum) - prob = ODEProblem(f, u0, tspan, grape.control_coefficients, saveat=dt) - sol = solve(prob) - sol.u -end - -""" - propagate() - -""" -function propagate(H0::Matrix{T}, ∂H_∂x::Vector{Matrix{T}}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, - γ, control_Hamiltonian::Vector{Matrix{T}}, control_coefficients::Vector{Vector{R}}, times) where {T <: Complex,R <: Real} - dim = size(H0)[1] - para_num = length(∂H_∂x) - H = Htot(H0, control_Hamiltonian, control_coefficients) - ρt = [Vector{ComplexF64}(undef, dim^2) for i in 1:length(times)] - ∂ρt_∂x = [[Vector{ComplexF64}(undef, dim^2) for i in 1:length(times)] for para in 1:para_num] - Δt = times[2] - times[1] - ρt[1] = evolute(H[1], Liouville_operator, γ, times, times[1]) * (ρ_initial |> vec) - for para in 1:para_num - ∂ρt_∂x[para][1] = -im * Δt * liouville_commu(∂H_∂x[para]) * ρt[1] - end - for t in 2:length(times) - expL = evolute(H[t], Liouville_operator, γ, times, times[t]) - ρt[t] = expL * ρt[t - 1] - for para in para_num - ∂ρt_∂x[para][t] = -im * Δt * liouville_commu(∂H_∂x[para]) * ρt[t] + expL * ∂ρt_∂x[para][t - 1] - end - end - ρt .|> vec2mat, ∂ρt_∂x .|> vec2mat -end - -""" - propagate!() - -""" -function propagate!(grape::GrapeControl) - grape.ρ, grape.∂ρ_∂x = propagate(grape.freeHamiltonian, grape.Hamiltonian_derivative, grape.ρ_initial, - grape.Liouville_operator, grape.γ, grape.control_Hamiltonian, - grape.control_coefficients, grape.times ) -end - -""" - propagate_analytical() -""" -function propagate_analytical(H0::Matrix{T}, ∂H_∂x::Vector{Matrix{T}}, ρ_initial::Matrix{T}, Liouville_operator::Vector{Matrix{T}}, - γ, control_Hamiltonian::Vector{Matrix{T}}, control_coefficients::Vector{Vector{R}}, times) where {T <: Complex,R <: Real} - - dim = size(H0)[1] - tnum = length(times) - para_num = length(∂H_∂x) - ctrl_num = length(control_Hamiltonian) - Δt = times[2] - times[1] - - H = Htot(H0, control_Hamiltonian, control_coefficients) - - ρt = [Vector{ComplexF64}(undef, dim^2) for i in 1:tnum] - ∂ρt_∂x = [[Vector{ComplexF64}(undef, dim^2) for para in 1:para_num] for i in 1:tnum] - δρt_δV = [[] for ctrl in 1:ctrl_num] - ∂xδρt_δV = [[[] for ctrl in 1:ctrl_num] for i in 1:para_num] - ∂H_L = [Matrix{ComplexF64}(undef, dim^2,dim^2) for i in 1:para_num] - Hc_L = [Matrix{ComplexF64}(undef, dim^2,dim^2) for i in 1:ctrl_num] - - ρt[1] = evolute(H[1], Liouville_operator, γ, times, times[1]) * (ρ_initial |> vec) - for pi in 1:para_num - ∂ρt_∂x[1][pi] = -im * Δt * ∂H_L[pi] * (ρ_initial |> vec) - ∂H_L[pi] = liouville_commu(∂H_∂x[pi]) - for ci in 1:ctrl_num - append!(δρt_δV[ci], [-im*Δt*Hc_L[ci]*ρt[1]]) - append!(∂xδρt_δV[pi][ci], [-im*Δt*Hc_L[ci]*∂ρt_∂x[1][pi]]) - end - end - - for cj in 1:ctrl_num - Hc_L[cj] = liouville_commu(control_Hamiltonian[cj]) - end - - for ti in 2:tnum - expL = evolute(H[ti], Liouville_operator, γ, times, times[ti]) - ρt[ti] = expL * ρt[ti-1] - for pk in 1:para_num - ∂ρt_∂x[ti][pk] = -im * Δt * ∂H_L[pk] * ρt[ti] + expL * ∂ρt_∂x[ti-1][pk] - for ck in 1:ctrl_num - for tk in 1:ti - δρt_δV_first = popfirst!(δρt_δV[ck]) - ∂xδρt_δV_first = popfirst!(∂xδρt_δV[pk][ck]) - δρt_δV_tp = expL * δρt_δV_first - ∂xδρt_δV_tp = -im * Δt * ∂H_L[pk] * expL * δρt_δV_first + expL * ∂xδρt_δV_first - append!(δρt_δV[ck], [δρt_δV_tp]) - append!(∂xδρt_δV[pk][ck], [∂xδρt_δV_tp]) - end - δρt_δV_last = -im * Δt * Hc_L[ck] * ρt[ti] - ∂xδρt_δV_last = -im * Δt * Hc_L[ck] * ∂ρt_∂x[ti][pk] - append!(δρt_δV[ck], [δρt_δV_last]) - append!(∂xδρt_δV[pk][ck], [∂xδρt_δV_last]) - end - end - end - - ρt_T = ρt[end] |> vec2mat - ∂ρt_T = [(∂ρt_∂x[end][para] |> vec2mat) for para in 1:para_num] - F_T = QFIM(ρt_T, ∂ρt_T) - - if para_num == 1 - Lx = SLD_eig(ρt_T, ∂ρt_T[1]) - anti_commu = 2*Lx[1]*Lx[1] - δF = [[0.0 for i in 1:tnum] for ctrl in 1:ctrl_num] - for tm in 1:tnum - for cm in 1:ctrl_num - ∂ρt_T_δV = δρt_δV[cm][tm] |> vec2mat - ∂xδρt_T_δV = ∂xδρt_δV[1][cm][tm] |> vec2mat - term1 = tr(∂xδρt_T_δV*Lx[1]) - term2 = tr(∂ρt_T_δV*anti_commu) - δF[cm][tm] = ((2*term1-0.5*term2) |> real) - end - end - - elseif para_num == 2 - F_det = F_T[1,1] * F_T[2,2] - F_T[1,2] * F_T[2,1] - δF = [[0.0 for i in 1:tnum] for ctrl in 1:ctrl_num] - for tm in 1:tnum - for cm in 1:ctrl_num - for pm in 1:para_num - ∂ρt_T_δV = δρt_δV[cm][tm] |> vec2mat - ∂xδρt_T_δV = ∂xδρt_δV[pm][cm][tm] |> vec2mat - term1 = tr(∂xδρt_T_δV * Lx[pm]) - anti_commu = 2 * Lx[pm] * Lx[pm] - term2 = tr(∂ρt_T_δV * anti_commu) - δF[cm][tm] = δF[cm][tm] - ((2*term1-0.5*term2) |> real) - end - δF[cm][tm] = δF[cm][tm] / F_det - end - end - - else - δF = [[0.0 for i in 1:tnum] for ctrl in 1:ctrl_num] - for tm in 1:tnum - for cm in 1:ctrl_num - for pm in 1:para_num - ∂ρt_T_δV = δρt_δV[cm][tm] |> vec2mat - ∂xδρt_T_δV = ∂xδρt_δV[pm][cm][tm] |> vec2mat - term1 = tr(∂xδρt_T_δV * Lx[pm]) - anti_commu = 2 * Lx[pm] * Lx[pm] - term2 = tr(∂ρt_T_δV * anti_commu) - δF[cm][tm] = δF[cm][tm] - ((2*term1-0.5*term2) |> real) / (F_T[pm][pm] * F_T[pm][pm]) - end - end - end - end - ρt |> vec2mat |> filterZeros!, ∂ρt_∂x |> vec2mat |> filterZeros!, δF -end - -""" - propagate_analitical!() -""" -function propagate_analitical!(grape::GrapeControl) - grape.ρ, grape.∂ρ_∂x, δF = propagate_analitical(grape.freeHamiltonian, grape.Hamiltonian_derivative, grape.ρ_initial, - grape.Liouville_operator, grape.γ, grape.control_Hamiltonian, - grape.control_coefficients, grape.times ) - δF -end - -""" - propagate_ODEAD!() -""" -function propagate_ODEAD!(grape::GrapeControl) - H(p) = Htot(grape.freeHamiltonian, grape.control_Hamiltonian, p) - dt = grape.times[2] - grape.times[1] - tspan = (grape.times[1], grape.times[end]) - u0 = grape.ρ_initial - Γ = grape.Liouville_operator - f(u, p, t) = -im * (H(p)[t2Num(tspan[1], dt, t)] * u + u * H(p)[t2Num(tspan[1], dt, t)]) + - ([grape.γ[i] * (Γ[i] * u * Γ[i]' - (Γ[i]' * Γ[i] * u + u * Γ[i]' * Γ[i] )) for i in 1:length(Γ)] |> sum) - p = grape.control_coefficients - prob = ODEProblem(f, u0, tspan, p, saveat=dt) - u = solve(prob).u - du = Zygote.jacobian(solve(remake(prob, u0=u, p), sensealg=QuadratureAdjoint())) - u, du -end - -""" - propagate_L_ODE!() -""" -function propagate_L_ODE!(grape::GrapeControl) - H = Htot(grape.freeHamiltonian, grape.control_Hamiltonian, grape.control_coefficients) - Δt = grape.times[2] - grape.times[1] - tspan = (grape.times[1], grape.times[end]) - u0 = grape.ρ_initial |> vec - evo(p, t) = evolute(p[t2Num(tspan[1], Δt, t)], grape.Liouville_operator, grape.γ, grape.times, t2Num(tspan[1], Δt, t)) - f(u, p, t) = evo(p, t) * u - prob = DiscreteProblem(f, u0, tspan, H,dt=Δt) - ρt = solve(prob).u - ∂ρt_∂x = Vector{Vector{Vector{eltype(u0)}}}(undef, 1) - for para in 1:length(grape.Hamiltonian_derivative) - devo(p, t) = -1.0im * Δt * liouville_commu(grape.Hamiltonian_derivative[para]) * evo(p, t) - du0 = devo(H, tspan[1]) * u0 - g(du, p, t) = evo(p, t) * du + devo(p, t) * ρt[t2Num(tspan[1], Δt, t)] - dprob = DiscreteProblem(g, du0, tspan, H,dt=Δt) - ∂ρt_∂x[para] = solve(dprob).u - end - - grape.ρ, grape.∂ρ_∂x = ρt |> vec2mat, ∂ρt_∂x |> vec2mat -end - diff --git a/src/QuanEstimation.jl b/src/QuanEstimation.jl index 50650fe..bca0012 100644 --- a/src/QuanEstimation.jl +++ b/src/QuanEstimation.jl @@ -1,25 +1,3 @@ module QuanEstimation -using LinearAlgebra -using Zygote -using DifferentialEquations -using JLD -using Random -using SharedArrays -using Base.Threads -using SparseArrays - -include("AsymptoticBound/CramerRao.jl") -include("Common/common.jl") -include("Common/utils.jl") -include("Control/GRAPE.jl") -include("Control/PSO.jl") -include("Dynamics/Adam.jl") -include("Dynamics/dynamcs.jl") -# include("QuanResources/") - -export sigmax, sigmay, sigmaz, sigmam -export GrapeControl, evolute, propagate!, QFI, CFI, gradient_CFI!,gradient_QFIM! -export Run, RunODE, RunMixed, RunADAM, RunPSO - -end +end \ No newline at end of file diff --git a/src/opt_target.jl b/src/opt_target.jl new file mode 100644 index 0000000..f7a5864 --- /dev/null +++ b/src/opt_target.jl @@ -0,0 +1,20 @@ +abstract type AbstractOpt end + +struct Opt <: AbstractOpt + opt_target::Symbol +end + +struct ControlOpt<:Opt end +struct StateOpt<:Opt end +struct MeasurementOpt<:Opt end +struct CompOpt<:Opt end +struct StateControlOpt<:CompOpt end +struct StateMeasurementOpt<:CompOpt end +struct StateControlMeasurementOpt<:CompOpt end + +ControlOpt() = ControlOpt(:Copt) +StateOpt() = StateOpt(:Sopt) +MeasurementOpt() = MeasurementOpt(:Mopt) +StateControlOpt() = StateControlOpt(:SCopt) +StateMeasurementOpt() = StateMeasurementOpt(:SMopt) +StateControlMeasurementOpt() = StateControlMeasurementOpt(:SCMOpt) \ No newline at end of file diff --git a/src/struct.jl b/src/struct.jl new file mode 100644 index 0000000..d10af0d --- /dev/null +++ b/src/struct.jl @@ -0,0 +1,12 @@ +abstract type AbstractSystem end + +mutable struct QuanEstSystem{T<:AbstractOpt, F<:AbstractObjective, A<:AbstractAlogrithm, D<:AbstractDynamics, O<:AbstractOutput} <:AbstractSystem + opt::T + objective::F + algorithm::A + dynamics::D + output::O +end + + +