Two-point statistics from Heracles

This notebook demonstrates how Heracles extracts the two-point statistics from a 3×2pt catalogue.

Important note

This notebook is only meant to give you an idea of how Heracles works.

It does not show everything that Heracles can do.

This is a toy, treat it is such!

Setup

Some required imports, nothing fancy.

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
import matplotlib as mpl

Now the Heracles imports:

• The top-level heracles module contains all general user-facing functionality.

• The heracles.healpy module contains mappers based on the healpy package.

• The heracles.notebook module contains a progress bar based on the ipywidgets package.

import heracles
import heracles.healpy
from heracles.notebook import Progress

If there is an import error on the last line of the previous block, it means you need to install the ipywidgets package.

Data set

The example uses a prepared example data set. The catalogue contains about 75 million galaxies over a contiguous 5% of the sky.

To download the example data, run the following cell.

import helpers

with Progress("example data") as progress:
    helpers.get_example_data(progress)

Basic parameters

This is the resolution parameter for measuring spectra from HEALPix maps. Here, we use nside = 1024 since that is the resolution at which the example data has been created. A value lmax of approximate 1.5x nside is fairly safe in terms of errors introduced by HEALPix.

nside = 1024
lmax = 1500

Catalogues

Heracles provides a flexible interface for loading catalogues from FITS files or arrays. It also provides a base class that can quickly be extended, e.g., to databases or more.

Here we use the FITS interface to read a catalogue from file. We could specify the columns to read, but the example catalogue does not have much more than we need here.

Importantly, catalogues are never read into memory all at once. Their page_size attribute determines how many rows are read at a time.

# load the FITS catalogue
catalog = heracles.FitsCatalog("catalog.fits")

Using the catalog.add_filter() method, we could add filters to the catalogue here, e.g., to strip rows with invalid values, or apply an extra mask.

Playground

The Catalog interface essentially provides an iterator over pages of rows. You can use it as such:

nrows = 0
for page in catalog:
    nrows += page.size

# no need to iterate to get the number of rows, really
assert nrows == catalog.size

print(f"there are {nrows:_} rows in your catalogue")

there are 76_889_615 rows in your catalogue

The page object is a mapping of column names to rows, with some additional features:

• The number of rows, page.size

• The names of columns, page.names

• A view of the underlying data, page.data

• Make a copy of the page, page.copy()

• Delete specific rows, page.delete(where)

• Return multiple columns at once, using page['a', 'b'] or page[['a', 'b']]

• Return columns while checking for invalid values, using page.get('a', 'b')

Visibility maps

The statistics of galaxy positions require knowledge of the a priori probability of detecting a galaxy at each point in the sky. We call this the visibility, and we use it in the form of a visibility maps, which is a full-sky maps of detection probabilities.

vmap = hp.read_map("vmap.fits.gz")

# fix UNSEEN pixels to zero and rescale to nside
vmap[vmap == hp.UNSEEN] = 0.0
vmap = hp.ud_grade(vmap, 2 * nside)

For a real survey, the visibility is a complicated function of position, observing conditions, selection, and tomographic bin. However, in this simulated example, the selection is the same for all positions and tomographic bins, and the visibility map is simply the footprint map of the survey.

# set visibility map of entire catalogue
catalog.visibility = vmap

Playground

Quick inspection of the visibility map:

hp.mollview(vmap, title="visibility", cmap="binary", bgcolor="none")
hp.graticule()
plt.show()

Tomographic binning

The catalog object can be used to read the entire catalogue. Of course, we would like to split our galaxies into individual tomographic bins. In the example data set, each galaxy has already been assigned a label for its tomographic bin, in the "BIN" column.

To perform the tomographic binning, we construct a dictionary that assigns an integer bin ID to a selection from the catalogue for each tomographic bin. This is done using the usual [] indexing syntax, which returns a new view of the FitsCatalog with the given selection applied.

catalogs = {i: catalog[f"BIN == {i}"] for i in range(1, 7)}

catalogs

