### SimulatedAnnealing Sampler Setup

Source: https://context7.com/fmagrini/bayes-bay/llms.txt

Implements simulated annealing with a decreasing temperature schedule during burn-in to escape local minima. Requires setting start temperature and cooling fraction.

```python
import bayesbay as bb

# Setup simulated annealing sampler
sa_sampler = bb.samplers.SimulatedAnnealing(
    temperature_start=10.0,
    cooling_fraction=0.8
)

inversion = bb.BayesianInversion(
    parameterization=parameterization,
    log_likelihood=log_likelihood,
    n_chains=10
)

# Burn-in is important for simulated annealing
inversion.run(
    sampler=sa_sampler,
    n_iterations=100000,
    burnin_iterations=30000,
    save_every=100,
    verbose=True
)
```

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

### Install BayesBay package

Source: https://github.com/fmagrini/bayes-bay/blob/main/README.md

Install the package in the current environment.

```console
$ python -m pip install .
```

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

### Install BayesBay from Source

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/installation.md

Build and install BayesBay from its source code. This is useful for developers who need to modify the code. The `-e` flag installs in editable mode.

```bash
git clone https://github.com/fmagrini/bayes-bridge
cd bayes-bridge
pip install -e .
```

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

### Set up development environment with venv

Source: https://github.com/fmagrini/bayes-bay/blob/main/README.md

Use these commands to create a virtual environment and install dependencies via pip.

```console
$ python -m venv bayesbay_dev
$ source bayesbay_dev/bin/activate
$ pip install -r envs/requirements_dev.txt
```

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

### ParallelTempering Sampler Setup

Source: https://context7.com/fmagrini/bayes-bay/llms.txt

Implements parallel tempering for improved exploration of multi-modal posteriors. Requires setting temperature parameters and swap frequency. More chains are recommended.

```python
import bayesbay as bb

# Setup parallel tempering sampler
pt_sampler = bb.samplers.ParallelTempering(
    temperature_max=5.0,
    chains_with_unit_temperature=0.4,
    swap_every=500
)

inversion = bb.BayesianInversion(
    parameterization=parameterization,
    log_likelihood=log_likelihood,
    n_chains=12
)

inversion.run(
    sampler=pt_sampler,
    n_iterations=200000,
    burnin_iterations=50000,
    save_every=200,
    verbose=True
)

# Check temperatures after sampling
for chain in inversion.chains:
    print(f"Chain {chain.id}: T={chain.temperature:.2f}")
```

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

### Inversion Output Statistics

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/12_transd_gaussian_mixture.ipynb

Example output showing the statistics for individual chains after running the inversion.

```text
Output:
Chain ID: 0
TEMPERATURE: 1
EXPLORED MODELS: 450000
ACCEPTANCE RATE: 65197/450000 (14.49 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(my_param_space): 2825/75051 (3.76%)
	DeathPerturbation(my_param_space): 2822/73838 (3.82%)
	NoisePerturbation(my_data): 15740/74627 (21.09%)
	ParamPerturbation(['my_param_space.mean', 'my_param_space.std', 'my_param_space.weight']): 43810/226484 (19.34%)
Chain ID: 1
TEMPERATURE: 1
EXPLORED MODELS: 450000
ACCEPTANCE RATE: 66472/450000 (14.77 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(my_param_space): 2908/75077 (3.87%)
	DeathPerturbation(my_param_space): 2906/74040 (3.92%)
	NoisePerturbation(my_data): 15748/75349 (20.90%)
	ParamPerturbation(['my_param_space.mean', 'my_param_space.std', 'my_param_space.weight']): 44910/225534 (19.91%)
Chain ID: 2
TEMPERATURE: 1
EXPLORED MODELS: 450000
ACCEPTANCE RATE: 59628/450000 (13.25 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(my_param_space): 2613/75539 (3.46%)
	DeathPerturbation(my_param_space): 2611/74536 (3.50%)
	NoisePerturbation(my_data): 15146/74962 (20.20%)
	ParamPerturbation(['my_param_space.mean', 'my_param_space.std', 'my_param_space.weight']): 39258/224963 (17.45%)
Chain ID: 3
TEMPERATURE: 1
EXPLORED MODELS: 450000
ACCEPTANCE RATE: 60787/450000 (13.51 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(my_param_space): 2439/74394 (3.28%)
	DeathPerturbation(my_param_space): 2438/74830 (3.26%)
	NoisePerturbation(my_data): 14820/75220 (19.70%)
	ParamPerturbation(['my_param_space.mean', 'my_param_space.std', 'my_param_space.weight']): 41090/225556 (18.22%)
Chain ID: 4
TEMPERATURE: 1
EXPLORED MODELS: 450000
ACCEPTANCE RATE: 59737/450000 (13.27 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(my_param_space): 2559/75227 (3.40%)
	DeathPerturbation(my_param_space): 2557/74293 (3.44%)
	NoisePerturbation(my_data): 15341/75074 (20.43%)
	ParamPerturbation(['my_param_space.mean', 'my_param_space.std', 'my_param_space.weight']): 39280/225406 (17.43%)
Chain ID: 5
TEMPERATURE: 1
EXPLORED MODELS: 450000
ACCEPTANCE RATE: 58972/450000 (13.10 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(my_param_space): 2478/74763 (3.31%)
	DeathPerturbation(my_param_space): 2475/74234 (3.33%)
	NoisePerturbation(my_data): 15324/75107 (20.40%)
	ParamPerturbation(['my_param_space.mean', 'my_param_space.std', 'my_param_space.weight']): 38695/225896 (17.13%)
Chain ID: 6
TEMPERATURE: 1
EXPLORED MODELS: 450000
ACCEPTANCE RATE: 66545/450000 (14.79 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(my_param_space): 2615/75104 (3.48%)
	DeathPerturbation(my_param_space): 2614/74131 (3.53%)
	NoisePerturbation(my_data): 15354/75098 (20.45%)
	ParamPerturbation(['my_param_space.mean', 'my_param_space.std', 'my_param_space.weight']): 45962/225667 (20.37%)
Chain ID: 7
TEMPERATURE: 1
EXPLORED MODELS: 450000
ACCEPTANCE RATE: 64728/450000 (14.38 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(my_param_space): 2779/74939 (3.71%)
	DeathPerturbation(my_param_space): 2778/73787 (3.76%)
	NoisePerturbation(my_data): 15585/75336 (20.69%)
	ParamPerturbation(['my_param_space.mean', 'my_param_space.std', 'my_param_space.weight']): 43586/225938 (19.29%)
Chain ID: 8
TEMPERATURE: 1
EXPLORED MODELS: 450000
ACCEPTANCE RATE: 58536/450000 (13.01 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(my_param_space): 2444/74991 (3.26%)
	DeathPerturbation(my_param_space): 2444/74991 (3.26%)
```

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

