micmac.likelihood.pixel module

class MicmacSampler(nside, lmax, nstokes, frequency_array, freq_inverse_noise, pos_special_freqs=[0, -1], freq_noise_c_ell=None, n_components=3, lmin=2, n_iter=8, limit_iter_cg=200, limit_iter_cg_eta=200, tolerance_CG=1e-08, atol_CG=1e-08, mask=None, save_CMB_chain_maps=False, save_eta_chain_maps=False, save_all_Bf_params=True, save_s_c_spectra=False, sample_r_Metropolis=True, sample_C_inv_Wishart=False, perturbation_eta_covariance=True, simultaneous_accept_rate=False, non_centered_moves=False, save_intermediary_centered_moves=False, limit_r_value=False, min_r_value=0, biased_version=False, classical_Gibbs=False, use_binning=False, bin_ell_distribution=None, acceptance_posdef=False, step_size_r=0.0001, covariance_Bf=None, use_scam_step_size=False, burn_in_scam=50, s_param_scam=5.76, epsilon_param_scam_r=1e-10, epsilon_param_scam_Bf=1e-11, scam_iteration_updates=50, indexes_free_Bf=False, number_iterations_sampling=100, number_iterations_done=0, seed=0, disable_chex=True, instrument_name='SO_SAT', spv_nodes_b=[])[source]

Bases: SamplingFunctions

Main MICMAC pixel sampling object to initialize and launch the Gibbs sampling in pixel domain.

The Gibbs sampling will always store Bf and r (or C) parameters

Parameters:
nside : int

nside of the input frequency maps

lmax : int

maximum multipole for the spherical harmonics transforms and harmonic domain objects,

nstokes : int

number of Stokes parameters

frequency_array : array[float]

array of frequencies, in GHz

freq_inverse_noise : array[float]

array of inverse noise for each frequency, in uK^-2

pos_special_freqs : list[int] (optional)

indexes of the special frequencies in the frequency array respectively for synchrotron and dust, default is [0,-1] for first and last frequencies

freq_noise_c_ell : array[float] of dimensions [frequencies, frequencies, lmax+1-lmin] or [frequencies, frequencies, lmax] (in which case it will be cut to lmax+1-lmin)

optional, noise power spectra for each frequency, in uK^2, dimensions

n_components : int (optional)

number of components for the mixing matrix, default 3

lmin : int (optional)

minimum multipole for the spherical harmonics transforms and harmonic domain objects, default 2

n_iter : int (optional)

number of iterations the spherical harmonics transforms (for map2alm transformations), default 8

limit_iter_cg : int (optional)

maximum number of iterations for the conjugate gradient for the CMB map sampling, default 200

limit_iter_cg_eta : int (optional)

maximum number of iterations for the conjugate gradient for eta maps sampling, default 200

tolerance_CG : float (optional)

tolerance for the conjugate gradient, default 1e-8

atol_CG : float (optional)

absolute tolerance for the conjugate gradient, default 1e-8

mask : None or array[float] of dimensions [n_pix] (optional)

mask to use in the sampling ; if not given, no mask is used, default None Note: the mask WILL NOT be applied to the input maps, it will be only used for the propagated noise covariance

save_CMB_chain_maps : bool (optional)

save the CMB chain maps, default False

save_eta_chain_maps : bool (optional)

save the eta chain maps, default False

sample_r_Metropolis : bool (optional)

sample r with a Metropolis-within-Gibbs step from the BB power spectrum of the reconstructed sample of the CMB map during the Gibbs iteration, default True Either sample_r_Metropolis or sample_C_inv_Wishart should True but not both

sample_C_inv_Wishart : bool (optional)

sample C_inv with Wishart distribution instead of simply r being sampled, default False Either sample_r_Metropolis or sample_C_inv_Wishart should True but not both

limit_r_value : bool (optional)

limit r value being sampled with the minmum r value given by min_r_value, default False

min_r_value : float (optional)

minimum r value accepted for r sample if limit_r_value is True, default 0

perturbation_eta_covariance : bool (optional)

