### Install Documentation Dependencies

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/README.md

Install the necessary Julia packages for building and serving the documentation locally.

```julia
julia -e 'using Pkg; Pkg.add(["Documenter", "DocumenterVitepress", "LiveServer"])'
```

--------------------------------

### Install NeuralEstimators.jl (Development)

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/README.md

Install the current development version of the NeuralEstimators.jl package from GitHub.

```julia
using Pkg; Pkg.add(url = "https://github.com/msainsburydale/NeuralEstimators.jl")
```

--------------------------------

### Activate and Instantiate Package Environment

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/CONTRIBUTING.md

Activates the package environment from the repository root and installs dependencies. This is a crucial first step for development.

```julia
using Pkg; Pkg.activate("."); Pkg.instantiate()
```

--------------------------------

### Activate and Instantiate Package Environment (Command Line)

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/CONTRIBUTING.md

Activates the package environment and instantiates dependencies using a Julia command-line execution. Useful for scripting or automated setups.

```bash
julia --project=. -e "using Pkg; Pkg.instantiate()"
```

--------------------------------

### Load GPU Acceleration Packages

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/advancedusage.md

Load the appropriate package for your GPU. Ensure the package is installed and a compatible GPU is available.

```julia
using CUDA, cuDNN
```

```julia
using AMDGPU
```

```julia
using Metal
```

```julia
using oneAPI
```

--------------------------------

### Install NeuralEstimators.jl (Stable)

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/README.md

Install the current stable release of the NeuralEstimators.jl package using Julia's package manager.

```julia
using Pkg; Pkg.add("NeuralEstimators")
```

--------------------------------

### Setup and Model Definition for Gaussian Process with Missing Data

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_missing_censored.md

This snippet sets up the necessary Julia packages and defines the prior distribution, spatial domain, and distance matrix for a Gaussian process model. It also defines structures for parameters and a Matern covariance function.

```julia
using NeuralEstimators, Flux
using Distributions: Uniform
using Distances, LinearAlgebra
using MLUtils: flatten
using SpecialFunctions: besselk, gamma
using Statistics: mean

# Prior and dimension of parameter vector
Π = (τ = Uniform(0, 1.0), ρ = Uniform(0, 0.4))
d = length(Π)

# Define the grid and compute the distance matrix
points = range(0, 1, 16)
S = expandgrid(points, points)
D = pairwise(Euclidean(), S, dims = 1)

# Collect model information for later use
ξ = (Π = Π, S = S, D = D)

# Struct for storing parameters and Cholesky factors
struct Parameters <: AbstractParameterSet
	θ
	L
end

# Matern covariance function
function matern(h, ρ, ν, σ² = one(typeof(h)))

    @assert h >= 0 "h should be non-negative"
    @assert ρ > 0 "ρ should be positive"
    @assert ν > 0 "ν should be positive"

    if h == 0
        σ²
    else
        d = h / ρ
        σ² * ((2^(1 - ν)) / gamma(ν)) * d^ν * besselk(ν, d)
    end
end

```

--------------------------------

### Ratio Estimation with NeuralEstimators.jl

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/index.md

Develops a neural estimator for parameters θ from data Z, using a specified prior and simulator, optimized for ratio estimation. This example demonstrates setting up, training, and assessing a RatioEstimator, including device specification.

```julia
using NeuralEstimators, Lux, Reactant

# Functions to sample from the prior and simulate data
d, n = 2, 100  # dimension of θ and number of replicates
sampler(K) = NamedMatrix(μ = randn(K), σ = rand(K))
simulator(θ::AbstractVector) = θ["μ"] .+ θ["σ"] .* sort(randn(n))
simulator(θ::AbstractMatrix) = reduce(hcat, map(simulator, eachcol(θ)))

# Neural network mapping n inputs into d outputs
network = Chain(Dense(n, 64, gelu), Dense(64, 64, gelu), Dense(64, d))

# Initialise a neural estimator
estimator = RatioEstimator(network, d; num_summaries = d)

# Train the estimator
estimator = train(estimator, sampler, simulator, device = reactant_device())

# Assess the estimator
θ_test = sampler(250)
Z_test = simulator(θ_test)
assessment = assess(estimator, θ_test, Z_test)
bias(assessment)
rmse(assessment)

# Apply to observed data (here, simulated as a stand-in)
θ = sampler(1)                   # ground truth (not known in practice)
Z = simulator(θ)                 # stand-in for real observations
sampleposterior(estimator, Z)    # approximate posterior sample
```

--------------------------------

### Install NeuralEstimators.jl and Dependencies

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

Installs essential Julia packages for neural inference, including Flux for neural networks, BSON for serialization, plotting libraries, and the NeuralEstimators.jl package itself from a specific Git revision.

```julia
# Install Julia packages
using Pkg
Pkg.add("Flux")
Pkg.add("BSON")
Pkg.add(["Plots", "StatsPlots", "CairoMakie"]) # Visualization 
Pkg.add(["Downloads", "RData"]) # Loading presimulated draws from 𝑝(𝜃, 𝐙)

Pkg.add(url = "https://github.com/msainsburydale/NeuralEstimators.jl", rev="workshop")

if install_cuda
  Pkg.add("CUDA")
end
```

--------------------------------

### Package dependencies for Neural Estimators

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Import necessary packages for using NeuralEstimators and statistical functions. Ensure these are installed before running the code.

