boxdm

boxdm

Author: Corentin Ravoux

Description : routines to extract a DM density box from a Saclay mocks output. The format of the output box corresponds to the one of a given Tomographic map. Tested on cori (NERSC)

BoxExtractor

BoxExtractor(pwd, box_dir, box_shape, size_cell, box_bound, master_file, interpolation_method='LINEAR')

Extract the underlying SaclayMocks dark-matter field onto a map grid.

Given the geometry of a tomographic map (or a set of lines of sight / catalog positions), samples the SaclayMocks Gaussian box(es), applies the growth factor and converts to a matter over-density, so a mock can be compared voxel-by-voxel with a reconstruction.

Store the box description and open a report log.

Parameters:

Name	Type	Description	Default
pwd	str	Output directory (and log location).	required
box_dir	str | list[str]	SaclayMocks box directory, or a list of directories for a multi-box footprint.	required
box_shape	tuple[int] | list	Box pixel shape(s).	required
size_cell	float	Box cell size (Mpc.h^-1).	required
box_bound	sequence	(ramin, ramax, decmin, decmax) footprint(s).	required
master_file	str	SaclayMocks master FITS file (growth table).	required
interpolation_method	str	"LINEAR" or "NEAREST".	'LINEAR'

Source code in lelantos/boxdm.py

def __init__(
    self,
    pwd,
    box_dir,
    box_shape,
    size_cell,
    box_bound,
    master_file,
    interpolation_method="LINEAR",
):
    """Store the box description and open a report log.

    Args:
        pwd (str): Output directory (and log location).
        box_dir (str | list[str]): SaclayMocks box directory, or a list of
            directories for a multi-box footprint.
        box_shape (tuple[int] | list): Box pixel shape(s).
        size_cell (float): Box cell size (Mpc.h^-1).
        box_bound (sequence): ``(ramin, ramax, decmin, decmax)`` footprint(s).
        master_file (str): SaclayMocks master FITS file (growth table).
        interpolation_method (str, optional): ``"LINEAR"`` or ``"NEAREST"``.
    """
    self.pwd = pwd
    self.size_cell = size_cell
    self.box_bound = box_bound
    self.box_shape = box_shape
    self.box_dir = box_dir
    self.master_file = master_file
    self.interpolation_method = interpolation_method
    self.log = utils.create_report_log(name=os.path.join(self.pwd, "Python_Report"))

create_box_array

create_box_array(X_tomo_min, X_tomo_max, Y_tomo_min, Y_tomo_max, Z_tomo_min, Z_tomo_max, shape_map)

Build the per-axis cartesian coordinate arrays of the output grid.

Parameters:

Name	Type	Description	Default
shape_map	tuple[int]	Output pixel shape.	required

Returns:

Name	Type	Description
tuple		(X_array, Y_array, Z_array) grid axes.

Source code in lelantos/boxdm.py

def create_box_array(
    self,
    X_tomo_min,
    X_tomo_max,
    Y_tomo_min,
    Y_tomo_max,
    Z_tomo_min,
    Z_tomo_max,
    shape_map,
):
    """Build the per-axis cartesian coordinate arrays of the output grid.

    Args:
        X_tomo_min, X_tomo_max, Y_tomo_min, Y_tomo_max, Z_tomo_min,
            Z_tomo_max (float): Cartesian bounds (Mpc.h^-1).
        shape_map (tuple[int]): Output pixel shape.

    Returns:
        tuple: ``(X_array, Y_array, Z_array)`` grid axes.
    """
    X_tomo_array = np.linspace(X_tomo_min, X_tomo_max, shape_map[0])
    Y_tomo_array = np.linspace(Y_tomo_min, Y_tomo_max, shape_map[1])
    Z_tomo_array = np.linspace(Z_tomo_min, Z_tomo_max, shape_map[2])
    return (X_tomo_array, Y_tomo_array, Z_tomo_array)

construct_DM_map

construct_DM_map(ra_array, dec_array, z_array, R_of_z, ra0_box, dec0_box, shape_map_output, Rmin, h, get_prop=None)

Sample the DM (and optional extra) fields from a single box.

Converts the output (RA, Dec, z) grid to SaclayMocks box coordinates and interpolates the density (and any requested extra fields) onto it.

Parameters:

Name	Type	Description	Default
ra_array, dec_array, z_array	array - like	Output grid sky coords.	required
R_of_z	callable	Redshift -> comoving distance.	required
ra0_box, dec0_box	float	Box footprint centre (degrees).	required
shape_map_output	tuple[int]	Output pixel shape (unused directly).	required
Rmin	float	Radial box lower bound (Mpc.h^-1).	required
h	float	Reduced Hubble constant.	required
get_prop	list[str]	Extra box fields to also sample.	None

Returns:

Name	Type	Description
tuple		(DM_map, Props_map) — the density map and stacked extra
		fields (or None).

Source code in lelantos/boxdm.py

def construct_DM_map(
    self,
    ra_array,
    dec_array,
    z_array,
    R_of_z,
    ra0_box,
    dec0_box,
    shape_map_output,
    Rmin,
    h,
    get_prop=None,
):
    """Sample the DM (and optional extra) fields from a single box.

    Converts the output (RA, Dec, z) grid to SaclayMocks box coordinates
    and interpolates the density (and any requested extra fields) onto it.

    Args:
        ra_array, dec_array, z_array (array-like): Output grid sky coords.
        R_of_z (callable): Redshift -> comoving distance.
        ra0_box, dec0_box (float): Box footprint centre (degrees).
        shape_map_output (tuple[int]): Output pixel shape (unused directly).
        Rmin (float): Radial box lower bound (Mpc.h^-1).
        h (float): Reduced Hubble constant.
        get_prop (list[str], optional): Extra box fields to also sample.

    Returns:
        tuple: ``(DM_map, Props_map)`` — the density map and stacked extra
        fields (or None).
    """
    self.log.add("Creation of (RA,DEC,R) data matrix")
    coords_ra_dec = np.moveaxis(
        np.array(
            np.meshgrid(ra_array, dec_array, h * R_of_z(z_array), indexing="ij")
        ),
        0,
        -1,
    )
    self.log.add("Conversion to (X,Y,Z) coordinates in the Saclay box")
    coords_box_saclay = np.zeros(coords_ra_dec.shape)
    (
        coords_box_saclay[:, :, :, 0],
        coords_box_saclay[:, :, :, 1],
        coords_box_saclay[:, :, :, 2],
    ) = utils.saclay_mock_sky_to_cartesian(
        coords_ra_dec[:, :, :, 0],
        coords_ra_dec[:, :, :, 1],
        coords_ra_dec[:, :, :, 2],
        ra0_box,
        dec0_box,
    )
    del coords_ra_dec
    self.log.add("Searching for the nearest pixels of the Saclay box")
    coords_pixels_box_saclay = np.zeros(coords_box_saclay.shape)
    (
        coords_pixels_box_saclay[:, :, :, 0],
        coords_pixels_box_saclay[:, :, :, 1],
        coords_pixels_box_saclay[:, :, :, 2],
    ) = utils.saclay_mock_coord_dm_map(
        coords_box_saclay[:, :, :, 0],
        coords_box_saclay[:, :, :, 1],
        coords_box_saclay[:, :, :, 2],
        Rmin,
        self.size_cell,
        self.box_shape,
        self.interpolation_method,
    )
    del coords_box_saclay
    self.log.add("Loading of the Saclay map")
    DM_mocks_map = utils.saclay_mock_get_box(self.box_dir, self.box_shape)
    self.log.add("Creation of the DM map")
    DM_map = utils.interpolate_map(
        self.interpolation_method, DM_mocks_map, coords_pixels_box_saclay
    )
    del DM_mocks_map
    if get_prop is not None:
        Props_map = []
        for i in range(len(get_prop)):
            self.log.add("Loading of the Saclay {} map".format(get_prop[i]))
            prop_mocks_map = utils.saclay_mock_get_box(
                self.box_dir, self.box_shape, name_box=get_prop[i]
            )
            self.log.add("Creation of the {} map".format(get_prop[i]))
            prop_map = utils.interpolate_map(
                self.interpolation_method, prop_mocks_map, coords_pixels_box_saclay
            )
            Props_map.append(prop_map)
            del prop_mocks_map
        Props_map = np.array(Props_map)
    else:
        Props_map = None
    del coords_pixels_box_saclay
    self.log.add("Multiplying by growth factor at redshift of the LOS")
    return (DM_map, Props_map)