{1: catalog.fits['BIN == 1'],
 2: catalog.fits['BIN == 2'],
 3: catalog.fits['BIN == 3'],
 4: catalog.fits['BIN == 4'],
 5: catalog.fits['BIN == 5'],
 6: catalog.fits['BIN == 6']}

Catalogue views (i.e., selections) have their own individual visibility, and inherit the visibility of the base catalogue by default.

Instead of the [] syntax, views can also be created using the catalog.where() method, where the visibility can be given as a parameter.

catalogs[3].visibility is vmap

True

Note that catalogs is a simple mapping of integer bin IDs to instances of type Catalog. You can create such a mapping in any way you like; the values do not have to come from the same FitsCatalog, or even the same type of catalogue.

Playground

Return the number of rows in tomographic bin 5.

# need to sum because FITS cannot tell the size a priori
nrows = sum(page.size for page in catalogs[5])

print(f"tomographic bin 5 contains {nrows:_} rows")

tomographic bin 5 contains 12_752_343 rows

Fields

To turn catalogues into spectra, Heracles requires a so-called mapper object that knows how to turn positions and values into spherical functions. Here, we construct a HealpixMapper instance with our desired nside and lmax parameters. When computing spherical harmonic coefficients, the HEALPix mapper will also remove the pixel window function, unless deconvolve=False is passed.

mapper = heracles.healpy.HealpixMapper(nside, lmax)

To specify the fields we wish to analyse, we construct a dictionary of keys and field definitions. Each field receives a mapper and a list of columns that it reads, plus potentially some other options.

For a standard 3×2pt analysis in harmonic space, we need

• A position field ("POS") for angular clustering and galaxy-galaxy lensing;

• A shear field ("SHE") for cosmic shear and galaxy-galaxy lensing.

When passing the column names, we could specify that the shear field should flip the sign of a column by adding a minus sign (e.g., "-E2" to flip the second shear component). However, we do not need to do that here.

Finally, we define the optional names of the masks ("VIS", "WHT") of the fields. These will be used further down below to compute mixing matrices.

fields = {
    "POS": heracles.Positions(
        mapper,
        "RA",
        "DEC",
        mask="VIS",
    ),
    "SHE": heracles.Shears(
        mapper,
        "RA",
        "DEC",
        "E1",
        "E2",
        "W",
        mask="WHT",
    ),
}

Mapping

The next step is to map the catalogues to position and shear fields for each tomographic bin. Heracles can map a set of catalogues all at once, using the map_catalogs() functions. We only need to pass in the fields and catalogues constructed previously.

If there is enough memory to hold all maps in memory at the same time, we can use parallel=True to construct maps for all entries in catalogs (i.e., tomographic bins) at the same time. Here, this has the advantage of reading the entire FITS file only once.

with Progress("mapping") as progress:
    data = heracles.map_catalogs(fields, catalogs, parallel=True, progress=progress)

/Users/ntessore/code/heracles-ec/heracles/heracles/fields.py:302: UserWarning: positions and visibility have different size
  warnings.warn("positions and visibility have different size")

The resulting data object is a dictionary with keys corresponding to each field ("POS", "SHE") and catalogue (0, 1, …). Results from Heracles are always of this form.

list(data.keys())[:6] + ["..."]

[('POS', 1), ('SHE', 1), ('POS', 2), ('SHE', 2), ('POS', 3), ('SHE', 3), '...']

Playground

Since the HEALPix mapper produces actual maps, we can take a quick look at the data for a tomographic bin ID, say 3.

i = 3

fig, ax = plt.subplots(2, 2, figsize=(12, 6))
fig.tight_layout()
plt.sca(ax[0, 0])
hp.cartview(
    data["POS", i],
    title="P",
    cmap="binary",
    min=-1.0,
    max=4.0,
    hold=True,
    lonra=[329.5, 53.5],
    latra=[-45.6, -15.7],
)
ax[0, 1].axis("off")
plt.sca(ax[1, 0])
hp.cartview(
    data["SHE", i][0],
    title="G1",
    cmap="RdGy",
    min=-0.5,
    max=0.5,
    hold=True,
    lonra=[329.5, 53.5],
    latra=[-45.6, -15.7],
)
plt.sca(ax[1, 1])
hp.cartview(
    data["SHE", i][1],
    title="G2",
    cmap="RdGy",
    min=-0.5,
    max=0.5,
    hold=True,
    lonra=[329.5, 53.5],
    latra=[-45.6, -15.7],
)
plt.show()

