surrDAMH 1.0

Python implementation of surrogate-accelerated Markov chain Monte Carlo methods for the numerical realization of the Bayesian inversion

The library is based on the delayed acceptance Metropolis-Hastings algorithm and accelerated using surrogate models and using parallelization. It provides samples from the posterior distribution π(u|y) ∝ f_η(y - G(u)) π₀(u), where y is a given vector of observations, G is an observation operator, f_η is a probability density function (pdf) of Gaussian observational noise, π₀(u) is Gaussian prior pdf.

Requirements

• numpy
• scipy
• pandas
• json
• mpi4py
• petsc4py (for "Darcy" example)
• MyFEM (for "Darcy" example)
  • github: https://github.com/Beremi/Simple_Python_PETSc_FEM/tree/SD
• pcdeflation (for the use of an own deflation basis in "Darcy" example)
  • make -C examples/solvers/pcdeflation clean
  • make -C examples/solvers/pcdeflation build
• cython (for pcdeflation build)

Instructions for running the sampling process

• problem_name
  • defines the name of the json configuration file that will be loaded: "conf/" + problem_name + ".json"
  • prepared toy examples: "simple", "simple_MPI", "Darcy" (in the form of json configuration files)
• N = number of sampling processes

run the sampling process:

python3 run.py problem_name N (oversubscribe)

• toy examples:

  • python3 run.py simple 4
  • python3 run.py simple_MPI 4
  • python3 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/problem_name.

run the visualization of obtained samples:

python3 run.py problem_name N visualize

• toy examples:

  • python3 run.py simple 4 visualize
  • python3 run.py simple_MPI 4 visualize
  • python3 run.py Darcy 4 visualize

• python3 run.py simple 4 visualize   <=>   python3 examples/visualization/simple.py 4

MPI processes

• 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)

spawned solvers:

• 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]) -> None
• get_observations(self) -> list[no_observations]

JSON configuration file example & comments

{

"problem_name": "Darcy",

"no_parameters": 4,

"no_observations": 3, length of the vector of the 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

},

"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": "polyharmonic",

"beta": 3

},

"surr_updater_parameters": { optional - arguments of class Surrogate_update that updates surrogate model

"kernel_type": "polyharmonic",

"beta": 3,

"max_centers": 500,

"expensive": false

}

}