Merge pull request #19 from foldfelis-QO/revise

Revise

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumStateBase"
 uuid = "73ce9c4f-35d1-4161-b9e6-26915895bfed"
 authors = ["JingYu Ning <[email protected]> and contributors"]
-version = "1.0.3"
+version = "1.1.0"
 
 [deps]
 ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
@@ -11,7 +11,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [compat]
-ClassicalOrthogonalPolynomials = "0.4"
+ClassicalOrthogonalPolynomials = "0.5"
 DataDeps = "0.7"
 KernelDensity = "0.6"
 julia = "1.6"

notebook/test_new_api.jl

@@ -0,0 +1,31 @@
+### A Pluto.jl notebook ###
+# v0.16.1
+
+using Markdown
+using InteractiveUtils
+
+# ╔═╡ 120da254-1c28-11ec-1938-a55c8cde5727
+using Pkg; Pkg.develop(path=".."); Pkg.activate("..")
+
+# ╔═╡ c99f9fc5-0ab1-444a-80f6-b5241edf2f8f
+begin
+	using QuantumStateBase
+	using Plots
+	using LinearAlgebra
+end
+
+# ╔═╡ db69bedc-2803-4268-9f56-34e4a88a324e
+state = SqueezedThermalState(Complex{Float64}, ξ(0.5, π/4), 0.5, )
+
+# ╔═╡ cd19d00d-bf2b-479d-8c88-493311d97303
+d = gaussian_state_sampler(Float64, state, 4096)
+
+# ╔═╡ b9ea40f4-1ed1-4862-8858-be07f2c04633
+scatter(d[1, :], d[2, :])
+
+# ╔═╡ Cell order:
+# ╟─120da254-1c28-11ec-1938-a55c8cde5727
+# ╠═c99f9fc5-0ab1-444a-80f6-b5241edf2f8f
+# ╠═db69bedc-2803-4268-9f56-34e4a88a324e
+# ╠═cd19d00d-bf2b-479d-8c88-493311d97303
+# ╠═b9ea40f4-1ed1-4862-8858-be07f2c04633

src/QuantumStateBase.jl

@@ -3,7 +3,7 @@ module QuantumStateBase
     using LinearAlgebra
     using ClassicalOrthogonalPolynomials
 
