Merge pull request #12 from meerkat321/methods_conflect

Implement value test for pdf

Project.toml

@@ -19,10 +19,9 @@ QuantumStateBase = "1.2"
 julia = "1.6"
 
 [extras]
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 QuantumStateBase = "73ce9c4f-35d1-4161-b9e6-26915895bfed"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 
 [targets]
-test = ["Test", "QuantumStateBase", "KernelDensity"]
-
+test = ["Test", "QuantumStateBase", "LsqFit"]

src/pdf.jl

@@ -3,14 +3,25 @@ export
     qpdf!
 
 """
-    qpdf([T=Float64], d::QuantumStateBHD, θ::Real, x::Real)
-    qpdf([T=Float64], d::QuantumStateBHD, θs::AbstractRange, xs::AbstractRange)
     qpdf([T=Float64], d::GaussianStateBHD, θ::Real, x::Real)
     qpdf([T=Float64], d::GaussianStateBHD, θs::AbstractRange, xs::AbstractRange)
+    qpdf([T=Float64], d::QuantumStateBHD, θ::Real, x::Real)
+    qpdf([T=Float64], d::QuantumStateBHD, θs::AbstractRange, xs::AbstractRange)
+    qpdf([T=Float64], ρ::AbstractArray, θ::Real, x::Real)
+    qpdf([T=Float64], ρ::AbstractArray, θs::AbstractRange, xs::AbstractRange)
+
+Quadrature prabability in intensity-to-measurement-phase quadrature coordinate.
+
+## Arguments
+
+* `T`: Data type, default as Float64.
+* `state`: State can be `GaussianStateBHD` distribution, `QuantumStateBHD` distribution or density matrix.
+* `θ`: Measurement phase, can be `Real` or `AbstractRange`.
+* `x`: Intensity in quadrature coordinate, can be `Real` or `AbstractRange`.
 
-Quadrature prabability at point (θ, x) or points (θs, xs)
+``p(\\rho, \\theta, x) = tr(\\hat{\\Pi}(\\theta, x) \\rho)``
 """
-qpdf(d, θ, x) = qpdf(Float64, d, θ, x)
+qpdf(state, θ, x) = qpdf(Float64, state, θ, x)
 
 function qpdf(T::Type{<:Real}, d::GaussianStateBHD, θ::Real, x::Real)
     μ = QuantumStateDistributions.mean(d, θ)
@@ -47,30 +58,20 @@ end
 
 qpdf(T::Type{<:Real}, d::QuantumStateBHD, θ, x) = qpdf(T, d.ρ, θ, x)
 
-"""
-    qpdf([T=Float64], ρ::AbstractArray, θ::Real, x::Real)
-    qpdf([T=Float64], ρ::AbstractArray, θs::AbstractRange, xs::AbstractRange)
-
-Quadrature prabability in intensity-to-measurement-phase quadrature coordinate.
-
-``p(\\rho, \\theta, x) = tr(\\hat{\\Pi}(\\theta, x) \\rho)``
-"""
-qpdf(ρ, θ, x) = qpdf(Float64, ρ, θ, x)
-
-function qpdf(T::Type{<:Real}, ρ, θ::Real, x::Real)
+function qpdf(T::Type{<:Real}, ρ::AbstractMatrix, θ::Real, x::Real)
     dim = size(ρ, 1)
     𝛑̂_res = Matrix{Complex{T}}(undef, dim, dim)
 
     return qpdf!(𝛑̂_res, ρ, θ, x)
 end
 
-function qpdf!(𝛑̂_res::AbstractMatrix, ρ::AbstractArray, θ::Real, x::Real)
+function qpdf!(𝛑̂_res::AbstractMatrix, ρ::AbstractMatrix, θ::Real, x::Real)
     dim = size(ρ, 1)
 
     return real_tr_mul(𝛑̂!(𝛑̂_res, θ, x, dim=dim), ρ)
 end
 
-function qpdf(T::Type{<:Real}, ρ, θs::AbstractRange, xs::AbstractRange)
+function qpdf(T::Type{<:Real}, ρ::AbstractMatrix, θs::AbstractRange, xs::AbstractRange)
     dim = size(ρ, 1)
     𝛑̂_res_vec = [Matrix{Complex{T}}(undef, dim, dim) for _ in 1:Threads.nthreads()]
     𝐩 = Matrix{T}(undef, length(θs), length(xs))
@@ -80,7 +81,7 @@ end
 
 function qpdf!(
     𝛑̂_res_vec::AbstractVector{Matrix{Complex{T}}}, 𝐩::Matrix{T},
-    ρ::AbstractArray, θs::AbstractRange, xs::AbstractRange
+    ρ::AbstractMatrix, θs::AbstractRange, xs::AbstractRange
 ) where {T}
     @sync for (j, x) in enumerate(xs)
         for (i, θ) in enumerate(θs)

