Merge pull request #26 from Eliassj/master

Merge ARSampling.jl (JuliaRegistries/General/151430)

Author: Eliassj

.github/workflows/Docs.yml

@@ -0,0 +1,62 @@
+name: Docs
+
+on:
+  push:
+    branches:
+      - master
+    paths:
+      - "docs/**"
+      - "src/**"
+      - "*.toml"
+    tags: ["*"]
+  pull_request:
+    branches:
+      - master
+    paths:
+      - "docs/**"
+      - "src/**"
+      - "*.toml"
+    types: [opened, synchronize, reopened]
+
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+
+jobs:
+  docs:
+    permissions:
+      actions: write
+      contents: write
+      pull-requests: read
+      statuses: write
+    name: Documentation
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: "1"
+      - name: Use Julia cache
+        uses: julia-actions/cache@v2
+      - name: Instantiate environment with development version of the package
+        run: |
+          julia --project=docs -e '
+            using Pkg
+            Pkg.develop(PackageSpec(path=pwd()))
+            Pkg.instantiate()'
+      - name: Run doctest
+        run: |
+          julia --project=docs -e '
+            using Documenter: DocMeta, doctest
+            using ARSampling
+            DocMeta.setdocmeta!(ARSampling, :DocTestSetup, :(using ARSampling); recursive=true)
+            doctest(ARSampling)'
+      - name: Generate and deploy documentation
+        run: julia --project=docs docs/make.jl
+        env:
+          JULIA_PKG_SERVER: ""
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
+          GKSwstype: "100" # https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988

.github/workflows/ReusableTest.yml

@@ -0,0 +1,52 @@
+name: Reusable test
+
+on:
+  workflow_call:
+    inputs:
+      version:
+        required: false
+        type: string
+        default: "1"
+      os:
+        required: false
+        type: string
+        default: ubuntu-latest
+      arch:
+        required: false
+        type: string
+        default: x64
+      allow_failure:
+        required: false
+        type: boolean
+        default: false
+      run_codecov:
+        required: false
+        type: boolean
+        default: false
+    secrets:
+      codecov_token:
+        required: true
+
+jobs:
+  test:
+    name: Julia ${{ inputs.version }} - ${{ inputs.os }} - ${{ inputs.arch }} - ${{ github.event_name }}
+    runs-on: ${{ inputs.os }}
+    continue-on-error: ${{ inputs.allow_failure }}
+
+    steps:
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ inputs.version }}
+          arch: ${{ inputs.arch }}
+      - name: Use Julia cache
+        uses: julia-actions/cache@v2
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+        if: ${{ inputs.run_codecov }}
+      - uses: codecov/codecov-action@v4
+        if: ${{ inputs.run_codecov }}
+        with:
+          file: lcov.info
+          token: ${{ secrets.codecov_token }}

.github/workflows/Test.yml