-    const DIM = 70
+    const DIM = 100
 
     function __init__()
         register(DataDep(
@@ -26,7 +26,7 @@ module QuantumStateBase
     include("operator.jl")
     include("state.jl")
 
-    # wigner
+    # # wigner
     include("wigner_util.jl")
     include("wigner.jl")
 

src/basis.jl

@@ -37,14 +37,16 @@ julia> state = FockState(1, dim=100);
 julia> state = FockState(1, rep=StateMatrix);
 ```
 """
-FockState(n; T=ComplexF64, dim=DIM, rep=StateVector) = FockState(T, n, dim, rep)
+FockState(T::Type{<:Number}, n; dim=DIM, rep=StateVector) = FockState(T, n, dim, rep)
+FockState(n; dim=DIM, rep=StateVector) = FockState(ComplexF64, n, dim, rep)
 
 """
     NumberState(n; T=ComplexF64, dim=DIM, rep=StateVector)
 
 Exactly the same with `FockState`.
 """
-NumberState(n; T=ComplexF64, dim=DIM, rep=StateVector) = FockState(T, n, dim, rep)
+NumberState(T::Type{<:Number}, n; dim=DIM, rep=StateVector) = FockState(T, n, dim, rep)
+NumberState(n; dim=DIM, rep=StateVector) = FockState(ComplexF64, n, dim, rep)
 
 """
     VacuumState(; T=ComplexF64, dim=DIM, rep=StateVector)
@@ -64,7 +66,8 @@ julia> state = VacuumState(dim=100);
 julia> state = VacuumState(rep=StateMatrix);
 ```
 """
-VacuumState(; T=ComplexF64, dim=DIM, rep=StateVector) = FockState(T, 0, dim, rep)
+VacuumState(T::Type{<:Number}; dim=DIM, rep=StateVector) = FockState(T, 0, dim, rep)
+VacuumState(; dim=DIM, rep=StateVector) = FockState(ComplexF64, 0, dim, rep)
 
 """
     SinglePhotonState(; T=ComplexF64, dim=DIM, rep=StateVector)
@@ -84,4 +87,5 @@ julia> state = SinglePhotonState(dim=100);
 julia> state = SinglePhotonState(rep=StateMatrix);
 ```
 """
-SinglePhotonState(; T=ComplexF64, dim=DIM, rep=StateVector) = FockState(T, 1, dim, rep)
+SinglePhotonState(T::Type{<:Number}; dim=DIM, rep=StateVector) = FockState(T, 1, dim, rep)
+SinglePhotonState(; dim=DIM, rep=StateVector) = FockState(ComplexF64, 1, dim, rep)

src/operator.jl

@@ -27,7 +27,8 @@ Creation operator in matrix representation
 
 ``\\hat{a}^{\\dagger}``
 """
-Creation(; dim=DIM) = diagm(-1 => sqrt.(1:dim-1))
+Creation(T::Type{<:Number}; dim=DIM) = diagm(-1 => sqrt.(T.(1:dim-1)))
+Creation(; dim=DIM) = Creation(ComplexF64, dim=dim)
 
 """
     create!(state::AbstractState)
@@ -44,16 +45,16 @@ julia> vec(state) == vec(SinglePhotonState())
 true
 ```
 """
-function create!(state::StateVector{<:Number})
+function create!(state::StateVector{T}) where {T}
     dim = state.dim
-    state.v = Creation(dim=dim) * state.v
+    state.v = Creation(T, dim=dim) * state.v
 
     return state
 end
 
-function create!(state::StateMatrix{<:Number})
+function create!(state::StateMatrix{T}) where {T}
     dim = state.dim
-    𝐜 = Creation(dim=dim)
+    𝐜 = Creation(T, dim=dim)
     state.𝛒 = 𝐜 * state.𝛒 * 𝐜'
 
     return state
@@ -86,7 +87,8 @@ Annihilation operator in matrix representation
 
 ``\\hat{a}``
 """
-Annihilation(; dim=DIM) = diagm(1 => sqrt.(1:dim-1))
+Annihilation(T::Type{<:Number}; dim=DIM) = diagm(1 => sqrt.(T.(1:dim-1)))
+Annihilation(; dim=DIM) = Annihilation(ComplexF64, dim=dim)
 
 """
     annihilate!(state::AbstractState)
@@ -103,16 +105,16 @@ julia> vec(state) == vec(VacuumState())
 true
 ```
 """
-function annihilate!(state::StateVector{<:Number})
+function annihilate!(state::StateVector{T}) where {T}
     dim = state.dim
-    state.v = Annihilation(dim=dim) * state.v
+    state.v = Annihilation(T, dim=dim) * state.v
 
     return state
 end
 
-function annihilate!(state::StateMatrix{<:Number})
+function annihilate!(state::StateMatrix{T}) where {T}
     dim = state.dim
-    𝐚 = Annihilation(dim=dim)
+    𝐚 = Annihilation(T, dim=dim)
     state.𝛒 = 𝐚 * state.𝛒 * 𝐚'
 
     return state
@@ -149,14 +151,16 @@ Vector in polar coordinate for complex plane.
 
 ``v = r e^{-i\\theta}``
 """
-struct ComplexVec{T <: Real}
+struct ComplexVec{T<:Real}
     r::T
     θ::T
 end
 
-Base.show(io::IO, complexvec::ComplexVec{T}) where {T} = print(io, "ComplexVec{$T}($(complexvec.r)exp(-$(complexvec.θ)im))")
+function Base.show(io::IO, complexvec::ComplexVec{T}) where {T}
+    print(io, "ComplexVec{$T}($(complexvec.r)exp(-$(complexvec.θ)im))")
+end
 
-z(complexvec::ComplexVec{<:Real}) = complexvec.r * exp(-im * complexvec.θ)
+z(complexvec::ComplexVec) = complexvec.r * exp(-im * complexvec.θ)
 
 """
     α(r::Real, θ::Real)
@@ -171,7 +175,8 @@ julia> α(1.5, π/4)
 ComplexVec{Float64}(1.5exp(-0.7853981633974483im))
 ```
 """
-α(r::T, θ::T) where {T} = ComplexVec{T}(r, θ)
+α(T::Type{<:Real}, r::Real, θ::Real) = ComplexVec{T}(T(r), T(θ))
+α(r, θ) = α(Float64, r, θ)
 
 """
     ξ(r::Real, θ::Real)
@@ -195,10 +200,12 @@ Displacement operator in matrix representation
 
 ``\\hat{D}(\\alpha) = exp(\\alpha \\hat{a}^{\\dagger} - \\alpha^{*} \\hat{a})``
 """
-function Displacement(α::ComplexVec{<:Real}; dim=DIM)
-    return exp(z(α) * Creation(dim=dim) - z(α)' * Annihilation(dim=dim))
+function Displacement(T::Type{<:Number}, α::ComplexVec; dim=DIM)
+    return exp(z(α) * Creation(T, dim=dim) - z(α)' * Annihilation(T, dim=dim))
 end
 
+Displacement(α::ComplexVec; dim=DIM) = Displacement(ComplexF64, α, dim=dim)
+
 """
     displace!(state::AbstractState)
 
@@ -214,16 +221,16 @@ julia> vec(state) == vec(CoherentState(α(5., π/4)))
 true
 ```
 """
-function displace!(state::StateVector{<:Number}, α::ComplexVec{<:Real})
+function displace!(state::StateVector{T}, α::ComplexVec) where {T}
     dim = state.dim
-    state.v = Displacement(α, dim=dim) * state.v
+    state.v = Displacement(T, α, dim=dim) * state.v
 
     return state
 end
 
-function displace!(state::StateMatrix{<:Number}, α::ComplexVec{<:Real})
+function displace!(state::StateMatrix{T}, α::ComplexVec) where {T}
     dim = state.dim
-    𝐝 = Displacement(α, dim=dim)
+    𝐝 = Displacement(T, α, dim=dim)
     state.𝛒 = 𝐝 * state.𝛒 * 𝐝'
 
     return state
@@ -240,10 +247,12 @@ Squeezing operator in matrix representation
 
 ``\\hat{S}(\\xi) = exp(\\frac{1}{2} (\\xi^{*} \\hat{a}^{2} - \\xi \\hat{a}^{\\dagger 2}))``
 """
-function Squeezing(ξ::ComplexVec{<:Real}; dim=DIM)
-    return exp(0.5 * z(ξ)' * Annihilation(dim=dim)^2 - 0.5 * z(ξ) * Creation(dim=dim)^2)
+function Squeezing(T::Type{<:Number}, ξ::ComplexVec; dim=DIM)
+    return exp(0.5 * z(ξ)' * Annihilation(T, dim=dim)^2 - 0.5 * z(ξ) * Creation(T, dim=dim)^2)
 end
 
+Squeezing(ξ::ComplexVec; dim=DIM) = Squeezing(ComplexF64, ξ, dim=dim)
+
 """
     squeeze!(state::AbstractState)
 
@@ -259,16 +268,16 @@ julia> vec(state) == vec(SqueezedState(ξ(0.5, π/4)))
 true
 ```
 """
-function squeeze!(state::StateVector{<:Number}, ξ::ComplexVec{<:Real})
+function squeeze!(state::StateVector{T}, ξ::ComplexVec) where {T}
     dim = state.dim
-    state.v = Squeezing(ξ, dim=dim) * state.v
+    state.v = Squeezing(T, ξ, dim=dim) * state.v
 
     return state
 end
 
-function squeeze!(state::StateMatrix{<:Number}, ξ::ComplexVec{<:Real})
+function squeeze!(state::StateMatrix{T}, ξ::ComplexVec) where {T}
     dim = state.dim
-    𝐬 = Squeezing(ξ, dim=dim)
+    𝐬 = Squeezing(T, ξ, dim=dim)
     state.𝛒 = 𝐬 * state.𝛒 * 𝐬'
 
     return state
@@ -278,56 +287,40 @@ end
 # measurement #
 ###############
 
-# ##### for arb. statein θ-x quadrature coordinate #####
+##### for arb. statein θ-x quadrature coordinate #####
 
 # |θ, x⟩ = ∑ₙ |n⟩ ⟨n|θ, x⟩ = ∑ₙ ψₙ(θ, x) |n⟩
 # ⟨n|θ, x⟩ = ψₙ(θ, x) = exp(im n θ) (2/π)^(1/4) exp(-x^2) Hₙ(√2 x)/√(2^n n!)
-# coeff_ψₙ = (2/π)^(1/4)/√(2^n n!)
-# ψₙ = coeff_ψₙ(n) exp(im n θ) exp(-x^2) Hₙ(√2 x)
-calc_coeff_ψₙ(n::BigInt) = (2/π)^(1/4) / sqrt(2^n * factorial(n))
-COEFF_ψₙ = [calc_coeff_ψₙ(big(n)) for n in 0:(DIM-1)]
-
-function extend_coeff_ψₙ!(n::Integer)
-    while length(COEFF_ψₙ)-1 < n
-        push!(COEFF_ψₙ, calc_coeff_ψₙ(big(length(COEFF_ψₙ))))
-    end
-end
-
-function coeff_ψₙ(n::Integer)
-    (n < length(COEFF_ψₙ)) && (return COEFF_ψₙ[n+1])
-
-    return calc_coeff_ψₙ(big(n))
-end
-
 function ψₙ(n::Integer, θ::Real, x::Real)
-    return coeff_ψₙ(n) * exp(im * n * θ - x^2) * hermiteh(n, sqrt(2)x)
+    return (2/π)^(1/4) * exp(im*n*θ - x^2) * hermiteh(n, sqrt(2)x) / sqrt(2^n * factorial(n))
 end
 
 function 𝛑̂!(result::Matrix{<:Complex}, θ::Real, x::Real; dim=DIM)
-    view(result, :, 1) .= ψₙ.(0:dim-1, θ, x)
+    view(result, :, 1) .= ψₙ.(big(0):big(dim-1), θ, x)
     result .= view(result, :, 1) * view(result, :, 1)'
 
     return result
 end
 
-function 𝛑̂(θ::Real, x::Real; dim=DIM, T=ComplexF64)
+function 𝛑̂(T::Type{<:Complex}, θ::Real, x::Real; dim=DIM)
     result = Matrix{T}(undef, dim, dim)
-    U = T.parameters[1]
 
-    return 𝛑̂!(result, U(θ), U(x), dim=dim)
+    return 𝛑̂!(result, θ, x, dim=dim)
 end
 
-# ##### for Gaussian state in θ-x quadrature coordinate #####
+𝛑̂(θ::Real, x::Real; dim=DIM) = 𝛑̂(ComplexF64, θ, x, dim=dim)
+
+##### for Gaussian state in θ-x quadrature coordinate #####
 
 # π̂ₓ = (â exp(-im θ) + ↠exp(im θ)) / 2
 
 tr_mul(𝐚, 𝐛) = sum(𝐚[i, :]' * 𝐛[:, i] for i in 1:size(𝐚, 1))
-create_μ(state::StateMatrix) = tr_mul(Creation(dim=state.dim), state.𝛒)
-create²_μ(state::StateMatrix) = tr_mul(Creation(dim=state.dim)^2, state.𝛒)
-annihilate_μ(state::StateMatrix) = tr_mul(Annihilation(dim=state.dim), state.𝛒)
-annihilate²_μ(state::StateMatrix) = tr_mul(Annihilation(dim=state.dim)^2, state.𝛒)
-create_annihilate_μ(state::StateMatrix) = tr_mul(
-    Creation(dim=state.dim) * Annihilation(dim=state.dim),
+create_μ(state::StateMatrix{T}) where {T} = tr_mul(Creation(T, dim=state.dim), state.𝛒)
+create²_μ(state::StateMatrix{T}) where {T} = tr_mul(Creation(T, dim=state.dim)^2, state.𝛒)
+annihilate_μ(state::StateMatrix{T}) where {T} = tr_mul(Annihilation(T, dim=state.dim), state.𝛒)
+annihilate²_μ(state::StateMatrix{T}) where {T} = tr_mul(Annihilation(T, dim=state.dim)^2, state.𝛒)
+create_annihilate_μ(state::StateMatrix{T}) where {T} = tr_mul(
+    Creation(T, dim=state.dim) * Annihilation(T, dim=state.dim),
     state.𝛒
 )
 

src/quadrature_pdf.jl

@@ -9,17 +9,19 @@ real_tr_mul(𝐚, 𝐛) = sum(real(𝐚[i, :]' * 𝐛[:, i]) for i in 1:size(
 
 Quadrature prabability at point (θ, x)
 """
-function q_pdf(state::AbstractState, θ::Real, x::Real; T=Float64)
-    𝛑̂_res = Matrix{complex(T)}(undef, state.dim, state.dim)
+function q_pdf(T::Type{<:Real}, state::AbstractState, θ::Real, x::Real)
+    𝛑̂_res = Matrix{Complex{T}}(undef, state.dim, state.dim)
 
     return q_pdf!(𝛑̂_res, state, θ, x)
 end
 
-function q_pdf!(𝛑̂_res::Matrix{Complex{T}}, state::StateMatrix, θ::Real, x::Real) where {T}
-    return real_tr_mul(𝛑̂!(𝛑̂_res, T(θ), T(x), dim=state.dim), state.𝛒)
+q_pdf(state::AbstractState, θ::Real, x::Real) = q_pdf(Float64, state, θ, x)
+
+function q_pdf!(𝛑̂_res::AbstractMatrix, state::StateMatrix, θ::Real, x::Real)
+    return real_tr_mul(𝛑̂!(𝛑̂_res, θ, x, dim=state.dim), state.𝛒)
 end
 
-function q_pdf!(𝛑̂_res::Matrix{Complex{T}}, state::StateVector, args...; kwargs...) where {T}
+function q_pdf!(𝛑̂_res::AbstractMatrix, state::StateVector, args...; kwargs...)
     return q_pdf!(𝛑̂_res, StateMatrix(state), args...; kwargs...)
 end
 
@@ -28,13 +30,15 @@ end
 
 Quadrature prabability at points (θs, xs)
 """
-function q_pdf(state::AbstractState, θs, xs; T=Float64)
-    𝛑̂_res_vec = [Matrix{complex(T)}(undef, state.dim, state.dim) for _ in 1:Threads.nthreads()]
+function q_pdf(T::Type{<:Real}, state::AbstractState, θs, xs)
+    𝛑̂_res_vec = [Matrix{Complex{T}}(undef, state.dim, state.dim) for _ in 1:Threads.nthreads()]
     𝐩 = Matrix{T}(undef, length(θs), length(xs))
 
     return q_pdf!(𝛑̂_res_vec, 𝐩, state, θs, xs)
 end
 
+q_pdf(state::AbstractState, θs, xs) = q_pdf(Float64, state, θs, xs)
+
 function q_pdf!(𝛑̂_res_vec::Vector{Matrix{Complex{T}}}, 𝐩::Matrix{T}, state::StateMatrix, θs, xs) where {T}
     @sync for (j, x) in enumerate(xs)
         for (i, θ) in enumerate(θs)