API

Base reconstruction class

Implementation of base reconstruction class.

class pyrecon.recon.BaseReconstruction(f=None, bias=None, los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, data_positions=None, randoms_positions=None, data_weights=None, randoms_weights=None, positions=None, position_type='pos', resampler='cic', decomposition=None, fft_plan='estimate', dtype='f8', mpiroot=None, mpicomm=<mpi4py.MPI.Intracomm object>, threshold_randoms=0.01, **kwargs)

Bases: BaseClass

Base template reconstruction class. Reconstruction algorithms should extend this class, by (at least) implementing:

A standard reconstruction would be:

# MyReconstruction is your reconstruction algorithm
recon = MyReconstruction(f=0.8, bias=2.0, nmesh=512, boxsize=1000., boxcenter=2000.)
recon.assign_data(positions_data, weights_data)
recon.assign_randoms(positions_randoms, weights_randoms)
recon.set_density_contrast()
recon.run()
positions_rec_data = recon.read_shifted_positions(positions_data)
# RecSym = remove large scale RSD from randoms
positions_rec_randoms = recon.read_shifted_positions(positions_randoms)
# Or RecIso
# positions_rec_randoms = recon.read_shifted_positions(positions_randoms, field='disp')
mesh_data

Mesh (3D grid) to assign (“paint”) galaxies to.

Type:

RealField

mesh_randoms

Mesh (3D grid) to assign (“paint”) randoms to.

Type:

RealField

f

Growth rate.

Type:

float

bias

Galaxy bias.

Type:

float

los

If None, local (varying) line-of-sight. Else line-of-sight (unit) 3-vector.

Type:

array

boxsize, boxcenter, cellsize, offset

Mesh properties.

Type:

array

Initialize BaseReconstruction.

Parameters:

  • f (float, default=None.) – Growth rate.

  • bias (float, default=None.) – Galaxy bias.

  • los (string, array_like, default=None) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis.

  • boxsize (array, float, default=None) – Physical size of the box along each axis, defaults to maximum extent taken by all input positions, times boxpad.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • boxpad (float, default=2.) – When boxsize is determined from input positions, take boxpad times the smallest box enclosing positions as boxsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize]? If False and input positions do not fit in the the box size, raise a ValueError.

  • positions (list, array, default=None) – Optionally, positions used to defined box size. Of shape (3, N) or (N, 3), depending on position_type.

  • data_positions (list, array, default=None) – Positions in the data catalog. Of shape (3, N) or (N, 3), depending on position_type. If provided, reconstruction will be run directly (and data_positions will be added to positions to define the box size).

  • randoms_positions (list, array, default=None) – Positions in the randoms catalog. See data_positions.

  • data_weights (array of shape (N,), default=None) – Optionally, weights in the data catalog.

  • randoms_weights (array of shape (N,), default=None) – Optionally, weights in the randoms catalog.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

  • fft_plan (string, default='estimate') – FFT planning. ‘measure’ may allow for faster FFTs, but is slower to set up than ‘estimate’.

  • resampler (string, ResampleWindow, default='tsc') – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’].

  • dtype (string, dtype, default='f8') – The data type to use for the mesh.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scattered across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

  • threshold_randoms (float, default=0.01) – mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

assign_data(positions, weights=None, **kwargs)

Assign (paint) data to mesh_data. This can be done slab-by-slab (e.g. to reduce memory footprint).

Parameters:

  • positions (array of shape (N, 3)) – Cartesian positions.

  • weights (array of shape (N,), default=None) – Weights; default to 1.

assign_randoms(positions, weights=None)

Same as assign_data(), but for random objects.

property beta

\(\beta\) parameter, as the ratio of the growth rate to the galaxy bias.

read_shifted_positions(positions, field='disp+rsd')

Read shifted positions i.e. the difference positions - self.read_shifts(positions, field=field). Output (and input) positions are wrapped if wrap.

Parameters:

  • positions (list, array) – Positions of shape (3, N) or (N, 3), depending on position_type.

  • field (string, default='disp+rsd') – Apply either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

  • kwargs (dict) – position_type, mpiroot can be provided to override default position_type, mpiroot.

Returns:

positions – Shifted positions, of same type as input.

Return type:

list, array

read_shifts(positions, field='disp+rsd')

Read displacement at input positions. To get shifted/reconstructed positions, given reconstruction instance recon:

positions_rec_data = positions_data - recon.read_shifts(positions_data)
# RecSym = remove large scale RSD from randoms
positions_rec_randoms = positions_randoms - recon.read_shifts(positions_randoms)
# Or RecIso
# positions_rec_randoms = positions_randoms - recon.read_shifts(positions_randoms, field='disp')

Or directly use read_shifted_positions() (which wraps output positions if wrap).

Parameters:

  • positions (array of shape (N, 3)) – Cartesian positions.

  • field (string, default='disp+rsd') – Either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

  • kwargs (dict) – position_type, mpiroot can be provided to override default position_type, mpiroot.

Returns:

shifts – Displacements.

Return type:

array of shape (N, 3)

run(*args, **kwargs)

Run reconstruction; to be implemented in your algorithm.

set_bias(bias)

Set bias. If callable, set real space mesh bias.

set_cosmo(f=None, bias=None, beta=None)

Set cosmology.