fill_DM_map

fill_DM_map(ra_array, dec_array, z_array, R_of_z, ra0_box, dec0_box, shape_map_output, Rmin, h, DM_map, i_box, get_prop=None, Props_map=None)

Fill the still-empty voxels of the DM map from one box of a set.

Only voxels not yet filled and inside box i_box's footprint are sampled, so several boxes can be stitched into one output map.

Parameters:

Name	Type	Description	Default
ra_array, dec_array, z_array	array - like	Output grid sky coords.	required
R_of_z	callable	Redshift -> comoving distance.	required
ra0_box, dec0_box	float	Box footprint centre (degrees).	required
shape_map_output	tuple[int]	Output pixel shape.	required
Rmin	float	Radial box lower bound (Mpc.h^-1).	required
h	float	Reduced Hubble constant.	required
DM_map	ndarray	Density map being filled (updated in place).	required
i_box	int	Index of the box in the multi-box set.	required
get_prop	list[str]	Extra box fields to also sample.	None
Props_map	ndarray	Extra-field maps being filled.	None

Returns:

Name	Type	Description
tuple		(DM_map, Props_map).

Source code in lelantos/boxdm.py

def fill_DM_map(
    self,
    ra_array,
    dec_array,
    z_array,
    R_of_z,
    ra0_box,
    dec0_box,
    shape_map_output,
    Rmin,
    h,
    DM_map,
    i_box,
    get_prop=None,
    Props_map=None,
):
    """Fill the still-empty voxels of the DM map from one box of a set.

    Only voxels not yet filled and inside box ``i_box``'s footprint are
    sampled, so several boxes can be stitched into one output map.

    Args:
        ra_array, dec_array, z_array (array-like): Output grid sky coords.
        R_of_z (callable): Redshift -> comoving distance.
        ra0_box, dec0_box (float): Box footprint centre (degrees).
        shape_map_output (tuple[int]): Output pixel shape.
        Rmin (float): Radial box lower bound (Mpc.h^-1).
        h (float): Reduced Hubble constant.
        DM_map (numpy.ndarray): Density map being filled (updated in place).
        i_box (int): Index of the box in the multi-box set.
        get_prop (list[str], optional): Extra box fields to also sample.
        Props_map (numpy.ndarray, optional): Extra-field maps being filled.

    Returns:
        tuple: ``(DM_map, Props_map)``.
    """
    self.log.add("Creation of (RA,DEC,R) data matrix")
    coords_ra_dec = np.moveaxis(
        np.array(
            np.meshgrid(ra_array, dec_array, h * R_of_z(z_array), indexing="ij")
        ),
        0,
        -1,
    )
    mask1 = DM_map[:, :, :] == None
    mask2 = coords_ra_dec[:, :, :, 0] >= self.box_bound[i_box][0]
    mask2 &= coords_ra_dec[:, :, :, 0] < self.box_bound[i_box][1]
    mask2 &= coords_ra_dec[:, :, :, 1] >= self.box_bound[i_box][2]
    mask2 &= coords_ra_dec[:, :, :, 1] < self.box_bound[i_box][3]
    if (len(mask1[mask1 == True]) == 0) | (len(mask2[mask2 == True]) == 0):
        self.log.add("No need of the Saclay box {}".format(i_box))
        return (DM_map, Props_map)
    mask = mask1 & mask2
    del mask1, mask2
    self.log.add(
        "Conversion to (X,Y,Z) coordinates in the Saclay box {}".format(i_box)
    )
    coords_box_saclay = np.zeros(coords_ra_dec.shape)
    (
        coords_box_saclay[:, :, :, 0],
        coords_box_saclay[:, :, :, 1],
        coords_box_saclay[:, :, :, 2],
    ) = utils.saclay_mock_sky_to_cartesian(
        coords_ra_dec[:, :, :, 0],
        coords_ra_dec[:, :, :, 1],
        coords_ra_dec[:, :, :, 2],
        ra0_box,
        dec0_box,
    )
    del coords_ra_dec
    self.log.add(
        "Searching for the nearest pixels of the Saclay box {}".format(i_box)
    )
    coords_pixels_box_saclay = np.zeros(coords_box_saclay.shape)
    (
        coords_pixels_box_saclay[:, :, :, 0],
        coords_pixels_box_saclay[:, :, :, 1],
        coords_pixels_box_saclay[:, :, :, 2],
    ) = utils.saclay_mock_coord_dm_map(
        coords_box_saclay[:, :, :, 0],
        coords_box_saclay[:, :, :, 1],
        coords_box_saclay[:, :, :, 2],
        Rmin,
        self.size_cell,
        self.box_shape[i_box],
        self.interpolation_method,
    )
    del coords_box_saclay
    self.log.add("Loading of the Saclay map {}".format(i_box))
    DM_mocks_map = utils.saclay_mock_get_box(
        self.box_dir[i_box], self.box_shape[i_box]
    )
    self.log.add("Creation of the DM map with box {}".format(i_box))
    DM_map[mask] = utils.interpolate_and_fill_map(
        self.interpolation_method, DM_mocks_map, coords_pixels_box_saclay[mask]
    )
    del DM_mocks_map
    if get_prop is not None:
        for i in range(len(get_prop)):
            self.log.add(
                "Loading of the Saclay {} map {}".format(get_prop[i], i_box)
            )
            prop_mocks_map = utils.saclay_mock_get_box(
                self.box_dir[i_box], self.box_shape[i_box], name_box=get_prop[i]
            )
            self.log.add(
                "Creation of the {} map with box {}".format(get_prop[i], i_box)
            )
            Props_map[i][mask] = utils.interpolate_and_fill_map(
                self.interpolation_method,
                prop_mocks_map,
                coords_pixels_box_saclay[mask],
            )
            del prop_mocks_map
    del coords_pixels_box_saclay, mask, get_prop
    return (DM_map, Props_map)

extract_box

extract_box(ra_array, dec_array, z_array, shape_map_output, rsd_box)

Extract the DM map from a single box (optionally with RSD field).

Parameters:

Name	Type	Description	Default
ra_array, dec_array, z_array	array - like	Output grid sky coords.	required
shape_map_output	tuple[int]	Output pixel shape.	required
rsd_box	bool	Also sample the eta_zz field for RSD.	required

Returns:

Name	Type	Description
tuple		(DM_map, Props_map, get_prop).

Source code in lelantos/boxdm.py

