process_data

multiprobe_framework.process_data

get_data_boss(path2boss, selection)

Get the BOSS data for the given selection.

Parameters:

Name	Type	Description	Default
path2boss		the path to the BOSS data	required
selection		the bin of galaxies either 'cmass' or 'lowz'	required

Returns:

Type	Description
	the right ascension, declination, redshift, and total weight of the galaxies

Source code in src/multiprobe_framework/process_data.py

def get_data_boss(path2boss, selection):
    """
    Get the BOSS data for the given selection.

    :param path2boss: the path to the BOSS data
    :param selection: the bin of galaxies either 'cmass' or 'lowz'

    :return: the right ascension, declination, redshift, and total weight of
             the galaxies
    """
    if selection == "lowz":
        z_min, z_max = 0.0, 0.4
    elif selection == "cmass":
        z_min, z_max = 0.4, 0.8
    else:
        raise ValueError("Choose either 'cmass' or 'lowz' for galaxy selections")

    all_data = {
        "ra": [],
        "dec": [],
        "z": [],
        "w_cp": [],
        "w_noz": [],
        "w_star": [],
        "w_seeing": [],
        "w_systot": [],
    }

    for region in ["North", "South"]:
        df = f"{path2boss}/galaxy_DR12v5_{selection.upper()}_{region}.fits.gz"
        with fits.open(df) as hdul:
            data = hdul[1].data
            useful = np.logical_and(data["Z"] >= z_min, data["Z"] <= z_max)

            all_data["ra"].append(data["RA"][useful])
            all_data["dec"].append(data["DEC"][useful])
            all_data["z"].append(data["Z"][useful])
            all_data["w_cp"].append(data["WEIGHT_CP"][useful])
            all_data["w_noz"].append(data["WEIGHT_NOZ"][useful])
            all_data["w_star"].append(data["WEIGHT_STAR"][useful])
            all_data["w_seeing"].append(data["WEIGHT_SEEING"][useful])
            all_data["w_systot"].append(data["WEIGHT_SYSTOT"][useful])

    concatenated_data = {key: np.concatenate(val) for key, val in all_data.items()}
    w_tot = (
        concatenated_data["w_star"]
        * concatenated_data["w_seeing"]
        * (concatenated_data["w_cp"] + concatenated_data["w_noz"] - 1.0)
    )

    return (
        concatenated_data["ra"],
        concatenated_data["dec"],
        concatenated_data["z"],
        w_tot,
    )

generate_boss_map(ra, dec=None, galaxy_weight=None, boss_binary_mask=None, nside=None, mode='regular', rot=True, seed_no=None)

Generates a galaxy overdensity map with options for rotation into galactic coordinates and creating a random map.

Parameters:

Name	Type	Description	Default
ra		the right ascension of the galaxies	required
dec		the declination of the galaxies	None
galaxy_weight		the weight of the galaxies	None
boss_binary_mask		the binary mask at the correct NSIDE	None
nside		the nside of the output map to be generated	None
mode		the mode of the map generation, either 'regular' or 'random'	'regular'
rot		whether to rotate the coordinates into galactic coordinates	True
seed_no		the seed number for random mode	None

Returns:

Type	Description
	the galaxy overdensity map generated from the input data

Source code in src/multiprobe_framework/process_data.py

def generate_boss_map(
    ra,
    dec=None,
    galaxy_weight=None,
    boss_binary_mask=None,
    nside=None,
    mode="regular",
    rot=True,
    seed_no=None,
):
    """
    Generates a galaxy overdensity map with options for rotation into galactic
    coordinates and creating a random map.

    :param ra: the right ascension of the galaxies
    :param dec: the declination of the galaxies
    :param galaxy_weight: the weight of the galaxies
    :param boss_binary_mask: the binary mask at the correct NSIDE
    :param nside: the nside of the output map to be generated
    :param mode: the mode of the map generation, either 'regular' or 'random'
    :param rot: whether to rotate the coordinates into galactic coordinates
    :param seed_no: the seed number for random mode

    :return: the galaxy overdensity map generated from the input data
    """

    if boss_binary_mask is None or nside is None or galaxy_weight is None:
        raise ValueError("boss_binary_mask, nside, and galaxy weight must be provided.")

    # Initialize deltag map and identify valid and masked regions based on
    # Pixel completeness following
    # https://arxiv.org/pdf/1809.07204 C_pix.0.8
    # ToDo verify this BOSS data handling against CGG Cls
    deltag = np.zeros(hp.nside2npix(nside))
    valid_idx = np.where(boss_binary_mask == 1)[0]
    masked_idx = np.where(boss_binary_mask == 0)[0]

    if mode == "regular":
        # Rotate coordinates if required and get pixel indices
        if rot and dec is not None:
            r = hp.Rotator(coord=["C", "G"])
            ra, dec = r(ra, dec, lonlat=True)
        pix = hp.ang2pix(nside, ra, dec, lonlat=True)
    elif mode == "random":
        if seed_no is None:
            raise ValueError("Seed number must be provided for random mode.")
        np.random.seed(seed_no)
        # Ensure random selection is only from valid indices
        pix = np.random.choice(valid_idx, size=ra.size)

    # Accumulate weights into the pixels
    for i, p in enumerate(pix):
        deltag[p] += galaxy_weight[i]

    # Normalize and mask the deltag map
    mean_valid = np.mean(deltag[valid_idx])
    deltag[valid_idx] = deltag[valid_idx] / mean_valid - 1 if mean_valid else 0
    deltag[masked_idx] = 0.0

    return deltag

