Merge pull request #86 from QuanEstimation/unit-test

Add PythonCall extension & Unit test

Author: hmyuuu

.github/workflows/CI.yml

@@ -0,0 +1,46 @@
+name: CI
+on:
+  push:
+    branches:
+      - unit-test
+    tags: '*'
+  pull_request:
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          # - '1.8'
+          - '1.9'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+        # include:
+        #   - os: windows-latest
+        #     version: '1'
+        #     arch: x64
+        #   - os: macOS-latest
+        #     version: '1'
+        #     arch: x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v1
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v2
+        with:
+          files: lcov.info

Project.toml

@@ -18,18 +18,24 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReinforcementLearning = "158674fc-8238-5cab-b5ba-03dfc80d1318"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
+[weakdeps]
+PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+
+[extensions]
+QuanEstimationPyExt = "PythonCall"
+
 [compat]
 BoundaryValueDiffEq = "2"
-Convex = "0.15"
 Distributions = "0.25"
-Flux = "0.12, 0.13"
+Flux = "0.12, 0.13, 0.14"
 Interpolations = "0.13, 0.14"
 IntervalSets = "0.5, 0.7"
 JLD2 = "0.4"
@@ -41,7 +47,6 @@ StableRNGs = "1"
 StatsBase = "0.33"
 Trapz = "2"
 Zygote = "0.6"
-julia = "1.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -1,6 +1,9 @@
 # QuanEstimation.jl
 