### Example Chain Statistics Output

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/42_transd_partition_mod.ipynb

This is an example output format for the statistics of a single chain after running the Bayesian inversion. It details acceptance rates for the overall chain and for specific perturbation types.

```text
Chain ID: 0
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 118583/400000 (29.65 %)
PARTIAL ACCEPTANCE RATES:
	BirthPerturbation(voronoi): 3968/56793 (6.99%)
	DeathPerturbation(voronoi): 3963/57442 (6.90%)
	NoisePerturbation(d_obs): 12926/57028 (22.67%)
	ParamPerturbation(voronoi.discretization): 8612/57154 (15.07%)
	ParamPerturbation(voronoi.y): 89114/171583 (51.94%)
```

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

### Install BayesBay using pip

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/installation.md

Install the BayesBay package directly from PyPI using pip. This is the recommended method for most users.

```bash
pip install bayesbay
```

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

### Example Bayesian Sampling Output

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/01_polyfit.ipynb

This is an example output from running the `print_statistics()` method on Bayesian inversion chains. It shows key metrics like explored models, acceptance rates, and partial acceptance rates for parameter perturbations.

```text
Output:
Chain ID: 0
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5655/50000 (11.31 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5655/50000 (11.31%)
Chain ID: 1
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5393/50000 (10.79 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5393/50000 (10.79%)
Chain ID: 2
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5542/50000 (11.08 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5542/50000 (11.08%)
Chain ID: 3
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5662/50000 (11.32 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5662/50000 (11.32%)
Chain ID: 4
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5580/50000 (11.16 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5580/50000 (11.16%)
Chain ID: 5
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5642/50000 (11.28 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5642/50000 (11.28%)
Chain ID: 6
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5392/50000 (10.78 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5392/50000 (10.78%)
Chain ID: 7
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5626/50000 (11.25 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5626/50000 (11.25%)
Chain ID: 8
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5714/50000 (11.43 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5714/50000 (11.43%)
Chain ID: 9
TEMPERATURE: 1
EXPLORED MODELS: 50000
ACCEPTANCE RATE: 5468/50000 (10.94 %)
PARTIAL ACCEPTANCE RATES:
	ParamPerturbation(['my_param_space.m0', 'my_param_space.m1', 'my_param_space.m2', 'my_param_space.m3']): 5468/50000 (10.94%)
```

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

### Install Package in Editable Mode

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/developer.md

Install the BayesBay package in editable mode for development.

```console
$ python -m pip install -e .
```

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

### Define True Model Parameters

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Sets up the ground truth values for the inversion example.

```python
m_true = 10
sigma_true = 0.5
d = np.random.normal(m_true, sigma_true)
```

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

### Setup Inversion Target

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/21_rayleigh.ipynb

Defines the target data and uncertainty parameters for the Bayesian inversion.

```python
target = bb.Target("rayleigh", 
                   d_obs, 
                   std_min=0.001, 
                   std_max=0.1, 
                   std_perturb_std=0.002)
```

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