Alms

Since we are working in harmonic space, the real-space maps we created so far are not what we are actually after. We therefore transform the data into harmonic space (i.e. \(a_{lm}\)) using the transform() function.

with Progress("transform") as progress:
    alms = heracles.transform(fields, data, progress=progress)

The resulting alms dictionary has the same keys as data, but contains the spherical harmonic coefficients of the maps.

list(alms.keys())[:6] + ["..."]

[('POS', 1), ('SHE', 1), ('POS', 2), ('SHE', 2), ('POS', 3), ('SHE', 3), '...']

Two-point statistics

We are now able to compute two-point statistics in the form of angular power spectra. To do so, we simply call the angular_power_spectra() function on alms. The function will automatically remove the noise bias from the spectra (unless debias=False is passed), and can optionally compute binned spectra (by passing the bins= and weights= parameters).

cls = heracles.angular_power_spectra(alms)

This computes the angular power spectra for all alms combinations. They are very many, arranged into the familiar dictionary format with entries such as ("POS", "SHE", 5, 4) for the position and shear cross-spectrum for bins 5 and 4 (in order).

list(cls.keys())[:6] + ["..."] + list(cls.keys())[-6:]

[('POS', 'POS', 1, 1),
 ('POS', 'SHE', 1, 1),
 ('POS', 'POS', 1, 2),
 ('POS', 'SHE', 1, 2),
 ('POS', 'POS', 1, 3),
 ('POS', 'SHE', 1, 3),
 '...',
 ('SHE', 'SHE', 5, 5),
 ('POS', 'SHE', 6, 5),
 ('SHE', 'SHE', 5, 6),
 ('POS', 'POS', 6, 6),
 ('POS', 'SHE', 6, 6),
 ('SHE', 'SHE', 6, 6)]

Mixing matrices

We usually do not have data over the entire sky, but are limited by the survey footprint and visibility. This affects the harmonic-space two-point statistics, and is modelled by the so-called mixing matrices.

Computing the mixing matrices requires additional computation for the visibility and shear weights in each tomographic bin. Square mixing matrices up to lmax require input spectra with 2x lmax, for which we create a new HEALPix mapper at twice the resolution.

mapper2 = heracles.healpy.HealpixMapper(2 * nside, 2 * lmax)

We now use the same catalogues to map the visibility and weights, using the map_catalogs() function as before.

Since each field gets its own mapper, we could also have computed all maps in one go.

# visibility maps are taken as-is from catalogue, so no columns
fields2 = {
    "VIS": heracles.Visibility(
        mapper2,
    ),
    "WHT": heracles.Weights(
        mapper2,
        "RA",
        "DEC",
        "W",
    ),
}

with Progress("mapping") as progress:
    data2 = heracles.map_catalogs(fields2, catalogs, parallel=True, progress=progress)

The data2 output has the same format as earlier, but now contains "VIS" and "WHT" maps as expected.

list(data2.keys())[:6] + ["..."]

[('VIS', 1), ('WHT', 1), ('VIS', 2), ('WHT', 2), ('VIS', 3), ('WHT', 3), '...']

Next we transform the maps for the mixing matrices …

with Progress("transform") as progress:
    alms2 = heracles.transform(fields2, data2, progress=progress)

… and compute their angular power spectra.

cls2 = heracles.angular_power_spectra(alms2)

With the angular power spectra of visibility and weight maps available, we can compute all mixing matrices with the mixing_matrices() function. It can optionally compute mixing matrices for binned spectra, by providing the bins= and weights= parameters.

with Progress("mixmats") as progress:
    mms = heracles.mixing_matrices(
        fields, cls2, l1max=lmax, l2max=lmax, progress=progress
    )

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

The mixing matrices are returned in a dictionary with names matching the angular power spectra:

• ("POS", "POS", i, j) contains the \(M^{TT,TT}\) mixing matrix;

• ("POS", "SHE", i, j) contains the \(M^{TE,TE} = M^{TB,TB}\) mixing matrix;

• ("SHE", "SHE", i, j) contains the stack of mixing matrices

  • \(M^{EE,EE} = M^{BB,BB}\),

  • \(M^{EE,BB} = M^{BB,EE}\),

  • \(M^{EB,EB}\).

For more details on mixing matrices in general, see, e.g., the paper by Brown, Castro & Taylor (2005).

list(mms.keys())[:6] + ["..."] + list(mms.keys())[-6:]

[('POS', 'POS', 1, 1),
 ('POS', 'SHE', 1, 1),
 ('POS', 'POS', 1, 2),
 ('POS', 'SHE', 1, 2),
 ('POS', 'POS', 1, 3),
 ('POS', 'SHE', 1, 3),
 '...',
 ('SHE', 'SHE', 5, 5),
 ('POS', 'SHE', 6, 5),
 ('SHE', 'SHE', 5, 6),
 ('POS', 'POS', 6, 6),
 ('POS', 'SHE', 6, 6),
 ('SHE', 'SHE', 6, 6)]

Playground

As before, the mapper has produced HEALPix maps, so we can have a look at the visibility and weight map for bin ID 3, say.

i = 3

fig, ax = plt.subplots(1, 2, figsize=(12, 8))
fig.tight_layout()
plt.sca(ax[0])
hp.cartview(
    data2["VIS", i],
    title="V",
    cmap="binary",
    hold=True,
    lonra=[329.5, 53.5],
    latra=[-45.6, -15.7],
)
plt.sca(ax[1])
hp.cartview(
    data2["WHT", i],
    title="W",
    cmap="binary",
    hold=True,
    lonra=[329.5, 53.5],
    latra=[-45.6, -15.7],
)
plt.show()

And here is the corresponding ("SHE", "SHE", i, i) mixing matrix \(M^{EE,EE}\):

plt.imshow(
    mms["SHE", "SHE", i, i][0], cmap="binary", norm=mpl.colors.LogNorm(vmin=1e-7)
)
plt.colorbar(pad=0.025, fraction=0.0465)
plt.show()

Theory

To model the expected angular power spectra, we require the mixing matrices and the expected full-sky angular power spectra from theory. Here, we use CAMB to compute the latter.

import camb
from camb.sources import SplinedSourceWindow

We set the CAMB cosmology to match the simulation that created the example data set.

# cosmology for the simulation
h = 0.7
Oc = 0.25
Ob = 0.05

# set up CAMB parameters for matter angular power spectrum
pars = camb.set_params(
    H0=100 * h, omch2=Oc * h**2, ombh2=Ob * h**2, NonLinear=camb.model.NonLinear_both
)
pars.Want_CMB = False
pars.min_l = 1
pars.set_for_lmax(2 * lmax, lens_potential_accuracy=1);

We also need the redshift distributions for the tomographic bins.

with np.load("nz.npz") as npz:
    z, nz = npz["z"], npz["nz"]

Given the redshift distributions, we can construct the CAMB source distributions for positions (counts) and shears (lensing).

sources = []
for i, nz_i in enumerate(nz):
    sources += [
        SplinedSourceWindow(source_type="counts", z=z, W=nz_i),
        SplinedSourceWindow(source_type="lensing", z=z, W=nz_i),
    ]
pars.SourceWindows = sources

Use the pars we constructed above to compute the full sky theory spectra up to lmax, setting raw_cl=True to return unscaled full-sky spectra.

results = camb.get_results(pars)
camb_cls = results.get_source_cls_dict(lmax=lmax, raw_cl=True)

This is the factor needed to convert from the convergence spectra returned by CAMB to shear E-mode spectra.

ell = np.arange(lmax + 1)
fl = -np.sqrt((ell + 2) * (ell + 1) * ell * (ell - 1))
fl /= np.clip(ell * (ell + 1), 1, None)

