API

FFT-based power spectrum estimator

Implementation of power spectrum estimator, following https://arxiv.org/abs/1704.02357. Apart from interface choices, differences w.r.t. original nbodykit’s implementation https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/convpower/fkp.py are:

  • real space positions are taken at mesh nodes, instead of 0.5 cell shift (matters only for ell > 0 in global line-of-sight)

  • normalization is computed with density obtained by paintaing data/randoms to mesh, instead of relying on \(\bar{n}_{i}\) column in the catalogs

  • FKP weights are treated as other weights

class pypower.fft_power.BasePowerSpectrumStatistics(edges, modes, power_nonorm, nmodes, wnorm=1.0, shotnoise_nonorm=0.0, power_zero_nonorm=None, power_direct_nonorm=None, attrs=None, mpicomm=None)

Bases: pypower.utils.BaseClass

Base template power statistic class. Specific power statistic should extend this class.

We recommend accessing power spectrum measurements through get_power(), or __call__() (accessed through my_power_statistic_instance()).

Initialize BasePowerSpectrumStatistics.

Parameters

  • edges (tuple of ndim arrays) – Edges used to bin power spectrum measurement.

  • modes (array) – Mean “wavevector” (e.g. \((k, \mu)\)) in each bin.

  • power_nonorm (array) – Power spectrum in each bin, without normalization.

  • nmodes (array) – Number of modes in each bin.

  • wnorm (float, default=1.) – Power spectrum normalization.

  • shotnoise_nonorm (float, default=0.) – Shot noise, without normalization.

  • power_zero_nonorm (float, array, default=0.) – Value of power spectrum at \(k = 0\).

  • power_direct_nonorm (array, default=0.) – Value of pair-count-based ‘direct’ power spectrum estimation, (e.g. PIP and angular upweights correction, eq. 26 of https://arxiv.org/abs/1912.08803), to be added to power_nonorm.

  • attrs (dict, default=None) – Dictionary of other attributes.

  • mpicomm (MPI communicator, default=None) – The MPI communicator, only used when saving (save() and save_txt()) statistics.

get_power(add_direct=True, remove_shotnoise=True, null_zero_mode=True, divide_wnorm=True, complex=True)

Return power spectrum, computed using various options.

Parameters

  • add_direct (bool, default=True) – Add pair-count-based ‘direct’ power spectrum measurement.

  • remove_shotnoise (bool, default=True) – Remove estimated shot noise.

  • null_zero_mode (bool, default=True) – Remove power spectrum at \(k = 0\) (if within edges).

  • divide_wnorm (bool, default=True) – Divide by estimated power spectrum normalization.

  • complex (bool, default=True) – Whether (True) to return the complex power spectrum, or (False) return its real part only.

  • Results –

  • ------- –

  • power (array) –

property k

Wavenumbers.

property kedges

Wavenumber edges.

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 ndim

Return binning dimensionality.

property power

Power spectrum, normalized and with shot noise removed.

rebin(factor=1)

Rebin power spectrum estimation in place, by factor(s) factor. A tuple must be provided in case ndim is greater than 1. Input factors must divide shape.

save(filename)

Save to filename.

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

Save power spectrum 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)
property shape

Return shape of binned power spectrum power.

property shotnoise

Normalized shot noise.

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.
classmethod sum(*others)

Sum input power spectra, weighted by their wnorm.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

property with_mpi

Whether to use MPI.

class pypower.fft_power.CatalogFFTPower(data_positions1, data_positions2=None, randoms_positions1=None, randoms_positions2=None, shifted_positions1=None, shifted_positions2=None, data_weights1=None, data_weights2=None, randoms_weights1=None, randoms_weights2=None, shifted_weights1=None, shifted_weights2=None, D1D2_twopoint_weights=None, D1R2_twopoint_weights=None, R1D2_twopoint_weights=None, D1S2_twopoint_weights=None, S1D2_twopoint_weights=None, edges=None, ells=(0, 2, 4), los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, dtype='f8', resampler='tsc', interlacing=2, position_type='xyz', weight_type='auto', weight_attrs=None, direct_engine='corrfunc', direct_limits=(0.0, 0.03333333333333333), direct_limit_type='degree', wnorm=None, shotnoise=None, mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Bases: pypower.fft_power.MeshFFTPower

Wrapper on MeshFFTPower to estimate power spectrum directly from positions and weights.

Initialize CatalogFFTPower, i.e. estimate power spectrum.

Warning

In case line-of-sight is not local, one can provide \(\mu\)-edges. In this case, integration over Legendre polynomials for multipoles is performed between the first and last \(\mu\)-edges. For example, with \(\mu\)-edges [0.2, 0.4, 0.8], integration is performed between \(\mu = 0.2\) and \(\mu = 0.8\). In all other cases, integration is performed between \(\mu = -1.0\) and \(\mu = 1.0\).

Note

When running with MPI, input positions and weights are assumed to be scatted on all MPI ranks of mpicomm. If this is not the case, use mpi.scatter_array().

Parameters

  • data_positions1 (list, array) – Positions in the first data catalog. Typically of shape (3, N) or (N, 3).

  • data_positions2 (list, array, default=None) – Optionally (for cross-correlation), positions in the second data catalog. See data_positions1.

  • randoms_positions1 (list, array, default=None) – Optionally, positions of the random catalog representing the first selection function. If no randoms are provided, selection function will be assumed uniform.

  • randoms_positions2 (list, array, default=None) – Optionally (for cross-correlation), positions in the second randoms catalog. See randoms_positions1.

  • shifted_positions1 (array, default=None) – Optionally, in case of BAO reconstruction, positions of the first shifted catalog.

  • shifted_positions2 (array, default=None) – Optionally, in case of BAO reconstruction, positions of the second shifted catalog.

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

  • data_weights2 (array of shape (N,), default=None) – Optionally (for cross-correlation), weights in the second data catalog.

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

  • randoms_weights2 (array of shape (N,), default=None) – Optionally (for cross-correlation), weights in the second randoms catalog.

  • shifted_weights1 (array, default=None) – Optionally, weights of the first shifted catalog. See data_weights1.

  • shifted_weights2 (array, default=None) – Optionally, weights of the second shifted catalog. See shifted_weights1.

  • edges (tuple, array, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu\)-edges (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). 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 find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’). For both \(k\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. edges[i] <= x < edges[i+1]. However, last \(\mu\)-bin is inclusive on both ends: edges[-2] <= mu <= edges[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\).

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • 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.

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

  • boxsize (float, default=None) – Physical size of the box, 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.

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

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

  • interlacing (bool, int, default=2) – Whether to use interlacing to reduce aliasing when painting the particles on the mesh. If positive int, the interlacing order (minimum: 2).

  • 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)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”inverse_bitwise”: each pair is weighted by \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\).

      Multiple bitwise weights can be provided as a list. Individual weights can additionally be provided as float arrays. In case of cross-correlations with floating weights, bitwise weights are automatically turned to IIP weights, i.e. \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1}))\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, “inverse_bitwise” if one of input weights is integer, else “product_individual”.

    In addition, angular upweights can be provided with D1D2_twopoint_weights, D1R2_twopoint_weights, etc. If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). Inverse probability weight is then computed as: \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0.

  • D1D2_twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles between first and second data catalogs. A WeightTwoPointEstimator instance (from pycorr) or any object with arrays sep (separations) and weight (weight at given separation) as attributes (i.e. to be accessed through twopoint_weights.sep, twopoint_weights.weight) or as keys (i.e. twopoint_weights['sep'], twopoint_weights['weight']) or as element (i.e. sep, weight = twopoint_weights).

  • D1R2_twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles between first data catalog and second randoms catalog. See D1D2_twopoint_weights.

  • R1D2_twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles between second data catalog and first randoms catalog. See D1D2_twopoint_weights.

  • D1S2_twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles between first data catalog and second shifted catalog. See D1D2_twopoint_weights.

  • S1D2_twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles between second data catalog and first shifted catalog. See D1D2_twopoint_weights.

  • direct_engine (string, default='corrfunc') – Engine for direct power spectrum computation (if input weights are bitwise weights), one of [“kdtree”, “corrfunc”].

  • direct_limits (tuple, default=(0., 2./60.)) – Limits of particle pair separations used in the direct power spectrum computation.

  • direct_limit_type (string, default='degree') – Type of direct_limits; i.e. are those angular limits (“degree”, “radian”), or 3D limits (“s”)?

  • wnorm (float, default=None) – Power spectrum normalization, to use instead of internal estimate obtained with normalization().

  • shotnoise (float, default=None) – Power spectrum shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation and in case of auto-correlation is obtained by dividing CatalogMesh.unnormalized_shotnoise() by power spectrum normalization.

  • 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.

property boxsize

Physical box size.

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_power.MeshFFTPower(mesh1, mesh2=None, edges=None, ells=(0, 2, 4), los=None, boxcenter=None, compensations=None, wnorm=None, shotnoise=None)

Bases: pypower.utils.BaseClass

Class that computes power spectrum from input mesh(es), using global or local line-of-sight, following https://arxiv.org/abs/1704.02357. In effect, this class merges nbodykit’s implementation of the global line-of-sight (periodic) algorithm of: https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/fftpower.py with the local line-of-sight algorithm of: https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/convpower/fkp.py

poles

Estimated power spectrum multipoles.

Type

PowerSpectrumMultipoles

wedges

Estimated power spectrum wedges (if relevant).

Type

PowerSpectrumWedges

Initialize MeshFFTPower, i.e. estimate power spectrum.

Warning

In case line-of-sight is not local, one can provide \(\mu\)-edges. In this case, integration over Legendre polynomials for multipoles is performed between the first and last \(\mu\)-edges. For example, with \(\mu\)-edges [0.2, 0.4, 0.8], integration is performed between \(\mu = 0.2\) and \(\mu = 0.8\). In all other cases, integration is performed between \(\mu = -1.0\) and \(\mu = 1.0\).

