Overview

Introduction

The hmpdf code computes one- and two-point PDFs, the covariance matrix of the one-point PDF, and the angular power spectrum/correlation function of cosmological fields.

The formalism is based on the halo model, thus, only fields that are mainly sourced by virialized matter make sense.

Currently, two fields are supported:

• thermal Sunyaev-Zel'dovich effect (Compton-y)
• weak lensing convergence

The code is interfaced both through C header files as well as through a convenient python wrapper.

Applications of the code can be found in

• Thiele, Hill, Smith 2019 (for the tSZ field)
• Thiele, Hill, Smith 2020 (for the weak lensing convergence field)

which we also request you cite if you use the code.

Building

Dependencies:

• FFTW3
• GSL (>= 2.4)
• CLASS (as of writing, the latest version has a bug that can occasionally cause segfaults; the fix is simple)

Compilation should be straightforward. You may have to edit the following parts of the Makefile:

• adapt PATHTOCLASS to your CLASS installation
• add FFTW3/GSL locations to LINKER if those are non-standard
• remove the OpenMP flag if you do not wish parallel execution

Then you can simply type make.

Afterwards, if you intend to use the python wrapper, type make python. This will run pip. It will also hardcode the location of libhmpdf.so into the python package, which has the advantage that you can recompile the C code without having to run pip again, and the disadvantage that you need to keep libhmpdf.so where it is.

It may also be convenient to copy libhmpdf.so into one of the default locations searched by the linker.

Usage

The C interface can be found in hmpdf.h. All exposed objects' names start with hmpdf_ to help you keep your name space clean.

Workflow is as follows:

1. allocate a new hmpdf_obj with hmpdf_new().
2. set all options (required and optional) with hmpdf_init(), this will also compute the data needed for all outputs. See hmpdf_configs_e for optional inputs.
3. get your output [hmpdf_get_op(), hmpdf_get_tp(), hmpdf_get_cov(), hmpdf_get_Cell(), hmpdf_get_Cphi(), hmpdf_get_map(), hmpdf_get_map_op()].
4. go to (3.) if you require any other outputs; go to (2.) if you want to re-run the code with different options.
5. free the memory associated with the hmpdf_obj with hmpdf_delete().

There is also a Python wrapper (hmpdf.py), with analogous workflow and simpler syntax.

Options

There is a number of optional arguments to hmpdf_init(), the most important ones to get you started are

• pixelization: hmpdf_pixel_side
• ell space filters: hmpdf_tophat_radius, hmpdf_gaussian_fwhm, hmpdf_custom_ell_filter (and hmpdf_custom_ell_filter_params)
• Gaussian noise: hmpdf_noise_pwr (and hmpdf_noise_pwr_params)
• multithreading: hmpdf_N_threads

See the documentation of hmpdf_configs_e for all options.

Runtime

The following were found with default settings on a laptop. These functions scale as hmpdf_N_M x hmpdf_N_z, this is omitted in the following:

• hmpdf_init(): a few seconds (depends strongly on the configuration)

  scales as hmpdf_N_theta to hmpdf_N_theta^2. Different if CLASS runtime starts to dominate. hmpdf_tsz takes longer than hmpdf_kappa.

• hmpdf_get_tp(): ~10 seconds on a single thread.

  scales as hmpdf_N_signal^2.

• hmpdf_get_cov(): ~20 minutes on 8 threads.

  scales as hmpdf_N_signal^2 x hmpdf_N_phi.

The simplified simulations (maps) scale as hmpdf_map_fsky / hmpdf_pixel_side^2.

Other functions are fast in comparison to hmpdf_init(). hmpdf_init(), hmpdf_get_cov(), hmpdf_get_map(), hmpdf_get_map_op() are parallelized in critical parts, while hmpdf_get_tp() does not get faster if hmpdf_N_threads is increased. The simplified simulations can easily become memory throughput-limited, in which case speed does not scale well with hmpdf_N_threads.

Error handling

All functions [except hmpdf_new()] return a non-zero int if an error occured. A traceback will be printed to stderr. You can remove the definition of the macro DEBUG in the Makefile, in which case no error handling at all will be happening (could give a marginal speed-up in some cases, but not recommended).

Thread safety

All functions are not threadsafe: making two calls on the same hmpdf_obj concurrently results in undefined behaviour.

Examples

Some examples are collected in examples/example.c and examples/example.py. You can compile the C code with

/* gcc -I../include -o example example.c -L.. -lhmpdf */

A simple calculation of the weak lensing convergence one-point PDF with all settings at default would look like this:

int example_kappa_onepoint(void)
{
    /* construct some binedges */
    int Nbins = 100; double kappamin = 0.0; double kappamax = 0.3;
    double binedges[Nbins+1];
    for (int ii=0; ii<=Nbins; ii++)
        binedges[ii] = kappamin + (double)(ii)*(kappamax-kappamin)/(double)(Nbins);
 
    /* get a new hmpdf_obj */
    hmpdf_obj *d = hmpdf_new();
    if (!(d))
        return -1;
 
    /* initialize with default settings */
    if (hmpdf_init(d, "example.ini", hmpdf_kappa,
                   1.0/* source redshift */))
        return -1;
 
    /* get the one-point PDF */
    double op[Nbins];
    if (hmpdf_get_op(d, Nbins, binedges, op,
                     1/* include two-halo term */,
                     0/* don't include noise */))
        return -1;
 
    /* free memory associated with the hmpdf_obj */
    if (hmpdf_delete(d))
        return -1;
 
    /* do something with the one-point PDF ... */
 
    return 0;
}

The same thing using the python wrapper is really easy:

def example_kappa_onepoint() :
    binedges = np.linspace(0.0, 0.3, num=101)
    with HMPDF() as d :
        d.init('example.ini', 'kappa', 1.0)
        op = d.get_op(binedges)
    # do something with the one-point PDF ...

Including the effect of an ell-space filter would look like this:

int example_ell_filter_use(void)
{
    hmpdf_obj *d = hmpdf_new();
    if (!(d))
        return -1;
 
    /* initialize with an ell-space filter which cuts off all
     * modes above ell = 5000
     */
    double ell_max = 5000.0;
    if (hmpdf_init(d, "example.ini", hmpdf_kappa, 1.0,
                   hmpdf_custom_ell_filter, &example_ell_filter,
                   hmpdf_custom_ell_filter_params, &ell_max))
        return -1;
 
    /* get your results ... */
 
    if (hmpdf_delete(d))
        return -1;
 
    return 0;
}

Here, we defined the function example_ell_filter conforming to the typedef hmpdf_ell_filter_f:

double example_ell_filter(double ell, void *p)
{
    double ell_max = *(double *)p;
    return (double)(ell < ell_max);
}

We can also do this in python (note, however, that due to the overhead in calling the python lambda this will be slower than in C):

def example_ell_filter_use() :
    example_ell_filter = lambda ell : float(ell < 5000.0)
    with HMPDF() as d :
        d.init('example.ini', 'kappa', 1.0,
               custom_ell_filter=example_ell_filter)
        # get your results ...

For a more involved example, in examples/WL_PDF_forecast/ you can find all ingredients to reproduce the Fisher forecast from the WL convergence PDF.