### Initialize and Print Model State

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Initializes the model state using parameterization and prints it to show the current numerical information. This is the starting point for MCMC iterations.

```python
state = parameterization.initialize()
print(state)
```

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

### VanillaSampler - Standard MCMC

Source: https://context7.com/fmagrini/bayes-bay/llms.txt

Uses the default VanillaSampler for standard reversible-jump MCMC. Suitable for well-behaved posteriors. Requires prior setup of parameterization and log_likelihood.

```python
import bayesbay as bb

# Setup inversion (using previous definitions)
# parameterization = ...
# log_likelihood = ...

inversion = bb.BayesianInversion(
    parameterization=parameterization,
    log_likelihood=log_likelihood,
    n_chains=8
)

# Use vanilla sampler (default)
sampler = bb.samplers.VanillaSampler()

inversion.run(
    sampler=sampler,
    n_iterations=100000,
    burnin_iterations=20000,
    save_every=100,
    verbose=True,
    print_every=5000,
    parallel_config={"n_jobs": 8}  # parallel chains
)

# Print chain statistics
for chain in inversion.chains:
    chain.print_statistics()
```

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

### Combine Parameter Spaces into a Parameterization

Source: https://context7.com/fmagrini/bayes-bay/llms.txt

Aggregate multiple ParameterSpace instances into a single Parameterization object to define the complete model parameterization. This example combines background and anomaly parameter spaces.

```python
import bayesbay as bb
import numpy as np

# Define multiple parameter spaces for a complex model
# Space 1: Background model parameters
background_velocity = bb.prior.UniformPrior(name="v_bg", vmin=3.0, vmax=4.5, perturb_std=0.1)
background_space = bb.parameterization.ParameterSpace(
    name="background",
    n_dimensions=1,
    parameters=[background_velocity]
)

# Space 2: Anomaly parameters (trans-dimensional)
anomaly_amplitude = bb.prior.GaussianPrior(name="amplitude", mean=0, std=0.5, perturb_std=0.1)
anomaly_space = bb.parameterization.ParameterSpace(
    name="anomalies",
    n_dimensions=None,
    n_dimensions_min=0,
    n_dimensions_max=10,
    parameters=[anomaly_amplitude]
)

# Combine into full parameterization
parameterization = bb.parameterization.Parameterization([background_space, anomaly_space])

# Initialize random state from prior
initial_state = parameterization.initialize()
print(f"Background velocity: {initial_state['background']['v_bg']}")
print(f"Number of anomalies: {initial_state['anomalies'].n_dimensions}")
```

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

### Set up development environment with mamba

Source: https://github.com/fmagrini/bayes-bay/blob/main/README.md

Use this command to create a development environment from the provided YAML file.

```console
$ mamba env create -f envs/environment_dev.yml
```

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

### Set up target and log-likelihood

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/02_hierarchical_polyfit.ipynb

Configures the data target with noise bounds and initializes the log-likelihood function.

```python
target = bb.Target("my_data", 
                   y_noisy, 
                   std_min=0, 
                   std_max=40, 
                   std_perturb_std=4,
                   noise_is_correlated=False)
```

```python
log_likelihood = bb.LogLikelihood(targets=target, fwd_functions=fwd_function)
```

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

### Import libraries and define constants

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/02_hierarchical_polyfit.ipynb

Initializes the environment by importing necessary libraries and setting up simulation constants.

```python
import bayesbay as bb
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(30)
```

```python
# dimensions and true coefficients
N_DIMS = 4
M0, M1, M2, M3 = 20, -10, -3, 1

# data and noise
N_DATA = 15
DATA_X = np.linspace(-5, 10, N_DATA)
DATA_NOISE_STD = 20
```

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

### Import libraries and define constants

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/31_sw_tomography.ipynb

Initializes the environment by importing necessary libraries for tomography, discretization, and Bayesian inference.

```python
import numpy as np
import matplotlib.pyplot as plt
import scipy

from seislib.tomography import SeismicTomography
import seislib.colormaps as scm

from bayesbay.discretization import Voronoi2D
from bayesbay.prior import UniformPrior
import bayesbay as bb

np.random.seed(10)
```

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

### Import libraries and define constants

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/11_gaussian_mixture.ipynb

Initializes the environment by importing necessary libraries and setting global constants for the Gaussian mixture.

```python
import bayesbay as bb
from math import sqrt, pi
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(30)
```

```python
MEANS = [140, 162, 177]  # Means of the Gaussians
STDS = [12, 5, 6]  # Standard deviations of the Gaussians
WEIGHTS = [0.4, 0.3, 0.3]  # Weights of each Gaussian in the mixture
N_SAMPLES = 10_000 # Number of samples to generate
```

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

### Define CustomPrior with User-Defined Functions

Source: https://context7.com/fmagrini/bayes-bay/llms.txt

Create a CustomPrior by providing your own log-prior and sampling functions for maximum flexibility. This example uses a log-uniform (Jeffreys) prior.

