beetroots.sampler package

Subpackages

• beetroots.sampler.saver package
  • Submodules
  • beetroots.sampler.saver.abstract_saver module
    • Saver
      • Saver.D
      • Saver.D_sampling
      • Saver.L
      • Saver.N
      • Saver.batch_size
      • Saver.check_need_to_save()
      • Saver.check_need_to_update_memory()
      • Saver.final_next_batch_size
      • Saver.freq_save
      • Saver.initialize_memory()
      • Saver.list_idx_sampling
      • Saver.memory
      • Saver.next_batch_size
      • Saver.results_path
      • Saver.save_additional()
      • Saver.save_forward_map_evals
      • Saver.save_to_file()
      • Saver.scaler
      • Saver.set_results_path()
      • Saver.t_last_init
      • Saver.t_last_save
      • Saver.update_memory()
  • beetroots.sampler.saver.my_saver module
    • MySaver
      • MySaver.D
      • MySaver.D_sampling
      • MySaver.L
      • MySaver.N
      • MySaver.batch_size
      • MySaver.final_next_batch_size
      • MySaver.freq_save
      • MySaver.initialize_memory()
      • MySaver.list_idx_sampling
      • MySaver.memory
      • MySaver.next_batch_size
      • MySaver.results_path
      • MySaver.save_forward_map_evals
      • MySaver.scaler
      • MySaver.t_last_init
      • MySaver.t_last_save
      • MySaver.update_memory()
  • Module contents
• beetroots.sampler.utils package
  • Submodules
  • beetroots.sampler.utils.mml module
    • EBayesMMLE
      • EBayesMMLE.N0
      • EBayesMMLE.N1
      • EBayesMMLE.dim
      • EBayesMMLE.exponent
      • EBayesMMLE.homogeneity
      • EBayesMMLE.mean_theta
      • EBayesMMLE.mean_theta_old
      • EBayesMMLE.scale
      • EBayesMMLE.sum_weights
      • EBayesMMLE.update()
      • EBayesMMLE.vmax
      • EBayesMMLE.vmin
    • EBayesMMLELogRate
      • EBayesMMLELogRate.N0
      • EBayesMMLELogRate.N1
      • EBayesMMLELogRate.dim
      • EBayesMMLELogRate.exponent
      • EBayesMMLELogRate.homogeneity
      • EBayesMMLELogRate.mean_theta
      • EBayesMMLELogRate.mean_theta_old
      • EBayesMMLELogRate.scale
      • EBayesMMLELogRate.sum_weights
      • EBayesMMLELogRate.vmax
      • EBayesMMLELogRate.vmin
    • EBayesMMLELogScale
      • EBayesMMLELogScale.N0
      • EBayesMMLELogScale.N1
      • EBayesMMLELogScale.dim
      • EBayesMMLELogScale.exponent
      • EBayesMMLELogScale.homogeneity
      • EBayesMMLELogScale.mean_theta
      • EBayesMMLELogScale.mean_theta_old
      • EBayesMMLELogScale.scale
      • EBayesMMLELogScale.sum_weights
      • EBayesMMLELogScale.vmax
      • EBayesMMLELogScale.vmin
    • EBayesMMLERate
      • EBayesMMLERate.N0
      • EBayesMMLERate.N1
      • EBayesMMLERate.dim
      • EBayesMMLERate.exponent
      • EBayesMMLERate.homogeneity
      • EBayesMMLERate.mean_theta
      • EBayesMMLERate.mean_theta_old
      • EBayesMMLERate.scale
      • EBayesMMLERate.sum_weights
      • EBayesMMLERate.vmax
      • EBayesMMLERate.vmin
    • EBayesMMLEScale
      • EBayesMMLEScale.N0
      • EBayesMMLEScale.N1
      • EBayesMMLEScale.dim
      • EBayesMMLEScale.exponent
      • EBayesMMLEScale.homogeneity
      • EBayesMMLEScale.mean_theta
      • EBayesMMLEScale.mean_theta_old
      • EBayesMMLEScale.scale
      • EBayesMMLEScale.sum_weights
      • EBayesMMLEScale.vmax
      • EBayesMMLEScale.vmin
    • pgd_rate()
    • pgd_rate_log()
    • pgd_scale()
    • pgd_scale_log()
  • beetroots.sampler.utils.my_sampler_params module
    • MySamplerParams
      • MySamplerParams.compute_correction_term
      • MySamplerParams.extreme_grad
      • MySamplerParams.history_weight
      • MySamplerParams.initial_step_size
      • MySamplerParams.is_stochastic
      • MySamplerParams.k_mtm
      • MySamplerParams.selection_probas
  • beetroots.sampler.utils.utils module
    • compute_nlpdf_spatial_proposal()
    • compute_sum_subsets_mean()
    • get_neighboring_pixels()
    • numba_logsumexp_stable()
    • sample_conditional_spatial_and_indicator_prior()
    • sample_conditional_spatial_prior()
    • sample_generalized_gaussian()
    • sample_smooth_indicator()
  • Module contents