test/pdf.jl

@@ -1,14 +1,101 @@
-@testset "qpdf" begin
+using LsqFit
+
+@testset "qpdf for GaussianStateBHD" begin
     r, θ, dim = 0.8, 2π, 100
     ρ = SqueezedState(r, θ, Matrix, dim=dim)
     d = GaussianStateBHD(ρ)
 
-    θs = LinRange(0, 2π, 101)
-    xs = LinRange(-3, 3, 101)
+    m = 10
+    n = 1000
+    θs = LinRange(0, 2π, m)
+    xs = LinRange(-3, 3, n)
+
+    g_pdf = Vector{Float64}(undef, n)
+
+    for i in 1:n
+        g_pdf[i] = qpdf(d, π/3, xs[i])
+    end
+
+    g_pdfs = qpdf(d, θs, xs)
+
+    moments = [QuantumStateDistributions.mean(d, π/3), QuantumStateDistributions.std(d, π/3)]
+    @. model(x, p) = 1 / (p[2] * √(2π)) * exp(-(x - p[1])^2 / (2 * p[2]^2))
+    p0 = [0.5, 0.5]
+
+    g_fit = curve_fit(model, xs, g_pdf, p0)
+    @test coef(g_fit) ≈ moments
+
+    for i in 1:m
+        g_fits = curve_fit(model, xs, g_pdfs[i, :], p0)
+        moments = [QuantumStateDistributions.mean(d, θs[i]), QuantumStateDistributions.std(d, θs[i])]
+
+        @test coef(g_fits) ≈ moments
+    end
+end
+
+@testset "qpdf for QuantumStateBHD" begin
+    r, θ, dim = 0.8, 2π, 100
+    ρ = SqueezedState(r, θ, Matrix, dim=dim)
+    q = QuantumStateBHD(ρ)
+    d = GaussianStateBHD(ρ)
+
+    m = 10
+    n = 1000
+    θs = LinRange(0, 2π, m)
+    xs = LinRange(-3, 3, n)
+
+    q_pdf = Vector{Float64}(undef, n)
+
+    for i in 1:n
+        q_pdf[i] = qpdf(q, π/3, xs[i])
+    end
+
+    q_pdfs = qpdf(q, θs, xs)
+
+    moments = [QuantumStateDistributions.mean(d, π/3), QuantumStateDistributions.std(d, π/3)]
+    @. model(x, p) = 1 / (p[2] * √(2π)) * exp(-(x - p[1])^2 / (2 * p[2]^2))
+    p0 = [0.5, 0.5]
+
+    q_fit = curve_fit(model, xs, q_pdf, p0)
+    @test coef(q_fit) ≈ moments
+
+    for i in 1:m
+        q_fits = curve_fit(model, xs, q_pdfs[i, :], p0)
+        moments = [QuantumStateDistributions.mean(d, θs[i]), QuantumStateDistributions.std(d, θs[i])]
+
+        @test coef(q_fits) ≈ moments
+    end
+end
+
+@testset "pdf for density matrix" begin
+    r, θ, dim = 0.8, 2π, 100
+    ρ = SqueezedState(r, θ, Matrix, dim=dim)
+    d = GaussianStateBHD(ρ)
+
+    m = 10
+    n = 1000
+    θs = LinRange(0, 2π, m)
+    xs = LinRange(-3, 3, n)
+
+    ρ_pdf = Vector{Float64}(undef, n)
+
+    for i in 1:n
+        ρ_pdf[i] = qpdf(ρ, π/3, xs[i])
+    end
+
+    ρ_pdfs = qpdf(ρ, θs, xs)
+
+    moments = [QuantumStateDistributions.mean(d, π/3), QuantumStateDistributions.std(d, π/3)]
+    @. model(x, p) = 1 / (p[2] * √(2π)) * exp(-(x - p[1])^2 / (2 * p[2]^2))
+    p0 = [0.5, 0.5]
+
+    ρ_fit = curve_fit(model, xs, ρ_pdf, p0)
+    @test coef(ρ_fit) ≈ moments
+
+    for i in 1:m
+        ρ_fits = curve_fit(model, xs, ρ_pdfs[i, :], p0)
+        moments = [QuantumStateDistributions.mean(d, θs[i]), QuantumStateDistributions.std(d, θs[i])]
 
-    @test size(qpdf(d, θs, xs)) == (101, 101)
-    @test length(qpdf(d, π/2, 1)) == 1
-    @test size(qpdf(Float32, d, θs, xs)) == (101, 101)
-    @test length(qpdf(Float32, d ,π/3, 2)) == 1
-    @test size(qpdf(d, LinRange(0, 2π, 10), LinRange(-3, 3, 10))) == (10, 10)
+        @test coef(ρ_fits) ≈ moments
+    end
 end