```julia
using NeuralEstimators
using Statistics: mean, std
```

--------------------------------

### Point Estimation with NeuralEstimators.jl

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/index.md

Develops a neural estimator for parameters θ from data Z, using a specified prior and simulator. This example demonstrates setting up, training, and assessing a PointEstimator.

```julia
using NeuralEstimators, Lux, Enzyme

# Functions to sample from the prior and simulate data
d, n = 2, 100  # dimension of θ and number of replicates
sampler(K) = NamedMatrix(μ = randn(K), σ = rand(K))
simulator(θ::AbstractVector) = θ["μ"] .+ θ["σ"] .* sort(randn(n))
simulator(θ::AbstractMatrix) = reduce(hcat, map(simulator, eachcol(θ)))

# Neural network mapping n inputs into d outputs
network = Chain(Dense(n, 64, gelu), Dense(64, 64, gelu), Dense(64, d))

# Initialise a neural estimator
estimator = PointEstimator(network, d; num_summaries = d)

# Train the estimator
estimator = train(estimator, sampler, simulator)

# Assess the estimator
θ_test = sampler(250)
Z_test = simulator(θ_test)
assessment = assess(estimator, θ_test, Z_test)
bias(assessment)
rmse(assessment)

# Apply to observed data (here, simulated as a stand-in)
θ = sampler(1)                   # ground truth (not known in practice)
Z = simulator(θ)                 # stand-in for real observations
estimate(estimator, Z)           # point estimate
```

--------------------------------

### Posterior Estimation with NeuralEstimators.jl

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/index.md

Develops a neural estimator for parameters θ from data Z, using a specified prior and simulator, with a Gaussian Mixture approximation for the posterior. This example demonstrates setting up, training, and assessing a PosteriorEstimator.

```julia
using NeuralEstimators, Lux, Enzyme

# Functions to sample from the prior and simulate data
d, n = 2, 100  # dimension of θ and number of replicates
sampler(K) = NamedMatrix(μ = randn(K), σ = rand(K))
simulator(θ::AbstractVector) = θ["μ"] .+ θ["σ"] .* sort(randn(n))
simulator(θ::AbstractMatrix) = reduce(hcat, map(simulator, eachcol(θ)))

# Neural network mapping n inputs into d outputs
network = Chain(Dense(n, 64, gelu), Dense(64, 64, gelu), Dense(64, d))

# Initialise a neural estimator
estimator = PosteriorEstimator(network, d; num_summaries = d, q = GaussianMixture)

# Train the estimator
estimator = train(estimator, sampler, simulator)

# Assess the estimator
θ_test = sampler(250)
Z_test = simulator(θ_test)
assessment = assess(estimator, θ_test, Z_test)
bias(assessment)
rmse(assessment)

# Apply to observed data (here, simulated as a stand-in)
θ = sampler(1)                   # ground truth (not known in practice)
Z = simulator(θ)                 # stand-in for real observations
sampleposterior(estimator, Z)    # approximate posterior sample
```

--------------------------------

### Quick Start: Neural Bayes Estimation (NBE)

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/README.md

Construct, train, and apply a Neural Bayes Estimator for parameters (μ, σ) from simulated normal data. Requires NeuralEstimators and Flux.jl. The simulator function generates data based on provided parameters.

```julia
using NeuralEstimators
using Flux

# Dimension of θ and number of replicates
d, n = 2, 100  

# Functions to sample from prior p(θ) and simulate data p(Z|θ)
sampler(K) = NamedMatrix(μ = randn(K), σ = rand(K))
simulator(θ::AbstractVector) = θ["μ"] .+ θ["σ"] .* sort(randn(n))
simulator(θ::AbstractMatrix) = reduce(hcat, map(simulator, eachcol(θ)))

# Neural network, an MLP mapping n inputs into d outputs
network = Chain(Dense(n, 64, gelu), Dense(64, 64, gelu), Dense(64, d))

# Initialise an NBE
estimator = PointEstimator(network)

# Train the estimator
estimator = train(estimator, sampler, simulator)

# Apply to observed data
Z = simulator(sampler(1))  # stand-in for real observations
estimate(estimator, Z)     # point estimate
```

--------------------------------

### Load Package Dependencies

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_spatiotemporal.md

Loads necessary packages for spatiotemporal data analysis and neural network operations. Ensure Julia is started with multi-threading support for parallel simulations.

```julia
using NeuralEstimators
using NeuralEstimators: getobs, numobs
using Flux
using Folds                    # parallel simulation (start Julia with --threads=auto)
using Distances
using Distributions: Uniform
using LinearAlgebra
using Statistics
using CairoMakie
```

--------------------------------

### Build and Serve Documentation Locally

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/README.md

Navigate to the docs directory and run this command to build the documentation and serve it locally for preview.

```bash
julia --project=. make.jl && julia -e 'using LiveServer; serve(dir="build/1")'
```

--------------------------------

### Constructing a Lux-based summary network

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Initialize a summary network using Lux, setting it to the identity function as expert summaries are precomputed.

```julia
summary_network = Lux.WrappedFunction(identity)
```

--------------------------------

### Run Tests with Project and Pkg

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/test/README.md

Run tests for the current project using the julia command with the --project flag and Pkg.test.