def extract_box(self, ra_array, dec_array, z_array, shape_map_output, rsd_box):
    """Extract the DM map from a single box (optionally with RSD field).

    Args:
        ra_array, dec_array, z_array (array-like): Output grid sky coords.
        shape_map_output (tuple[int]): Output pixel shape.
        rsd_box (bool): Also sample the ``eta_zz`` field for RSD.

    Returns:
        tuple: ``(DM_map, Props_map, get_prop)``.
    """
    (
        R0,
        z0,
        R_of_z,
        z_of_R,
        Rmin,
        Rmax,
        h,
    ) = utils.saclay_mock_box_cosmo_parameters(self.box_shape, self.size_cell)
    ra0_box, dec0_box = utils.saclay_mock_center_of_the_box(self.box_bound)
    if rsd_box:
        get_prop = ["eta_zz"]
    else:
        get_prop = None
    DM_map, Props_map = self.construct_DM_map(
        ra_array,
        dec_array,
        z_array,
        R_of_z,
        ra0_box,
        dec0_box,
        shape_map_output,
        Rmin,
        h,
        get_prop=get_prop,
    )
    return (DM_map, Props_map, get_prop)

extract_box_multiple

extract_box_multiple(ra_array, dec_array, z_array, shape_map_output, rsd_box)

Extract the DM map by stitching several boxes over the footprint.

Parameters:

Name	Type	Description	Default
ra_array, dec_array, z_array	array - like	Output grid sky coords.	required
shape_map_output	tuple[int]	Output pixel shape.	required
rsd_box	bool	Also sample the eta_zz and vz fields for RSD.	required

Returns:

Name	Type	Description
tuple		(DM_map, Props_map, get_prop).

Source code in lelantos/boxdm.py

def extract_box_multiple(
    self, ra_array, dec_array, z_array, shape_map_output, rsd_box
):
    """Extract the DM map by stitching several boxes over the footprint.

    Args:
        ra_array, dec_array, z_array (array-like): Output grid sky coords.
        shape_map_output (tuple[int]): Output pixel shape.
        rsd_box (bool): Also sample the ``eta_zz`` and ``vz`` fields for RSD.

    Returns:
        tuple: ``(DM_map, Props_map, get_prop)``.
    """
    DM_map = np.full(tuple(shape_map_output), None)
    if rsd_box:
        get_prop = ["eta_zz", "vz"]
        Props_map = np.array(
            [np.full(tuple(shape_map_output), None) for i in range(len(get_prop))]
        )
    else:
        Props_map, get_prop = None, None
    for i_box in range(len(self.box_dir)):
        (
            R0,
            z0,
            R_of_z,
            z_of_R,
            Rmin,
            Rmax,
            h,
        ) = utils.saclay_mock_box_cosmo_parameters(
            self.box_shape[i_box], self.size_cell
        )
        ra0_box, dec0_box = utils.saclay_mock_center_of_the_box(
            self.box_bound[i_box]
        )
        DM_map, Props_map = self.fill_DM_map(
            ra_array,
            dec_array,
            z_array,
            R_of_z,
            ra0_box,
            dec0_box,
            shape_map_output,
            Rmin,
            h,
            DM_map,
            i_box,
            get_prop=get_prop,
            Props_map=Props_map,
        )
    if len(DM_map[DM_map == None] != 0):
        self.log.add("WARNING : Not enough boxes to fill the Dark Matter map")
    self.log.add("Multiplying by growth factor at redshift of the LOS")
    return (DM_map, Props_map, get_prop)

multiply_by_growth

multiply_by_growth(DM_map, z_array)

Scale each radial slice of the map by the growth factor at its z.

Parameters:

Name	Type	Description	Default
DM_map	ndarray	Gaussian density map.	required
z_array	array - like	Redshift of each radial slice.	required

Returns:

Type	Description
	numpy.ndarray: The growth-scaled map.

Source code in lelantos/boxdm.py

def multiply_by_growth(self, DM_map, z_array):
    """Scale each radial slice of the map by the growth factor at its z.

    Args:
        DM_map (numpy.ndarray): Gaussian density map.
        z_array (array-like): Redshift of each radial slice.

    Returns:
        numpy.ndarray: The growth-scaled map.
    """
    Z = fitsio.FITS(self.master_file)[2]["Z"][:]
    G = fitsio.FITS(self.master_file)[2]["G"][:]
    G_of_Z = interpolate.interp1d(Z, G)
    for i in range(len(z_array)):
        DM_map[:, :, i] = DM_map[:, :, i] * G_of_Z(z_array[i])
    return DM_map

multiply_los_by_growth

multiply_los_by_growth(DM_los, z_array)

Scale line-of-sight values by the growth factor at each redshift.

Parameters:

Name	Type	Description	Default
DM_los	ndarray	Gaussian density along a line of sight.	required
z_array	array - like	Redshift of each pixel.	required

Returns:

Type	Description
	numpy.ndarray: The growth-scaled line of sight.

Source code in lelantos/boxdm.py

def multiply_los_by_growth(self, DM_los, z_array):
    """Scale line-of-sight values by the growth factor at each redshift.

    Args:
        DM_los (numpy.ndarray): Gaussian density along a line of sight.
        z_array (array-like): Redshift of each pixel.

    Returns:
        numpy.ndarray: The growth-scaled line of sight.
    """
    Z = fitsio.FITS(self.master_file)[2]["Z"][:]
    G = fitsio.FITS(self.master_file)[2]["G"][:]
    G_of_Z = interpolate.interp1d(Z, G)
    DM_los = DM_los * G_of_Z(z_array)
    return DM_los

convert_to_matter_field

convert_to_matter_field(gaussian_array)

Log-normal transform a Gaussian field into a matter over-density.

Parameters:

Name	Type	Description	Default
gaussian_array	ndarray	Gaussian field values.	required

Returns:

Type	Description
	numpy.ndarray: The matter over-density exp(g - sigma^2/2) - 1.

Source code in lelantos/boxdm.py

def convert_to_matter_field(self, gaussian_array):
    """Log-normal transform a Gaussian field into a matter over-density.

    Args:
        gaussian_array (numpy.ndarray): Gaussian field values.

    Returns:
        numpy.ndarray: The matter over-density ``exp(g - sigma^2/2) - 1``.
    """
    sigma_l = np.std(gaussian_array)
    density_matter = np.exp((gaussian_array - (sigma_l**2 / 2)).astype(float)) - 1
    return density_matter

compute_rsd

compute_rsd(DM_map, Props_map, z_array, Cosmo)

Apply a redshift-space distortion correction to the DM map.

Parameters:

Name	Type	Description	Default
DM_map	ndarray	Density map.	required
Props_map	ndarray	Extra fields (uses the velocity term).	required
z_array	array - like	Redshift of each radial slice.	required
Cosmo		Object exposing hubble(z).	required

Returns:

Type	Description
	numpy.ndarray: The RSD-corrected map.

Source code in lelantos/boxdm.py

def compute_rsd(self, DM_map, Props_map, z_array, Cosmo):
    """Apply a redshift-space distortion correction to the DM map.

    Args:
        DM_map (numpy.ndarray): Density map.
        Props_map (numpy.ndarray): Extra fields (uses the velocity term).
        z_array (array-like): Redshift of each radial slice.
        Cosmo: Object exposing ``hubble(z)``.

    Returns:
        numpy.ndarray: The RSD-corrected map.
    """
    for i in range(len(z_array)):
        DM_map[:, :, i] = DM_map[:, :, i] - (
            ((1 + z_array[i]) * Props_map[0][:, :, i]) / Cosmo.hubble(z_array[i])
        )
    return DM_map

construct_DM_LOS

construct_DM_LOS(ra_array, dec_array, z_array, R_of_z, ra0_box, dec0_box, Rmin, h)

Sample the DM field along lines of sight from a single box.

Parameters:

Name	Type	Description	Default
ra_array, dec_array, z_array	array - like	Line-of-sight sky coords.	required
R_of_z	callable	Redshift -> comoving distance.	required
ra0_box, dec0_box	float	Box footprint centre (degrees).	required
Rmin	float	Radial box lower bound (Mpc.h^-1).	required
h	float	Reduced Hubble constant.	required

Returns:

Type	Description
	numpy.ndarray: DM field value at each line-of-sight pixel.

Source code in lelantos/boxdm.py

def construct_DM_LOS(
    self, ra_array, dec_array, z_array, R_of_z, ra0_box, dec0_box, Rmin, h
):
    """Sample the DM field along lines of sight from a single box.

    Args:
        ra_array, dec_array, z_array (array-like): Line-of-sight sky coords.
        R_of_z (callable): Redshift -> comoving distance.
        ra0_box, dec0_box (float): Box footprint centre (degrees).
        Rmin (float): Radial box lower bound (Mpc.h^-1).
        h (float): Reduced Hubble constant.

    Returns:
        numpy.ndarray: DM field value at each line-of-sight pixel.
    """
    self.log.add("Creation of (RA,DEC,R) data matrix")
    R_array = h * R_of_z(z_array)
    self.log.add("Conversion to (X,Y,Z) coordinates in the Saclay box")
    X, Y, Z = (
        np.zeros(len(ra_array)),
        np.zeros(len(dec_array)),
        np.zeros(len(R_array)),
    )
    X, Y, Z = utils.saclay_mock_sky_to_cartesian(
        ra_array, dec_array, R_array, ra0_box, dec0_box
    )
    self.log.add("Searching for the nearest pixels of the Saclay box")
    i, j, k = (
        np.zeros(len(R_array)).astype(int),
        np.zeros(len(R_array)).astype(int),
        np.zeros(len(R_array)).astype(int),
    )
    i, j, k = utils.saclay_mock_coord_dm_map(
        X, Y, Z, Rmin, self.size_cell, self.box_shape, self.interpolation_method
    )
    del X, Y, Z
    self.log.add("Loading of the Saclay map")
    DM_mocks_map = utils.saclay_mock_get_box(self.box_dir, self.box_shape)
    self.log.add("Creation of the DM los")
    DM_LOS = np.zeros(len(R_array))
    DM_LOS[:] = DM_mocks_map[i[:], j[:], k[:]]
    del DM_mocks_map, i, j, k
    return DM_LOS

fill_DM_LOS

fill_DM_LOS(ra_array, dec_array, z_array, R_of_z, ra0_box, dec0_box, Rmin, h, DM_LOS, i_box)

Fill still-empty line-of-sight pixels from one box of a set.

Parameters:

Name	Type	Description	Default
ra_array, dec_array, z_array	array - like	Line-of-sight sky coords.	required
R_of_z	callable	Redshift -> comoving distance.	required
ra0_box, dec0_box	float	Box footprint centre (degrees).	required
Rmin	float	Radial box lower bound (Mpc.h^-1).	required
h	float	Reduced Hubble constant.	required
DM_LOS	ndarray	Line-of-sight values (updated in place).	required
i_box	int	Index of the box in the multi-box set.	required

Returns:

Type	Description
	numpy.ndarray: The updated line-of-sight array.

Source code in lelantos/boxdm.py

def fill_DM_LOS(
    self,
    ra_array,
    dec_array,
    z_array,
    R_of_z,
    ra0_box,
    dec0_box,
    Rmin,
    h,
    DM_LOS,
    i_box,
):
    """Fill still-empty line-of-sight pixels from one box of a set.

    Args:
        ra_array, dec_array, z_array (array-like): Line-of-sight sky coords.
        R_of_z (callable): Redshift -> comoving distance.
        ra0_box, dec0_box (float): Box footprint centre (degrees).
        Rmin (float): Radial box lower bound (Mpc.h^-1).
        h (float): Reduced Hubble constant.
        DM_LOS (numpy.ndarray): Line-of-sight values (updated in place).
        i_box (int): Index of the box in the multi-box set.

    Returns:
        numpy.ndarray: The updated line-of-sight array.
    """
    self.log.add("Creation of (RA,DEC,R) data matrix")
    R_array = h * R_of_z(z_array)
    mask1 = DM_LOS[:] == None
    mask2 = ra_array[:] >= self.box_bound[i_box][0]
    mask2 &= ra_array[:] < self.box_bound[i_box][1]
    mask2 &= dec_array[:] >= self.box_bound[i_box][2]
    mask2 &= dec_array[:] < self.box_bound[i_box][3]
    if (len(mask1[mask1 == True]) == 0) | (len(mask2[mask2 == True]) == 0):
        self.log.add("No need of the Saclay box {}".format(i_box))
        return DM_LOS
    mask = mask1 & mask2
    del mask1, mask2
    self.log.add(
        "Conversion to (X,Y,Z) coordinates in the Saclay box {}".format(i_box)
    )
    X, Y, Z = (
        np.zeros(len(ra_array)),
        np.zeros(len(dec_array)),
        np.zeros(len(R_array)),
    )
    X, Y, Z = utils.saclay_mock_sky_to_cartesian(
        ra_array, dec_array, R_array, ra0_box, dec0_box
    )
    self.log.add(
        "Searching for the nearest pixels of the Saclay box {}".format(i_box)
    )
    i, j, k = (
        np.zeros(len(R_array)).astype(int),
        np.zeros(len(R_array)).astype(int),
        np.zeros(len(R_array)).astype(int),
    )
    i, j, k = utils.saclay_mock_coord_dm_map(
        X,
        Y,
        Z,
        Rmin,
        self.size_cell,
        self.box_shape[i_box],
        self.interpolation_method,
    )
    del X, Y, Z
    self.log.add("Loading of the Saclay map {}".format(i_box))
    DM_mocks_map = utils.saclay_mock_get_box(
        self.box_dir[i_box], self.box_shape[i_box]
    )
    self.log.add("Creation of the DM los with box {}".format(i_box))
    DM_LOS[:][mask] = DM_mocks_map[i[:][mask], j[:][mask], k[:][mask]]
    del DM_mocks_map, i, j, k
    return DM_LOS

extract_los

extract_los(ra_array, dec_array, z_array)

Extract DM lines of sight from a single box.

Parameters:

Name	Type	Description	Default
ra_array, dec_array, z_array	array - like	Line-of-sight sky coords.	required

Returns:

Type	Description
	numpy.ndarray: DM field along the lines of sight.

Source code in lelantos/boxdm.py

def extract_los(self, ra_array, dec_array, z_array):
    """Extract DM lines of sight from a single box.

    Args:
        ra_array, dec_array, z_array (array-like): Line-of-sight sky coords.

    Returns:
        numpy.ndarray: DM field along the lines of sight.
    """
    (
        R0,
        z0,
        R_of_z,
        z_of_R,
        Rmin,
        Rmax,
        h,
    ) = utils.saclay_mock_box_cosmo_parameters(self.box_shape, self.size_cell)
    ra0_box, dec0_box = utils.saclay_mock_center_of_the_box(self.box_bound)
    los = self.construct_DM_LOS(
        ra_array, dec_array, z_array, R_of_z, ra0_box, dec0_box, Rmin, h
    )
    return los

extract_los_multiple

extract_los_multiple(ra_array, dec_array, z_array)

Extract DM lines of sight by stitching several boxes.

Parameters:

Name	Type	Description	Default
ra_array, dec_array, z_array	array - like	Line-of-sight sky coords.	required

Returns:

Type	Description
	numpy.ndarray: DM field along the lines of sight.

Source code in lelantos/boxdm.py