```python
import bayesbay as bb
import numpy as np

# Define custom log-prior function (log-uniform/Jeffreys prior)
def log_prior_jeffreys(value, position=None):
    if value <= 0:
        return -np.inf
    return -np.log(value)

# Define custom sampling function
def sample_jeffreys(position=None):
    # Sample from log-uniform between 0.1 and 10
    return np.exp(np.random.uniform(np.log(0.1), np.log(10)))

# Create custom prior
jeffreys_prior = bb.prior.CustomPrior(
    name="scale_parameter",
    log_prior=log_prior_jeffreys,
    sample=sample_jeffreys,
    perturb_std=0.3
)

# Use in parameterization
param_space = bb.parameterization.ParameterSpace(
    name="custom_space",
    n_dimensions=1,
    parameters=[jeffreys_prior]
)
```

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

### Configure Pre-commit Hook for Formatting and Docs

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/developer.md

Set up a pre-commit hook to automatically format code with black and build documentation for notebooks.

```sh
#!/bin/sh

black src/bayesbay

if git diff --cached --name-status | grep -q 'docs/source/tutorials/*.ipynb'; then
     cd docs
     make html
fi
```

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

### Import libraries and define constants

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/42_transd_partition_mod.ipynb

Initializes the environment by importing necessary modules and setting a random seed for reproducibility.

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from bayesbay.discretization import Voronoi1D
from bayesbay.prior import UniformPrior
from bayesbay.parameterization import Parameterization
from bayesbay import Target, LogLikelihood, BayesianInversion
np.random.seed(30)
```

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

### Initialize and Run Bayesian Inversion

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Initializes the BayesianInversion instance with a log-likelihood and parameterization, then executes the sampling process and prints statistics for each chain.

```python
inversion = bb.BayesianInversion(
    log_likelihood=log_likelihood, 
    parameterization=parameterization, 
    n_chains=10
)
inversion.run(
    n_iterations=50_000, 
    burnin_iterations=5_000, 
    save_every=50, 
    verbose=False, 
)
for chain in inversion.chains:
    chain.print_statistics()
```

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

### Import BayesBay and Dependencies

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Initializes the environment with required libraries for BayesBay and plotting.

```python
import bayesbay as bb

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm

plt.rc('text', usetex=True)
plt.rc('font', family='serif')
```

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

### Initialize and run Bayesian inversion

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/11_gaussian_mixture.ipynb

Configures the BayesianInversion object with log-likelihood, parameterization, and chain count, then executes the sampling process.

```python
inversion = bb.BayesianInversion(
    log_likelihood=log_likelihood, 
    parameterization=parameterization, 
    n_chains=20
)
inversion.run(
    n_iterations=250_000, 
    burnin_iterations=50_000, 
    save_every=100, 
    verbose=False, 
)
```

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

### Plotting KDE with ArviZ

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/11_gaussian_mixture.ipynb

This code snippet visualizes kernel density estimates (KDE) with specified HDI probabilities and contour styles using ArviZ. It requires ArviZ and Matplotlib to be installed. The plot is displayed using `plt.show()`.

```python
f'{key}_3': inferred_value[:,2]},
        marginals=True,
        reference_values={f'{key}_1': true_value[0], 
                          f'{key}_2': true_value[1], 
                          f'{key}_3': true_value[2]},
        reference_values_kwargs={'color': 'yellow',
                                 'ms': 10},
        kind='kde',
        kde_kwargs={
            'hdi_probs': [0.3, 0.6, 0.9],  # Plot 30%, 60% and 90% HDI contours
            'contourf_kwargs': {'cmap': 'Blues'},
            },
        ax=axes,
        textsize=10
        )
    fig.suptitle(key.upper())
    plt.tight_layout()
    plt.show()
```

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

### Import BayesBay libraries and define constants

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/41_simple_partition_mod.ipynb

Initializes the environment by importing necessary modules and setting the random seed for reproducibility.

```python
import numpy as np
import matplotlib.pyplot as plt
from bayesbay.discretization import Voronoi1D
from bayesbay.prior import UniformPrior
from bayesbay.parameterization import Parameterization
from bayesbay import Target, LogLikelihood, BayesianInversion
np.random.seed(30)
```

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

### Import Libraries and Define Constants

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/03_transd_polyfit.ipynb

Initializes the environment by importing necessary libraries and setting constants for the true model, dimension range, and synthetic data generation.

```python
import numpy as np
import matplotlib
import bayesbay as bb
```

```python
# True model
M_TRUE = np.array([20, -10, -3, 1])

# Constants for dimensions
N_DIMS_MIN = 1
N_DIMS_MAX = 15

# Constants for data
N_DATA = 15
DATA_NOISE_STD = 20
DATA_X = np.linspace(-5, 10, N_DATA)
```

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

### Setting up Bayesian Inversion

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/22_rayleigh_love.ipynb

This code configures the BayesianInversion object with the defined parameterization, log-likelihood function, and number of chains. It then runs the inversion process with specified iterations and burn-in periods.

```python
log_likelihood = bb.LogLikelihood(
    targets=[target_rayleigh, target_love], 
    fwd_functions=[forward_rayleigh, forward_love]
)