Parameters

  • mesh1 (CatalogMesh, RealField, ComplexField) – First mesh. If RealField, assumed to be \(1 + \delta\) or \(\bar{n} (1 + \delta)\). In case of ComplexField, assumed to be the FFT of \(\delta\) (or \(1 + \delta\)), i.e. unit density.

  • mesh2 (CatalogMesh, RealField, ComplexField, default=None) – In case of cross-correlation, second mesh, with same size and physical extent (boxsize and boxcenter) that mesh1.

  • edges (tuple, array, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu\)-edges (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). 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 find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’). For both \(k\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. edges[i] <= x < edges[i+1]. However, last \(\mu\)-bin is inclusive on both ends: edges[-2] <= mu <= edges[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\).

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • los (string, array, default=None) – 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.

  • boxcenter (float, array, default=None) – Box center; defaults to 0. Used only if provided mesh1 and mesh2 are not CatalogMesh.

  • compensations (list, tuple, string, default=None) – Compensations to apply to mesh to (optionally) correct for particle-mesh assignment scheme; e.g. ‘cic’ (resp. ‘cic-sn’) for cic assignment scheme, with (resp. without) interlacing. In case mesh2 is not None (cross-correlation), provide a list (or tuple) of two such strings (for mesh1 and mesh2, respectively). Used only if provided mesh1 or mesh2 are not CatalogMesh.

  • wnorm (float, default=None) – Power spectrum normalization, to use instead of internal estimate obtained with normalization().

  • shotnoise (float, default=None) – Power spectrum shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation or both mesh1 and mesh2 are pmesh.pm.RealField, and in case of auto-correlation is obtained by dividing CatalogMesh.unnormalized_shotnoise() of mesh1 by power spectrum normalization.

property boxsize

Physical box size.

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_power.MetaPowerSpectrumStatistics(name, bases, class_dict)

Bases: pypower.utils.BaseMetaClass

Metaclass to return correct power spectrum statistic.

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 pypower.fft_power.PowerSpectrumMultipoles(edges, modes, power_nonorm, nmodes, ells, **kwargs)

Bases: pypower.fft_power.BasePowerSpectrumStatistics

Power spectrum multipoles binned in \(k\).

Initialize PowerSpectrumMultipoles.

Parameters

  • edges (tuple of ndim arrays) – Edges used to bin power spectrum measurement.

  • modes (array) – Mean “wavevector” (e.g. \((k, \mu)\)) in each bin.

  • power_nonorm (array) – Power spectrum in each bin, without normalization.

  • nmodes (array) – Number of modes in each bin.

  • ells (tuple, list.) – Multipole orders.

  • kwargs (dict) – Other arguments for BasePowerSpectrumStatistics.

get_power(add_direct=True, remove_shotnoise=True, null_zero_mode=True, divide_wnorm=True, complex=True)

Return power spectrum, computed using various options.

Parameters

  • add_direct (bool, default=True) – Add direct power spectrum measurement.

  • remove_shotnoise (bool, default=True) – Remove estimated shot noise.

  • null_zero_mode (bool, default=True) – Remove power spectrum at \(k = 0\) (if within edges).

  • divide_wnorm (bool, default=True) – Divide by estimated power spectrum normalization.

  • complex (bool, default=True) – Whether (True) to return the complex power spectrum, or (False) return its real part if even multipoles, imaginary part if odd multipole.

  • Results –

  • ------- –

  • power (array) –

property k

Wavenumbers.

property kavg

Mode-weighted average wavenumber = k.

property kedges

Wavenumber edges.

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 ndim

Return binning dimensionality.

property power

Power spectrum, normalized and with shot noise removed.

rebin(factor=1)

Rebin power spectrum estimation in place, by factor(s) factor. A tuple must be provided in case ndim is greater than 1. Input factors must divide shape.

save(filename)

Save to filename.

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

Save power spectrum 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)
property shape

Return shape of binned power spectrum power.

property shotnoise

Normalized shot noise.

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.
classmethod sum(*others)

Sum input power spectra, weighted by their wnorm.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

property with_mpi

Whether to use MPI.

class pypower.fft_power.PowerSpectrumStatistics(*args, statistic='wedge', **kwargs)

Bases: pypower.utils.BaseClass

Entry point to power spectrum statistics.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_power.PowerSpectrumWedges(edges, modes, power_nonorm, nmodes, wnorm=1.0, shotnoise_nonorm=0.0, power_zero_nonorm=None, power_direct_nonorm=None, attrs=None, mpicomm=None)

Bases: pypower.fft_power.BasePowerSpectrumStatistics

Power spectrum binned in \((k, \mu)\).

Initialize BasePowerSpectrumStatistics.

Parameters

  • edges (tuple of ndim arrays) – Edges used to bin power spectrum measurement.

  • modes (array) – Mean “wavevector” (e.g. \((k, \mu)\)) in each bin.

  • power_nonorm (array) – Power spectrum in each bin, without normalization.

  • nmodes (array) – Number of modes in each bin.

  • wnorm (float, default=1.) – Power spectrum normalization.

  • shotnoise_nonorm (float, default=0.) – Shot noise, without normalization.

  • power_zero_nonorm (float, array, default=0.) – Value of power spectrum at \(k = 0\).

  • power_direct_nonorm (array, default=0.) – Value of pair-count-based ‘direct’ power spectrum estimation, (e.g. PIP and angular upweights correction, eq. 26 of https://arxiv.org/abs/1912.08803), to be added to power_nonorm.

  • attrs (dict, default=None) – Dictionary of other attributes.

  • mpicomm (MPI communicator, default=None) – The MPI communicator, only used when saving (save() and save_txt()) statistics.

get_power(add_direct=True, remove_shotnoise=True, null_zero_mode=True, divide_wnorm=True, complex=True)

Return power spectrum, computed using various options.

Parameters

  • add_direct (bool, default=True) – Add pair-count-based ‘direct’ power spectrum measurement.

  • remove_shotnoise (bool, default=True) – Remove estimated shot noise.

  • null_zero_mode (bool, default=True) – Remove power spectrum at \(k = 0\) (if within edges).

  • divide_wnorm (bool, default=True) – Divide by estimated power spectrum normalization.

  • complex (bool, default=True) – Whether (True) to return the complex power spectrum, or (False) return its real part only.

  • Results –

  • ------- –

  • power (array) –

property k

Wavenumbers.

property kavg

Mode-weighted average wavenumber.

property kedges

Wavenumber edges.

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 mu

Cosine angle to line-of-sight.

property muavg

Mode-weighted average \(\mu\).

property muedges

\(\mu\)-edges.

property ndim

Return binning dimensionality.

property power

Power spectrum, normalized and with shot noise removed.

rebin(factor=1)

Rebin power spectrum estimation in place, by factor(s) factor. A tuple must be provided in case ndim is greater than 1. Input factors must divide shape.

save(filename)

Save to filename.

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

Save power spectrum 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)
property shape

Return shape of binned power spectrum power.

property shotnoise

Normalized shot noise.

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.
classmethod sum(*others)

Sum input power spectra, weighted by their wnorm.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

property with_mpi

Whether to use MPI.

pypower.fft_power.find_unique_edges(x, x0, xmin=0.0, xmax=inf, mpicomm=mpi4py.MPI.COMM_WORLD)

Construct unique edges for distribution of Cartesian distances corresponding to coordinates x. Taken from https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/fftpower.py.

Parameters

  • x (list of ndim arrays) – List of ndim (broadcastable) coordinate arrays.

  • x0 (array_like of shape (ndim, )) – 3-vector of fundamental coordinate separation.

  • xmin (float, default=0.) – Minimum separation.

  • xmax (float, default=np.inf) – Maximum separation.

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

Returns

edges – Edges, starting at 0, such that each bin contains a unique value of Cartesian distances.

Return type

array

pypower.fft_power.get_power_statistic(statistic='wedge')

Return BasePowerSpectrumStatistics subclass corresponding to statistic (either ‘wedge’ or ‘multipole’).

pypower.fft_power.get_real_Ylm(ell, m, modules=None)

Return a function that computes the real spherical harmonic of order (ell, m). Adapted from https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/convpower/fkp.py.

Note

Faster evaluation will be achieved if sympy and numexpr are available. Else, fallback to numpy and scipy’s functions.

Parameters

  • ell (int) – The degree of the harmonic.

  • m (int) – The order of the harmonic; abs(m) <= ell.

  • modules (str, default=None) – If ‘sympy’, use sympy + numexpr to speed up calculation. If ‘scipy’, use scipy. If None, defaults to sympy if installed, else scipy.

Returns

Ylm – A function that takes 3 arguments: (xhat, yhat, zhat) unit-normalized Cartesian coordinates and returns the specified Ylm.

Return type

callable

References

https://en.wikipedia.org/wiki/Spherical_harmonics#Real_form

pypower.fft_power.normalization(mesh1, mesh2=None, uniform=False, resampler='cic', cellsize=10.0, fields=None)

Return DESI-like normalization, summing over mesh cells:

\[A = 1/dV \frac{\alpha_{2} \sum_{i} n_{d,1}^{i} n_{r,2}^{i} + \alpha_{1} \sum_{i} n_{d,2}^{i} n_{r,1}^{i}}{2}\]

\(n_{d,1}^{i}\) and \(n_{r,1}^{i}\) are the first data and randoms density meshes, as obtained by painting data \(w_{d}\) and random weights \(w_{r}\) on the same mesh (of cell volume \(dV\)), using the cic assignment scheme. The sum then runs over the mesh cells. \(\alpha_{1} = \sum_{i} w_{d,1}^{i} / \sum_{i} w_{r,1}^{i}\) and \(\alpha_{2} = \sum_{i} w_{d,2}^{i} / \sum_{i} w_{r,2}^{i}\) where the sum of weights is performed over the catalogs. If no randoms are provided, density is supposed to be uniform and mesh1 and mesh2 are assumed to occupy the same physical volume.