def extract_los_multiple(self, ra_array, dec_array, z_array):
    """Extract DM lines of sight by stitching several boxes.

    Args:
        ra_array, dec_array, z_array (array-like): Line-of-sight sky coords.

    Returns:
        numpy.ndarray: DM field along the lines of sight.
    """
    los = np.full(z_array.shape, None)
    for i_box in range(len(self.box_dir)):
        (
            R0,
            z0,
            R_of_z,
            z_of_R,
            Rmin,
            Rmax,
            h,
        ) = utils.saclay_mock_box_cosmo_parameters(
            self.box_shape[i_box], self.size_cell
        )
        ra0_box, dec0_box = utils.saclay_mock_center_of_the_box(
            self.box_bound[i_box]
        )
        los = self.fill_DM_LOS(
            ra_array,
            dec_array,
            z_array,
            R_of_z,
            ra0_box,
            dec0_box,
            Rmin,
            h,
            los,
            i_box,
        )
    if len(los[los == None] != 0):
        self.log.add("Warning : Not enough boxes to fill the Dark Matter map")
    return los

extract_delta_multiple

extract_delta_multiple(ra, dec, z)

Extract DM values at catalog (QSO) positions from several boxes.

Parameters:

Name	Type	Description	Default
ra, dec, z	array - like	Catalog sky coordinates.	required

Returns:

Type	Description
	numpy.ndarray: DM field at each catalog position.

Source code in lelantos/boxdm.py

def extract_delta_multiple(self, ra, dec, z):
    """Extract DM values at catalog (QSO) positions from several boxes.

    Args:
        ra, dec, z (array-like): Catalog sky coordinates.

    Returns:
        numpy.ndarray: DM field at each catalog position.
    """
    delta_quasars = np.full(z.shape, None)
    for i_box in range(len(self.box_dir)):
        (
            R0,
            z0,
            R_of_z,
            z_of_R,
            Rmin,
            Rmax,
            h,
        ) = utils.saclay_mock_box_cosmo_parameters(
            self.box_shape[i_box], self.size_cell
        )
        ra0_box, dec0_box = utils.saclay_mock_center_of_the_box(
            self.box_bound[i_box]
        )
        delta_quasars = self.fill_DM_LOS(
            ra, dec, z, R_of_z, ra0_box, dec0_box, Rmin, h, delta_quasars, i_box
        )
    if len(delta_quasars[delta_quasars == None] != 0):
        self.log.add("Warning : Not enough boxes to fill the Dark Matter map")
    return delta_quasars

extract_delta

extract_delta(ra, dec, z)

Extract DM values at catalog (QSO) positions from a single box.

Parameters:

Name	Type	Description	Default
ra, dec, z	array - like	Catalog sky coordinates.	required

Returns:

Type	Description
	numpy.ndarray: DM field at each catalog position.

Source code in lelantos/boxdm.py

def extract_delta(self, ra, dec, z):
    """Extract DM values at catalog (QSO) positions from a single box.

    Args:
        ra, dec, z (array-like): Catalog sky coordinates.

    Returns:
        numpy.ndarray: DM field at each catalog position.
    """
    (
        R0,
        z0,
        R_of_z,
        z_of_R,
        Rmin,
        Rmax,
        h,
    ) = utils.saclay_mock_box_cosmo_parameters(self.box_shape, self.size_cell)
    ra0_box, dec0_box = utils.saclay_mock_center_of_the_box(self.box_bound)
    delta_quasars = self.construct_DM_LOS(
        ra, dec, z, R_of_z, ra0_box, dec0_box, Rmin, h
    )
    return delta_quasars

create_box

create_box(map_property_file, name, rsd_box=False, growth_multiplication=True, matter_field=True, gaussian_smoothing=None, shape_map_output=None)

Build the DM map matching a tomographic map and write it to disk.

Reads the map property file to get the geometry/cosmology, samples the box(es), applies the growth factor, log-normal transform and optional Gaussian smoothing, and writes the resulting map (plus RSD fields).

Parameters:

Name	Type	Description	Default
map_property_file	str	Property file of the reference map.	required
name	str	Output map file.	required
rsd_box	bool	Also extract RSD fields.	False
growth_multiplication	bool	Apply the growth factor.	True
matter_field	bool	Apply the log-normal transform.	True
gaussian_smoothing	float	Smoothing scale (Mpc.h^-1).	None
shape_map_output	tuple[int]	Override the output shape.	None

Source code in lelantos/boxdm.py

def create_box(
    self,
    map_property_file,
    name,
    rsd_box=False,
    growth_multiplication=True,
    matter_field=True,
    gaussian_smoothing=None,
    shape_map_output=None,
):
    """Build the DM map matching a tomographic map and write it to disk.

    Reads the map property file to get the geometry/cosmology, samples the
    box(es), applies the growth factor, log-normal transform and optional
    Gaussian smoothing, and writes the resulting map (plus RSD fields).

    Args:
        map_property_file (str): Property file of the reference map.
        name (str): Output map file.
        rsd_box (bool, optional): Also extract RSD fields.
        growth_multiplication (bool, optional): Apply the growth factor.
        matter_field (bool, optional): Apply the log-normal transform.
        gaussian_smoothing (float, optional): Smoothing scale (Mpc.h^-1).
        shape_map_output (tuple[int], optional): Override the output shape.
    """
    property_file = tomographic_objects.MapPixelProperty(name=map_property_file)
    property_file.read()
    if shape_map_output is None:
        shape_map_output = property_file.shape
    X_tomo_min, Y_tomo_min, Z_tomo_min = property_file.boundary_cartesian_coord[0]
    X_tomo_max, Y_tomo_max, Z_tomo_max = property_file.boundary_cartesian_coord[1]
    X_tomo_array, Y_tomo_array, Z_tomo_array = self.create_box_array(
        X_tomo_min,
        X_tomo_max,
        Y_tomo_min,
        Y_tomo_max,
        Z_tomo_min,
        Z_tomo_max,
        shape_map_output,
    )
    suplementary_parameters = utils.return_suplementary_parameters(
        property_file.coordinate_transform, property=property_file
    )
    (rcomov, distang, inv_rcomov, inv_distang) = utils.get_cosmo_function(
        property_file.Omega_m
    )
    ra_array, dec_array, z_array = utils.convert_cartesian_to_sky(
        X_tomo_array,
        Y_tomo_array,
        Z_tomo_array,
        property_file.coordinate_transform,
        inv_rcomov=inv_rcomov,
        inv_distang=inv_distang,
        distang=distang,
        suplementary_parameters=suplementary_parameters,
    )
    ra_array = np.degrees(ra_array)
    dec_array = np.degrees(dec_array)
    del (
        X_tomo_array,
        Y_tomo_array,
        Z_tomo_array,
        suplementary_parameters,
        rcomov,
        distang,
        inv_rcomov,
        inv_distang,
    )
    if type(self.box_dir) == list:
        (dm_map, prop_maps, prop) = self.extract_box_multiple(
            ra_array, dec_array, z_array, shape_map_output, rsd_box
        )
    else:
        (dm_map, prop_maps, prop) = self.extract_box(
            ra_array, dec_array, z_array, shape_map_output, rsd_box
        )
    del ra_array, dec_array
    if growth_multiplication:
        dm_map = self.multiply_by_growth(dm_map, z_array)
    del z_array
    if matter_field:
        dm_map = self.convert_to_matter_field(dm_map)
    if gaussian_smoothing is not None:
        gaussian_smoothing_pix = gaussian_smoothing * utils.pixel_per_mpc(
            property_file.size, shape_map_output
        )
        dm_map = utils.gaussian_smoothing(dm_map, gaussian_smoothing_pix)
    dm_map_object = tomographic_objects.TomographicMap(
        map_array=np.array(dm_map), name=name
    )
    del dm_map
    dm_map_object.write()
    del dm_map_object
    if rsd_box:
        for i in range(len(prop)):
            prop_map_object = tomographic_objects.TomographicMap(
                map_array=prop_maps[i], name=f"{name}_{prop[i]}"
            )
            prop_map_object.write()
        del prop_map_object
    del prop_maps