approach to compute difference between CMB noise component for eta log proba instead of repeating the CG for each Bf sampling, default True

simultaneous_accept_rate : bool (optional)

use the simultaneous accept rate for the patches of the Bf sampling, default False

biased_version : bool (optional)

use the biased version of the likelihood, so no computation of the correction term, default False

classical_Gibbs : bool (optional)

sampling only for s_c and the CMB covariance, and neither Bf or eta, default False

use_binning : bool (optional)

use binning for the sampling of inverse Wishart CMB covariance, if False bin_ell_distribution will not be used, default False

bin_ell_distribution : array[int] (optional)

binning distribution for the sampling of inverse Wishart CMB covariance, default None

acceptance_posdef : accept only positive definite matrices C sampling, bool

step_size_r : float (optional)

step size for the Metropolis-Hastings sampling of r, default 1e-4

covariance_Bf : None or array[float] of dimensions [(n_frequencies-len(pos_special_freqs))*(n_components-1), (n_frequencies-len(pos_special_freqs))*(n_components-1)] (optional)

covariance for the Metropolis-Hastings sampling of Bf ; will be repeated if multiresoltion case, default None

use_scam_step_size : bool (optional)

use the SCAM step size for the Metropolis-Hastings sampling of Bf and r (Haario et al. 2005), default False

burn_in_scam : int (optional)

number of burn-in iterations before using adaptive step-size (SCAM), not used if use_scam_steps_size is False, default 50

s_param_scam : float (optional)

s parameter for the SCAM step size (see Haario et al. 2001, Haario et al. 2005), default (2.4)**2

epsilon_param_scam_r : float (optional)

epsilon parameter for the SCAM step size for r (see Haario et al. 2001, Haario et al. 2005), default 1e-10

epsilon_param_scam_Bf : float (optional)

epsilon parameter for the SCAM step size for Bf (see Haario et al. 2001, Haario et al. 2005), default 1e-11

scam_iteration_updates : int (optional)

number of iterations for which the SCAM step size will be updated (the variance is updated successively for every scam_iteration_updates iterations), default 100

indexes_free_Bf : bool or array[int] (optional)

indexes of the free Bf parameters to actually sample and leave the rest of the indices fixed, array of integers, default False to sample all Bf

number_iterations_sampling : int (optional)

maximum number of iterations for the sampling, default 100

number_iterations_done : int (optional)

number of iterations already accomplished, in case the chain is resuming from a previous run, usually set by exterior routines, default 0

seed : int or array[jnp.uint32] (optional)

seed for the JAX PRNG random number generator to start the chain or array of a previously computed seed, default 0

disable_chex : bool (optional)

disable chex tests (to improve speed)

instrument_name : str (optional)

name of the instrument as expected by cmbdb or given as ‘customized_instrument’ if redefined by user, default ‘SO_SAT’ see https://github.com/dpole/cmbdb/blob/master/cmbdb/experiments.yaml

fwhm : float (optional)

FWHM of the beam in arcmin, default None (no beam) ; not implemented yet

spv_nodes_b : list[dictionaries] (optional)

tree for the spatial variability, to generate from a yaml file, default [] in principle set up by get_nodes_b

__init__(nside, lmax, nstokes, frequency_array, freq_inverse_noise, pos_special_freqs=[0, -1], freq_noise_c_ell=None, n_components=3, lmin=2, n_iter=8, limit_iter_cg=200, limit_iter_cg_eta=200, tolerance_CG=1e-08, atol_CG=1e-08, mask=None, save_CMB_chain_maps=False, save_eta_chain_maps=False, save_all_Bf_params=True, save_s_c_spectra=False, sample_r_Metropolis=True, sample_C_inv_Wishart=False, perturbation_eta_covariance=True, simultaneous_accept_rate=False, non_centered_moves=False, save_intermediary_centered_moves=False, limit_r_value=False, min_r_value=0, biased_version=False, classical_Gibbs=False, use_binning=False, bin_ell_distribution=None, acceptance_posdef=False, step_size_r=0.0001, covariance_Bf=None, use_scam_step_size=False, burn_in_scam=50, s_param_scam=5.76, epsilon_param_scam_r=1e-10, epsilon_param_scam_Bf=1e-11, scam_iteration_updates=50, indexes_free_Bf=False, number_iterations_sampling=100, number_iterations_done=0, seed=0, disable_chex=True, instrument_name='SO_SAT', spv_nodes_b=[])[source]

