API

panco2.panco2

class panco2.PressureProfileFitter(sz_map_file, hdu_data, hdu_rms, z, M_500, coords_center=None, map_size=None, cosmo=FlatLambdaCDM(H0=70.0 km / (Mpc s), Om0=0.3, Tcmb0=0.0 K, Neff=3.04, m_nu=None, Ob0=None))

The main class of panco2, that manages data loading, model definition, and MCMC sampling.

Parameters

sz_map_file (str) – Path to a FITS file containing the SZ map and noise RMS
hdu_data (int) – Index of the FITS extension in which the SZ map is stored
hdu_rms (int) – Index of the FITS extension in which the RMS map is stored
z (float) – Cluster’s redshift
M_500 (float) – A guess of the cluster’s mass [Msun]. This is used to build the starting point of the MCMC.
coords_center (astropy.coordinates.SkyCoord, optional) – Coordinate to consider as the center of the map. If None (default), the center of the FITS map is used.
map_size (float, optional) – The size of the map to be considered [arcmin]. If None (default), the entire FITS map is used.
cosmo (astropy.cosmology.Cosmology, optional) – The cosmology to assume for distance computations. Defaults to flat LCDM with h=0.7, Om0=0.3.

add_covmat(covmat=None, inv_covmat=None)

Adds in a (inverse) noise covariance matrix to be used in the log-likelihood computation, in order to account for spatial correlations in the noise. The panco2.noise_covariance module offers ways to compute such matrices from different types of inputs: noise power spectrum, noise maps, etc – see documentation.

Parameters

covmat (ndarray, optional) – Noise covariance matrix [data units**2]. shape=(n_pix_map**2, n_pix_map**2).
inv_covmat (ndarray, optional) – Inverse noise covariance matrix [data units**2]. shape=(n_pix_map**2, n_pix_map**2).

Raises

Exception – If neither the covariance matrix or inverted covariance matrix is provided.

Notes

If inv_covmat is provided, it is taken as-is for the log-likelihood computation. If it is not, and covmat is given, the code will try to compute the Moore-Penrose pseudo-inverse matrix using scipy.linalg.pinv, and check that the matrix product C @ C-1 is close to identity.

add_filtering(beam_fwhm=0.0, pad=0, ell=None, tf=None)

Initialize convolution by a beam and transfer function in the model computation.

Parameters

beam_fwhm (float, optional) – The FWHM of your gaussian beam [arcsec]. Can be set to 0 if you don’t want beam convolution (e.g. if the beam is already in the transfer function).
tf (array) – Filtering measurement.
ell (array | tuple of arrays) – Multipole moments at which the filtering was measured. Can be a tuple of ellx and elly, see Notes.
pad (int) – Padding to be added to the sides of the map before convolution [pixels].

Notes

The code can deal with 1D or 2D transfer functions, depending on the inputs:
- For a 1D transfer function (isotropic filtering), ell and tf should be 1D arrays of same shape
- For a 2D transfer function, ell should be a tuple of two 1D arrays (ell_x and ell_y), and tf should be a 2D array with shape as the outer product (ell_x x ell_y)

add_integ_Y(Y, dY, r)

Add a constraint on the integrated Compton parameter to the fit.

Parameters

Y (float) – Integrated Compton parameter value [kpc2]
dY (float) – 1 sigma uncertainty on Y [kpc2]
r (float) – Radius within which the Compton parameter was integrated [kpc]

add_mask(mask)

Add a mask to discard part of the data in the model fitting. The panco2.masks module offers ways to create this mask for some standard usecases – see documentation.

Parameters

mask (np.ndarray) – Boolean mask, where True means the pixel is masked. shape=(n_pix_map, n_pix_map).

Raises

Exception – If the mask is not an array of the right shape.
Warning – If more than half the map pixels are masked.

add_point_sources(coords, beam_fwhm)

Add point sources to the model.

Parameters

coords (list of astropy.coordinates.SkyCoord) – The positions of each source in the sky
beam_fwhm (float [arcsec]) – The FWHM of the Gaussian model to be used for point sources

define_model(r_bins, zero_level=True)

Define options of the model to be used to fit the map.

Parameters

r_bins (array) – List of radial bins [kpc].
zero_level (bool, optional) – If True (default), a constant zero level in the map is also fitted.