create_LOS

create_LOS(property_file, pixel_name, name, growth_multiplication=True, matter_field=True)

Build DM values along the lines of sight of a pixel file.

Parameters:

Name	Type	Description	Default
property_file	str	Map property/pixel geometry file.	required
pixel_name	str	Pixel (line-of-sight) file.	required
name	str	Output pixel file.	required
growth_multiplication	bool	Apply the growth factor.	True
matter_field	bool	Apply the log-normal transform.	True

Source code in lelantos/boxdm.py

def create_LOS(
    self,
    property_file,
    pixel_name,
    name,
    growth_multiplication=True,
    matter_field=True,
):
    """Build DM values along the lines of sight of a pixel file.

    Args:
        property_file (str): Map property/pixel geometry file.
        pixel_name (str): Pixel (line-of-sight) file.
        name (str): Output pixel file.
        growth_multiplication (bool, optional): Apply the growth factor.
        matter_field (bool, optional): Apply the log-normal transform.
    """
    pixel = tomographic_objects.Pixel(name=pixel_name)
    pixel.read()
    property_file = tomographic_objects.MapPixelProperty(name=property_file)
    property_file.read()
    suplementary_parameters = utils.return_suplementary_parameters(
        property_file.coordinate_transform, property=property_file
    )
    (rcomov, distang, inv_rcomov, inv_distang) = utils.get_cosmo_function(
        property_file.Omega_m
    )
    X_array, Y_array, Z_array = (
        np.array(pixel.pixel_array[:, 0]),
        np.array(pixel.pixel_array[:, 1]),
        np.array(pixel.pixel_array[:, 2]),
    )
    ra_array, dec_array, z_array = utils.convert_cartesian_to_sky(
        X_array,
        Y_array,
        Z_array,
        property_file.coordinate_transform,
        inv_rcomov=inv_rcomov,
        inv_distang=inv_distang,
        distang=distang,
        suplementary_parameters=suplementary_parameters,
    )
    ra_array = np.degrees(ra_array)
    dec_array = np.degrees(dec_array)
    if type(self.box_dir) == list:
        los = self.extract_los(ra_array, dec_array, z_array)
    else:
        los = self.extract_los_multiple(ra_array, dec_array, z_array)
    self.log.add("Multiplying by growth factor at redshift of the LOS")
    if growth_multiplication:
        los = self.multiply_los_by_growth(los, z_array)
    if matter_field:
        los = self.convert_to_matter_field(los)
    self.log.add("Mean delta extracted : {}".format(np.mean(los)))
    pixel_out = tomographic_objects.Pixel(name=name, pixel_array=los)
    pixel_out.read()

create_catalog

create_catalog(cat_name, type_catalog, name, growth_multiplication=True, matter_field=True)

Extract DM over-density at the positions of a catalog.

Parameters:

Name	Type	Description	Default
cat_name	str	Catalog file.	required
type_catalog	str	Catalog type ("qso"/"void"/"galaxy").	required
name	str	Output name (unused for the return value).	required
growth_multiplication	bool	Apply the growth factor.	True
matter_field	bool	Apply the log-normal transform.	True

Returns:

Type	Description
	numpy.ndarray: DM over-density at each catalog position.

Source code in lelantos/boxdm.py

def create_catalog(
    self,
    cat_name,
    type_catalog,
    name,
    growth_multiplication=True,
    matter_field=True,
):
    """Extract DM over-density at the positions of a catalog.

    Args:
        cat_name (str): Catalog file.
        type_catalog (str): Catalog type (``"qso"``/``"void"``/``"galaxy"``).
        name (str): Output name (unused for the return value).
        growth_multiplication (bool, optional): Apply the growth factor.
        matter_field (bool, optional): Apply the log-normal transform.

    Returns:
        numpy.ndarray: DM over-density at each catalog position.
    """
    catalog = tomographic_objects.Catalog.init_catalog_from_fits(
        cat_name, type_catalog
    )
    ra, dec, z = catalog.coord[:, 0], catalog.coord[:, 1], catalog.coord[:, 2]
    ra[ra > 180] = ra[ra > 180] - 360
    if type(self.box_dir) == list:
        delta_quasars = self.create_delta_catalog_several_boxes(ra, dec, z)
    else:
        delta_quasars = self.create_delta_catalog(ra, dec, z)
    self.log.add("Multiplying by growth factor at redshift of the LOS")
    if growth_multiplication:
        delta_quasars = self.multiply_los_by_growth(delta_quasars, z)
    if matter_field:
        delta_quasars = self.convert_to_matter_field(delta_quasars)
    self.log.add("Mean delta extracted : {}".format(np.mean(delta_quasars)))
    return delta_quasars

BoxPlot

BoxPlot(pwd, map_name, box_name, property_file, limit=None)

Bases: object

Compare a reconstructed map with its underlying DM box via 2D plots.

Store the map / box files and plot limits.

Parameters:

Name	Type	Description	Default
pwd	str	Output directory.	required
map_name	str	Reconstructed map file.	required
box_name	str	DM box map file.	required
property_file	str	Shared map property file.	required
limit	sequence	(xmin, xmax, ymin, ymax) plot limits.	None

Source code in lelantos/boxdm.py

def __init__(self, pwd, map_name, box_name, property_file, limit=None):
    """Store the map / box files and plot limits.

    Args:
        pwd (str): Output directory.
        map_name (str): Reconstructed map file.
        box_name (str): DM box map file.
        property_file (str): Shared map property file.
        limit (sequence, optional): ``(xmin, xmax, ymin, ymax)`` plot limits.
    """
    self.pwd = pwd
    self.map_name = map_name
    self.box_name = box_name
    self.property_file = property_file
    self.limit = limit

contourplot staticmethod

contourplot(x, y, ncont=10, colors=None, pltclf=True, binsx=100, binsy=100, log=False)

Draw contour lines of the 2D histogram of (x, y).

Parameters:

Name	Type	Description	Default
x, y	array - like	Point coordinates.	required
ncont	int	Number of contour levels.	10
colors	optional	Contour colour(s).	None
pltclf	bool	Clear the figure first.	True
binsx, binsy	int	Histogram bin counts.	required
log	bool	Log-scale the histogram counts.	False

Source code in lelantos/boxdm.py

@staticmethod
def contourplot(
    x, y, ncont=10, colors=None, pltclf=True, binsx=100, binsy=100, log=False
):
    """Draw contour lines of the 2D histogram of ``(x, y)``.

    Args:
        x, y (array-like): Point coordinates.
        ncont (int, optional): Number of contour levels.
        colors (optional): Contour colour(s).
        pltclf (bool, optional): Clear the figure first.
        binsx, binsy (int, optional): Histogram bin counts.
        log (bool, optional): Log-scale the histogram counts.
    """
    H, xedges, yedges = np.histogram2d(x, y, bins=(binsx, binsy))
    H = np.rot90(H)
    H = np.flipud(H)
    xcenters = (xedges[:-1] + xedges[1:]) / 2.0
    ycenters = (yedges[:-1] + yedges[1:]) / 2.0
    if log is True:
        H[(H > 0)] = np.log(H[(H > 0)])
    if pltclf:
        plt.clf()
    plt.contour(xcenters, ycenters, H, ncont, colors=colors, linewidths=2)