Parameters

  • mesh1 (CatalogMesh, RealField, ComplexField) – First mesh. If RealField, density is assumed to be uniform, mesh1.csum()/np.prod(mesh1.pm.BoxSize). If ComplexField, assumed to be the FFT of \(\delta\) (or \(1 + \delta\)), i.e. unit density.

  • mesh2 (CatalogMesh, RealField, ComplexField, default=None) – Second mesh, for cross-correlations.

  • uniform (bool, default=False) – Whether to assume uniform selection function (only revelant when both mesh1 and mesh2 are CatalogMesh).

  • resampler (string, ResampleWindow, default='cic') – Particle-mesh assignment scheme. Choices are [‘ngp’, ‘cic’, ‘tsc’, ‘pcs’].

  • cellsize (array, float) – Physical size of mesh cells used to paint mesh1 and mesh2 (if instance of CatalogMesh).

Returns

norm – Normalization.

Return type

float

pypower.fft_power.normalization_from_nbar(nbar, weights=None, data_weights=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Return BOSS/eBOSS-like normalization, summing over \(\bar{n}_{i}\) and weight columns, i.e.:

\[\alpha \sum_{i=0}^{N} w_{i} \bar{n}_{i}\]
Parameters

  • nbar (array of shape (N,)) – \(\bar{n}_{i}\) (comoving density) column.

  • weights (array of shape (N,), default=None) – Weights, if any.

  • data_weights (array of shape (N,), default=None) – Data weights, to normalize randoms weights by alpha = sum(data_weights)/sum(weights).

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

Returns

norm – Normalization.

Return type

float

pypower.fft_power.project_to_basis(y3d, edges, los=(0, 0, 1), ells=None, antisymmetric=False, exclude_zero=False)

Project a 3D statistic on to the specified basis. The basis will be one of:

  • 2D \((x, \mu)\) bins: \(\mu\) is the cosine of the angle to the line-of-sight

  • 2D \((x, \ell)\) bins: \(\ell\) is the multipole number, which specifies the Legendre polynomial when weighting different \(\mu\) bins.

Adapted from https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/fftpower.py.

Notes

In single precision (float32/complex64) nbodykit’s implementation is fairly imprecise due to incorrect binning of \(x\) and \(\mu\) modes. Here we cast mesh coordinates to the maximum precision of input edges, which makes computation much more accurate in single precision.

Notes

We deliberately set to 0 the number of modes beyond Nyquist, as it is unclear whether to count Nyquist as \(\mu\) or \(-\mu\) (it should probably be half weight for both). Our safe choice ensures consistent results between hermitian compressed and their associated uncompressed fields.

Notes

The 2D \((x, \ell)\) bins will be computed only if ells is specified. See return types for further details. For both \(x\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. mode mode falls in bin i if edges[i] <= mode < edges[i+1]. However, last \(\mu\)-bin is inclusive on both ends: edges[-2] <= mu <= edges[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\).

Warning

Integration over Legendre polynomials for multipoles is performed between the first and last \(\mu\)-edges, e.g. with \(\mu\)-edges [0.2, 0.4, 0.8], integration is performed between \(\mu = 0.2\) and \(\mu = 0.8\).

Parameters

  • y3d (RealField or ComplexField) – The 3D array holding the statistic to be projected to the specified basis.

  • edges (list of arrays, (2,)) – List of arrays specifying the edges of the desired \(x\) bins and \(\mu\) bins; assumed sorted.

  • los (array_like, default=(0, 0, 1)) – The line-of-sight direction to use, which \(\mu\) is defined with respect to.

  • ells (tuple of ints, default=None) – If provided, a list of integers specifying multipole numbers to project the 2D \((x, \mu)\) bins on to.

Returns

  • result (tuple) –

    The 2D binned results; a tuple of (xmean2d, mumean2d, y2d, n2d), where:

    • xmean2darray_like, (nx, nmu)

      The mean \(x\) value in each 2D bin

    • mumean2darray_like, (nx, nmu)

      The mean \(\mu\) value in each 2D bin

    • y2darray_like, (nx, nmu)

      The mean y3d value in each 2D bin

    • n2darray_like, (nx, nmu)

      The number of values averaged in each 2D bin

  • result_poles (tuple or None) – The multipole results; if ells supplied it is a tuple of (xmean1d, poles, n1d), where:

    • xmean1darray_like, (nx,)

      The mean \(x\) value in each 1D multipole bin

    • polesarray_like, (nell, nx)

      The mean multipoles value in each 1D bin

    • n1darray_like, (nx,)

      The number of values averaged in each 1D bin

Painting catalog on mesh

Implementation of methods to paint a catalog on mesh; workhorse is CatalogMesh.

pypower.mesh.ArrayMesh(array, boxsize, mpiroot=0, mpicomm=mpi4py.MPI.COMM_WORLD)

Turn numpy array into pmesh.pm.RealField.

Parameters

  • array (array) – Mesh numpy array gathered on mpiroot.

  • boxsize (array) – Physical box size.

  • mpiroot (int, default=0) – MPI rank where input array is gathered.

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

Returns

mesh

Return type

pmesh.pm.RealField

class pypower.mesh.CatalogMesh(data_positions, data_weights=None, randoms_positions=None, randoms_weights=None, shifted_positions=None, shifted_weights=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, dtype='f8', resampler='tsc', interlacing=2, position_type='xyz', copy=False, mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Bases: pypower.utils.BaseClass

Class to paint catalog of positions and weights to mesh.

Initialize CatalogMesh.

Note

When running with MPI, input positions and weights are assumed to be scatted on all MPI ranks of mpicomm. If this is not the case, use mpi.scatter_array().

Parameters

  • data_positions (list, array) – Positions in the data catalog. Typically of shape (3, N) or (N, 3).

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

  • randoms_positions (list, array) – Positions in the randoms catalog. Typically of shape (3, N) or (N, 3).

  • randoms_weights (array of shape (N,), default=None) – Randoms weights.

  • shifted_positions (array, default=None) – Optionally, in case of BAO reconstruction, positions of the shifted catalog.

  • shifted_weights (array, default=None) – Optionally, in case of BAO reconstruction, weigths of the shifted catalog.

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

  • boxsize (float, default=None) – Physical size of the box. If not provided, see positions.

  • boxcenter (array, float, default=None) – Box center. If not provided, see 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.

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

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

  • dtype (string, dtype, default='f8') – The data type to use for the mesh. Input positions and weights are cast to the corresponding (real) precision.

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

  • interlacing (bool, int, default=2) – Whether to use interlacing to reduce aliasing when painting the particles on the mesh. If positive int, the interlacing order (minimum: 2).

  • 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)

  • copy (bool, default=False) – If False, avoids copy of positions and weights if they are of (real) type dtype, mpiroot is None, and position_type is “pos” (for positions). Setting to True is only useful if one wants to modify positions or weights that have been passed as input while keeping those attached to the current mesh instance the same.

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

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

clone(data_positions=None, data_weights=None, randoms_positions=None, randoms_weights=None, shifted_positions=None, shifted_weights=None, boxsize=None, cellsize=None, nmesh=None, boxcenter=None, dtype=None, resampler=None, interlacing=None, position_type='xyz', mpicomm=None)

Clone current instance, i.e. copy and set new positions and weights. Arguments ‘boxsize’, ‘nmesh’, ‘boxcenter’, ‘dtype’, ‘resampler’, ‘interlacing’, ‘mpicomm’, if None, are overriden by those of the current instance.

property compensation

Return dictionary specifying compensation scheme for particle-mesh resampling.

save(filename)

Save to filename.

to_mesh(field=None, dtype=None, compensate=False)

Paint positions/weights to mesh.

Parameters

  • field (string, default=None) –

    Field to paint to mesh, one of:

    • ”data”: data positions and weights

    • ”shifted”: shifted positions and weights (available only if shifted positions are provided)

    • ”randoms”: randoms positions and weights

    • ”data-normalized_shifted”: shifted positions and weights, renormalized (by alpha)

      such that their sum is same as data weights

    • ”data-normalized_randoms”: randoms positions and weights, renormalized (by alpha)

      such that their sum is same as data weights

    • ”fkp”: FKP field, i.e. data - alpha * (shifted if provided else randoms)

    • None: defaults to “data” if no shifted/randoms, else “fkp”

  • dtype (string, dtype, default='f8') – The data type of the mesh when painting, to override current dtype.

  • compensate (bool, default=False) – Wether to apply compensation for particle-mesh assignment scheme.

Returns

out – Mesh, with values in “weights” units (not normalized as density).

Return type

RealField

unnormalized_shotnoise()

Return unnormalized shotnoise, as:

\[\sum_{i=1}^{N_{g}} w_{i,g}^{2} + \alpha^{2} \sum_{i=1}^{N_{r}} w_{i,r}^{2}\]

Where the sum runs over data (and optionally) shifted/randoms weights.

property with_mpi

Whether to use MPI.

property with_randoms

Whether randoms positions have been provided.

property with_shifted

Whether “shifted” positions have been provided (e.g. for reconstruction).

Approximate window

Implementation of (approximate) window function estimation and convolution. Typically, the window function will be estimated through CatalogSmoothWindow, and window function matrices using PowerSpectrumSmoothWindowMatrix, following https://arxiv.org/abs/2106.06324.

class pypower.smooth_window.CatalogSmoothWindow(randoms_positions1=None, randoms_positions2=None, randoms_weights1=None, randoms_weights2=None, edges=None, projs=None, power_ref=None, los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, dtype=None, resampler=None, interlacing=None, position_type='xyz', weight_type='auto', weight_attrs=None, wnorm=None, shotnoise=None, mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Bases: pypower.fft_power.MeshFFTPower

Wrapper on MeshFFTPower to estimate window function from input random positions and weigths.

Initialize CatalogSmoothWindow, i.e. estimate power spectrum window.

Parameters

  • randoms_positions1 (list, array, default=None) – Positions in the first randoms catalog. Typically of shape (3, N) or (N, 3).

  • randoms_positions2 (list, array, default=None) – Optionally (for cross-correlation), positions in the second randoms catalog. See randoms_positions1.

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

  • randoms_weights2 (array of shape (N,), default=None) – Optionally (for cross-correlation), weights in the second randoms catalog.

  • edges (tuple, array, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu-edges\) (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). 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 find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’).

  • projs (list, default=None) – List of Projection instances or (multipole, wide-angle order) tuples. If None, and power_ref is provided, the list of projections is set to be able to compute window convolution of theory power spectrum multipoles of orders power_ref.ells. Namely, for maximum theory and output multipole \(\ell_{\mathrm{max}}\), window function multipoles will be computed at wide-angle order 0 up to \(2 \ell_{\mathrm{max}}\) (maximum order yielded by the product of any theory and output Legendre polynomial, see e.g. eq. C8 of https://arxiv.org/pdf/2106.06324.pdf). In addition, if chosen line-of-sight is local (either ‘firstpoint’ or ‘endpoint’), odd poles of the window function will be computed at wide-angle order 1, up to \(2 \ell_{\mathrm{max}} + 1\) (maximum non-zero odd pole of wide angle order 1 generated by even poles up to \(\ell_{\mathrm{max}}\) is \(\ell_{\mathrm{max}} + 1\)). Finally, if any of power_ref.ells is odd, all (even and odd) poles will be computed at wide-angle orders 0 and 1, up to \(2 \ell_{\mathrm{max}}\) and \(2 \ell_{\mathrm{max}} + 1\), respectively.

  • power_ref (PowerSpectrumMultipoless, default=None) – “Reference” power spectrum estimation, e.g. of the actual data. It is used to set default values for projs, los, boxsize, boxcenter, nmesh, interlacing, resampler and wnorm if those are None.

  • los (string, array, default=None) – 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. If None, defaults to line-of-sight used in estimation of power_ref.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis. If None, defaults to the value used in estimation of power_ref.

  • boxsize (float, default=None) – Physical size of the box, defaults to maximum extent taken by all input positions, times boxpad. If None, defaults to the value used in estimation of power_ref.

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

  • 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.

  • dtype (string, dtype, default=None) – The data type to use for input positions and weights and the mesh. If None, defaults to the value used in estimation of power_ref if provided, else ‘f8’.

  • resampler (string, ResampleWindow, default=None) – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’]. If None, defaults to the value used in estimation of power_ref.

  • interlacing (bool, int, default=None) – Whether to use interlacing to reduce aliasing when painting the particles on the mesh. If positive int, the interlacing order (minimum: 2). If None, defaults to the value used in estimation of power_ref.

  • 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)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, else “product_individual”.

    If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of weights if the denominator is zero (defaulting to 0).

  • wnorm (float, default=None) – Window function normalization. If None, defaults to the value used in estimation of power_ref, rescaled to the input random weights — which yields a correct normalization of the window function for the power spectrum estimation power_ref. If power_ref provided, use internal estimate obtained with normalization() — which is wrong (the normalization poles.wnorm can be reset a posteriori using the above recipe).

  • shotnoise (float, default=None) – Window function shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation and in case of auto-correlation is obtained by dividing CatalogMesh.unnormalized_shotnoise() by window function normalization.

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

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

property boxsize

Physical box size.

classmethod concatenate(*others, **kwargs)

Concatenate poles. Same argument as PowerSpectrumSmoothWindow.concatenate_x().

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.smooth_window.CorrelationFunctionSmoothWindow(sep, corr, projs)

Bases: pypower.utils.BaseClass

Correlation window function multipoles.

Initialize CorrelationFunctionSmoothWindow.

Parameters

  • modes (array) – Mean separation.

  • corr (array) – Mean correlation.

  • projs (list) – List of Projection instances or (multipole, wide-angle order) tuples.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.smooth_window.CorrelationFunctionSmoothWindowMatrix(sep, projsin, projsout=None, window=None, sum_wa=True, default_zero=False, attrs=None)

Bases: pypower.wide_angle.BaseMatrix

Class computing matrix for window product in configuration space.

projmatrix

Array of shape (len(self.projsout), len(self.projsin), len(self.x)).

Type

array

Initialize CorrelationFunctionSmoothWindowMatrix.

Parameters

  • sep (array) – Input (and ouput) separations.

  • projsin (list) – Input projections.

  • projsout (list, default=None) – Output projections. Defaults to propose_out(projsin, sum_wa=sum_wa).

  • window (CorrelationFunctionSmoothWindow, PowerSpectrumSmoothWindow) – Window function to convolve power spectrum with. If a PowerSpectrumSmoothWindow instance is provided, it is transformed to configuration space.

  • sum_wa (bool, default=True) – Whether to perform summation over output wide-angle orders. Always set to True except for debugging purposes.

  • default_zero (bool, default=False) – If a given projection is not provided in window function, set to 0. Else an IndexError is raised.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

static propose_out(projsin, sum_wa=True)

Propose output projections given proposed input projections projsin. If sum_wa is True (typically always the case), return projections with Projection.wa_order set to None (all wide-angle orders have been summed).

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

resum_input_odd_wide_angle(**kwargs)

Resum odd wide-angle orders. By default, line-of-sight is chosen as that save in attrs (attrs['los_type']). To override, use input kwargs which will be passed to CorrelationFunctionOddWideAngleMatrix.

run()

Set up transform, i.e. compute matrix:

\[W_{\ell,\ell^{\prime}}^{(n,n^{\prime})}(s) = \delta_{n n^{\prime}} \sum_{L} C_{\ell \ell^{\prime} L} Q_{L}^{(n)}(s)\]

with \(\ell\) multipole order corresponding to projout.ell and \(\ell^{\prime}\) to projin.ell, \(n\) wide angle order corresponding to projout.wa_order and \(n^{\prime}\) to projin.wa_order. If sum_wa is True, or output projout.wa_order is None, sum over \(n\) (always the case except for debugging purposes). For example, see q. D5 and D6 of arXiv:1810.05051.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

class pypower.smooth_window.PowerSpectrumSmoothWindow(edges, modes, power_nonorm, nmodes, projs, **kwargs)

Bases: pypower.fft_power.BasePowerSpectrumStatistics

Power spectrum window function multipoles.

Initialize PowerSpectrumSmoothWindow.

Parameters

  • edges (tuple of ndim arrays) – Edges used to bin window function measurement.

  • modes (array) – Mean “wavenumber” (\(k\)) in each bin.

  • power_nonorm (array) – Power spectrum in each bin, without normalization.

  • nmodes (array) – Number of modes in each bin.

  • projs (list) – List of Projection instances or (multipole, wide-angle order) tuples.

  • kwargs (dict) – Other arguments for BasePowerSpectrumStatistics.

classmethod concatenate_proj(*others)

Concatenate input window functions, along projections.

Parameters

others (list of PowerSpectrumSmoothWindow) – List of window functions to be concatenated.

Returns

new

Return type

PowerSpectrumSmoothWindow

classmethod concatenate_x(*others, select='nmodes', frac_nyq=None)

Concatenate input window functions, along k-coordinates. k-edges and k-coordinates are taken from the first provided window; then expanded by those of other provided windows, if they cover a wider range. Eventually, for each bin (low, high), loop through all windows and select that with the largest number of modes in the given bin, if exists (bins are declared equal when exact floating point matching of (low, high)). In case two windows have the same number of modes in the same bin, the first provided one is selected. Therefore, different results may be obtained when changing the order of input windows.

Note

Typically, you will want to input windows with decreasing box size (largest box size first).

Parameters

  • others (list of PowerSpectrumSmoothWindow) – List of window functions to be concatenated.

  • select (string, default='nmodes') – How to select input windows for each k (if several); ‘nmodes’: select window with highest number of modes.

  • frac_nyq (float, tuple, default=None) – Optionally, fraction of Nyquist frequency where to cut input windows (e.g. 0.8). If a float, the same for all input windows; else a tuple or a list of such fraction for each input window.

Returns

new

Return type

PowerSpectrumSmoothWindow

classmethod from_power(power, wa_order=0, **kwargs)

Build window function from input PowerSpectrumMultipoless.

Parameters

  • power (PowerSpectrumMultipoless) – Power spectrum measurement to convert into PowerSpectrumSmoothWindow.

  • wa_order (int, default=0) – Wide-angle order used for input power spectrum measurement.

Returns

window

Return type

PowerSpectrumSmoothWindow

get_power(add_direct=True, remove_shotnoise=True, null_zero_mode=False, divide_wnorm=True, complex=True)

Return power spectrum, computed using various options.

Parameters

  • add_direct (bool, default=True) – Add direct power spectrum measurement.

  • remove_shotnoise (bool, default=True) – Remove estimated shot noise.

  • null_zero_mode (bool, default=True) – Remove power spectrum at \(k = 0\) (if within edges).

  • divide_wnorm (bool, default=True) – Divide by estimated power spectrum normalization.

  • complex (bool, default=True) – Whether (True) to return the complex power spectrum, or (False) return its real part if even multipoles, imaginary part if odd multipole.

  • Results –

  • ------- –

  • power (array) –

property k

Wavenumbers.

property kavg

Mode-weighted average wavenumber = k.

property kedges

Wavenumber edges.

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 ndim

Return binning dimensionality.

property power

Power spectrum, normalized and with shot noise removed.

rebin(factor=1)

Rebin power spectrum estimation in place, by factor(s) factor. A tuple must be provided in case ndim is greater than 1. Input factors must divide shape.

save(filename)

Save to filename.

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

Save power spectrum 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)
property shape

