pyilc

pyilc is a pure-python implementation of the needlet internal linear combination (NILC) algorithm for CMB component separation. Harmonic-space ILC is also implemented in the code. For details, see McCarthy & Hill (2023) arXiv:2307.01043.

diffusive_inpaint

This repository also includes an inpainting code, diffusive_inpaint, that diffusively inpaints a masked region with the mean of the unmasked neighboring pixels. This README is intended for pyilc, not diffusive_inpaint; in the diffusive_inpaint sub-directory, we include a sample .py file diffusive_inpaint/diffusive_inpaint_example.py, which should make clear how to use diffusive_inpaint/diffusive_inpaint.py.

Requirements

pyilc requires python3, numpy, matplotlib, and healpy (and all of their requirements). Recently, support for maps in CAR pixelization has been added; to use this functionality (which is still in 'beta'), you will need pixell.

Using the code

pyilc is public; if you use it in a publication, please cite the paper arXiv:2307.01043 and (optionally) arXiv:2308.16260.

Additionally, if you use NILC you should cite the original NILC reference, https://ui.adsabs.harvard.edu/abs/arXiv:0807.0773.

If you use an ILC that deprojects some component, please cite the constrained ILC references https://ui.adsabs.harvard.edu/abs/2009ApJ...694..222C/abstract and https://ui.adsabs.harvard.edu/abs/arXiv:1006.5599.

If using a moment-based deprojection, please cite https://ui.adsabs.harvard.edu/abs/arXiv:1701.00274.

See here for slides with an overview of how to use pyilc, presented to the SO FG AWG (Simons Observatory Foreground Analysis Working Group).

Basic usage

To run pyilc, a yaml input file is required. pyilc is run by running the main.py file using python with this input file as an argument:

python pyilc/main.py sample_input.yml

We have included a sample input file pyilc_input_example_general.yml, which serves as documentation of the different input options. The main ones are described here.

We go into detail about the input and output structure below. In general, an input file will contain a list of input frequency maps in units of K_CMB on which the NILC is to be performed, along with a path specifying what directory the output should be saved in and a prefix and suffix with which to save the products. The output products that are saved are the needlet coefficients of the input maps, the (maps of) covariance and inverse frequency-frequency covariance matrices, and the ILC maps (and ILC weights if requested in the input file). Before performing NILC, the code will check whether these products already exist in the specified output directory with the specified prefix (in the case of the input needlet coefficients and the covariance and inverse covariance matrices) and the specified prefix AND suffix (in the case of the weights and final ILC map):

• If the map or weights exist, the code will not compute anything as the products already exist.
• If the covariance/inverse covariance matrices exist, the code will load these and use these to compute the final ILC weights and map, then save the weights/map with the specified prefix AND suffix.
• If no covariance products exist, the code will compute these, save them with the specified prefix, then use them to compute the weights and map and save them with the specified prefix AND suffix.

As the covariance matrix computation/inversion is often the most computationally expensive step in performing NILC, this allows several versions of ILC maps (different deprojections, SEDs, etc.) to be computed with the same input without recomputing the covariances.

Examples

There are some examples of usage in the notebooks/ folder, where we explicitly run a .yaml file from the input/ folder which will reproduce the ILC calculation undeprojected y map in arXiv:2307.01043 . Note that it first downloads the preprocessed input single-frequency maps from https://users.flatironinstitute.org/~fmccarthy/ymaps_PR4_McCH23/inpainted_input_maps/ (the preprocessing is done with diffusive_inpaint.py and is described in arXiv:2307.01043 ). There is also an example notebook which performs HILC on the Planck data to create a temperature map. The HILC calculation is much faster than the NILC calculation, and this can be run very quickly.

Input structure

Specifying the input maps

The maps on which the (N)ILC will be performed should all be saved separately at some location /path/to/location/map_freq_X.fits. These files should be included in the .yml file as a list of strings:

freq_map_files: ['/path/to/location/map_freq_X.fits',
                 '/path/to/location/map_freq_Y.fits',...]

Note that there is another input parameter N_freqs that is required in the input file, and which must be equal to the length of the freq_map_files list, or else an error will be thrown. The maps must be listed in order of decreasing resolution, as specified by the user-specified input beams, as described below.