Parameters:

  • f (float) – Growth rate. If None and beta` is provided, set f as beta * bias; else growth rate is left unchanged.

  • bias (float) – Bias. If None, bias is left unchanged.

  • beta (float) – \(\beta\) parameter. If not None, overrides f as beta * bias.

set_density_contrast(threshold_randoms=None, smoothing_radius=15.0, check=False, kw_weights=None)

Set \(\delta\) field mesh_delta from data and randoms fields mesh_data and mesh_randoms.

Note

This method follows Julian’s reconstruction code.

Parameters:

  • threshold_randoms (float, default=0.01) – If provided, override value given at initialization. mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

  • smoothing_radius (float, default=15) – Smoothing scale.

  • check (bool, default=False) – If True, run some tests (printed in logger) to assess whether enough randoms have been used.

  • kw_weights (dict, default=None) – Optionally, for optimal weights, dictionary with keys ‘P0’ (fiducial power spectrum value), and ‘nbar’ (comoving density at given comoving distance). If ‘nbar’ not provided, mesh_randoms rescaled to the data density is used instead.

set_f(f)

Set growth rate. If callable, set real space mesh f.

set_los(los=None)

Set line-of-sight.

Parameters:

los (string, array) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

set_threshold_randoms(threshold_randoms=0.01)

Set threshold on randoms to define the density contrast.

exception pyrecon.recon.ReconstructionError

Bases: Exception

Error raised when issue with reconstruction.

with_traceback()

Exception.with_traceback(tb) – set self.__traceback__ to tb and return self.

pyrecon.recon.format_positions_weights_wrapper(func)

Method wrapper applying _format_positions and _format_weigths on input.

Multi grid reconstruction

Re-implementation of Martin J. White’s reconstruction code.

class pyrecon.multigrid.MultiGridReconstruction(*args, mpicomm=<mpi4py.MPI.Intracomm object>, **kwargs)

Bases: OriginalMultiGridReconstruction

Any update / test / improvement upon original algorithm.

Initialize BaseReconstruction.

Parameters:

  • f (float, default=None.) – Growth rate.

  • bias (float, default=None.) – Galaxy bias.

  • los (string, array_like, default=None) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis.

  • boxsize (array, float, default=None) – Physical size of the box along each axis, defaults to maximum extent taken by all input positions, times boxpad.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • boxpad (float, default=2.) – When boxsize is determined from input positions, take boxpad times the smallest box enclosing positions as boxsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize]? If False and input positions do not fit in the the box size, raise a ValueError.

  • positions (list, array, default=None) – Optionally, positions used to defined box size. Of shape (3, N) or (N, 3), depending on position_type.

  • data_positions (list, array, default=None) – Positions in the data catalog. Of shape (3, N) or (N, 3), depending on position_type. If provided, reconstruction will be run directly (and data_positions will be added to positions to define the box size).

  • randoms_positions (list, array, default=None) – Positions in the randoms catalog. See data_positions.

  • data_weights (array of shape (N,), default=None) – Optionally, weights in the data catalog.

  • randoms_weights (array of shape (N,), default=None) – Optionally, weights in the randoms catalog.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

  • fft_plan (string, default='estimate') – FFT planning. ‘measure’ may allow for faster FFTs, but is slower to set up than ‘estimate’.

  • resampler (string, ResampleWindow, default='tsc') – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’].

  • dtype (string, dtype, default='f8') – The data type to use for the mesh.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scattered across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

  • threshold_randoms (float, default=0.01) – mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

assign_data(positions, weights=None, **kwargs)

Assign (paint) data to mesh_data. This can be done slab-by-slab (e.g. to reduce memory footprint).

Parameters:

  • positions (array of shape (N, 3)) – Cartesian positions.

  • weights (array of shape (N,), default=None) – Weights; default to 1.

assign_randoms(positions, weights=None)

Same as assign_data(), but for random objects.

property beta

\(\beta\) parameter, as the ratio of the growth rate to the galaxy bias.

read_shifted_positions(positions, field='disp+rsd')

Read shifted positions i.e. the difference positions - self.read_shifts(positions, field=field). Output (and input) positions are wrapped if wrap.

Parameters:

  • positions (list, array) – Positions of shape (3, N) or (N, 3), depending on position_type.

  • field (string, default='disp+rsd') – Apply either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

  • kwargs (dict) – position_type, mpiroot can be provided to override default position_type, mpiroot.

Returns:

positions – Shifted positions, of same type as input.

Return type:

list, array

read_shifts(positions, field='disp+rsd')

Read displacement at input positions by deriving the computed displacement potential mesh_phi (finite difference scheme). See BaseReconstruction.read_shifts() for input parameters.

run(jacobi_damping_factor=0.4, jacobi_niterations=5, vcycle_niterations=6)

Run reconstruction, i.e. set displacement potential attr:mesh_phi from mesh_delta. Default parameter values are the same as in Martin’s code.

Parameters:

  • jacobi_damping_factor (float, default=0.4) – Damping factor for Jacobi iterations.

  • jacobi_niterations (int, default=5) – Number of Jacobi iterations.

  • vcycle_niterations (int, default=6) – Number of V-cycle calls.

set_bias(bias)

Set bias. If callable, set real space mesh bias.

set_cosmo(f=None, bias=None, beta=None)

Set cosmology.

Parameters:

  • f (float) – Growth rate. If None and beta` is provided, set f as beta * bias; else growth rate is left unchanged.

  • bias (float) – Bias. If None, bias is left unchanged.

  • beta (float) – \(\beta\) parameter. If not None, overrides f as beta * bias.

set_density_contrast(*args, **kwargs)

See BaseReconstruction.set_density_contrast.

set_f(f)

Set growth rate. If callable, set real space mesh f.

set_los(los=None)

Set line-of-sight.

Parameters:

los (string, array) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

set_threshold_randoms(threshold_randoms=0.01)

Set threshold on randoms to define the density contrast.

class pyrecon.multigrid.OriginalMultiGridReconstruction(*args, mpicomm=<mpi4py.MPI.Intracomm object>, **kwargs)

Bases: BaseReconstruction

ctypes-based implementation for Martin J. White’s reconstruction code, using full multigrid V-cycle based on damped Jacobi iteration. We re-implemented https://github.com/martinjameswhite/recon_code/blob/master/multigrid.cpp, allowing for non-cubic (rectangular) mesh. Numerical agreement in the Zeldovich displacements between original code and this re-implementation is numerical precision (absolute and relative difference of 1e-10). To test this, change float to double and increase precision in io.cpp/write_data in the original code.

Initialize BaseReconstruction.

Parameters:

  • f (float, default=None.) – Growth rate.

  • bias (float, default=None.) – Galaxy bias.

  • los (string, array_like, default=None) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis.

  • boxsize (array, float, default=None) – Physical size of the box along each axis, defaults to maximum extent taken by all input positions, times boxpad.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • boxpad (float, default=2.) – When boxsize is determined from input positions, take boxpad times the smallest box enclosing positions as boxsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize]? If False and input positions do not fit in the the box size, raise a ValueError.

  • positions (list, array, default=None) – Optionally, positions used to defined box size. Of shape (3, N) or (N, 3), depending on position_type.

  • data_positions (list, array, default=None) – Positions in the data catalog. Of shape (3, N) or (N, 3), depending on position_type. If provided, reconstruction will be run directly (and data_positions will be added to positions to define the box size).

  • randoms_positions (list, array, default=None) – Positions in the randoms catalog. See data_positions.

  • data_weights (array of shape (N,), default=None) – Optionally, weights in the data catalog.

  • randoms_weights (array of shape (N,), default=None) – Optionally, weights in the randoms catalog.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

  • fft_plan (string, default='estimate') – FFT planning. ‘measure’ may allow for faster FFTs, but is slower to set up than ‘estimate’.

  • resampler (string, ResampleWindow, default='tsc') – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’].

  • dtype (string, dtype, default='f8') – The data type to use for the mesh.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scattered across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

  • threshold_randoms (float, default=0.01) – mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

assign_data(positions, weights=None, **kwargs)

Assign (paint) data to mesh_data. This can be done slab-by-slab (e.g. to reduce memory footprint).

Parameters:

  • positions (array of shape (N, 3)) – Cartesian positions.

  • weights (array of shape (N,), default=None) – Weights; default to 1.

assign_randoms(positions, weights=None)

Same as assign_data(), but for random objects.

property beta

\(\beta\) parameter, as the ratio of the growth rate to the galaxy bias.

read_shifted_positions(positions, field='disp+rsd')

Read shifted positions i.e. the difference positions - self.read_shifts(positions, field=field). Output (and input) positions are wrapped if wrap.