```bash
julia --project=. -e "using Pkg; Pkg.test()"
```

--------------------------------

### Apply Estimator to Observed Data (Point)

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Applies a trained point estimator to observed data to get a point estimate.

```julia
θ = sampler(1)                 # ground truth (not known in practice)
Z = simulator(θ)               # stand-in for real observations

estimate(estimator, Z)             # point estimate
```

--------------------------------

### Deep-Learning Backend: Lux with Reactant

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_temporal.md

Import the Lux deep-learning library and configure Reactant for GPU usage. For optimal performance, pass `device = reactant_device()` to the `train()` function.

```julia
using Lux

# NB: for the most computationally efficient setup, pass device = reactant_device() to train()
using Reactant
Reactant.set_default_backend("gpu")
```

--------------------------------

### Set GPU Acceleration Flag

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

This code snippet sets a boolean flag to control GPU acceleration. It is set to false for the workshop to avoid lengthy installation times.

```python
# Set to true to enable GPU acceleration (requires an NVIDIA GPU and CUDA;
# we use the CPU in this workshop due to Colab installation time)
install_cuda = false
```

--------------------------------

### Intel GPU Support

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_replicated.md

Load oneAPI package for Intel GPU acceleration.

```julia
using oneAPI
```

--------------------------------

### Prepare Data for Estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_gridded.md

Prepare ground truth parameters and simulate data to be used as a stand-in for observed data.

```julia
θ = Parameters(Matrix([0.1]"))   # ground truth (not known in practice)
Z = simulator(θ)                 # stand-in for real data
```

--------------------------------

### Apply Estimator to Observed Data

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_replicated.md

Apply a trained estimator to observed data to obtain point estimates or posterior samples. This example uses simulated data as a stand-in for real observed data.

```julia
θ = sampler(1)                      # ground truth (not known in practice)
Z = simulator(θ, 100)               # stand-in for real data
```

--------------------------------

### Constructing a Flux-based summary network

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Initialize a summary network using Flux, setting it to the identity function as expert summaries are precomputed.

```julia
summary_network = identity
```

--------------------------------

### SimpleChains.jl Integration via Lux.jl

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/advancedusage.md

Integrate NeuralEstimators.jl with SimpleChains.jl through Lux.jl for CPU-optimized small networks. Requires defining a Lux network and converting it.

```julia
using NeuralEstimators, Lux, SimpleChains

# Define Lux network and convert it to SimpleChains
network  = Lux.Chain(...)
adaptor  = ToSimpleChainsAdaptor(...) # declare input size
network  = adaptor(network)

# Then proceed with Lux workflow... 
```

--------------------------------

### Run Tests with Pkg.test

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/test/README.md

Execute the test suite for the NeuralEstimators package using the Pkg.test function.

```julia
Pkg.test("NeuralEstimators")
```

--------------------------------

### Visualize Model Realizations (Julia)

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

Generates and plots heatmap visualizations for model realizations across a range of parameter values to illustrate phase transitions. This code requires the Plots package.

```julia
# True parameter values
theta = [0.2 0.45 0.7 0.881 1.0 1.2 1.4]

# Simulations
Z = simulate(theta)

# Create plots for each theta value
plots_array = []
for k in 1:size(theta, 2)
    z = Z[:, :, 1, k]
    p = Plots.heatmap(z', c=:magma,
                title="θ = $(theta[k])",
                colorbar=false,
                aspect_ratio=:equal,
                axis=nothing,
                border=:none,
                framestyle=:none)
    push!(plots_array, p)
end

# Combine into single plot
Plots.plot(plots_array..., layout=(1, 7), size=(2000, 400))
```

--------------------------------

### Create Ratio Estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Initializes a RatioEstimator with a specified summary network and number of summaries.

```julia
estimator = RatioEstimator(summary_network, d; num_summaries = num_expert_summaries + num_learned_summaries)
```

--------------------------------

### Sampler Constructor for Training

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_irregularspatial.md

Define a constructor `sampler` that accepts an integer `K` and samples parameters `θ`, spatial configurations `S` from a Matérn cluster process, and returns a `Parameters` object. This is used during training.

```julia
function sampler(K::Integer)
	θ = 0.5 * rand(1, K)

	# Sample spatial configurations from a Matérn cluster process on [0, 1]²
	n = rand(200:300, K)
	λ = rand(Uniform(10, 50), K)
	S = [maternclusterprocess(λ = λ[k], μ = n[k]/λ[k]) for k ∈ 1:K]

	Parameters(θ, S)
end
```

--------------------------------

### Load Julia Packages

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

This snippet loads core packages for Flux, BSON, NeuralEstimators, plotting with StatsPlots and CairoMakie, and data handling with Downloads and RData. It also includes conditional loading for CUDA if a GPU is to be used.

```julia
# Load Julia packages
using Flux
using BSON: @load
using NeuralEstimators
using Plots, StatsPlots, CairoMakie
using CairoMakie: plot
using Downloads, RData
using Downloads: download
using Statistics

# If using a GPU, load CUDA and check that the GPU is functional
if install_cuda
  using CUDA
  CUDA.functional()
end
```

--------------------------------

### Train Neural Estimator with On-the-fly Simulation

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_spatiotemporal.md

Alternatively, train the estimator by simulating data on-the-fly during the training process. This is useful when pre-simulation is not feasible or desired.