Main MICMAC pixel sampling object to initialize and launch the Gibbs sampling in pixel domain.

The Gibbs sampling will always store Bf and r (or C) parameters

Parameters:
nside : int

nside of the input frequency maps

lmax : int

maximum multipole for the spherical harmonics transforms and harmonic domain objects,

nstokes : int

number of Stokes parameters

frequency_array : array[float]

array of frequencies, in GHz

freq_inverse_noise : array[float]

array of inverse noise for each frequency, in uK^-2

pos_special_freqs : list[int] (optional)

indexes of the special frequencies in the frequency array respectively for synchrotron and dust, default is [0,-1] for first and last frequencies

freq_noise_c_ell : array[float] of dimensions [frequencies, frequencies, lmax+1-lmin] or [frequencies, frequencies, lmax] (in which case it will be cut to lmax+1-lmin)

optional, noise power spectra for each frequency, in uK^2, dimensions

n_components : int (optional)

number of components for the mixing matrix, default 3

lmin : int (optional)

minimum multipole for the spherical harmonics transforms and harmonic domain objects, default 2

n_iter : int (optional)

number of iterations the spherical harmonics transforms (for map2alm transformations), default 8

limit_iter_cg : int (optional)

maximum number of iterations for the conjugate gradient for the CMB map sampling, default 200

limit_iter_cg_eta : int (optional)

maximum number of iterations for the conjugate gradient for eta maps sampling, default 200

tolerance_CG : float (optional)

tolerance for the conjugate gradient, default 1e-8

atol_CG : float (optional)

absolute tolerance for the conjugate gradient, default 1e-8

mask : None or array[float] of dimensions [n_pix] (optional)

mask to use in the sampling ; if not given, no mask is used, default None Note: the mask WILL NOT be applied to the input maps, it will be only used for the propagated noise covariance

save_CMB_chain_maps : bool (optional)

save the CMB chain maps, default False

save_eta_chain_maps : bool (optional)

save the eta chain maps, default False

sample_r_Metropolis : bool (optional)

sample r with a Metropolis-within-Gibbs step from the BB power spectrum of the reconstructed sample of the CMB map during the Gibbs iteration, default True Either sample_r_Metropolis or sample_C_inv_Wishart should True but not both

sample_C_inv_Wishart : bool (optional)

sample C_inv with Wishart distribution instead of simply r being sampled, default False Either sample_r_Metropolis or sample_C_inv_Wishart should True but not both

limit_r_value : bool (optional)

limit r value being sampled with the minmum r value given by min_r_value, default False

min_r_value : float (optional)

minimum r value accepted for r sample if limit_r_value is True, default 0

perturbation_eta_covariance : bool (optional)

approach to compute difference between CMB noise component for eta log proba instead of repeating the CG for each Bf sampling, default True

simultaneous_accept_rate : bool (optional)

use the simultaneous accept rate for the patches of the Bf sampling, default False

biased_version : bool (optional)

use the biased version of the likelihood, so no computation of the correction term, default False

classical_Gibbs : bool (optional)

sampling only for s_c and the CMB covariance, and neither Bf or eta, default False

use_binning : bool (optional)

use binning for the sampling of inverse Wishart CMB covariance, if False bin_ell_distribution will not be used, default False

bin_ell_distribution : array[int] (optional)

binning distribution for the sampling of inverse Wishart CMB covariance, default None

acceptance_posdef : accept only positive definite matrices C sampling, bool

step_size_r : float (optional)

step size for the Metropolis-Hastings sampling of r, default 1e-4

covariance_Bf : None or array[float] of dimensions [(n_frequencies-len(pos_special_freqs))*(n_components-1), (n_frequencies-len(pos_special_freqs))*(n_components-1)] (optional)

