Python implementation of surrogate-accelerated Markov chain Monte Carlo methods for Bayesian inversion
Provides samples from the posterior distribution π(u|y) ∝ fη(y - G(u)) π0(u), where y is a given vector of observations, G is an observation operator, fη is probability density function (pdf) of Gaussian observational noise, π0(u) is Gaussian prior pdf.
- numpy
- scipy
- pandas
- json
- mpi4py
- petsc4py (for "Darcy" example)
- MyFEM (for "Darcy" example)
- pcdeflation (for the use of an own deflation basis in "Darcy" example)
make -C examples/solvers/pcdeflation cleanmake -C examples/solvers/pcdeflation build
- cython (for pcdeflation build)
- conf_name
- prepared toy examples: "simple", "simple_MPI", "Darcy"
- loads configuration file "conf/" + conf_name + ".json"
- N = number of sampling processes
-
toy examples:
python3 run.py simple 4python3 run.py simple_MPI 4python3 run.py Darcy 4
-
use "oversubscribe" if there are not enough slots available, for example:
python3 run.py simple 4 oversubscribe
-
python3 run.py simple 4<=>mpirun -n 4 python3 -m mpi4py surrDAMH/process_SAMPLER.py : -n 1 python3 -m mpi4py surrDAMH/process_SOLVER.py simple : -n 1 python3 -m mpi4py surrDAMH/process_COLLECTOR.py -
python3 run.py simple 4 oversubscribe<=>mpirun -n 4 --oversubscribe python3 -m mpi4py surrDAMH/process_SAMPLER.py : -n 1 --oversubscribe python3 -m mpi4py surrDAMH/process_SOLVER.py simple : -n 1 --oversubscribe python3 -m mpi4py surrDAMH/process_COLLECTOR.py
Obtained samples are saved into saved_samples/conf_name.
-
toy examples:
python3 run.py simple 4 visualizepython3 run.py simple_MPI 4 visualizepython3 run.py Darcy 4 visualize
-
python3 run.py simple 4 visualize<=>python3 examples/visualization/simple.py 4
- process_SAMPLER.py: N sampling processes based on the Metropolis-Hastings (MH) and the delayed-acceptance MH algorithm
- process_SOLVER.py: solvers parent process that spawns new MPI processes (number and type specified in the JSON configuration file)
- process_COLLECTOR.py: collects snapshots and uses them for the construction and updates of the surrogate model ("poly" or "rbf" ... polynomial or radial basis functions based surrogate model)
- provide evaluations of an observation operator
observations = G(parameters)
- e.g. wrapper of an external library (see the "Darcy" example)
Path and constructor arguments of the spawed solver (wrapper) class are specified in the JSON configuration file. The following methods are required (executed by all MPI processes in the spawned communicator):
set_parameters(self, parameters: list[no_parameters]) -> Noneget_observations(self) -> list[no_observations]
{
"no_parameters": 4, *length of the vector of parameters
"no_observations": 3, length of the vector of observations, repetitive observations are not supported
"problem_parameters": {
"prior_mean": 0.0, scalar or vector of length no_parameters
"prior_std": 1.0, scalar or vector or covariance matrix
"observations": [9.62638828, -5.90755323, -3.71883564], vector of observations
"noise_std": [0.2, 0.1, 0.1] standard deviation or covariance matrix of observational noise
"saved_samples_name": "Darcy", optional - use "saved_samples/saved_samples_name/" instead of saved_samples/conf_name/
"visualization_filename": "Darcy.py", optional - use a specific visualization script
"visualization_path": "examples/visualization/", optional - folder of the visualization script
},
"paths_to_append": [ optional
"/home/simona/GIT/Simple_Python_PETSc_FEM"
],
"solver_module_path": "examples/solvers/FEM_interfaces.py",
"solver_module_name": "FEM_interfaces",
"solver_init_name": "FEM", spawned solver class
"solver_parameters": { spawned solver class constructor arguments (in the form of a dictionary)
},
"solver_parent_parameters": { optional
"maxprocs": 4 number of MPI processes of each spawned solver
},
"no_solvers": 3, number of separate solvers to be spawned
"samplers_list": [ dictionary of sampling algoritms run by all sampling processes in consecutive order
{
"type": "MH", MH (does not use surrogate model) or DAMH (uses surrogate model)
"max_samples": 200, maximal length of the chain
"time_limit": 60, maximal sampling time in seconds
"proposal_std": 0.2, standard deviation or covariance matrix of Gaussian proposal distribution
"surrogate_is_updated": true update surrogate while sampling
},
{
"type": "DAMH",
"max_samples": 2000,
"time_limit": 60,
"proposal_std": 0.2,
"surrogate_is_updated": true
},
{
"type": "DAMH",
"max_samples": 10000,
"time_limit": 60,
"proposal_std": 0.2,
"surrogate_is_updated": false
}
],
"surrogate_type": "rbf", "poly" or "rbf" (polynomial or radial basis functions based surrogate model)
"surr_solver_parameters": { optional - arguments of class Surrogate_apply that evaluates surrogate model
"kernel_type": 1
},
"surr_updater_parameters": { optional - arguments of class Surrogate_update that updates surrogate model
"kernel_type": 1,
"no_keep": 500,
"expensive": false
}
}