inversion = bb.BayesianInversion(
    parameterization=parameterization, 
    log_likelihood=log_likelihood,
    n_chains=20
)

inversion.run(
    sampler=None, 
    n_iterations=350_000, 
    burnin_iterations=75_000, 
    save_every=150, 
    verbose=False
)
```

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

### Initialize Vs prior and Voronoi discretization

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/51_sw_rf_joint.ipynb

Sets up a custom initialization function for the shear wave velocity (Vs) prior and defines a 1D Voronoi tessellation for parameterizing the subsurface. The Voronoi discretization controls the number and boundaries of layers.

```python
def initialize_vs(param, positions=None):
    vmin, vmax = param.get_vmin_vmax(positions)
    sorted_vals = np.sort(np.random.uniform(vmin, vmax, positions.size))
    return sorted_vals

vs = bb.prior.UniformPrior(name="vs", 
                           vmin=[2.2, 2.8, 3.3, 4], 
                           vmax=[3.9, 4.6, 4.8, 5], 
                           position=[0, 20, 60, 150],
                           perturb_std=0.15)
vs.set_custom_initialize(initialize_vs)

voronoi = Voronoi1D(
    name="voronoi", 
    vmin=0,
    vmax=150,
    perturb_std=10,
    n_dimensions=None, 
    n_dimensions_min=4,
    n_dimensions_max=15,
    parameters=[vs], 
    birth_from='neighbour')
parameterization = bb.parameterization.Parameterization(voronoi)
```

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

### Configuring and Running the Bayesian Inversion

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/31_sw_tomography.ipynb

Initializes and runs the Bayesian inversion process. This involves specifying the parameterization, log-likelihood function, number of chains, and the number of iterations for sampling and burn-in.

```python
inversion = bb.BayesianInversion(
    parameterization=parameterization, 
    log_likelihood=log_likelihood,
    n_chains=20
)
inversion.run(
    sampler=None, 
    n_iterations=75_000, 
    burnin_iterations=25_000, 
    save_every=250, 
    verbose=False,
    print_every=25_000
)
```

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

### Initialize Joint LogLikelihood for Multiple Datasets

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Demonstrates how to initialize a LogLikelihood instance for joint inversion of multiple datasets. Each dataset (target) is associated with its own forward function.

```python
log_likelihood = bb.likelihood.LogLikelihood(
            targets=[target_1, target_2], 
            fwd_functions=[fwd_func_1, fwd_func_2]
            )
```

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

### Configure Log Likelihood

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/41_simple_partition_mod.ipynb

Sets up the target observation and links it to the forward function for likelihood evaluation.

```python
target = Target("d_obs", 
                d_obs, 
                std_min=0, 
                std_max=20, 
                std_perturb_std=2)

log_likelihood = LogLikelihood(targets=target, fwd_functions=fwd_function)
```

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

### Configure Voronoi parameterization

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/42_transd_partition_mod.ipynb

Sets up the Voronoi1D discretization with a uniform prior on the parameter values and defined dimension constraints.

```python
y = UniformPrior('y', vmin=-35, vmax=35, perturb_std=3.5)
voronoi = Voronoi1D(
    name="voronoi", 
    vmin=0,
    vmax=10,
    perturb_std=0.75,
    n_dimensions=None, 
    n_dimensions_min=2,
    n_dimensions_max=40,
    parameters=[y], 
    birth_from='prior'
)
parameterization = Parameterization(voronoi)
```

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

### Initialize ParameterSpace for Trans-dimensional Inference

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/12_transd_gaussian_mixture.ipynb

Initialize a ParameterSpace with an unknown number of dimensions to enable trans-dimensional inference. Set minimum and maximum dimensions for the parameter space.

```python
param_space = bb.parameterization.ParameterSpace(
    name='my_param_space', 
    n_dimensions=None, # Trans-dimensional setting
    n_dimensions_min=1, # Minimum number of dimensions (i.e., Gaussians in the mixture)
    n_dimensions_max=7, # Maximum number of dimensions (i.e., Gaussians in the mixture)
    parameters=[mean, std, weight], 
)
parameterization = bb.parameterization.Parameterization(param_space)
```

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

### Configure and Run Bayesian Inversion

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/12_transd_gaussian_mixture.ipynb

Initializes the BayesianInversion object and executes the sampling process with specified iterations and burn-in phases.

```python
inversion = bb.BayesianInversion(
    log_likelihood=log_likelihood, 
    parameterization=parameterization, 
    n_chains=20
)
inversion.run(
    sampler=bb.samplers.SimulatedAnnealing(temperature_start=5),
    n_iterations=450_000, 
    burnin_iterations=150_000, 
    save_every=100, 
    verbose=False, 
)
```

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

### Initialize Parameterization

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Aggregates all parameter spaces into a single parameterization instance for the inversion.

```python
parameterization = bb.parameterization.Parameterization(param_space)
print(parameterization)
```

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

### Initialize LogLikelihood

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/11_gaussian_mixture.ipynb

Combines the target and forward function into a LogLikelihood instance for Bayesian sampling.

```python
log_likelihood = bb.LogLikelihood(targets=target, fwd_functions=fwd_function)
```

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

### Setting up the Log Likelihood Function

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/31_sw_tomography.ipynb

Configures the log-likelihood function for the Bayesian inversion. It defines the observed data ('d_obs'), its associated uncertainties, and links it to the forward model function.

```python
target = bb.likelihood.Target('d_obs', 
                              d_obs, 
                              std_min=0, 
                              std_max=0.01, 
                              std_perturb_std=0.001,
                              noise_is_correlated=False)
