QuickCW package

Submodules

QuickCW.CWFastLikelihoodNumba module

C 2021 Matthew Digman and Bence Becsy numba version of fast likelihood

class QuickCW.CWFastLikelihoodNumba.CWInfo(*args, **kwargs)

Bases: CWInfo

simple jitclass to store the various parmeters in a way that can be accessed quickly from a numba environment

Parameters

• Npsr – Number of pulsars

• params_in – Array of all parameters in the same order as in par_names

• par_names – List of parameter names - must follow certain naming conventions so we can identify parameters

• par_names_cw_ext – Subset of all parameters that describe the CW signal and are projection parameters (previously called extrinsic parameters)

• par_names_cw_int – Subset of all parameters that describe the CW signal and are shape parameters (previously called intrinsic parameters)

class_type = jitclass.CWInfo#7f5a8350e7f0<Npsr:int64,cw_p_dists:array(float64, 1d, A),cw_p_phases:array(float64, 1d, A),rn_gammas:array(float64, 1d, A),rn_log10_As:array(float64, 1d, A),cos_gwtheta:float64,cos_inc:float64,gwphi:float64,log10_fgw:float64,log10_h:float64,log10_mc:float64,phase0:float64,psi:float64,gwb_gamma:float64,gwb_log10_A:float64,idx_dists:array(int64, 1d, A),idx_phases:array(int64, 1d, A),idx_rn_gammas:array(int64, 1d, A),idx_rn_log10_As:array(int64, 1d, A),idx_cos_gwtheta:int64,idx_cos_inc:int64,idx_gwphi:int64,idx_log10_fgw:int64,idx_log10_h:int64,idx_log10_mc:int64,idx_phase0:int64,idx_psi:int64,idx_gwb_gamma:int64,idx_gwb_log10_A:int64,idx_rn:array(int64, 1d, A),idx_gwb:array(int64, 1d, A),idx_int:array(int64, 1d, A),idx_cw_ext:array(int64, 1d, A),idx_cw_int:array(int64, 1d, A)>

class QuickCW.CWFastLikelihoodNumba.FastLikeInfo(*args, **kwargs)

Bases: FastLikeInfo

simple jitclass to store the various elements of fast likelihood calculation in a way that can be accessed quickly from a numba environment

Parameters

• logdet_base – Contibution to logdet without CW signal

• pos – (number of pulsars, 3) array holding 3d unit vector pointing towards each pulsar given by psr.pos

• pdist – (number of pulsars, 2) array holding distance and error of distance for each pulsar in kpc given by psr.pdist

• toas – List of arrays of TOAs for each pulsar

• Nvecs – List of N vectors

• Nrs – Residuals times inverse sqareroot N vectors

• max_toa – Maximum TOA over all pulsars

• x0 – CWInfo object

• Npsr – Number of pulsars

• isqrNvecs – Inverse squareroot of N vectors

• TNvs – T vectros times inverse squareroot N vectors

• dotTNrs – Precalculated dot product of Nrs and TNvs

• chol_Sigmas – List of Cholesky decompositions of Sigma matrices

• phiinvs – List of phiinv matrices

• includeCW – Switch if we want to include the contribution of the CW signal or not [True]

• prior_recovery – If True, we return constant likelihood to be used for prior recovery diagnostic test [False]

class_type = jitclass.FastLikeInfo#7f5a8330dcd0<resres:float64,logdet:float64,resres_array:array(float64, 1d, A),logdet_array:array(float64, 1d, A),logdet_base:float64,logdet_base_orig:float64,pos:array(float64, 2d, C),pdist:array(float64, 2d, C),toas:ListType[array(float64, 1d, C)],Npsr:int64,max_toa:float64,phiinvs:ListType[array(float64, 1d, C)],dotTNrs:ListType[array(float64, 1d, C)],Nvecs:ListType[array(float64, 1d, C)],isqrNvecs:ListType[array(float64, 1d, C)],Nrs:ListType[array(float64, 1d, C)],TNvs:ListType[array(float64, 2d, F)],chol_Sigmas:ListType[array(float64, 2d, F)],MMs:array(float64, 3d, C),NN:array(float64, 2d, C),cos_gwtheta:float64,gwphi:float64,log10_fgw:float64,log10_mc:float64,cw_p_dists:array(float64, 1d, A),gwb_gamma:float64,gwb_log10_A:float64,rn_gammas:array(float64, 1d, A),rn_log10_As:array(float64, 1d, A),includeCW:bool,prior_recovery:bool>

class QuickCW.CWFastLikelihoodNumba.FastLikeMaster(psrs, pta, params, x0, includeCW=True, prior_recovery=False)

Bases: object

class to store pta things so they do not have to be recomputed when red noise is recomputed

get Class for generating the fast CW likelihood.

Parameters

• pta – enterprise pta object.

• params – Dictionary of noise parameters.

• x0 – CWInfo object, which is partially redundant with params but better handled by numba

• includeCW – Switch if we want to include the contribution of the CW signal or not [True]

• prior_recovery – If True, we return constant likelihood to be used for prior recovery diagnostic test [False]

get_new_FastLike(x0, params)

recompute_FastLike(FLI, x0, params, chol_update=False, mask=None)

QuickCW.CWFastLikelihoodNumba.cholupdate(L_in, diag_diff)

Jitted routine to update the Cholesky of a matrix instead of recomputing it. In principle could be faster, but right now it is not, so it’s not used.

Parameters

• L_in – Previous Cholesky that we want to update

• diag_diff – Differences in the diagonal of the original matrix

Return L

Updated Cholesky

QuickCW.CWFastLikelihoodNumba.cholupdate_loop(chol_Sigmas, pls_temp, old_phiinvs, Npsr)

Jitted loop over Sigma matrices to update their Cholesky. Currently not faster than recalculating from scratch, so not used.

Parameters

• chol_Sigmas – List of Cholesky decompositions of Sigma matrices

• pls_temp – New list of tuples returned by pta.get_phiinv

• old_phiinvs – Old list of phiinv matrices

• Npsr – Number of pulsars (also number of Sigma matrices)

Return chol_Sigmas

List of updated Choleskies

Return logdet_array

New array containing logdet values

Return new_phiinvs

List of new phiinv matrices

QuickCW.CWFastLikelihoodNumba.create_Sigma(phiinv_loc, TNT, Sigma)

create just the upper triangle of the Sigma matrix with phiinv_loc added to the diagonal, lower triangle will be garbage

Parameters

• phiinv_loc – phiinv matrix

• TNT – Precomputed matrix product of TNT

• Sigma – Initialized mtrix with the right shape to hold Sigma (allows for just overwriting instead of creating new array each time we update this)

Return Sigma

Upper trinagle Sigma matrix - lower triangle is junk

QuickCW.CWFastLikelihoodNumba.get_lnlikelihood_helper(x0, resres, logdet, pos, pdist, NN, MMs, includeCW=True, prior_recovery=False)

jittable helper for calculating the log likelihood in CWFastLikelihood

Parameters