covariance for the Metropolis-Hastings sampling of Bf ; will be repeated if multiresoltion case, default None

use_scam_step_size : bool (optional)

use the SCAM step size for the Metropolis-Hastings sampling of Bf and r (Haario et al. 2005), default False

burn_in_scam : int (optional)

number of burn-in iterations before using adaptive step-size (SCAM), not used if use_scam_steps_size is False, default 50

s_param_scam : float (optional)

s parameter for the SCAM step size (see Haario et al. 2001, Haario et al. 2005), default (2.4)**2

epsilon_param_scam_r : float (optional)

epsilon parameter for the SCAM step size for r (see Haario et al. 2001, Haario et al. 2005), default 1e-10

epsilon_param_scam_Bf : float (optional)

epsilon parameter for the SCAM step size for Bf (see Haario et al. 2001, Haario et al. 2005), default 1e-11

scam_iteration_updates : int (optional)

number of iterations for which the SCAM step size will be updated (the variance is updated successively for every scam_iteration_updates iterations), default 100

indexes_free_Bf : bool or array[int] (optional)

indexes of the free Bf parameters to actually sample and leave the rest of the indices fixed, array of integers, default False to sample all Bf

number_iterations_sampling : int (optional)

maximum number of iterations for the sampling, default 100

number_iterations_done : int (optional)

number of iterations already accomplished, in case the chain is resuming from a previous run, usually set by exterior routines, default 0

seed : int or array[jnp.uint32] (optional)

seed for the JAX PRNG random number generator to start the chain or array of a previously computed seed, default 0

disable_chex : bool (optional)

disable chex tests (to improve speed)

instrument_name : str (optional)

name of the instrument as expected by cmbdb or given as ‘customized_instrument’ if redefined by user, default ‘SO_SAT’ see https://github.com/dpole/cmbdb/blob/master/cmbdb/experiments.yaml

fwhm : float (optional)

FWHM of the beam in arcmin, default None (no beam) ; not implemented yet

spv_nodes_b : list[dictionaries] (optional)

tree for the spatial variability, to generate from a yaml file, default [] in principle set up by get_nodes_b

property all_samples_s_c

Returns all the CMB sampled maps from the initial WF and fluctuation maps

generate_CMB(return_spectra=True)[source]

Returns CMB spectra of scalar modes only and tensor modes only (with r=1) Both CMB spectra are either returned in the usual form [number_correlations,lmax+1], or in the red_cov form if return_spectra == False

generate_input_freq_maps_from_fgs(freq_maps_fgs, r_true=0, return_only_freq_maps=True, return_only_maps=False)[source]

Generate input frequency maps (CMB+foregrounds) from the input frequency foregrounds maps, return either the full frequency maps, the full frequency and CMB maps alone, or the full frequency and CMB maps with the theoretical reduced covariance matrices for the CMB scalar and tensor modes

Parameters:
freq_maps_fgs : array[float] of dimensions [n_frequencies,nstokes,n_pix]

input frequency foregrounds maps

r_true : float, optional

input tensor-to-scalar ratio r to generate the CMB, default to 0

return_only_freq_maps : bool (optional)

return only the full frequency maps, bool

return_only_maps : bool (optional)

return only the full frequency and CMB maps alone, bool

Returns:

  • input_freq_maps (array[float] of dimensions [n_frequencies,nstokes,n_pix]) – input frequency maps

  • input_cmb_maps (array[float] of dimensions [nstokes,n_pix]) – input CMB maps

  • theoretical_red_cov_r0_total (array[float] of dimensions [lmax+1-lmin,nstokes,nstokes]) – theoretical reduced covariance matrix for the CMB scalar modes

  • theoretical_red_cov_r1_tensor (array[float] of dimensions [lmax+1-lmin,nstokes,nstokes]) – theoretical reduced covariance matrix for the CMB tensor modes