@@ -0,0 +1,35 @@
+name: Test
+
+on:
+  push:
+    branches:
+      - master
+    tags: ["*"]
+  workflow_dispatch:
+
+jobs:
+  test:
+    uses: ./.github/workflows/ReusableTest.yml
+    with:
+      os: ${{ matrix.os }}
+      version: ${{ matrix.version }}
+      arch: ${{ matrix.arch }}
+      allow_failure: ${{ matrix.allow_failure }}
+      run_codecov: ${{ matrix.version == '1' && matrix.os == 'ubuntu-latest' }}
+    secrets:
+      codecov_token: ${{ secrets.CODECOV_TOKEN }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - "lts"
+          - "1"
+        os:
+          - ubuntu-latest
+          - macOS-latest
+
+          - windows-latest
+
+        arch:
+          - x64
+        allow_failure: [false]

Project.toml

@@ -3,14 +3,27 @@ uuid = "c75e803d-635f-53bd-ab7d-544e482d8c75"
 version = "0.2.1"
 
 [deps]
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
+[weakdeps]
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+
 [compat]
-ForwardDiff = "0.10.1"
-StatsBase = "0.24, 0.34"
-julia = "0.7.0, 1"
+Compat = "4.18.1"
+DifferentiationInterface = "0.7"
+ForwardDiff = "1"
+Makie = "0.24.9"
+SpecialFunctions = "2"
+StatsBase = "0.34.10"
+julia = "1"
+
+[extensions]
+MakieExtension = "Makie"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -12,43 +12,62 @@ This package is useful for efficient sampling from log-concave univariate densit
 ## Examples
 
 ```julia
-using AdaptiveRejectionSampling
-using Plots
+using AdaptiveRejectionSampling: Objective, ARSampler, sample!, eval_hull, abscissae
+using SpecialFunctions: gamma
+using CairoMakie
 ```
 
-
 ### Sampling from a shifted normal distribution
 
 
 ```julia
-# Define function to be sampled
-μ, σ = 1.0, 2.0
-f(x) = exp(-0.5(x - μ)^2 / σ^2) / sqrt(2pi * σ^2) 
-support = (-Inf, Inf)
+const μ, σ = 1.0, 2.0
 
-# Build the sampler and simulate 10,000 samples
-sampler = RejectionSampler(f, support, max_segments = 5)
-@time sim = run_sampler!(sampler, 10000);
+function bench()
+    # Define function to be sampled
+    f(x) = log(exp(-0.5(x - μ)^2 / σ^2) / sqrt(2pi * σ^2)) 
+    support = (-10., 10.)
+    
+    # Build the sampler and simulate 10,000 samples
+    obj = Objective(f)
+    sampler = ARSampler(obj, support)
+    b = @be deepcopy(sampler) sample!(_, 10000, true, 10);
+return sampler, b
+end
+s, b = bench();
+b
 ```
 
-      0.010434 seconds (192.15 k allocations: 3.173 MiB)
-    
+      Benchmark: 126 samples with 1 evaluation
+      min    616.951 μs (17 allocs: 80.523 KiB)
+      median 641.532 μs (17 allocs: 80.523 KiB)
+      mean   749.812 μs (17.10 allocs: 80.625 KiB, 0.75% gc time)
+      max    13.744 ms (20 allocs: 83.711 KiB, 94.06% gc time)
 
 Let's verify the result
 
 
 ```julia
 # Plot the results and compare to target distribution
+f(x) = log(exp(-0.5(x - μ)^2 / σ^2) / sqrt(2pi * σ^2)) 
 x = range(-10.0, 10.0, length=100)
-envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
-target = [f(xi) for xi in x]
-
-histogram(sim, normalize = true, label = "Histogram")
-plot!(x, [target envelop], width = 2, label = ["Normal(μ, σ)" "Envelop"])
+envelop = eval_hull.(s.upper_hull, x)
+target = f.(x)
+samples = sample!(s, 10000, max_segments = 5)
+
+fig, ax, p = hist(samples, bins=25, normalization = :pdf, label = "Samples");
+lines!(ax, -10..10, x -> exp(f(x)), label = "Target", color = :orange)
+lines!(ax, x, exp.(envelop), label = "Envelop", color = :black, alpha = 0.8)
+absc = ARS.abscissae(s.upper_hull)
+scatter!(ax, absc, exp.(eval_hull.(s.upper_hull, absc)), label = "Abscissae", color = :red)
+axislegend(ax)
+fig
+# histogram(sim, normalize = true, label = "Histogram")
+# plot!(x, [target envelop], width = 2, label = ["Normal(μ, σ)" "Envelop"])
 ```
 
 
-![](img/example1.png)
+![](img/example1.svg)
 
 
 ### Let's try a Gamma
@@ -59,24 +78,28 @@ plot!(x, [target envelop], width = 2, label = ["Normal(μ, σ)" "Envelop"])
 f(x) = β^α * x^(α-1) * exp(-β*x) / gamma(α)
 support = (0.0, Inf)
 
-# Build the sampler and simulate 10,000 samples
-sampler = RejectionSampler(f, support)
-@time sim = run_sampler!(sampler, 10000) 
-
-# Plot the results and compare to target distribution
-x = range(0.0, 5.0, length=100)
-envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
-target = [f(xi) for xi in x]
-
-histogram(sim, normalize = true, label = "Histogram")
-plot!(x, [target envelop], width = 2, label = ["Gamma(α, β)" "Envelop"])
+obj = Objective(x -> log(f(x)))
+sam = ARSampler(obj, support)
+@time sim = sample!(sam, 10000, max_segments = 5)
+
+# Verify result
+x = range(0.0, 8.0, length=100)
+envelop = eval_hull.(sam.upper_hull, x)
+target = f.(x)
+mx = maximum(sim)
+fig, ax, p = hist(sim, bins=50, normalization = :pdf, label = "Samples");
+lines!(ax, 0..mx, y -> f(y), label = "Target", color = :orange)
+lines!(ax, 0..mx, y -> exp(eval_hull(sam.upper_hull, y)), label = "Envelop", color = :black, alpha = 0.8)
+absc = abscissae(sam.upper_hull)
+scatter!(ax, absc, exp.(eval_hull.(sam.upper_hull, absc)), label = "Abscissae", color = :red)
+axislegend(ax)
+fig
 ```
 
-      0.007299 seconds (182.00 k allocations: 3.027 MiB)
-    
+      0.000772 seconds (100 allocations: 85.211 KiB)
 
 
-![](img/example2.png)
+![](img/example2.svg)
 
 ### Truncated distributions and unknown normalization constant
 

ext/MakieExtension.jl

@@ -0,0 +1,81 @@
+module MakieExtension
+using Makie
+import Makie: mixin_colormap_attributes
+using AdaptiveRejectionSampling
+import AdaptiveRejectionSampling: hullplot, hullplot!, eval_hull, abscissae, ARSampler
+
+@recipe Hullplot begin
+    abscissae = true
+    target = true
+    upper_hull = true
+    lower_hull = true
+    target_linewidth = @inherit linewidth
+    upper_hull_linewidth = @inherit linewidth
+    lower_hull_linewidth = @inherit linewidth
+    target_label = "Target"
+    upper_hull_label = "Upper hull"
+    lower_hull_label = "Lower hull"
+    abscissae_label = "Abscissae"
+end
+
+
+function Makie.plot!(p::Hullplot{<:Tuple{<:Any, <:ARSampler}})
+
+    sam = p[2][]
+    upper_hull = sam.upper_hull
+    lower_hull = sam.lower_hull
+    obj = sam.objective
+    xs = p[1]
+
+    i = 1
+    if p.lower_hull[]
+        lines!(
+            p,
+            xs,
+            x -> exp(eval_hull(lower_hull, x)),
+            color=Cycled(i),
+            linewidth = p.lower_hull_linewidth,
+            label = p.lower_hull_label
+        )
+        i += 1
+    end
+    if p.upper_hull[]
+        lines!(
+            p,
+            xs,
+            x -> exp(eval_hull(upper_hull, x)),
+            linewidth = p.upper_hull_linewidth,
+            color=Cycled(i),
+            label = p.upper_hull_label
+        )
+        i += 1
+    end
+    if p.target[]
+        lines!(
+            p,
+            xs,
+            x -> exp(obj.f(x)),
+            color=Cycled(i),
+            label = p.target_label,
+            linewidth = p.target_linewidth
+        )
+        i += 1
+    end
+    if p.abscissae[]
+        absc = abscissae(upper_hull)
+        scatter!(
+            p,
+            absc,
+            exp.(eval_hull.(upper_hull, absc)),
+            color = Cycled(i),
+            label = p.abscissae_label
+        )
+    end
+
+end
+
+function Makie.get_plots(plot::Hullplot)
+    return plot.plots
+end
+
+end # End of module