generate_kids_map(nside, e1, e2, w, idx, rotate=False, return_w2_sigma2=False, m=0.0, seed_no=None)

Generate a shear component maps from the given ellipticities and weights of the galaxies.

Parameters:

Name	Type	Description	Default
nside		the nside of the map to be generated	required
e1		the first component of the galaxy ellipticity	required
e2		the second component of the galaxy ellipticity	required
w		the weight of the galaxies	required
idx		the pixel indices of the galaxies	required
rotate		whether to rotate the ellipticities randomly for a noise realization	False
return_w2_sigma2		whether to return the weight squared sigma squared map	False
m		the multiplicative bias	0.0
seed_no		the seed number for random mode	None

Returns:

Type	Description
	the galaxy ellipticity map and the weight map

Source code in src/multiprobe_framework/process_data.py

def generate_kids_map(
    nside, e1, e2, w, idx, rotate=False, return_w2_sigma2=False, m=0.0, seed_no=None
):
    """
    Generate a shear component maps from the given ellipticities and weights of
    the galaxies.

    :param nside: the nside of the map to be generated
    :param e1: the first component of the galaxy ellipticity
    :param e2: the second component of the galaxy ellipticity
    :param w: the weight of the galaxies
    :param idx: the pixel indices of the galaxies
    :param rotate: whether to rotate the ellipticities randomly for a noise realization
    :param return_w2_sigma2: whether to return the weight squared sigma squared map
    :param m: the multiplicative bias
    :param seed_no: the seed number for random mode

    :return: the galaxy ellipticity map and the weight map
    """

    npix = hp.nside2npix(nside)
    e1 /= 1 + m
    e2 /= 1 + m

    if rotate:
        np.random.seed(seed_no)
        alpha = np.pi * np.random.rand(len(e1))
        e = np.sqrt(e1**2 + e2**2)
        e1 = np.cos(2.0 * alpha) * e
        e2 = np.sin(2.0 * alpha) * e

    e1_map = np.bincount(idx, weights=w * e1, minlength=npix)
    e2_map = np.bincount(idx, weights=w * e2, minlength=npix)
    w_map = np.bincount(idx, weights=w, minlength=npix)

    good_pixel = w_map > 0
    e1_map[good_pixel] /= w_map[good_pixel]
    e2_map[good_pixel] /= w_map[good_pixel]

    if return_w2_sigma2:
        w2_sigma2_map = np.bincount(
            idx, weights=w**2 * (e1**2 + e2**2) / 2, minlength=npix
        )
        return e1_map, e2_map, w_map, w2_sigma2_map

    return e1_map, e2_map, w_map

get_data_kids(path_kids, tomo_bin, nside)

Get the KiDS data for the given redshift range.

Parameters:

Name	Type	Description	Default
path_kids		the path to the KiDS data	required

Returns:

Type	Description

Source code in src/multiprobe_framework/process_data.py

def get_data_kids(path_kids, tomo_bin, nside):
    """
    Get the KiDS data for the given redshift range.

    :param path_kids: the path to the KiDS data

    :return:
    """

    # mbias bestfits following the KiDS1000xtSZ pipeline:
    # (https://github.com/tilmantroester/KiDS-1000xtSZ)
    m_bias = [-0.009, -0.011, -0.015, 0.002, 0.007]

    # The cuts in redshift
    Z_CUTS = [(0.1, 0.3), (0.3, 0.5), (0.5, 0.7), (0.7, 0.9), (0.9, 1.2)]

    # for this tomographic bin get the correct data

    z_cut = Z_CUTS[tomo_bin - 1]
    path2kids = os.path.join(
        path_kids, f"KiDS-1000_z{z_cut[0]}-{z_cut[1]}nside_{nside}_galactic.npz"
    )
    with np.load(path2kids) as data:
        w = data["w"]
        e1 = data["e1"]
        e2 = data["e2"]
        pixel_idx = data["pixel_idx"]

    return e1, e2, w, m_bias[tomo_bin - 1], pixel_idx