```julia
# Alternatively, simulate on-the-fly:
# estimator = train(estimator, sampler, simulator; K = K)
```

--------------------------------

### Download and Load Pre-trained Model

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

Download a pre-trained neural network model from a URL and load its state into the current model. This is useful for reusing models trained elsewhere, especially on GPUs.

```python
# Specify which estimator we are downloading
estimator_file = "NBE.bson"

# Download the file
download(workshop_url * estimator_file, estimator_file)

# Load the trained neural network
model_state = Flux.state(NBE)
@load estimator_file model_state
NBE = Flux.loadmodel!(NBE, model_state)
```

--------------------------------

### Define Network Parameters

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_temporal.md

Sets up parameters for neural network construction, including observation dimension and number of parameters.

```julia
p = 2    # dimension of each observation
d = 6    # number of parameters (entries of A and innovation standard deviations)
num_summaries = 3d
```

--------------------------------

### Define sampler and simulator functions

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Define functions to sample parameters from a prior distribution and simulate data. The simulator returns summary statistics directly.

```julia
d = 2              # dimension of θ
num_summaries = 3  # number of expert summaries

# Function to sample from the prior
sampler(K) = NamedMatrix(μ = randn(K), σ = rand(K))

# Function to simulate data
function simulator(θ::AbstractVector, n::Integer)
    Z = θ["μ"] .+ θ["σ"] .* randn(n)
    S = [mean(Z), std(Z), log(n)]
    return S
end
simulator(θ::AbstractVector, n) = simulator(θ, rand(n))
simulator(θ::AbstractMatrix, n) = reduce(hcat, simulator.(eachcol(θ), Ref(n)))
```

--------------------------------

### Sample from Prior Distribution (Julia)

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

Defines a function to sample parameters from a uniform prior distribution. The function takes the number of samples as input and returns parameters as a matrix.

```julia
# Sampling from the prior distribution
# N: number of samples to draw from the prior
function prior(N)
    theta = rand(N) .* (1.5 - 0.05) .+ 0.05  # uniform(0.05, 1.5)
    return theta'  # return as matrix
end
```

--------------------------------

### Sample VAR Process Parameters

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_temporal.md

Samples parameters for a VAR process, ensuring stationarity by rejecting draws with a spectral radius >= 1. The output is a NamedMatrix containing coefficients and innovation standard deviations.

```julia
function sampler(K)
    p = 2   # dimension of the VAR process
    A = Matrix{Float64}(undef, p^2, 0)
    while size(A, 2) < K
        candidates = 0.5 .* randn(p^2, 2K)   # oversample then filter
        valid = [ρ(reshape(candidates[:, k], p, p)) < 1 for k in 1:2K]
        A = hcat(A, candidates[:, valid])
    end
    NamedMatrix(
        A₁₁ = A[1, 1:K],
        A₁₂ = A[2, 1:K],
        A₂₁ = A[3, 1:K],
        A₂₂ = A[4, 1:K],
        σ₁  = abs.(0.5 .* randn(K)),
        σ₂  = abs.(0.5 .* randn(K)),
    )
end

ρ(A) = maximum(abs.(eigvals(A)))   # spectral radius
```

--------------------------------

### Package Dependencies

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_replicated.md

Load necessary packages for neural estimation, flux, distributions, and plotting.

```julia
using NeuralEstimators
using Flux             
using Distributions: InverseGamma
using CairoMakie
```

--------------------------------

### Metal M-Series GPU Support

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_replicated.md

Load Metal package for Apple Metal M-Series GPU acceleration.

```julia
using Metal
```

--------------------------------

### Lux.jl Backend Integration (Idiomatic)

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/advancedusage.md

Leverage NeuralEstimators.jl with Lux.jl using its idiomatic functional style. This involves explicit management of parameters and states.

```julia
using NeuralEstimators, Lux, Random, Optimisers

network    = Lux.Chain(...)
estimator  = PointEstimator(network)

# Initialize the parameters/states
rng        = Random.default_rng()
ps, st     = Lux.setup(rng, estimator)

# Training
optimiser  = Adam(5e-4)
trainstate = Lux.Training.TrainState(estimator, ps, st, optimiser)
trainstate = train(trainstate, sampler, simulator)
ps         = trainstate.parameters
st         = trainstate.states

assess(estimator, θ_test, Z_test, ps, st)
estimate(estimator, Z, ps, st)
```

--------------------------------

### Construct a Unet Neural Estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_gridded_nonstationary.md

Initializes a Unet network and wraps it in a PointEstimator. This is the first step in setting up a neural estimator.

```julia
network = Unet(1)
estimator = PointEstimator(network)
```

--------------------------------

### Constructing a RatioEstimator with expert summaries

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Create a RatioEstimator using expert summaries. This estimator is useful for likelihood ratio estimation tasks.

```julia
estimator = RatioEstimator(summary_network, d; num_summaries = num_summaries)
```

--------------------------------

### Download Pre-simulated Draws

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

Downloads pre-simulated draws from the joint distribution p(θ, Z) using a provided URL. These draws are essential for training and evaluating neural networks.

```python
# Download the draws from p(θ, Z)
workshop_url = "https://raw.githubusercontent.com/msainsburydale/NeuralIncompleteData/main/workshop/"
download(workshop_url * "parameter_draws.rds", "parameter_draws.rds")
download(workshop_url * "data_draws.rds", "data_draws.rds")

# Load the draws into Julia
theta = load("parameter_draws.rds")
Z = load("data_draws.rds")
size(Z)
```