densityplot staticmethod

densityplot(x, y, binsx=200, binsy=200, scaleperdeg2=False, maxdensity=0, mindensity=0, pltclf=True, logscale=False)

Draw a normalised 2D density image of (x, y).

Parameters:

Name	Type	Description	Default
x, y	array - like	Point coordinates.	required
binsx, binsy	int	Histogram bin counts.	required
scaleperdeg2	bool	Normalise to a density per deg^2.	False
maxdensity, mindensity	float	Clip the density range.	required
pltclf	bool	Clear the figure first.	True
logscale	bool	Use a logarithmic colour scale.	False

Source code in lelantos/boxdm.py

@staticmethod
def densityplot(
    x,
    y,
    binsx=200,
    binsy=200,
    scaleperdeg2=False,
    maxdensity=0,
    mindensity=0,
    pltclf=True,
    logscale=False,
):
    """Draw a normalised 2D density image of ``(x, y)``.

    Args:
        x, y (array-like): Point coordinates.
        binsx, binsy (int, optional): Histogram bin counts.
        scaleperdeg2 (bool, optional): Normalise to a density per deg^2.
        maxdensity, mindensity (float, optional): Clip the density range.
        pltclf (bool, optional): Clear the figure first.
        logscale (bool, optional): Use a logarithmic colour scale.
    """
    H, xedges, yedges = np.histogram2d(x, y, bins=(binsx, binsy))
    if scaleperdeg2 == True:
        delta_ra = abs(xedges[1] - xedges[0])
        delta_dec = abs(yedges[1] - yedges[0])
        cosdelta = np.cos(yedges[:-1] * np.pi / 180.0)
        for k in range(binsy):
            H[:, k] = H[:, k] / (delta_ra * delta_dec * cosdelta[k])
    H = np.rot90(H)
    H = np.flipud(H)
    if maxdensity > 0:
        wcut = np.where((H > maxdensity))
        H[wcut] = maxdensity
    if mindensity > 0:
        wcut = np.where((H < mindensity))
        H[wcut] = 0
    Hmasked = np.ma.masked_where(H == 0, H)
    if pltclf:
        plt.clf()
    if logscale is False:
        extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
        Hmaskednorm = Hmasked / np.max(Hmasked)
        plt.imshow(
            Hmaskednorm,
            interpolation="bilinear",
            cmap="jet",
            origin="lower",
            extent=extent,
            norm=Normalize(vmin=0.0, vmax=1.0),
        )
    else:
        plt.pcolormesh(xedges, yedges, Hmasked, norm=LogNorm(), cmap=plt.cm.jet)
    cbar = plt.colorbar()
    cbar.set_label("Normalized density")
    if scaleperdeg2 == True:
        plt.title(r"density / deg$^2$")

gaussian_plot staticmethod

gaussian_plot(X, Y, binsx=200, binsy=200, ncont=4)

Fit a 2D Gaussian to the (X, Y) histogram and overlay contours.

Parameters:

Name	Type	Description	Default
X, Y	array - like	Point coordinates.	required
binsx, binsy	int	Histogram bin counts.	required
ncont	int | sequence	Contour levels (scaled by the fit amplitude if not an int).	4

Returns:

Type	Description
	numpy.ndarray: The fitted Gaussian parameters.

Source code in lelantos/boxdm.py

@staticmethod
def gaussian_plot(X, Y, binsx=200, binsy=200, ncont=4):
    """Fit a 2D Gaussian to the ``(X, Y)`` histogram and overlay contours.

    Args:
        X, Y (array-like): Point coordinates.
        binsx, binsy (int, optional): Histogram bin counts.
        ncont (int | sequence, optional): Contour levels (scaled by the fit
            amplitude if not an int).

    Returns:
        numpy.ndarray: The fitted Gaussian parameters.
    """
    data, xedges, yedges = np.histogram2d(X, Y, bins=(binsx, binsy))
    xcenters = (xedges[:-1] + xedges[1:]) / 2.0
    ycenters = (yedges[:-1] + yedges[1:]) / 2.0
    Fitter = utils.gaussian_fitter_2d(inpdata=data)
    p, success = Fitter.FitGauss2D()
    x, y = np.indices((binsx, binsy), dtype=np.float)
    gauss = Fitter.Gaussian2D(*p)
    data_fitted = gauss(y, x)
    if type(ncont) != int:
        ncont = np.array(ncont) * p[0]
    plt.contour(
        xcenters,
        ycenters,
        data_fitted.reshape(binsx, binsy),
        ncont,
        linewidths=2,
        colors="w",
    )
    return p

plot_gaussian_line staticmethod

plot_gaussian_line(gaussian_fit, binsx, limit)

Plot the major-axis line of a 2D Gaussian fit.

Parameters:

Name	Type	Description	Default
gaussian_fit	sequence	Fitted Gaussian parameters.	required
binsx	int	Histogram bin count used for the fit.	required
limit	sequence	(xmin, xmax, ymin, ymax) plot limits.	required

Source code in lelantos/boxdm.py

@staticmethod
def plot_gaussian_line(gaussian_fit, binsx, limit):
    """Plot the major-axis line of a 2D Gaussian fit.

    Args:
        gaussian_fit (sequence): Fitted Gaussian parameters.
        binsx (int): Histogram bin count used for the fit.
        limit (sequence): ``(xmin, xmax, ymin, ymax)`` plot limits.
    """
    xcenter = (gaussian_fit[1] - binsx) * ((limit[1] - limit[0]) / binsx)
    ycenter = (gaussian_fit[2] - binsx) * ((limit[3] - limit[2]) / binsx)
    xcenter = 0
    ycenter = 0
    angle = gaussian_fit[5]
    pente = 1 / np.tan(np.radians(angle))
    x_arrray = np.linspace(limit[0], limit[1], 100)
    y_array = ycenter + pente * (x_arrray - xcenter)
    plt.plot(x_arrray, y_array)

plot_scatter_los_DM staticmethod

plot_scatter_los_DM(pixel_dm_name, pixel_name, gaussian_smoothing=None)

2D-histogram the reconstructed vs DM delta along lines of sight.

Parameters:

Name	Type	Description	Default
pixel_dm_name	str	DM pixel file.	required
pixel_name	str	Reconstructed pixel file.	required
gaussian_smoothing	float	Unused placeholder.	None

Source code in lelantos/boxdm.py

@staticmethod
def plot_scatter_los_DM(self, pixel_dm_name, pixel_name, gaussian_smoothing=None):
    """2D-histogram the reconstructed vs DM delta along lines of sight.

    Args:
        pixel_dm_name (str): DM pixel file.
        pixel_name (str): Reconstructed pixel file.
        gaussian_smoothing (float, optional): Unused placeholder.
    """
    pixel_dm = tomographic_objects.Pixel(name=pixel_dm_name)
    pixel_dm.read()
    pixel = tomographic_objects.Pixel(name=pixel_name)
    pixel.read()
    delta_pixel_tomo = pixel.pixel_array[:, 4]
    sigma_pixel_tomo = pixel.pixel_array[:, 3]
    mask = sigma_pixel_tomo < 0.3
    delta_pixel_tomo = delta_pixel_tomo[mask]
    pixel_DM = pixel_dm.pixel_array[mask]
    plt.hist2d(delta_pixel_tomo, pixel_DM, 100)
    print(np.corrcoef(delta_pixel_tomo, pixel_DM))

plot_scatter_box_DM