• x0 – CWInfo object

• resres – Array holding the (residual|residual) inner product

• logdet – Log determinant piece of the likelihood

• pos – (number of pulsars, 3) array holding 3d unit vector pointing towards each pulsar given by psr.pos

• pdist – (number of pulsars, 2) array holding distance and error of distance for each pulsar in kpc given by psr.pdist

• NN – N matrices holding (filter|residual) type inner products

• MMs – M matrices holding (filter|filter) type inner products

• includeCW – Switch if we want to include the contribution of the CW signal or not [True]

• prior_recovery – If True, we return constant likelihood to be used for prior recovery diagnostic test [False]

Return log_L

Log likelihood value

QuickCW.CWFastLikelihoodNumba.isclose(a, b, rtol=1e-05, atol=1e-08)

check if close in same way as np.isclose

Parameters

• a – First float to use in comparison

• b – Second float to use in comparison

• rtol – Realtive tolerance - default:1e-5

• atol – Absolute tolerance - default:1e-8

Returns

True/False indicating if a and b are close to each other

QuickCW.CWFastLikelihoodNumba.logdet_Sigma_helper(chol_Sigma)

get logdet sigma from cholesky of Sigma

Parameters

chol_Sigma – Cholesky decomposition of Sigma matrix

Return 2*res

Contribution to logdet from this Sigma matrix

QuickCW.CWFastLikelihoodNumba.update_intrinsic_params(x0, isqrNvecs, Nrs, pos, pdist, toas, NN, MMs, SigmaTNrProds, invchol_Sigma_TNs, idxs, dist_only=True)

Calculate inner products N=(res|S), M=(S|S)

Parameters

• x0 – CWInfo object

• isqrNvecs – Inverse squareroot of N vectors

• Nrs – Residuals times inverse sqareroot N vectors

• pos – (number of pulsars, 3) array holding 3d unit vector pointing towards each pulsar given by psr.pos

• pdist – (number of pulsars, 2) array holding distance and error of distance for each pulsar in kpc given by psr.pdist

• toas – List of arrays of TOAs for each pulsar

• NN – N matrices holding (filter|residual) type inner products

• MMs – M matrices holding (filter|filter) type inner products

• SigmaTNrProds –

  —

• invchol_Sigma_TNs –

  —

• idxs – Indices of pulsar for which we want to update things

• dist_only – Option to skip parts of calculation whn only updating pulsar distances - not currently used

QuickCW.CWFastLikelihoodNumba.update_intrinsic_params2(x0, isqrNvecs, Nrs, pos, pdist, toas, NN, MMs, TNvs, chol_Sigmas, idxs, resres_array, dotTNrs)

Calculate inner products N=(res|S), M=(S|S)

Parameters

• x0 – CWInfo object

• isqrNvecs – Inverse squareroot of N vectors

• Nrs – Residuals times inverse sqareroot N vectors

• pos – (number of pulsars, 3) array holding 3d unit vector pointing towards each pulsar given by psr.pos

• pdist – (number of pulsars, 2) array holding distance and error of distance for each pulsar in kpc given by psr.pdist

• toas – List of arrays of TOAs for each pulsar

• NN – N matrices holding (filter|residual) type inner products

• MMs – M matrices holding (filter|filter) type inner products

• TNvs – T vectros times inverse squareroot N vectors

• chol_Sigmas – List of Cholesky decompositions of Sigma matrices

• idxs – Indices of pulsar for which we want to update things

• resres_array – Array containing contributions to (res|res)

• dotTNrs – Precalculated dot product of Nrs and TNvs

QuickCW.CWFastPrior module

C 2021 Bence Becsy MCMC for CW fast likelihood (w/ Neil Cornish and Matthew Digman) Helpers for Getting Prior Draws

class QuickCW.CWFastPrior.FastPrior(pta, psrs, par_names_cw_ext)

Bases: object

helper class to set up information about priors

Parameters

• pta – enterprise PTA object

• psrs – list of enterprise pulsar objects

• par_names_cw_ext – list of projection parameters (previously called extrinsic parameters)

get_lnprior(x0)

wrapper to get ln prior

get_sample(idx)

wrapper to quickly return random prior draw for the (idx)th parameter

class QuickCW.CWFastPrior.FastPriorInfo(*args, **kwargs)

Bases: FastPriorInfo

simple jitclass to store the various elements of fast prior calculation in a way that can be accessed quickly from a numba environment

Parameters

• uniform_par_ids – Indices of parameters with uniform prior

• uniform_lows – Lower prior bounds of uniform prior parameters

• uniform_highs – Upper prior bounds of uniform prior parameters

• lin_exp_par_ids – Indices of parameters with linear exponential prior

• lin_exp_lows – Lower prior bounds of linear exponential prior parameters

• lin_exp_highs – Upper prior bounds of linear exponential prior parameters

• normal_par_ids – Indices of parameters with normal prior

• normal_mus – Means of normal prior parameters

• normal_sigs – Standard deviations of normal prior parameters

• dm_par_ids – Indices of parameters with DM distance prior

• dm_dists – Mean distances of DM distance prior parameters

• dm_errs – Errors of DM distance prior parameters

• px_par_ids – Indices of parameters with parallax distance prior

• px_mus – Means of parallax distance prior parameters

• px_errs – Errors of parallax distance prior parameters

• cut_par_ids – Indices of parameters where an extra cut might be needed to stay in prior range (like distance)

• cut_lows – Lower bounds of these parameters

• cut_highs – Upper bounds of these parameters

• cw_ext_par_ids – Indices of projection parameters (previously called extrinsic parameters)

• cw_ext_lows – Lower bounds of extrinsic parameters

• cw_ext_highs – Upper bounds of extrinsic parameters

• global_common – Part of the log prior independent of the parameter values

class_type = jitclass.FastPriorInfo#7f5a83192f40<uniform_par_ids:array(int64, 1d, A),uniform_lows:array(float64, 1d, A),uniform_highs:array(float64, 1d, A),lin_exp_par_ids:array(int64, 1d, A),lin_exp_lows:array(float64, 1d, A),lin_exp_highs:array(float64, 1d, A),normal_par_ids:array(int64, 1d, A),normal_mus:array(float64, 1d, A),normal_sigs:array(float64, 1d, A),dm_par_ids:array(int64, 1d, A),dm_dists:array(float64, 1d, A),dm_errs:array(float64, 1d, A),px_par_ids:array(int64, 1d, A),px_mus:array(float64, 1d, A),px_errs:array(float64, 1d, A),cut_par_ids:array(int64, 1d, A),cut_lows:array(float64, 1d, A),cut_highs:array(float64, 1d, A),cw_ext_par_ids:array(int64, 1d, A),cw_ext_lows:array(float64, 1d, A),cw_ext_highs:array(float64, 1d, A),global_common:float64>

QuickCW.CWFastPrior.get_FastPriorInfo(pta, psrs, par_names_cw_ext)

get FastPriorInfo object from pta

Parameters

• pta – enterprise PTA object

