Python API documentation

pymaster contains three basic classes:

and a number of functions

pymaster also comes with a flat-sky version with most of the same functionality:

Many of these function accept or return sets of power spectra (arrays with one element per angular multipole) or bandpowers (binned versions of power spectra). In all cases, these are returned and provided as 2D arrays with shape [n_cls][nl], where n_cls is the number of power spectra and nl is either the number of multipoles or bandpowers. In all cases, n_cls should correspond with the spins of the two fields being correlated, and the ordering is as follows:

  • Two spin-0 fields: n_cls=1, [C_T1T2]
  • One spin-0 field and one spin>0 field: n_cls=2, [C_TE,C_TB]
  • Two spin>0 fields: n_cls=4, [C_E1E2,C_E1B2,C_E2B1,C_B1B2]

All sky maps accepted and returned by these functions are in the form of HEALPix maps exclusively with RING ordering.

class pymaster.field.NmtField(mask, maps, spin=None, templates=None, beam=None, purify_e=False, purify_b=False, n_iter_mask_purify=3, tol_pinv=1e-10, wcs=None, n_iter=3, lmax_sht=-1, masked_on_input=False, lite=False)

An NmtField object contains all the information describing the fields to correlate, including their observed maps, masks and contaminant templates.

Parameters:
  • mask – array containing a map corresponding to the field’s mask. Should be 1-dimensional for a HEALPix map or 2-dimensional for a map with rectangular pixelization.
  • maps – array containing the observed maps for this field. Should be at least 2-dimensional. The first dimension corresponds to the number of maps, which should be 1 for a spin-0 field and 2 otherwise. The other dimensions should be [npix] for HEALPix maps or [ny,nx] for maps with rectangular pixels. For a spin>0 field, the two maps to pass should be the usual Q/U Stokes parameters for polarization, or e1/e2 (gamma1/gamma2 etc.) in the case of cosmic shear. It is important to note that NaMaster uses the same polarization convention as HEALPix (i.e. with the x-coordinate growing with increasing colatitude theta). It is however more common for galaxy ellipticities to be provided using the IAU convention (i.e. x grows with declination). In this case, the sign of the e2/gamma2 map should be swapped before using it to create an NmtField. See more here . If None, this field will only contain a mask but no maps. The field can then be used to compute a mode-coupling matrix, for instance, but not actual power spectra.
  • spin – field’s spin. If None it will be set to 0 if there is a single map on input, and will default to 2 if there are 2 maps.
  • templates – array containing a set of contaminant templates for this field. This array should have shape [ntemp][nmap]…, where ntemp is the number of templates, nmap should be 1 for spin-0 fields and 2 otherwise. The other dimensions should be [npix] for HEALPix maps or [ny,nx] for maps with rectangular pixels. The best-fit contribution from each contaminant is automatically removed from the maps unless templates=None.
  • beam – spherical harmonic transform of the instrumental beam (assumed to be rotationally symmetric - i.e. no m dependence). If None, no beam will be corrected for. Otherwise, this array should have at least as many elements as the maximum multipole sampled by the maps + 1 (e.g. if a HEALPix map, it should contain 3*nside elements, corresponding to multipoles from 0 to 3*nside-1).
  • purify_e – use pure E-modes?
  • purify_b – use pure B-modes?
  • n_iter_mask_purify – number of iterations used to compute an accurate SHT of the mask when using E/B purification.
  • tol_pinv – when computing the pseudo-inverse of the contaminant covariance matrix, all eigenvalues below tol_pinv * max_eval will be treated as singular values, where max_eval is the largest eigenvalue. Only relevant if passing contaminant templates that are likely to be highly correlated.
  • wcs – a WCS object if using rectangular pixels (see http://docs.astropy.org/en/stable/wcs/index.html).
  • n_iter – number of iterations when computing a_lms.
  • lmax_sht – maximum multipole up to which map power spectra will be computed. If negative or zero, the maximum multipole given the map resolution will be used (e.g. 3 * nside - 1 for HEALPix maps).
  • masked_on_input – set to True if input maps and templates are already multiplied by the masks. Note that this is not advisable if you’re using purification.
  • lite – set to True if you want to only store the bare minimum necessary to run a standard pseudo-Cl with deprojection and purification, but you don’t care about deprojection bias. This will reduce the memory taken up by the resulting object.
get_alms()

Returns a 2D array ([nmap][nlm]) corresponding to the observed harmonic coefficients of this field.

Returns:2D array of alms
get_maps()

Returns a 2D array ([nmap][npix]) corresponding to the observed maps for this field. If the field was initialized with contaminant templates, the maps returned by this function have their best-fit contribution from these contaminants removed.

Returns:2D array of maps
get_templates()

Returns a 3D array ([ntemp][nmap][npix]) corresponding to the contaminant templates passed when initializing this field.

Returns:3D array of maps
class pymaster.field.NmtFieldFlat(lx, ly, mask, maps, spin=None, templates=None, beam=None, purify_e=False, purify_b=False, tol_pinv=1e-10, masked_on_input=False, lite=False)

An NmtFieldFlat object contains all the information describing the flat-sky fields to correlate, including their observed maps, masks and contaminant templates.

Parameters:
  • lx,ly (float) – size of the patch in the x and y directions (in radians)
  • mask – 2D array (nx,ny) containing a HEALPix map corresponding to the field’s mask.
  • maps – 2 2D arrays (nmaps,nx,ny) containing the observed maps for this field. The first dimension corresponds to the number of maps, which should be 1 for a spin-0 field and 2 otherwise. If None, this field will only contain a mask but no maps. The field can then be used to compute a mode-coupling matrix, for instance, but not actual power spectra.
  • spin – field’s spin. If None it will be set to 0 if there is a single map on input, and will default to 2 if there are 2 maps.
  • templates – array of maps (ntemp,nmaps,nx,ny) containing a set of contaminant templates for this field. This array should have shape [ntemp][nmap][nx][ny], where ntemp is the number of templates, nmap should be 1 for spin-0 fields and 2 for spin-2 fields, and nx,ny define the patch. The best-fit contribution from each contaminant is automatically removed from the maps unless templates=None
  • beam – 2D array (2,nl) defining the FT of the instrumental beam (assumed to be rotationally symmetric). beam[0] should contain the values of l for which de beam is defined, with beam[1] containing the beam values. If None, no beam will be corrected for.
  • purify_e – use pure E-modes?
  • purify_b – use pure B-modes?
  • tol_pinv – when computing the pseudo-inverse of the contaminant covariance matrix, all eigenvalues below tol_pinv * max_eval will be treated as singular values, where max_eval is the largest eigenvalue. Only relevant if passing contaminant templates that are likely to be highly correlated.
  • masked_on_input – set to True if input maps and templates are already multiplied by the masks. Note that this is not advisable if you’re using purification.
  • lite – set to True if you want to only store the bare minimum necessary to run a standard pseudo-Cl with deprojection and purification, but you don’t care about deprojection bias. This will reduce the memory taken up by the resulting object.
get_maps()

Returns a 3D array ([nmap][ny][nx]) corresponding to the observed maps for this field. If the field was initialized with contaminant templates, the maps returned by this function have their best-fit contribution from these contaminants removed.

Returns:3D array of flat-sky maps
get_templates()

Returns a 4D array ([ntemp][nmap][ny][nx]) corresponding to the contaminant templates passed when initializing this field.

Returns:4D array of flat-sky maps
class pymaster.bins.NmtBin(nside=None, bpws=None, ells=None, weights=None, nlb=None, lmax=None, is_Dell=False, f_ell=None)

An NmtBin object defines the set of bandpowers used in the computation of the pseudo-Cl estimator. The definition of bandpowers is described in Section 3.6 of the scientific documentation. We provide several convenience constructors that cover a range of common use cases requiring fewer parameters (see NmtBin.from_nside_linear(), NmtBin.from_lmax_linear() and Nmt.from_edges()).

Parameters:
  • nside (int) – HEALPix nside resolution parameter of the maps you intend to correlate. The maximum multipole considered for bandpowers will be 3*nside-1, unless lmax is set.
  • ells (array-like) – array of integers corresponding to different multipoles
  • bpws (array-like) – array of integers that assign the multipoles in ells to different bandpowers. All negative values will be ignored.
  • weights (array-like) – array of floats corresponding to the weights associated to each multipole in ells. The sum of weights within each bandpower is normalized to 1.
  • nlb (int) – integer value corresponding to a constant bandpower width. I.e. the bandpowers will be defined as consecutive sets of nlb multipoles from l=2 to l=lmax (see below) with equal weights. If this argument is provided, the values of ells, bpws and weights are ignored.
  • lmax (int) – integer value corresponding to the maximum multipole used by these bandpowers. If None, it will be set to 3*nside-1. In any case the actual maximum multipole will be chosen as the minimum of lmax, 3*nside-1 and the maximum element of ells (e.g. if you are using CAR maps and don’t care about nside, you can pass whatever lmax you want and e.g. nside=lmax).
  • is_Dell (boolean) – if True, the output of all pseudo-Cl computations carried out using this bandpower scheme (e.g. from pymaster.workspaces.NmtWorkspace.decouple_cell()) will be multiplied by ell * (ell + 1) / 2 * PI, where ell is the multipole order (no prefactor otherwise).
  • f_ell (array-like) – if present, this is array represents an ell-dependent function that will be multiplied by all pseudo-Cl computations carried out using this bandpower scheme. If not None, the value of is_Dell is ignored.
bin_cell(cls_in)

Bins a power spectrum into bandpowers. This is carried out as a weighted average over the multipoles in each bandpower.

Parameters:cls_in (array-like) – 2D array of power spectra
Returns:array of bandpowers
classmethod from_edges(ell_ini, ell_end, is_Dell=False)

Convenience constructor for general equal-weight bands. All ells in the interval [ell_ini, ell_end) will be binned with equal weights across the band.

Parameters:
  • ell_ini (int) – array containing the lower edges of each bandpower.
  • ell_end (int) – array containing the upper edges of each bandpower.
  • is_Dell (boolean) – if True, the output of all pseudo-Cl computations carried out using this bandpower scheme (e.g. from pymaster.workspaces.NmtWorkspace.decouple_cell()) will be multiplied by ell * (ell + 1) / 2 * PI, where ell is the multipole order (no prefactor otherwise).
classmethod from_lmax_linear(lmax, nlb, is_Dell=False)

Convenience constructor for generic linear binning.

Parameters:
  • lmax (int) – integer value corresponding to the maximum multipole used by these bandpowers.
  • nlb (int) – integer value corresponding to a constant bandpower width. I.e. the bandpowers will be defined as consecutive sets of nlb multipoles from l=2 to l=lmax with equal weights.
  • is_Dell (boolean) – if True, the output of all pseudo-Cl computations carried out using this bandpower scheme (e.g. from pymaster.workspaces.NmtWorkspace.decouple_cell()) will be multiplied by ell * (ell + 1) / 2 * PI, where ell is the multipole order (no prefactor otherwise).
classmethod from_nside_linear(nside, nlb, is_Dell=False)

Convenience constructor for HEALPix maps with linear binning.

Parameters:
  • nside (int) – HEALPix nside resolution parameter of the maps you intend to correlate. The maximum multipole considered for bandpowers will be 3*nside-1.
  • nlb (int) – integer value corresponding to a constant bandpower width. I.e. the bandpowers will be defined as consecutive sets of nlb multipoles from l=2 to l=lmax with equal weights.
  • is_Dell (boolean) – if True, the output of all pseudo-Cl computations carried out using this bandpower scheme (e.g. from pymaster.workspaces.NmtWorkspace.decouple_cell()) will be multiplied by ell * (ell + 1) / 2 * PI, where ell is the multipole order (no prefactor otherwise).
get_effective_ells()

Returns an array with the effective multipole associated to each bandpower. These are computed as a weighted average of the multipoles within each bandpower.

Returns:effective multipoles for each bandpower
get_ell_list(ibin)

Returns an array with the multipoles in the ibin-th bandpower

Parameters:ibin (int) – bandpower index
Returns:multipoles associated with bandpower ibin
get_ell_max(b)

Returns the maximum ell value used by bandpower with index b

Returns:maximum ell value
get_ell_min(b)

Returns the minimum ell value used by bandpower with index b

Returns:minimum ell value
get_n_bands()

Returns the number of bandpowers stored in this object

Returns:number of bandpowers
get_nell_list()

Returns an array with the number of multipoles in each bandpower stored in this object

Returns:number of multipoles per bandpower
get_weight_list(ibin)

Returns an array with the weights associated to each multipole in the ibin-th bandpower

Parameters:ibin (int) – bandpower index
Returns:weights associated to multipoles in bandpower ibin
unbin_cell(cls_in)

Un-bins a set of bandpowers into a power spectrum. This is simply done by assigning a constant value for every multipole in each bandpower (corresponding to the value of that bandpower).

Parameters:cls_in (array-like) – array of bandpowers
Returns:array of power spectra
class pymaster.bins.NmtBinFlat(l0, lf)

An NmtBinFlat object defines the set of bandpowers used in the computation of the pseudo-Cl estimator. The definition of bandpowers is described in Section 3.6 of the scientific documentation. Note that currently pymaster only supports top-hat bandpowers for flat-sky power spectra.

Parameters:
  • l0 (array-like) – array of floats corresponding to the lower bound of each bandpower.
  • lf (array-like) – array of floats corresponding to the upper bound of each bandpower. lf should have the same shape as l0
bin_cell(ells, cls_in)

Bins a power spectrum into bandpowers. This is carried out as a weighted average over the multipoles in each bandpower.

Parameters:
  • ells (array-like) – multipole values at which the input power spectra are defined
  • cls_in (array-like) – 2D array of input power spectra
Returns:

array of bandpowers
get_effective_ells()

Returns an array with the effective multipole associated to each bandpower. These are computed as a weighted average of the multipoles within each bandpower.

Returns:effective multipoles for each bandpower
get_n_bands()

Returns the number of bandpowers stored in this object

Returns:number of bandpowers
unbin_cell(cls_in, ells)

Un-bins a set of bandpowers into a power spectrum. This is simply done by assigning a constant value for every multipole in each bandpower (corresponding to the value of that bandpower).

Parameters:
  • cls_in (array-like) – array of bandpowers
  • ells (array-like) – array of multipoles at which the power spectra should be intepolated
Returns:

array of power spectra
class pymaster.workspaces.NmtWorkspace

NmtWorkspace objects are used to compute and store the coupling matrix associated with an incomplete sky coverage, and used in the MASTER algorithm. When initialized, this object is practically empty. The information describing the coupling matrix must be computed or read from a file afterwards.

compute_coupling_matrix(fl1, fl2, bins, is_teb=False, n_iter=3, lmax_mask=-1, l_toeplitz=-1, l_exact=-1, dl_band=-1)

Computes coupling matrix associated with the cross-power spectrum of two NmtFields and an NmtBin binning scheme. Note that the mode coupling matrix will only contain ells up to the maximum multipole included in the NmtBin bandpowers.

Parameters:
  • fl1,fl2 (NmtField) – fields to correlate
  • bin (NmtBin) – binning scheme
  • is_teb (boolean) – if true, all mode-coupling matrices (0-0,0-2,2-2) will be computed at the same time. In this case, fl1 must be a spin-0 field and fl1 must be spin-2.
  • n_iter – number of iterations when computing a_lms.
  • lmax_mask – maximum multipole for masks. If smaller than the maximum multipoles of the fields, it will be set to that.
  • l_toeplitz – if a positive number, the Toeplitz approximation described in Louis et al. 2020 (arXiv:2010.14344) will be used. In that case, this quantity corresponds to ell_toeplitz in Fig. 3 of that paper.
  • l_exact – if l_toeplitz>0, this quantity corresponds to ell_exact in Fig. 3 of Louis et al. 2020. Ignored if l_toeplitz<=0.
  • dl_band – if l_toeplitz>0, this quantity corresponds to Delta ell_band in Fig. 3 of Louis et al. 2020. Ignored if l_toeplitz<=0.
couple_cell(cl_in)

Convolves a set of input power spectra with a coupling matrix (see Eq. 6 of the C API documentation).

Parameters:cl_in – set of input power spectra. The number of power spectra must correspond to the spins of the two fields that this NmtWorkspace object was initialized with (i.e. 1 for two spin-0 fields, 2 for one spin-0 and one spin-2 field and 4 for two spin-2 fields).
Returns:coupled power spectrum
decouple_cell(cl_in, cl_bias=None, cl_noise=None)

Decouples a set of pseudo-Cl power spectra into a set of bandpowers by inverting the binned coupling matrix (se Eq. 4 of the C API documentation).

Parameters:
  • cl_in – set of input power spectra. The number of power spectra must correspond to the spins of the two fields that this NmtWorkspace object was initialized with (i.e. 1 for two spin-0 fields, 2 for one spin-0 and one spin-2 field, 4 for two spin-2 fields and 7 if this NmtWorkspace was created using is_teb=True).
  • cl_bias – bias to the power spectrum associated to contaminant residuals (optional). This can be computed through pymaster.deprojection_bias().
  • cl_noise – noise bias (i.e. angular power spectrum of masked noise realizations).
Returns:

set of decoupled bandpowers
get_bandpower_windows()

Get bandpower window functions. Convolve the theory power spectra with these as an alternative to the combination decouple_cell(couple_cell(.

Returns:bandpower windows with shape [n_cls, n_bpws, n_cls, lmax+1].
get_coupling_matrix()

Returns the currently stored mode-coupling matrix.

Returns:mode-coupling matrix. The matrix will have shape [nrows,nrows], with nrows = n_cls * n_ells, where n_cls is the number of power spectra (1, 2 or 4 for spin0-0, spin0-2 and spin2-2 correlations) and n_ells = lmax + 1 (normally lmax = 3 * nside - 1). The assumed ordering of power spectra is such that the l-th element of the i-th power spectrum be stored with index l * n_cls + i.
read_from(fname)

Reads the contents of an NmtWorkspace object from a FITS file.

Parameters:fname (str) – input file name
update_coupling_matrix(new_matrix)

Updates the stored mode-coupling matrix.

The new matrix (new_matrix) must have shape [nrows,nrows], with nrows = n_cls * n_ells, where n_cls is the number of power spectra (1, 2 or 4 for spin0-0, spin0-2 and spin2-2 correlations) and n_ells = lmax + 1 (normally lmax = 3 * nside - 1). The assumed ordering of power spectra is such that the l-th element of the i-th power spectrum be stored with index l * n_cls + i.

Parameters:new_matrix – matrix that will replace the mode-coupling matrix.
write_to(fname)

Writes the contents of an NmtWorkspace object to a FITS file.

Parameters:fname (str) – output file name
class pymaster.workspaces.NmtWorkspaceFlat

NmtWorkspaceFlat objects are used to compute and store the coupling matrix associated with an incomplete sky coverage, and used in the flat-sky version of the MASTER algorithm. When initialized, this object is practically empty. The information describing the coupling matrix must be computed or read from a file afterwards.

compute_coupling_matrix(fl1, fl2, bins, ell_cut_x=[1.0, -1.0], ell_cut_y=[1.0, -1.0], is_teb=False)

Computes coupling matrix associated with the cross-power spectrum of two NmtFieldFlats and an NmtBinFlat binning scheme.

Parameters:
  • fl1,fl2 (NmtFieldFlat) – fields to correlate
  • bin (NmtBinFlat) – binning scheme
  • ell_cut_x (float(2)) – remove all modes with ell_x in the interval [ell_cut_x[0],ell_cut_x[1]] from the calculation.
  • ell_cut_y (float(2)) – remove all modes with ell_y in the interval [ell_cut_y[0],ell_cut_y[1]] from the calculation.
  • is_teb (boolean) – if true, all mode-coupling matrices (0-0,0-2,2-2) will be computed at the same time. In this case, fl1 must be a spin-0 field and fl1 must be spin-2.
couple_cell(ells, cl_in)

Convolves a set of input power spectra with a coupling matrix (see Eq. 6 of the C API documentation).

Parameters:
  • ells – list of multipoles on which the input power spectra are defined
  • cl_in – set of input power spectra. The number of power spectra must correspond to the spins of the two fields that this NmtWorkspaceFlat object was initialized with (i.e. 1 for two spin-0 fields, 2 for one spin-0 and one spin-2 field and 4 for two spin-2 fields).
Returns:

coupled power spectrum. The coupled power spectra are returned at the multipoles returned by calling get_ell_sampling() for any of the fields that were used to generate the workspace.
decouple_cell(cl_in, cl_bias=None, cl_noise=None)

Decouples a set of pseudo-Cl power spectra into a set of bandpowers by inverting the binned coupling matrix (se Eq. 4 of the C API documentation).

Parameters:
  • cl_in – set of input power spectra. The number of power spectra must correspond to the spins of the two fields that this NmtWorkspaceFlat object was initialized with (i.e. 1 for two spin-0 fields, 2 for one spin-0 and one spin-2 field, 4 for two spin-2 fields and 7 if this NmtWorkspaceFlat was created using is_teb=True). These power spectra must be defined at the multipoles returned by get_ell_sampling() for any of the fields used to create the workspace.
  • cl_bias – bias to the power spectrum associated to contaminant residuals (optional). This can be computed through pymaster.deprojection_bias_flat().
  • cl_noise – noise bias (i.e. angular power spectrum of masked noise realizations).
Returns:

set of decoupled bandpowers
read_from(fname)

Reads the contents of an NmtWorkspaceFlat object from a FITS file.

Parameters:fname (str) – input file name
write_to(fname)

Writes the contents of an NmtWorkspaceFlat object to a FITS file.

Parameters:fname (str) – output file name
pymaster.workspaces.compute_coupled_cell(f1, f2)

Computes the full-sky angular power spectra of two masked fields (f1 and f2) without aiming to deconvolve the mode-coupling matrix. Effectively, this is equivalent to calling the usual HEALPix anafast routine on the masked and contaminant-cleaned maps.

Parameters:f1,f2 (NmtField) – fields to correlate
Returns:array of coupled power spectra
pymaster.workspaces.compute_coupled_cell_flat(f1, f2, b, ell_cut_x=[1.0, -1.0], ell_cut_y=[1.0, -1.0])

Computes the angular power spectra of two masked flat-sky fields (f1 and f2) without aiming to deconvolve the mode-coupling matrix. Effectively, this is equivalent to computing the map FFTs and averaging over rings of wavenumber. The returned power spectrum is defined at the multipoles returned by the method get_ell_sampling() of either f1 or f2.

Parameters:
  • f1,f2 (NmtFieldFlat) – fields to correlate
  • b (NmtBinFlat) – binning scheme defining output bandpower
  • ell_cut_x (float(2)) – remove all modes with ell_x in the interval [ell_cut_x[0],ell_cut_x[1]] from the calculation.
  • ell_cut_y (float(2)) – remove all modes with ell_y in the interval [ell_cut_y[0],ell_cut_y[1]] from the calculation.
Returns:

array of coupled power spectra
pymaster.workspaces.compute_full_master(f1, f2, b, cl_noise=None, cl_guess=None, workspace=None, n_iter=3, lmax_mask=-1, l_toeplitz=-1, l_exact=-1, dl_band=-1)

Computes the full MASTER estimate of the power spectrum of two fields (f1 and f2). This is equivalent to successively calling:

  • pymaster.NmtWorkspace.compute_coupling_matrix()
  • pymaster.deprojection_bias()
  • pymaster.compute_coupled_cell()
  • pymaster.NmtWorkspace.decouple_cell()
Parameters:
  • f1,f2 (NmtField) – fields to correlate
  • b (NmtBin) – binning scheme defining output bandpower
  • cl_noise – noise bias (i.e. angular power spectrum of masked noise realizations) (optional).
  • cl_guess – set of power spectra corresponding to a best-guess of the true power spectra of f1 and f2. Needed only to compute the contaminant cleaning bias (optional).
  • workspace (NmtWorkspace) – object containing the mode-coupling matrix associated with an incomplete sky coverage. If provided, the function will skip the computation of the mode-coupling matrix and use the information encoded in this object.
  • n_iter – number of iterations when computing a_lms.
  • lmax_mask – maximum multipole for masks. If smaller than the maximum multipoles of the fields, it will be set to that.
  • l_toeplitz – if a positive number, the Toeplitz approximation described in Louis et al. 2020 (arXiv:2010.14344) will be used. In that case, this quantity corresponds to ell_toeplitz in Fig. 3 of that paper.
  • l_exact – if l_toeplitz>0, this quantity corresponds to ell_exact in Fig. 3 of Louis et al. 2020. Ignored if l_toeplitz<=0.
  • dl_band – if l_toeplitz>0, this quantity corresponds to Delta ell_band in Fig. 3 of Louis et al. 2020. Ignored if l_toeplitz<=0.
Returns:

set of decoupled bandpowers
pymaster.workspaces.compute_full_master_flat(f1, f2, b, cl_noise=None, cl_guess=None, ells_guess=None, workspace=None, ell_cut_x=[1.0, -1.0], ell_cut_y=[1.0, -1.0])

Computes the full MASTER estimate of the power spectrum of two flat-sky fields (f1 and f2). This is equivalent to successively calling:

  • pymaster.NmtWorkspaceFlat.compute_coupling_matrix()
  • pymaster.deprojection_bias_flat()
  • pymaster.compute_coupled_cell_flat()
  • pymaster.NmtWorkspaceFlat.decouple_cell()
Parameters:
  • f1,f2 (NmtFieldFlat) – fields to correlate
  • b (NmtBinFlat) – binning scheme defining output bandpower
  • cl_noise – noise bias (i.e. angular power spectrum of masked noise realizations) (optional). This power spectrum should correspond to the bandpowers defined by b.
  • cl_guess – set of power spectra corresponding to a best-guess of the true power spectra of f1 and f2. Needed only to compute the contaminant cleaning bias (optional).
  • ells_guess – multipoles at which cl_guess is defined.
  • workspace (NmtWorkspaceFlat) – object containing the mode-coupling matrix associated with an incomplete sky coverage. If provided, the function will skip the computation of the mode-coupling matrix and use the information encoded in this object.
  • nell_rebin (int) – number of sub-intervals into which the base k-intervals will be sub-sampled to compute the coupling matrix
  • ell_cut_x (float(2)) – remove all modes with ell_x in the interval [ell_cut_x[0],ell_cut_x[1]] from the calculation.
  • ell_cut_y (float(2)) – remove all modes with ell_y in the interval [ell_cut_y[0],ell_cut_y[1]] from the calculation.
Returns:

set of decoupled bandpowers
pymaster.workspaces.deprojection_bias(f1, f2, cls_guess, n_iter=3)

Computes the bias associated to contaminant removal to the cross-pseudo-Cl of two fields.

Parameters:
  • f1,f2 (NmtField) – fields to correlate
  • cls_guess – set of power spectra corresponding to a best-guess of the true power spectra of f1 and f2.
  • n_iter – number of iterations when computing a_lms.
Returns:

deprojection bias power spectra.
pymaster.workspaces.deprojection_bias_flat(f1, f2, b, ells, cls_guess, ell_cut_x=[1.0, -1.0], ell_cut_y=[1.0, -1.0])

Computes the bias associated to contaminant removal to the cross-pseudo-Cl of two flat-sky fields. The returned power spectrum is defined at the multipoles returned by the method get_ell_sampling() of either f1 or f2.

Parameters:
  • f1,f2 (NmtFieldFlat) – fields to correlate
  • b (NmtBinFlat) – binning scheme defining output bandpower
  • ells – list of multipoles on which the proposal power spectra are defined
  • cls_guess – set of power spectra corresponding to a best-guess of the true power spectra of f1 and f2.
  • ell_cut_x (float(2)) – remove all modes with ell_x in the interval [ell_cut_x[0],ell_cut_x[1]] from the calculation.
  • ell_cut_y (float(2)) – remove all modes with ell_y in the interval [ell_cut_y[0],ell_cut_y[1]] from the calculation.
Returns:

deprojection bias power spectra.
pymaster.workspaces.uncorr_noise_deprojection_bias(f1, map_var, n_iter=3)

Computes the bias associated to contaminant removal in the presence of uncorrelated inhomogeneous noise to the auto-pseudo-Cl of a given field f1.

Parameters:
  • f1 (NmtField) – fields to correlate
  • map_cls_guess – array containing a HEALPix map corresponding to the local noise variance (in one sterad).
  • n_iter – number of iterations when computing a_lms.
Returns:

deprojection bias power spectra.
class pymaster.covariance.NmtCovarianceWorkspace

NmtCovarianceWorkspace objects are used to compute and store the coupling coefficients needed to calculate the Gaussian covariance matrix under the approximations in astro-ph/0307515 and arXiv/1609.09730. When initialized, this object is practically empty. The information describing the coupling coefficients must be computed or read from a file afterwards.

compute_coupling_coefficients(fla1, fla2, flb1=None, flb2=None, lmax=None, n_iter=3, l_toeplitz=-1, l_exact=-1, dl_band=-1)

Computes coupling coefficients of the Gaussian covariance between the power spectra of two pairs of NmtField objects (fla1, fla2, flb1 and flb2). Note that you can reuse this workspace for the covariance of power spectra between any pairs of fields as long as the fields have the same masks as those passed to this function, and as long as the binning scheme used are also the same.

Parameters:
  • fla1,fla2 (NmtField) – fields contributing to the first power spectrum whose covariance you want to compute.
  • flb1,flb2 (NmtField) – fields contributing to the second power spectrum whose covariance you want to compute. If None, fla1,fla2 will be used.
  • n_iter – number of iterations when computing a_lms.
  • l_toeplitz – if a positive number, the Toeplitz approximation described in Louis et al. 2020 (arXiv:2010.14344) will be used. In that case, this quantity corresponds to ell_toeplitz in Fig. 3 of that paper.
  • l_exact – if l_toeplitz>0, this quantity corresponds to ell_exact in Fig. 3 of Louis et al. 2020. Ignored if l_toeplitz<=0.
  • dl_band – if l_toeplitz>0, this quantity corresponds to Delta ell_band in Fig. 3 of Louis et al. 2020. Ignored if l_toeplitz<=0.
read_from(fname)

Reads the contents of an NmtCovarianceWorkspace object from a FITS file.

Parameters:fname (str) – input file name
write_to(fname)

Writes the contents of an NmtCovarianceWorkspace object to a FITS file.

Parameters:fname (str) – output file name
class pymaster.covariance.NmtCovarianceWorkspaceFlat

NmtCovarianceWorkspaceFlat objects are used to compute and store the coupling coefficients needed to calculate the Gaussian covariance matrix under a flat-sky version the Efstathiou approximations in astro-ph/0307515 and arXiv/1609.09730. When initialized, this object is practically empty. The information describing the coupling coefficients must be computed or read from a file afterwards.

compute_coupling_coefficients(fla1, fla2, bin_a, flb1=None, flb2=None, bin_b=None)

Computes coupling coefficients of the Gaussian covariance between the power spectra of two pairs of NmtFieldFlat objects (fla1, fla2, flb1 and flb2). Note that you can reuse this workspace for the covariance of power spectra between any pairs of fields as long as the fields have the same masks as those passed to this function, and as long as the binning scheme used are also the same.

Parameters:
  • fla1,fla2 (NmtFieldFlat) – fields contributing to the first power spectrum whose covariance you want to compute.
  • bin_a (NmtBinFlat) – binning scheme for the first power spectrum.
  • flb1,flb2 (NmtFieldFlat) – fields contributing to the second power spectrum whose covariance you want to compute. If None, fla1,fla2 will be used.
  • bin_b (NmtBinFlat) – binning scheme for the second power spectrum. If none, bin_a will be used.
read_from(fname)

Reads the contents of an NmtCovarianceWorkspaceFlat object from a FITS file.

Parameters:fname (str) – input file name
write_to(fname)

Writes the contents of an NmtCovarianceWorkspaceFlat object to a FITS file.

Parameters:fname (str) – output file name
pymaster.covariance.gaussian_covariance(cw, spin_a1, spin_a2, spin_b1, spin_b2, cla1b1, cla1b2, cla2b1, cla2b2, wa, wb=None, coupled=False)

Computes Gaussian covariance matrix for power spectra using the information precomputed in cw (a NmtCovarianceWorkspace object). cw should have been initialized using four NmtField objects (let’s call them a1, a2, b1 and b2), corresponding to the two pairs of fields whose power spectra we want the covariance of. These power spectra should have been computed using two NmtWorkspace objects, wa and wb, which must be passed as arguments of this function (the power spectrum for fields a1 and a2 was computed with wa, and that of b1 and b2 with wb). Using the same notation, clXnYm should be a prediction for the power spectrum between fields Xn and Ym. These predicted input power spectra should be defined for all ells <=3*nside (where nside is the HEALPix resolution parameter of the fields that were correlated).

Parameters:
  • cw (NmtCovarianceWorkspace) – workspaces containing the precomputed coupling coefficients.
  • spin_Xn (int) – spin for the n-th field in pair X.
  • clXnYm – prediction for the cross-power spectrum between fields Xn and Ym.
  • wX (NmtWorkspace) – workspace containing the mode-coupling matrix pair X. If wb is None, the code will assume wb=wa.
  • coupled – if True, the covariance matrix of the mode-coupled pseudo-Cls will be computed. Otherwise it’ll be the covariance of decoupled bandpowers.
pymaster.covariance.gaussian_covariance_flat(cw, spin_a1, spin_a2, spin_b1, spin_b2, larr, cla1b1, cla1b2, cla2b1, cla2b2, wa, wb=None)

Computes Gaussian covariance matrix for flat-sky power spectra using the information precomputed in cw (a NmtCovarianceWorkspaceFlat object). cw should have been initialized using four NmtFieldFlat objects (let’s call them a1, a2, b1 and b2), corresponding to the two pairs of fields whose power spectra we want the covariance of. These power spectra should have been computed using two NmtWorkspaceFlat objects, wa and wb, which must be passed as arguments of this function (the power spectrum for fields a1 and a2 was computed with wa, and that of b1 and b2 with wb). Using the same notation, clXnYm should be a prediction for the power spectrum between fields Xn and Ym. These predicted input power spectra should be defined in a sufficiently well sampled range of ells given the map properties from which the power spectra were computed. The values of ell at which they are sampled are given by larr.

Please note that, while the method used to estimate these covariance is sufficiently accurate in a large number of scenarios, it is based on a numerical approximation, and its accuracy should be assessed if in doubt. In particular, we discourage users from using it to compute any covariance matrix involving B-mode components.

Parameters:
  • cw (NmtCovarianceWorkspaceFlat) – workspaces containing the precomputed coupling coefficients.
  • spin_Xn (int) – spin for the n-th field in pair X.
  • larr – values of ell at which the following power spectra are computed.
  • clXnYm – prediction for the cross-power spectrum between fields Xn and Ym.
  • wX (NmtWorkspaceFlat) – workspace containing the mode-coupling matrix pair X. If wb is None, the code will assume wb=wa.
class pymaster.utils.NmtWCSTranslator(wcs, axes)

This class takes care of interpreting a WCS object in terms of a Clenshaw-Curtis grid.

Parameters:
pymaster.utils.mask_apodization(mask_in, aposize, apotype='C1')

Apodizes a mask with an given apodization scale using different methods. A given pixel is determined to be “masked” if its value is 0.

Parameters:
  • mask_in – input mask, provided as an array of floats corresponding to a HEALPix map in RING order.
  • aposize – apodization scale in degrees.
  • apotype – apodization type. Three methods implemented: “C1”, “C2” and “Smooth”. See the description of the C-function nmt_apodize_mask in the C API documentation for a full description of these methods.
Returns:

apodized mask as a HEALPix map
pymaster.utils.mask_apodization_flat(mask_in, lx, ly, aposize, apotype='C1')

Apodizes a flat-sky mask with an given apodization scale using different methods. A given pixel is determined to be “masked” if its value is 0.

Parameters:
  • mask_in – input mask, provided as a 2D array (ny,nx) of floats.
  • lx (float) – patch size in the x-axis (in radians)
  • ly (float) – patch size in the y-axis (in radians)
  • aposize – apodization scale in degrees.
  • apotype – apodization type. Three methods implemented: “C1”, “C2” and “Smooth”. See the description of the C-function nmt_apodize_mask in the C API documentation for a full description of these methods.
Returns:

apodized mask as a 2D array (ny,nx)
pymaster.utils.synfast_flat(nx, ny, lx, ly, cls, spin_arr, beam=None, seed=-1)

Generates a flat-sky Gaussian random field according to a given power spectrum. This function is the flat-sky equivalent of healpy’s synfast.

Parameters:
  • nx (int) – number of pixels in the x-axis
  • ny (int) – number of pixels in the y-axis
  • lx (float) – patch size in the x-axis (in radians)
  • ly (float) – patch size in the y-axis (in radians)
  • cls (array-like) – array containing power spectra. Shape should be [n_cls][n_ell], where n_cls is the number of power spectra needed to define all the fields. This should be n_cls = n_maps * (n_maps + 1) / 2, where n_maps is the total number of maps required (1 for each spin-0 field, 2 for each spin>0 field). Power spectra must be provided only for the upper-triangular part in row-major order (e.g. if n_maps is 3, there will be 6 power spectra ordered as [1-1,1-2,1-3,2-2,2-3,3-3].
  • spin_arr (array-like) – array containing the spins of all the fields to generate.
  • array-like (beam) – 2D array containing the instrumental beam of each field to simulate (the output map(s) will be convolved with it)
  • seed (int) – RNG seed. If negative, will use a random seed.
Returns:

a number of arrays (1 for each spin-0 field, 2 for each spin-2 field) of size (ny,nx) containing the simulated maps.
pymaster.utils.synfast_spherical(nside, cls, spin_arr, beam=None, seed=-1, wcs=None)

Generates a full-sky Gaussian random field according to a given power spectrum. This function should produce outputs similar to healpy’s synfast.

Parameters:
  • nside (int) – HEALpix resolution parameter. If you want rectangular pixel maps, ignore this parameter and pass a WCS object as wcs (see below).
  • cls (array-like) – array containing power spectra. Shape should be [n_cls][n_ell], where n_cls is the number of power spectra needed to define all the fields. This should be n_cls = n_maps * (n_maps + 1) / 2, where n_maps is the total number of maps required (1 for each spin-0 field, 2 for each spin>0 field). Power spectra must be provided only for the upper-triangular part in row-major order (e.g. if n_maps is 3, there will be 6 power spectra ordered as [1-1,1-2,1-3,2-2,2-3,3-3].
  • spin_arr (array-like) – array containing the spins of all the fields to generate.
  • array-like (beam) – 2D array containing the instrumental beam of each field to simulate (the output map(s) will be convolved with it)
  • seed (int) – RNG seed. If negative, will use a random seed.
  • wcs – a WCS object (see http://docs.astropy.org/en/stable/wcs/index.html).
Returns:

a number of full-sky maps (1 for each spin-0 field, 2 for each spin-2 field).