perform_Gibbs_sampling(input_freq_maps, c_ell_approx, CMB_c_ell, init_params_mixing_matrix, initial_guess_r=1e-08, initial_wiener_filter_term=None, initial_fluctuation_maps=None, theoretical_r0_total=None, theoretical_r1_tensor=None, **dictionnary_additional_parameters)[source]
Perform sampling steps with:

  1. The sampling of eta by computing eta = x + C_approx^(1/2) N_c^{-1/2} y ; where x is band-limited

  2. A CG for the Wiener filter (WF) and fluctuation variables s_c: (s_c - s_{c,WF})^t (C^{-1} + N_c^{-1}) (s_c - s_{c,WF})

  3. The c_ell sampling, either by parametrizing it by r or by sampling an inverse Wishart distribution

  4. Mixing matrix Bf sampling with: -(d - B_c s_c)^t N^{-1} B_f (B_f^t N^{-1} B_f)^{-1} B_f^t N^{-1} (d - B_c s_c) + eta^t (Id + C_{approx}^{1/2} N_c^{-1} C_{approx}^{1/2}) eta

The results of the chain will be stored in the class attributes, depending if the save options are put to True or False:

  • self.all_samples_eta (if self.save_eta_chain_maps is True)

  • self.all_samples_wiener_filter_maps (if self.save_CMB_chain_maps is True)

  • self.all_samples_fluctuation_maps (if self.save_CMB_chain_maps is True)

  • self.all_samples_r (if self.sample_r_Metropolis is True)

  • self.all_samples_CMB_c_ell (if self.sample_C_inv_Wishart is True)

  • self.all_params_mixing_matrix_samples (always)

This same function can be used to continue a chain from a previous run, by giving the number of iterations already done in the MicmacSampler object, giving the chains to the attributes of the object, and giving the last iteration results as initial guesses.

Parameters:
input_freq_maps : array[float] of dimensions [frequencies, nstokes, n_pix]

input frequency maps

c_ell_approx : array[float] of dimensions [number_correlations, lmax+1]

approximate CMB power spectra for the latent parameter eta defining the ad-hoc correction term

CMB_c_ell : array[float] of dimensions [number_correlations, lmax+1]

CMB power spectra, where number_correlations is the number of auto- and cross-correlations relevant considering the number of Stokes parameters

init_params_mixing_matrix : array[float] of dimensions [len_params]

initial parameters for the mixing matrix elements Bf; expected to be given flattened as [Bf_s1, Bf_s2, …, Bf_sn, Bf_d1, …, Bf_dn]

initial_guess_r : float (optional)

initial guess for r, default 1e-8

initial_wiener_filter_term : array[float] of dimensions [nstokes, n_pix] or empty (optional)

initial guess for the Wiener filter term, default empty array

initial_fluctuation_maps : array[float] of dimensions [nstokes, n_pix] or empty (optional)

initial guess for the fluctuation maps, default empty array

theoretical_r0_total : array[float] of dimensions [number_correlations, lmax+1-lmin] (optional)

theoretical reduced covariance matrix for the CMB scalar modes, default empty array

theoretical_r1_tensor : array[float] of dimensions [number_correlations, lmax+1-lmin] (optional)

theoretical reduced covariance matrix for the CMB tensor modes, default empty array

dictionnary_additional_parameters : dictionary

additional parameters to give to the function, currently only the ones related to the SCAM step size

Notes

The formalism relies on the ability to have an inverse for C_approx (even though it is never computed effectively in the code), and may lead to numerical instabilities if the C_approx matrix is not well-conditioned.

update_one_sample(one_sample)[source]

Update the samples with one sample to add

update_samples(all_samples)[source]

Update the samples with new samples to add

Parameters:
all_samples : dictionary

dictionary of all the samples to update

update_variable(all_samples, new_samples_to_add)[source]

Update the samples with new samples to add by stacking them

Parameters:
all_samples : array[float] of dimensions [n_samples,n_pix]

previous samples to update

new_samples_to_add : array[float] of dimensions [n_samples,n_pix]

new samples to add

Returns:

all_samples – updated samples

Return type:

array[float] of dimensions [n_samples+n_samples,n_pix]