From ada53cd855a072f8ff1bbde61102da8d4669a63d Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Sat, 15 Jan 2022 12:28:58 -1000 Subject: [PATCH 1/9] switch logvol and logwt calculations to use quadratic trapezoid rule --- src/step.jl | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/step.jl b/src/step.jl index 8041bb17..48e8d88d 100644 --- a/src/step.jl +++ b/src/step.jl @@ -10,11 +10,14 @@ function step(rng, model, sampler::Nested; kwargs...) v_dead = vs[:, idx_dead] # update weight using trapezoidal rule - logvol = log1mexp(-1 / sampler.nactive) - logwt = logl_dead + logvol + # logvol = log1mexp(-1 / sampler.nactive) + logvol = -sampler.dlv + # logdvol = log(exp(0.5) * sum(exp(logvol + sample.dlv, -logvol)))) + logdvol = logvol - log(sampler.nactive) - log(2) + # logwt = logl_dead + logvol + logwt = logl_dead + logdvol # sample a new live point without bounds - point = rand(rng, eltype(us), sampler.ndims) bound = Bounds.fit(Bounds.NoBounds, us) proposal = Proposals.Rejection() u, v, ll, nc = proposal(rng, v_dead, logl_dead, bound, model) @@ -29,7 +32,8 @@ function step(rng, model, sampler::Nested; kwargs...) logz = logwt h = logl_dead - logz logzerr = sqrt(h / sampler.nactive) - logvol -= 1 / sampler.nactive + # # logvol -= 1 / sampler.nactive + # logvol -= sampler.dlv sample = (u = u_dead, v = v_dead, logwt = logwt, logl = logl_dead) state = (it = 1, ncall = ncall, us = us, vs = vs, logl = logl, logl_dead = logl_dead, @@ -91,8 +95,13 @@ function step(rng, model, sampler, state; kwargs...) ncall = state.ncall + nc since_update += nc - # update weight - logwt = state.logvol + logl_dead + # update weight using trapezoidal rule + # logvol = log1mexp(-1 / sampler.nactive) + logvol = state.logvol - sampler.dlv + # logdvol = log(exp(0.5) * sum(exp(logvol + sample.dlv, -logvol)))) + logdvol = logvol - log(sampler.nactive) - log(2) + # logwt = state.logvol + logl_dead + logwt = logaddexp(state.logl_dead, logl_dead) + logdvol # update evidence and information logz = logaddexp(state.logz, logwt) From 0965663f9eace71cd587b9a12e67f7f875da5620 Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Sat, 15 Jan 2022 13:04:49 -1000 Subject: [PATCH 2/9] incorporate quad. trap. rule to logzerr calculation --- src/step.jl | 20 ++++++-------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/src/step.jl b/src/step.jl index 48e8d88d..0e7d8cab 100644 --- a/src/step.jl +++ b/src/step.jl @@ -10,11 +10,8 @@ function step(rng, model, sampler::Nested; kwargs...) v_dead = vs[:, idx_dead] # update weight using trapezoidal rule - # logvol = log1mexp(-1 / sampler.nactive) logvol = -sampler.dlv - # logdvol = log(exp(0.5) * sum(exp(logvol + sample.dlv, -logvol)))) logdvol = logvol - log(sampler.nactive) - log(2) - # logwt = logl_dead + logvol logwt = logl_dead + logdvol # sample a new live point without bounds @@ -30,10 +27,8 @@ function step(rng, model, sampler::Nested; kwargs...) # update evidence and information logz = logwt - h = logl_dead - logz - logzerr = sqrt(h / sampler.nactive) - # # logvol -= 1 / sampler.nactive - # logvol -= sampler.dlv + h = exp(logl_dead - logz + logdvol) * logl_dead - logz + logzerr = sqrt(h * sampler.dlv) sample = (u = u_dead, v = v_dead, logwt = logwt, logl = logl_dead) state = (it = 1, ncall = ncall, us = us, vs = vs, logl = logl, logl_dead = logl_dead, @@ -96,19 +91,16 @@ function step(rng, model, sampler, state; kwargs...) since_update += nc # update weight using trapezoidal rule - # logvol = log1mexp(-1 / sampler.nactive) logvol = state.logvol - sampler.dlv - # logdvol = log(exp(0.5) * sum(exp(logvol + sample.dlv, -logvol)))) logdvol = logvol - log(sampler.nactive) - log(2) - # logwt = state.logvol + logl_dead logwt = logaddexp(state.logl_dead, logl_dead) + logdvol # update evidence and information logz = logaddexp(state.logz, logwt) - h = (exp(logwt - logz) * logl_dead + - exp(state.logz - logz) * (state.h + state.logz) - logz) - logzerr = h ≥ 0 ? sqrt(h / sampler.nactive) : NaN - logvol = state.logvol - 1 / sampler.nactive + logzterm = exp(state.logl_dead - logz + logdvol) * state.logl_dead + + exp(logl_dead - logz + logdvol) * logl_dead + h = logzterm + exp(state.logz - logz) * (state.h + state.logz) - logz + logzerr = h ≥ 0 ? sqrt(state.logzerr^2 + (h - state.h) * sampler.dlv) : NaN ## prepare returns sample = (u = u_dead, v = v_dead, logwt = logwt, logl = logl_dead) From 03a906cbf20ce9f5815fdb4d5151f0d3b4653d96 Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Sat, 15 Jan 2022 13:09:57 -1000 Subject: [PATCH 3/9] relax constraint for pathological h --- src/step.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/step.jl b/src/step.jl index 0e7d8cab..b3ec5c75 100644 --- a/src/step.jl +++ b/src/step.jl @@ -100,7 +100,7 @@ function step(rng, model, sampler, state; kwargs...) logzterm = exp(state.logl_dead - logz + logdvol) * state.logl_dead + exp(logl_dead - logz + logdvol) * logl_dead h = logzterm + exp(state.logz - logz) * (state.h + state.logz) - logz - logzerr = h ≥ 0 ? sqrt(state.logzerr^2 + (h - state.h) * sampler.dlv) : NaN + logzerr = sqrt(state.logzerr^2 + (h - state.h) * sampler.dlv) ## prepare returns sample = (u = u_dead, v = v_dead, logwt = logwt, logl = logl_dead) From e57f330718f9bf66cef0c4402c8f483b2bfdd793 Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Sat, 15 Jan 2022 14:40:49 -1000 Subject: [PATCH 4/9] add quadratic trapezoidal rule to add_live_points --- src/step.jl | 44 +++++++++++++++++++++++++++----------------- 1 file changed, 27 insertions(+), 17 deletions(-) diff --git a/src/step.jl b/src/step.jl index b3ec5c75..4f443242 100644 --- a/src/step.jl +++ b/src/step.jl @@ -31,8 +31,8 @@ function step(rng, model, sampler::Nested; kwargs...) logzerr = sqrt(h * sampler.dlv) sample = (u = u_dead, v = v_dead, logwt = logwt, logl = logl_dead) - state = (it = 1, ncall = ncall, us = us, vs = vs, logl = logl, logl_dead = logl_dead, - logz = logz, logzerr = logzerr, h = h, logvol = logvol, + state = (it = 1, ncall = ncall, us = us, vs = vs, logl = logl, + logl_dead = logl_dead, logz = logz, logzerr = logzerr, h = h, logvol = logvol, since_update = since_update, has_bounds = false, active_bound = nothing) return sample, state @@ -104,8 +104,8 @@ function step(rng, model, sampler, state; kwargs...) ## prepare returns sample = (u = u_dead, v = v_dead, logwt = logwt, logl = logl_dead) - state = (it = it, ncall = ncall, us = state.us, vs = state.vs, logl = state.logl, logl_dead = logl_dead, - logz = logz, logzerr = logzerr, h = h, logvol = logvol, + state = (it = it, ncall = ncall, us = state.us, vs = state.vs, logl = state.logl, + logl_dead = logl_dead, logz = logz, logzerr = logzerr, h = h, logvol = logvol, since_update = since_update, has_bounds = has_bounds, active_bound = active_bound) return sample, state @@ -208,36 +208,46 @@ end # add remaining live points to `samples` function add_live_points(samples, model, sampler, state) - logvol = -state.it / sampler.nactive - log(sampler.nactive) - + prev_logvol = state.logvol prev_logz = state.logz prev_h = state.h + prev_logzerr = state.logzerr + prev_logl_dead = state.logl_dead - local logl, logz, h, logzerr + local logz, h, logzerr, logl_dead, logvol N = length(samples) @inbounds for (i, idx) in enumerate(eachindex(state.logl)) # get new point u = state.us[:, idx] v = state.vs[:, idx] - logl = state.logl[idx] + logl_dead = state.logl[idx] + + # update weight using trapezoidal rule + logvol = state.logvol + log1p(-i / (sampler.nactive + 1)) + dlv = prev_logvol - logvol + logdvol = logvol + log(exp(dlv) - 1) - log(2) + logwt = logaddexp(prev_logl_dead, logl_dead) + logdvol - # update sampler - logwt = logvol + logl + # update evidence and information logz = logaddexp(prev_logz, logwt) - h = (exp(logwt - logz) * logl + - exp(prev_logz - logz) * (prev_h + prev_logz) - logz) - logzerr = sqrt(h / sampler.nactive) + logzterm = exp(prev_logl_dead - logz + logdvol) * prev_logl_dead + + exp(logl_dead - logz + logdvol) * logl_dead + h = logzterm + exp(prev_logz - logz) * (prev_h + prev_logz) - logz + logzerr = sqrt(max(zero(h), prev_logzerr^2 + (h - prev_h) * dlv)) + prev_logvol = logvol prev_logz = logz prev_h = h + prev_logzerr = logzerr + prev_logl_dead = logl_dead - sample = (u = u, v = v, logwt = logwt, logl = logl) + sample = (u = u, v = v, logwt = logwt, logl = logl_dead) save!!(samples, sample, N + i, model, sampler) end - state = (it = state.it + sampler.nactive, us = state.us, vs = state.vs, logl = logl, - logz = logz, logzerr = logzerr, logvol = logvol, - since_update = state.since_update, has_bounds = state.has_bounds, active_bound = state.active_bound) + state = (it = state.it + sampler.nactive, us = state.us, vs = state.vs, logl = state.logl, + logl_dead = logl_dead, logz = logz, logzerr = logzerr, h = h, logvol = logvol, + since_update = state.since_update, has_bounds = state.has_bounds, active_bound = state.active_bound) return samples, state end From 2f993f6bca9354bc757f53f2f4615fead94639a4 Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Sat, 15 Jan 2022 16:24:27 -1000 Subject: [PATCH 5/9] rerun benchmarks with new integrator --- bench/Manifest.toml | 6 ++++++ bench/Project.toml | 1 + bench/sampling.jl | 22 ++++++++++------------ bench/sampling_results.csv | 10 +++++----- 4 files changed, 22 insertions(+), 17 deletions(-) diff --git a/bench/Manifest.toml b/bench/Manifest.toml index db05500b..264c4d95 100644 --- a/bench/Manifest.toml +++ b/bench/Manifest.toml @@ -665,6 +665,12 @@ git-tree-sha1 = "39c9f91521de844bad65049efd4f9223e7ed43f9" uuid = "171d559e-b47b-412a-8079-5efa626c420e" version = "0.1.14" +[[deps.StableRNGs]] +deps = ["Random", "Test"] +git-tree-sha1 = "3be7d49667040add7ee151fefaf1f8c04c8c8276" +uuid = "860ef19b-820b-49d6-a774-d7a799459cd3" +version = "1.0.0" + [[deps.Static]] deps = ["IfElse"] git-tree-sha1 = "b4912cd034cdf968e06ca5f943bb54b17b97793a" diff --git a/bench/Project.toml b/bench/Project.toml index e77475ec..46d4a7fe 100644 --- a/bench/Project.toml +++ b/bench/Project.toml @@ -5,5 +5,6 @@ NestedSamplers = "41ceaf6f-1696-4a54-9b49-2e7a9ec3782e" ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" diff --git a/bench/sampling.jl b/bench/sampling.jl index f4b26a64..f069b3b3 100644 --- a/bench/sampling.jl +++ b/bench/sampling.jl @@ -2,27 +2,25 @@ using BenchmarkTools using CSV using NestedSamplers -using PythonCall +using ProgressLogging +using StableRNGs using Statistics using StatsBase +rng = StableRNG(112358) rows = [] -dims = [2, 4, 8, 16, 32] -for D in dims - +dims = 2 .^ (1:5) +@progress for D in dims model, true_lnZ = Models.CorrelatedGaussian(D) - splr = Nested(D, 50D; proposal=Proposals.Slice(), bounds=Bounds.Ellipsoid) + splr = Nested(D, 50D; proposal=Proposals.Slice(slices=5), bounds=Bounds.Ellipsoid) # run once to extract values from state, also precompile - ch, state = sample(model, splr; dlogz=0.01, chain_type=Array) - lnZ = state.logz - lnZstd = state.logzerr + ch, state = sample(rng, model, splr; dlogz=0.01, chain_type=Array) + t = @belapsed sample($rng, $model, $splr; dlogz=0.01, chain_type=Array) - tt = @belapsed sample($model, $splr; dlogz=0.01, chain_type=Array) + dlnZ = state.logz - true_lnZ - dlnZ = abs(true_lnZ - lnZ) - - row = (; library="NestedSamplers.jl", D, t=median(tt), lnZ, lnZstd, dlnZ) + row = (; library="NestedSamplers.jl", D, t=t, lnZ=state.logz, lnZstd=state.logzerr, dlnZ) @info "$row" push!(rows, row) end diff --git a/bench/sampling_results.csv b/bench/sampling_results.csv index 501c3bad..bebe228a 100644 --- a/bench/sampling_results.csv +++ b/bench/sampling_results.csv @@ -1,6 +1,6 @@ library,D,t,lnZ,lnZstd,dlnZ -NestedSamplers.jl,2,0.129553849,-3.7464515508498475,0.1241802259632353,0.012655386209909025 -NestedSamplers.jl,4,0.95785889,-6.233094208166701,0.1334147665362398,0.045180940536692304 -NestedSamplers.jl,8,4.539512133,-10.575110892409716,0.13712909419237015,0.12434640904846894 -NestedSamplers.jl,16,17.816551031,-18.752425012323915,0.13595239531892792,0.3202050261597691 -NestedSamplers.jl,32,150.735562133,-34.75372182270856,0.14125000020140985,0.8322655713364 +NestedSamplers.jl,2,0.201037605,-3.776611302476892,0.1259604503086115,-0.017504365417135404 +NestedSamplers.jl,4,0.954831354,-6.23589557023661,0.1343478011972437,-0.0479823026066013 +NestedSamplers.jl,8,4.369661922,-10.752241751423023,0.13698920791533648,-0.30147726806177566 +NestedSamplers.jl,16,28.311188809,-18.63596853727137,0.1415574543503746,-0.20374855110722478 +NestedSamplers.jl,32,249.568858136,-34.867406263641186,0.1436942445755534,-0.945950012269023 From 8e82f7c42f929d460e76f2bc2be7ef297b08cc18 Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Sat, 15 Jan 2022 20:28:36 -1000 Subject: [PATCH 6/9] add dynesty benchmark to set expectations --- bench/CondaPkg.toml | 2 ++ bench/sampling.jl | 25 +++++++++++++++++++++++-- bench/sampling_results.csv | 12 +++++++----- docs/Project.toml | 1 - docs/src/benchmarks.md | 16 +++++++++++----- test/models.jl | 2 +- 6 files changed, 44 insertions(+), 14 deletions(-) create mode 100644 bench/CondaPkg.toml diff --git a/bench/CondaPkg.toml b/bench/CondaPkg.toml new file mode 100644 index 00000000..c4c0f35f --- /dev/null +++ b/bench/CondaPkg.toml @@ -0,0 +1,2 @@ +[deps] +dynesty = "" diff --git a/bench/sampling.jl b/bench/sampling.jl index f069b3b3..b33182ed 100644 --- a/bench/sampling.jl +++ b/bench/sampling.jl @@ -3,15 +3,18 @@ using BenchmarkTools using CSV using NestedSamplers using ProgressLogging +using PythonCall using StableRNGs using Statistics using StatsBase +dy = pyimport("dynesty") + rng = StableRNG(112358) rows = [] -dims = 2 .^ (1:5) -@progress for D in dims +dims = 2 .^ (1:4) +@progress name="NestedSamplers.jl" for D in dims model, true_lnZ = Models.CorrelatedGaussian(D) splr = Nested(D, 50D; proposal=Proposals.Slice(slices=5), bounds=Bounds.Ellipsoid) # run once to extract values from state, also precompile @@ -25,6 +28,24 @@ dims = 2 .^ (1:5) push!(rows, row) end +@progress name="dynesty" for D in dims[begin:end-1] + model, true_lnZ = Models.CorrelatedGaussian(D) + splr = dy.NestedSampler( + model.prior_transform_and_loglikelihood.loglikelihood, + model.prior_transform_and_loglikelihood.prior_transform, + D; nlive=50D, bound="single", sample="slice", slices=5 + ) + t = @elapsed splr.run_nested(dlogz=0.01) + res = splr.results + lnZ = PyArray(res["logz"], copy=false)[end] + lnZstd = PyArray(res["logzerr"], copy=false)[end] + dlnZ = lnZ - true_lnZ + + row = (; library="dynesty", D, t, lnZ, lnZstd, dlnZ) + @info "$row" + push!(rows, row) +end + path = joinpath(@__DIR__, "sampling_results.csv") CSV.write(path, rows) @info "output saved to $path" diff --git a/bench/sampling_results.csv b/bench/sampling_results.csv index bebe228a..fff55ce2 100644 --- a/bench/sampling_results.csv +++ b/bench/sampling_results.csv @@ -1,6 +1,8 @@ library,D,t,lnZ,lnZstd,dlnZ -NestedSamplers.jl,2,0.201037605,-3.776611302476892,0.1259604503086115,-0.017504365417135404 -NestedSamplers.jl,4,0.954831354,-6.23589557023661,0.1343478011972437,-0.0479823026066013 -NestedSamplers.jl,8,4.369661922,-10.752241751423023,0.13698920791533648,-0.30147726806177566 -NestedSamplers.jl,16,28.311188809,-18.63596853727137,0.1415574543503746,-0.20374855110722478 -NestedSamplers.jl,32,249.568858136,-34.867406263641186,0.1436942445755534,-0.945950012269023 +NestedSamplers.jl,2,0.018802904,-3.776611302476892,0.1259604503086115,-0.017504365417135404 +NestedSamplers.jl,4,0.112532803,-6.173820528139972,0.1334784546731703,0.014092739490036976 +NestedSamplers.jl,8,1.136714685,-10.527414660722169,0.13686799596023871,-0.07665017736092139 +NestedSamplers.jl,16,11.023224754,-19.108389493075787,0.14231162869377487,-0.6761695069116413 +dynesty,2,13.67880283,-3.733850115101514,0.17078345040236537,0.025256821958242526 +dynesty,4,87.696596209,-6.076658702382752,0.1860338795913046,0.11125456524725674 +dynesty,8,600.509202249,-10.214865171723309,0.1932735763541268,0.23589931163793842 diff --git a/docs/Project.toml b/docs/Project.toml index 37811a52..37497d4f 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,7 +7,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" diff --git a/docs/src/benchmarks.md b/docs/src/benchmarks.md index f4aad572..a68f804c 100644 --- a/docs/src/benchmarks.md +++ b/docs/src/benchmarks.md @@ -65,11 +65,14 @@ This benchmark uses [`Models.CorrelatedGaussian`](@ref) and simply measures the ### Timing ```@example sample-benchmark -using CSV, DataFrames, NestedSamplers, Plots # hide +using CSV, DataFrames, NestedSamplers, StatsPlots # hide benchdir = joinpath(dirname(pathof(NestedSamplers)), "..", "bench") # hide results = DataFrame(CSV.File(joinpath(benchdir, "sampling_results.csv"))) # hide -plot(results.D, results.t, label="NestedSamplers.jl", marker=:o, yscale=:log10, # hide - ylabel="runtime (s)", xlabel="prior dimension", leg=:topleft, ylims=(1e-2, 1e4)) # hide +groups = groupby(results, :library) # hide +@df groups[1] plot(:D, :t, label="NestedSamplers.jl", marker=:o, yscale=:log10, # hide + ylabel="runtime (s)", xlabel="prior dimension", leg=:topleft, # hide + ylims=(1e-2, 1e4), markerstrokecolor=:auto) # hide +@df groups[2] plot!(:D, :t, label="dynesty", marker=:triangle, markerstrokecolor=:auto) # hide ``` ### Accuracy @@ -77,7 +80,10 @@ plot(results.D, results.t, label="NestedSamplers.jl", marker=:o, yscale=:log10, The following shows the Bayesian evidence estmiate as compared to the true value ```@example sample-benchmark -plot(results.D, results.dlnZ, yerr=results.lnZstd, label="NestedSamplers.jl", # hide - marker=:o, ylabel="ΔlnZ", xlabel="prior dimension", leg=:topleft) # hide +@df groups[1] plot(:D, :dlnZ, yerr=:lnZstd, label="NestedSamplers.jl", # hide + marker=:o, ylabel="ΔlnZ", xlabel="prior dimension", # hide + leg=:topleft, markerstrokecolor=:auto) # hide +@df groups[2] plot!(:D, :dlnZ, yerr=:lnZstd, label="dynesty", # hide + marker=:o, markerstrokecolor=:auto) # hide hline!([0.0], c=:black, ls=:dash, alpha=0.7, label="") # hide ``` diff --git a/test/models.jl b/test/models.jl index d35159fa..abc01c11 100644 --- a/test/models.jl +++ b/test/models.jl @@ -12,7 +12,7 @@ const test_props = [ const MAXZSCORES = Dict(zip( Iterators.product(test_bounds, test_props), - [4, 3, 9, 8, 6, 3, 5, 7, 4, 4] # rwalk is bad... + [4, 3, 9, 8, 6, 3, 5, 12, 4, 4] # rwalk is bad... )) function test_logz(measured, actual, error, bound, proposal) From e0245e6cc16eebc58f33c247d312c91b9585fed3 Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Sun, 16 Jan 2022 08:16:11 -1000 Subject: [PATCH 7/9] run full dynesty set --- bench/sampling.jl | 4 ++-- bench/sampling_results.csv | 15 ++++++++------- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/bench/sampling.jl b/bench/sampling.jl index b33182ed..c1919080 100644 --- a/bench/sampling.jl +++ b/bench/sampling.jl @@ -19,7 +19,7 @@ dims = 2 .^ (1:4) splr = Nested(D, 50D; proposal=Proposals.Slice(slices=5), bounds=Bounds.Ellipsoid) # run once to extract values from state, also precompile ch, state = sample(rng, model, splr; dlogz=0.01, chain_type=Array) - t = @belapsed sample($rng, $model, $splr; dlogz=0.01, chain_type=Array) + t = @elapsed sample(rng, model, splr; dlogz=0.01, chain_type=Array) dlnZ = state.logz - true_lnZ @@ -28,7 +28,7 @@ dims = 2 .^ (1:4) push!(rows, row) end -@progress name="dynesty" for D in dims[begin:end-1] +@progress name="dynesty" for D in dims model, true_lnZ = Models.CorrelatedGaussian(D) splr = dy.NestedSampler( model.prior_transform_and_loglikelihood.loglikelihood, diff --git a/bench/sampling_results.csv b/bench/sampling_results.csv index fff55ce2..cc260fef 100644 --- a/bench/sampling_results.csv +++ b/bench/sampling_results.csv @@ -1,8 +1,9 @@ library,D,t,lnZ,lnZstd,dlnZ -NestedSamplers.jl,2,0.018802904,-3.776611302476892,0.1259604503086115,-0.017504365417135404 -NestedSamplers.jl,4,0.112532803,-6.173820528139972,0.1334784546731703,0.014092739490036976 -NestedSamplers.jl,8,1.136714685,-10.527414660722169,0.13686799596023871,-0.07665017736092139 -NestedSamplers.jl,16,11.023224754,-19.108389493075787,0.14231162869377487,-0.6761695069116413 -dynesty,2,13.67880283,-3.733850115101514,0.17078345040236537,0.025256821958242526 -dynesty,4,87.696596209,-6.076658702382752,0.1860338795913046,0.11125456524725674 -dynesty,8,600.509202249,-10.214865171723309,0.1932735763541268,0.23589931163793842 +NestedSamplers.jl,2,0.040124793,-3.776611302476892,0.1259604503086115,-0.017504365417135404 +NestedSamplers.jl,4,0.264282685,-5.905758259445603,0.1307813273766663,0.28215500818440553 +NestedSamplers.jl,8,1.248085493,-10.509974669347862,0.1365781040925085,-0.05921018598661476 +NestedSamplers.jl,16,10.817955738,-17.873550793433207,0.14071085448722837,0.5586691927309388 +dynesty,2,13.998151217,-3.6544882444005133,0.1757184239472305,0.10461869265924317 +dynesty,4,89.716491323,-6.21727948417353,0.18737407387137883,-0.02936621654352134 +dynesty,8,614.704773432,-10.486069868517422,0.1955724257483754,-0.03530538515617465 +dynesty,16,4002.620846303,-18.38592241896873,0.19906552168471234,0.04629756719541689 From 8e7fc6e42fdbd69f7af5b371cab591a98615d1a5 Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Mon, 25 Jul 2022 22:06:37 -1000 Subject: [PATCH 8/9] update CI --- .github/workflows/ci.yml | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ff504ca4..56739269 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -19,7 +19,7 @@ jobs: fail-fast: false matrix: version: - - '1.3' + - '1.6' - '1' - 'nightly' os: @@ -34,16 +34,7 @@ jobs: with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: actions/cache@v1 - env: - cache-name: cache-artifacts - with: - path: ~/.julia/artifacts - key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} - restore-keys: | - ${{ runner.os }}-test-${{ env.cache-name }}- - ${{ runner.os }}-test- - ${{ runner.os }}- + - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 From 5667ae0012372c267064e5cea9b9b2f863b95003 Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Mon, 25 Jul 2022 22:24:50 -1000 Subject: [PATCH 9/9] switch to literate --- docs/Project.toml | 2 + docs/make.jl | 20 ++++++--- .../correlated.md => examples/correlated.jl | 42 +++++++----------- .../examples/eggbox.md => examples/eggbox.jl | 44 +++++++------------ .../examples/shells.md => examples/shells.jl | 44 +++++++------------ 5 files changed, 64 insertions(+), 88 deletions(-) rename docs/src/examples/correlated.md => examples/correlated.jl (70%) rename docs/src/examples/eggbox.md => examples/eggbox.jl (71%) rename docs/src/examples/shells.md => examples/shells.jl (72%) diff --git a/docs/Project.toml b/docs/Project.toml index 37497d4f..aefbd4f1 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,9 +4,11 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" +NestedSamplers = "41ceaf6f-1696-4a54-9b49-2e7a9ec3782e" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" diff --git a/docs/make.jl b/docs/make.jl index 44cd73bf..57fb4d6e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,15 @@ using Documenter using NestedSamplers +using Literate: markdown + +# preprocess tutorials using Literate +srcdir = abspath(joinpath(@__DIR__, "..", "examples")) +outdir = abspath(joinpath(@__DIR__, "src")) + +examples = map(Iterators.filter(endswith(".jl"), readdir(srcdir; join=true))) do file + markdown(file, outdir) + return replace(basename(file), ".jl" => ".md") +end DocMeta.setdocmeta!( NestedSamplers, @@ -12,11 +22,7 @@ makedocs( sitename = "NestedSamplers.jl", pages = [ "Home" => "index.md", - "Examples" => [ - "Gaussian Shells" => "examples/shells.md", - "Correlated Gaussian" => "examples/correlated.md", - "Eggbox" => "examples/eggbox.md", - ], + "Examples" => examples, "Benchmarks" => "benchmarks.md", "API/Reference" => "api.md" ], @@ -30,4 +36,8 @@ makedocs( ] ) +# clean up markdown files generated by Literate +map(f -> rm(joinpath(outdir, f)), examples) + +# CI only: deploy docs deploydocs(repo = "github.com/TuringLang/NestedSamplers.jl.git", push_preview=true, devbranch="main") diff --git a/docs/src/examples/correlated.md b/examples/correlated.jl similarity index 70% rename from docs/src/examples/correlated.md rename to examples/correlated.jl index 040533a8..8826b19f 100644 --- a/docs/src/examples/correlated.md +++ b/examples/correlated.jl @@ -1,3 +1,4 @@ +#= # Correlated Gaussian This example will explore a highly-correlated Gaussian using [`Models.CorrelatedGaussian`](@ref). This model uses a conjuage Gaussian prior, see the docstring for the mathematical definition. @@ -8,28 +9,22 @@ For this example, you'll need to add the following packages ```julia julia>]add Distributions MCMCChains Measurements NestedSamplers StatsBase StatsPlots ``` - -```@setup correlated +=# using AbstractMCMC using Random AbstractMCMC.setprogress!(false) Random.seed!(8452) -``` -## Define model +# ## Define model -```@example correlated using NestedSamplers -# set up a 4-dimensional Gaussian +## set up a 4-dimensional Gaussian D = 4 -model, logz = Models.CorrelatedGaussian(D) -nothing; # hide -``` +model, logz = Models.CorrelatedGaussian(D); -let's take a look at a couple of parameters to see what the likelihood surface looks like +# let's take a look at a couple of parameters to see what the likelihood surface looks like -```@example correlated using StatsPlots θ1 = range(-1, 1, length=1000) @@ -44,42 +39,35 @@ heatmap( xlabel="θ1", ylabel="θ2" ) -``` -## Sample +# ## Sample -```@example correlated using MCMCChains using StatsBase -# using single Ellipsoid for bounds -# using Gibbs-style slicing for proposing new points +## using single Ellipsoid for bounds +## using Gibbs-style slicing for proposing new points sampler = Nested(D, 50D; bounds=Bounds.Ellipsoid, proposal=Proposals.Slice() ) names = ["θ_$i" for i in 1:D] chain, state = sample(model, sampler; dlogz=0.01, param_names=names) -# resample chain using statistical weights +## resample chain using statistical weights chain_resampled = sample(chain, Weights(vec(chain[:weights])), length(chain)); -nothing # hide -``` -## Results +# ## Results -```@example correlated chain_resampled -``` -```@example correlated +# + corner(chain_resampled) -``` -```@example correlated +# + using Measurements logz_est = state.logz ± state.logzerr diff = logz_est - logz println("logz: $logz") println("estimate: $logz_est") println("diff: $diff") -nothing # hide -``` diff --git a/docs/src/examples/eggbox.md b/examples/eggbox.jl similarity index 71% rename from docs/src/examples/eggbox.md rename to examples/eggbox.jl index 23ad14f2..51fb2ae8 100644 --- a/docs/src/examples/eggbox.md +++ b/examples/eggbox.jl @@ -1,3 +1,4 @@ +#= # Eggbox This example will explore the classic eggbox function using [`Models.Eggbox`](@ref). @@ -8,26 +9,20 @@ For this example, you'll need to add the following packages ```julia julia>]add Distributions MCMCChains Measurements NestedSamplers StatsBase StatsPlots ``` - -```@setup eggbox +=# using AbstractMCMC using Random AbstractMCMC.setprogress!(false) Random.seed!(8452) -``` -## Define model +# ## Define model -```@example eggbox using NestedSamplers -model, logz = Models.Eggbox() -nothing; # hide -``` +model, logz = Models.Eggbox(); -let's take a look at a couple of parameters to see what the log-likelihood surface looks like +# let's take a look at a couple of parameters to see what the log-likelihood surface looks like -```@example eggbox using StatsPlots x = range(0, 1, length=1000) @@ -41,47 +36,40 @@ heatmap( xlabel="x", ylabel="y", ) -``` -## Sample +# ## Sample -```@example eggbox using MCMCChains using StatsBase -# using multi-ellipsoid for bounds -# using default rejection sampler for proposals +## using multi-ellipsoid for bounds +## using default rejection sampler for proposals sampler = Nested(2, 500) chain, state = sample(model, sampler; dlogz=0.01, param_names=["x", "y"]) -# resample chain using statistical weights +## resample chain using statistical weights chain_resampled = sample(chain, Weights(vec(chain[:weights])), length(chain)); -nothing # hide -``` -## Results +# ## Results -```@example eggbox chain_resampled -``` -```@example eggbox +# + marginalkde(chain[:x], chain[:y]) plot!(xlims=(0, 1), ylims=(0, 1), sp=2) plot!(xlims=(0, 1), sp=1) plot!(ylims=(0, 1), sp=3) -``` -```@example eggbox +# + density(chain_resampled, xlims=(0, 1)) vline!(0.1:0.2:0.9, c=:black, ls=:dash, sp=1) vline!(0.1:0.2:0.9, c=:black, ls=:dash, sp=2) -``` -```@example eggbox +# + using Measurements logz_est = state.logz ± state.logzerr diff = logz_est - logz println("logz: $logz") println("estimate: $logz_est") println("diff: $diff") -nothing # hide -``` diff --git a/docs/src/examples/shells.md b/examples/shells.jl similarity index 72% rename from docs/src/examples/shells.md rename to examples/shells.jl index ae432cb1..92286c85 100644 --- a/docs/src/examples/shells.md +++ b/examples/shells.jl @@ -1,3 +1,4 @@ +#= # Gaussian Shells This example will explore the classic Gaussian shells model using [`Models.GaussianShells`](@ref). @@ -8,26 +9,20 @@ For this example, you'll need to add the following packages ```julia julia>]add Distributions MCMCChains Measurements NestedSamplers StatsBase StatsPlots ``` - -```@setup shells +=# using AbstractMCMC using Random AbstractMCMC.setprogress!(false) Random.seed!(8452) -``` -## Define model +# ## Define model -```@example shells using NestedSamplers -model, logz = Models.GaussianShells() -nothing; # hide -``` +model, logz = Models.GaussianShells(); -let's take a look at a couple of parameters to see what the likelihood surface looks like +# let's take a look at a couple of parameters to see what the likelihood surface looks like -```@example shells using StatsPlots x = range(-6, 6, length=1000) @@ -41,47 +36,40 @@ heatmap( xlabel="x", ylabel="y", ) -``` -## Sample +# ## Sample -```@example shells using MCMCChains using StatsBase -# using multi-ellipsoid for bounds -# using default rejection sampler for proposals +## using multi-ellipsoid for bounds +## using default rejection sampler for proposals sampler = Nested(2, 1000) chain, state = sample(model, sampler; dlogz=0.05, param_names=["x", "y"]) -# resample chain using statistical weights +## resample chain using statistical weights chain_resampled = sample(chain, Weights(vec(chain[:weights])), length(chain)); -nothing # hide -``` -## Results +# ## Results -```@example shells chain_resampled -``` -```@example shells +# + marginalkde(chain[:x], chain[:y]) plot!(xlims=(-6, 6), ylims=(-2.5, 2.5), sp=2) plot!(xlims=(-6, 6), sp=1) plot!(ylims=(-2.5, 2.5), sp=3) -``` -```@example shells +# + density(chain_resampled) vline!([-5.5, -1.5, 1.5, 5.5], c=:black, ls=:dash, sp=1) vline!([-2, 2], c=:black, ls=:dash, sp=2) -``` -```@example shells +# + using Measurements logz_est = state.logz ± state.logzerr diff = logz_est - logz println("logz: $logz") println("estimate: $logz_est") println("diff: $diff") -nothing # hide -```