• psrs – List of enterprise pulsar objects

• par_names_cw_ext – List of parameter names which are projection parameters (previously called extrinsic parameters)

QuickCW.CWFastPrior.get_lnprior(x0, FPI)

wrapper to get lnprior from jitted helper

Parameters

• x0 – Array of parameter values for which we want to calculate the prior

• FPI – FastPriorInfo object

Returns

log prior

QuickCW.CWFastPrior.get_lnprior_array(samples, FPI)

wrapper to get lnprior from jitted helper

Parameters

• samples – 2d array of parameter sets for which we want to calculate the prior

• FPI – FastPriorInfo object

Returns

Array of log prior values

QuickCW.CWFastPrior.get_lnprior_helper(x0, uniform_par_ids, uniform_lows, uniform_highs, lin_exp_par_ids, lin_exp_lows, lin_exp_highs, normal_par_ids, normal_mus, normal_sigs, dm_par_ids, dm_dists, dm_errs, px_par_ids, px_mus, px_errs, global_common)

jittable helper for calculating the log prior

Parameters

• x0 – Array of parameters for which we want to calculate the prior

• uniform_par_ids – Indices of parameters with uniform prior

• uniform_lows – Lower prior bounds of uniform prior parameters

• uniform_highs – Upper prior bounds of uniform prior parameters

• lin_exp_par_ids – Indices of parameters with linear exponential prior

• lin_exp_lows – Lower prior bounds of linear exponential prior parameters

• lin_exp_highs – Upper prior bounds of linear exponential prior parameters

• normal_par_ids – Indices of parameters with normal prior

• normal_mus – Means of normal prior parameters

• normal_sigs – Standard deviations of normal prior parameters

• dm_par_ids – Indices of parameters with DM distance prior

• dm_dists – Mean distances of DM distance prior parameters

• dm_errs – Errors of DM distance prior parameters

• px_par_ids – Indices of parameters with parallax distance prior

• px_mus – Means of parallax distance prior parameters

• px_errs – Errors of parallax distance prior parameters

• global_common – Part of the log prior independent of the parameter values

Return log_prior

Log prior

QuickCW.CWFastPrior.get_lnprior_helper_array(x0s, uniform_par_ids, uniform_lows, uniform_highs, lin_exp_par_ids, lin_exp_lows, lin_exp_highs, normal_par_ids, normal_mus, normal_sigs, dm_par_ids, dm_dists, dm_errs, px_par_ids, px_mus, px_errs, global_common)

jittable helper for calculating the log prior for an array of points same as get_lnprior_helper, but can be called on an array of parameters

Parameters

• x0s – 2D array of multiple parameter sets for which we want to calculate the prior

• uniform_par_ids – Indices of parameters with uniform prior

• uniform_lows – Lower prior bounds of uniform prior parameters

• uniform_highs – Upper prior bounds of uniform prior parameters

• lin_exp_par_ids – Indices of parameters with linear exponential prior

• lin_exp_lows – Lower prior bounds of linear exponential prior parameters

• lin_exp_highs – Upper prior bounds of linear exponential prior parameters

• normal_par_ids – Indices of parameters with normal prior

• normal_mus – Means of normal prior parameters

• normal_sigs – Standard deviations of normal prior parameters

• dm_par_ids – Indices of parameters with DM distance prior

• dm_dists – Mean distances of DM distance prior parameters

• dm_errs – Errors of DM distance prior parameters

• px_par_ids – Indices of parameters with parallax distance prior

• px_mus – Means of parallax distance prior parameters

• px_errs – Errors of parallax distance prior parameters

• global_common – Part of the log prior independent of the parameter values

Return log_priors

Array of log prior values

QuickCW.CWFastPrior.get_sample_full(n_par, FPI)

helper to get a full prior draw sample

Parameters

• n_par – Number of parameters

• FPI – FastPriorInfo object

Return new_point

Array with parameters drawn from their priors

QuickCW.CWFastPrior.get_sample_helper(idx, uniform_par_ids, uniform_lows, uniform_highs, lin_exp_par_ids, lin_exp_lows, lin_exp_highs, normal_par_ids, normal_mus, normal_sigs, dm_par_ids, dm_dists, dm_errs, px_par_ids, px_mus, px_errs)

jittable helper for prior draws - same as get_sample_helper_full but only for a single parameter

Parameters

• idx – Index of parameters for which we want a prior draw

• uniform_par_ids – Indices of parameters with uniform prior

• uniform_lows – Lower prior bounds of uniform prior parameters

• uniform_highs – Upper prior bounds of uniform prior parameters

• lin_exp_par_ids – Indices of parameters with linear exponential prior

• lin_exp_lows – Lower prior bounds of linear exponential prior parameters

• lin_exp_highs – Upper prior bounds of linear exponential prior parameters

• normal_par_ids – Indices of parameters with normal prior

• normal_mus – Means of normal prior parameters

• normal_sigs – Standard deviations of normal prior parameters

• dm_par_ids – Indices of parameters with DM distance prior

• dm_dists – Mean distances of DM distance prior parameters

• dm_errs – Errors of DM distance prior parameters

• px_par_ids – Indices of parameters with parallax distance prior

• px_mus – Means of parallax distance prior parameters

• px_errs – Errors of parallax distance prior parameters

Returns

Prior draw of the parameter of interest

QuickCW.CWFastPrior.get_sample_helper_full(n_par, uniform_par_ids, uniform_lows, uniform_highs, lin_exp_par_ids, lin_exp_lows, lin_exp_highs, normal_par_ids, normal_mus, normal_sigs, dm_par_ids, dm_dists, dm_errs, px_par_ids, px_mus, px_errs)

jittable helper for prior draws

Parameters

• n_par – Number of parameters in likelihood/prior

• uniform_par_ids – Indices of parameters with uniform prior

• uniform_lows – Lower prior bounds of uniform prior parameters

• uniform_highs – Upper prior bounds of uniform prior parameters

• lin_exp_par_ids – Indices of parameters with linear exponential prior

• lin_exp_lows – Lower prior bounds of linear exponential prior parameters

• lin_exp_highs – Upper prior bounds of linear exponential prior parameters

• normal_par_ids – Indices of parameters with normal prior

• normal_mus – Means of normal prior parameters

• normal_sigs – Standard deviations of normal prior parameters

• dm_par_ids – Indices of parameters with DM distance prior

• dm_dists – Mean distances of DM distance prior parameters

• dm_errs – Errors of DM distance prior parameters

• px_par_ids – Indices of parameters with parallax distance prior

• px_mus – Means of parallax distance prior parameters

• px_errs – Errors of parallax distance prior parameters

Return res

Array of prior draws

QuickCW.CWFastPrior.get_sample_idxs(old_point, idx_choose, FPI)

get just some indexes drawn from a prior

Parameters

• old_point – Old array of parameters

• idx_choose – List of indices for which we want to drawn from the prior

• FPI – FastPriorInfo object

Return new_point

New array of parameters

QuickCW.OutputUtils module