Now we can compute the theory spectra for our observations, using the CAMB results and the mixing matrices we computed earlier. We store everything in dictionary using the same format as before.

theory = {}
for i in range(1, 7):
    for j in range(i, 7):
        # get the full-sky spectra; B-mode is assumed zero
        cl_pp = camb_cls[f"W{2 * i - 1}xW{2 * j - 1}"]
        cl_pe = fl * camb_cls[f"W{2 * i - 1}xW{2 * j}"]
        cl_pb = np.zeros_like(cl_pe)
        cl_ep = fl * camb_cls[f"W{2 * i}xW{2 * j - 1}"]
        cl_bp = np.zeros_like(cl_ep)
        cl_ee = fl**2 * camb_cls[f"W{2 * i}xW{2 * j}"]
        cl_bb = np.zeros_like(cl_ee)
        cl_eb = np.zeros_like(cl_ee)
        cl_be = np.zeros_like(cl_ee)

        # all mixing matrix combinations
        theory["POS", "POS", i, j] = mms["POS", "POS", i, j] @ cl_pp
        theory["POS", "SHE", i, j] = np.array(
            [
                mms["POS", "SHE", i, j] @ cl_pe,
                mms["POS", "SHE", i, j] @ cl_pb,
            ]
        )
        theory["POS", "SHE", j, i] = np.array(
            [
                mms["POS", "SHE", j, i] @ cl_ep,
                mms["POS", "SHE", j, i] @ cl_bp,
            ]
        )
        theory["SHE", "SHE", i, j] = np.array(
            [
                [
                    mms["SHE", "SHE", i, j][0] @ cl_ee
                    + mms["SHE", "SHE", i, j][1] @ cl_bb,
                    mms["SHE", "SHE", i, j][2] @ cl_eb,
                ],
                [
                    mms["SHE", "SHE", i, j][2] @ cl_be,
                    mms["SHE", "SHE", i, j][0] @ cl_bb
                    + mms["SHE", "SHE", i, j][1] @ cl_ee,
                ],
            ]
        )

Results

Finally, we can plot the results for positions and shears.

ell = np.arange(lmax + 1)

fig, ax = plt.subplots(6, 6, figsize=(6, 6), sharex=True, sharey=True)

for i in range(1, 7):
    for j in range(1, i):
        ax[j - 1, i - 1].axis("off")
    for j in range(i, 7):
        ax[j - 1, i - 1].plot(
            ell[1:], cls["POS", "POS", i, j][1:], c="C0", lw=1.5, zorder=3.0, alpha=0.5
        )
        ax[j - 1, i - 1].plot(
            ell[1:], theory["POS", "POS", i, j][1:], c="C0", lw=1.0, zorder=4.0
        )
        ax[j - 1, i - 1].axhline(0.0, c="k", lw=0.8, zorder=-1)
        ax[j - 1, i - 1].tick_params(axis="both", which="both", direction="in")

ax[0, 0].set_xscale("log")
ax[0, 0].set_xlim(1 / 3, lmax * 3)
ax[0, 0].xaxis.get_major_locator().set_params(numticks=99)
ax[0, 0].xaxis.get_minor_locator().set_params(
    numticks=99, subs=np.arange(0.1, 1.0, 0.1)
)
ax[0, 0].set_yscale(
    "symlog", linthresh=1e-7, linscale=0.45, subs=np.arange(0.1, 1.0, 0.1)
)
ax[0, 0].set_ylim(-2e-7, 2e-6)

fig.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.0, wspace=0.0, hspace=0.0)

fig.supxlabel("angular mode $\\ell$", y=-0.05, va="top")
fig.supylabel("galaxy clustering $C_\\ell$", x=-0.1, ha="right")

plt.show()

fig, ax = plt.subplots(6, 6, figsize=(6, 6), sharex=True, sharey=True)