define_priors(P_bins=None, conv=<scipy.stats._distn_infrastructure.rv_continuous_frozen object>, zero=<scipy.stats._distn_infrastructure.rv_continuous_frozen object>, ps_fluxes=[])

Parameters

P_bins (ss.distributions or list of ss.distributions) – Priors on the pressure bins in binned profile mode. If only one distribution, all bins will have the same prior. If a list, the length should be the number of bins. [keV cm-3]
conv (ss.distributions instance) – Prior on the conversion coefficient. Defaults to N(0, 1). [map units]
zero (ss.distributions instance) – Prior on the zero level of the map. Defaults to N(0, 1). [map units]
ps_fluxes (list of ss.distributions) – Priors on the flux of each point source [map units]

Raises

Exception – If the list of priors don’t have the same size as the list of parameters (for the point source fluxes or the pressure bins).

Notes

To create distributions, use the scipy.stats module.

Examples

Normal distribution with mean 0 and spread 1:

import scipy.stats as ss

prior_on_parameter = ss.norm(0.0, 1.0)
Uniform distribution between 0 and 1:

import scipy.stats as ss

prior_on_parameter = ss.uniform(0.0, 1.0)

dump_to_file(file_name)

Writes a PressureProfileFitter object as a dill serialized file.

Parameters: file_name (str) – Name of the file to write

panco2.results

panco2.results.get_best_fit(chains_clean, ppf)

Get best-fit from Markov chains.

Parameters

chains_clean (pd.DataFrame) – Dataframe of the chains, produced by clean_chains()
ppf (panco2.PressureProfileFitter) – The main PressureProfileFitter object created for the fit. Used to compute the number of degrees of freedom.

Returns

int – The index of the best-fit vector in chains_clean
float – The value of reduced chi2 for the best-fit
dict – The best-fitting parameters

panco2.results.load_chains(out_chains_file, burn, discard, clip_percent=0, verbose=False)

Loads raw Markov chains, cleans them, and arranges them in a convenient data frame.

Parameters

out_chains_file (str) – Path to the npz file containing the raw chains produced by panco2.PressureProfileFitter.run_mcmc
burn (int) – Length to discard as burn-in.
discard (int) – Thinning length. See Notes.
clip_percent (float) – Percentile considered as extreme to weed-out bad chains. See Notes.
verbose (bool) – Wether or not you want results of the cleaning to be printed.

Returns

chains_clean – Cleaned chains in a DataFrame. Keys are the parameter names, as well as lnprior (log-prior value), lnlike (log-likelihood), chain (which walker the sample belongs to) and step (which step of the walk the sample was obtained at).

Return type

pd.DataFrame

Notes

Cleaning is done as follows:

Discard burn burn-in steps and keep one sample every discard steps
Find chains that are always out of the [clip_percent, 100 - clip_percent] interval for every parameter and discard them.

panco2.results.mcmc_corner_plot(chains_clean, per_chain=False, show_probs=True, filename=None, ppf=None)

Plots 1D and 2D posterior sampled by the MCMC.

Parameters

chains_clean (pd.DataFrame) – Markov chains, output of load_chains
per_chain (bool) – If True, plots one distribution per chain. This can make very cluttered plots.
show_probs (bool) – Wether or not the log likelihood and log prior chains should be shown along with the parameters
filename (str or None) – Filename to save the plot
ppf (PressureProfileFitter instance or None) – The main panco2.PressureProfileFitter instance. If provided, the marginal prior distribution will be shown for comparison with the posterior.

Return type

fig

panco2.results.mcmc_matrices_plot(chains_clean, ppf, filename=None)

Produces an overly complicated plot of the correlation and covariance matrices of Markov chains, where:

The lower corner (diagonal not included) shows the correlation coefficient between parameters;
The upper corner (diagonal included) shows the absolute covariance between parameters.

Parameters

chains_clean (pd.DataFrame) – Markov chains, output of load_chains
ppf (PressureProfileFitter instance) – The main panco2.PressureProfileFitter instance.
filename (str or None) – Filename to save the plot

Return type

fig, ax

panco2.results.mcmc_trace_plot(chains_clean, show_probs=True, filename=None)

Create MCMC trace plots, i.e. chains evolution with steps.

Parameters

chains_clean (pd.DataFrame) – Markov chains, output of load_chains
show_probs (bool) – Wether or not the log likelihood and log prior chains should be shown along with the parameters
filename (str or None) – Filename to save the plot

Return type