C 2021 Bence Becsy MCMC for CW fast likelihood (w/ Neil Cornish and Matthew Digman) Helpers to print outputs

QuickCW.OutputUtils.output_hdf5_end(chain_params, samples, log_likelihood, acc_fraction, fisher_diag, par_names, verbosity)

output to hdf5 file at end of body of loop

Parameters

• chain_params – ChainParams object

• samples – Array with posterior samples

• log_likelihood – Array with log likelihood values corresponding to samples

• acc_fraction – Array with acceptance rates

• diag (fisher) – Array with the digonal elements of fisher matrix

• par_names – List of parameter names

• verbosity – Parameter indicating how much info to print (higher value means more prints)

QuickCW.OutputUtils.output_hdf5_loop(itrn, chain_params, samples, log_likelihood, acc_fraction, fisher_diag, par_names, N, verbosity)

output to hdf5 at loop iteration

Parameters

• itrn – Overall index (as opposed to the block index itri or the index within saved values itrb)

• chain_params – ChainParams object

• samples – Array with posterior samples

• log_likelihood – Array with log likelihood values corresponding to samples

• acc_fraction – Array with acceptance rates

• diag (fisher) – Array with the digonal elements of fisher matrix

• par_names – List of parameter names

• N – Total number of samples

• verbosity – Parameter indicating how much info to print (higher value means more prints)

QuickCW.OutputUtils.print_acceptance_progress(itrn, N, n_int_block, a_yes, a_no, t_itr, ti_loop, tf1_loop, Ts, verbosity)

print the acceptance fraction

Parameters

• itrn – Overall index (as opposed to the block index itri or the index within saved values itrb)

• N – Total number of iterations we asked for at startup

• n_int_block – Number of iterations within a block

• a_yes – Array holding number of accepted steps

• a_no – Array holding number of rejected steps

• t_itr – Time just before calling this function got from time.perf_counter()

• ti_loop – Time after initialization got from time.perf_counter()

• tf1_loop – Time at start of this block got from time.perf_counter()

• Ts – Temperature ladder

• verbosity – Parameter indicating how much info to print (higher value means more prints)

QuickCW.PulsarDistPriors module

QuickCW.PulsarDistPriors.DMDistParameter(dist=0, err=1, size=None)

Class factory for DM distance parameters with a pdf that is flat for dist+-err and a half Gaussian beyond that

Parameters

• dist – mean distance

• err – distance error

• size – length for vector parameter

Returns

DMDist parameter class

QuickCW.PulsarDistPriors.DMDistPrior(value, dist, err)

Prior function for DMDist parameters.

Parameters

• value – point where we want the prior evaluated

• dist – mean distance

• err – distance error

Returns

prior value

QuickCW.PulsarDistPriors.DMDistSampler(dist, err, size=None)

Sampling function for DMDist parameters.

Parameters

• dist – mean distance

• err – distance error

• size – length for vector parameter

Returns

random draw from prior (float or ndarray with lenght size)

QuickCW.PulsarDistPriors.PXDistParameter(dist=0, err=1, size=None)

Class factory for PX distance parameters with a pdf of inverse Gaussian (since parallax is Gaussian)

Parameters

• dist – mean distance

• err – distance error

• size – length for vector parameter

Returns

PXDist parameter class

QuickCW.PulsarDistPriors.PXDistPrior(value, dist, err)

Prior function for PXDist parameters.

Parameters

• value – point where we want the prior evaluated

• dist – mean distance

• err – distance error

Returns

prior value

QuickCW.PulsarDistPriors.PXDistSampler(dist, err, size=None)

Sampling function for PXDist parameters.

Parameters

• dist – mean distance

• err – distance error

• size – length for vector parameter

Returns

random draw from prior (float or ndarray with lenght size)

QuickCW.QuickCW module

C 2021 Bence Becsy MCMC for CW fast likelihood (w/ Neil Cornish and Matthew Digman)

QuickCW.QuickCW.QuickCW(chain_params, psrs, noise_json=None, use_legacy_equad=False, include_ecorr=True, amplitude_prior='UL', gwb_gamma_prior=None, psr_distance_file=None, backend_selection=True)

Set up all essential objects for QuickCW to do MCMC iterations

Parameters

• chain_params – ChainParams object

• psrs – enterprise pulsar objects

• noise_json – JSON file with noise dictionary [None]

• use_legacy_equad – Option to use old convention for equad [False]

• include_ecorr – Option to include ECORR white noise [True]

• amplitude_prior – Prior to use on CW amplitude; ‘UL’ indicates uniform in amplitude prior (used for upper limits); ‘detection’ indicates uniform in log amplitude prior (used for Bayes factor calculation/detection)

• gwb_gamma_prior – Option to specify prior range on GWB spectral index gamma; None means we use the default np.array([0,7]) [None]

• psr_distance_file – File containing parallax and DM distance information for pulsars; If None, we use Gaussian prior with pulsar distance and error from psr objects [None]

• backend_selection – Option to use an enterprise Selection based on backend; Usually use True for real data False for simulated data [True]

Return pta

enterprise PTA object

Return mcc

MCMCChain onject

QuickCW.QuickCW.get_default_args(func)

Gets default arguments from a python function

QuickCW.QuickCW.per_pulsar_prior(enterprise_pulsar: enterprise.pulsar.Pulsar, pulsar_distances: dict, cw_delay_args: Optional[dict] = None, CWSignal_args: Optional[dict] = None)

Creates a CW signal applying distance priors to individual pulsars based on DM or PX in the pulsar_distances dict

Parameters

• enterprise_pulsar – enterprise pulsar object

• pulsar_distances – dictionary containing pulsar distance info

• cw_delay_args – arguments to be passed on to deterministic.cw_delay

• CWSignal_args – arguments to be passed on to deterministic.CWSignal

Return cw

enterprise signal object with the CW model

QuickCW.QuickCorrectionUtils module

C 2021 Bence Becsy MCMC for CW fast likelihood (w/ Neil Cornish and Matthew Digman) utils for correcting parameters to nominal ranges

QuickCW.QuickCorrectionUtils.check_merged(log10_fgw, log10_mc, max_toa)

check the maximum toa is not such that the source has already merged, and if so draw new parameters to avoid starting from nan likelihood

Parameters

• log10_fgw – log base 10 of GW frequency

• log10_mc – log base 10 of chirp mass

• max_toa – Latest TOA in any pulsar in the array

Returns

True if source already merged, False otherwise

QuickCW.QuickCorrectionUtils.correct_extrinsic(sample, x0)

correct extrinsic parameters for phases and cosines

Parameters

• sample – Array with parameters

• x0 – CWInfo object

Return sample

Corrected sample array

QuickCW.QuickCorrectionUtils.correct_extrinsic_array(samples, x0)

correct extrinsic parameters for phases and cosines

same as correct_extrinsic but handles mutiple samples

Parameters

• samples – 2D array with parameters for multiple samples

• x0 – CWInfo object

Return samples

Corrected samples array