Beams

There is also additional metadata about the input maps that must be included in the input file. In particular the beams with which the maps are convolved must be specified. There are two options: Gaussian beams, or more general, one-dimensional ($\ell$-dependent) beams. Gaussian beams are specified as follows:

beam_type: 'Gaussians'  
beam_FWHM_arcmin: [FWHM_X,FWHM_Y,...]

where FWHM_X is the FWHM of the beam, in arcminutes, of the map at /path/to/location/map_freq_X.fits. The FWHMS must be listed in the same order as the maps they correspond to in freq_map_files. Also, the FWHMS must be decreasing, i.e., the maps must be read in in order from lowest-resolution to highest-resolution, or else an error will be thrown (note this might mean that they are not read in in a monotonically increasing frequency order, e.g., if the maps come from different instruments).

Alternatively, 1D beams can be specified as follows:

beam_type: '1DBeams'
beam_files:['/path/to/location/beam_freq_X.txt',
            '/path/to/location/beam_freq_Y.txt',...]

where '/path/to/location/beam_freq_X.txt' contains an array of shape (LMAX,2) where the first column specifies the $\ell$ and the second column specifies the beam at $\ell$. LMAX should be at least as high as the ELLMAX (maximum multipole) at which the NILC is being performed (this is a user-specified parameter in the input file -- see below).

The maps should all be in units of CMB thermodynamic temperature, $\mathrm{K}_{\mathrm{CMB}}$.

Frequency coverage

In order to calculate the SED of a given component, the code needs to know what frequencies the maps correspond to. There are two options: delta-function passbands and performing passband integration. The former can be indicated by including the following in the input file:

bandpass_type: 'DeltaBandpasses'
freqs_delta_ghz: [freq_X, freq_Y, ...]

where freq_X refers to the frequency (in GHz) of the appropriate input map (the order should be the same as that of freq_map_files). To perform passband integration, the following should be indicated in the input file:

bandpass_type: 'ActualBandpasses'
freq_bp_files: ['/path/to/location/bandpass_freq_X.txt',
                '/path/to/location/bandpass_freq_Y.txt',,...]
               

where '/path/to/location/bandpass_freq_X.txt' is a file containing the passband (first column = frequency in GHz, second column = transmission in the same convention as used in Planck) of the appropriate input map (again, the order should be the same as that of freq_map_files).

Specifying the resolution of the output ILC map

The maximum multipole (ELLMAX) used in all spherical harmonic transformations should be specified with

ELLMAX: ELLMAX

The N_side of the output map should be specified with

N_side: N_side

This should be the same N_side as the highest-resolution input map.

Optionally, the FWHM (in arcminutes) at which NILC should be performed (and with which the output map is convolved) can be specified according to

perform_ILC_at_beam: FWHM

If this is not specified, the ILC will be performed on maps which are all convolved to the same beam as that of the highest-resolution input map.

Specifying the type of ILC to perform

Currently, pyilc can perform two types of ILC: NILC with Gaussian needlets and harmonic ILC.

Gaussian NILC

For Gaussian NILC with N_scales needlet scales, constructed from Gaussians with FWHMs of FWHM_1,..., the following needs to be included in the input file:

wavelet_type: 'GaussianNeedlets' 
N_scales: N_scales
GH_FWHM_arcmin: [FWHM_1,FWHM_2,...,]

GH_FWHM_arcmin should be a list of the FWHMs in arcminutes of the Gaussians which are used to define the needlet scales. There must be N_scales-1 entries. They must be in decreasing order.

Real-space filter size

The size of the real-space domains used for calculating the covariances is defined adaptively in the code by specifying a threshold for the ILC bias and calculating the area required such that the number of modes within the area will be enough to unsure the fractional ILC bias is lower than this threshold. This threshold can be modified by specifying

ILC_bias_tol: b_tol

If unspecified, this is set by default to 0.01.

HILC

For harmonic ILC (HILC), the following needs to be included in the input file

wavelet_type: 'TopHatHarmonic'
BinSize: 20

BinSize is **one integer