--------------------------------

### Train Neural Estimator with Pre-simulated Data

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_spatiotemporal.md

Train a neural estimator using pre-simulated training and validation datasets. This approach can be faster if the simulator is computationally expensive.

```julia
K = 2500 # size of the training set
θ_train = sampler(K)
θ_val   = sampler(K)
Z_train = simulator(θ_train);
Z_val   = simulator(θ_val);
estimator = train(estimator, θ_train, θ_val, Z_train, Z_val)
```

--------------------------------

### Save and Load Lux Estimator State

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/advancedusage.md

Save and load the parameters and states of a Lux-based neural estimator using BSON. Initialize an estimator with the same architecture before loading the parameters and states.

```julia
using Lux
using BSON: @save, @load

# Save
@save "estimator.bson" parameters=estimator.ps states=estimator.st

# Load (initialise an estimator with the same architecture, then load the parameters/states)
@load "estimator.bson" parameters states
estimator = Lux.setparam(estimator, parameters)
```

--------------------------------

### Create Posterior Estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Initializes a PosteriorEstimator with a specified summary network, distribution, and number of summaries.

```julia
estimator = PosteriorEstimator(summary_network, d; q = GaussianMixture, num_summaries = num_expert_summaries + num_learned_summaries)
```

--------------------------------

### Load Pretrained Network

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

Optionally load a pretrained network for the NPE. This involves downloading the model state and loading it into the NPE object.

```python
# Optional: Load pretrained network
estimator_file = "NPE.bson"
download(workshop_url * estimator_file, estimator_file)
model_state = Flux.state(NPE)
@load estimator_file model_state
NPE = Flux.loadmodel!(NPE, model_state);
```

--------------------------------

### LuxEstimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/API/estimators.md

A convenience wrapper for Lux.jl-based estimators, bundling the estimator with its parameters and states.

```APIDOC
## LuxEstimator

`LuxEstimator` bundles a Lux-based estimator together with its parameters and states for a unified, backend-agnostic API.
```

--------------------------------

### Sampling Parameters and Simulating Data for Gaussian Process

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_missing_censored.md

This snippet includes functions for sampling parameters from the prior distribution and for simulating data marginally from the Gaussian process model. It handles the computation of Cholesky factors and reshaping the simulated data.

```julia
# Constructor for above struct
function sample(K::Integer, ξ)

	# Sample parameters from the prior
	Π = ξ.Π
	tau = rand(Π.τ, K)
	rho = rand(Π.ρ, K)
	nu = 1 # fixed smoothness

	# Compute Cholesky factors
	L = map(1:K) do k
		C = matern.(UpperTriangular(ξ.D), rho[k], nu, 1)
		L = cholesky(Symmetric(C)).L
		convert(Array, L)
	end
	L = Base.stack(L)

	# Concatenate into matrix
	theta = permutedims(hcat(tau, rho))

	Parameters(theta, L)
end

# Marginal simulation from the data model
function simulate(parameters::Parameters, m::Integer)

	K = size(parameters, 2)
	tau = parameters.θ[1, :]
	L = parameters.L
	G = isqrt(size(L, 1)) # side-length of grid

	Z = map(1:K) do k
		z = simulategaussian(L[:, :, k], m)
		z = z + tau[k] * randn(size(z)...)
		z = reshape(z, G, G, 1, :)
		z
	end

	return Z
end

```

--------------------------------

### Sampler for Gaussian Process Parameters

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_gridded.md

A constructor function `sampler` that generates `K` samples from the prior distribution of the range parameter `θ` and wraps them in a `NamedMatrix`.

```julia
function sampler(K::Integer)
    θ = 0.5 * rand(K)               # K samples from p(θ) = Unif(0, 0.5)
    Parameters(NamedMatrix(θ = θ))  # Wrap as a named matrix and pass to matrix constructor
end
```

--------------------------------

### NVIDIA GPU Support

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_replicated.md

Load CUDA and cuDNN packages for NVIDIA GPU acceleration.

```julia
using CUDA, cuDNN
```

--------------------------------

### Train Estimator with Simulation On-the-fly

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_temporal.md

Trains an estimator using simulation on-the-fly. This method is useful when generating training data is computationally intensive and can be done in parallel.

```julia
estimator = train(estimator, sampler, simulator; K = 10000, simulator_args = T)
```

--------------------------------

### Parameters Constructor from Data

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_irregularspatial.md

Define a constructor for `Parameters` that accepts pre-computed parameter matrix `θ` and spatial locations `S`. It calculates covariance matrices, Cholesky factors, and spatial graphs.

```julia
function Parameters(θ::Matrix, S)
	# Covariance matrices and Cholesky factors
	L = Folds.map(axes(θ, 2)) do k
		D = pairwise(Euclidean(), S[k], dims = 1)
		Σ = Symmetric(exp.(-D ./ θ[k]))
		cholesky(Σ).L
	end

	# Spatial graphs
	g = spatialgraph.(S)

	Parameters(θ, L, g, S)
end
```

--------------------------------

### Include New Julia Files

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/CONTRIBUTING.md