QuickCW.QuickCorrectionUtils.correct_intrinsic(sample, x0, freq_bounds, cut_par_ids, cut_lows, cut_highs)

correct intrinsic parameters for phases and cosines

Parameters

• sample – Array with parameters

• x0 – CWInfo object

• freq_bounds – Lower and upper prior bounds of GW frequency

• cut_par_ids – Indices of parameters needing extra check

• cut_lows – Lower bounds of parameters needing extra check

• cut_highs – Upper bounds of parameters needing extra check

Return sample

Corrected sample array

QuickCW.QuickCorrectionUtils.reflect_cosines(cos_in, angle_in, rotfac=3.141592653589793, modfac=6.283185307179586)

helper to reflect cosines of coordinates around poles to get them between -1 and 1,

which requires also rotating the signal by rotfac each time, then mod the angle by modfac

Parameters

• cos_in – cosine of polar angle

• angle_in – azimuthal angle

• rotfac – angle to get to the other side of pole [np.pi]

• modfac – angle to change azimuthal angle by before getting back to the same point [2*np.pi]

Return cos_in

Reflected cosine of polar angle

Return angle_in

Reflected azimuthal angle

QuickCW.QuickCorrectionUtils.reflect_cosines_array(cos_ins, angle_ins, rotfac=3.141592653589793, modfac=6.283185307179586)

helper to reflect cosines of coordinates around poles to get them between -1 and 1,

which requires also rotating the signal by rotfac each time, then mod the angle by modfac same as reflect_cosines but handles arrays

Parameters

• cos_in – cosine of polar angle

• angle_in – azimuthal angle

• rotfac – angle to get to the other side of pole [np.pi]

• modfac – angle to change azimuthal angle by before getting back to the same point [2*np.pi]

Return cos_in

Reflected cosine of polar angle

Return angle_in

Reflected azimuthal angle

QuickCW.QuickCorrectionUtils.reflect_into_range(x, x_low, x_high)

reflect an arbitrary parameter into a nominal range

Parameters

• x – Value of parameter to reflect in range

• x_low – Lower bound of parameter

• x_high – Upper bound of prameter

Return res

Value of parameter reflected into range

QuickCW.QuickFisherHelpers module

C 2021 Bence Becsy MCMC for CW fast likelihood (w/ Neil Cornish and Matthew Digman) Helpers to get fisher matrices

QuickCW.QuickFisherHelpers.fisher_rn_mm_pp_diagonal_helper(params, x0_swap, FLI_swap, flm, par_names, epsilon_gammas, epsilon_log10_As, Npsr, get_intrinsic_diag=True, start_safe=False, get_gwb=True)

helper to get the mm and pp values needed to calculate the diagonal fisher eigenvectors for the red noise parameters

Parameters

• params –

  –

• x0_swap – CWInfo object

• FLI_swap – FastLikeInfo object

• flm – FastLikeMaster object

• par_names – List of parameter names

• epsilon_gammas – Perturbation values for gamma parameters

• epsilon_log10_As – Perturbation values for log10 amplitudes

• Npsr – Number of pulsars

• get_intrinsic_diag –

  –

• start_safe –

  –

• get_gwb –

  –

Return pp1s

–

Return mm1s

–

Return nn1s

–

Return helper_tuple0

–

Return pps_gwb

–

Return mms_gwb

–

Return nns_gwb

–

QuickCW.QuickFisherHelpers.fisher_synthetic_FLI_helper(helper_tuple_RR, x0_swap, FLI_swap, params_old, dist_mode=False, phase_mode=False)

helper to construct synthetic likelihoods for each pulsar from input MMs and NNs one by one

Parameters

• helper_tuple_RR –

  –

• x0_swap – CWInfo object

• FLI_swap – FastLikelihoodInfo object

• params_old – Old parameters to set FLI back to once done

• dist_mode – Specify if pulsar distances are being perturbed [False]

• phase_mode – Specify if phases are being perturbed [False]

Return rrs

–

QuickCW.QuickFisherHelpers.get_FLI_mem(FLI_swap)

store everything needed to reset FLI for non-red noise updates

Parameters

FLI_swap – FastLikelihoodInfo object

Return MMs0

list of M matrices

Return NN0

list of N vectors

Return resres_array0

array with resres values

Return logdet_array0

array with logdet contributions

Return logdet_base_old

base value of logdet

QuickCW.QuickFisherHelpers.get_fisher_diagonal(samples_fisher, par_names, x0_swap, flm, FLI_swap, get_intrinsic_diag=True, start_safe=False)

get the diagonal elements of the fisher matrix for all parameters

Parameters

• samples_fisher –

  –

• par_names – List of parameter names

• x0_swap – CWInfo object

• flm – FastLikeMaster object

• FLI_swap – FastLikeInfo object

• get_intrinsic_diag –

  –

• start_safe –

  –

Return 1/np.sqrt(fisher_diag)

–

Return pp2s

–

Return mm2s

–

Return nn2s

–

Return epsilons

–

Return helper_tuple0

–

Return pps

–

Return mms

–

Return nns

–

QuickCW.QuickFisherHelpers.get_fisher_eigenvectors(params, par_names, par_names_to_perturb, pta, epsilon=0.0001)

get fisher eigenvectors for a generic set of parameters the slow way

Parameters

• params – Parameter values where fisher is to be calculated

• par_names – List of parameter names

• par_names_to_perturb – Subset of par_names for which we want to calculate the fisher

• pta – enterprise PTA object

• epsilon – Perturbation values [1e-4]

Returns

Matrix with fisher eigenvectors

QuickCW.QuickFisherHelpers.get_fisher_eigenvectors_common(params, x0_swap, FLI_swap, diagonal_data, epsilon=0.0001, default_all=False)

update just the 4x4 block of common eigenvectors

Parameters

• params – Parameters where the fisher is to be calculated

• x0_swap – CWInfo object

• FLI_swap – FastLikeInfo object

• diagonal_data –

  –

• epsilon – Perturbation values [1e-4]

• default_all –

  –

Returns

Matrix with fisher eigenvectors

QuickCW.QuickFisherHelpers.get_fisher_rn_block_eigenvectors(params, par_names, x0_swap, flm, FLI_swap, diagonal_data_loc)

get the diagonal elements of the fisher matrix with the needed local stabilizations

Parameters

• params –

  –

• par_names – List of parameter names

• x0_swap – CWInfo object

• flm – FastLikeMaster object

• FLI_swap – FastLikeInfo object

• diagonal_data_loc –

  –

Return eig_rn

Matrix with RN eigenvectors

QuickCW.QuickFisherHelpers.get_fishers(samples, par_names, x0_swap, flm, FLI_swap, get_diag=True, get_common=True, get_rn_block=True, get_intrinsic_diag=True, start_safe=False)

get all the red noise eigenvectors in a block, and if get_diag is True also get all the diagonal fisher matrix elements

Parameters

• samples – Array of samples

• par_names – List of parameter names

• x0_swap – CWInfo object

• flm – FastLikeMaster object

• FLI_swap – FastLikeInfo object