Return shape of binned power spectrum power.

property shotnoise

Normalized shot noise.

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.
classmethod sum(*others)

Sum input power spectra, weighted by their wnorm.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

to_real(**kwargs)

Transform power spectrum window function to configuration space.

Parameters

kwargs (dict) – Arguments for power_to_correlation_window().

Returns

window

Return type

CorrelationFunctionSmoothWindow

property with_mpi

Whether to use MPI.

class pypower.smooth_window.PowerSpectrumSmoothWindowMatrix(kout, projsin, projsout=None, k=None, kin_rebin=1, kin_lim=(0.0001, 1.0), sep=None, window=None, xy=1, q=0, sum_wa=True, default_zero=False, attrs=None)

Bases: pypower.wide_angle.BaseMatrix

Class computing matrix for window convolution in Fourier space.

Initialize PowerSpectrumSmoothWindowMatrix.

Parameters

  • kout (array) – Output wavenumbers.

  • projsin (list) – Input projections.

  • projsout (list, default=None) – Output projections. Defaults to propose_out(projsin, sum_wa=sum_wa).

  • k (array, default=None) – Wavenumber for Hankel transforms; must be log-spaced. If None, use sep and xy instead to determine k.

  • kin_rebin (tuple, default=1) – To save some memory, rebin along input k-coordinates by this factor.

  • kin_lim (tuple, default=(1e-4, 1.)) – To save some memory, pre-cut input k-coordinates to these limits.

  • sep (array, default=None) – Separations for Hankel transforms; must be log-spaced. If None, use k and xy instead to determine sep.

  • window (CorrelationFunctionSmoothWindow, PowerSpectrumSmoothWindow) – Window function to convolve power spectrum with. If a PowerSpectrumSmoothWindow instance is provided, it is transformed to configuration space.

  • xy (float, default=1) – If one of k or sep is None, set it following e.g. xy/sep[::-1].

  • q (int, default=0) – Power-law tilt to regularize Hankel transforms.

  • sum_wa (bool, default=True) – Whether to perform summation over output wide-angle orders. Always set to True except for debugging purposes.

  • default_zero (bool, default=False) – If a given projection is not provided in window function, set to 0. Else an IndexError is raised.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