Parameters:

  • positions (list, array) – Positions of shape (3, N) or (N, 3), depending on position_type.

  • field (string, default='disp+rsd') – Apply either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

  • kwargs (dict) – position_type, mpiroot can be provided to override default position_type, mpiroot.

Returns:

positions – Shifted positions, of same type as input.

Return type:

list, array

read_shifts(positions, field='disp+rsd')

Read displacement at input positions by deriving the computed displacement potential mesh_phi (finite difference scheme). See BaseReconstruction.read_shifts() for input parameters.

run(jacobi_damping_factor=0.4, jacobi_niterations=5, vcycle_niterations=6)

Run reconstruction, i.e. set displacement potential attr:mesh_phi from mesh_delta. Default parameter values are the same as in Martin’s code.

Parameters:

  • jacobi_damping_factor (float, default=0.4) – Damping factor for Jacobi iterations.

  • jacobi_niterations (int, default=5) – Number of Jacobi iterations.

  • vcycle_niterations (int, default=6) – Number of V-cycle calls.

set_bias(bias)

Set bias. If callable, set real space mesh bias.

set_cosmo(f=None, bias=None, beta=None)

Set cosmology.

Parameters:

  • f (float) – Growth rate. If None and beta` is provided, set f as beta * bias; else growth rate is left unchanged.

  • bias (float) – Bias. If None, bias is left unchanged.

  • beta (float) – \(\beta\) parameter. If not None, overrides f as beta * bias.

set_density_contrast(ran_min=0.75, smoothing_radius=15.0, **kwargs)

Set \(\delta\) field mesh_delta from data and randoms fields mesh_data and mesh_randoms.

Note

This method follows Martin’s reconstruction code: we decided to change the ran_min prescription. At least ran_min should depend on random weights. See also Martin’s notes below.

Parameters:

  • ran_min (float, default=0.75) – mesh_randoms points below this threshold have their density contrast set to 0.

  • smoothing_radius (float, default=15) – Smoothing scale, see RealMesh.smooth_gaussian().

  • kwargs (dict) – Optional arguments for RealMesh.smooth_gaussian().

set_f(f)

Set growth rate. If callable, set real space mesh f.

set_los(los=None)

Set line-of-sight.

Parameters:

los (string, array) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

set_threshold_randoms(threshold_randoms=0.01)

Set threshold on randoms to define the density contrast.

Iterative FFT reconstruction

Implementation of Burden et al. 2015 (https://arxiv.org/abs/1504.02591) algorithm.

class pyrecon.iterative_fft.IterativeFFTReconstruction(f=None, bias=None, los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, data_positions=None, randoms_positions=None, data_weights=None, randoms_weights=None, positions=None, position_type='pos', resampler='cic', decomposition=None, fft_plan='estimate', dtype='f8', mpiroot=None, mpicomm=<mpi4py.MPI.Intracomm object>, threshold_randoms=0.01, **kwargs)

Bases: BaseReconstruction

Implementation of Burden et al. 2015 (https://arxiv.org/abs/1504.02591) field-level (as opposed to IterativeFFTParticleReconstruction) algorithm.

Initialize BaseReconstruction.

Parameters:

  • f (float, default=None.) – Growth rate.

  • bias (float, default=None.) – Galaxy bias.

  • los (string, array_like, default=None) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis.

  • boxsize (array, float, default=None) – Physical size of the box along each axis, defaults to maximum extent taken by all input positions, times boxpad.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • boxpad (float, default=2.) – When boxsize is determined from input positions, take boxpad times the smallest box enclosing positions as boxsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize]? If False and input positions do not fit in the the box size, raise a ValueError.

  • positions (list, array, default=None) – Optionally, positions used to defined box size. Of shape (3, N) or (N, 3), depending on position_type.

  • data_positions (list, array, default=None) – Positions in the data catalog. Of shape (3, N) or (N, 3), depending on position_type. If provided, reconstruction will be run directly (and data_positions will be added to positions to define the box size).

  • randoms_positions (list, array, default=None) – Positions in the randoms catalog. See data_positions.

  • data_weights (array of shape (N,), default=None) – Optionally, weights in the data catalog.

  • randoms_weights (array of shape (N,), default=None) – Optionally, weights in the randoms catalog.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

  • fft_plan (string, default='estimate') – FFT planning. ‘measure’ may allow for faster FFTs, but is slower to set up than ‘estimate’.

  • resampler (string, ResampleWindow, default='tsc') – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’].

  • dtype (string, dtype, default='f8') – The data type to use for the mesh.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scattered across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

  • threshold_randoms (float, default=0.01) – mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

assign_data(positions, weights=None, **kwargs)

Assign (paint) data to mesh_data. This can be done slab-by-slab (e.g. to reduce memory footprint).

Parameters:

  • positions (array of shape (N, 3)) – Cartesian positions.

  • weights (array of shape (N,), default=None) – Weights; default to 1.

assign_randoms(positions, weights=None)

Same as assign_data(), but for random objects.

property beta

\(\beta\) parameter, as the ratio of the growth rate to the galaxy bias.

read_shifted_positions(positions, field='disp+rsd')

Read shifted positions i.e. the difference positions - self.read_shifts(positions, field=field). Output (and input) positions are wrapped if wrap.

Parameters:

  • positions (list, array) – Positions of shape (3, N) or (N, 3), depending on position_type.

  • field (string, default='disp+rsd') – Apply either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

  • kwargs (dict) – position_type, mpiroot can be provided to override default position_type, mpiroot.

Returns:

positions – Shifted positions, of same type as input.

Return type:

list, array

read_shifts(positions, field='disp+rsd')

Read displacement at input positions. To get shifted/reconstructed positions, given reconstruction instance recon:

positions_rec_data = positions_data - recon.read_shifts(positions_data)
# RecSym = remove large scale RSD from randoms
positions_rec_randoms = positions_randoms - recon.read_shifts(positions_randoms)
# Or RecIso
# positions_rec_randoms = positions_randoms - recon.read_shifts(positions_randoms, field='disp')

Or directly use read_shifted_positions() (which wraps output positions if wrap).

Parameters:

  • positions (array of shape (N, 3)) – Cartesian positions.

  • field (string, default='disp+rsd') – Either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

  • kwargs (dict) – position_type, mpiroot can be provided to override default position_type, mpiroot.

Returns:

shifts – Displacements.

Return type:

array of shape (N, 3)

run(niterations=3)

Run reconstruction, i.e. compute Zeldovich displacement fields mesh_psi.

Parameters:

niterations (int, default=3) – Number of iterations.

set_bias(bias)

Set bias. If callable, set real space mesh bias.

set_cosmo(f=None, bias=None, beta=None)

Set cosmology.

Parameters:

  • f (float) – Growth rate. If None and beta` is provided, set f as beta * bias; else growth rate is left unchanged.

  • bias (float) – Bias. If None, bias is left unchanged.

  • beta (float) – \(\beta\) parameter. If not None, overrides f as beta * bias.

