Add test gaussian

.gitignore

@@ -10,3 +10,4 @@ coverage
 docs/build/
 env
 node_modules
+example/*.pdf

Project.toml

@@ -1,9 +1,6 @@
 name = "InclusiveScans"
 uuid = "ef74772d-374f-5edb-a43a-0dbb31824c9d"
-authors = [
-    "Uwe Hernandez Acosta",
-    " Simeon Ehrig",
-]
+authors = ["Uwe Hernandez Acosta", "Simeon Ehrig"]
 version = "0.1.0"
 
 [compat]

examples/ConstrainedGaussians/Project.toml

@@ -0,0 +1,18 @@
+name = "ConstrainedGaussians"
+uuid = "72cb5a96-92ad-43ae-9774-2d852a9e5056"
+authors = ["Uwe Hernandez Acosta <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
+
+[compat]
+Distributions = "0.25.116"
+DomainSets = "0.7.14"
+QuadGK = "2.11.1"
+Random = "1.11.0"
+StatsPlots = "0.15.7"

examples/ConstrainedGaussians/src/ConstrainedGaussians.jl

@@ -0,0 +1,22 @@
+module ConstrainedGaussians
+
+
+using Distributions
+using DomainSets
+using Random
+using QuadGK
+using StatsPlots
+#using InclusiveScans
+
+export ConstrainedGaussian1D
+export target_dist, filter, norm, normalize
+export endpoints
+export generate_CPU
+export plot_compare
+
+include("interface.jl")
+include("impl.jl")
+include("generate.jl")
+include("plot_compare.jl")
+
+end

examples/ConstrainedGaussians/src/generate.jl

@@ -0,0 +1,36 @@
+
+function generate_batch_CPU(rng::AbstractRNG, batch_size::Int, dist::ConstrainedGaussian1D)
+    rand_args = rand(rng, Uniform(DomainSets.endpoints(dist.domain)...), batch_size)
+
+    # opt this block out to be evaluated on a single element
+    # and broadcast over everything
+    # -------
+    rand_vals = target_dist.(dist, rand_args)
+    rel_rand_vals = normalize(dist, rand_vals)
+    rand_probs = rand(rng, batch_size)
+    mask = filter.(rel_rand_vals, rand_probs)
+    # -------
+
+    accepted_args = rand_args[mask] # this one with InclusiveScan.jl
+
+    return accepted_args
+end
+
+function generate_CPU(
+    rng::AbstractRNG,
+    dist::ConstrainedGaussian1D{T},
+    N::Int;
+    batch_size::Int = 1000,
+) where {T}
+    accepted_args = T[]
+    sizehint!(accepted_args, N + batch_size - 1)
+
+    nrun = 0
+    while nrun <= N
+        accepted_args_batch = generate_batch_CPU(rng, batch_size, dist)
+        append!(accepted_args, accepted_args_batch)
+        nrun += length(accepted_args_batch)
+    end
+
+    return accepted_args
+end

examples/ConstrainedGaussians/src/impl.jl

@@ -0,0 +1,39 @@
+
+# constrained Gaussian
+
+struct ConstrainedGaussian1D{T<:Real,DIST,DOM<:AbstractInterval} <: AbstractTestDist{T}
+    mu::T
+    sig::T
+    dist::DIST
+    domain::DOM
+
+    function ConstrainedGaussian1D(
+        mu::T,
+        sig::T,
+        domain::DOM = Interval(-5 * sig, 5 * sig),
+    ) where {T<:Real,DOM<:AbstractInterval}
+        dist = Normal(mu, sig)
+        return new{T,typeof(dist),DOM}(mu, sig, dist, domain)
+    end
+end
+
+function ConstrainedGaussian1D(mu::T, sig::T, domain::D) where {T<:Real,D<:Tuple}
+    return ConstrainedGaussian1D(mu, sig, Interval(domain...))
+end
+
+endpoints(dist::ConstrainedGaussian1D) = DomainSets.endpoints(dist.domain)
+
+function norm(d::ConstrainedGaussian1D{T}) where {T}
+    mu = d.mu
+    dom = d.domain
+
+    if mu in dom
+        return target_dist(d, mu)
+    end
+
+    if mu <= minimum(dom)
+        return target_dist(d, minimum(dom))
+    end
+
+    return target_dist(d, maximum(dom))
+end

examples/ConstrainedGaussians/src/interface.jl

@@ -0,0 +1,16 @@
+
+
+# general abstract test-distribution
+
+abstract type AbstractTestDist{T} end
+
+Base.broadcastable(dist::AbstractTestDist) = Ref(dist)
+target_dist(d::AbstractTestDist{T}, x) where {T} = x in d.domain ? pdf(d.dist, x) : zero(T)
+
+function normalize(d, x)
+    return x / norm(d)
+end
+
+function filter(rel_rand_val, rand_prob)
+    return rel_rand_val >= rand_prob
+end

examples/ConstrainedGaussians/src/plot_compare.jl

@@ -0,0 +1,24 @@
+function plot_compare(samples, dist; kwargs...)
+    P = histogram(
+        samples;
+        label = "samples",
+        xlabel = "x",
+        ylabel = "normalized event count",
+        nbins = 100,
+        normalize = :pdf,
+        opacity = 0.5,
+        kwargs...,
+    )
+
+    tot_weight, _ = quadgk(x -> target_dist(dist, x), DomainSets.endpoints(dist.domain)...)
+    plot!(
+        P,
+        range(endpoints(dist)...; length = 100),
+        x -> target_dist(dist, x) / tot_weight;
+        label = "normalized target dist.",
+        line = (2, :black, :dash),
+        alpha = 0.5,
+    )
+
+    return P
+end

examples/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+ConstrainedGaussians = "72cb5a96-92ad-43ae-9774-2d852a9e5056"
+InclusiveScans = "ef74772d-374f-5edb-a43a-0dbb31824c9d"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

examples/README.md

@@ -0,0 +1,91 @@
+# Example: Constrained Gaussian
+
+This repository demonstrates how to generate samples from a **constrained Gaussian distribution**,
+a non-trivial extension of the Gaussian distribution confined to a specified interval.
+The distribution and sampling functionality (CPU-based) are implemented in the `ConstrainedGaussians.jl` package.
+
+## Features
+
+- sampling of constrained Gaussian distributions.
+- Comparison of generated samples with the target distribution using histograms.
+
+---
+
+## Setup Instructions
+
+### 1. Clone the Repository
+
+Ensure you have cloned the repository containing this example.
+
+```bash
+git clone https://github.com/QEDjl-project/InclusiveScans.jl
+cd <repository-directory>/example
+```
+
+### 2. Initialize the Example Environment
+
+Since `InclusiveScans.jl` and `ConstrainedGaussians.jl` are not registered Julia packages,
+you need to initialize the environment manually. Run the following command inside the `example` directory:
+
+```bash
+julia --project=. init.jl
+```
+
+This script will:
+
+- Link the source code of `InclusiveScans.jl` and `ConstrainedGaussians.jl` to the example environment.
+- Install all other necessary dependencies.
+- Instantiate the example environment
+
+---
+
+## Running the Example
+
+### 1. Execute the Example Script
+
+To generate samples and visualize the constrained Gaussian distribution, execute the following command inside the `example` directory:
+
+```bash
+julia --project example.jl
+```
+
+### 2. Output
+
+- The script generates **1 million samples** from a constrained Gaussian distribution.
+- A histogram comparing the generated samples with the target distribution is plotted.
+- The histogram is saved as `example_plot_compare.pdf` in the `example` directory.
+
+### 3. Interactive Exploration (Optional)
+
+For interactive experimentation, open the Jupyter notebook `example.iypnb`:
+
+1. Ensure you have the `IJulia` package installed:
+
+   ```julia
+   using Pkg
+   Pkg.add("IJulia")
+   ```
+
+2. Start Jupyter and open the notebook:
+
+   ```bash
+   jupyter notebook
+   ```
+
+This allows you to explore the constrained Gaussian example interactively using Julia's kernel.
+
+---
+
+## Requirements
+
+Ensure you have the following tools installed:
+
+- **Julia**: Version 1.6 or later is recommended.
+- **Jupyter** (optional): For interactive notebooks.
+
+---
+
+## Troubleshooting
+
+- **Dependency Issues**: If you encounter problems during initialization, verify that the `init.jl` script executed successfully and installed all required packages.
+- **Missing `IJulia`**: Ensure the `IJulia` package is installed to enable Julia kernels in Jupyter.