log_likelihood = bb.likelihood.LogLikelihood(targets=target, fwd_functions=forward)
```

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

### Create virtualenv Virtual Environment

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/installation.md

Create a new virtual environment using virtualenv, specifying the Python version. Replace `<path-to-new-env>` and `<env-name>` accordingly.

```bash
virtualenv <path-to-new-env>/<env-name> -p=3.10
```

```bash
virtualenv ~/my_envs/bb_env -p=3.10
```

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

### Configure Voronoi Parameterization

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/21_rayleigh.ipynb

Initializes the Voronoi1D parameterization object for the inversion.

```python
voronoi = Voronoi1D(
    name="voronoi", 
    vmin=0,
    vmax=150,
    perturb_std=10,
    n_dimensions=None, 
    n_dimensions_min=4,
    n_dimensions_max=15,
    parameters=[vs], 
    birth_from='neighbour'
)
parameterization = bb.parameterization.Parameterization(voronoi)
```

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

### Implement 1D Voronoi Discretization

Source: https://context7.com/fmagrini/bayes-bay/llms.txt

Sets up a 1D Voronoi tessellation for trans-dimensional problems where the number of layers is unknown.

```python
import bayesbay as bb
import numpy as np

# Define 1D Voronoi discretization for Earth structure
velocity_prior = bb.prior.UniformPrior(
    name="vs",
    vmin=2.5,
    vmax=5.0,
    perturb_std=0.15
)

voronoi = bb.discretization.Voronoi1D(
    name="earth_structure",
    vmin=0,           # top of model (surface)
    vmax=100,         # bottom of model (100 km depth)
    perturb_std=5.0,  # site position perturbation
    n_dimensions=None,
    n_dimensions_min=2,
    n_dimensions_max=30,
    parameters=[velocity_prior],
    birth_from="neighbour"  # or "prior"
)

# Use in parameterization
parameterization = bb.parameterization.Parameterization(voronoi)
state = parameterization.initialize()

# Access discretization and parameter values
depths = state["earth_structure"]["discretization"]
velocities = state["earth_structure"]["vs"]
print(f"Number of cells: {len(depths)}")
print(f"Cell positions: {depths}")
print(f"Velocities: {velocities}")

# After inversion, analyze results
# results = inversion.get_results()
# samples_sites = results["earth_structure.discretization"]
# samples_vs = results["earth_structure.vs"]

# Compute statistics on interpolated profiles
interp_depths = np.linspace(0, 100, 200)
# stats = bb.discretization.Voronoi1D.get_tessellation_statistics(
#     samples_sites, samples_vs, interp_depths
# )
# print(f"Mean velocity profile shape: {stats['mean'].shape}")

# Plotting utilities
# bb.discretization.Voronoi1D.plot_tessellation_density(
#     samples_sites, samples_vs, position_bins=100
# )
```

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

### Run Bayesian Inversion

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/51_sw_rf_joint.ipynb

Configure and execute the inversion process using a specified sampler and iteration parameters.

```python
inversion = bb.BayesianInversion(
    parameterization=parameterization, 
    log_likelihood=log_likelihood,
    n_chains=24
)
inversion.run(
    sampler=bb.samplers.SimulatedAnnealing(temperature_start=3), 
    n_iterations=600_000, 
    burnin_iterations=150_000, 
    save_every=150,
    verbose=False
)
```

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

### Configure Data Target

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/11_gaussian_mixture.ipynb

Creates a Target object to handle observed data and noise parameters, including bounds for standard deviation.

```python
target = bb.Target('my_data', 
                   data_obs, 
                   std_min=0, 
                   std_max=0.01, 
                   std_perturb_std=0.001,
                   noise_is_correlated=False)
