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, interpolation_positions=None)

Utility class for Voronoi tessellation in 2D

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, shapely.geometry.MultiPolygon], 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 (int, 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 2 and 100. 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 2 and 100. 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 = 2, n_dimensions_max = 100, 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_min)
```
parameters (List[Prior], 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 nearest-neighbour lookup at every perturbation of the discretization, stored in each state’s cache under the key "kdtree". Use this when the forward function needs distances, multiple nearest neighbours, or query points that vary between iterations; when it only interpolates the tessellation onto a fixed set of points, prefer interpolation_positions, which is much faster
interpolation_positions (np.ndarray of shape (m, 2), optional) – fixed positions onto which the tessellation is interpolated during the sampling; when given, every state carries in its cache the index of the Voronoi cell each position falls in, kept up to date through exact incremental updates, so that forward functions can interpolate the tessellation through get_interpolated_values() (see also get_nearest_site_indices() and set_interpolation_positions()). Configure these positions before starting the first sampling run

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

Notes

Subclasses with polygonal domains draw the newborn site by rejection sampling from the position prior restricted to the polygon. Redrawing is valid here because the birth-position proposal is state independent and exactly equals that restricted prior; unlike a move proposal centred on the current site, it introduces no state-dependent normalization.

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_interpolated_values(ps_state, param)

interpolates the values that the given free parameter takes in each Voronoi cell onto the positions registered through set_interpolation_positions() (nearest-neighbour interpolation, i.e., piecewise constant within each Voronoi cell), given the discretization found in the passed parameter space state

This is equivalent to (but, during the sampling, much faster than) querying the registered positions against a nearest-neighbour structure built on the Voronoi sites:

indices = self.get_nearest_site_indices(ps_state)
interp_values = values[indices]

Parameters:

ps_state (ParameterSpaceState) – the parameter space state holding the discretization
param (Union[str, np.ndarray]) – either the name of a free parameter associated with this discretization or an array of values, one per Voronoi cell

Returns:

the parameter values interpolated onto the registered positions

Return type:

np.ndarray of shape (m,)

get_kdtree(ps_state)

Return the state’s site KD-tree, building it on demand when needed.

The tree uses the coordinates appropriate to the tessellation geometry: planar coordinates for Voronoi2D and 3-D unit vectors for Voronoi2DSphere. Spherical longitude-latitude query positions must therefore be converted with Voronoi2DSphere.lonlat_to_xyz() before calling query. The accessor supports initialized, sampled, nested-birth, and user-created parameter-space states.

get_nearest_site_indices(ps_state)

returns, for each of the positions registered through set_interpolation_positions(), the index of the Voronoi cell it belongs to (i.e., of the nearest Voronoi site), given the discretization found in the passed parameter space state

The indices are read from the cache of the given state, where they are maintained through exact incremental updates at every perturbation of the discretization. If they are not present in the cache (e.g., for a state that was not generated by this class, such as a user-provided starting state), they are computed from scratch and stored in it.

Parameters:: ps_state (ParameterSpaceState) – the parameter space state holding the discretization
Returns:: the index of the Voronoi cell each registered position belongs to
Return type:: np.ndarray of shape (m,)

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) – a list of arrays containing the Voronoi-site positions
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 “percentiles”

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]

initialize_interpolation(ps_state)

computes from scratch the nearest-site assignment of the positions registered through set_interpolation_positions() and stores it in the cache of the given parameter space state

This is done automatically at the initialization of each Markov chain and, incrementally, at every perturbation; calling this method is only needed for states created by other means (e.g. a user-provided starting state).

Parameters:: ps_state (ParameterSpaceState) – the parameter space state whose cache is to be filled
Returns:: the same instance, with "interp_nearest" and "interp_affinity" stored in its cache
Return type:: 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
interp_positions ((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, clip_polygon=None, **kwargs)

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, optional) – when given, the Voronoi nuclei are displayed, styled by these keyword arguments (passed to matplotlib.pyplot.plot). By default (None), the nuclei are not displayed
voronoi_plot_2d_kwargs (dict) – keyword arguments passed to scipy.spatial.voronoi_plot_2d, used to plot the Voronoi interfaces. When clip_polygon is given, only its entries line_colors and line_width are used, to style the edges of the clipped cells
clip_polygon (Union[np.ndarray, shapely.geometry.Polygon, shapely.geometry.MultiPolygon], optional) – region of interest to which the displayed tessellation is clipped. When given, each Voronoi cell is intersected with this polygon before being drawn
kwargs (dict, optional) – when clip_polygon is given, additional keyword arguments passed to matplotlib.axes.Axes.fill (e.g., transform, alpha, zorder); ignored otherwise

Returns:

ax (matplotlib.axes.Axes) – The axes containing the tessellation.
cbar (Union[matplotlib.colorbar.Colorbar, None]) – The colourbar, or None when param_values is omitted.

Notes

Voronoi2D cell edges are straight segments in Cartesian coordinates, so they are drawn exactly and require no boundary-spacing setting. Spherical Voronoi edges are curved and therefore the plotting methods of Voronoi2DSphere sample those curves internally at a fixed angular spacing suitable for plotting.

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()

Draw a site uniformly from the rectangular or polygonal position prior.

For a polygon, rejection sampling redraws from the state-independent bounding-box proposal until the site lies inside. The result is exactly the uniform prior restricted to the polygon, so it is suitable as the birth proposal without a Hastings correction.

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)

set_interpolation_positions(positions)

registers a fixed set of positions onto which the tessellation is interpolated during the sampling

Once registered, every parameter space state carries in its cache:

"interp_nearest": the index of the Voronoi site each registered position belongs to, i.e. the nearest site
"interp_affinity": the affinity between each registered position and its Voronoi site, i.e. the negative squared Euclidean distance (Voronoi2D) or the cosine of the angular distance (Voronoi2DSphere)

Both entries are kept up to date at every perturbation of the discretization through exact incremental updates (a perturbation only affects the assignments of the positions within the perturbed Voronoi cells), which is much faster than a nearest-neighbour query of all the registered positions at every Markov chain iteration. Registering a new set of positions increments an internal version; existing state caches are then recomputed automatically on first use.

Configure the interpolation positions before starting the first sampling run. Process-based workers receive a serialized copy of the discretization configuration, so later changes to the original object do not propagate to workers that are already running.

Parameters:: positions (np.ndarray of shape (m, 2)) – the interpolation positions

back to top