LensPix: Fast MPI full sky transforms for HEALPix

NOTE:LensPix is an old code, the modern LensPyx python code should be much faster and more accurate.

NEW: (Apr 14) Healpix/MPI compatibility updates; new sample LensReconExample.f90 code; AsymmCouplings changes for use in E/B separation
(Sept 12) functions for read/write of lensing potential and deflection (spin) maps, a_lm rotation, etc.
(Dec 11) functions for transforming arrays of maps and alms, packed alm and array formats made public, HealpixMap_MakeSeparableNGfromG function for simulating local non-Gaussianity.

The code requires HEALPix, with features:

• Simulate lensed temperature and polarized CMB maps with efficient MPI parallelization
• Cubic interpolation for fast and accurate lensed simulation
• Non-lensing routines for spherical harmonic transforms on a cluster using MPI
• Compute gradient maps (spin 1, and also up to spin 3)
• Compute harmonic coefficients of maps of spin fields
• Various basic visualisation routines for producing ppm and matlab files without IDL (as above)
• Convenient F90 wrapper code for some HEALPix structures and routines (HealpixObj)
• Calculating mask Alm coupling matrices, cut-sky supported modes and harmonic E/B separation
• Transforms of arrays of maps and calculating cross-spectra
• Includes latest updates WeightMixer cross power spectrum code and functions

Download the code

A version-tracked (git) copy of the code is available on a shared repository site; after downloading please email if you would instead like to have access to the repository and hence latest builds.

There is no HEALPix-quality documentation, though details of the theory are given in the appendix of my original paper and the splining scheme in Appendix E.4 of arXiv:0801.0554. A discussion of the lensing potential-E correlation is in Lewis, Challinor & Hanson (arXiv:1101.2234). The download includes a sample program SimLens.f90 which should give you a good idea of how things work using the HealpixObj wrapper routines. If you prefer lower-level HEALPix-style routines you can also use spin_alm_tools.f90 directly. The code can be run on single CPUs without MPI, though for high resolution fast lensed map making you will likely need more than 2GB of memory. Rather than wasting lots of time documenting things now (which probably no one will read), please ask any questions on CosmoCoffee and I will respond there so everyone can read it.

Note for intel compilers: you may need to compile with -heap-arrays to avoid getting segmentation faults using some routines with large l_max.

To run the demo SimLens program you will need to edit the Makefile and job file, and also change params.ini to have the correct path to your HEALPix data files. Then run with "qsub job" or "./simlens params.ini". The program will read in some unlensed Cls from a text file, simulate unlensed a_lms, compute the gradient of the potential, then use this to construct a lensed map, outputting the text lensed C_ls in that realization.

The way MPI is implemented is that you call a HealpixInit(H,..) function with an nside and lmax you want to use before anything else. All but the main thread stay in this routine until HealpixFree(H) is called by the main thread. The main thread calls transforms like HealpixMap2Alm(H,....), which are automatically farmed out over the other threads. This is demonstrated in the SimLens demo program. There are numerous ways things could be improved.

There are almost certainly bugs, please let me know if you find out what they are, as well as any improvements. The implementation is not complete, in that there are no fits routines for storing spin fields or lensing maps, etc, etc. Using interpolation for the pixel re-mapping might well save a lot of time and memory, and probably worth implementing if you want to do lots of simulations.

For an accurate calculation of the lensing potential, unlensed and lensed CMB power spectra see CAMB. For parameter estimation including CMB lensing see CosmoMC.

Accuracy

The cubic interpolation method works by simulating an equicylindricallly-pixellized high-resolution map, then doing cubic interpolation (ref) on this grid for the deflection angles in each Healpix pixel. The high-res pixels are equally spaced in θ, but the number of pixels in φ depends on the number of MPI threads (for each MPI slice it chooses the number of pixels in φ so that the maximum angular spacing is similar to or smaller than the θ-spacing, and the number of φ-pixels is of the form 2ⁿ3^m so that the FFT is fast). To this extent the result is not the invariant under changing the number of threads. Using a simpler cubic interpolation scheme would speed up the code significantly, however it appears to make it much less stable - the result does then not converge as you increase interp_factor.

For a given l_max using nside=2048 and interp_factor=1.5 seems to be give power spectra accurate to about 0.1% at l<2000 (or nside=1024 and interp_factor=3; for accurate lensing basically needs an effective high-res nside of about 3000). Errors on BB power spectrum are larger. Note that all spectra are affected by choice of l_max; BB especially requires high l_max for accurate results (and an input power spectrum with non-linear corrections).

Comparison with CAMB

The expectation over simulations of the simulated lensed C_l should agree with the lensed power spectrum from CAMB (using accurate full-sky method of astro-ph/0502425). For this to be true either l_max needs to be much larger than the range of interest, so all relevant power is included, or l_max must be the same in both cases (and of course pixellization/interpolation errors must be negligible). For l_max=2250 the smoothed average fractional difference from 50 simulated spectra does indeed agree well (~0.1%) on all spectra (with usual caveats about BB).