set_density_contrast(threshold_randoms=None, smoothing_radius=15.0, check=False, kw_weights=None)

Set \(\delta\) field mesh_delta from data and randoms fields mesh_data and mesh_randoms.

Note

This method follows Julian’s reconstruction code.

Parameters:

  • threshold_randoms (float, default=0.01) – If provided, override value given at initialization. mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

  • smoothing_radius (float, default=15) – Smoothing scale.

  • check (bool, default=False) – If True, run some tests (printed in logger) to assess whether enough randoms have been used.

  • kw_weights (dict, default=None) – Optionally, for optimal weights, dictionary with keys ‘P0’ (fiducial power spectrum value), and ‘nbar’ (comoving density at given comoving distance). If ‘nbar’ not provided, mesh_randoms rescaled to the data density is used instead.

set_f(f)

Set growth rate. If callable, set real space mesh f.

set_los(los=None)

Set line-of-sight.

Parameters:

los (string, array) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

set_threshold_randoms(threshold_randoms=0.01)

Set threshold on randoms to define the density contrast.

Iterative FFT particle reconstruction

Re-implementation of Bautista et al. 2018 (https://arxiv.org/pdf/1712.08064.pdf) algorithm.

class pyrecon.iterative_fft_particle.IterativeFFTParticleReconstruction(f=None, bias=None, los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, data_positions=None, randoms_positions=None, data_weights=None, randoms_weights=None, positions=None, position_type='pos', resampler='cic', decomposition=None, fft_plan='estimate', dtype='f8', mpiroot=None, mpicomm=<mpi4py.MPI.Intracomm object>, threshold_randoms=0.01, **kwargs)

Bases: OriginalIterativeFFTParticleReconstruction

Any update / test / improvement upon original algorithm.

Initialize BaseReconstruction.

Parameters:

  • f (float, default=None.) – Growth rate.

  • bias (float, default=None.) – Galaxy bias.

  • los (string, array_like, default=None) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis.

  • boxsize (array, float, default=None) – Physical size of the box along each axis, defaults to maximum extent taken by all input positions, times boxpad.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • boxpad (float, default=2.) – When boxsize is determined from input positions, take boxpad times the smallest box enclosing positions as boxsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize]? If False and input positions do not fit in the the box size, raise a ValueError.

  • positions (list, array, default=None) – Optionally, positions used to defined box size. Of shape (3, N) or (N, 3), depending on position_type.

  • data_positions (list, array, default=None) – Positions in the data catalog. Of shape (3, N) or (N, 3), depending on position_type. If provided, reconstruction will be run directly (and data_positions will be added to positions to define the box size).

  • randoms_positions (list, array, default=None) – Positions in the randoms catalog. See data_positions.

  • data_weights (array of shape (N,), default=None) – Optionally, weights in the data catalog.

  • randoms_weights (array of shape (N,), default=None) – Optionally, weights in the randoms catalog.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

  • fft_plan (string, default='estimate') – FFT planning. ‘measure’ may allow for faster FFTs, but is slower to set up than ‘estimate’.

  • resampler (string, ResampleWindow, default='tsc') – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’].

  • dtype (string, dtype, default='f8') – The data type to use for the mesh.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scattered across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

  • threshold_randoms (float, default=0.01) – mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

assign_data(positions, weights=None)

Assign (paint) data to mesh_data. Keeps track of input positions (for run()) and weights (for set_density_contrast()). See BaseReconstruction.assign_data() for parameters.

assign_randoms(positions, weights=None)

Same as assign_data(), but for random objects.

property beta

\(\beta\) parameter, as the ratio of the growth rate to the galaxy bias.

read_shifted_positions(positions, field='disp+rsd')

Read shifted positions i.e. the difference positions - self.read_shifts(positions, field=field). Output (and input) positions are wrapped if wrap.

Parameters:

  • positions (array of shape (N, 3), string) – Cartesian positions. Pass string ‘data’ to get the shift positions for the input data positions passed to assign_data(). Note that in this case, shifts are read at the reconstructed data real-space positions.

  • field (string, default='disp+rsd') – Apply either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

Returns:

positions – Shifted positions.

Return type:

array of shape (N, 3)

read_shifts(positions, field='disp+rsd')

Read displacement at input positions.

Note

Data shifts are read at the reconstructed real-space positions, while random shifts are read at the redshift-space positions, is that consistent?

Parameters:

  • positions (array of shape (N, 3), string) – Cartesian positions. Pass string ‘data’ to get the displacements for the input data positions passed to assign_data(). Note that in this case, shifts are read at the reconstructed data real-space positions.

  • field (string, default='disp+rsd') – Either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

Returns:

shifts – Displacements.

Return type:

array of shape (N, 3)

run(niterations=3)

Run reconstruction, i.e. compute reconstructed data real-space positions (_positions_rec_data) and Zeldovich displacements fields mesh_psi.

Parameters:

niterations (int) – Number of iterations.

set_bias(bias)

Set bias. If callable, set real space mesh bias.

set_cosmo(f=None, bias=None, beta=None)

Set cosmology.