New .jl files must be included in the main module file. Files in subdirectories like Estimators/ and ApproximateDistributions/ are auto-included alphabetically.

```julia
include("Estimators/Estimators.jl")
include("Estimators/Ensemble.jl")
include("Estimators/PointEstimator.jl")
include("Estimators/PosteriorEstimator.jl")
include("Estimators/QuantileEstimator.jl")
include("Estimators/RatioEstimator.jl")

include("ApproximateDistributions/ApproximateDistributions.jl")
include("ApproximateDistributions/GaussianMixture.jl")
include("ApproximateDistributions/NormalisingFlow.jl")

include("train.jl")
include("assess.jl")
include("DataParameters.jl")
include("Architectures.jl")
include("inference.jl")
include("summarystatistics.jl")
include("losses.jl")
include("missingdata.jl")
include("utility.jl")
```

--------------------------------

### GPU Performance with CUDA and Reactant

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/advancedusage.md

Ensure correct loading order of CUDA.jl and Reactant.jl for GPU acceleration. Load CUDA.jl/cuDNN.jl before Reactant.jl when using both.

```julia
using CUDA, cuDNN
using Reactant
Reactant.set_default_backend("gpu")
```

--------------------------------

### Lux package dependency

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Import the Lux package for neural network operations when using the Lux backend.

```julia
using Lux, Enzyme
```

--------------------------------

### Spatio-temporal Field Simulator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_spatiotemporal.md

Simulates a spatio-temporal field over a regular grid. It can optionally return the first differences of the field. Supports different input formats for parameters and time steps.

```julia
function simulator(θ::AbstractVector, T::Integer; grid_size = 10, first_difference::Bool = true)
    ρ = θ["θ₁"]
    φ = θ["θ₂"]

    # Regular spatial grid
    locs = [(i/grid_size, j/grid_size) for i in 1:grid_size for j in 1:grid_size]
    n_spatial = length(locs)

    # Spatial covariance matrix and its Cholesky factor
    Σ = [matern15(norm(collect(locs[i]) .- collect(locs[j])), ρ) for i in 1:n_spatial, j in 1:n_spatial]
    L = cholesky(Σ + 1e-6 * I).L

    # Dynamic system
    Z = zeros(n_spatial, T)
    Z[:, 1] = L * randn(n_spatial)
    for t in 2:T
        Z[:, t] = φ .* Z[:, t-1] .+ sqrt(1 - φ^2) .* (L * randn(n_spatial))
    end

    # Reshape into format required by our chosen architecture
    Z = reshape(Z, grid_size, grid_size, 1, T)
    return first_difference ? firstdifference(Z) : Z
end
simulator(θ::AbstractVector, T; kwargs...) = simulator(θ, rand(T); kwargs...)
simulator(θ::AbstractMatrix, T = 10:30; kwargs...) = Folds.map(ϑ -> simulator(ϑ, T; kwargs...), eachcol(θ))
```

--------------------------------

### Constructing a PosteriorEstimator with expert summaries

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Create a PosteriorEstimator using expert summaries. This estimator is suitable for approximating posterior distributions, here using a GaussianMixture.

```julia
estimator = PosteriorEstimator(summary_network, d; num_summaries = num_summaries, q = GaussianMixture)
```

--------------------------------

### Initialize Posterior Estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

Initialize the estimator using PosteriorEstimator, which is designed for full posterior inference. This replaces the PointEstimator used for neural point estimation.

```python
NPE = PosteriorEstimator(network, d);
```

--------------------------------

### Instantiate DeepSet Network

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_gridded.md

Create a DeepSet object by combining the inner (ψ) and outer (ϕ) network components.

```julia
# DeepSet object
network = DeepSet(ψ, ϕ)
```

--------------------------------

### Train a Neural Estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_gridded_nonstationary.md

Trains the initialized estimator using training and validation data generated by a sampler. Requires a simulator function.

```julia
K = 2500
θ_train = sampler(K)
θ_val   = sampler(K)
estimator = train(estimator, θ_train, θ_val, simulator)
```

--------------------------------

### Clone NeuralEstimators.jl Repository

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/README.md

Use this command to download the package source code from GitHub.

```bash
git clone https://github.com/msainsburydale/NeuralEstimators.jl.git
```

--------------------------------

### Parameters Constructor with Cholesky Factor Computation

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_gridded.md

Constructs `Parameters` by calculating the Cholesky factor of the covariance matrix for a given parameter `θ`. Assumes a regular grid for spatial locations.

```julia
function Parameters(θ::AbstractMatrix; grid_dim = 16)
	# Spatial locations: regular grid over the unit square
	pts = range(0, 1, length = grid_dim)
	S = expandgrid(pts, pts)

	# Pairwise distances, covariance matrices, and Cholesky factors
	D = pairwise(Euclidean(), S, dims = 1)
	K = size(θ, 2)
	L = Folds.map(1:K) do k
		Σ = exp.(-D ./ θ[k])
		cholesky(Symmetric(Σ)).L
	end

	Parameters(θ, L)
end
```

--------------------------------

### Matrix Parameter Data Simulator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_replicated.md

Simulate data for multiple parameter vectors stored in a matrix (θ), with a default range for the number of replicates (m).

```julia
simulator(θ::AbstractMatrix, m = 10:100) = [simulator(ϑ, m) for ϑ in eachcol(θ)]
```

--------------------------------