prepare_catalog_kids(catalog_filename, column_names={'x': 'x', 'y': 'y', 'e1': 'e1', 'e2': 'e2', 'w': 'w', 'm': 'm', 'c1': 'c1', 'c2': 'c2'}, c_correction='data', m_correction='catalog', z_min=None, z_max=None, selections=[('weight', 'gt', 0.0)], hdu_idx=1)

Prepare lensing catalogues.

This function a a few of the other functions related to handling the KiDS data are shamelessly copied from the excellent public repository (https://github.com/tilmantroester/KiDS-1000xtSZ) associated with the paper: https://arxiv.org/abs/2109.04458

Parameters:

Name	Type	Description	Default
catalog_filename		the path to the catalog file	required
column_names		the column names in the catalog	{'x': 'x', 'y': 'y', 'e1': 'e1', 'e2': 'e2', 'w': 'w', 'm': 'm', 'c1': 'c1', 'c2': 'c2'}
c_correction		the c correction to be applied	'data'
m_correction		the m correction to be applied	'catalog'
z_min		the minimum redshift	None
z_max		the maximum redshift	None
selections		the selections to be applied	[('weight', 'gt', 0.0)]
hdu_idx		the index of the HDU in the FITS file	1

Returns:

Type	Description
	the right ascension, declination, ellipticity components, weight, and multiplicative bias of the galaxies

Source code in src/multiprobe_framework/process_data.py

def prepare_catalog_kids(
    catalog_filename,
    column_names={
        "x": "x",
        "y": "y",
        "e1": "e1",
        "e2": "e2",
        "w": "w",
        "m": "m",
        "c1": "c1",
        "c2": "c2",
    },
    c_correction="data",
    m_correction="catalog",
    z_min=None,
    z_max=None,
    selections=[("weight", "gt", 0.0)],
    hdu_idx=1,
):
    """
    Prepare lensing catalogues.

    This function a a few of the other functions related to handling the KiDS
    data are shamelessly copied from the excellent public repository
    (https://github.com/tilmantroester/KiDS-1000xtSZ) associated with the
    paper: https://arxiv.org/abs/2109.04458

    :param catalog_filename: the path to the catalog file
    :param column_names: the column names in the catalog
    :param c_correction: the c correction to be applied
    :param m_correction: the m correction to be applied
    :param z_min: the minimum redshift
    :param z_max: the maximum redshift
    :param selections: the selections to be applied
    :param hdu_idx: the index of the HDU in the FITS file

    :return: the right ascension, declination, ellipticity components, weight,
             and multiplicative bias of the galaxies
    """
    hdu_idx = 1
    hdu = fits.open(catalog_filename)

    mask = np.ones(hdu[hdu_idx].data.size, dtype=bool)
    # Apply selections to catalog
    for col, op, val in selections:
        if op == "eq":
            mask = np.logical_and(mask, hdu[hdu_idx].data[col] == val)
        elif op == "neq":
            mask = np.logical_and(mask, hdu[hdu_idx].data[col] != val)
        elif op == "gt":
            mask = np.logical_and(mask, hdu[hdu_idx].data[col] > val)
        elif op == "ge":
            mask = np.logical_and(mask, hdu[hdu_idx].data[col] >= val)
        elif op == "lt":
            mask = np.logical_and(mask, hdu[hdu_idx].data[col] < val)
        elif op == "le":
            mask = np.logical_and(mask, hdu[hdu_idx].data[col] <= val)
        else:
            raise ValueError(f"Operator {op} not supported.")

    # For convenience, redshift cuts can be applied through the arguments
    # z_min and z_max as well.
    if z_min is not None and z_max is not None:
        z = hdu[hdu_idx].data[column_names["z"]]
        mask = np.logical_and(mask, np.logical_and(z > z_min, z <= z_max))

    ra = hdu[hdu_idx].data[column_names["x"]][mask]
    dec = hdu[hdu_idx].data[column_names["y"]][mask]
    w = hdu[hdu_idx].data[column_names["w"]][mask]
    e1 = hdu[hdu_idx].data[column_names["e1"]][mask]
    e2 = hdu[hdu_idx].data[column_names["e2"]][mask]

    # Apply c correction
    if c_correction == "catalog":
        # Use c correction supplied by the catalog
        c1 = hdu[hdu_idx].data[column_names["c1"]][mask]
        c2 = hdu[hdu_idx].data[column_names["c2"]][mask]
        c1_mask = c1 > -99
        c2_mask = c2 > -99
        e1[c1_mask] -= c1[c1_mask]
        e2[c2_mask] -= c2[c2_mask]
    elif c_correction == "data":
        # Calculate c correction from the weighted ellipticity average
        c1 = np.sum(w * e1) / np.sum(w)
        c2 = np.sum(w * e2) / np.sum(w)
        e1 -= c1
        e2 -= c2

    # Apply m correction
    if m_correction == "catalog":
        m = hdu[hdu_idx].data[column_names["m"]][mask]
    else:
        m = np.zeros_like(w)

    hdu.close()

    return ra, dec, e1, e2, w, m

get_nmt_cl_data(spec, maps1, maps2, mask1, mask2, config_dict, path2wsp, noise_realization=None, instrument_beam=None, no_pixwin=False, nside=2048)

Return the NaMaster computed Cl for a set of maps and masks.

Maps should be specified at lists, with the spin of the map implicit in the number of maps passed in the list. Noise realization: remove an average noise cl from the computed fullsky Cl.

Source code in src/multiprobe_framework/process_data.py

def get_nmt_cl_data(
    spec,
    maps1,
    maps2,
    mask1,
    mask2,
    config_dict,
    path2wsp,
    noise_realization=None,
    instrument_beam=None,
    no_pixwin=False,
    nside=2048,
):
    """
    Return the NaMaster computed Cl for a set of maps and masks.

    Maps should be specified at lists, with the spin of the map implicit in the
    number of maps passed in the list. Noise realization: remove an average
    noise cl from the computed fullsky Cl.
    """

    beam = {}
    beam[0] = hp.pixwin(nside=nside, pol=True)[0]
    beam[2] = hp.pixwin(nside=nside, pol=True)[1]

    if no_pixwin:
        beam[0] = None
        beam[2] = None

    spin_1 = utils.get_spin(maps1)
    spin_2 = utils.get_spin(maps2)

    f_0 = nmt.NmtField(mask1, maps1, beam=beam[spin_1], n_iter=3)
    f_1 = nmt.NmtField(mask2, maps2, beam=beam[spin_2], n_iter=3)

    if instrument_beam is not None:
        if spec[0] == "t":
            ins_beam_shape = len(instrument_beam)
            instrument_beam = np.pad(
                instrument_beam, (0, 3 * nside - ins_beam_shape), "constant"
            )  # pad to correct size
            beam = (
                beam[spin_1] * instrument_beam
                if beam[spin_1] is not None
                else instrument_beam
            )
            f_0 = nmt.NmtField(mask1, maps1, beam=beam)

        if spec[1] == "t":
            ins_beam_shape = len(instrument_beam)
            instrument_beam = np.pad(
                instrument_beam, (0, 3 * nside - ins_beam_shape), "constant"
            )  # pad to correct size
            beam = (
                beam[spin_1] * instrument_beam
                if beam[spin_1] is not None
                else instrument_beam
            )
            f_1 = nmt.NmtField(mask2, maps2, beam=beam)

    # set up the Namaster workspace
    w = nmt.NmtWorkspace()
    w.read_from(path2wsp + f"/pymaster_workspace_{spec}.fits")

    cl_coupled = nmt.compute_coupled_cell(f_0, f_1)

    if noise_realization is not None:
        cl_compute = w.decouple_cell(cl_coupled, cl_noise=noise_realization)

    else:
        cl_compute = w.decouple_cell(cl_coupled)

    if (spin_1 == 0) & (spin_2 == 0):
        print(cl_compute[0])

        assert cl_compute[0].shape[0] == config_dict[f"{spec}_n_bin"]
        return cl_compute[0]

    elif (spin_1 == 0) & (spin_2 == 2) or (spin_1 == 2) & (spin_2 == 0):
        assert cl_compute[0].shape[0] == config_dict[f"{spec}_n_bin"]

        return cl_compute[0]

    elif (spin_1 == 2) & (spin_2 == 2):
        cl_compute_ee = cl_compute[0]
        cl_compute_bb = cl_compute[3]

        # assert cl_compute_ee.shape[0] == config_dict[f"{spec}_n_bin"]

        return cl_compute_ee, cl_compute_bb