propose_out(sum_wa=True)

Propose output projections given proposed input projections projsin. If sum_wa is True (typically always the case), return projections with Projection.wa_order set to None (all wide-angle orders have been summed).

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

resum_input_odd_wide_angle(**kwargs)

Resum odd wide-angle orders. By default, line-of-sight is chosen as that save in attrs (attrs['los_type']). To override, use input kwargs which will be passed to PowerSpectrumOddWideAngleMatrix.

run()

Set matrix. Provided arXiv:2106.06324 eq. 2.5:

\[W_{\ell\ell^{\prime}}^{(n)}(k) = \frac{2}{\pi} i^{\ell} (-i)^{\ell^{\prime}} \int ds s^{2} j_{\ell}(ks) j_{\ell^{\prime}}(k^{\prime}s) \sum_{L} C_{\ell \ell^{\prime} L} Q_{L}^{(n)}(s)\]

with \(\ell\) corresponding to projout.ell and \(\ell^{\prime}\) to projin.ell, \(k\) to kout and \(k^{\prime}\) to kin. \(n\) is the wide-angle order proj.wa_order. Yet, to avoid bothering with complex values, we only work with the imaginary part of odd power spectra (input and output). In addition, we include the \(dk^{\prime} k^{\prime 2}\) volume element (arXiv:2106.06324 eq. 2.7). Hence we actually implement:

\[W_{\ell,\ell^{\prime}}^{(n)}(k) = dk^{\prime} k^{\prime 2} \frac{2}{\pi} (-1)^{\ell/2} (-1)^{\ell^{\prime}/2} \int ds s^{2} j_{\ell}(ks) j_{\ell^{\prime}}(k^{\prime}s) \sum_{L} C_{\ell \ell^{\prime} L} Q_{L}^{(n)}(s)\]

Note that we do not include \(k^{-n}\) as this factor is included in PowerSpectrumOddWideAngleMatrix.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

pypower.smooth_window.power_to_correlation_window(fourier_window, sep=None, k=None, smooth=None)

Compute correlation window function by taking Hankel transforms of input power spectrum window function.

Parameters

  • fourier_window (PowerSpectrumSmoothWindow) – Power spectrum window function.

  • sep (array, default=None) – Separations \(s\) where to compute Hankel transform; defaults to inverse of fourier_window wavenumbers.

  • k (array, default=None) – Wavenumbers where to interpolate the window function. If provided, \(k\)-space volume element will be computed as \(4 \pi dk k^{2}\). Else, defaults to fourier_window.k and fourier_window.volume.

  • smooth (float, array, default=None) – If not None, if float, radius of Gaussian smoothing. Else, smoothing kernel, should be the same size as used k (see above).

Returns

window – Correlation window function.

Return type

CorrelationFunctionSmoothWindow

pypower.smooth_window.weights_trapz(x)

Return weights for trapezoidal integration.

pypower.smooth_window.wigner3j_square(ellout, ellin, prefactor=True)

Return the coefficients corresponding to the product of two Legendre polynomials, corresponding to \(C_{\ell \ell^{\prime} L}\) of e.g. eq. 2.2 of https://arxiv.org/pdf/2106.06324.pdf, with \(\ell\) corresponding to ellout and \(\ell^{\prime}\) to ellin.

Parameters

  • ellout (int) – Output order.

  • ellin (int) – Input order.

  • prefactor (bool, default=True) – Whether to include prefactor \((2 \ell + 1)/(2 L + 1)\) for window convolution.

Returns

  • ells (list) – List of mulipole orders \(L\).

  • coeffs (list) – List of corresponding window coefficients.

Accurate window

Implementation of window function estimation, following https://github.com/cosmodesi/GC_derivations, and https://fr.overleaf.com/read/hpgbwqzmtcxn.

class pypower.fft_window.CatalogFFTWindow(randoms_positions1=None, randoms_positions2=None, randoms_weights1=None, randoms_weights2=None, edgesin=None, projsin=None, edges=None, ells=None, power_ref=None, los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, dtype=None, resampler=None, interlacing=None, position_type='xyz', weight_type='auto', weight_attrs=None, wnorm=None, shotnoise=None, edgesin_type='smooth', mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Bases: pypower.fft_window.MeshFFTWindow

Wrapper on MeshFFTWindow to estimate window function from input random positions and weigths.

Initialize CatalogFFTWindow, i.e. estimate power spectrum window matrix.

Parameters

  • randoms_positions1 (list, array, default=None) – Positions in the first randoms catalog. Typically of shape (3, N) or (N, 3).

  • randoms_positions2 (list, array, default=None) – Optionally (for cross-correlation), positions in the second randoms catalog. See randoms_positions1.

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

  • randoms_weights2 (array of shape (N,), default=None) – Optionally (for cross-correlation), weights in the second randoms catalog.

  • edgesin (dict, array, list) – An array of \(k\)-edges which defines the theory \(k\)-binning; corresponding derivatives will be computed using get_correlation_function_tophat_derivative(); or a dictionary of such array for each theory projection. Else a list of derivatives (callable) of theory correlation function w.r.t. each theory basis vector, e.g. each in \(k\)-bin; or a dictionary of such list for each theory projection.

  • projsin (list, default=None) – List of Projection instances or (multipole, wide-angle order) tuples. If None, and power_ref is provided, the list of projections is set to be able to compute window convolution of theory power spectrum multipoles of orders power_ref.ells.

  • power_ref (PowerSpectrumMultipoles, default=None) – “Reference” power spectrum estimation, e.g. of the actual data. It is used to set default values for edges, ells, los, boxsize, boxcenter, nmesh, interlacing, resampler and wnorm if those are None.

  • edges (tuple, array, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu\)-edges (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). 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 find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’). For both \(k\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. bins[i] <= x < bins[i+1]. However, last \(\mu\)-bin is inclusive on both ends: bins[-2] <= mu <= bins[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\). If None, defaults to the edges used in estimation of power_ref.

  • ells (list, tuple, default=(0, 2, 4)) – Output multipole orders.

  • los (string, array, default=None) – 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. If None, defaults to line-of-sight used in estimation of power_ref.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis. If None, defaults to the value used in estimation of power_ref.

  • boxsize (float, default=None) – Physical size of the box, defaults to maximum extent taken by all input positions, times boxpad. If None, defaults to the value used in estimation of power_ref.

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

  • 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.

  • dtype (string, dtype, default=None) – The data type to use for input positions and weights and the mesh. If None, defaults to the value used in estimation of power_ref if provided, else ‘f8’.

  • resampler (string, ResampleWindow, default=None) – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’]. If None, defaults to the value used in estimation of power_ref.

  • interlacing (bool, int, default=None) – Whether to use interlacing to reduce aliasing when painting the particles on the mesh. If positive int, the interlacing order (minimum: 2). If None, defaults to the value used in estimation of power_ref.

  • 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)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, else “product_individual”.

    If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of weights if the denominator is zero (defaulting to 0).

  • wnorm (float, default=None) – Window function normalization. If None, defaults to the value used in estimation of power_ref, rescaled to the input random weights — which yields a correct normalization of the window function for the power spectrum estimation power_ref. If power_ref provided, use internal estimate obtained with normalization() — which is wrong (the normalization poles.wnorm can be reset a posteriori using the above recipe).

  • shotnoise (float, default=None) – Window function shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation and in case of auto-correlation is obtained by dividing CatalogMesh.unnormalized_shotnoise() by window function normalization.

  • edgesin_type (str, default='smooth') – Technique to transpose edgesin to Fourier space. ‘smooth’ uses get_correlation_function_tophat_derivative(); ‘fourier-grid’ paints edgesin on the Fourier mesh, then takes the FFT.

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

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

property boxsize

Physical box size.

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_window.MeshFFTWindow(mesh1=None, mesh2=None, edgesin=None, projsin=None, power_ref=None, edges=None, ells=None, los=None, periodic=False, boxcenter=None, compensations=None, wnorm=None, shotnoise=None, edgesin_type='smooth', **kwargs)