### AMD ROCm GPU Support

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_replicated.md

Load AMDGPU package for AMD ROCm GPU acceleration.

```julia
using AMDGPU
```

--------------------------------

### Sample Prior and Simulate Censored Data

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_missing_censored.md

Defines functions to sample parameters from a prior distribution and simulate censored data based on a given parameter set and sample size. Includes transformations to uniform margins.

```julia
# Libraries used throughout this example
using NeuralEstimators, Flux
using Folds
using CUDA # GPU if it is available
using LinearAlgebra: Symmetric, cholesky
using Distributions: cdf, Uniform, Normal, quantile
using CairoMakie   

# Sampling θ from the prior distribution
function sample(K)
	ρ = rand(Uniform(-0.99, 0.99), K)
	δ = rand(Uniform(0.0, 1.0), K)
	θ = vcat(ρ', δ')
	return θ 
end

# Marginal simulation of Z | θ
function simulate(θ, m) 
	Z = Folds.map(1:size(θ, 2)) do k
		ρ = θ[1, k]
		δ = θ[2, k]
		Σ = [1 ρ; ρ 1]
		L = cholesky(Symmetric(Σ)).L
		X = L * randn(2, m)                 # Standard Gaussian margins
		X = -log.(1 .- cdf.(Normal(), X))   # Transform to unit exponential margins
		R = -log.(1 .- rand(1, m))         
		Y = δ .* R .+ (1 - δ) .* X          
		Z = F.(Y; δ = δ)                     # Transform to uniform margins
	end
	return Z
end

# Marginal distribution function; see Huser and Wadsworth (2019)
function F(y; δ)
	if δ == 0.5 
        u = 1 .- exp.(- 2 .* y) .* (1 .+ 2 .* y) 
    else 
        u = 1 .- (δ ./ (2 .* δ .- 1)) .* exp.(- y ./ δ) .+ ((1 .- δ) ./ (2 * δ .- 1)) .* exp.( - y ./ (1 - δ)) 
    end
	return u
end
```

--------------------------------

### Training the neural estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Train the estimator using the provided sampler and simulator functions. The training process adapts to a range of sample sizes specified in `n_training`.

```julia
# Sample sizes used during training
n_training = 30:1000

estimator = train(estimator, sampler, simulator; simulator_args = (n_training,))
```

--------------------------------

### Constructing a PointEstimator with expert summaries

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Create a PointEstimator by providing the summary network and specifying the number of expert summaries. This estimator provides a single point estimate.

```julia
estimator = PointEstimator(summary_network, d; num_summaries = num_summaries)
```

--------------------------------

### samplesize

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/API/architectures.md

User-defined summary statistic for sample size.

```APIDOC
## samplesize

### Description
A user-defined summary statistic function, often useful as a summary in `DeepSet` objects.

### Usage
Corresponds to a summary statistic that can be used in user-defined summaries.
```

--------------------------------

### Download and Load Sea Ice Extent Data

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/tutorials/introduction.ipynb

Downloads a sea ice extent dataset from a specified URL and loads it into a Julia environment. This snippet requires the `workshop_url` and `dest_file` variables to be defined, and the `download` and `load` functions to be available.

```python
# Name of the file
dest_file = "sea_ice_extent_complete.rds"

# Download the file
download(workshop_url * dest_file, dest_file)

# Read the .rds object into Julia
sea_ice_extent = load(dest_file)

# Inspect the object
println("Size: ", size(sea_ice_extent))
println("Type: ", typeof(sea_ice_extent))
```

--------------------------------

### Construct Neural Bayes Estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_missing_censored.md

Initializes a DeepSet network and a PointEstimator, then trains the neural Bayes estimator using simulated data. Requires `Chain`, `Conv`, `Dense`, `relu`, `flatten`, `DeepSet`, `PointEstimator`, `train`, `sample`, `simulate`, and `simulator_args`.

```julia
# Construct DeepSet object
ψ = Chain(
	Conv((10, 10), 1 => 16,  relu),
	Conv((5, 5),  16 => 32,  relu),
	Conv((3, 3),  32 => 64, relu),
	flatten
	)
ϕ = Chain(
	Dense(64, 256, relu),
	Dense(256, d, exp)
	)
network = DeepSet(ψ, ϕ)

# Initialise point estimator
θ̂ = PointEstimator(network)

# Train neural Bayes estimator
H = 50
θ̂ = train(θ̂, sample, simulate, simulator_args = H, ξ = ξ, K = 1000, epochs = 10)
```

--------------------------------

### Univariate Data Simulator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_replicated.md

Simulate univariate data for a single parameter vector (θ) with a fixed number of replicates (m).

```julia
simulator(θ::AbstractVector, m::Integer) = θ["μ"] .+ θ["σ"] .* randn(1, m)
```

--------------------------------

### Train the Estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Trains the estimator using the provided sampler and simulator functions.

```julia
estimator = train(estimator, sampler, simulator)
```

--------------------------------

### Save and Load Flux Estimator State

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/advancedusage.md

Save and load the model state of a Flux-based neural estimator using BSON. Initialize an estimator with the same architecture before loading the state.

```julia
using Flux
using BSON: @save, @load

# Save
model_state = Flux.state(estimator)
@save "estimator.bson" model_state

# Load (initialise an estimator with the same architecture, then load the state)
@load "estimator.bson" model_state
Flux.loadmodel!(estimator, model_state)
```