fig

panco2.results.plot_data_model_residuals(ppf, par_vec=None, par_dic=None, smooth=0.0, fig=None, axs=None, lims=None, cbar_label=None, cbar_fact=1.0, cbar_hztl=False, cmap='RdBu_r', mask_color='0.75', snr_contours=None, separate_ps_model=False, filename_plot=None, filename_maps=None)

Plot data, model, and residuals maps.

Parameters

ppf (PressureProfileFitter instance) – The main panco2.PressureProfileFitter instance.
par_vec (list or None) – Vector in the parameter space to use to compute the model map to be shown. Either one of par_vec or par_dic must be provided.
par_dic (dict or None) – Same as par_vec, but in a dictionary. Either one of par_vec or par_dic must be provided.
smooth (float [pixels]) – Width (sigma) of a gaussian kernel to be used to smooth the maps, for visual purposes.
fig (plt.Figure or None) – If provided, an existing figure can be used
axs (list of plt.Axis or None) – If provided, existing axes can be used. The length of the list must be consistent with the number of maps to show
lims (tuple or None) – Limits of the color maps, in data units
cbar_label (str or None) – Label for the colorbar
cbar_fact (float) – A factor by which all maps are to be multiplied before plotting. Useful for very small units like Compton-y or Jy/beam
cbar_hztl (bool) – Makes the colorbar horizontal instead of vertical.
cmap (str or mpl.colors.Colorbar) – The color map to use. Always make pretty plots!
mask_color (matplotlib color) – Color to fill the masked pixels with.
snr_contours (None or array-like) – Signal-to-noise levels to be overplotted as contours, in numbers of sigma.
separate_ps_model (bool) – If your model fits both SZ and point sources, makes the model and residuals for SZ/PS/SZ+PS (i.e. 7 total plots)
filename_plot (str or None) – Filename to save the plot
filename_maps (str or None) – Filename to save the maps in FITS format

Return type

fig, axs

Raises

Exception – If neither par_vec or par_dic are provided.

panco2.results.plot_data_model_residuals_1d(ppf, par_vec=None, par_dic=None, chains_clean=None, fig=None, ax=None, y_label=None, y_fact=1.0, x_log=False, plot_beam=True, filename=None)

Parameters

ppf (PressureProfileFitter instance) – The main panco2.PressureProfileFitter instance.
par_vec (list or None) – Vector in the parameter space to use to compute the model map to be shown. Either one of par_vec or par_dic must be provided.
par_dic (dict or None) – Same as par_vec, but in a dictionary. Either one of par_vec or par_dic must be provided.
chains_clean (pd.DataFrame) – Markov chains, output of load_chains
fig (plt.Figure or None) – If provided, an existing figure can be used
ax (plt.Axis or None) – If provided, existing axes can be used. The length of the list must be consistent with the number of maps to show
y_label (str or None) – Label for the y axis
y_fact (float) – A factor by which all profiles are to be multiplied before plotting. Useful for very small units like Compton-y or Jy/beam
x_log (bool) – Log-scale for the x axis
plot_beam (bool) – Overplot the beam profile.
filename (str or None) – Filename to save the plot

Return type

fig, ax

Raises

Exception – if neither par_dic, par_vec, or chains_clean are provided.

Notes

If par_dic or par_vec are provided, only one line will be plotted for the model (corresponding to one position in the parameter space). If chains_clean is provided, confidence intervals will be drawn, giving more information, but making the process slower.

panco2.results.plot_profile(chains_clean, ppf, r_range, P_compare=None, kwargs_compare={}, label=None, filename=None, **kwargs)

Plots the pressure profile recovered by PANCO2 from the Markov chains.

Parameters

chains_clean (pd.DataFrame) – Markov chains, output of load_chains
ppf (PressureProfileFitter instance) – The main panco2.PressureProfileFitter instance.
r_range (np.array [kpc]) – The radial range on which to show the profile.
P_compare (np.array [keV cm-3] or None) – A profile to compare panco2’s results with, such as the true profile when using simulated maps. Must be the same shape as r_range. If provided, a bottom panel will be added showing the ratio between results and P_compare.
kwargs_compare (dict) – keyword arguments to pass plt.plot for the compared profile (should include label, and any other style info you want)
label (str or None) – Label of the curve for legend purposes
filename (str or None) – Filename to save the plot
**kwargs (dict) – Other options to pass to plt.plot