Parameters:

  • f (float) – Growth rate. If None and beta` is provided, set f as beta * bias; else growth rate is left unchanged.

  • bias (float) – Bias. If None, bias is left unchanged.

  • beta (float) – \(\beta\) parameter. If not None, overrides f as beta * bias.

set_density_contrast(threshold_randoms=None, smoothing_radius=15.0, check=False, kw_weights=None)

Set \(\delta\) field mesh_delta from data and randoms fields mesh_data and mesh_randoms.

Note

This method follows Julian’s reconstruction code. mesh_data and mesh_randoms fields are assumed to be smoothed already.

Parameters:

  • threshold_randoms (float, default=0.01) – If provided, override value given at initialization. mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

  • smoothing_radius (float, default=15) – Smoothing scale, see RealMesh.smooth_gaussian().

  • check (bool, default=False) – If True, run some tests (printed in logger) to assess whether enough randoms have been used.

set_f(f)

Set growth rate. If callable, set real space mesh f.

set_los(los=None)

Set line-of-sight.

Parameters:

los (string, array) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

set_threshold_randoms(threshold_randoms=0.01)

Set threshold on randoms to define the density contrast.

class pyrecon.iterative_fft_particle.OriginalIterativeFFTParticleReconstruction(f=None, bias=None, los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, data_positions=None, randoms_positions=None, data_weights=None, randoms_weights=None, positions=None, position_type='pos', resampler='cic', decomposition=None, fft_plan='estimate', dtype='f8', mpiroot=None, mpicomm=<mpi4py.MPI.Intracomm object>, threshold_randoms=0.01, **kwargs)

Bases: BaseReconstruction

Exact re-implementation of Bautista et al. 2018 (https://arxiv.org/pdf/1712.08064.pdf) algorithm at https://github.com/julianbautista/eboss_clustering/blob/master/python/recon.py. Numerical agreement in the Zeldovich displacements between original codes and this re-implementation is machine precision (absolute and relative difference of 1e-12).

Initialize BaseReconstruction.

Parameters:

  • f (float, default=None.) – Growth rate.

  • bias (float, default=None.) – Galaxy bias.

  • los (string, array_like, default=None) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis.

  • boxsize (array, float, default=None) – Physical size of the box along each axis, defaults to maximum extent taken by all input positions, times boxpad.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • boxpad (float, default=2.) – When boxsize is determined from input positions, take boxpad times the smallest box enclosing positions as boxsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize]? If False and input positions do not fit in the the box size, raise a ValueError.

  • positions (list, array, default=None) – Optionally, positions used to defined box size. Of shape (3, N) or (N, 3), depending on position_type.

  • data_positions (list, array, default=None) – Positions in the data catalog. Of shape (3, N) or (N, 3), depending on position_type. If provided, reconstruction will be run directly (and data_positions will be added to positions to define the box size).

  • randoms_positions (list, array, default=None) – Positions in the randoms catalog. See data_positions.

  • data_weights (array of shape (N,), default=None) – Optionally, weights in the data catalog.

  • randoms_weights (array of shape (N,), default=None) – Optionally, weights in the randoms catalog.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

  • fft_plan (string, default='estimate') – FFT planning. ‘measure’ may allow for faster FFTs, but is slower to set up than ‘estimate’.

  • resampler (string, ResampleWindow, default='tsc') – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’].

  • dtype (string, dtype, default='f8') – The data type to use for the mesh.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scattered across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

  • threshold_randoms (float, default=0.01) – mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

assign_data(positions, weights=None)

Assign (paint) data to mesh_data. Keeps track of input positions (for run()) and weights (for set_density_contrast()). See BaseReconstruction.assign_data() for parameters.

assign_randoms(positions, weights=None)

Same as assign_data(), but for random objects.

property beta

\(\beta\) parameter, as the ratio of the growth rate to the galaxy bias.

read_shifted_positions(positions, field='disp+rsd')

Read shifted positions i.e. the difference positions - self.read_shifts(positions, field=field). Output (and input) positions are wrapped if wrap.

Parameters:

  • positions (array of shape (N, 3), string) – Cartesian positions. Pass string ‘data’ to get the shift positions for the input data positions passed to assign_data(). Note that in this case, shifts are read at the reconstructed data real-space positions.

  • field (string, default='disp+rsd') – Apply either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

Returns:

positions – Shifted positions.

Return type:

array of shape (N, 3)

read_shifts(positions, field='disp+rsd')

Read displacement at input positions.

Note

Data shifts are read at the reconstructed real-space positions, while random shifts are read at the redshift-space positions, is that consistent?

Parameters:

  • positions (array of shape (N, 3), string) – Cartesian positions. Pass string ‘data’ to get the displacements for the input data positions passed to assign_data(). Note that in this case, shifts are read at the reconstructed data real-space positions.

  • field (string, default='disp+rsd') – Either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

Returns:

shifts – Displacements.

Return type:

array of shape (N, 3)

run(niterations=3)

Run reconstruction, i.e. compute reconstructed data real-space positions (_positions_rec_data) and Zeldovich displacements fields mesh_psi.

Parameters:

niterations (int) – Number of iterations.

set_bias(bias)

Set bias. If callable, set real space mesh bias.

set_cosmo(f=None, bias=None, beta=None)

Set cosmology.

Parameters:

  • f (float) – Growth rate. If None and beta` is provided, set f as beta * bias; else growth rate is left unchanged.

  • bias (float) – Bias. If None, bias is left unchanged.

  • beta (float) – \(\beta\) parameter. If not None, overrides f as beta * bias.

set_density_contrast(threshold_randoms=None, smoothing_radius=15.0, check=False, kw_weights=None)

Set \(\delta\) field mesh_delta from data and randoms fields mesh_data and mesh_randoms.

Note

This method follows Julian’s reconstruction code. mesh_data and mesh_randoms fields are assumed to be smoothed already.

Parameters:

  • threshold_randoms (float, default=0.01) – If provided, override value given at initialization. mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

  • smoothing_radius (float, default=15) – Smoothing scale, see RealMesh.smooth_gaussian().

  • check (bool, default=False) – If True, run some tests (printed in logger) to assess whether enough randoms have been used.

set_f(f)

Set growth rate. If callable, set real space mesh f.

set_los(los=None)

Set line-of-sight.

Parameters:

los (string, array) – If los is None, use local (varying) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

set_threshold_randoms(threshold_randoms=0.01)

Set threshold on randoms to define the density contrast.

Plane-parallel FFT reconstruction