plot_scatter_box_DM(name, gaussian_smoothing=None, cut_redshift_coef=None, cut_pixel=None, binsx=200, binsy=200, ncont=4, rotate=False)

Correlate the reconstructed map with the DM box, voxel by voxel.

Produces density, Gaussian-fit and 2D-histogram figures of the reconstructed vs DM over-density, with a linear fit and correlation.

Parameters:

Name	Type	Description	Default
name	str	Output figure base name.	required
gaussian_smoothing	float	Smoothing scale (Mpc.h^-1).	None
cut_redshift_coef	int	Keep only 1/coef of the z range.	None
cut_pixel	int	Keep only the first cut_pixel slices.	None
binsx, binsy	int	Histogram bin counts.	required
ncont	int	Number of contour levels.	4
rotate	bool	Swap the two axes.	False

Returns:

Name	Type	Description
tuple		(poly_fit, gaussian_fit, correlation_matrix).

Source code in lelantos/boxdm.py

def plot_scatter_box_DM(
    self,
    name,
    gaussian_smoothing=None,
    cut_redshift_coef=None,
    cut_pixel=None,
    binsx=200,
    binsy=200,
    ncont=4,
    rotate=False,
):
    """Correlate the reconstructed map with the DM box, voxel by voxel.

    Produces density, Gaussian-fit and 2D-histogram figures of the
    reconstructed vs DM over-density, with a linear fit and correlation.

    Args:
        name (str): Output figure base name.
        gaussian_smoothing (float, optional): Smoothing scale (Mpc.h^-1).
        cut_redshift_coef (int, optional): Keep only 1/coef of the z range.
        cut_pixel (int, optional): Keep only the first ``cut_pixel`` slices.
        binsx, binsy (int, optional): Histogram bin counts.
        ncont (int, optional): Number of contour levels.
        rotate (bool, optional): Swap the two axes.

    Returns:
        tuple: ``(poly_fit, gaussian_fit, correlation_matrix)``.
    """
    (list_map, list_box) = self.get_list_boxes(
        gaussian_smoothing=gaussian_smoothing,
        cut_redshift_coef=cut_redshift_coef,
        cut_pixel=cut_pixel,
    )
    if gaussian_smoothing is not None:
        if cut_redshift_coef is not None:
            name = name + "_smoothing{}_cutpart{}".format(
                gaussian_smoothing, cut_redshift_coef
            )
        elif cut_pixel is not None:
            name = name + "_smoothing{}_cutpixel{}".format(
                gaussian_smoothing, cut_pixel
            )
        else:
            name = name + "_smoothing{}".format(gaussian_smoothing)
    else:
        if cut_redshift_coef is not None:
            name = name + "_cutpart{}".format(cut_redshift_coef)
        elif cut_pixel is not None:
            name = name + "_cutpixel{}".format(cut_pixel)
        else:
            name = name

    #### Analysis ####

    poly_fit = np.polyfit(list_map, list_box, 1)
    x_poly = np.linspace(-0.2, 0.2, 100)
    y_poly = poly_fit[0] * x_poly + poly_fit[1]
    correlation_matrix = np.corrcoef(list_map, list_box)

    #### Plots ####
    xlabel = r"$\delta_{Fmap}$"
    ylabel = r"$\delta_{m}$"
    if rotate:
        list_box, list_map = list_map, list_box

    plt.figure()
    self.densityplot(
        list_map,
        list_box,
        binsx=binsx,
        binsy=binsy,
        scaleperdeg2=False,
        maxdensity=0,
        mindensity=0,
        pltclf=True,
        logscale=False,
    )
    self.contourplot(
        list_map,
        list_box,
        ncont=ncont,
        colors="w",
        pltclf=False,
        binsx=binsx,
        binsy=binsy,
        log=False,
    )
    plt.xlim(self.limit[0:2])
    plt.ylim(self.limit[2:4])
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    if rotate:
        plt.xlabel(ylabel)
        plt.ylabel(xlabel)
    plt.savefig(name + "_contour.pdf", format="pdf", dpi=100)

    plt.figure()
    self.densityplot(
        list_map,
        list_box,
        binsx=binsx,
        binsy=binsy,
        scaleperdeg2=False,
        maxdensity=0,
        mindensity=0,
        pltclf=True,
        logscale=False,
    )
    gaussian_fit = self.gaussian_plot(
        list_map, list_box, binsx=binsx, binsy=binsy, ncont=ncont
    )
    plt.xlim(self.limit[0:2])
    plt.ylim(self.limit[2:4])
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    if rotate:
        plt.xlabel(ylabel)
        plt.ylabel(xlabel)
    plt.savefig(name + "_gaussian_fit.pdf", format="pdf")

    plt.figure()
    plt.hist2d(list_map, list_box, 300)
    plt.plot(x_poly, y_poly, "r-")
    plt.xlim(self.limit[0:2])
    plt.ylim(self.limit[2:4])
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    if rotate:
        plt.xlabel(ylabel)
        plt.ylabel(xlabel)
    plt.savefig(name + "_histo2d.pdf", format="pdf")

    del list_box, list_map

    return (poly_fit, gaussian_fit, correlation_matrix)

get_list_boxes

get_list_boxes(gaussian_smoothing=None, cut_redshift_coef=None, cut_pixel=None)

Flatten the reconstructed map and DM box into paired 1D arrays.

Parameters:

Name	Type	Description	Default
gaussian_smoothing	float	Smoothing scale (Mpc.h^-1).	None
cut_redshift_coef	int	Keep only 1/coef of the z range.	None
cut_pixel	int	Keep only the first cut_pixel slices.	None

Returns:

Name	Type	Description
tuple		(list_map, list_box) — flattened, aligned voxel values.

Source code in lelantos/boxdm.py

def get_list_boxes(
    self, gaussian_smoothing=None, cut_redshift_coef=None, cut_pixel=None
):
    """Flatten the reconstructed map and DM box into paired 1D arrays.

    Args:
        gaussian_smoothing (float, optional): Smoothing scale (Mpc.h^-1).
        cut_redshift_coef (int, optional): Keep only 1/coef of the z range.
        cut_pixel (int, optional): Keep only the first ``cut_pixel`` slices.

    Returns:
        tuple: ``(list_map, list_box)`` — flattened, aligned voxel values.
    """
    tomography_map = tomographic_objects.TomographicMap.init_from_property_files(
        self.property_file, name=self.map_name
    )
    tomography_map.read()
    box = tomographic_objects.TomographicMap.init_from_property_files(
        self.property_file, name=self.box_name
    )
    box.read()
    max_x, max_y, max_z = (
        tomography_map.shape[0],
        tomography_map.shape[1],
        tomography_map.shape[2],
    )
    if cut_redshift_coef is not None:
        max_z = tomography_map.shape[2] // cut_redshift_coef
    if cut_pixel is not None:
        max_z = cut_pixel
    min_x, min_y, min_z = 0, 0, 0
    length_list = (max_x - min_x) * (max_y - min_y) * (max_z - min_z)
    box_DM = box.map_array[min_x:max_x, min_y:max_y, min_z:max_z]
    if gaussian_smoothing is not None:
        gaussian_smoothing_pix = gaussian_smoothing * utils.pixel_per_mpc(
            box.size, box.shape
        )
        box_DM = utils.gaussian_smoothing(box_DM, gaussian_smoothing_pix)
    list_box = box_DM.reshape(length_list)
    del box_DM, box
    map_3D = tomography_map.map_array[min_x:max_x, min_y:max_y, min_z:max_z]
    list_map = map_3D.reshape(length_list)
    del map_3D, tomography_map
    return (list_map, list_box)