for i in range(1, 7):
    for j in range(1, i):
        ax[j - 1, i - 1].axis("off")
    for j in range(i, 7):
        ax[j - 1, i - 1].plot(
            ell[2:],
            cls["SHE", "SHE", i, j][0, 0, 2:],
            c="C0",
            lw=1.5,
            zorder=3.0,
            alpha=0.5,
        )
        ax[j - 1, i - 1].plot(
            ell[2:], theory["SHE", "SHE", i, j][0, 0, 2:], c="C0", lw=1.0, zorder=4.0
        )
        ax[j - 1, i - 1].plot(
            ell[2:],
            cls["SHE", "SHE", i, j][1, 1, 2:],
            c="C1",
            lw=1.5,
            zorder=1.0,
            alpha=0.5,
        )
        ax[j - 1, i - 1].plot(
            ell[2:], theory["SHE", "SHE", i, j][1, 1, 2:], c="C1", lw=1.0, zorder=2.0
        )
        ax[j - 1, i - 1].axhline(0.0, c="k", lw=0.8, zorder=-1)
        ax[j - 1, i - 1].tick_params(axis="both", which="both", direction="in")

ax[0, 0].set_xscale("log")
ax[0, 0].set_xlim(1 / 3, lmax * 3)
ax[0, 0].xaxis.get_major_locator().set_params(numticks=99)
ax[0, 0].xaxis.get_minor_locator().set_params(
    numticks=99, subs=np.arange(0.1, 1.0, 0.1)
)
ax[0, 0].set_yscale(
    "symlog", linthresh=1e-10, linscale=0.45, subs=np.arange(0.1, 1.0, 0.1)
)
ax[0, 0].set_ylim(-3e-10, 5e-9)

fig.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.0, wspace=0.0, hspace=0.0)

fig.supxlabel("angular mode $\\ell$", y=-0.05, va="top")
fig.supylabel("cosmic shear $C_\\ell$", x=-0.1, ha="right")

plt.show()

fig, ax = plt.subplots(6, 6, figsize=(6, 6), sharex=True, sharey=True)

for i in range(1, 7):
    for j in range(1, 7):
        ax[j - 1, i - 1].plot(
            ell[2:],
            cls["POS", "SHE", i, j][0, 2:],
            c="C0",
            lw=1.5,
            zorder=3.0,
            alpha=0.5,
        )
        ax[j - 1, i - 1].plot(
            ell[2:], theory["POS", "SHE", i, j][0, 2:], c="C0", lw=1.0, zorder=4.0
        )
        ax[j - 1, i - 1].plot(
            ell[2:],
            cls["POS", "SHE", i, j][1, 2:],
            c="C1",
            lw=1.5,
            zorder=1.0,
            alpha=0.5,
        )
        ax[j - 1, i - 1].plot(
            ell[2:], theory["POS", "SHE", i, j][1, 2:], c="C1", lw=1.0, zorder=2.0
        )
        ax[j - 1, i - 1].axhline(0.0, c="k", lw=0.8, zorder=-1)
        ax[j - 1, i - 1].tick_params(axis="both", which="both", direction="in")

ax[0, 0].set_xscale("log")
ax[0, 0].set_xlim(1 / 3, lmax * 3)
ax[0, 0].xaxis.get_major_locator().set_params(numticks=99)
ax[0, 0].xaxis.get_minor_locator().set_params(
    numticks=99, subs=np.arange(0.1, 1.0, 0.1)
)
ax[0, 0].set_yscale(
    "symlog", linthresh=1e-9, linscale=0.45, subs=np.arange(0.1, 1.0, 0.1)
)
ax[0, 0].set_ylim(-8e-8, 4e-8)

fig.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.0, wspace=0.0, hspace=0.0)

fig.supxlabel("angular mode $\\ell$", y=-0.05, va="top")
fig.supylabel("galaxy--galaxy lensing $C_\\ell$", x=-0.1, ha="right")

plt.show()

Output

Additionally, we can write the data we produced to file.

heracles.write_maps("example-data_maps.fits", data, clobber=True)
heracles.write("example-mask-spectra.fits", cls2, clobber=True)
heracles.write("example-theory.fits", theory, clobber=True)
heracles.write("example-spectra.fits", cls, clobber=True)
heracles.write("example-mixmats.fits", mms, clobber=True)