Return type

fig, axs

panco2.noise_covariance

panco2.noise_covariance.covmat_from_noise_maps(noise_maps, method='lw')

Computes the pixel noise covariance matrix and its inverse from random noise realizations.

Parameters

noise_maps (ndarray) – The noise realizations, shape=(n_maps, nx, nx)
method (str, optional) – How to compute the covariance. Must be either “np”, in which case the covariance is the sample covariance computed via numpy.cov, or “lw” for the Ledoit-Wolf covariance estimated using sklearn.covariance.ledoit_wolf() (default).

Returns

ndarray – The covariance matrix, shape=(nx^2, nx^2)
ndarray – The inverted covariance matrix, shape=(nx^2, nx^2)

Raises

Exception – If method is not set to “np” or “lw”

panco2.noise_covariance.covmat_from_powspec(ell, C_ell, n_pix, pix_size, n_maps=1000, method='lw', return_maps=False)

Computes the pixel noise covariance matrix and its inverse from a noise power spectrum by generating many random maps with the input power spectrum.

Parameters

ell (ndarray) – Multipole numbers
C_ell (ndarray) – Power spectrum values
n_pix (int) – Number of pixels on each dimension of the map (i.e. the map is n_pix * n_pix)
pix_size (float) – Size of one pixel in the map [arcsec]
n_maps (int, optional) – Number of maps to be created. Defaults to 1000.
method (str, optional) – How to compute the covariance. Must be either “np”, in which case the covariance is the sample covariance computed via numpy.cov, or “lw” for the Ledoit-Wolf covariance estimated using sklearn.covariance.ledoit_wolf() (default).
return_maps (bool, optional) – If True, this function also returns the noise realizations used to compute the covariance, default is False.

Returns

ndarray – The covariance matrix, shape=(nx^2, nx^2)
ndarray – The inverted covariance matrix, shape=(nx^2, nx^2)
ndarray – If return_maps is True, the noise maps used to compute the covariance, shape=(n_maps, nx, nx)

panco2.noise_covariance.covmats_from_noise_map(noise_map, pix_size, n_maps=1000, method='lw', return_maps=False)

Computes the pixel noise covariance matrix and its inverse from a noise map by generating many random realizations with the same power spectrum.

Parameters

noise_map (ndarray) – The noise realizations, shape=(nx, nx)
pix_size (float) – The size of one pixel in the map [arcsec]
n_maps (int, optional) – Number of maps to be created. Defaults to 1000.
method (str, optional) – How to compute the covariance. Must be either “np”, in which case the covariance is the sample covariance computed via numpy.cov, or “lw” for the Ledoit-Wolf covariance estimated using sklearn.covariance.ledoit_wolf() (default).
return_maps (bool, optional) – If True, this function also returns the noise realizations used to compute the covariance, default is False.

Returns

ndarray – The covariance matrix, shape=(nx^2, nx^2)
ndarray – The inverted covariance matrix, shape=(nx^2, nx^2)
ndarray – If return_maps is True, the noise maps used to compute the covariance, shape=(n_maps, nx, nx)

panco2.masks

panco2.masks.cut_mask_from_fits(ppf, filename, hdu_mask)

Crops a mask from a FITS file.

Parameters

ppf (panco2.PressureProfileFitter instance) – The panco2 object to create a mask for. This is used to read in the WCS and map size.
filename (str) – Path to the FITS file containing the mask to be cropped.
hdu_mask (int) – Index of the extension of the FITS file containing the mask.

Returns

Boolean mask with the same shape as the data. Warning: this does not check the mask – if the FITS mask is 0 for masked pixels and 1 elsewhere, you’ll need to invert it to use it in panco2.PressureProfileFitter.add_mask().

Return type

np.ndarray

panco2.masks.mask_holes(ppf, coords, radius)

Creates a boolean mask by punching circular holes at given sky positions.

Parameters

ppf (panco2.PressureProfileFitter instance) – The panco2 object to create a mask for. This is used to read in the WCS and map size.
coords (List of astropy.coordinates.SkyCoord instances) – Coordinates of the centers of the holes.
radius (float) – Radius of each hole, in pixels.

Returns

Boolean mask with circular holes at given positions. The convention is the right one to use as input for panco2.PressureProfileFitter.add_mask(), i.e. the mask is True in the pixels to be masked and False elsewhere.

Return type

np.ndarray