Bases: pypower.fft_power.MeshFFTPower

Class that computes window function from input mesh(es), using global or local line-of-sight, see:

poles

Window matrix.

Type

PowerSpectrumFFTWindowMatrix

Initialize MeshFFTWindow.

Parameters

  • mesh1 (CatalogMesh, RealField, default=None) – First mesh.

  • mesh2 (CatalogMesh, RealField, default=None) – In case of cross-correlation, second mesh, with same size and physical extent (boxsize and boxcenter) that mesh1.

  • edgesin (dict, array, list) – An array of \(k\)-edges which defines the theory \(k\)-binning; corresponding derivatives will be computed (see edgesin_type); or a dictionary of such array for each theory projection. Else a list of derivatives (callable) of theory correlation function w.r.t. each theory basis vector, e.g. each in \(k\)-bin; or a dictionary of such list for each theory projection. If periodic is True, this should correspond to the derivatives of theory power spectrum (instead of correlation function) w.r.t. each theory basis vector, e.g. each in \(k\) bin.

  • projsin (list, default=None) – List of Projection instances or (multipole, wide-angle order) tuples. If None, and power_ref is provided, the list of projections is set to be able to compute window convolution of theory power spectrum multipoles of orders power_ref.ells.

  • power_ref (CatalogFFTPower, MeshFFTPower, PowerSpectrumWedges, PowerSpectrumMultipoles, default=None) – “Reference” power spectrum estimation, e.g. of the actual data. It is used to set default values for edges, ells, los, boxcenter, compensations and wnorm if those are None.

  • edges (tuple, array, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu\)-edges (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). 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 find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’). For both \(k\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. edges[i] <= x < edges[i+1]. However, last \(\mu\)-bin is inclusive on both ends: edges[-2] <= mu <= edges[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\). If None, defaults to the edges used in estimation of power_ref.

  • ells (list, tuple, default=(0, 2, 4)) – Output multipole orders. If None, defaults to the multipoles used in estimation of power_ref.

  • los (string, array, default=None) – 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. If None, defaults to the line-of-sight used in estimation of power_ref.

  • periodic (bool, default=False) – If True, selection function is assumed uniform, periodic. In this case, mesh1 may be None; in this case nmesh and boxsize default to that of power_ref, else may be set with kwargs.

  • boxcenter (float, array, default=None) – Box center; defaults to 0. Used only if provided mesh1 and mesh2 are not CatalogMesh. If None, defaults to the value used in estimation of power_ref.

  • compensations (list, tuple, string, default=None) – Compensations to apply to mesh to (optionally) correct for particle-mesh assignment scheme; e.g. ‘cic’ (resp. ‘cic-sn’) for cic assignment scheme, with (resp. without) interlacing. In case mesh2 is not None (cross-correlation), provide a list (or tuple) of two such strings (for mesh1 and mesh2, respectively). Used only if provided mesh1 or mesh2 are not CatalogMesh.

  • wnorm (float, default=None) – Window function normalization. If None, defaults to the value used in estimation of power_ref, rescaled to the input random weights — which yields a correct normalization of the window function for the power spectrum estimation power_ref. If power_ref provided, use internal estimate obtained with normalization() — which is wrong (the normalization poles.wnorm can be reset a posteriori using the above recipe).

  • shotnoise (float, default=None) – Window function shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation or both mesh1 and mesh2 are pmesh.pm.RealField, and in case of auto-correlation is obtained by dividing CatalogMesh.unnormalized_shotnoise() of mesh1 by window function normalization.

  • edgesin_type (str, default='smooth') – Technique to transpose edgesin to Fourier space, relevant only if periodic is False. ‘smooth’ uses get_correlation_function_tophat_derivative(); ‘fourier-grid’ paints edgesin on the Fourier mesh (akin to the periodic case), then takes the FFT.

  • kwargs (dict) – Arguments for ParticleMesh in case mesh1 is not provided (as may be the case if periodic is True), typically boxsize, nmesh, mpicomm.

property boxsize

Physical box size.

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_window.PowerSpectrumFFTWindowMatrix(matrix, xin, xout, projsin, projsout, nmodes, wnorm=1.0, attrs=None, mpicomm=None)

Bases: pypower.wide_angle.BaseMatrix

Window matrix, relating “theory” input to “observed” output.

Initialize PowerSpectrumFFTWindowMatrix.

Parameters

  • matrix (array) – 2D array representing window matrix.

  • xin (array, list) – List of input “theory” coordinates. If single array, assumed to be the same for all input projections projsin.

  • xout (list) – List of output “theory” coordinates. If single array, assumed to be the same for all output projections projsout.

  • projsin (list) – List of input “theory” projections.

  • projsout (list) – List of output “observed” projections.

  • nmodes (array) – Number of modes in each bin.

  • wnorm (float, default=1.) – Window function normalization.

  • attrs (dict, default=None) – Dictionary of other attributes.

  • mpicomm (MPI communicator, default=None) – The MPI communicator, only used when saving (save()) matrix.

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

classmethod from_power(power, xin, projin=(0, 0), **kwargs)

Create window function from input PowerSpectrumMultipoles.

Parameters
Returns

matrix

Return type

PowerSpectrumFFTWindowMatrix

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

resum_input_odd_wide_angle(**kwargs)

Resum odd wide-angle orders. Input kwargs will be passed to PowerSpectrumOddWideAngleMatrix.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

pypower.fft_window.get_correlation_function_tophat_derivative(kedges, ell=0, k=None, **kwargs)

Return a list of callable corresponding to the derivative of the correlation function w.r.t. \(k\)-bins.

Parameters

  • kedges (array) – \(k\)-edges of the \(k\)-bins.

  • ell (int, default=0) – Multipole order.

  • k (array, default=None) – If None, calculation will be analytic, which will work if ell in [0, 2, 4], or sympy package is installed (such analytic integration with sympy may take several seconds). If not None, this is the \(k\) log-spaced array for numerical FFTlog integration.

  • kwargs (dict) – If k is not None, other arguments for fftlog.PowerToCorrelation.

Returns

toret – List of callables, taking configuration-space separation s as input.

Return type

list

Direct power spectrum estimator

Implementation of direct estimation of power spectrum multipoles, i.e. summing over particle pairs. This should be mostly used to sum over pairs at small separations, otherwise the calculation will be prohibitive.

class pypower.direct_power.BaseDirectPowerEngine(modes, positions1, positions2=None, weights1=None, weights2=None, ells=(0, 2, 4), limits=(0.0, 0.03333333333333333), limit_type='degree', position_type='xyz', weight_type='auto', weight_attrs=None, twopoint_weights=None, los='firstpoint', boxsize=None, dtype='f8', mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD, **kwargs)

Bases: pypower.utils.BaseClass

Direct power spectrum measurement, summing over particle pairs.

Initialize BaseDirectPowerEngine.

Parameters

  • modes (array) – Wavenumbers at which to compute power spectrum.

  • positions1 (list, array) – Positions in the first data catalog. Typically of shape (3, N) or (N, 3).

  • positions2 (list, array, default=None) – Optionally, for cross-power spectrum, positions in the second catalog. See positions1.

  • weights1 (array, list, default=None) – Weights of the first catalog. Not required if weight_type is either None or “auto”. See weight_type.

  • weights2 (array, list, default=None) – Optionally, for cross-pair counts, weights in the second catalog. See weights1.

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • limits (tuple, default=(0., 2./60.)) – Limits of particle pair separations.

  • limit_type (string, default='degree') – Type of limits; i.e. are those angular limits (“degree”, “radian”), or 3D limits (“s”)?

  • 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)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”inverse_bitwise”: each pair is weighted by \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\).

      Multiple bitwise weights can be provided as a list. Individual weights can additionally be provided as float arrays. In case of cross-correlations with floating weights, bitwise weights are automatically turned to IIP weights, i.e. \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1}))\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, “inverse_bitwise” if one of input weights is integer, else “product_individual”.

    In addition, angular upweights can be provided with twopoint_weights. If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). Inverse probability weight is then computed as: \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0.

  • twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles. A WeightTwoPointEstimator instance (from pycorr) or any object with arrays sep (separations) and weight (weight at given separation) as attributes (i.e. to be accessed through twopoint_weights.sep, twopoint_weights.weight) or as keys (i.e. twopoint_weights['sep'], twopoint_weights['weight']) or as element (i.e. sep, weight = twopoint_weights).

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

  • dtype (string, dtype, default='f8') – The data type to use for input positions and weights.

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

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

run()

Method that computes the power spectrum and set power_nonorm, to be implemented in your new engine.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.direct_power.CorrfuncDirectPowerEngine(modes, positions1, positions2=None, weights1=None, weights2=None, ells=(0, 2, 4), limits=(0.0, 0.03333333333333333), limit_type='degree', position_type='xyz', weight_type='auto', weight_attrs=None, twopoint_weights=None, los='firstpoint', boxsize=None, dtype='f8', mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD, **kwargs)

Bases: pypower.direct_power.BaseDirectPowerEngine

Direct power spectrum measurement, using Corrfunc.

Initialize BaseDirectPowerEngine.