• get_diag –

  –

• get_common –

  –

• get_rn_block –

  –

• get_intrinsic_diag –

  –

• start_safe –

  –

Return eig_rn

Matrix with RN eigenvectors

Return fisher_diag

Diagonal fisher

Return eig_common

Matrix of common parameter fishers

QuickCW.QuickFisherHelpers.params_perturb_helper(params, x0_swap, FLI_swap, flm, par_names, idxs_targ, epsilons, dist_mode=False, phase_mode=False, mask=None)

helper to perturb the specified parameters by factors of epsilon en masse

Parameters

• params – array of parameter values

• x0_swap – CWInfo object

• FLI_swap – FastLikelihoodInfo object

• flm – FastLikeMaster object

• par_names – List of parameter names

• idxs_targ – Indices of parameters to be perturbed

• epsilon – Amount of perturbation

• dist_mode – Specify if pulsar distances are being perturbed [False]

• phase_mode – Specify if phases are being perturbed [False]

• mask – Mask to specify which pulsars need to be updated; if None, all pulsars are updated [None]

Return paramsPP

Perturbed parameters

Return FLI_memp

FastLikeInfo object for the perturbed parameters

QuickCW.QuickFisherHelpers.safe_reset_swap(FLI_swap, x0_swap, params_old, FLI_mem0)

safely reset everything back to the initial values as input for self consistency in future calculations

Parameters

• FLI_swap – FastLikeInfo object to be reset

• x0_swap – CWInfo object

• params_old – Parameters to which we want to reset

• FLI_mem0 – Parts of FLI object saved to memory

QuickCW.QuickMCMCUtils module

C 2021 Bence Becsy MCMC for CW fast likelihood (w/ Neil Cornish and Matthew Digman) Helpers for MCMC; extrinsic blocks and parallel tempering

class QuickCW.QuickMCMCUtils.ChainParams(T_max: float, n_chain: int, n_block_status_update: int, n_int_block: int = 1000, n_update_fisher: int = 100000, save_every_n: int = 10000, fisher_eig_downsample: int = 10, T_ladder: ~typing.Optional[list] = None, includeCW: bool = True, prior_recovery: bool = False, verbosity: int = 1, freq_bounds: ~numpy.ndarray = array([   nan, 1.e-07]), gwb_comps: int = 14, cos_gwtheta_bounds: ~numpy.ndarray = array([-1, 1]), gwphi_bounds: ~numpy.ndarray = array([0., 6.28318531]), de_history_size: int = 5000, thin_de: int = 10000, log_fishers: bool = False, log_mean_likelihoods: bool = True, savefile: ~typing.Optional[str] = None, thin: int = 100, samples_precision: type = <class 'numpy.float32'>, save_first_n_chains: int = 1, prior_draw_prob: float = 0.1, de_prob: float = 0.6, fisher_prob: float = 0.3, rn_emp_dist_file: ~typing.Optional[str] = None, dist_jump_weight: float = 0.2, rn_jump_weight: float = 0.3, gwb_jump_weight: float = 0.1, common_jump_weight: float = 0.2, all_jump_weight: float = 0.2, fix_rn: bool = False, zero_rn: bool = False, fix_gwb: bool = False, zero_gwb: bool = False)

Bases: object

store basic parameters the govern the evolution of the mcmc chain

Parameters

• T_max – Maximum temperature of PT ladder

• n_chain – Number of PT chains

• n_block_status_update – Number of blocks between status updates

• n_int_block – Number of iterations in a block [1_000]

• n_update_fisher – Number of iterations between Fisher updates [100_000]

• save_every_n – Number of iterations between saving intermediate results (needs to be intiger multiple of n_int_block) [10_000]

• fisher_eig_downsample – Multiplier for how much less to do more expensive updates to fisher eigendirections for red noise and common parameters compared to diagonal elements [10]

• T_ladder – Temperature ladder; if None, geometrically spaced ladder is made with n_chain chains reaching T_max [None]

• includeCW – If False, we are not including the CW in the likelihood (good for testing) [True]

• prior_recovery – If True, likelihood is set to a constant (good for testing the prior recovery of the MCMC) [False]

• verbosity – Parameter indicating how much info to print (higher value means more prints) [1]

• freq_bounds – Lower and upper prior bounds on the GW frequency of the CW; np.nan lower bound is automatically turned into one over the observation time [[np.nan, 1.e-07]]

• gwb_comps – Number of frequency components to model in the GWB [14]

• cos_gwtheta_bounds – Prior bounds on the cosine of the GW theta sky location parameter (useful e.g. for targeted searches) [[-1,1]]

• gwphi_bounds – Prior bounds on the the GW phi sky location parameter (useful e.g. for targeted searches) [[0,2*np.pi]]

• de_history_size – Size of the differential evolution buffer

• thin_de – How much to thin samples for the DE buffer

• log_fishers –

  –

• log_mean_likelihoods –

  –

• savefile – File name to save the results to, if None, no results are saved [None]

• thin – How much to thin the samples by for saving [100]

• samples_precision – Precision to use for the saved samples [np.single]

• save_first_n_chains – Number of PT chains to save [1]

• prior_draw_prob – Probability of prior draws [0.1]

• de_prob – Probability of DE jumps [0.6]

• fisher_prob – Probability of fisher updates [0.3]

• rn_emp_dist_file – Filename with empirical distribution to use for per psr RN, if None, do not do empirical distribution jumps [None]

• dist_jump_weight – Weight if jumps changing pulsar distances [0.2]

• rn_jump_weight – Weight of jumps changing RN parameters [0.3]

• gwb_jump_weight – Weight of jumps changing GWB parameters [0.1]

• common_jump_weight – Weight of jumps changing common CW shape parameters (sky location, frequency, chirp mass) [0.2]

• all_jump_weight – Weight of jumps changing all parameters [0.2]

• fix_rn – If True, we fix per psr RN parameters to the value it starts at [False]

• zero_rn – If True, we fix per psr RN amplitude to a very low value effectively turning it off [False]

• fix_gwb – If True, we fix GWB parameters to the value it starts at [False]

• zero_gwb – If True, we fix GWB amplitude to a very low value effectively turning it off [False]

class QuickCW.QuickMCMCUtils.MCMCChain(chain_params, psrs, pta, max_toa, noisedict, ti)

Bases: object

store the miscellaneous objects needed to manage the mcmc chain

Parameters

• chain_params – ChainParams object

• psrs – List of enterprise pulsar objects

• pta – enterprise PTA object

• max_toa – Latest TOA in any pulsar in the array

• noisedict – Noise dictionary

• ti – Time after initialization got from time.perf_counter()

advance_N_blocks(N_blocks)

advance the state of the MCMC system by N_blocks of size

advance_block()

advance the state of the mcmc chain by 1 entire block, updating fisher matrices and differential evolution as necessary

do_status_update(itrn, N_blocks)

print a status update

output_and_wrap_state(itrn, N_blocks, do_output=True)

wrap the samples around to the first element and save the old ones to the hdf5 file

