bayesbay.discretization.Voronoi2D

class bayesbay.discretization.Voronoi2D(name, vmin=None, vmax=None, polygon=None, perturb_std=1, n_dimensions=None, n_dimensions_min=2, n_dimensions_max=100, n_dimensions_init_range=0.3, parameters=None, birth_from='neighbour', compute_kdtree=False)

Utility class for Voronoi tessellation in 1D

Parameters:
  • name (str) – name attributed to the Voronoi tessellation, for display and storing purposes

  • vmin (Union[Number, np.ndarray]) – minimum/maximum value bounding each dimension. Ignored when polygon is not None

  • vmax (Union[Number, np.ndarray]) – minimum/maximum value bounding each dimension. Ignored when polygon is not None

  • polygon (Union[np.ndarray, shapely.geometry.Polygon], optional) – polygon defining the domain of the Voronoi tessellation; Voronoi sites outside this polygon are not allowed

  • perturb_std (Union[Number, np.ndarray]) – standard deviation of the Gaussians used to randomly perturb the Voronoi sites in each dimension.

  • n_dimensions (Number, optional) – number of dimensions. None (default) results in a trans-dimensional discretization, with the dimensionality of the parameter space allowed to vary in the range n_dimensions_min-n_dimensions_max

  • n_dimensions_min (Number, optional) – minimum and maximum number of dimensions, by default 1 and 10. These parameters are ignored if n_dimensions is not None, i.e. if the discretization is not trans-dimensional

  • n_dimensions_max (Number, optional) – minimum and maximum number of dimensions, by default 1 and 10. These parameters are ignored if n_dimensions is not None, i.e. if the discretization is not trans-dimensional

  • n_dimensions_init_range (Number, optional) –

    percentage of the range n_dimensions_min - n_dimensions_max used to initialize the number of dimensions (0.3. by default). For example, if n_dimensions_min = 1, n_dimensions_max = 10, and n_dimensions_init_range = 0.5, the maximum number of dimensions at the initialization is:

    int((n_dimensions_max - n_dimensions_min) * n_dimensions_init_range + n_dimensions_max)
    

  • parameters (List[Parameter], optional) – a list of free parameters, by default None

  • birth_from ({"prior", "neighbour"}, optional) – whether to initialize the free parameters associated with the newborn Voronoi cell by randomly drawing from their prior or by perturbing the value found in the nearest Voronoi cell (default)

  • compute_kdtree (bool) – whether to compute a kd-tree for quick nearest-neighbor lookup at every perturbation of the discretization

Reference Details

is_leaf
name
parameters

all the free parameters defined in this parameter space

perturbation_funcs

the list of perturbation functions allowed in the parameter space linked to the Voronoi discretization. Each function takes in a state (see State) and returns a new state along with the corresponding partial acceptance probability,