```

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

### Execute Bayesian inversion

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/02_hierarchical_polyfit.ipynb

Initializes and runs the Bayesian inversion process with specified chain and iteration parameters.

```python
inversion = bb.BayesianInversion(
    log_likelihood=log_likelihood, 
    parameterization=parameterization, 
    n_chains=10
)
inversion.run(
    n_iterations=50_000, 
    burnin_iterations=20_000, 
    save_every=100, 
    verbose=False, 
)
```

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

### Define prior probability distributions

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/02_hierarchical_polyfit.ipynb

Sets up Gaussian and Uniform priors for the model parameters.

```python
m0 = bb.prior.GaussianPrior(name="m0", mean=20, std=1, perturb_std=0.5)
m1 = bb.prior.UniformPrior(name="m1", vmin=-13, vmax=-7, perturb_std=0.4)
m2 = bb.prior.UniformPrior(name="m2", vmin=-10, vmax=4, perturb_std=0.5)
m3 = bb.prior.GaussianPrior(name="m3", mean=1, std=0.1, perturb_std=0.1)
```

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

### Configure Parallel Jobs for Sampling

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Demonstrates how to manually specify the number of parallel jobs using the parallel_config argument in the run method.

```python
inversion.run(
    n_iterations=50_000, 
    burnin_iterations=5_000, 
    save_every=50, 
    verbose=False, 
    parallel_config=dict(n_jobs=5)
)
```

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

### Run Bayesian Inversion with BayesianInversion

Source: https://context7.com/fmagrini/bayes-bay/llms.txt

Configures and executes an MCMC inversion process using defined parameterizations and likelihood functions.

```python
import bayesbay as bb
import numpy as np

# Define observed data and noise
m_true = 10
sigma = 0.5
d_obs = np.random.normal(m_true, sigma)

# Define prior for the parameter
m_prior = bb.prior.UniformPrior(name="m", vmin=5, vmax=15, perturb_std=0.4)

# Create parameter space with fixed dimensionality
param_space = bb.parameterization.ParameterSpace(
    name="model",
    n_dimensions=1,
    parameters=[m_prior]
)
parameterization = bb.parameterization.Parameterization(param_space)

# Define forward function
def fwd_function(state: bb.State) -> np.ndarray:
    return state["model"]["m"]

# Create target with known noise
target = bb.likelihood.Target(
    name="data",
    dobs=d_obs,
    covariance_mat_inv=1/sigma**2
)

# Setup log likelihood
log_likelihood = bb.likelihood.LogLikelihood(
    targets=target,
    fwd_functions=fwd_function
)

# Run Bayesian inversion with 10 chains
inversion = bb.BayesianInversion(
    parameterization=parameterization,
    log_likelihood=log_likelihood,
    n_chains=10
)

inversion.run(
    n_iterations=50000,
    burnin_iterations=5000,
    save_every=50,
    verbose=True,
    print_every=1000
)

# Get results
results = inversion.get_results(concatenate_chains=True)
m_samples = np.array(results["model.m"])
print(f"Mean estimate: {np.mean(m_samples):.3f}")
print(f"Std estimate: {np.std(m_samples):.3f}")
```

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

### Define ParameterSpace and Parameterization

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/11_gaussian_mixture.ipynb

Initializes a parameter space with specific dimensions and parameters, then aggregates it into a parameterization object.

```python
param_space = bb.parameterization.ParameterSpace(
    name="my_param_space", 
    n_dimensions=3, 
    parameters=[mean, std, weight], 
)
parameterization = bb.parameterization.Parameterization(param_space)
```

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

### Configure Noise Standard Deviation Parameters

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/51_sw_rf_joint.ipynb

Sets up the grid layout and defines standard deviation bounds for surface wave and receiver function noise.

```python
gs_bottom = gridspec.GridSpecFromSubplotSpec(
    2, 2, subplot_spec=gs_main[1], wspace=0.05, hspace=0.05
)

std_min_sw, std_max_sw = 0.001, 0.1
std_min_rf, std_max_rf = 0.001, 0.1
```

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

### Verify Prior Distribution

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Prints the bounds and probability density values for the defined prior to verify its configuration.

```python
print(f'a = {m_prior.vmin}')
print(f'b = {m_prior.vmax}')
print(f'p(m=5) = {round(np.exp(m_prior.log_prior(5)), 10)}')
print(f'p(m=10) = {round(np.exp(m_prior.log_prior(10)), 10)}')
print(f'p(m=15) = {round(np.exp(m_prior.log_prior(10)), 10)}')
print(f'p(m=20) = {round(np.exp(m_prior.log_prior(20)), 10)}')
```

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

### Define forward model for surface waves

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/51_sw_rf_joint.ipynb

Creates a partial function for forward modeling of surface-wave dispersion curves (Rayleigh and Love waves) using the PhaseDispersion class from the disba library. It takes the current state of the Voronoi discretization as input.

```python
def forward_sw(state, wave='rayleigh', mode=0):
    voronoi = state["voronoi"]
    voronoi_sites = voronoi["discretization"]
    thickness = Voronoi1D.compute_cell_extents(voronoi_sites)
    vs = voronoi["vs"]
    vp = vs * VP_VS
    rho = 0.32 * vp + 0.77
    pd = PhaseDispersion(thickness, vp, vs, rho)
    d_pred = pd(PERIODS, mode=mode, wave=wave).velocity
    return d_pred

