From 86c988e87ad30e37f62a535301f72955e1185d78 Mon Sep 17 00:00:00 2001 From: JingYu Ning Date: Thu, 29 Jul 2021 15:19:10 +0800 Subject: [PATCH] fallback to vanilla rejection method --- Project.toml | 4 +- docs/src/get_started.md | 2 +- src/sampler.jl | 90 +++++++++-------------------------------- test/sampler.jl | 26 +++--------- 4 files changed, 27 insertions(+), 95 deletions(-) diff --git a/Project.toml b/Project.toml index 1661d54..2f0805f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuantumStateBase" uuid = "73ce9c4f-35d1-4161-b9e6-26915895bfed" authors = ["JingYu Ning and contributors"] -version = "1.0.0" +version = "1.0.1" [deps] ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" @@ -12,8 +12,8 @@ Mmap = "a63ad114-7e13-5084-954f-fe012c677804" [compat] ClassicalOrthogonalPolynomials = "0.4" -KernelDensity = "0.6" DataDeps = "0.7" +KernelDensity = "0.6" julia = "1.6" [extras] diff --git a/docs/src/get_started.md b/docs/src/get_started.md index 8969e0f..72ee1eb 100644 --- a/docs/src/get_started.md +++ b/docs/src/get_started.md @@ -247,7 +247,7 @@ julia> heatmap(θs, xs, ps') ## Quantum state sampler Here, we can sample points from quadrature probability density function of the quantum state. -The sampler is implemented by special adaptive rejection method. +The sampler is implemented by rejection method. ```julia julia> points = rand(state, 4096); diff --git a/src/sampler.jl b/src/sampler.jl index f125ab8..5b87b59 100644 --- a/src/sampler.jl +++ b/src/sampler.jl @@ -64,101 +64,49 @@ function ranged_rand(range::Tuple{T, T}) where {T <: Number} return range[1] + (range[2]-range[1]) * rand(T) end -function reject!(new_points, θ_range, gen_point, p, g, c) - Threads.@threads for i in 1:size(new_points, 2) - new_points[:, i] .= gen_point() - while ( - !(θ_range[1] ≤ new_points[1, i] ≤ θ_range[2]) || - p(Threads.threadid(), new_points[:, i]...) / g(new_points[:, i]...) < c - ) - new_points[:, i] .= gen_point() - end - end - - return new_points -end - """ state_sampler( state::AbstractState, n::Integer; - warm_up_n=128, batch_size=64, c=0.9, θ_range=(0., 2π), x_range=(-10., 10.), - show_log=false + c=1, θ_range=(0., 2π), x_range=(-10., 10.), ) Random points sampled from quadrature probability density function of Gaussian `state`. * `state`: Quantum state. * `n`: N points. -* `warm_up_n`: N points sampled from uniform random and accepted by rejection method. -* `batch_size`: Adapt `g` for every `batch_size` points. * `θ_range`: Sampling range of θ. * `x_range`: Sampling range of x. """ function state_sampler( state::AbstractState, n::Integer; - warm_up_n=128, batch_size=64, c=0.9, θ_range=(0., 2π), x_range=(-10., 10.), - show_log=false + c=1, θ_range=(0., 2π), x_range=(-10., 10.), ) sampled_points = Matrix{Float64}(undef, 2, n) - 𝛑̂_res_vec = [Matrix{complex(Float64)}(undef, state.dim, state.dim) for _ in 1:Threads.nthreads()] + 𝛑̂_res_vec = [ + Matrix{complex(Float64)}(undef, state.dim, state.dim) + for _ in 1:Threads.nthreads() + ] - return state_sampler!( - sampled_points, 𝛑̂_res_vec, - state, - warm_up_n, batch_size, c, θ_range, x_range, - show_log - ) + return state_sampler!(sampled_points, 𝛑̂_res_vec, state, c, θ_range, x_range) end +# rejection method: \frac{f(x)}{c * g(x)} ≥ u +# here, say g(x) is a uniform distribution, and c = 1 function state_sampler!( sampled_points::AbstractMatrix{T}, 𝛑̂_res_vec::Vector{Matrix{Complex{T}}}, - state::StateMatrix, - warm_up_n::Integer, batch_size::Integer, c::Real, θ_range, x_range, - show_log::Bool + state::StateMatrix, c::Real, θ_range, x_range ) where {T} n = size(sampled_points, 2) - p = (thread_id, θ, x) -> q_pdf!(𝛑̂_res_vec[thread_id], state, θ, x) - - show_log && @info "Warm up" - warm_up_n = n < warm_up_n ? n : warm_up_n - warm_up_points = view(sampled_points, :, 1:warm_up_n) - gen_rand_point = () -> [ranged_rand(θ_range), ranged_rand(x_range)] - kde_result = kde((ranged_rand(n, θ_range), ranged_rand(n, x_range))) - g = (θ, x) -> pdf(kde_result, θ, x) - reject!(warm_up_points, θ_range, gen_rand_point, p, g, c) - - show_log && @info "Start to generate data" - batch = div(n-warm_up_n, batch_size) - for b in 1:batch - ref_range = 1:(warm_up_n+(b-1)*batch_size) - ref_points = view(sampled_points, :, ref_range) - new_range = (warm_up_n+(b-1)*batch_size+1):(warm_up_n+b*batch_size) - new_points = view(sampled_points, :, new_range) - - h = KernelDensity.default_bandwidth((ref_points[1, :], ref_points[2, :])) - gen_point_from_g = () -> ref_points[:, rand(ref_range)] + randn(2)./h - kde_result = kde((ref_points[1, :], ref_points[2, :]), bandwidth=h) - g = (θ, x) -> pdf(kde_result, θ, x) - reject!(new_points, θ_range, gen_point_from_g, p, g, c) - - show_log && @info "progress: $b/$(batch+1)" - end - rem_n = rem(n-warm_up_n, batch_size) - if rem_n > 0 - ref_range = 1:(n-rem_n) - ref_points = view(sampled_points, :, ref_range) - new_range = (n-rem_n+1):n - new_points = view(sampled_points, :, new_range) - - h = KernelDensity.default_bandwidth((ref_points[1, :], ref_points[2, :])) - gen_point_from_g = () -> ref_points[:, rand(ref_range)] + randn(2)./h - kde_result = kde((ref_points[1, :], ref_points[2, :]), bandwidth=h) - g = (θ, x) -> pdf(kde_result, θ, x) - reject!(new_points, θ_range, gen_point_from_g, p, g, c) - end - show_log && @info "progress: $(batch+1)/$(batch+1)" + sampled_points[1, :] .= sort!(ranged_rand(n, θ_range)) + + Threads.@threads for i in 1:n + p = (thread_id, x) -> q_pdf!(𝛑̂_res_vec[thread_id], state, sampled_points[1, i], x) - sampled_points .= sampled_points[:, sortperm(sampled_points[1, :])] + sampled_points[2, i] = ranged_rand(x_range) + while (p(Threads.threadid(), sampled_points[2, i]) < c*rand()) + sampled_points[2, i] = ranged_rand(x_range) + end + end return sampled_points end diff --git a/test/sampler.jl b/test/sampler.jl index 92af596..40eed42 100644 --- a/test/sampler.jl +++ b/test/sampler.jl @@ -48,20 +48,10 @@ end single_point_pdf = (θ, x) -> q_pdf(state, θ, x) @test single_point_pdf.(θs, xs') ≈ ground_truth_pdf - n = 4096 - data = state_sampler(state, n) + n = 10000 + @time data = state_sampler(state, n) sampled_pdf = pdf(kde((data[1, :], data[2, :])), θs, xs) - @test sum(abs.(sampled_pdf .- ground_truth_pdf)) / n < 5e-2 - - n = 4096-1 - data = state_sampler(state, n) - sampled_pdf = pdf(kde((data[1, :], data[2, :])), θs, xs) - @test sum(abs.(sampled_pdf .- ground_truth_pdf)) / n < 5e-2 - - n = 4100 - data = state_sampler(state, n, warm_up_n=100, batch_size=97) - sampled_pdf = pdf(kde((data[1, :], data[2, :])), θs, xs) - @test sum(abs.(sampled_pdf .- ground_truth_pdf)) / n < 5e-2 + @test sum(abs.(sampled_pdf .- ground_truth_pdf)) / n < 5e-4 end @testset "wrapping" begin @@ -76,14 +66,8 @@ end state = SinglePhotonState() @test size(rand(state)) == (2, 1) - @test size(rand(state, 4096, show_log=false)) == (2, 4096) - @test size(rand(state, 4100, warm_up_n=97, show_log=false)) == (2, 4100) - @test size(rand(state, 4100, warm_up_n=97, batch_size=100, show_log=false)) == (2, 4100) - @test size(rand(state, 100, warm_up_n=200, batch_size=100, show_log=false)) == (2, 100) + @test size(rand(state, 4096)) == (2, 4096) state = SinglePhotonState(rep=StateMatrix) @test size(rand(state)) == (2, 1) - @test size(rand(state, 4096, show_log=false)) == (2, 4096) - @test size(rand(state, 4100, warm_up_n=97, show_log=false)) == (2, 4100) - @test size(rand(state, 4100, warm_up_n=97, batch_size=100, show_log=false)) == (2, 4100) - @test size(rand(state, 100, warm_up_n=200, batch_size=100, show_log=false)) == (2, 100) + @test size(rand(state, 4096)) == (2, 4096) end