Submodules

beetroots.sampler.abstract_sampler module

class beetroots.sampler.abstract_sampler.Sampler(D: int, L: int, N: int, rng: ~numpy.random._generator.Generator = Generator(PCG64) at 0x7FC163622580)

Bases: ABC

abstract class designed to sample from a distribution defined with the Posterior class.

Parameters:

• D (int) – total number of physical parameters to reconstruct

• L (int) – number of observables per component \(n\)

• N (int) – total number of pixels to reconstruct

• rng (np.random.Generator) – numpy random generator (for sampling reproducibility), by default np.random.default_rng(42)

D

total number of physical parameters to reconstruct

Type:

int

ESS_OPTIM = 1000

number of random reproduced observations to draw to evaluate the model checking p-value for optimization procedures

L

number of observables per component \(n\)

Type:

int

N

total number of pixels to reconstruct

Type:

int

generate_random_start_Theta(posterior: Posterior)

generates a random element of the hypercube defined by the lower and upper bounds with uniform distribution

Parameters:

posterior (Posterior) – contains the lower and upper bounds of the hypercube

Returns:

Theta – random element of the hypercube defined by the lower and upper bounds with uniform distribution

Return type:

np.array of shape (N, D)

abstract generate_random_start_Theta_1pix(Theta: ndarray, posterior: Posterior, idx_pix: ndarray) → ndarray

draws a random vectors for components \(n\) (e.g., a pixel \(\theta_n\)).

Parameters:

• Theta (np.ndarray) – current iterate

• posterior (Posterior) – contains the lower and upper bounds of the hypercube

• idx_pix (np.ndarray) – indices of the pixels

Returns:

random element of the hypercube defined by the lower and upper bounds with uniform distribution

Return type:

np.array of shape (n_pix, self.k_mtm, D)

get_rng_state()

extracts the current state and inc of the random generator

Returns:

• rng_state_array (np.array of ints) – current state of the random generator

• rng_inc_array (np.array of ints) – current inc of the random generator

Note

1. state and inc are very large integers, and thus need to be converted to hex format (later to an array of ints) to be saved in a .h5 file https://docs.python.org/3/library/stdtypes.html#int.to_bytes

2. need 32 bytes in length: otherwise, the inverse operation int.from_bytes(state_array,byteorder) does not coincide with the original int value

rng

numpy random generator (for sampling reproducibility)

Type:

np.random.Generator

abstract sample(posterior: Posterior, saver: Saver, max_iter: int, Theta_0: ndarray | None, T_BI: int | None, **kwargs) → None

main method of the class, runs the sampler on the input posterior distribution

Parameters:

• posterior (Posterior) – probability distribution to be sampled

• saver (Saver) – object responsible for progressively saving the Markov chain data during the run

• max_iter (int) – total duration of a Markov chain

• Theta_0 (Optional[np.ndarray], optional) – starting point

• T_BI (int, optional) – duration of the Burn-in phase

sample_regu_hyperparams(posterior: Posterior, regu_weights_optimizer: EBayesMMLE, t: int, current_Theta: ndarray) → ndarray

updates the spatial regularization weight vector (denoted \(\tau\) in Palud et al. [2023]) using the marginal likelihood optimizer defined in Vidal et al. [2020]

Parameters:

• posterior (Posterior) – posterior distribution

• regu_weights_optimizer (EBayesMMLE) – marginal likelihood optimizer

• t (int) – time step

• current_Theta (np.ndarray) – current iterate

Returns:

updated spatial regularization weight vector

Return type:

np.ndarray

set_rng_state(state_bytes, inc_bytes)

sets the state and inc of the random generator to a specific value. Enables to re-launch a sampling when the rng state is regularly saved (see MySaver class).

Parameters:

• state_bytes (bytes) – state to be set

• inc_bytes (bytes) – inc to be set

beetroots.sampler.my_sampler module

Contains a class of sampler used in the Meudon PDR code Bayesian inversion problems

class beetroots.sampler.my_sampler.MySampler(my_sampler_params: ~beetroots.sampler.utils.my_sampler_params.MySamplerParams, D: int, L: int, N: int, rng: ~numpy.random._generator.Generator = Generator(PCG64) at 0x7FC1608BFAC0)

Bases: Sampler

Defines the sampler proposed in Palud et al. [2023] that randomly combines two transition kernels :

1. a independent MTM-chromatic Gibbs transition kernel

2. a position-dependent MALA transition kernel with the RMSProp pre-conditioner

Parameters:

• my_sampler_params (MySamplerParams) – contains the main parameters of the algorithm

• D (int) – total number of physical parameters to reconstruct