--------------------------------

### Package Dependencies for Neural Estimators

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_irregularspatial.md

Load necessary packages for neural estimation, GNNs, distributions, plotting, distances, parallel computation, and linear algebra.

```julia
using NeuralEstimators
using Flux
using GraphNeuralNetworks
using Distributions: Uniform
using CairoMakie
using Distances
using Folds                    # parallel simulation (start Julia with --threads=auto)
using LinearAlgebra            # Cholesky factorisation
using Statistics: mean
```

--------------------------------

### Random Size Data Simulator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_replicated.md

Simulate data with a random number of replicates (m) drawn from a range for a single parameter vector (θ).

```julia
simulator(θ::AbstractVector, m) = simulator(θ, rand(m))
```

--------------------------------

### Create Point Estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Initializes a PointEstimator with a specified summary network and number of summaries.

```julia
estimator = PointEstimator(summary_network, d; num_summaries = num_expert_summaries + num_learned_summaries)
```

--------------------------------

### Construct Neural Network Summary

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Defines the summary network and specifies the total number of summaries (expert + learned).

```julia
num_expert_summaries  = 2
num_learned_summaries = 3d
summary_network = Chain(Dense(n, 64, gelu), Dense(64, 64, gelu), Dense(64, num_learned_summaries))
```

--------------------------------

### logsamplesize

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/API/architectures.md

User-defined summary statistic for the log of sample size.

```APIDOC
## logsamplesize

### Description
A user-defined summary statistic function for the log of sample size, often useful as a summary in `DeepSet` objects.

### Usage
Corresponds to a summary statistic that can be used in user-defined summaries.
```

--------------------------------

### Compress

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/API/architectures.md

An output layer for compressing outputs to satisfy constraints.

```APIDOC
## Compress

### Description
An output layer that can be used at the end of a neural network to ensure the outputs satisfy certain constraints. This type should be incorporated as a separate layer in the final stage of a `Chain`.

### Usage
Can be used as an output layer to enforce output constraints.
```

--------------------------------

### Simulate and Visualize Spatio-Temporal Dependence

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_spatiotemporal.md

Simulates spatio-temporal data with varying temporal dependence and visualizes the results. This is useful for understanding the effect of the temporal range parameter (θ₂).

```julia
# Fix θ₁, vary θ₂
θ = NamedMatrix(
    θ₁ = [0.3, 0.3],
    θ₂ = [0.1, 0.8]
)

# Simulate several time steps
Z = simulator(θ, 5; grid_size = 32, first_difference = false)

labels = [
    "Low temporal dependence (θ₂ = 0.1)",
    "High temporal dependence (θ₂ = 0.8)"
]

fig = Figure(size = (1350, 600));

for k in 1:2
    Z_k = Z[k]
    T = size(Z_k, 4)

    for t in 1:T
        ax = Axis(
            fig[k, t],
            title = t == 1 ? labels[k] : "t = $t",
            aspect = DataAspect()
        )

        hidedecorations!(ax)

        hm = heatmap!(
            ax,
            Z_k[:, :, 1, t],
            colormap = :balance,
            colorrange = (-2.5, 2.5)
        )

        if t == T
            Colorbar(fig[k, T + 1], hm)
        end
    end
end

display(fig)
```

--------------------------------

### Construct and initialize DeepSet estimator

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_missing_censored.md

Builds a neural network using the DeepSet architecture, where τ is incorporated as an input to the outer network. Initializes a PointEstimator with this network.

```julia
# Construct neural network based on DeepSet architecture
ψ = Chain(Dense(n * 2, w, relu),Dense(w, w, relu))    
ϕ = Chain(Dense(w + 1, w, relu), final_layer)
network = DeepSet(ψ, ϕ)

# Initialise the estimator
estimator = PointEstimator(network)
```

--------------------------------

### Apply Estimator to Observed Data (Posterior/Ratio)

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Applies a trained posterior or ratio estimator to observed data to sample from the posterior distribution.

```julia
θ = sampler(1)                 # ground truth (not known in practice)
Z = simulator(θ)               # stand-in for real observations

sampleposterior(estimator, Z)      # posterior sample
```

--------------------------------

### GNNSummary

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/API/architectures.md

A module for graph neural network summaries.

```APIDOC
## GNNSummary

### Description
A module for graph neural network summaries, useful when constructing neural estimators.

### Usage
This module can be used as a component in neural network construction.
```

--------------------------------

### Define Sampling and Simulation Functions

Source: https://github.com/msainsburydale/neuralestimators.jl/blob/main/docs/src/examples/data_expert_summaries.md

Defines functions for sampling from the prior and simulating data, including expert summary statistics computation.

```julia
d = 2    # dimension of θ
n = 100  # number of replicates (fixed in this example)

# Function to sample from the prior
sampler(K) = NamedMatrix(μ = randn(K), σ = rand(K))

# Functions to simulate data
function simulator(θ::AbstractVector)
    Z = θ["μ"] .+ θ["σ"] .* sort(randn(n))
    return Z
end
function simulator(θ::AbstractMatrix)
    Z = reduce(hcat, map(simulator, eachcol(θ)))
    S = vcat(mean(Z, dims = 1), std(Z, dims = 1))
    return DataAndSummaries(Z, S)
end
```