update_de_history(itrn)

update de history array

update_fishers(itrn)

handle fisher matrix update logic

update_fishers_partial(itrn1, itrn2)

handle fisher matrix update logic, itrn1 and itrn2 must be ranges over which no changes in the total set of intrinsic parameters occur so that we can skip the intrinsic update

validate_consistent(itrb, full_validate=False)

QuickCW.QuickMCMCUtils.add_rn_eig_starting_point(samples, par_names, x0_swap, flm, FLI_swap, chain_params, Npsr, FPI)

add a fisher eig jump to the starting point of each chain based only on the fisher matrix at the first point

Parameters

• samples – Samples to perturb

• par_names – List of parameter names

• x0_swap – CWInfo object

• flm – FastLikeMaster object

• FLI_swap – FastLikeInfo object

• chain_params – ChainParams object

• Npsr – Number of pulsars

• FPI – FastPriorInfo object

Return samples

Perturbed samples

QuickCW.QuickMCMCUtils.do_extrinsic_block(n_chain, samples, itrb, Ts, x0s, FLIs, FPI, n_par_tot, log_likelihood, n_int_block, fisher_diag, a_yes, a_no)

do blocks of just the extrinsic parameters, which should be very fast

Parameters

• n_chain – Number of PT chains

• samples – Array holding posterior samples

• itrb – Index within saved values (as opposed to block index itri or overall index itrn)

• Ts – List of PT temperatures

• x0s – List of CWInfo objects

• FLIs – List of FastLikeInfo objects

• FPI – FastPriorInfo object

• n_par_tot – Number of total parameters

• log_likelihood – Array holding log likelihood values

• n_int_block – Number of iterations to be done in a block

• fisher_diag – Diagonal fisher

• a_yes – Array to hold number of accepted steps

• a_no – Array to hold number of rejected steps

QuickCW.QuickMCMCUtils.do_pt_swap(n_chain, samples, itrb, Ts, a_yes, a_no, x0s, FLIs, log_likelihood, fisher_diag)

do the parallel tempering swap

Parameters

• n_chain – Number of PT chains

• samples – Array holding posterior samples

• itrb – Index within saved values (as opposed to block index itri or overall index itrn)

• Ts – List of PT temperatures

• a_yes – Array to hold number of accepted steps

• a_no – Array to hold number of rejected steps

• x0s – List of CWInfo objects

• FLIs – List of FastLikeInfo objects

• log_likelihood – Array holding log likelihood values

• fisher_diag – Diagonal fisher

QuickCW.QuickMCMCUtils.get_param_names(pta)

get the name Lists for various parameters

Parameters

pta – enterprise PTA object

Return par_names

List of parameter names

Return par_names_cw

List of parameter names describing the CW signal

Return par_names_cw_int

List of parameter names which are projection parameters (previously called extrinsic parameters)

Return par_names_cw_ext

List of parameter names which are shape parameters (previously called intrinsic parameters)

Return par_names_noise

List of noise parameter names

QuickCW.QuickMCMCUtils.initialize_de_buffer(sample0, n_par_tot, par_names, chain_params, x0_swap, FPI, eig_rn)

set up differential evolution

Parameters

• sample0 – Parameter values to start the RN from

• n_par_tot – Number of total parameters

• par_names – List of parameter names

• chain_params – ChainParams object

• x0_swap – CWInfo object

• FPI – FastPriorInfo object

• eig_rn – RN eigenvectors

Return de_history

Array holding initial differenctial evolution buffer

QuickCW.QuickMCMCUtils.initialize_sample_helper(chain_params, n_par_tot, Npsr, max_toa, par_names, par_names_cw_ext, par_names_cw_int, FPI, pta, noisedict, rn_emp_dist)

initialize starting samples for each chain to a random point

Parameters

• chain_params – ChainParams object

• n_par_tot – Total number of parameters

• Npsr – Number of pulsars

• max_toa – Latest TOA in any pulsar in the array

• par_names – List of parameter names

• par_names_cw_ext – List of parameter names which are projection parameters (previously called extrinsic parameters)

• par_names_cw_int – List of parameter names which are shape parameters (previously called intrinsic parameters)

• FPI – FastPriorInfo object

• pta – enterprise PTA object

• noisedict – Noise dictionary

• rn_emp_dist – RN empirical distributions

Return samples

Array to hold posterior samples initialized for the first sample

QuickCW.QuickMCMCUtils.uniform(low=0.0, high=1.0, size=None)

Draw samples from a uniform distribution.

Samples are uniformly distributed over the half-open interval [low, high) (includes low, but excludes high). In other words, any value within the given interval is equally likely to be drawn by uniform.

Note

New code should use the ~numpy.random.Generator.uniform method of a ~numpy.random.Generator instance instead; please see the random-quick-start.

Parameters

lowfloat or array_like of floats, optional

Lower boundary of the output interval. All values generated will be greater than or equal to low. The default value is 0.

highfloat or array_like of floats

Upper boundary of the output interval. All values generated will be less than or equal to high. The high limit may be included in the returned array of floats due to floating-point rounding in the equation low + (high-low) * random_sample(). The default value is 1.0.

sizeint or tuple of ints, optional

Output shape. If the given shape is, e.g., (m, n, k), then m * n * k samples are drawn. If size is None (default), a single value is returned if low and high are both scalars. Otherwise, np.broadcast(low, high).size samples are drawn.

Returns

outndarray or scalar

Drawn samples from the parameterized uniform distribution.

See Also

randint : Discrete uniform distribution, yielding integers. random_integers : Discrete uniform distribution over the closed

interval [low, high].

random_sample : Floats uniformly distributed over [0, 1). random : Alias for random_sample. rand : Convenience function that accepts dimensions as input, e.g.,

rand(2,2) would generate a 2-by-2 array of floats, uniformly distributed over [0, 1).

random.Generator.uniform: which should be used for new code.

Notes

The probability density function of the uniform distribution is

\[p(x) = \frac{1}{b - a}\]

anywhere within the interval [a, b), and zero elsewhere.

When high == low, values of low will be returned. If high < low, the results are officially undefined and may eventually raise an error, i.e. do not rely on this function to behave when passed arguments satisfying that inequality condition. The high limit may be included in the returned array of floats due to floating-point rounding in the equation low + (high-low) * random_sample(). For example:

>>> x = np.float32(5*0.99999999)
>>> x
5.0

Examples

Draw samples from the distribution:

>>> s = np.random.uniform(-1,0,1000)

All values are within the given interval:

>>> np.all(s >= -1)
True
>>> np.all(s < 0)
True

Display the histogram of the samples, along with the probability density function:

>>> import matplotlib.pyplot as plt
>>> count, bins, ignored = plt.hist(s, 15, density=True)
>>> plt.plot(bins, np.ones_like(bins), linewidth=2, color='r')
>>> plt.show()

QuickCW.QuickMTHelpers module

C 2021 Bence Becsy MCMC for CW fast likelihood (w/ Neil Cornish and Matthew Digman)

