CosmoSIS module for window function Cl bias

window.py

"""
Module to apply window function matrix to Cls

Compatible with matrices generated by https://github.com/ntessore/clwin.

Configuration:

Add an entry of the form `m_<spectrum> = <matrix.fits>` for each spectrum that
needs to be modified.

Additionally, set `use_fsky = <bool>` to apply rescaling by the `f_sky` value
from the matrix FITS file.

Example:

    use_fsky = true
    m_shear_cl = survey_window_spin2.fits

"""
from __future__ import print_function

from builtins import range
from builtins import object
from cosmosis.datablock import option_section, names
import numpy as np
from astropy.io import fits

def setup(options):
    mat = {}

    use_fsky = options.get_bool(option_section, 'use_fsky', False)

    print('apply f_sky correction:', use_fsky)

    print('loading matrices for spectra:')
    for _, key in options.keys(option_section):
        if key.startswith('m_'):
            spec = key[2:]
            fil = options.get_string(option_section, key)
            print('    -  %s: %s' % (spec, fil))
            mat[spec] = fits.getdata(fil)
            if use_fsky:
                fsky = fits.getval(fil, 'FSKY')
                mat[spec] *= 1./fsky

    return mat

def execute(block, config):
    print('applying matrices to spectra:')
    for spec in config:
        mat = config[spec]

        # lmax for input and output spectrum
        lm0 = mat.shape[0]-1
        lm1 = mat.shape[1]-1
        
        # ell values (integers) for input and output spectrum
        el0 = np.arange(lm0+1, dtype=float)
        el1 = np.arange(lm1+1, dtype=float)

        print('- %s: [0, %d] -> [0, %d]' % (spec, lm1, lm0))

        # input and output section names
        input_section = spec
        output_section = spec + '_w'

        # number of bins in spectrum
        nba = block[input_section, 'nbin_a']
        nbb = block[input_section, 'nbin_b']

        # input ell values (non-integer)
        el2 = block[input_section, 'ell']

        # output ell values (integer)
        block[output_section, 'ell'] = el0

        for b1 in range(1, nba+1):
            for b2 in range(1, nbb+1):
                name = 'bin_{}_{}'.format(b1, b2)
                if not block.has_value(input_section, name):
                    continue

                # input power spectrum (non-integer ells)
                cl2 = block[input_section, name]

                # interpolate integer input ells
                cl1 = np.interp(el1, el2, cl2)

                # apply matrix
                cl0 = np.dot(mat, cl1)

                # output power spectrum (integer ells)
                block[output_section, name] = cl0

    return 0