-[![Dev](https://img.shields.io/badge/docs-stable-blue.svg)](https://quanestimation.github.io/QuanEstimation/)
+[![][docs-img]][docs-url]
+[![][action-img]][action-url]
+[![][codecov-img]][codecov-url]
+
 
 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, and the optimizations with respect to probe states, controls or measurements, as well as comprehensive optimizations in quantum metrology. Futhermore, QuanEstimation can generate not only optimal quantum parameter estimation schemes, but also adaptive measurement schemes.
 
@@ -19,3 +22,12 @@ If you use QuanEstimation in your research, please cite the following paper:
 [1] M. Zhang, H.-M. Yu, H. Yuan, X. Wang, R. Demkowicz-Dobrzański, and J. Liu, 
 QuanEstimation: An open-source toolkit for quantum parameter estimation, 
 [Phys. Rev. Research **4**, 043057 (2022).](https://doi.org/10.1103/PhysRevResearch.4.043057)
+
+
+
+[action-img]: https://github.com/QuanEstimation/QuanEstimation.jl/actions/workflows/CI.yml/badge.svg
+[action-url]: https://github.com/QuanEstimation/QuanEstimation.jl/actions
+[codecov-img]: https://codecov.io/gh/QuanEstimation/QuanEstimation.jl/branch/master/graph/badge.svg
+[codecov-url]: https://codecov.io/gh/QuanEstimation/QuanEstimation.jl?branch=test-codecov
+[docs-img]: https://img.shields.io/badge/docs-stable-blue.svg
+[docs-url]: https://quanestimation.github.io/QuanEstimation/

examples/CramerRao_bounds.jl

@@ -21,7 +21,7 @@ M = [M1, M2]
 # time length for the evolution
 tspan = range(0., 50., length=2000)
 # dynamics
-rho, drho = QuanEstimation.expm(tspan, rho0, H0, dH, decay)
+rho, drho = QuanEstimation.expm(tspan, rho0, H0, dH, decay=decay)
 # calculation of the CFI and QFI
 Im, F = Float64[], Float64[]
 for ti in 2:length(tspan)

ext/QuanEstimationPyExt.jl

@@ -0,0 +1,45 @@
+module QuanEstimationPyExt
+using PythonCall
+using QuanEstimation
+np = pyimport("numpy")
+# QuanEstimation.Htot(H0::PyArray, Hc::PyList, ctrl) = QuanEstimation.Htot(Matrix(H0), [Matrix(hc) for hc in Hc], ctrl)
+QuanEstimation.ControlOpt(ctrl::PyList, ctrl_bound::PyList, seed::Py) = QuanEstimation.ControlOpt(pyconvert(Vector{Vector{Float64}}, ctrl), pyconvert(Vector{Float64}, ctrl_bound), seed)
+QuanEstimation.Lindblad(
+	H0::PyArray,
+	dH::PyList,
+	Hc::PyList,
+	ctrl::PyList,
+	ρ0::PyArray,
+	tspan::PyArray,
+	decay_opt::PyList,
+	γ::PyList;
+	kwargs...,
+) = QuanEstimation.Lindblad(
+	pyconvert(Matrix{ComplexF64}, H0),
+	pyconvert(Vector{Matrix}, dH),
+	pyconvert(Vector{Matrix}, Hc),
+	pyconvert(Vector{Vector{Float64}}, ctrl),
+	pyconvert(Matrix{ComplexF64}, ρ0),
+	pyconvert(Vector, tspan),
+	pyconvert(Vector{Matrix}, decay_opt),
+	pyconvert(Vector, γ);
+	kwargs...,
+)
+QuanEstimation.Lindblad(
+	H0::PyArray,
+	dH::PyList,
+	ψ0::PyList,
+	tspan::PyArray,
+	decay_opt::PyList,
+	γ::PyList;
+	kwargs...,
+) = QuanEstimation.Lindblad(
+    pyconvert(Matrix{ComplexF64}, H0),
+	pyconvert(Vector{Matrix}, dH),
+	pyconvert(Vector{ComplexF64}, ψ0),
+	pyconvert(Vector, tspan),
+	pyconvert(Vector{Matrix}, decay_opt),
+	pyconvert(Vector, γ);
+	kwargs...,
+)
+end

src/Algorithm/DDPG.jl

@@ -3,8 +3,6 @@ using Flux: glorot_normal, glorot_uniform
 using StableRNGs
 using IntervalSets
 
-Base.copyto!(dest::ReinforcementLearning.NeuralNetworkApproximator, src::ReinforcementLearning.NeuralNetworkApproximator) =
-           Flux.loadparams!(dest.model, Flux.params(src))
 
 function Base.rsplit( v, l::Int)
     u = reshape(v,l,length(v)÷l)

src/ObjectiveFunc/AsymptoticBound/AnalogCramerRao.jl

@@ -1,3 +1,4 @@
+using Convex
 function decomposition(A)
     C = bunchkaufman(A; check = false)
     R = sqrt(Array(C.D)) * C.U'C.P

src/Parameterization/Lindblad/LindbladDynamics.jl

@@ -51,18 +51,18 @@ function liouville_dissip_py(A::Array{T}) where {T<:Complex}
 end
 
 function dissipation(
-    Γ::AbstractVector,
+    Γ::V,
     γ::Vector{R},
     t::Int = 1,
-) where {T<:Complex,R<:Real}
+) where {V<:AbstractVector,R<:Real}
     [γ[i] * liouville_dissip(Γ[i]) for i = 1:length(Γ)] |> sum
 end
 
 function dissipation(
-    Γ::AbstractVector,
+    Γ::V,
     γ::Vector{Vector{R}},
     t::Int = 1,
-) where {T<:Complex,R<:Real}
+) where {V<:AbstractVector,R<:Real}
     [γ[i][t] * liouville_dissip(Γ[i]) for i = 1:length(Γ)] |> sum
 end
 
@@ -81,20 +81,21 @@ function liouvillian(
     -1.0im * freepart + dissp
 end
 
-function Htot(H0::Matrix{T}, Hc::Vector{Matrix{T}}, ctrl) where {T<:Complex,R}
+function Htot(H0::T, Hc::V, ctrl) where {T<:AbstractArray, V<:AbstractVector}
     [H0] .+ ([ctrl[i] .* [Hc[i]] for i = 1:length(ctrl)] |> sum)
 end
 
 function Htot(
-    H0::Matrix{T},
-    Hc::Vector{Matrix{T}},
+    H0::T,
+    Hc::V,
     ctrl::Vector{R},
-) where {T<:Complex,R<:Real}
-    H0 + ([ctrl[i] * Hc[i] for i = 1:length(ctrl)] |> sum)
+) where {T<:AbstractArray,V<:AbstractVector, R<:Real}
+    H0 + ([ctrl[i] .* [Hc[i]] for i = 1:length(ctrl)] |> sum)
 end
 
-function Htot(H0::Vector{Matrix{T}}, Hc::Vector{Matrix{T}}, ctrl) where {T<:Complex}
-    H0 + ([ctrl[i] .* [Hc[i]] for i = 1:length(ctrl)] |> sum)
+
+function Htot(H0::V1, Hc::V2, ctrl) where {V1,V2<:AbstractVector}
+    H0 + ([ctrl[i] * Hc[i] for i = 1:length(ctrl)] |> sum)
 end
 
 function expL(H, decay_opt, γ, dt, tj = 1)
@@ -118,7 +119,7 @@ When applied to the case of single parameter.
 function expm(
     tspan::AbstractVector,
     ρ0::AbstractMatrix,
-    H0::AbstractMatrix,
+    H0::AbstractVecOrMat,
     dH::AbstractMatrix;
     decay::Union{AbstractVector, Missing}=missing,
     Hc::Union{AbstractVector, Missing}=missing,
@@ -197,7 +198,7 @@ The dynamics of a density matrix is of the form  ``\partial_t\rho=-i[H,\rho]+\su
 function expm(
     tspan::AbstractVector,
     ρ0::AbstractMatrix,
-    H0::AbstractMatrix,
+    H0::AbstractVecOrMat,
     dH::AbstractVector;
     decay::Union{AbstractVector, Missing}=missing,
     Hc::Union{AbstractVector, Missing}=missing,
@@ -241,7 +242,8 @@ function expm(
     para_num = length(dH)
     ctrl_num = length(Hc)
     ctrl_interval = ((length(tspan) - 1) / length(ctrl0[1])) |> Int
-    ctrl = [repeat(ctrl0[i], 1, ctrl_interval) |> transpose |> vec for i = 1:ctrl_num]
+    ## TODO reconstruct repeat feature
+    ctrl = [repeat(ctrl0[i], outer=ctrl_interval) for i = 1:ctrl_num]
 
     H = Htot(H0, Hc, ctrl)
     dH_L = [liouville_commu(dH[i]) for i = 1:para_num]
@@ -271,7 +273,7 @@ end
 function expm_py(
     tspan,
     ρ0::AbstractMatrix,
-    H0::AbstractMatrix,
+    H0::AbstractVecOrMat,
     dH::AbstractMatrix,
     Hc::AbstractVector,
     decay_opt::AbstractVector,
@@ -304,7 +306,7 @@ end
 function expm_py(
     tspan,
     ρ0::AbstractMatrix,
-    H0::AbstractMatrix,
+    H0::AbstractVecOrMat,
     dH::AbstractVector,
     decay_opt::AbstractVector,
     γ,
@@ -353,7 +355,7 @@ expm(tspan, ρ0, H0, dH, decay, Hc, ctrl) =
 function secondorder_derivative(
     tspan::AbstractVector,
     ρ0::AbstractMatrix,
-    H0::AbstractMatrix,
+    H0::AbstractVecOrMat,
     dH::AbstractVector,
     dH_∂x::AbstractVector,
     decay_opt::Vector{Matrix{T}},
@@ -399,7 +401,7 @@ When applied to the case of single parameter.
 function ode(
     tspan::AbstractVector,
     ρ0::AbstractMatrix,
-    H0::AbstractMatrix,
+    H0::AbstractVecOrMat,
     dH::AbstractMatrix;
     decay::Union{AbstractVector, Missing}=missing,
     Hc::Union{AbstractVector, Missing}=missing,
@@ -480,7 +482,7 @@ The dynamics of a density matrix is of the form  ``\partial_t\rho=-i[H,\rho]+\su
 function ode(
     tspan::AbstractVector,
     ρ0::AbstractMatrix,
-    H0::AbstractMatrix,
+    H0::AbstractVecOrMat,
     dH::AbstractVector;
     decay::Union{AbstractVector, Missing}=missing,
     Hc::Union{AbstractVector, Missing}=missing,
@@ -559,7 +561,7 @@ ode(tspan, ρ0, H0, dH, decay, Hc, ctrl) =
 function ode_py(
     tspan,
     ρ0::AbstractMatrix,
-    H0::AbstractMatrix,
+    H0::AbstractVecOrMat,
     dH::AbstractMatrix,
     Hc::AbstractVector,
     Γ::AbstractVector,
@@ -589,7 +591,7 @@ end
 function ode_py(
     tspan,
     ρ0::AbstractMatrix,
-    H0::AbstractMatrix,
+    H0::AbstractVecOrMat,
     dH::AbstractVector,
     Γ::AbstractVector,
     γ,

src/Parameterization/Lindblad/LindbladWrapper.jl

@@ -175,7 +175,7 @@ function _ini_measurement!(opt::Mopt_LinearComb, dim::Int; eps=GLOBAL_EPS)
 	end
 
 	if ismissing(B)
-		opt.B = [rand(opt.rng, length(POVM_basis)) for _ in 1:M_num]
+		opt.B = [rand(opt.rng, length(opt.POVM_basis)) for _ in 1:M_num]
 	end
 end
 

src/Parameterization/Parameterization.jl

@@ -21,12 +21,12 @@ abstract type Ode <: AbstractDynamicsProblemType end
 # check if the dynamics are with noise
 isNoisy(::noiseless) = false
 isNoisy(::noisy) = true
-isNoisy(dynamics::AbstractDynamics) = dynamics.noise_type |> eval |> isNoisy
+isNoisy(dynamics::AbstractDynamics) = dynamics.noise_type |> isNoisy
 
 # check if the dynamics are in control
 isCtrl(::free) = false
 isCtrl(::controlled) = true
-isCtrl(dynamics::AbstractDynamics) = dynamics.ctrl_type |> eval |> isCtrl
+isCtrl(dynamics::AbstractDynamics) = dynamics.ctrl_type |> isCtrl
 
 include("Lindblad/Lindblad.jl")
 include("Kraus/Kraus.jl")

src/QuanEstimation.jl

@@ -1,4 +1,12 @@
 module QuanEstimation
+export ControlOpt, ControlMeasurementOpt, CMopt, StateMeasurementOpt, SMopt, StateControlMeasurementOpt, SCMopt, opt_target, Htot
+export isCtrl, isNoisy
+export Lindblad
+export Output
+export autoGRAPE
+export QFIM_obj
+export QuanEstSystem
+export solve
 using Random
 using LinearAlgebra
 using Zygote
@@ -7,7 +15,6 @@ using DelimitedFiles
 using StatsBase
 using Flux
 using ReinforcementLearning
-using Convex
 using SCS
 using BoundaryValueDiffEq
 using Trapz