\[\begin{split}\underbrace{\alpha_{p}}_{\begin{array}{c} \text{Partial} \\ \text{acceptance} \\ \text{probability} \end{array}} = \underbrace{\frac{p\left({\bf m'}\right)}{p\left({\bf m}\right)}}_{\text{Prior ratio}} \underbrace{\frac{q\left({\bf m} \mid {\bf m'}\right)}{q\left({\bf m'} \mid {\bf m}\right)}}_{\text{Proposal ratio}} \underbrace{\lvert \mathbf{J} \rvert}_{\begin{array}{c} \text{Jacobian} \\ \text{determinant} \end{array}},\end{split}\]
perturbation_weights

a list of perturbation weights, corresponding to each of the perturbation_funcs() that determines the probability of each of them to be chosen during each step

The weights are not normalized and have the following default values:

  • Birth/Death perturbations: 1

  • Parameter values perturbation: 3

  • Voronoi site perturbation: 1

trans_d

indicates whether the current configuration allows changes in dimensionality

add_hyper_params(hyper_params)

Sets the attributes from the given dict and checks for errors

Parameters:

hyper_params (Dict[str, Union[Number, np.ndarray]]) – dictionary of attributes to be set

birth(old_ps_state)

creates a new Voronoi cell, initializes all free parameters associated with it, and returns the pertubed state along with the log of the corresponding partial acceptance probability,

\[\begin{split}\underbrace{\alpha_{p}}_{\begin{array}{c} \text{Partial} \\ \text{acceptance} \\ \text{probability} \end{array}} = \underbrace{\frac{p\left({\bf m'}\right)}{p\left({\bf m}\right)}}_{\text{Prior ratio}} \underbrace{\frac{q\left({\bf m} \mid {\bf m'}\right)}{q\left({\bf m'} \mid {\bf m}\right)}}_{\text{Proposal ratio}} \underbrace{\lvert \mathbf{J} \rvert}_{\begin{array}{c} \text{Jacobian} \\ \text{determinant} \end{array}}.\end{split}\]

In this case, the prior probability of the model \({\bf m}\) is

\[p({\bf m}) = p({\bf c} \mid k) p(k) \prod_i{p({\bf v}_i \mid {\bf c})} ,\]

where \(k\) denotes the number of Voronoi cells, each entry of the vector \({\bf c}\) corresponds to the position of a Voronoi site, and each \(i\)th free parameter \({\bf v}\) has the same dimensionality as \({\bf c}\).

Following [1], \(p({\bf c} \mid k) = \frac{k! \left(N - k \right)!}{N!}\). If we then assume that \(p(k) = \frac{1}{\Delta k}\), where \(\Delta k = k_{max} - k_{min}\), the prior ratio reads

\[\frac{p({\bf m'})}{p({\bf m})} = \frac{(k+1) \prod_i p(v_i^{k+1})}{(N-k)},\]

where \(p(v_i^{k+1})\) denotes the prior probability of the newly born \(i\)th parameter, which may be dependent on \({\bf c}\). The proposal ratio reads

\[\frac{q({\bf m} \mid {\bf m'})}{q({\bf m'} \mid {\bf m})} = \frac{(N-k)}{(k+1) \prod_i q_{v_i}^{k+1}},\]

where \(q_{v_i}^{k+1}\) denotes the proposal probability for the newly born \(i\)th parameter in the new dimension. It is easy to show that, in the case of a birth from neighbor [1] or a birth from prior [2] (see birth_from), \(\lvert \mathbf{J} \rvert = 1\) and \(\alpha_{p} = \frac{p({\bf m'})}{p({\bf m})} \frac{q({\bf m} \mid {\bf m'})}{q({\bf m'} \mid {\bf m})}\). It follows that

\[\alpha_{p} = \frac{(k+1) \prod_i p(v_i^{k+1})}{(N-k)} \frac{(N-k)}{(k+1) \prod_i q_{v_i}^{k+1}} = \frac{\prod_i p(v_i^{k+1})}{\prod_i{q_{v_i}^{k+1}}}.\]

In the case of a birth from prior, \(q_{v_i}^{k+1} = p(v_i^{k+1})\) and

\[\alpha_{p} = \frac{\prod_i p(v_i^{k+1})}{\prod_i{p(v_i^{k+1})}} = 1.\]

In the case of a birth from neighbor, \(q_{v_i}^{k+1} = \frac{1}{\theta \sqrt{2 \pi}} \exp \lbrace -\frac{\left( v_i^{k+1} - v_i \right)^2}{2\theta^2} \rbrace\), where the newly born value, \(v_i^{k+1}\), is generated by perturbing the original value, \(v_i\), of the \(i\)th parameter. This is achieved through a random deviate from the normal distribution \(\mathcal{N}(v_i, \theta)\), with \(\theta\) denoting the standard deviation of the Gaussian used to carry out the perturbation (see, for example, bayesbay.prior.UniformPrior.perturb_std) . The partial acceptance probability is then computed numerically.

Parameters:

old_ps_state (ParameterSpaceState) – current parameter space state

Returns:

  • ParameterSpaceState – new parameter space state

  • Number – log of the partial acceptance probability, \(log(\alpha_{p}) = \log(\frac{\prod_i p(v_i^{k+1})}{\prod_i{q_{v_i}^{k+1}}})\)

References

death(old_ps_state)

removes a new Voronoi cell and returns the pertubed state along with the log of the corresponding partial acceptance probability,

\[\begin{split}\underbrace{\alpha_{p}}_{\begin{array}{c} \text{Partial} \\ \text{acceptance} \\ \text{probability} \end{array}} = \underbrace{\frac{p\left({\bf m'}\right)}{p\left({\bf m}\right)}}_{\text{Prior ratio}} \underbrace{\frac{q\left({\bf m} \mid {\bf m'}\right)}{q\left({\bf m'} \mid {\bf m}\right)}}_{\text{Proposal ratio}} \underbrace{\lvert \mathbf{J} \rvert}_{\begin{array}{c} \text{Jacobian} \\ \text{determinant} \end{array}}.\end{split}\]

It is straightforward to show that this equals the reciprocal of the partial acceptance probability obtained in the case of a birth perturbation (see birth()), i.e.,

\[\alpha_{p} = \frac{\prod_i{q_{v_i}^{k+1}}}{\prod_i p(v_i^{k+1})}.\]
Parameters:

old_ps_state (ParameterSpaceState) – current parameter space state

Returns:

  • ParameterSpaceState – new parameter space state

  • Number – log of the partial acceptance probability, \(log(\alpha_{p}) = -\log(\frac{\prod_i p(v_i^{k+1})}{\prod_i{q_{v_i}^{k+1}}})\)

get_hyper_param(hyper_param, position=None)

Retrieves the value corresponding to the specified attribute, which may be a function of position

Parameters:
  • hyper_param (str) – the name of the attribute

  • position (Union[np.ndarray, Number], optional) – the position (in the discretization domain) associated with the value, None by default

Returns:

value corresponding to the specified attribute

Return type:

Union[Number, np.ndarray]

get_perturb_std(*args)

get the standard deviation of the Gaussian used to perturb the discretization

get_perturb_std_birth(position=None)

get the standard deviation of the Gaussian used to perturb the parameter at birth, which may be dependent on the position in the discretization domain

Parameters:

position (Union[Number, np.ndarray], optional) – the position in the discretization domain at which the standard deviation of the Gaussian used to perturb the parameter at birth is returned. Default is None

Returns:

standard deviation of the Gaussian used to perturb the parameter at birth, possibly at the specified position

Return type:

Number

static get_tessellation_statistics(samples_voronoi_cells, samples_param_values, interp_positions, percentiles=(10, 90))

get the mean, median, std and percentiles of the given ensemble

Parameters:
  • samples_voronoi_cells (list) – either a list of Voronoi-cell extents or of Voronoi-site positions (see input_type)

  • samples_param_values (list) – a list of parameter values to draw statistics from

  • interp_positions (np.ndarray) – points to interpolate

  • percentiles (tuple, optional) – percentiles to calculate, by default (10, 90)

Returns:

a dictionary with these keys: “mean”, “median”, “std” and “percentile”

Return type:

dict

has_hyper_param(hyper_param)

Whether or not the Prior instance has the specified attribute

Parameters:

hyper_param (str) – the name of the attribute

initialize(position=None)

initializes the parameter space linked to the Voronoi tessellation

Returns:

an initial parameter space state, or a list of parameter space states

Return type:

Union[ParameterSpaceState, List[ParameterSpaceState]

static interpolate_tessellation(voronoi_sites, param_values, interp_positions)

nearest neighbour interpolation based on Voronoi-site positions and values associated with them.

Parameters:
  • voronoi_sites ((n, 2) np.ndarray) – the positions of the Voronoi sites

  • param_values ((n,) np.ndarray) – the parameter values associated with each Voronoi cell

  • query_points ((m, 2) np.ndarray) – the positions where interpolation is performed

Returns:

interpolated values

Return type:

np.ndarray

log_prior(*args)

BayesBay implements the grid trick, which calculates the prior probability of a Voronoi discretization through the combinatorial formula \({N \choose k}^{-1}\), with k denoting the number of Voronoi sites and N the number of possible positions allowed for the sites [3].

References

log_prob_initialize_discretization(ps_state)

The log of the partial acceptance probability of the birth of the discretization. This includes only the discretization but not the parameter values.

\(\frac{p(k')}{p(k)} \frac{p(c'|k')}{p(c|k)} \frac{q(c'|m)}{q(c|m')}\)

Parameters:

ps_state (ParameterSpaceState) – the newly-born parameter space state

Returns:

the log of the partial acceptance probability of the birth of the discretization

Return type:

Number

nearest_neighbour(discretization, query_point)

returns the index of the nearest neighbour of a given query point in the discretization

Parameters:
  • discretization (np.ndarray) – the discretization

  • query_point (Union[Number, np.ndarray]) – the query point

Returns:

the index of the nearest neighbour point

Return type:

int

perturb_value(old_ps_state, isite)

perturbs the value of one Voronoi site and calculates the log of the partial acceptance probability

\[\begin{split}\underbrace{\alpha_{p}}_{\begin{array}{c} \text{Partial} \\ \text{acceptance} \\ \text{probability} \end{array}} = \underbrace{\frac{p\left({\bf m'}\right)}{p\left({\bf m}\right)}}_{\text{Prior ratio}} \underbrace{\frac{q\left({\bf m} \mid {\bf m'}\right)}{q\left({\bf m'} \mid {\bf m}\right)}}_{\text{Proposal ratio}} \underbrace{\lvert \mathbf{J} \rvert}_{\begin{array}{c} \text{Jacobian} \\ \text{determinant} \end{array}}.\end{split}\]
Parameters:
  • old_ps_state (ParameterSpaceState) – the current parameter space state

  • isite (Number) – the index of the Voronoi site to be perturbed

Returns:

the new parameter space state and its associated partial acceptance probability excluding log likelihood ratio

Return type:

Tuple[ParameterSpaceState, Number]

static plot_tessellation(voronoi_sites, param_values=None, ax=None, cmap='viridis', norm=None, vmin=None, vmax=None, voronoi_sites_kwargs=None, voronoi_plot_2d_kwargs=None)

display the Voronoi tessellation

Parameters:
  • voronoi_sites (np.ndarray of shape (m, 2)) – 2D Voronoi-site positions

  • param_values (np.ndarray, optional) – parameter values associated with each Voronoi cell. These could represent the physical property inferred in each cell of the discretized medium

  • ax (matplotlib.axes.Axes, optional) – an optional Axes object to plot on

  • cmap (Union[str, matplotlib.colors.Colormap]) – the Colormap instance or registered colormap name used to map scalar data to colors

  • norm (Union[str, matplotlib.colors.Normalize]) – the normalization method used to scale scalar data to the [0, 1] range before mapping to colors using cmap. By default, a linear scaling is used, mapping the lowest value to 0 and the highest to 1.

  • vmin (Number) – minimum and maximum values used to create the colormap

  • vmax (Number) – minimum and maximum values used to create the colormap

  • voronoi_sites_kwargs (dict) – keyword arguments passed to matplotlib.pyplot.plot, used to plot the voronoi nuclei

  • voronoi_plot_2d_kwargs (dict) – keyword arguments passed to scipy.spatial.voronoi_plot_2d, used to plot the Voronoi interfaces

Returns:

ax – The Axes object containing the 2D histogram

Return type:

matplotlib.axes.Axes

sample(*args)

sample a random ParameterSpaceState instance, including the number of dimensions, the Voronoi sites, and the parameter values

sample_discretization()

sample a parameter space state that only contains the discretization

sample_site()

draws a Voronoi-site position at random within the discretization domain

set_custom_initialize(initialize_func)

sets a custom initialization function

Parameters:

initialize_func (Callable[["Prior", np.ndarray], np.ndarray]) – The function to use for initialization. This function should take a Prior instance and optionally an array of positions as input arguments, and produce an array of values as output.

Examples

def my_init(
    param: bb.prior.Prior,
    position: np.ndarray
) -> np.ndarray:
    print("This is my custom init!")
    return np.ones(len(position))

my_param.set_custom_initialize(my_init)

back to top