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
See also
- classmethod load_from_file(file_name)
Load a PressureProfileFitter object from a dill serialized file.
- Parameters
file_name (str) – Name of the file dump
- Returns
The panco2.PressureProfileFitter object
- Return type
See also
- run_mcmc(n_chains, max_steps, n_threads, n_check=1000.0, max_delta_tau=0.01, min_autocorr_times=100, out_chains_file='./chains.npz', plot_convergence=None, progress=True, verbose_monitor=True)
Runs MCMC sampling of the posterior distribution.
- Parameters
n_chains (int) – Number of emcee walkers
max_steps (int) – Maximum number of steps in the Markov chains. The final number of points can be lower if convergence is accepted before max_steps is reached – see Notes.
n_threads (int) – Number of parallel threads to use
n_check (int) – Number of steps between two convergence checks
max_delta_tau (float) – Maximum relative difference of the autocorrelation length between two convergence checks to end sampling
min_autocorr_times (int) – Minimum ratio between the chains length and the autocorrelation length to end samlling
out_chains_file (str) – Path to a .npz file in which the chains will be stored. Default: “./chains.npz”.
plot_convergence (str or None) – Filename to save a plot of the autocorrelation function evolution and convergence test. If None (default), the plot is not produced.
progress (bool) – If True (default), displays a progressbar showing the sampling progression
verbose_monitor (bool) – If True (default), prints the result of each convergence check
- Returns
chains – Markov chains. Each key is a parameter, and the values are 2D arrays of shape (n_chains, n_steps).
- Return type
dict
Notes
The convergence check are performed every n_check steps. Convergence is accepted if:
In the last two checks, the mean autocorrelation time had a relative variation smaller than max_delta_tau compared to the previous step.
The chain is longer than min_autocorr_times times the current mean autocorrelation time.
- write_sim_map(par_vec, out_file, filter_noise=False, corr_noise=False)
Write a FITS map with given parameter values
- Parameters
par_vec (list or array) – Vector in the parameter space
out_file (str) – Path to a FITS file to which the map will be written. If the file already exists, it will be overwritten.
filter_noise (bool) – If True, convolve the noise realization by your filtering kernel. Defaults to False.
corr_noise (bool) – If True, the random noise realization will be drawn using the noise covariance matrix. Defaults to False.
Notes
The created FITS file contains the following extensions:
HDU 0: primary, contains header and no data.
HDU 1: TOTAL, contains the model map (SZ+PS+noise) and header.
HDU 2: SZ, contains the SZ model map and header.
HDU 3: PS, contains the PS model map and header.
HDU 4: NOISE, contains the noise map realization and header.
HDU 5: RMS, contains the noise RMS map and header.
All maps are cropped identically to the data used for the fit, and the headers are adjusted accordingly.
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