Parameters

  • modes (array) – Wavenumbers at which to compute power spectrum.

  • positions1 (list, array) – Positions in the first data catalog. Typically of shape (3, N) or (N, 3).

  • positions2 (list, array, default=None) – Optionally, for cross-power spectrum, positions in the second catalog. See positions1.

  • weights1 (array, list, default=None) – Weights of the first catalog. Not required if weight_type is either None or “auto”. See weight_type.

  • weights2 (array, list, default=None) – Optionally, for cross-pair counts, weights in the second catalog. See weights1.

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • limits (tuple, default=(0., 2./60.)) – Limits of particle pair separations.

  • limit_type (string, default='degree') – Type of limits; i.e. are those angular limits (“degree”, “radian”), or 3D limits (“s”)?

  • 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)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”inverse_bitwise”: each pair is weighted by \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\).

      Multiple bitwise weights can be provided as a list. Individual weights can additionally be provided as float arrays. In case of cross-correlations with floating weights, bitwise weights are automatically turned to IIP weights, i.e. \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1}))\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, “inverse_bitwise” if one of input weights is integer, else “product_individual”.

    In addition, angular upweights can be provided with twopoint_weights. If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). Inverse probability weight is then computed as: \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0.

  • twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles. A WeightTwoPointEstimator instance (from pycorr) or any object with arrays sep (separations) and weight (weight at given separation) as attributes (i.e. to be accessed through twopoint_weights.sep, twopoint_weights.weight) or as keys (i.e. twopoint_weights['sep'], twopoint_weights['weight']) or as element (i.e. sep, weight = twopoint_weights).

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

  • dtype (string, dtype, default='f8') – The data type to use for input positions and weights.

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

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

run()

Method that computes the power spectrum and set power_nonorm, to be implemented in your new engine.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.direct_power.DirectPower(*args, engine='corrfunc', **kwargs)

Bases: object

Entry point to direct power engines.

Parameters

  • engine (string, default='kdtree') – Name of direct power engine, one of [‘kdtree’].

  • args (list) – Arguments for direct power engine, see BaseDirectPowerEngine.

  • kwargs (dict) – Arguments for direct power engine, see BaseDirectPowerEngine.

Returns

engine

Return type

BaseDirectPowerEngine

class pypower.direct_power.KDTreeDirectPowerEngine(modes, positions1, positions2=None, weights1=None, weights2=None, ells=(0, 2, 4), limits=(0.0, 0.03333333333333333), limit_type='degree', position_type='xyz', weight_type='auto', weight_attrs=None, twopoint_weights=None, los='firstpoint', boxsize=None, dtype='f8', mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD, **kwargs)

Bases: pypower.direct_power.BaseDirectPowerEngine

Direct power spectrum measurement, summing over particle pairs, identified with KDTree.

Initialize BaseDirectPowerEngine.

Parameters

  • modes (array) – Wavenumbers at which to compute power spectrum.

  • positions1 (list, array) – Positions in the first data catalog. Typically of shape (3, N) or (N, 3).

  • positions2 (list, array, default=None) – Optionally, for cross-power spectrum, positions in the second catalog. See positions1.

  • weights1 (array, list, default=None) – Weights of the first catalog. Not required if weight_type is either None or “auto”. See weight_type.

  • weights2 (array, list, default=None) – Optionally, for cross-pair counts, weights in the second catalog. See weights1.

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • limits (tuple, default=(0., 2./60.)) – Limits of particle pair separations.

  • limit_type (string, default='degree') – Type of limits; i.e. are those angular limits (“degree”, “radian”), or 3D limits (“s”)?

  • 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)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”inverse_bitwise”: each pair is weighted by \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\).

      Multiple bitwise weights can be provided as a list. Individual weights can additionally be provided as float arrays. In case of cross-correlations with floating weights, bitwise weights are automatically turned to IIP weights, i.e. \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1}))\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, “inverse_bitwise” if one of input weights is integer, else “product_individual”.

    In addition, angular upweights can be provided with twopoint_weights. If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). Inverse probability weight is then computed as: \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0.

  • twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles. A WeightTwoPointEstimator instance (from pycorr) or any object with arrays sep (separations) and weight (weight at given separation) as attributes (i.e. to be accessed through twopoint_weights.sep, twopoint_weights.weight) or as keys (i.e. twopoint_weights['sep'], twopoint_weights['weight']) or as element (i.e. sep, weight = twopoint_weights).

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

  • dtype (string, dtype, default='f8') – The data type to use for input positions and weights.

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

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

run()

Method that computes the power spectrum and set power_nonorm, to be implemented in your new engine.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.direct_power.MetaDirectPower(name, bases, class_dict)

Bases: pypower.utils.BaseMetaClass

Metaclass to return correct direct power engine.

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 pypower.direct_power.RegisteredDirectPowerEngine(name, bases, class_dict)

Bases: pypower.utils.BaseMetaClass

Metaclass registering BaseDirectPowerEngine-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

pypower.direct_power.get_default_nrealizations(weights)

Return default number of realizations given input bitwise weights = the number of bits in input weights plus one.

pypower.direct_power.get_direct_power_engine(engine='corrfunc')

Return BaseDirectPowerEngine-subclass corresponding to input engine name.

Parameters

engine (string, default='kdtree') – Name of direct power engine, one of [‘kdtree’].

Returns

engine – Direct power engine class.

Return type

type

pypower.direct_power.get_inverse_probability_weight(*weights, noffset=1, nrealizations=None, default_value=0.0, dtype='f8')

Return inverse probability weight given input bitwise weights. Inverse probability weight is computed as: \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2} \& ...))\). If denominator is 0, weight is set to default_value.

Parameters

  • weights (int arrays) – Bitwise weights.

  • noffset (int, default=1) – The offset to be added to the bitwise counts in the denominator (defaults to 1).

  • nrealizations (int, default=None) – Number of realizations (defaults to the number of bits in input weights plus one).

  • default_value (float, default=0.) – Default weight value, if the denominator is zero (defaults to 0).

  • dtype (string, np.dtype) – Type for output weight.

Returns

weight – IIP weight.

Return type

array

Wide-angle corrections

Implementation of odd wide-angle matrices: - CorrelationFunctionOddWideAngleMatrix for correlation function - PowerSpectrumOddWideAngleMatrix for power spectrum, following https://arxiv.org/abs/2106.06324.

class pypower.wide_angle.BaseMatrix(value, xin, xout, projsin, projsout, weightsin=None, weightsout=None, attrs=None)

Bases: pypower.utils.BaseClass

Base class to represent a linear transform of the theory model, from input projections projsin to output projections projsout.

matrix

2D array representing linear transform. First axis is input, second is output.

Type

array

xin

List of input “theory” coordinates.

Type

list

xout

List of output “theory” coordinates.

Type

list

projsin

List of input “theory” projections.

Type

list

projsout

List of output “observed” projections.

Type

list

Initialize BaseMatrix.

Parameters

  • value (array) – 2D array representing linear transform.

  • xin (array, list) – List of input “theory” coordinates. If single array, assumed to be the same for all input projections projsin.

  • xout (list) – List of output “theory” coordinates. If single array, assumed to be the same for all output projections projsout.

  • projsin (list) – List of input “theory” projections.

  • projsout (list) – List of output “observed” projections.

  • weightsin (array, list, default=None) – Optionally, list of weights to apply when rebinning input “theory” coordinates.

  • weightsout (array, list, default=None) – Optionally, list of weights to apply when rebinning output “observed” coordinates.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

class pypower.wide_angle.CorrelationFunctionOddWideAngleMatrix(sep, projsin, projsout=None, wa_orders=1, los='firstpoint', attrs=None)

Bases: pypower.wide_angle.BaseMatrix

Class computing matrix for odd wide-angle expansion of the correlation function.

Initialize CorrelationFunctionOddWideAngleMatrix.

Parameters

  • k (array) – Input (and ouput) separations.

  • projsin (list) – Input projections.

  • projsout (list, default=None) – Output projections. Defaults to propose_out(projsin, wa_orders=wa_orders). If output projections have Projection.wa_order None, wide-angle orders are summed over.

  • wa_orders (int, list) – Wide-angle expansion orders. So far order 1 only is supported.

  • los (string) –

    Choice of line-of-sight, either:

    • ’firstpoint’: the separation vector starts at the end of the line-of-sight

    • ’endpoint’: the separation vector ends at the end of the line-of-sight.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

static propose_out(projsin, wa_orders=1)

Propose output projections (i.e. multipoles at wide-angle order > 0) that can be computed given proposed input projections projsin.

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

run()

Set matrix:

\[M_{\ell\ell^{\prime}}^{(n,n^{\prime})}(s) = - \frac{\ell \left(\ell - 1\right)}{2 \ell \left(2 \ell - 1\right)} \delta_{\ell,\ell - 1} \delta_{n^{\prime},0} + \frac{\left(\ell + 1\right) \left(\ell + 2\right)}{2 \ell \left(2 \ell + 3\right)} \delta_{\ell,\ell + 1} \delta_{n^{\prime},0}\]

if \(\ell\) is odd and \(n = 1\), else:

\[M_{\ell\ell^{\prime}}^{(0,n^{\prime})}(s) = \delta_{\ell,ell^{\prime}} \delta_{n^{\prime},0}\]

with \(\ell\) multipole order corresponding to projout.ell and \(\ell^{\prime}\) to projin.ell, \(n\) wide angle order corresponding to projout.wa_order and \(n^{\prime}\) to projin.wa_order. If output projout.wa_order is None, sum over \(n\) (correct only if no window convolution is accounted for).

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

class pypower.wide_angle.PowerSpectrumOddWideAngleMatrix(k, projsin, projsout=None, d=1.0, wa_orders=1, los='firstpoint', attrs=None)

Bases: pypower.wide_angle.BaseMatrix

Class computing matrix for odd wide-angle expansion of the power spectrum. Adapted from https://github.com/fbeutler/pk_tools/blob/master/wide_angle_tools.py

Initialize PowerSpectrumOddWideAngleMatrix.

Parameters

  • k (array) – Input (and ouput) wavenumbers.

  • projsin (list) – Input projections.

  • projsout (list, default=None) – Output projections. Defaults to propose_out(projsin, wa_orders=wa_orders). If output projections have Projection.wa_order None, wide-angle orders are summed over.

  • d (float, default=1) – Distance at the effective redshift. Use \(1\) if already included in window functions.

  • wa_orders (int, list) – Wide-angle expansion orders. So far order 1 only is supported.

  • los (string) –

    Choice of line-of-sight, either:

    • ’firstpoint’: the separation vector starts at the end of the line-of-sight

    • ’endpoint’: the separation vector ends at the end of the line-of-sight.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns

matrix – New matrix, of same type as others[0].

Return type

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

propose_out(wa_orders=1)

Propose output projections (i.e. multipoles at wide-angle order > 0) that can be computed given proposed input projections projsin.

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

run()

Set matrix:

\[M_{\ell\ell^{\prime}}^{(n,n^{\prime})}(k) = - \frac{\ell \left(\ell - 1\right)}{2 \ell \left(2 \ell - 1\right) d} \delta_{\ell,\ell - 1} \delta_{n^{\prime},0} \left[\frac{\ell - 1}{k} - \partial_{k} \right] - \frac{\left(\ell + 1\right) \left(\ell + 2\right)}{2 \ell \left(2 \ell + 3\right) d} \delta_{\ell,\ell + 1} \delta_{n^{\prime},0} \left[ \frac{\ell + 2}{k} + \partial_{k} \right]\]

if \(\ell\) is odd and \(n = 1\), else:

\[M_{\ell\ell^{\prime}}^{(0,n^{\prime})}(k) = \delta_{\ell,ell^{\prime}} \delta_{n^{\prime},0}\]

with \(\ell\) multipole order corresponding to projout.ell and \(\ell^{\prime}\) to projin.ell, \(n\) wide angle order corresponding to projout.wa_order and \(n^{\prime}\) to projin.wa_order. If output projout.wa_order is None, sum over \(n\) (correct only if no window convolution is accounted for). Derivatives \(\partial_{k}\) are computed with finite differences, see arXiv:2106.06324 eq. 3.3.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

class pypower.wide_angle.Projection(ell, wa_order='default', default_wa_order=0)

Bases: pypower.utils.BaseClass

Class representing a “projection”, i.e. multipole and wide-angle expansion order.

ell

Multipole order.

Type

int

wa_order

Wide-angle order.

Type

int, None

Initialize Projection.

Parameters

  • ell (int) – Multipole order.

  • wa_order (int, None, default='default') – Wide-angle order. If ‘default’, defaults to default_wa_order.

  • default_wa_order (int, default=0) – Default wide-angle order to use if wa_order is ‘default’.

clone(**kwargs)

Clone current projection, optionally updating ell and wa_order (using kwargs).

latex(inline=False)

Return latex string for current projection. If inline is True, add surrounding dollar $ signs.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

pypower.wide_angle.odd_wide_angle_coefficients(ell, wa_order=1, los='firstpoint')

Compute coefficients of odd wide-angle expansion, i.e.:

\[- \frac{\ell \left(\ell - 1\right)}{2 \ell \left(2 \ell - 1\right)}, \frac{\left(\ell + 1\right) \left(\ell + 2\right)}{2 \ell \left(2 \ell + 3\right)}\]

For the first point line-of-sight. See https://fr.overleaf.com/read/hpgbwqzmtcxn. A minus sign is applied on both factors if los is ‘endpoint’.

Parameters

  • ell (int) – (Odd) multipole order.

  • wa_order (int, default=1) – Wide-angle expansion order. So far only order 1 is supported.

  • los (string) –

    Choice of line-of-sight, either:

    • ’firstpoint’: the separation vector starts at the end of the line-of-sight

    • ’endpoint’: the separation vector ends at the end of the line-of-sight.

Returns

  • ells (list) – List of multipole orders of correlation function.

  • coeffs (list) – List of coefficients to apply to correlation function multipoles corresponding to output ells.

Utilities

A few utilities.

class pypower.utils.BaseClass

Bases: object

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

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

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

Bases: type

Metaclass 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

pypower.utils.cartesian_to_sky(positions, wrap=True, degree=True)

Transform cartesian coordinates into distance, RA, Dec.

Parameters

  • positions (array of shape (3, N), list of 3 arrays) – Positions 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

rdd – Right ascension, declination and distance.

Return type

list of 3 arrays

pypower.utils.distance(positions)

Return cartesian distance, taking coordinates along position first axis.

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

Print exception with a logger.

pypower.utils.mkdir(dirname)

Try to create dirnm and catch OSError.

pypower.utils.pack_bitarrays(*arrays, dtype=<class 'numpy.uint64'>)

Pack bit arrays into a list of integer arrays. Inverse operation is unpack_bitarray(), i.e. unpack_bitarrays(pack_bitarrays(*arrays, dtype=dtype))``is ``arrays, whatever integer dtype is.

Parameters

  • arrays (bool arrays) – Arrays of integers or booleans whose elements should be packed to bits.

  • dtype (string, dtype) – Type of output integer arrays.

Returns

arrays – List of integer arrays of type dtype, representing input boolean arrays.

Return type

list

pypower.utils.popcount(*arrays)

Return number of 1 bits in each value of input array. Inspired from https://github.com/numpy/numpy/issues/16325.

pypower.utils.rebin(array, new_shape, statistic=<function sum>)

Bin an array in all axes based on the target shape, by summing or averaging. Number of output dimensions must match number of input dimensions and new axes must divide old ones.

Taken from https://stackoverflow.com/questions/8090229/resize-with-averaging-or-rebin-a-numpy-2d-array and https://nbodykit.readthedocs.io/en/latest/_modules/nbodykit/binned_statistic.html#BinnedStatistic.reindex.

Example

>>> m = np.arange(0,100,1).reshape((10,10))
>>> n = rebin(m, new_shape=(5,5), statistic=np.sum)
>>> print(n)
[[ 22 30 38 46 54]

[102 110 118 126 134] [182 190 198 206 214] [262 270 278 286 294] [342 350 358 366 374]]

pypower.utils.reformat_bitarrays(*arrays, dtype=<class 'numpy.uint64'>, copy=True)

Reformat input integer arrays into list of arrays of type dtype. If, e.g. 6 arrays of type np.uint8 are input, and dtype is np.uint32, a list of 2 arrays is returned.

Parameters

  • arrays (integer arrays) – Arrays of integers to reformat.

  • dtype (string, dtype) – Type of output integer arrays.

  • copy (bool, default=True) – If False, avoids copy of input arrays if dtype is uint8.

Returns

arrays – List of integer arrays of type dtype, representing input integer arrays.

Return type

list

pypower.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().

pypower.utils.sky_to_cartesian(rdd, degree=True, dtype=None)

Transform distance, RA, Dec into cartesian coordinates.

Parameters

  • rdd (array of shape (3, N), list of 3 arrays) – Right ascension, declination and distance.

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

Returns

positions – Positions x, y, z in cartesian coordinates.

Return type

list of 3 arrays

pypower.utils.unpack_bitarrays(*arrays)

Unpack integer arrays into a bit array. Inverse operation is pack_bitarray(), i.e. pack_bitarrays(unpack_bitarrays(*arrays), dtype=arrays.dtype)``is ``arrays.

Parameters

arrays (integer arrays) – Arrays of integers whose elements should be unpacked to bits.

Returns

arrays – List of boolean arrays of type np.uint8, representing input integer arrays.

Return type

list

MPI

pypower.mpi.domain_decompose(mpicomm, smoothing, positions1, weights1=None, positions2=None, weights2=None, boxsize=None, domain_factor=None)

Adapted from https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/pair_counters/domain.py. Decompose positions and weights on a grid of MPI processes. Requires mpi4py and pmesh.

Parameters

  • mpicomm (MPI communicator) – The MPI communicator.

  • smoothing (float) – The maximum Cartesian separation implied by the user’s binning.

  • positions1 (array of shape (N, 3)) – Positions in the first catalog.

  • positions2 (array of shape (N, 3), default=None) – Optionally, for cross-pair counts, positions in the second catalog. See positions1.

  • weights1 (list, array, default=None) – Optionally, weights of the first catalog.

  • weights2 (list, array, default=None) – Optionally, weights in the second catalog.

  • boxsize (array, default=None) – For periodic wrapping, the 3 side-lengths of the periodic cube.

  • domain_factor (int, default=None) – Multiply the size of the MPI mesh by this factor. If None, defaults to 2 in case boxsize is None, else (periodic wrapping) 1.

Returns

(positions1, weights1), (positions2, weights2) – The (decomposed) set of positions and weights.

Return type

arrays

pypower.mpi.gather_array(data, root=0, mpicomm=mpi4py.MPI.COMM_WORLD)

Taken from https://github.com/bccp/nbodykit/blob/master/nbodykit/utils.py Gather the input data array from all ranks to the specified root. This uses Gatherv, which avoids mpi4py pickling, and also avoids the 2 GB mpi4py limit for bytes using a custom datatype

Parameters

  • data (array_like) – The data on each rank to gather.

  • root (int, Ellipsis, default=0) – The rank number to gather the data to. If root is Ellipsis or None, broadcast the result to all ranks.

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

Returns

recvbuffer – The gathered data on root, and None otherwise.

Return type

array_like, None

pypower.mpi.local_size(size, mpicomm=mpi4py.MPI.COMM_WORLD)

Divide global size into local (process) size.

Parameters

  • size (int) – Global size.

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

Returns

localsize – Local size. Sum of local sizes over all processes equals global size.

Return type

int

pypower.mpi.scatter_array(data, counts=None, root=0, mpicomm=mpi4py.MPI.COMM_WORLD)

Taken from https://github.com/bccp/nbodykit/blob/master/nbodykit/utils.py Scatter the input data array across all ranks, assuming data is initially only on root (and None on other ranks). This uses Scatterv, which avoids mpi4py pickling, and also avoids the 2 GB mpi4py limit for bytes using a custom datatype

Parameters

  • data (array_like or None) – On root, this gives the data to split and scatter.

  • counts (list of int) – List of the lengths of data to send to each rank.

  • root (int, default=0) – The rank number that initially has the data.

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

Returns

recvbuffer – The chunk of data that each rank gets.

Return type

array_like