• L (int) – number of observables per component \(n\)

• N (int) – total number of pixels to reconstruct

• rng (numpy.random.Generator, optional) – random number generator (for reproducibility), by default np.random.default_rng(42)

D

total number of physical parameters to reconstruct

Type:

int

L

number of observables per component \(n\)

Type:

int

N

total number of pixels to reconstruct

Type:

int

alpha

weight of past values of \(v\) in the exponential decay (cf RMSProp preconditioner), denoted \(\alpha \in ]0,1[\) in the article

Type:

float

compute_correction_term

wether or not to use the correction term (denoted \(\gamma\) in the article) during the sampling (only used if is_stochastic=True)

Type:

bool

compute_derivatives_2nd_order

wether to compute the expensive second order derivatives terms. Only true when the sampler runs a Markov chain (2nd order never used in optimization) and when the correction term denoted \(\gamma\) is to be computed.

Type:

bool

current

contains all the data about the current iterate (including the evaluations of the forward map and derivatives, etc.)

Type:

dict

eps0

step size used in the Position-dependent MALA transition kernel, denoted \(\epsilon > 0\) in the article

Type:

float

generate_new_sample_mtm(t: int, posterior: Posterior)

generates a new sample using the MTM transition kernel

Parameters:

• t (int) – current iteration index

• posterior (Posterior) – target posterior distribution to sample from

Returns:

• accepted (bool) – wether or not the candidate was accepted

• log_proba_accept (float) – log of the acceptance proba

generate_new_sample_pmala_rmsprop(t: int, posterior: Posterior)

generates a new sample using the position-dependent MALA transition kernel

Parameters:

• t (int) – current iteration index

• posterior (Posterior) – negative log posterior class

Returns:

• accepted (bool) – wether or not the candidate was accepted

• log_proba_accept (float) – log of the acceptance proba

generate_random_start_Theta_1pix(Theta: ndarray, posterior: Posterior, idx_pix: ndarray) → ndarray

draws a random vectors for components \(n\) (e.g., a pixel \(\theta_n\)). The distribution used to draw these vectors is:

• the smooth indicator prior

• a combination of the smooth indicator prior and of a Gaussian mixture defined with the set of all combinations of neighbors of component \(n\)

Parameters:

• Theta (np.ndarray) – current iterate

• posterior (Posterior) – contains the lower and upper bounds of the hypercube

• idx_pix (np.ndarray) – indices of the pixels

Returns:

random element of the hypercube defined by the lower and upper bounds with uniform distribution

Return type:

np.array of shape (n_pix, self.k_mtm, D)

:raises ValueError : if posterior.prior_indicator is None:

k_mtm

number of candidates in the MTM kernel, denoted \(K\) in the article

Type:

int

lambda_

limit value that avoids division by zero when computing the RMSProp preconditioner, denoted \(\eta > 0\) in the article

Type:

float

rng

random number generator (for reproducibility)

Type:

numpy.random.Generator

sample(posterior: Posterior, saver: Saver, max_iter: int, Theta_0: ndarray | None = None, disable_progress_bar: bool = False, regu_spatial_N0: int | float = inf, regu_spatial_scale: float = 1.0, regu_spatial_vmin: float = 1e-08, regu_spatial_vmax: float = 100000000.0, T_BI: int = 0) → None

main method of the class, runs the sampler

Parameters:

• posterior (Posterior) – probability distribution to be sampled

• saver (Saver) – object responsible for progressively saving the Markov chain data during the run

• max_iter (int) – total duration of a Markov chain

• Theta_0 (Optional[np.ndarray], optional) – starting point, by default None

• disable_progress_bar (bool, optional) – wether to disable the progress bar, by default False

• regu_spatial_N0 (Union[int, float], optional) – number of iterations defining the initial update phase (for spatial regularization weight optimization). np.infty means that the optimization phase never starts, and that the weight optimization is not applied. by default np.infty

• regu_spatial_scale (Optional[float], optional) – scale parameter involved in the definition of the projected gradient step size (for spatial regularization weight optimization). by default 1.0

• regu_spatial_vmin (Optional[float], optional) – lower limit of the admissible interval (for spatial regularization weight optimization), by default 1e-8

• regu_spatial_vmax (Optional[float], optional) – upper limit of the admissible interval (for spatial regularization weight optimization), by default 1e8

• T_BI (int, optional) – duration of the Burn-in phase, by default 0

selection_probas

vector of selection probabilities for the MTM and PMALA kernels, respectively, i.e., \([p_{MTM}, 1 - p_{MTM}]\)

Type:

np.ndarray

stochastic

if True, the algorithm performs sampling, and optimization otherwise

Type:

bool

v

RMSProp gradient variance vector, denoted \(v\) in the article

Type:

np.ndarray

Module contents

Contains the definition of the sampler and of a saver class