Implementation of Burden et al. 2015 (https://arxiv.org/abs/1504.02591) algorithm.

class pyrecon.plane_parallel_fft.PlaneParallelFFTReconstruction(los=None, **kwargs)

Bases: BaseReconstruction

Implementation of Eisenstein et al. 2007 (https://arxiv.org/pdf/astro-ph/0604362.pdf) algorithm. Section 3, paragraph starting with ‘Restoring in full the …’

Initialize IterativeFFTReconstruction.

Parameters:

  • los (string, array, default=None) – May be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • kwargs (dict) – See BaseReconstruction for parameters.

assign_data(positions, weights=None, **kwargs)

Assign (paint) data to mesh_data. This can be done slab-by-slab (e.g. to reduce memory footprint).

Parameters:

  • positions (array of shape (N, 3)) – Cartesian positions.

  • weights (array of shape (N,), default=None) – Weights; default to 1.

assign_randoms(positions, weights=None)

Same as assign_data(), but for random objects.

property beta

\(\beta\) parameter, as the ratio of the growth rate to the galaxy bias.

read_shifted_positions(positions, field='disp+rsd')

Read shifted positions i.e. the difference positions - self.read_shifts(positions, field=field). Output (and input) positions are wrapped if wrap.

Parameters:

  • positions (list, array) – Positions of shape (3, N) or (N, 3), depending on position_type.

  • field (string, default='disp+rsd') – Apply either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

  • kwargs (dict) – position_type, mpiroot can be provided to override default position_type, mpiroot.

Returns:

positions – Shifted positions, of same type as input.

Return type:

list, array

read_shifts(positions, field='disp+rsd')

Read displacement at input positions. To get shifted/reconstructed positions, given reconstruction instance recon:

positions_rec_data = positions_data - recon.read_shifts(positions_data)
# RecSym = remove large scale RSD from randoms
positions_rec_randoms = positions_randoms - recon.read_shifts(positions_randoms)
# Or RecIso
# positions_rec_randoms = positions_randoms - recon.read_shifts(positions_randoms, field='disp')

Or directly use read_shifted_positions() (which wraps output positions if wrap).

Parameters:

  • positions (array of shape (N, 3)) – Cartesian positions.

  • field (string, default='disp+rsd') – Either ‘disp’ (Zeldovich displacement), ‘rsd’ (RSD displacement), or ‘disp+rsd’ (Zeldovich + RSD displacement).

  • kwargs (dict) – position_type, mpiroot can be provided to override default position_type, mpiroot.

Returns:

shifts – Displacements.

Return type:

array of shape (N, 3)

run()

Run reconstruction, i.e. compute Zeldovich displacement fields mesh_psi.

set_bias(bias)

Set bias. If callable, set real space mesh bias.

set_cosmo(f=None, bias=None, beta=None)

Set cosmology.

Parameters:

  • f (float) – Growth rate. If None and beta` is provided, set f as beta * bias; else growth rate is left unchanged.

  • bias (float) – Bias. If None, bias is left unchanged.

  • beta (float) – \(\beta\) parameter. If not None, overrides f as beta * bias.

set_density_contrast(threshold_randoms=None, smoothing_radius=15.0, check=False, kw_weights=None)

Set \(\delta\) field mesh_delta from data and randoms fields mesh_data and mesh_randoms.

Note

This method follows Julian’s reconstruction code.

Parameters:

  • threshold_randoms (float, default=0.01) – If provided, override value given at initialization. mesh_randoms points below this threshold times mean random weights have their density contrast set to 0. For a more consistent thresholding, pass e.g. (‘noise’, 0.01) to set the threshold has \(0.01 \sum w^2 / \sum w\) where \(\sum w\), \(\sum w^2\) are the (total) sum of random weights and squared random weights.

  • smoothing_radius (float, default=15) – Smoothing scale.

  • check (bool, default=False) – If True, run some tests (printed in logger) to assess whether enough randoms have been used.

  • kw_weights (dict, default=None) – Optionally, for optimal weights, dictionary with keys ‘P0’ (fiducial power spectrum value), and ‘nbar’ (comoving density at given comoving distance). If ‘nbar’ not provided, mesh_randoms rescaled to the data density is used instead.

set_f(f)

Set growth rate. If callable, set real space mesh f.

set_los(los=None)

Set line-of-sight.

Parameters:

los (string, array_like) – May be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

set_threshold_randoms(threshold_randoms=0.01)

Set threshold on randoms to define the density contrast.

Metrics

Routines to estimate reconstruction efficiency:

This requires the following packages:

class pyrecon.metrics.BasePowerRatio

Bases: BaseClass

Base template class to compute power ratios. Specific statistic should extend this class.

get_ratio(complex=False, **kwargs)

Return power spectrum ratio, computed using various options.

Parameters:

  • complex (bool, default=False) – Whether (True) to return the ratio of complex power spectra, or (False) return the ratio of their real part only.

  • kwargs (dict) – Optionally, arguments for BasePowerSpectrumStatistics.get_power().

  • Results

  • -------

  • ratio (array)

modeavg(axis=0, method=None)

Return average of modes for input axis.

Parameters:

  • axis (int, default=0) – Axis.

  • method (str, default=None) – If None, return average separation from modes. If ‘mid’, return bin mid-points.

Returns:

modeavg – 1D array of size shape[axis].

Return type:

array

property ratio

Power spectrum ratio.

rebin(factor=1)

Rebin statistic, by factor(s) factor. Input factors must divide shape.

save(filename)

Save to filename.

save_txt(filename, fmt='%.12e', delimiter=' ', header=None, comments='# ', **kwargs)

Save power spectrum ratio as txt file.

Warning

Attributes are not all saved, hence there is load_txt() method.

Parameters:

  • filename (str) – File name.

  • fmt (str, default='%.12e') – Format for floating types.

  • delimiter (str, default=' ') – String or character separating columns.

  • header (str, list, default=None) – String that will be written at the beginning of the file. If multiple lines, provide a list of one-line strings.

  • comments (str, default=' #') – String that will be prepended to the header string.

  • kwargs (dict) – Arguments for get_power().

select(*xlims)

Restrict statistic to provided coordinate limits in place.

For example:

statistic.select((0, 0.3)) # restrict first axis to (0, 0.3)
statistic.select(None, (0, 0.2)) # restrict second axis to (0, 0.2)
slice(*slices)

Slice statistics in place. If slice step is not 1, use rebin(). For example:

statistic.slice(slice(0, 10, 2), slice(0, 6, 3)) # rebin by factor 2 (resp. 3) along axis 0 (resp. 1), up to index 10 (resp. 6)
statistic[:10:2, :6:3] # same as above, but return new instance.
class pyrecon.metrics.MeshFFTCorrelator(mesh_reconstructed, mesh_initial, edges=None, los=None, ells=None, compensations=None)

Bases: BasePowerRatio

Estimate correlation between two meshes (reconstructed and initial fields), i.e.:

\[r(k) = \frac{P_{\mathrm{rec},\mathrm{init}}}{\sqrt{P_{\mathrm{rec}}P_{\mathrm{init}}}}\]

Initialize MeshFFTCorrelation.

Parameters:

  • mesh_reconstructed (CatalogMesh, RealField) – Mesh with reconstructed density field. If RealField, should be \(1 + \delta\) or \(\bar{n} (1 + \delta)\).

  • mesh_initial (CatalogMesh, RealField) – Mesh with initial density field (before structure formation). If RealField, should be \(1 + \delta\) or \(\bar{n} (1 + \delta)\).

  • edges (tuple, array, default=None) – \(k\)-edges. One can also provide \(\mu-edges\) (hence a tuple (kedges, muedges)). kedges may be a dictionary, with keys ‘min’ (minimum \(k\), defaults to 0), ‘max’ (maximum \(k\), defaults to np.pi/(boxsize/nmesh)), ‘step’ (if not provided pypower.fft_power.find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’).

  • los (string, array, default='firstpoint') – If los is ‘firstpoint’ (resp. ‘endpoint’), use local (varying) first point (resp. end point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • ells (tuple, list, default=None) – Multipoles to compute, which are used to compute wedges (using PowerSpectrumMultipoles.to_wedges()) when los is local (firstpoint, endpoint). In this case, if None, defaults to (0, 2, 4).

  • compensations (list, tuple, default=None) – Compensations to apply to mesh to (optionally) correct for particle-mesh assignment scheme; e.g. ‘cic-sn’ (resp. ‘cic’) for cic assignment scheme, with (resp. without) interlacing. Provide a list (or tuple) of two such strings (for mesh_reconstructed and mesh_initial`, respectively). Used only if provided mesh_reconstructed or mesh_initial are not CatalogMesh.

get_ratio(complex=False, **kwargs)

Return correlation, computed using various options.

Parameters:

  • complex (bool, default=False) – Whether (True) to return the ratio of complex power spectra, or (False) return the ratio of their real part only.

  • kwargs (dict) – Optionally, arguments for BasePowerSpectrumStatistics.get_power().

  • Results

  • -------

  • ratio (array)

modeavg(axis=0, method=None)

Return average of modes for input axis.

Parameters:

  • axis (int, default=0) – Axis.

  • method (str, default=None) – If None, return average separation from modes. If ‘mid’, return bin mid-points.

Returns:

modeavg – 1D array of size shape[axis].

Return type:

array

property ratio

Power spectrum ratio.

rebin(factor=1)

Rebin statistic, by factor(s) factor. Input factors must divide shape.

save(filename)

Save to filename.

save_txt(filename, fmt='%.12e', delimiter=' ', header=None, comments='# ', **kwargs)

Save power spectrum ratio as txt file.

Warning

Attributes are not all saved, hence there is load_txt() method.

Parameters:

  • filename (str) – File name.

  • fmt (str, default='%.12e') – Format for floating types.

  • delimiter (str, default=' ') – String or character separating columns.

  • header (str, list, default=None) – String that will be written at the beginning of the file. If multiple lines, provide a list of one-line strings.

  • comments (str, default=' #') – String that will be prepended to the header string.

  • kwargs (dict) – Arguments for get_power().

select(*xlims)

Restrict statistic to provided coordinate limits in place.

For example:

statistic.select((0, 0.3)) # restrict first axis to (0, 0.3)
statistic.select(None, (0, 0.2)) # restrict second axis to (0, 0.2)
slice(*slices)

Slice statistics in place. If slice step is not 1, use rebin(). For example:

statistic.slice(slice(0, 10, 2), slice(0, 6, 3)) # rebin by factor 2 (resp. 3) along axis 0 (resp. 1), up to index 10 (resp. 6)
statistic[:10:2, :6:3] # same as above, but return new instance.
to_propagator(growth=1.0)

Return propagator, using computed power spectra.

Parameters:

growth (float, default=1.) – Growth factor (and galaxy bias) to turn initial field to the linearly-evolved galaxy density field at the redshift of interest.

Returns:

propagator

Return type:

MeshFFTPropagator

to_transfer(growth=1.0)

Return transfer function, using computed power spectra.

Parameters:

growth (float, default=1.) – Growth factor (and galaxy bias) to turn initial field to the linearly-evolved galaxy density field at the redshift of interest.

Returns:

transfer

Return type:

MeshFFTTransfer

class pyrecon.metrics.MeshFFTPropagator(mesh_reconstructed, mesh_initial, edges=None, los=None, ells=None, compensations=None, growth=1.0)

Bases: BasePowerRatio

Estimate propagator, i.e.:

\[g(k) = \frac{P_{\mathrm{rec},\mathrm{init}}}{G(z) b(z) P_{\mathrm{init}}}\]

Initialize MeshFFTPropagator.

Parameters:

  • mesh_reconstructed (CatalogMesh, RealField) – Mesh with reconstructed density field.

  • mesh_initial (CatalogMesh, RealField) – Mesh with initial density field (before structure formation).

  • edges (tuple, array, default=None) – \(k\)-edges. One can also provide \(\mu-edges\) (hence a tuple (kedges, muedges)). kedges may be a dictionary, with keys ‘min’ (minimum \(k\), defaults to 0), ‘max’ (maximum \(k\), defaults to np.pi/(boxsize/nmesh)), ‘step’ (if not provided pypower.fft_power.find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’).

  • los (string, array, default='firstpoint') – If los is ‘firstpoint’ (resp. ‘endpoint’), use local (varying) first point (resp. end point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • ells (tuple, list, default=None) – Multipoles to compute, which are used to compute wedges (using PowerSpectrumMultipoles.to_wedges()) when los is local (firstpoint, endpoint). In this case, if None, defaults to (0, 2, 4).

  • compensations (list, tuple, default=None) – Compensations to apply to mesh to (optionally) correct for particle-mesh assignment scheme; e.g. ‘cic-sn’ (resp. ‘cic’) for cic assignment scheme, with (resp. without) interlacing. Provide a list (or tuple) of two such strings (for mesh_reconstructed and mesh_initial`, respectively). Used only if provided mesh_reconstructed or mesh_initial are not CatalogMesh.

  • growth (float, default=1.) – Growth factor (and galaxy bias) to turn initial field to the linearly-evolved galaxy density field at the redshift of interest.

get_ratio(complex=False, **kwargs)

Return propagator, computed using various options.

Parameters:

  • complex (bool, default=False) – Whether (True) to return the ratio of complex power spectra, or (False) return the ratio of their real part only.

  • kwargs (dict) – Optionally, arguments for BasePowerSpectrumStatistics.get_power().

  • Results

  • -------

  • ratio (array)

modeavg(axis=0, method=None)

Return average of modes for input axis.

Parameters:

  • axis (int, default=0) – Axis.

  • method (str, default=None) – If None, return average separation from modes. If ‘mid’, return bin mid-points.

Returns:

modeavg – 1D array of size shape[axis].

Return type:

array

property ratio

Power spectrum ratio.

rebin(factor=1)

Rebin statistic, by factor(s) factor. Input factors must divide shape.

save(filename)

Save to filename.

save_txt(filename, fmt='%.12e', delimiter=' ', header=None, comments='# ', **kwargs)

Save power spectrum ratio as txt file.

Warning

Attributes are not all saved, hence there is load_txt() method.

Parameters:

  • filename (str) – File name.

  • fmt (str, default='%.12e') – Format for floating types.

  • delimiter (str, default=' ') – String or character separating columns.

  • header (str, list, default=None) – String that will be written at the beginning of the file. If multiple lines, provide a list of one-line strings.

  • comments (str, default=' #') – String that will be prepended to the header string.

  • kwargs (dict) – Arguments for get_power().

select(*xlims)

Restrict statistic to provided coordinate limits in place.

For example:

statistic.select((0, 0.3)) # restrict first axis to (0, 0.3)
statistic.select(None, (0, 0.2)) # restrict second axis to (0, 0.2)
slice(*slices)

Slice statistics in place. If slice step is not 1, use rebin(). For example:

statistic.slice(slice(0, 10, 2), slice(0, 6, 3)) # rebin by factor 2 (resp. 3) along axis 0 (resp. 1), up to index 10 (resp. 6)
statistic[:10:2, :6:3] # same as above, but return new instance.
class pyrecon.metrics.MeshFFTTransfer(mesh_reconstructed, mesh_initial, edges=None, los=None, ells=None, compensations=None, growth=1.0)

Bases: BasePowerRatio

Estimate transfer function, i.e.:

\[t(k) = \sqrt{\frac{P_{\mathrm{rec}}}{G^{2}(z) b^{2}(z) P_{\mathrm{init}}}}\]

Initialize MeshFFTTransfer.

Parameters:

  • mesh_reconstructed (CatalogMesh, RealField) – Mesh with reconstructed density field.

  • mesh_initial (CatalogMesh, RealField) – Mesh with initial density field (before structure formation).

  • edges (tuple, array, default=None) – \(k\)-edges. One can also provide \(\mu-edges\) (hence a tuple (kedges, muedges)). kedges may be a dictionary, with keys ‘min’ (minimum \(k\), defaults to 0), ‘max’ (maximum \(k\), defaults to np.pi/(boxsize/nmesh)), ‘step’ (if not provided pypower.fft_power.find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’).

  • los (string, array, default='firstpoint') – If los is ‘firstpoint’ (resp. ‘endpoint’), use local (varying) first point (resp. end point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • ells (tuple, list, default=None) – Multipoles to compute, which are used to compute wedges (using PowerSpectrumMultipoles.to_wedges()) when los is local (firstpoint, endpoint). In this case, if None, defaults to (0, 2, 4).

  • compensations (list, tuple, default=None) – Compensations to apply to mesh to (optionally) correct for particle-mesh assignment scheme; e.g. ‘cic-sn’ (resp. ‘cic’) for cic assignment scheme, with (resp. without) interlacing. Provide a list (or tuple) of two such strings (for mesh_reconstructed and mesh_initial`, respectively). Used only if provided mesh_reconstructed or mesh_initial are not CatalogMesh.

  • growth (float, default=1.) – Growth factor (and galaxy bias) to turn initial field to the linearly-evolved galaxy density field at the redshift of interest.

get_ratio(complex=False, **kwargs)

Return transfer function, computed using various options.

Parameters:

  • complex (bool, default=False) – Whether (True) to return the ratio of complex power spectra, or (False) return the ratio of their real part only.

  • kwargs (dict) – Optionally, arguments for BasePowerSpectrumStatistics.get_power().

  • Results

  • -------

  • ratio (array)

modeavg(axis=0, method=None)

Return average of modes for input axis.

Parameters:

  • axis (int, default=0) – Axis.

  • method (str, default=None) – If None, return average separation from modes. If ‘mid’, return bin mid-points.

Returns:

modeavg – 1D array of size shape[axis].

Return type:

array

property ratio

Power spectrum ratio.

rebin(factor=1)

Rebin statistic, by factor(s) factor. Input factors must divide shape.

save(filename)

Save to filename.

save_txt(filename, fmt='%.12e', delimiter=' ', header=None, comments='# ', **kwargs)

Save power spectrum ratio as txt file.

Warning

Attributes are not all saved, hence there is load_txt() method.

Parameters:

  • filename (str) – File name.

  • fmt (str, default='%.12e') – Format for floating types.

  • delimiter (str, default=' ') – String or character separating columns.

  • header (str, list, default=None) – String that will be written at the beginning of the file. If multiple lines, provide a list of one-line strings.

  • comments (str, default=' #') – String that will be prepended to the header string.

  • kwargs (dict) – Arguments for get_power().

select(*xlims)

Restrict statistic to provided coordinate limits in place.

For example:

statistic.select((0, 0.3)) # restrict first axis to (0, 0.3)
statistic.select(None, (0, 0.2)) # restrict second axis to (0, 0.2)
slice(*slices)

Slice statistics in place. If slice step is not 1, use rebin(). For example:

statistic.slice(slice(0, 10, 2), slice(0, 6, 3)) # rebin by factor 2 (resp. 3) along axis 0 (resp. 1), up to index 10 (resp. 6)
statistic[:10:2, :6:3] # same as above, but return new instance.
class pyrecon.metrics.MetaBasePowerRatio(name, bases, class_dict)

Bases: BaseMetaClass

Metaclass adding to-wedges transforms, properties and methods to BasePowerRatio-derived classes.

mro()

Return a type’s method resolution order.

set_logger()

Add attributes for logging:

  • logger

  • methods log_debug, log_info, log_warning, log_error, log_critical

Utilities

class pyrecon.utils.BaseClass

Bases: object

Base class that implements copy(). To be used throughout this package.

class pyrecon.utils.BaseMetaClass(name, bases, class_dict)

Bases: type

Meta class to add logging attributes to BaseClass derived classes.

mro()

Return a type’s method resolution order.

set_logger()

Add attributes for logging:

  • logger

  • methods log_debug, log_info, log_warning, log_error, log_critical

class pyrecon.utils.DistanceToRedshift(distance, zmax=100.0, nz=2048, interp_order=3)

Bases: object

Class that holds a conversion distance -> redshift.

Initialize DistanceToRedshift. Creates an array of redshift -> distance in log(redshift) and instantiates a spline interpolator distance -> redshift.

Parameters:

  • distance (callable) – Callable that provides distance as a function of redshift (array).

  • zmax (float, default=100.) – Maximum redshift for redshift <-> distance mapping.

  • nz (int, default=2048) – Number of points for redshift <-> distance mapping.

  • interp_order (int, default=3) – Interpolation order, e.g. 1 for linear interpolation, 3 for cubic splines.

class pyrecon.utils.MemoryMonitor(pid=None)

Bases: object

Class that monitors memory usage and clock, useful to check for memory leaks.

>>> with MemoryMonitor() as mem:
        '''do something'''
        mem()
        '''do something else'''

Initalize MemoryMonitor and register current memory usage.

Parameters:

pid (int, default=None) – Process identifier. If None, use the identifier of the current process.

pyrecon.utils.cartesian_to_sky(position, wrap=True, degree=True)

Transform Cartesian coordinates into distance, RA, Dec.

Parameters:

  • position (array of shape (N, 3)) – Position in Cartesian coordinates.

  • wrap (bool, default=True) – Whether to wrap RA in \([0, 2 \pi]\).

  • degree (bool, default=True) – Whether RA, Dec are in degrees (True) or radians (False).

Returns:

  • dist (array) – Distance.

  • ra (array) – Right Ascension.

  • dec (array) – Declination.

pyrecon.utils.distance(position)

Return cartesian distance, taking coordinates along position last axis.

pyrecon.utils.exception_handler(exc_type, exc_value, exc_traceback)

Print exception with a logger.

pyrecon.utils.mkdir(dirname)

Try to create dirname and catch OSError.

pyrecon.utils.safe_divide(x, y, inplace=False)

Divide x by y after replacing 0 in y by 1. If inplace is True, x and y modified in-place.

pyrecon.utils.setup_logging(level=20, stream=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>, filename=None, filemode='w', **kwargs)

Set up logging.

Parameters:

  • level (string, int, default=logging.INFO) – Logging level.

  • stream (_io.TextIOWrapper, default=sys.stdout) – Where to stream.

  • filename (string, default=None) – If not None stream to file name.

  • filemode (string, default='w') – Mode to open file, only used if filename is not None.

  • kwargs (dict) – Other arguments for logging.basicConfig().

pyrecon.utils.sky_to_cartesian(dist, ra, dec, degree=True, dtype=None)

Transform distance, RA, Dec into Cartesian coordinates.

Parameters:

  • dist (array of shape (N,)) – Distance.

  • ra (array of shape (N,)) – Right Ascension.

  • dec (array of shape (N,)) – Declination.

  • degree (default=True) – Whether RA, Dec are in degrees (True) or radians (False).

  • dtype (numpy.dtype, default=None) – numpy.dtype for returned array.

Returns:

position – Position in Cartesian coordinates.

Return type:

array of shape (N, 3)