Power Spectrum Estimation

The galsim.pse module was developed largely by Joe Zuntz and tweaked by assorted GalSim developers. This development and testing of this module took place in a separate (private) repository before the code was moved into the GalSim repository, but there are some demonstrations of the performance of this code in devel/modules/lensing_engine.pdf.

class galsim.pse.PowerSpectrumEstimator(N=100, sky_size_deg=10.0, nbin=15)[source]

Class for estimating the shear power spectrum from gridded shears.

This class stores all the data used in power spectrum estimation that is fixed with the geometry of the problem - the binning and spin weighting factors.

The only public method is estimate(), which is called with 2D g1 and g2 arrays on a square grid. It assumes the flat sky approximation (where ell and k are interchangeable), and rebins the observed ell modes into a user-defined number of logarithimic bins in ell. Given that the grid parameters are precomputed and stored when the PowerSpectrumEstimator is initialized, computation of the PS for multiple sets of shears corresponding to the same grid setup can proceed more rapidly than if everything had to be recomputed each time.

Below is an example of how to use this code (relying on GalSim to provide the arrays of g1 and g2, though that is by no means required, and assuming that the user is sitting in the examples/ directory):

>>> grid_size = 10.  # Define the total grid extent, in degrees
>>> ngrid = 100      # Define the number of grid points in each dimension: (ngrid x ngrid)
>>> n_ell = 15       # Choose the number of logarithmic bins in ell or k for outputs
>>>
>>> # Define a lookup-table for the power spectrum as a function of k based on the outputs
>>> # of iCosmo (see demo11.py for more description of how this was generated).
>>> my_tab = galsim.LookupTable(file='data/cosmo-fid.zmed1.00.out')
>>>
>>> # Generate a galsim.PowerSpectrum with this P(k), noting the units.
>>> my_ps = galsim.PowerSpectrum(my_tab, units=galsim.radians)
>>>
>>> # Build a grid of shear values with the desired parameters.
>>> g1, g2 = my_ps.buildGrid(grid_spacing=grid_size/ngrid, ngrid=ngrid,
...                          units=galsim.degrees)
>>>
>>> # Initialize a PowerSpectrumEstimator with the chosen grid geometry and number of ell
>>> # bins. Note that these values are actually the default, so we didn't technically have
>>> # to specifythem.
>>> my_pse = galsim.pse.PowerSpectrumEstimator(ngrid, grid_size, n_ell)
>>>
>>> # Estimate the power based on this set of g1, g2.  If we get another set of shears for
>>> # the same grid geometry, we can reuse the same PowerSpectrumEstimator object.
>>> ell, P_e, P_b, P_eb = my_pse.estimate(g1, g2)

The output NumPy arrays ell, P_e, P_b, and P_eb contain the effective ell value, the E-mode auto-power spectrum, the B-mode auto-power spectrum, and the EB cross-power spectrum. The units are inverse radians for ell, and radians^2 for the output power spectra.

Some important notes:

  1. Power spectrum estimation requires a weight function which decides how the averaging is done across ell within each bin. By default that weighting is flat in ell using an analytic calculation of the area in ell space, but this is easy to change with the _bin_power function. (Note this area averaged bin weighting is only approximate for the higher frequency bins in which the lower ell edge is greater than pi * ngrid / grid_size, due to the annular ell region being cut off by the square grid edges beyond this value.) A keyword allows for weighting by the power itself.

  2. This is the power spectrum of the gridded data, not the underlying field - we do not account for the effects of the finite grid (basically, ignoring all the reasons why power spectrum estimation is hard - see devel/modules/lensing_engine.pdf in the GalSim repository). Users must account for the contribution of noise in g1, g2 and any masking.

  3. The binning is currently fixed as uniform in log(ell).

  4. The code for this class uses the notation of the GREAT10 handbook (Kitching et al. 2011, http://dx.doi.org/10.1214/11-AOAS484), equations 17-21.

galsim.pse.PowerSpectrumEstimator.__init__

estimate(g1, g2, weight_EE=False, weight_BB=False, theory_func=None)[source]

Compute the EE, BB, and EB power spectra of two 2D arrays g1 and g2.

For example usage, see the docstring for the PowerSpectrumEstimator class.

Parameters:

  • g1 – The shear component g1 as a square 2D NumPy array.

  • g2 – The shear component g2 as a square 2D NumPy array.

  • weight_EE – If True, then the E auto-power spectrum is re-computed weighting by the power within each logarithmically-spaced ell bin. [default: False]

  • weight_BB – If True, then the B auto-power spectrum is re-computed weighting by the power within each logarithmically-spaced ell bin. [default: False]

  • theory_func – An optional callable function that can be used to get an idealized value of power at each point on the grid, and then see what results it gives for our chosen ell binning. [default: None]