forward_rayleigh = partial(forward_sw, wave='rayleigh', mode=0)
forward_love = partial(forward_sw, wave='love', mode=0)
```

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

### Create venv Virtual Environment

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/installation.md

Use this command to create a new virtual environment with venv. Replace `<path-to-new-env>` and `<env-name>` with your desired paths and environment name.

```bash
python -m venv <path-to-new-env>/bb_env
```

```bash
python -m venv ~/my_envs/bb_env
```

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

### Activate venv Virtual Environment

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/installation.md

Activate the virtual environment created with venv. This command needs to be run in your shell.

```bash
source <path-to-new-env>/<env-name>/bin/activate
```

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

### Save and Load Data to State Cache

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Demonstrates how to save arbitrary data to the State's cache for later retrieval and how to load it back. Useful for reusing computationally intensive results across MCMC iterations.

```python
# Save
my_result = np.random.randn() # Any data type works
state.save_to_cache(name='my_result_name', value=my_result)

# Load
my_result = state.load_from_cache(name='my_result_name')
```

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

### Initialize Shear Wave Velocity Prior

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/21_rayleigh.ipynb

Defines a uniform prior for shear wave velocity (vs) with specified bounds and positions. Includes a custom initialization function to sort random uniform values within the parameter bounds.

```python
def initialize_vs(param, positions=None):
    vmin, vmax = param.get_vmin_vmax(positions)
    sorted_vals = np.sort(np.random.uniform(vmin, vmax, positions.size))
    return sorted_vals

vs = bb.prior.UniformPrior(name="vs", 
                                    vmin=[2.2, 2.8, 3.3, 4], 
                                    vmax=[3.9, 4.6, 4.8, 5], 
                                    position=[0, 20, 60, 150],
                                    perturb_std=0.15)
vs.set_custom_initialize(initialize_vs)
```

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

### Define Parameter Space and Parameterization

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/03_transd_polyfit.ipynb

Configures the parameter space with a range of allowed dimensions and initializes the parameterization object.

```python
# Define parameterization
param_space = bb.parameterization.ParameterSpace(
    name="my_param_space", 
    n_dimensions_min=N_DIMS_MIN, 
    n_dimensions_max=N_DIMS_MAX, 
    parameters=[coefficients]
)
parameterization = bb.parameterization.Parameterization(param_space)
```

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

### Create conda/mamba Virtual Environment

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/installation.md

Create a new virtual environment using conda or mamba, specifying the Python version. Replace `<env-name>` with your desired environment name.

```bash
conda create -n <env-name> python=3.10
```

```bash
conda create -n bb_env python=3.10
```

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

### Creating Partial Forward Functions

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/22_rayleigh_love.ipynb

This uses functools.partial to create specialized forward functions for Rayleigh and Love waves, fixing the 'wave' and 'mode' arguments of the general forward_sw function.

```python
from functools import partial

forward_rayleigh = partial(forward_sw, wave='rayleigh', mode=0)
forward_love = partial(forward_sw, wave='love', mode=0)
```

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

### Visualizing Sampled Voronoi Tessellations

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/31_sw_tomography.ipynb

Displays multiple sampled Voronoi tessellations from the inversion results. Each panel shows a different realization of the tessellation and the corresponding velocity distribution, allowing for assessment of model uncertainty.

```python
fig, axes = plt.subplots(2, 3, figsize=(10, 6))
random_indexes = np.random.choice(range(len(results['voronoi.vel'])), size=9, replace=False)
for ipanel, (ax, irandom) in enumerate(zip(axes.ravel(), random_indexes)):
    voronoi_sites = results['voronoi.discretization'][irandom]
    velocity = results['voronoi.vel'][irandom]
    ax, cbar = Voronoi2D.plot_tessellation(
        voronoi_sites, 
        velocity, 
        ax=ax, 
        cmap=scm.roma,
        voronoi_sites_kwargs=dict(markersize=0)
    )
    ax.tick_params(labelleft=False, labelbottom=False)
    ax.set_xlabel('')
    ax.set_ylabel('')
    cbar.set_label('Velocity [km/s]')
    if ipanel in [0, 3]:
        ax.tick_params(labelleft=True)
        ax.set_ylabel('y')
    if ipanel not in [2, 5]:
        cbar.set_ticklabels('')
        cbar.set_label('')
    if ipanel>2:
        ax.set_xlabel('x')
        ax.tick_params(labelbottom=True)
plt.tight_layout()
plt.show()
```

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

### Initialize LogLikelihood Instance

Source: https://github.com/fmagrini/bayes-bay/blob/main/docs/source/tutorials/00_quickstart.ipynb

Initializes a LogLikelihood instance using the defined Target and forward function(s). This object is central to the Bayesian inference process, quantifying the probability of the data given the model.

```python
log_likelihood = bb.likelihood.LogLikelihood(targets=target, fwd_functions=fwd_function)
print(log_likelihood)
```