rewrite integrator

.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

bench/CondaPkg.toml

@@ -0,0 +1,2 @@
+[deps]
+dynesty = ""

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"

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"

bench/sampling.jl

@@ -2,27 +2,46 @@
 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, 4, 8, 16, 32]
-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(), 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 = @elapsed sample(rng, model, splr; dlogz=0.01, chain_type=Array)
+
+    dlnZ = state.logz - true_lnZ
 
-    tt = @belapsed sample($model, $splr; dlogz=0.01, chain_type=Array)
+    row = (; library="NestedSamplers.jl", D, t=t, lnZ=state.logz, lnZstd=state.logzerr, dlnZ)
+    @info "$row"
+    push!(rows, row)
+end
 
-    dlnZ = abs(true_lnZ - lnZ)
+@progress name="dynesty" for D in dims
+    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="NestedSamplers.jl", D, t=median(tt), lnZ, lnZstd, dlnZ)
+    row = (; library="dynesty", D, t, lnZ, lnZstd, dlnZ)
     @info "$row"
     push!(rows, row)
 end

bench/sampling_results.csv

@@ -1,6 +1,9 @@
 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.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

docs/Project.toml

@@ -4,10 +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"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+NestedSamplers = "41ceaf6f-1696-4a54-9b49-2e7a9ec3782e"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"

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")

docs/src/benchmarks.md

@@ -65,19 +65,25 @@ 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
 
 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
 ```

docs/src/examples/correlated.md → 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
-```

docs/src/examples/eggbox.md → 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
-```