QuickCW.QuickMTHelpers.add_rn_eig_jump(scale_eig0, scale_eig1, new_point, rn_base, idx_rn, Npsr, all_eigs=False)

add a fisher eigenvalue jump to the red noise parameters in place

Parameters

• scale_eig0 – Amount to scale jump in gamma values by

• scale_eig1 – Amount to scale in log10_A values by

• new_point – Parameter values to add RN jump to

• rn_base – RN values to jump from (usually justa slice of new_point)

• idx_rn – Indices of new_point containing RN parameters

• Npsr – Number of pulsars

• all_eigs – If True, perturb all pulsars’ RN, if False, pick randomly [False]

Return new_point

Perturbed parameter values

QuickCW.QuickMTHelpers.do_intrinsic_update_mt(mcc, itrb)

do the intrinsic update using the multiple try mcmc algorithm

Parameters

• mcc – MCMCChain onject

• itrb – Index within saved values (as opposed to block index itri or overall index itrn)

Return mcc.FLI_swap

FastLikeInfo object

QuickCW.QuickMTHelpers.do_mt_step(mcc, j, itrb, new_point, samples_current, FLI_mem_save, recompute_rn, log_proposal_ratio)

compute the multiple tries and chose a sample

Parameters

• mcc – MCMCChain onject

• j – Index of PT chain

• itrb – Index within saved values (as opposed to block index itri or overall index itrn)

• new_point – Proposed new point (with new shape parameters)

• samples_current – Current point in parameter space

• FLI_mem_save – Parts of FLI object saved to memory

• recompute_rn – If True, recompute everything needed to go to new RN parameters

• log_proposal_ratio – Log of the proposal ratio needed to calculate acceptance probability

Return log_acc_ratio

Log of acceptance probability

Return chosen_trial

Index of chosen trial

Return sample_choose

Parameters of the chosen trial

Return log_Ls[chosen_trial]

Log likelihood of the chosen trial

QuickCW.QuickMTHelpers.get_mt_weights(x0_extras, FLI_use, Ts, log_posterior_old, tries, log_prior_news)

Helper function to quickly return multiple tries and their likelihoods fo MTMCMC

Parameters

• x0_extras – List of extra CWInfo objects for parallelizing multiple try

• FLI_use – FastLikeInfo object

• Ts – List of PT temperatures

• log_posterior_old – Log posterior at old parameters

• tries – Parameters at a set of multiple tries for which we want to calculate the weights

• log_prior_news – Log prior values at propose new points

Return mt_weights

Multiple try weights

Return log_Ls

Log likelihoods

Return log_mt_norm_shift

Amount to shift the multiple try weights (helps with using floating point precision efficiently)

QuickCW.QuickMTHelpers.get_ref_mt_weights(x0_extras, FLI_use, Ts, log_posterior_old, chosen_trial, ref_tries, log_prior_refs)

Helper function to quickly return multiple tries and their likelihoods fo MTMCMC

Parameters

• x0_extras – List of extra CWInfo objects for parallelizing multiple try

• FLI_use – FastLikeInfo object

• Ts – List of PT temperatures

• log_posterior_old – Log posterior at old parameters

• chosen_trial – Index of chosen trial

• ref_tries – Parameters at a set of reference multiple tries for which we want to calculate the weights

• log_prior_refs – Log prior values at reference points

Return ref_mt_weights

Reference point multiply try weights

Return log_ref_mt_norm_shift

Amount to shift the reference point multiple try weights (helps with using floating point precision efficiently)

QuickCW.QuickMTHelpers.set_params(sample_set, jumps, fisher_mask, random_draws_from_prior, x0)

assign parameters to tries for multiple try mcmc

Parameters

• sample_set – Samples to start from

• jumps – Precaluclated fisher jumps to use

• fisher_mask – Mask determining which projection parameters to do fisher jump vs prior draw in

• random_draws_from_prior – Precalculated prior draws to use

• x0 – CWInfo object

Return ref_tries

2D array holding samples at multiple trials

QuickCW.QuickMTHelpers.uniform(low=0.0, high=1.0, size=None)

Draw samples from a uniform distribution.

Samples are uniformly distributed over the half-open interval [low, high) (includes low, but excludes high). In other words, any value within the given interval is equally likely to be drawn by uniform.

Note

New code should use the ~numpy.random.Generator.uniform method of a ~numpy.random.Generator instance instead; please see the random-quick-start.

Parameters

lowfloat or array_like of floats, optional

Lower boundary of the output interval. All values generated will be greater than or equal to low. The default value is 0.

highfloat or array_like of floats

Upper boundary of the output interval. All values generated will be less than or equal to high. The high limit may be included in the returned array of floats due to floating-point rounding in the equation low + (high-low) * random_sample(). The default value is 1.0.

sizeint or tuple of ints, optional

Output shape. If the given shape is, e.g., (m, n, k), then m * n * k samples are drawn. If size is None (default), a single value is returned if low and high are both scalars. Otherwise, np.broadcast(low, high).size samples are drawn.

Returns

outndarray or scalar

Drawn samples from the parameterized uniform distribution.

See Also

randint : Discrete uniform distribution, yielding integers. random_integers : Discrete uniform distribution over the closed

interval [low, high].

random_sample : Floats uniformly distributed over [0, 1). random : Alias for random_sample. rand : Convenience function that accepts dimensions as input, e.g.,

rand(2,2) would generate a 2-by-2 array of floats, uniformly distributed over [0, 1).

random.Generator.uniform: which should be used for new code.

Notes

The probability density function of the uniform distribution is

\[p(x) = \frac{1}{b - a}\]

anywhere within the interval [a, b), and zero elsewhere.

When high == low, values of low will be returned. If high < low, the results are officially undefined and may eventually raise an error, i.e. do not rely on this function to behave when passed arguments satisfying that inequality condition. The high limit may be included in the returned array of floats due to floating-point rounding in the equation low + (high-low) * random_sample(). For example:

>>> x = np.float32(5*0.99999999)
>>> x
5.0

Examples

Draw samples from the distribution:

>>> s = np.random.uniform(-1,0,1000)

All values are within the given interval:

>>> np.all(s >= -1)
True
>>> np.all(s < 0)
True

Display the histogram of the samples, along with the probability density function:

>>> import matplotlib.pyplot as plt
>>> count, bins, ignored = plt.hist(s, 15, density=True)
>>> plt.plot(bins, np.ones_like(bins), linewidth=2, color='r')
>>> plt.show()

QuickCW.const_mcmc module

store mcmc constants for global reference

QuickCW.lapack_wrappers module

C 2021 Matthew Digman various jit compatible interfaces to cython lapack functions

QuickCW.lapack_wrappers.solve_triangular(x, y, lower_a=True, trans_a=True, unitdiag=False, overwrite_b=False)

solve x*B=y

Parameters

• x – triangular matrix (must be either type of contiguous)

• y – vector (must be fortran ordered)

Return B

Solution to x*B=y

QuickCW.test_fisher_convergence module