Hybrid-Basis Inference for Large-Scale Galaxy Clustering

arXiv eprint GitHub release (latest by date) Documentation status Build status Licence

Harmonia is a Python package that combines clustering statistics decomposed in spherical and Cartesian Fourier bases for large-scale galaxy clustering likelihood analysis.

Installation

We recommend that you install nbodykit first by following these instructions.

After that, you can install Harmonia simply using pip:

pip install harmoniacosmo

Note that only here does the name harmoniacosmo appear because unfortunately on PyPI the project name harmonia has already been taken.

Documentation

  • Recipes (under construction): tutorials in the format of integrated notebooks showcasing the use of Harmonia will be gradually added, so look out for any updates! For now, application offers some scripts that illustrate the use of Harmonia.

  • API Reference: more detailed documentation of classes and functions.

Attribution

If you would like to acknowledge this work, please cite Wang et al. (2020). You may use the following BibTeX record:

@article{Wang_2020b,
    author={Wang, M.~S. and Avila, S. and Bianchi, D. and Crittenden, R. and Percival, W.~J.},
    title={Hybrid-basis inference for large-scale galaxy clustering: combining spherical and {Cartesian} {Fourier} analyses},
    year={2020},
    eprint={2007.14962},
    archivePrefix={arXiv},
    primaryClass={astro-ph.CO},
}

Licence

Harmonia is made freely available under the GPL v3.0 licence.

Recipes

Tutorials in the format of integrated notebooks showcasing the use of Harmonia will be gradually added, so look out for any updates! For now, application offers some scripts that illustrate the use of Harmonia.

API Reference

Algorithms (algorithms)

Provide algorithms for Fourier basis evaluation and integration, spectrum discretisation and structured array manipulation for cosmological data.

Structured arrays (arrays)

Provide structured arrays for cosmological data.

DataArray()

Abstract data array with save and load methods.

SphericalArray(disc)

Structured array for spherically decomposed cosmological data.

CartesianArray(orders, wavenumbers[, …])

Structured array for Cartesian decomposition of cosmological data.


exception harmonia.algorithms.arrays.IndexingError[source]

Exception raised for unsupported slicing or indexing in __getitem__ methods.

class harmonia.algorithms.arrays.DataArray[source]

Abstract data array with save and load methods.

save(output_file, file_extension)[source]

Save the structured array.

Parameters
  • output_file (str or pathlib.Path) – Output file path.

  • extension ({‘.pkl’, ‘.npz’}) – Output file extension.

classmethod load(input_file)[source]

Load the structured array from a .npz or .pkl file.

Parameters

input_file (str or pathlib.Path) – Input file path.

class harmonia.algorithms.arrays.SphericalArray(disc)[source]

Structured array for spherically decomposed cosmological data.

Array is initialised with a discrete spectrum of Fourier modes and consists of three fields: the ‘index’ field of \((\ell, m_\ell, n_\ell)\) triplets, the ‘wavenumber’ field of discrete \(k_{\ell n}\), and the ‘coefficient’ field of spherically decomposed data.

Parameters

disc (DiscreteSpectrum) – Discrete spectrum associated with the structured array.

Variables
  • array (numpy.ndarray) – Structured NumPy array.

  • size (int) – Total number of elements in the array. This should equal the sum of disc.mode_counts.

  • attrs (dict) – Attributes of the structured array inherited from disc.

See also

DiscreteSpectrum

__getitem__(key)[source]

Get the ‘coefficient’ field value(s).

The access key can be an integer, a slice expression, a tuple of index triplet or a string, e.g. [-1], [:], [(0, 0, 1)] or 'degree_0'.

Parameters

key (int, slice, tuple(int, int, int) or str) – ‘coefficient’ field access key.

Returns

‘coefficient’ field data entry.

Return type

complex

__setitem__(key, value)[source]

Set the ‘coefficient’ field value(s).

Parameters
  • key (int, tuple of int or slice) – ‘coefficient’ field access key.

  • value (complex) – ‘coefficient’ field data entry.

vectorise(pivot, collapse=None)[source]

Returrn a data vector from the ‘coefficient’ field.

Vectorisation is performed by pivoting in either of the following orders of precedence—

  • ‘natural’: ordered by \((\ell, m, n)\);

  • ‘spectral’: ordered by \((k_{\ell n}, m)\).

Subarrays of equivalent \((\ell, n)\) may be further collapsed over spherical order \(m\) by simple averaging or averaging in quadrature.

Parameters
  • pivot ({‘natural’, ‘spectral’}) – Pivot order for vectorisation.

  • collapse ({None, ‘mean’, ‘qaudrature’}, optional) – If not None (default), subarrays are collapsed over equivalent spherical order \(m\) by averaging (‘mean’) or averaging in quadrature (‘qaudrature’).

Returns

vectorised_data – Vectorised coefficient data.

Return type

numpy.ndarray

class harmonia.algorithms.arrays.CartesianArray(orders, wavenumbers, mode_counts=None, shot_noise=None)[source]

Structured array for Cartesian decomposition of cosmological data.

Array is initialised with three fields: the ‘order’ field of the Legendre multipoles, the ‘wavenumber’ field of \(k\)-bin centres and the ‘power’ field for power spectrum multipole measurements.

Parameters
  • orders (list of tuple of int) – Orders of the power spectrum multipole.

  • wavenumbers (float, array_like) – Wavenumbers of the multipole data.

  • mode_counts (list of int or None, optional) – Mode counts in wavenumber bins (default is None).

  • shot_noise (float or None, optional) – Shot noise level (default is None).

Variables
  • array (numpy.ndarray) – Structured NumPy array.

  • size (int) – Total number of elements in the array. This should equal the product of the numbers of wavenumbers and multipoles.

  • attrs (dict) – Initialisation parameters as attributes.

__getitem__(key)[source]

Access the ‘power’ field.

The access key can be an integer positional index, a slice, a tuple of (order, wavenumber) or a string, e.g. [-1], [:], [(0, 0.04)] or 'power_0'.

Parameters

key (int, slice, tuple(int, float) or str) – ‘power’ field access key.

Returns

‘power’ field data entry.

Return type

float

__setitem__(key, value)[source]

Set the ‘power’ field value(s).

Parameters
  • key (int, tuple(int, float) or slice) – ‘power’ field access key.

  • value (float) – ‘power’ field data entry.

vectorise(pivot)[source]

Return a data vector from the ‘power’ field.

Vectorisation is performed by pivoting in either of the following orders of precedence—

  • ‘order’: ordered by multipole order \(\ell\);

  • ‘wavenumber’: ordered by wavenumber \(k\).

Parameters

pivot ({‘order’, ‘wavenumber’}) – Pivot order for vectorisation.

Returns

vectorised_data – Vectorised power spectrum data.

Return type

numpy.ndarray

Fourier bases (bases)

Compute quantities related to Fourier basis functions.

spherical_harmonic(ell, m, theta, phi[, conj])

Evaluate the spherical harmonic function.

spherical_besselj(ell, x[, derivative])

Evaluate the spherical Bessel function of the first kind or its derivative.

spherical_besselj_root(ell, nmax[, only, …])

Find positive zero(s) \(u_{\ell n}\) of the spherical Bessel function \(j_{\ell}\) of the first kind of order \(\ell\), or its derivative \(j'_{\ell}\), up to a maximum number \(n_\textrm{max}\).


harmonia.algorithms.bases.spherical_harmonic(ell, m, theta, phi, conj=False)[source]

Evaluate the spherical harmonic function.

This returns \(Y_{\ell m}(\theta, \phi)\) of degree \(\ell \geqslant 0\) and order \(\vert {m} \vert \leqslant \ell\) at the polar angle \(\theta \in [0, \pi]\) and the azimuthal angle \(\phi \in [0, 2\pi]\).

Parameters
  • ell (int, array_like) – Degree of the spherical harmonic function, ell >= 0.

  • m (int, array_like) – Order of the spherical harmonic function, |m| <= ell.

  • theta (float, array_like) – Polar angle in the interval [0, pi].

  • phi (float, array_like) – Azimuthal angle in the interval [0, 2*pi].

  • conj (bool, optional) – If True (default is False), return the complex conjugate.

Returns

\(Y_{\ell m}\) value at theta and phi.

Return type

complex, array_like

harmonia.algorithms.bases.spherical_besselj(ell, x, derivative=False)[source]

Evaluate the spherical Bessel function of the first kind or its derivative.

This returns \(j_{\ell}(x)\) or \(j'_{\ell}(x)\) of order \(\ell \geqslant 0\).

Parameters
  • ell (int, array_like) – Order of the spherical harmonic function, ell >= 0.

  • x (float, array_like) – Argument of the spherical Bessel function.

  • derivative (bool, optional) – If True (default is False), evaluate the derivative instead.

Returns

\(j_{\ell}\) or \(j'_{\ell}\) value at x.

Return type

float, array_like

harmonia.algorithms.bases.spherical_besselj_root(ell, nmax, only=True, derivative=False)[source]

Find positive zero(s) \(u_{\ell n}\) of the spherical Bessel function \(j_{\ell}\) of the first kind of order \(\ell\), or its derivative \(j'_{\ell}\), up to a maximum number \(n_\textrm{max}\).

Solving for roots of the spherical Bessel function relies on the identity \(j_\ell(x) = \sqrt{\pi/(2x)} J_{\ell + 1/2}(x)\), where \(J_\ell(x)\) is the Bessel funcion of the first kind. Solving for roots of the derivative function employs the interval bisection method, with the interval ansatz \(\ell + 1 \leqslant x \leqslant n_\textrm{max} \operatorname{max}\{4, \ell\}\).

Parameters
  • ell (int) – Order of the spherical Bessel function, ell >= 0.

  • nmax (int) – Maximum number of positive zeros to be found, nmax >= 1.

  • only (bool, optional) – If True (default), return the maximal root only.

  • derivative (bool, optional) – If True (default is False), compute the zero(s) of the derivative function instead.

Returns

u_ell – Positive zero(s) for order ell (in ascending order).

Return type

float, array_like

Spectrum discretisation (discretisation)

Discretise the Fourier spectrum under boundary conditions.

DiscreteSpectrum(radius, condition, highcut)

Discrete Fourier spectrum under radial boundary conditions.


class harmonia.algorithms.discretisation.DiscreteSpectrum(radius, condition, highcut, maxdeg=None, lowcut=0.0, mindeg=0, comm=None)[source]

Discrete Fourier spectrum under radial boundary conditions.

The spectral modes are indexed by tuple \((\ell, n)\), where \(\ell\) is the spherical degree associated with the spherical harmonic and Bessel functions, and \(n\) is the spherical depth associated with roots of the spherical Bessel function.

The discrete wavenumbers are

\[k_{\ell n} = \frac{u_{\ell n}}{R} \,, \quad \text{where} \quad \{ u_{\ell n}: n = 1, \dots, n_{\textrm{max},\ell} \}_{\ell}\]

are zeros of the spherical Bessel functions of order \(\ell\) if the boundary condition is Dirichlet, or zeros of their derivatives if the boundary condition is Neumann. The maximum spherical depth \(n_{\textrm{max},\ell}\) corresponds to the largest wavenumber allowed in the cutoff range.

The normalisation coefficients derived from completeness relations are

\[\begin{split}\kappa_{\ell n} = \begin{cases} (2/R^3) j_{\ell+1}^{-2}(u_{\ell n}) \,, \quad \text{for Dirichlet boundary conditions;} \\ (2/R^3) j_{\ell}^{-2}(u_{\ell n}) \left[ 1 - \ell (\ell + 1) / u_{\ell n}^2 \right]^{-1} \,, \quad \text{for Neumann boundary conditions.} \end{cases}\end{split}\]
Parameters
  • radius (float) – Boundary radius.

  • condition ({‘dirichlet’, ‘neumann’}) – Either Dirichlet or Neumann boundary condition.

  • highcut (float) – Upper cutoff wavenumber of the spectrum.

  • maxdeg (int or None, optional) – Maximum spherical degree (default is None).

  • lowcut (float, optional) – Lower cutoff wavenumber of the spectrum (default is 0.).

  • mindeg (int, optional) – Minimum spherical degree (default is 0).

  • comm (mpi4py.MPI.Comm) – MPI communicator.

Variables
  • degrees (list of int) – Spherical degrees associated with the discrete spectrum.

  • depths (list of int) – Spherical depths associated with each spherical degree.

  • mode_counts (list of int) – Total number of spectral modes associated with each spherical degree counting spherical order multiplicities.

  • roots (dict) – Spherical Bessel roots associated with the discrete spectrum as a dictionary accessed by integer spherical degrees.

  • attrs (dict) – Discrete spectrum attributes including the following keys: ‘min_wavenumber’ and ‘max_wavenumber’ for minimum and maximum wavenumbers, ‘boundary_radius’ and ‘bounded_volume’ for the boundary radius and bounded volume, and ‘boundary_condition’ for the imposed boundary condition.

property wavenumbers

Discrete mode wavenumbers \(k_{\ell n}\).

Returns

Mode wavenumbers as a dictionary accessed by doublet tuples \((\ell, n)\).

Return type

dict

property normalisations

Normalisation coefficients \(\kappa_{\ell n}\).

Returns

Normalisation coefficients as a dictionary accessed by doublet tuples \((\ell, n)\).

Return type

dict

Numerical integration (integration)

Integrate numerically in specified coordinate systems.

Warning

Quadrature integration of spherical functions may suffer from poor convergence.

Spherical integrals

angular_integral(angular_func)

Compute the full spherical angular integral.

radial_integral(radial_func, rmax)

Compute the radial integral up to the given maximum radius.

pixelated_angular_integral(angular_func, nside)

Compute the full spherical angular integral with pixelation.


harmonia.algorithms.integration.angular_integral(angular_func)[source]

Compute the full spherical angular integral.

Notes

Arguments \((\theta, \phi)\) of angular_func must be in radians in the domain \([0, \pi] \times [0, 2\pi]\).

Parameters

angular_func (callable) – Angular function to be integrated.

Returns

Angular integral value.

Return type

complex

harmonia.algorithms.integration.radial_integral(radial_func, rmax)[source]

Compute the radial integral up to the given maximum radius.

Parameters
  • radial_func (callable) – Radial function to be integrated.

  • rmax (float) – Upper radial limit rmax > 0.

Returns

integral – Radial integral value.

Return type

float

harmonia.algorithms.integration.pixelated_angular_integral(angular_func, nside)[source]

Compute the full spherical angular integral with pixelation.

Notes

Arguments \((\theta, \phi)\) of angular_func must be in radians in the domain \([0, \pi] \times [0, 2\pi]\).

Parameters
  • angular_func (callable) – Angular function to be integrated.

  • nside (int) – ‘NSIDE’ parameter for healpy pixelation.

Returns

Angular integral value.

Return type

complex

Cosmological modelling (cosmology)

Provide general and fiducial cosmological models and compute derived model quantities.

Cosmological models are implemented by nbodykit (see documentation).

BaseModel(*args, **kwargs)

Fixed (fiducial) cosmological model built with a parameter dictionary read from an external file, e.g.


Background geometry (geometry)

Compute geometrical quantities in a cosmological background.

redshift_from_distance(cosmo[, lg_num_sample])

Invert redshift-to-distance relationship of a cosmological model to redshift-from-distance.

differential_AP_distortion(fiducial_z_to_r, …)

Compute the differential Alcock–Paczynski distortion between a fiducial and a cosmological redshift-to-distance conversion as a fuction of redshift.


harmonia.cosmology.geometry.redshift_from_distance(cosmo, lg_num_sample=5)[source]

Invert redshift-to-distance relationship of a cosmological model to redshift-from-distance.

Notes

This is useful when the Alcock–Paczynski effect needs to be included in modelling. Only valid for redshift between 1.e-3 and 100.

Parameters
  • cosmo (nbodykit.cosmology.cosmology.Cosmology) – Cosmological model.

  • lg_num_sample (float, optional) – Base-10 logarithm of the number of redshift points to sample the comoving distance as a function of redshift (default is 5, i.e. 100000 samle points).

Returns

Redshift-from-distance function.

Return type

callable

harmonia.cosmology.geometry.differential_AP_distortion(fiducial_z_to_r, variable_z_to_r, max_redshift=10.0, lg_num_sample=5)[source]

Compute the differential Alcock–Paczynski distortion between a fiducial and a cosmological redshift-to-distance conversion as a fuction of redshift.

Parameters
  • fiducial_z_to_r (callable) – Fiducial redshift-to-distance conversion.

  • variable_z_to_r (callable) – Variable redshift-to-distance conversion.

  • max_redshift (float, optional) – Maximum redshift to sample for interpolation (default is 10).

  • lg_num_sample (float, optional) – Base-10 logarithm of the number of redshift points to sample the differential distortion as a function of redshift (default is 5, i.e. 100000 samle points).

Returns

Differential distortion as a fuction of redshift.

Return type

callable

Scale dependence (scale_dependence)

Compute scale-dependence modifications to galaxy clustering from local primordial non-Gausianity \(f_\textrm{NL}\).

The scale-independent linear bias is modified as

\[b_1(z) \mapsto b_1(z) + f_\textrm{NL} [b_1(z) - p] A(k,z) \,,\]

where \(p\) is a tracer-dependent parameter and the scale-dependence modification kernel

\[A(k,z) = 1.27 \left( \frac{H_0}{\mathrm{c}} \right)^2 \frac{3\Omega_\textrm{m,0} \delta_\textrm{c}}{k^2 D(z) T(k)}\]

relates to the cosmological model and its growth factor \(D(z)\) (normalised to unity at the current epoch; hence the numerical factor 1.27) and transfer function \(T(k)\) (normalised to unity as \(k \to 0\)).

scale_dependence_modification(cosmo, redshift)

Return the scale-dependence modification kernel \(A(k,z)\) for a given cosmological model.

scale_dependent_bias(b_1, f_nl, cosmo[, …])

Return the scale-dependent bias modulated by local primordial non-Gaussianity.

modified_power_spectrum(b_1, f_nl, cosmo[, …])

Return the tracer power spectrum modified by primordial non-Gaussianity.


harmonia.cosmology.scale_dependence.scale_dependence_modification(cosmo, redshift)[source]

Return the scale-dependence modification kernel \(A(k,z)\) for a given cosmological model.

Parameters
Returns

Scale-dependence modification kernel as a function of wavenumber (in \(h\)/Mpc).

Return type

callable

harmonia.cosmology.scale_dependence.scale_dependent_bias(b_1, f_nl, cosmo, redshift=0.0, tracer_p=1.0)[source]

Return the scale-dependent bias modulated by local primordial non-Gaussianity.

Parameters
  • b_1 (float) – Scale-independent linear bias at input redshift.

  • f_nl (float) – Local primordial on-Gaussianity.

  • cosmo (nbodykit.cosmology.cosmology.Cosmology) – Cosmological model.

  • redshift (float, optional) – Redshift at which quantities are evaluated (default is 0.).

  • tracer_p (float, optional) – Tracer-dependent parameter \(p\) (default is 1.).

Returns

Scale-dependent bias as a function of wavenumber (in \(h\)/Mpc).

Return type

callable

harmonia.cosmology.scale_dependence.modified_power_spectrum(b_1, f_nl, cosmo, redshift=0.0, tracer_p=1.0, nbar=None, contrast=None)[source]

Return the tracer power spectrum modified by primordial non-Gaussianity.

Parameters
  • b_1 (float) – Scale-independent linear bias at input redshift.

  • f_nl (float) – Local primordial non-Gaussianity.

  • cosmo (nbodykit.cosmology.cosmology.Cosmology) – Cosmological model.

  • redshift (float, optional) – Redshift at which quantities are evaluated (default is 0.).

  • tracer_p (float, optional) – Tracer-dependent parameter \(p\) (default is 1.).

  • nbar (float or None, optional) – If not None (default), add 1/nbar as shot noise to the resulting power spectrum.

  • contrast (float or None, optional) – If not None (default), add additional 1/(contrast*nbar) as shot noise to the resulting power spectrum.

Returns

Tracer power spectrum modified by primordial non-Gaussianity as a function of wavenumber (in \(h\)/Mpc).

Return type

callable

Catalogue mapper (mapper)

Build and process discrete catalogues into Fourier-space maps.

Note

Unless otherwise specified, the length dimension in the module is in units of Mpc/\(h\).

Catalogue maker (catalogue_maker)

Make discrete catalogues from observed or simulated realisations.

SourceCatalogue(*args, **kwargs)

Catalogue from external sources.

RandomCatalogue(*args, **kwargs)

Uniform random catalogue.

SphericalFKPCatalogue(radius[, …])

FKP-style paired catalogues in a spherical domain.

spherical_indicator(cartesian_position, …)

Indicate whether an object lies within the a spherical domain centred at the origin.


harmonia.mapper.catalogue_maker.spherical_indicator(cartesian_position, bounding_radius)[source]

Indicate whether an object lies within the a spherical domain centred at the origin.

Parameters
  • cartesian_position (float, array_like) – Object position in Cartesian coordinates.

  • bounding_radius (float) – Bounding radius of the spherical domain.

Returns

True if the objection position lies within the spherical domain.

Return type

bool numpy.ndarray

class harmonia.mapper.catalogue_maker.SphericalFKPCatalogue(radius, data_catalogue=None, source_file=None, source_kwargs=None, random_catalogue=None, contrast=None, mask=None, selection=None, weight=None, random_seed=None, apply_selection_as_veto=False)[source]

FKP-style paired catalogues in a spherical domain.

Parameters
  • data_catalogue (SourceCatalogue or None, optional) – Data catalogue. Cannot be None (default) unless source_file and source_kwargs are provided as an alternative to read the catalogue source file.

  • source_file (str or None, optional) – Catalogue source file path. This can only be None (default) if data_catalogue is already provided.

  • source_kwargs (dict or None, optional) – Parameters to pass to SourceCatalogue. This can only be None (default) if data_catalogue is already provided.

  • random_catalogue (nbodykit.base.catalog.CatalogSource or None, optional) – Random catalogue (default is None).

  • contrast (float or None, optional) – Mean density contrast compared to the data catalogue used to generate a random catalogue (default is None). Ignored if random_catalogue is provided.

  • mask (callable or None, optional) – Any veto mask function to be applied to both the data and random catalogues. Must be a function of 3-d Cartesian coordinates only, assuming the origin is at the centre of the catalogues.

  • selection (callable or None, optional) – Any selection function (normalised to unity) to be applied to both the data and random catalogues. Must be a function of 3-d Cartesian coordinates only, assuming the origin is at the centre of the catalogues.

  • weight (callable or None, optional) – Any weight function to be applied to both the data and random catalogues. Must be a function of 3-d Cartesian coordinates only, assuming the origin is at the centre of the catalogues.

  • random_seed (int or None, optional) – Random seed of the random catalogue (default is None).

  • apply_selection_as_veto (bool, optional) – If True (default is False), selection returning boolean values is applied to the ‘Selection’ column of the catalogues instead of the ‘NZ’ column (see nbodykit documentation for this peculiarity.)

Variables

Notes

If selection is provided, it must be normalised to the overall mean number density of the data catalogue.

Map transform (map_transform)

Transform discrete catalogues to Fourier-space map data.

SphericalMap(catalogues, disc[, initialise])

Discretised spherical Fourier map from catalogue sources.

CartesianMap(catalogues, orders[, kmin, …])

Compressed Cartesian Fourier-space map from catalogue sources.


class harmonia.mapper.map_transform.SphericalMap(catalogues, disc, initialise=True)[source]

Discretised spherical Fourier map from catalogue sources.

Parameters
  • catalogues (SphericalFKPCatalogue) – FKP-style paired catalogues in a spherical domain.

  • disc (DiscreteSpectrum) – Discrete spectrum for the map.

  • initialise (bool, optional) – If True (default), map transform is performed upon creation.

Variables
  • catalogues (SphericalFKPCatalogue) – FKP-style paired catalogues associated with the map.

  • disc (DiscreteSpectrum) – Discrete spectrum associated with the map.

  • attrs (dict) – Attributes inherited upon creation.

property density_contrast

Spherical Fourier coefficients for the density contrast between data and random catalogues.

Notes

When the spherical map is initialsed upon creation or this is directly or indirectly accessed without initialising the spherical map upon creation, discrete spherical Fourier transform is performed by direct summation over all selected objects in the paired catalogues.

Calling this method stores the transformed number densities of data and random catalogues internally for further processing, e.g. for mode_power(). Computational redundancy is reduced by employing parity relations between spherical harmonics of the same degree but opposite orders.

Returns

Density contrast coefficients of the spherical map.

Return type

SphericalArray

property mode_power

Spherical Fourier mode power suitably normalised.

In the simplest case of a full-sky statistically isotropic map without any mode coupling, this is equivalent to the Cartesian power spectrum at the same wavenumbers.

Returns

Spherical mode power at discrete mode wavenumbers with mode counts and indices.

Return type

dict

class harmonia.mapper.map_transform.CartesianMap(catalogues, orders, kmin=None, kmax=None, dk=None, num_mesh=256, resampler='tsc', interlaced=True)[source]

Compressed Cartesian Fourier-space map from catalogue sources.

Parameters
  • catalogues (SphericalFKPCatalogue or FKPCatalog) – FKP-style paired catalogues in a spherical domain.

  • orders (sequence of int) – Orders of the power spectrum multipoles.

  • kmin (float, optional) – Minimum wavenumber of the compressed map (in \(h\)/Mpc) (default is None).

  • kmax (float or None, optional) – Maximum wavenumber of the compressed map (in \(h\)/Mpc) (default is None).

  • dk (float or None, optional) – Wavenumber bin width (in \(h\)/Mpc) (default is None).

  • num_mesh (int, optional) – Mesh number per dimension for interpolating the discrete catalogues on a grid.

  • resampler ({‘cic’, ‘tsc’, ‘pcs’}, optional) – Grid assignment scheme (default is ‘tsc’) for catalogue interpolation.

  • interlaced (bool, optional) – If True (default), use interlacing for aliasing mitigation.

Variables

Map reader (reader)

Model Fourier-space map statistics and construct cosmological likelihoods.

Note

Unless otherwise specified, the length dimension in the module is in units of Mpc/\(h\).

Spherical couplings (couplings)

Compute spherical Fourier couplings coefficients for given survey spefications and cosmological models.

SphericalCoefficientWarning

Warning issued for poorly determined spherical coefficient.

Couplings(disc[, survey_specs, cosmo_specs, …])

Angular, radial and RSD coupling coefficients computed for given survey and cosmological specifications.

Kernels

Coupling kernels are integrands without the coordinate Jacobian, which may include the following factors:

  • angular mask \(M(\hat{\mathbf{r}})\);

  • radial selection \(\phi(r)\);

  • radial weight \(w(r)\) and its derivative \(w'(r)\);

  • clustering evolution, \(D(z)\), which is the linear growth factor normalised to unity at the \(z = 0\) epoch;

  • linear bias + clustering evolution \(G(z, k) = b(z, k) D(z) / b(z_*, k) D(z_*)\) normalised to unity at a fiducial epoch \(z_*\), where \(b(z, k)\) is the scale-dependent linear bias;

  • linear growth rate + clustering evolution \(F(z) = f(z) D(z) / f(z_*) D(z_*)\) normalised to unity at a fiducial epoch \(z_*\), where \(f(z)\) is the linear growth rate;

  • differential Alcock–Paczynski (AP) distortion

    \[\gamma(z) = \frac{ \operatorname{d}\!\breve{\chi}(z) }{ \operatorname{d}\!\chi(z) } \,,\]

    where \(\breve{r} = \breve{\chi}(z)\) is the distance converted from redshift in a fiducial cosmology rather than from the true comoving distance–redshift correspondence \(z = \chi^{-1}(r)\).

Couplings

Coupling coefficients are computed by integrating the angular, radial and RSD coupling kernels

\[\begin{split}\begin{align*} M_{\mu\nu} &= \int \operatorname{d}^2\!\hat{\mathbf{r}} \, Y_{\ell_\mu m_\mu}^*(\hat{\mathbf{r}}) M(\hat{\mathbf{r}}) Y_{\ell_\nu m_\nu}(\hat{\mathbf{r}}) \,, \\ \Phi_{\mu\nu} &= \kappa_{\ell_\nu n_\nu} \int \operatorname{d}\!r \, r^2 w(\breve{r}) j_{\ell_\mu}(k_{\ell_\mu n_\mu} \breve{r}) j_{\ell_\nu}(k_{\ell_\nu n_\nu} r) G(z(r), k_{\ell_\nu n_\nu}) \phi(r) \,, \\ \Upsilon_{\mu\nu} &= \frac{\kappa_{\ell_\nu n_\nu}}{k_{\ell_\nu n_\nu}} \int \operatorname{d}\!r \, r^2 \frac{\operatorname{d}\!}{\operatorname{d}\!\breve{r}} \left[ w(\breve{r}) j_{\ell_\mu}(k_{\ell_\mu n_\mu} \breve{r}) \right] j'_{\ell_\nu}(k_{\ell_\nu n_\nu} r) \gamma(z(r)) F(z(r)) \phi(r) \,, \end{align*}\end{split}\]

over the spherical volume element, where \(k_{\ell n}\) is the discrete wavenumber.

When there is no angular masking (i.e. \(M(\hat{\mathbf{r}})\) is constant), the coupling coefficients reduce to \(M_{\mu\nu} = \delta_{\mu\nu}\); if in addition radial selection, weighting and evolutionary effects are all absent and there is no AP correction, then \(M_{\mu\nu} \Phi_{\mu\nu} = \delta_{\mu\nu}\).


exception harmonia.reader.couplings.SphericalCoefficientWarning[source]

Warning issued for poorly determined spherical coefficient.

class harmonia.reader.couplings.Couplings(disc, survey_specs=None, cosmo_specs=None, initialise=True, external_angular_couplings=None, pixelate=None, comm=None)[source]

Angular, radial and RSD coupling coefficients computed for given survey and cosmological specifications.

Notes

Survey specification functions must be given in spherical coordinates and may include the following:

  • ‘mask’ for angular mask;

  • ‘selection’, ‘weight’ and ‘weight_derivative’ for radial selection, weighting and weight derivative.

Cosmological specification functions must be given in redshift and/or the radial coordinate and may include the following:

  • ‘z_from_r’ and ‘chi_of_z’ for cosmological comoving distance-to-redshift conversion and fiducial redshift-to-comoving distance conversion;

  • ‘bias_evolution’, ‘growth_evolution’ for linear bias and linear growth rate, normalised to unity at a fiducial epoch;

  • ‘clustering_evolution’ for clustering evolution as a function of the redshift normalised to unity at the current epoch;

  • ‘differential_AP_distortion’ for differential AP distortion as a function of the redshift.

Parameters
  • disc (DiscreteSpectrum) – Discrete spectrum associated with the couplings.

  • survey_specs (dict{str: callable or None} or None, optional) – Survey specification functions as detailed above.

  • cosmo_specs (dict{str: callable or None} or None, optional) – Cosmological specification functions as detailed above.

  • initialise (bool, optional) – If True, compile all coupling coefficients upon creation.

  • external_angular_couplings (dict{tuple(tuple, tuple): complex} or None, optional) – Pre-compute angular couplings (default is None).

  • pixelate (int or None, optional) – If not None (default), this sets the ‘NSIDE’ parameter of healpy pixelation for angular coupling integration.

  • comm (mpi4py.MPI.Comm or None, optional) – MPI communicator (default is None).

Variables
  • disc (DiscreteSpectrum) – Discrete spectrum associated with the couplings.

  • couplings (dict{str: dict}) – Directory for all coupling coefficients of different coupling types.

__getitem__(key)[source]

Access coupling coefficient by key.

Notes

Only accessible if initialised with compile_couplings(). The access key is a tuple specifying the coupling type and the pair of triplet/double indices, e.g. ['angular', (0, 0), (3, -1, 1)] for \(M_{0,0,3,-1}\).

Parameters

key (tuple) – Coefficient access key. Must contain elements of types specified above.

Returns

Coupling coefficient.

Return type

float or complex

Raises

AttributeError – If the coupling coefficients have not been initialised.

load_angular_couplings(angular_couplings)[source]

Load pre-computed angular coupling coefficients which are independent of the cosmological model.

Parameters

angular_couplings (dict{tuple(tuple, tuple): complex}) – Pre-compute angular couplings.

Raises

ValueError – If the number of entries in angular_couplings do not match the class instance.

compile_couplings(pixelate=None)[source]

Compile all coupling coefficients.

Parameters

pixelate (int or None, optional) – If not None (default), this sets the ‘NSIDE’ parameter of healpy pixelation for angular coupling integration.

compute_coefficient(mu, nu, coupling_type, pixelate=None)[source]

Compute coupling coefficients for given triplet indices.

Parameters
  • mu, nu (tuple(int, int) or tuple(int, int, int)) – Coefficient triplet or reduced doublet index.

  • coupling_type ({‘angular’, ‘radial’, ‘RSD’}) – Coupling type.

  • pixelate (int or None, optional) – If not None (default), this sets the ‘NSIDE’ parameter of healpy pixelation for computing angular coupling coefficients.

Returns

cofficient – Coupling coefficient of the specified type.

Return type

float or complex numpy.ndarray

save(output_file)[source]

Save compiled couplings as a .npz file.

Parameters

output_file (str or pathlib.Path) – Output file path.

classmethod load(input_file, comm=None)[source]

Load compiled couplings from a .npz file.

Parameters

input_file (str or pathlib.Path) – Input file path.

Cosmological likelihoods (likelihoods)

Construct likelihoods from Spherical Fourier coefficients and Cartesian power spectrum multipoles of random fields for cosmological parameter inference.

Probability distributions

complex_normal_pdf(data_vector, …[, …])

Compute the complex normal probability density function or its natural logarithm given the zero-centred data vector and its covariance matrix.

multivariate_normal_pdf(data_vector, …[, …])

Compute the multivariate normal probability density function or its natural logarithm given the data vector, its mean vector and covariance matrix.

modified_student_pdf(data_vector, …[, ret_log])

Compute the multivariate modified Student probability density function or its natural logarithm given the data vector, its mean vector and covariance matrix.

Likelihood inference

spherical_covariance(pivot, spherical_model, …)

Compute the parametrised covariance matrix of spherical Fourier coefficients.

cartesian_moments(pivot, orders, cartesian_model)

Compute the parametrised mean and covariance of Cartesian power spectrum multipoles.

LogLikelihood([spherical_data, …])

Construct the logarithmic likelihood function from cosmological data.


harmonia.reader.likelihoods.chi_square(data_vector, covariance_matrix)[source]

Calculate chi-square from zero-centred data vector and its covariance matrix.

Parameters
  • data_vector (complex, array_like) – Zero-centred data vector.

  • covariance_matrix (complex, array_like) – Covariance matrix.

Returns

chi_sq – Chi-square value.

Return type

float numpy.ndarray

harmonia.reader.likelihoods.complex_normal_pdf(data_vector, covariance_matrix, ret_log=True, downscale=None)[source]

Compute the complex normal probability density function or its natural logarithm given the zero-centred data vector and its covariance matrix.

Parameters
  • data_vector (complex, array_like) – Zero-centred data vector.

  • covariance_matrix (complex, array_like) – Covariance matrix.

  • ret_log (bool, optional) – If True (default), return logarithmic probability density.

  • downscale (float or None, optional) – If not None (default), the data vector and covariance matrix are simultaneous downscaled to avoid numerical issue.

Returns

(Logarithmic) probability density.

Return type

float, array_like

harmonia.reader.likelihoods.multivariate_normal_pdf(data_vector, expectation_vector, covariance_matrix, ret_log=True)[source]

Compute the multivariate normal probability density function or its natural logarithm given the data vector, its mean vector and covariance matrix.

Parameters
  • data_vector (float, array_like) – Data vector.

  • expectation_vector (float, array_like) – Mean vector.

  • covariance_matrix (float, array_like) – Covariance matrix.

  • ret_log (bool, optional) – If True (default), return logarithmic probability density.

Returns

density – Logarithmic probability density value.

Return type

float, array_like

harmonia.reader.likelihoods.modified_student_pdf(data_vector, expectation_vector, covariance_matrix, degree, ret_log=True)[source]

Compute the multivariate modified Student probability density function or its natural logarithm given the data vector, its mean vector and covariance matrix.

Parameters
  • data_vector (float, array_like) – Data vector.

  • expectation_vector (float, array_like) – Mean vector.

  • covariance_matrix (float, array_like) – Covariance matrix.

  • degree (int) – The degree number. This could be the number of empirical covariance matrices used to obtain the estimated covariance_matrix.

  • ret_log (bool, optional) – If True (default), return logarithmic probability density.

Returns

(Logarithmic) probability density value.

Return type

float, array_like

exception harmonia.reader.likelihoods.LikelihoodWarning[source]

Likelihood evaluation warning.

harmonia.reader.likelihoods.spherical_covariance(pivot, spherical_model, **kwargs)[source]

Compute the parametrised covariance matrix of spherical Fourier coefficients.

Parameters
  • pivot ({‘natural’, ‘spectral’}) – Pivot order for vectorisation.

  • spherical_model (SphericalCorrelator) – Spherical correlator base model.

  • **kwargs – Parameters (other than pivot) to be passed to correlator_matrix() of spherical_correlator.

Returns

covariance_matrix – Covariance matrix.

Return type

complex numpy.ndarray

harmonia.reader.likelihoods.cartesian_moments(pivot, orders, cartesian_model, covariance_estimator=None, mode_counts=None, **kwargs)[source]

Compute the parametrised mean and covariance of Cartesian power spectrum multipoles.

Parameters
  • pivot ({‘order’, ‘wavenumber’}) – Pivot order for vectorisation.

  • orders (list of int) – Orders of the power spectrum multipoles.

  • cartesian_model (CartesianMultipoles) – Cartesian power multipoles base model.

  • covariance_estimator (CovarianceEstimator or None, optional) – Cartesian power multipole covariance estimator. Its wavenumbers must match wavenumbers associated with cartesian_model. If None, no correlation between power spectrum multipoles is assumed but mode_counts must be provided for calculating the power spectrum variance.

  • mode_counts (int, array_like or None, optional) – Number of independent modes for each power spectrum measurement (default is None) used to calculate the power spectrum variance. Ignored if covariance_estimator is provided.

  • **kwargs – Parameters (other than orders) to be passed to convolved_power_multipoles() of cartesian_model.

Returns

  • expectation (float numpy.ndarray) – Power spectrum expectation at the specified wavenumbers.

  • covariance (float numpy.ndarray) – Power spectrum variance at the specified wavenumbers.

class harmonia.reader.likelihoods.LogLikelihood(spherical_data=None, cartesian_data=None, covariance_estimator=None, mode_counts=None, base_spherical_model=None, base_cartesian_model=None, spherical_pivot='natural', cartesian_pivot='order', nbar=None, contrast=None, tracer_p=1.0, comm=None)[source]

Construct the logarithmic likelihood function from cosmological data.

Parameters
  • spherical_data (SphericalArray or None, optional) – Spherical Fourier coefficient data (default is None).

  • cartesian_data (CartesianArray or None, optional) – Spherical Fourier coefficient data (default is None).

  • covariance_estimator (CovarianceEstimator or None, optional) – Cartesian multipole covariance estimator (default is None).

  • mode_counts (int, array_like or None, optional) – Number of independent modes for each Cartesian data point (default is None) as an alternative to covariance_estimator. Ignored if covariance_estimator is provided.

  • base_spherical_model (SphericalCorrelator or None, optional) – Baseline spherical correlator model (default is None).

  • base_cartesian_model (CartesianMultipoles or None, optional) – Baseline Cartesian multipole model (default is None).

  • spherical_pivot ({‘natural’, ‘spectral’}, optional) – Pivot order for spherical map data vectorisation (default is ‘natural’).

  • cartesian_pivot ({‘order’, ‘wavenumber’}, optional) – Pivot order for Cartesian map data vectorisation (default is ‘order’).

  • nbar (float or None, optional) – Mean particle number density (in cubic \(h\)/Mpc). If None (default), shot noise is neglected.

  • contrast (float or None, optional) – If not None (default), this adds additional shot noise level 1 / (contrast * nbar) due to a FKP-style random catalogue.

  • tracer_p (float, optional) – Tracer-dependent parameter for bias modulation by f_nl (default is 1.).

  • comm (mpi4py.MPI.Comm or None, optional) – MPI communicator (default is None).

Variables
  • attrs (dict) – Directory holding input parameters not corresponding to any of the following attributes.

  • spherical_data (SphericalArray or None) – Spherical Fourier coefficient data.

  • cartesian_data (CartesianArray or None) – Spherical Fourier coefficient data.

  • covariance_estimator (CovarianceEstimator or None) – Cartesian multipole covariance estimator.

  • base_spherical_model (SphericalCorrelator or None) – Baseline spherical correlator model.

  • base_cartesian_model (CartesianMultipoles or None) – Baseline Cartesian multipole model.

spherical_map_likelihood(b_1, f_nl, exclude_degrees=(), compression_matrix=None, **kwargs)[source]

Evaluate the spherical map logarithmic likelihood.

Parameters
  • b_1 (float) – Scale-independent linear bias.

  • f_nl (float or None) – Local primordial non-Gaussianity.

  • exclude_degrees (tuple of int, optional) – If not empty (default), modes whose spherical degree match one of its elements are removed from the likelihood.

  • compression_matrix (numpy.ndarray or None, optional) – If not None (default), both the data vector and the model covariance matrix are processed for data compression. This must be compatible with exclude_degrees, i.e. it accounts for elements removed from the data vector and covariance matrix by exclude_degrees.

  • **kwargs – Additional parameters to be passed to correlator_matrix() of base_spherical_model.

Returns

log_likelihood – Logarithmic likelihood.

Return type

float

cartesian_map_likelihood(b_1, f_nl, orders=None, num_samples=None, **kwargs)[source]

Evaluate the Cartesian map logarithmic likelihood.

Parameters
  • b_1 (float) – Scale-independent linear bias of the tracer particles.

  • f_nl (float or None) – Local primordial non-Gaussianity.

  • orders (list of int or None, optional) – Orders of the power spectrum multipoles. If None (default), only the monopole is used.

  • num_samples (int or None, optional) – If None (default), the normal distribution is used without correction for covariance estimation uncertainty; otherwise it is passed as degree to modified_student_pdf() for covariance estimation uncertainty correction 1.

  • **kwargs – Additional parameters to be passed to convolved_power_multipoles() of base_cartesian_model.

Returns

log_likelihood – Logarithmic likelihood.

Return type

float

1

Sellentin E. & Heavens A. F., 2016. MNRAS 456(1), L132–L136. [arXiv: 1511.05969]

Two-point correlator models (models)

Compute Fourier-space two-point correlator models.

Cartesian multipoles

The standard Kaiser model of plane-parallel power spectrum multipoles is implemented,

\[\begin{split}\begin{align*} P_0(k, z) &= \left[ b_1(z)^2 + \frac{2}{3} b(z, k) f(z) + \frac{1}{5} f(z)^2 \right] P_\mathrm{m}(k, z) \,, \\ P_2(k, z) &= \left[ \frac{4}{3} b(z, k) f(z) + \frac{4}{7} f(z)^2 \right] P_\mathrm{m}(k, z) \,, \\ P_4(k, z) &= \frac{8}{35} f(z)^2 P_\mathrm{m}(k, z) \,, \end{align*}\end{split}\]

with shot noise, window convolution and the integral constraint taken into account 1 2; here \(f(z)\) is the linear growth rate, \(b(z, k) = b_1(z) + f_\textrm{NL} \Delta b(k, z)\) is the scale-dependent linear bias including the modification \(\Delta b\) due to local primordial non-Gaussianity \(f_\textrm{NL}\) (see scale_dependence), and \(P_{\textrm{m}}\) is the matter power spectrum.

1

Wilson M. J. et al., 2017. MNRAS 464(3), 3121–3130. [arXiv: 1511.07799]

2

Beutler F. et al., 2017. MNRAS 466(2), 2242–2260. [arXiv: 1607.03150]

Spherical correlator

Spherical 2-point correlators are computed from coupling coefficients (see couplings) as a sum of the signal part

\[\left\langle \delta_\mu \delta_\nu \right\rangle = \sum_\sigma M_{\mu\sigma} M^*_{\nu\sigma} \left[ b_*(k_\sigma) \Phi_{\mu\sigma} + f_* \Upsilon_{\mu\sigma} \right] \left[ b_*(k_\sigma) \Phi_{\nu\sigma} + f_* \Upsilon_{\nu\sigma} \right] \kappa_\sigma^{-1} P_{\textrm{m}*}(k_\sigma) \,,\]

and the shot noise part

\[\left\langle \epsilon_\mu \epsilon_\nu \right\rangle = \frac{1}{\bar{n}} M_{\mu\nu} \int \operatorname{d}\!r r^2 (w^2\phi)(r) j_\mu(r) j_\nu(r) \,,\]

where \(*\) denotes quantities computed at the fiducial epoch \(z_*\), \(\kappa\) denotes the normalisation coefficients (see discretisation), \(\bar{n}\) is the mean number density and \(j_\mu(r) \equiv j_{\ell_\mu}(k_{\ell_\mu n_\mu} r)\).

CartesianMultipoles(wavenumbers, redshift[, …])

Power spectrum multipoles predicted for a given cosmological model and survey specifications.

SphericalCorrelator(disc, redshift[, cosmo, …])

Spherical two-point correlator predicted for a given cosmological model and survey specifications.


class harmonia.reader.models.CartesianMultipoles(wavenumbers, redshift, cosmo=None, power_spectrum=None, growth_rate=None, mask_multipoles=None, window_multipoles=None)[source]

Power spectrum multipoles predicted for a given cosmological model and survey specifications.

Parameters
  • wavenumbers (float, array_like) – Wavenumbers at which the model is evaluated.

  • redshift (float) – Redshift at which the model is evaluated.

  • cosmo (nbodykit.cosmology.cosmology.Cosmology or None, optional) – Baseline cosmological model used to produce the transfer function and power spectrum and to compute the linear growth rate. This can be subsequently updated when calling convolved_power_multipoles(). If None (default) and not subsequently updated, primordial non-Gaussianity modifications cannot be computed.

  • power_spectrum (callable or None, optional) – Baseline linear matter power spectrum model at redshift. Ignored when cosmo is provided. This cannot be None (default) unless it is subsequently updated when calling convolved_power_multipoles().

  • growth_rate (float or None, optional) – Baseline linear growth rate at redshift. If None (default), this is set by cosmo if it is provided; otherwise this is set to 0. This can be subsequently updated when calling two_point_correlator() or correlator_matrix().

  • mask_multipoles (numpy.ndarray or None, optional) – Survey mask multipoles given at sampled separations (default is None). Orders and sampled separations must be sorted.

  • window_multipoles (CartesianArray or None, optional) – Survey window multipoles given at sampled wavenumbers (default is None). If provided, integral constraint corrections are applied. Orders and sampled wavenumbers must be sorted.

Variables
  • mask_multipoles (numpy.ndarray or None) – Survey mask multipoles given at sampled separations.

  • window_multipoles (CartesianArray or None) – Survey window multipoles given at sampled wavenumbers.

  • cosmo (nbodykit.cosmology.cosmology.Cosmology or None) – Cosmological model.

  • matter_power_spectrum (callable) – Matter power spectrum at redshift.

  • growth_rate (float) – Linear growth rate at redshift.

  • attrs (dict) – Any other attributes inherited from input parameters.

convolved_power_multipoles(orders, b_1=None, f_nl=None, nbar=None, contrast=None, tracer_p=1.0, update_model_kwargs=None)[source]

Compute the convolved power spectrum multipoles.

Parameters
  • orders (list of int) – Orders of the power spectrum multipoles. Values only allowed from the set {0, 2, 4}.

  • b_1 (float or None, optional) – Scale-independent linear bias at input redshift. If None (default), no tracer bias is assumed relative to the matter power spectrum.

  • f_nl (float or None, optional) – Local primordial non-Gaussianity (default is None).

  • nbar (float or None, optional) – Mean particle number density (in cubic \(h\)/Mpc). If None (default), shot noise is neglected.

  • contrast (float or None, optional) – If not None (default), this adds additional shot noise 1 / (contrast * nbar) from a FKP-style random catalogue.

  • tracer_p (float, optional) – Tracer-dependent parameter for bias modulation by f_nl (default is 1.).

  • update_model_kwargs (dict or None, optional) – Parameters cosmo, power_spectrum and growth_rate passed as keyword arguments to update the baseline cosmological model.

Returns

convolved_power – Convolved power spectrum multipoles.

Return type

CartesianArray

static kaiser_factors(ell, b, f)[source]

Return the standard Kaiser power spectrum multipole as a multipole of the matter power spectrum.

Notes

Linear bias b and growth rate f must be specified at the same redshift.

Parameters
  • ell (int) – Order of the multipole in {0, 2, 4}.

  • b (float, array_like) – Linear bias of the tracer particles.

  • f (float or None) – Linear growth rate.

Returns

Standard Kaiser factor to multiply by the matter power spectrum.

Return type

float, array_like

class harmonia.reader.models.SphericalCorrelator(disc, redshift, cosmo=None, power_spectrum=None, growth_rate=None, couplings=None, coupling_disc=None, survey_specs=None, cosmo_specs=None, ini_shot_noise=True, comm=None)[source]

Spherical two-point correlator predicted for a given cosmological model and survey specifications.

Notes

Summation is performed over spherical modes to compute the correlator model. Owing to mode couplings, sums need to be truncated at higher mode wavenumbers than wavenumbers under consideration to ensure convergence. We suggest the user trial discrete spectra with different upper wavenumber cutoffs for the couplings and use radialised_power() to gauge convergence. In the future, an automated facility for this purpose may be provided.

Parameters
  • disc (DiscreteSpectrum) – Discrete spectrum associated with the correlator model.

  • redshift (float) – Redshift at which the model is evaluated.

  • cosmo (nbodykit.cosmology.cosmology.Cosmology or None, optional) – Baseline cosmological model used to produce the transfer function (and the power spectrum and linear growth rate if these are not set but required in model evaluation). This can be subsequently updated when calling two_point_correlator() or correlator_matrix(). If None (default) and not subsequently updated, primordial non-Gaussianity modifications cannot be computed.

  • power_spectrum (callable or None, optional) – Baseline linear matter power spectrum model at redshift. Ignored when cosmo is provided; otherwise this cannot be None (default) unless it is subsequently updated when calling two_point_correlator() or correlator_matrix().

  • growth_rate (float or None, optional) – Baseline linear growth rate at redshift. If None (default), this is set by cosmo (if provided); otherwise this is set to zero. This can be subsequently updated when calling two_point_correlator() or correlator_matrix().

  • couplings (Couplings or None, optional) – Baseline coupling coefficients consistent with the underlying cosmological model cosmo. If None (default), this is compiled from survey_specs and cosmo_specs if either is provided; otherwise all couplings are assumed to be trivial (i.e. angular and radial couplings are Kronecker deltas). This can be subsequently updated when calling two_point_correlator() or correlator_matrix().

  • coupling_disc (Couplings or None, optional) – Discrete spectrum for the baseline coupling coefficients (default is None). This must be provided if couplings is None and needs to be compiled. To ensure sum convergence in the correlator model, the upper wavenumber cutoff for this spectrum may need to be higher than that of the correlator model, i.e. disc.

  • survey_specs (dict{str: callable or None} or None, optional) – Survey specification functions to be passed as survey_specs to Couplings when couplings are compiled. Also used in shot noise calculations.

  • cosmo_specs (dict{str: callable, bool or None} or None, optional) – Baseline cosmological specification functions to be passed as cosmo_specs to Couplings when couplings are compiled. If not None (default), it must be a dictionary holding keys listed in Couplings: if callable values are passed to the keys 'chi_of_z', 'clustering_evolution', 'growth_evolution' or 'differential_AP_distortion', they should be consistent with the current cosmo; otherwise True can be passed here and their values are then derived from cosmo (which must then be set), or unspecified keys assume None values; also note some keys are linked and values must be simultaneously provided. This can be subsequently updated when calling two_point_correlator() or correlator_matrix().

  • ini_shot_noise (bool, optional) – If True (default), shot noise integrals are evaluated upon initialisation.

  • comm (mpi4py.MPI.Comm or None, optional) – MPI communicator. If None (default), no multiprocessing is invoked.

Variables

See also

Couplings

More details related to couplings and especially the cosmo_specs parameter.

two_point_correlator(mu, nu, b_1, f_nl=None, nbar=None, contrast=None, tracer_p=1.0, update_model_kwargs=None)[source]

Compute two-point correlator for given triplet indices.

Parameters
  • mu, nu (tuple(int, int, int)) – Coefficient triplet index.

  • b_1 (float) – Scale-independent linear bias at input redshift.

  • f_nl (float or None, optional) – Local primordial non-Gaussianity (default is None).

  • nbar (float or None, optional) – Mean particle number density (in cubic \(h\)/Mpc). If None (default), shot noise is neglected.

  • contrast (float or None, optional) – If not None (default), this adds additional shot noise 1 / (contrast * nbar) from a FKP-style random catalogue.

  • tracer_p (float, optional) – Tracer-dependent parameter for bias modulation by f_nl (default is 1.).

  • update_model_kwargs (dict or None, optional) – Parameters cosmo, power_spectrum, growth_rate and cosmo_specs passed as keyword arguments to update the baseline cosmological model. If cosmo_specs is passed (including None values), radial and RSD couplings will be updated.

Returns

Spherical Fourier 2-point function value for given triplet indices. If cosmo_specs is passed, radial and RSD couplings will be updated.

Return type

complex

correlator_matrix(pivot, b_1=None, f_nl=None, nbar=None, contrast=None, tracer_p=1.0, diagonal=False, shot_noise_only=False, update_model_kwargs=None, report_progress=False)[source]

Compute two-point correlator matrix for some vetorisation of all spectrum modes.

Parameters
  • pivot ({‘natural’, ‘spectral’}) – Pivot order for vectorisation.

  • b_1 (float) – Scale-independent linear bias at input redshift.

  • f_nl (float or None, optional) – Local primordial non-Gaussianity (default is None).

  • nbar (float or None, optional) – Mean particle number density (in cubic \(h\)/Mpc). If None (default), shot noise is neglected.

  • contrast (float or None, optional) – If not None (default), this adds additional shot noise 1 / (contrast * nbar) from a FKP-style random catalogue.

  • tracer_p (float, optional) – Tracer-dependent parameter for bias modulation by f_nl (default is 1.).

  • diagonal (bool, optional) – If True (default is False), return only the diagonal matrix part.

  • shot_noise_only (bool, optional) – If True (default is False), return only the shot noise correlator matrix.

  • update_model_kwargs (dict or None, optional) – Parameters cosmo, power_spectrum, growth_rate and cosmo_specs passed as keyword arguments to update the baseline cosmological model. If cosmo_specs is passed (including None values), radial and RSD couplings will be updated.

  • report_progress (bool, optional) – If True (default is False), progress status will be reported.

Returns

Two-point correlator matrix vectorised by the given pivot.

Return type

complex numpy.ndarray

See also

SphericalArray

radialised_power(b_1=None, f_nl=None, nbar=None, contrast=None, tracer_p=1.0, shot_noise_only=False)[source]

Compute the radialised spherical mode power.

Notes

This relies on correlator_matrix() with diagonal=True. Results are only meaningful in the radialisation limit when all coupling coefficients are trivial (i.e. Kronecker deltas) so that spherical modes are mutually independent 3. Mode power is averaged over equivalent spherical orders and suitably normalised so that it matches the Cartesian power spectrum for a full-sky statistically isotripoc map.

Parameters
  • b_1 (float) – Scale-independent linear bias at input redshift.

  • f_nl (float or None, optional) – Local primordial non-Gaussianity (default is None).

  • nbar (float or None, optional) – Mean particle number density (in cubic \(h\)/Mpc). If None (default), shot noise is neglected.

  • contrast (float or None, optional) – If not None (default), this adds additional shot noise 1 / (contrast * nbar) from a FKP-style random catalogue.

  • tracer_p (float, optional) – Tracer-dependent parameter for bias modulation by f_nl (default is 1.).

  • shot_noise_only (bool, optional) – If True (default is False), return only the shot noise power.

Returns

Radialised spherical mode powers at wavenumbers with mode indices.

Return type

dict

See also

mode_power

3

Rassat A. & Refregier A., 2012. A&A 540, A115. [arXiv: 1112.3100]

Survey factory (surveyor)

Produce from survey definition and specifications the quantities needed for map analysis, e.g. mask map, selection function, mask/window function multipoles, covariance estimates and data compressor.

Note

Unless otherwise specified, the length dimension in the module is in units of Mpc/\(h\).

Coordinates (coordinates)

Handle the survey coordinate systems.

cartesian_to_spherical(cartesian_coords)

Convert 3-d Cartesian coordinate arrays to spherical coordinate arrays.

spherical_to_cartesian(spherical_coords)

Convert 3-d spherical coordinate arrays to Cartesian coordinate arrays.

sky_to_spherical(sky_coords[, z_to_r])

Convert 3-d (or 2-d) sky coordinate arrays (Z, DEC, RA) (or (DEC, RA)) to spherical coordinate arrays.

spherical_to_sky(spherical_coords[, z_from_r])

Convert 3-d spherical coordinate arrays to sky coordinate arrays.

to_box_coords(native_coord_system[, …])

Convert a function defined for a native 3-d curvilinear coordinate system to an equivalent function defined for a box in Cartesian coordinates.


harmonia.surveyor.coordinates.cartesian_to_spherical(cartesian_coords)[source]

Convert 3-d Cartesian coordinate arrays to spherical coordinate arrays.

The coordinate transformation is given by

\[r = \sqrt{x^2 + y^2 + z^2} \,, \quad \theta = \arccos(z/r) \,, \quad \phi = \arctan(y/x) \,,\]

where the image of \(\arccos\) is \([0, \pi]\), and \(\arctan\) has an extended image set \([0, 2\pi]\).

Parameters

cartesian_coords (float, array_like) – Cartesian coordinates.

Returns

spherical_coords – Spherical coordinates.

Return type

float numpy.ndarray

harmonia.surveyor.coordinates.spherical_to_cartesian(spherical_coords)[source]

Convert 3-d spherical coordinate arrays to Cartesian coordinate arrays.

The coordinate transformation is given by

\[x = r \sin\theta \cos\phi \,, \quad y = r \sin\theta \sin\phi \,, \quad z = r \cos\theta \,.\]
Parameters

spherical_coords (float, array_like) – Spherical coordinates.

Returns

cartesian_coords – Cartesian coordinates.

Return type

float numpy.ndarray

harmonia.surveyor.coordinates.sky_to_spherical(sky_coords, z_to_r=None)[source]

Convert 3-d (or 2-d) sky coordinate arrays (Z, DEC, RA) (or (DEC, RA)) to spherical coordinate arrays.

The spherical surface coordinate transformation is given by

\[\theta = \frac{\pi}{180} (90 - \delta) \,, \quad \phi = \frac{\pi}{180} \alpha \,.\]

where \(\delta\) is the declination (DEC) and \(\alpha\) the right ascension (RA) both given in degrees.

The radial coordinate transformation from redshifts is given by cosmo.comoving_distance() if the input coordinates are 3-d and a redshift-to-distance conversion function is provided.

Parameters
  • sky_coords (float, array_like) – Sky coordinates (2-d or 3-d).

  • z_to_r (callable or None, optional) – Redshift-to-distance conversion (default is None).

Returns

spherical_coords – Spherical surface coordinates.

Return type

float numpy.ndarray

harmonia.surveyor.coordinates.spherical_to_sky(spherical_coords, z_from_r=None)[source]

Convert 3-d spherical coordinate arrays to sky coordinate arrays.

The spherical surface coordinate transformation is given by

\[\delta = 90 - \frac{180}{\pi} \theta \,, \quad \alpha = \frac{180}{\pi} \phi \,.\]

where \(\delta\) is the declination (DEC) and \(\alpha\) the right ascension (RA) both given in degrees.

The radial coordinate transformation is given by a distance-to-redshift conversion if it provided.

Parameters
  • sky_coords (float, array_like) – Sky coordinates.

  • z_from_r (callable or None, optional) – Distance-to-redshift conversion (default is None).

Returns

sky_coords – Sky coordinates.

Return type

float numpy.ndarray

harmonia.surveyor.coordinates.to_box_coords(native_coord_system, box_centre=None, conversion_kwargs=None)[source]

Convert a function defined for a native 3-d curvilinear coordinate system to an equivalent function defined for a box in Cartesian coordinates.

Parameters
  • native_coord_system ({‘null’, ‘spherical’, ‘sky’}, optional) – Native coordinate system of the function. If ‘cartesian’ (with the implicit assumption that the origin is at the box centre), the function is unconverted unless box_shift is also provided.

  • box_centre (float, array_like or Nonw) – If provided, translate the box coordinates so that the (positive) box_centre is moved to(0, 0, 0). This is needed to when the wrapped function needs to accept box coordinates with the origin placed at its corner whilst the native_coord_system has the origin at the centre of the box (e.g. ‘spherical’).

  • conversion_kwargs (dict or None, optional) – Additional parameters to use in conversion, e.g. z_from_r if native_coords is ‘sky’ (default is None).

Returns

The original function now accepting Cartesian coordinates.

Return type

callable

Definition (definition)

Produce survey definition and specifications for processing and analysing cosmological data.

generate_mask_by_sky_fraction(coord_system)

Generate mask function given the sky fraction corresponding to a spherical cap (rather than a wedge).

generate_mask_from_map(coord_system[, …])

Generate mask function from a veto mask map.

generate_selection_by_distribution(…[, …])

Generate selection function based on a probability distribution suitably rescale horizontally and vectically.

generate_selection_from_samples(…)

Generate selection function interpolated from normalised selection function samples at sample coodinates.

generate_selection_samples(…[, …])

Generate selection function samples from redshift samples.


harmonia.surveyor.definition.generate_mask_by_sky_fraction(coord_system, sky_fraction=1.0, box_shift=None)[source]

Generate mask function given the sky fraction corresponding to a spherical cap (rather than a wedge).

Parameters
  • coord_system ({‘cartesian’, ‘spherical’, ‘sky’}) – Coordinate system of the generated mask function. If ‘sky’, the returned function can accept either 2-d or 3-d coordinates.

  • sky_fraction (float, optional) – Sky coverage, 0 < sky_fraction <= 1.

  • box_shift (float, array_like or None, optional) – Passed as box_centre to to_box_coords() so that the resulting mask function accepts Cartesian coordinates with origin at a box corner. Ignored unless coord_system is ‘cartesian’.

Returns

mask_function – A mask function accepting coord_system as its coordinate system and returning boolean values.

Return type

callable

harmonia.surveyor.definition.generate_mask_from_map(coord_system, mask_map=None, nside=None, mask_map_file=None, box_shift=None, ret_nside=False)[source]

Generate mask function from a veto mask map.

Parameters
  • coord_system ({‘cartesian’, ‘spherical’, ‘sky’}) – Coordinate system of the generated mask function. If ‘sky’, the returned function can accept either 2-d or 3-d coordinates.

  • mask_map (numpy.ndarray or None, optional) – A veto array corresponding to pixels in a healpy mask map with parameter nside (default is None). Ignored if a healpy-generated .fits file for the mask map is provided.

  • nside (int or None, optional) – ‘NSIDE’ parameter of the healpy mask map (default is None). Ignored if a healpy-generated .fits file for the mask map is provided.

  • mask_map_file (str or pathlib.Path) – healpy-generated .fits file for the mask map. Must contain the ‘NSIDE’ parameter in its header.

  • box_shift (float, array_like or None, optional) – Passed as box_centre to to_box_coords() so that the resulting mask function accepts Cartesian coordinates with origin at a box corner. Ignored unless coord_system is ‘cartesian’.

  • ret_nside (bool, optional) – If True (default is False), return the nside parameter from mask_map read in.

Returns

mask_function – A mask function accepting coord_system as its coordinate system and returning boolean values.

Return type

callable

harmonia.surveyor.definition.generate_selection_by_distribution(coord_scale, distribution, peak, scale=None, location=None, shape=None)[source]

Generate selection function based on a probability distribution suitably rescale horizontally and vectically.

Notes

If the generated selection function detects 3-d input coordinate arrays, it would assume the coordinates are the Cartesian; otherwise it assumes the coordinates are the radial coordinate.

Parameters
  • coord_scale (float) – Coordinate scale of the selection function, e.g. the maximum comoving radius of a catalogue.

  • peak (float) – Peak value of the selection function used for renormalising the probability density function by its maximum value.

  • distribution ({‘gaussian’, ‘gamma’}) – Distribution to use, either ‘gaussian’ or ‘gamma’.

  • scale, location, shape (float or None, optional) – Scale, location and shape parameter of the distribution. For a normal distribution, scale is the standard deviation and location is the mean in units of coord_scale; for a gamma distribution, see scipy.stats.gamma.

Returns

Normalised selection function of the same coordinate as the sample coordinates.

Return type

callable

harmonia.surveyor.definition.generate_selection_by_cut(low_end, high_end)[source]

Generate selection function at a constant value by a cut.

Notes

If the generated selection function detects 3-d input coordinate arrays, it would assume the coordinates are the Cartesian; otherwise it assumes the coordinates are the radial coordinate.

Parameters

low_end, high_end (float) – Low and high end of the cut in some selection coordinate.

Returns

Constant cut selection function.

Return type

callable

harmonia.surveyor.definition.generate_selection_samples(selection_coordinate, coord_scale, redshift_samples, cosmo, sky_fraction=1.0, bins=None)[source]

Generate selection function samples from redshift samples.

The selection function is normalised by the mean number density inferred from the samples and the sky fraction.

Parameters
  • selection_coordinate ({‘z’, ‘r’}) – Coordinate of the generated selection function, either redshift ‘z’ or comoving distance ‘r’.

  • coord_scale (float) – Coordinate scale used to rescale the selection coordinates, e.g. the maximum comoving radius of the catalogue to which the selection function is applied.

  • redshift_samples (float, array_like) – Redshift samples.

  • cosmo (nbodykit.cosmology.cosmology.Cosmology) – Cosmological model providing the redshift-to-distance conversion and thus normalisation of the selection function by volume.

  • sky_fraction (float, optional) – The sky coverage of redshift samples as a fraction used to normalise the selection function (default is 1.), 0 < sky_fraction <= 1.

  • bins (int, str or None, optional) – Number of bins or binning scheme used to generate the selection function samples (see numpy.histogram_bin_edges()). If None (default), Scott’s rule for binning is used.

Returns

samples, coords – Normalised (dimensionless) selection function samples and corresponding coordinates.

Return type

numpy.ndarray

harmonia.surveyor.definition.generate_selection_from_samples(sample_selections, sample_coords)[source]

Generate selection function interpolated from normalised selection function samples at sample coodinates.

Notes

If the generated selection function detects 3-d input coordinate arrays, it would assume the coordinates are the Cartesian; otherwise it assumes the coordinates are the radial coordinate.

Parameters

sample_selections, sample_coords (float numpy.ndarray) – Selection function samples and coordinates.

Returns

Normalised selection function of the same coordinate as the sample coordinates.

Return type

callable

Synthesis (synthesis)

Systhesise FKP-style paired random catalogues for given survey specifications.

This can be used for determining survey window functions and estimating correlation induced by geometric filtering.

SyntheticCatalogue(number_density, boxsize)

Synthetic catalogue for probing survey geometry.

CovarianceEstimator(realisations[, reference])

Covariance matrix estimator for power spectrum multipoles at fiducial values.

generate_compression_matrix(…[, …])

Generate a compression matrix for spherical modes.


class harmonia.surveyor.synthesis.SyntheticCatalogue(number_density, boxsize, contrast=None, expansion=1.0, sphericalise=None, mask=None, selection=None, weight=None, apply_selection_as_veto=False)[source]

Synthetic catalogue for probing survey geometry.

Parameters
  • mean_density (float) – Mean number density (in cubic \(h\)/Mpc).

  • boxsize (float, array_like) – Catalogue boxsize (in Mpc/\(h\)) as a scalar or a triple of scalars.

  • contrast (float or None, optional) – Contrast of the number density of the secondary random catalogue to that of the primary. If None (default), the secondary catalogue is not produced. This cannot be None if sphericalise is provided.

  • expansion (float, optional) – Expansion factor of the catalogue box for grid assignment (default is 1.). Ignored if sphericalise is set.

  • sphericalise (float or None, optional) – If not None (default), this is passed to the radius parameter of SphericalFKPCatalogue for instantiating the synthetic catalogue.

  • mask, selection, weight (callable or None, optional) – Survey mask, selection or weight function (default is None). Must be given as a function of Cartesian coordinate only assuming the origin is at the centre of the catalogues.

  • apply_selection_as_veto (bool, optional) – If True (default is False), selection returning boolean values is applied to the ‘Selection’ column of the catalogues instead of the ‘NZ’ column (see nbodykit documentation for this peculiarity.)

Variables
compute_spherical_coefficients(disc)[source]

Compute the spherical Fourier density contrast coefficients of the synthetic catalogue.

Parameters

disc (DiscreteSpectrum) – Discrete spectrum.

Returns

Density contrast coefficients.

Return type

SphericalArray

See also

SphericalMap

compute_power(orders, kmin=0.0001, kmax=None, dk=None, num_mesh=256, resampler='tsc', interlaced=True)[source]

Compute the power spectrum multipoles of the synthetic catalogue.

Parameters
  • orders (sequence of int) – Orders of the power spectrum multipoles.

  • kmin (float, optional) – Minimum wavenumber of the map (in \(h\)/Mpc) (default is 1.e-4).

  • kmax (float or None, optional) – Maximum wavenumber of the map (in \(h\)/Mpc) (default is None). If None, kmax is the Nyquist wavenumber determined from num_mesh.

  • dk (float or None, optional) – Wavenumber bin width (in \(h\)/Mpc) (default is None).

  • num_mesh (int, optional) – Mesh number per dimension for interpolating the discrete catalogues on a grid.

  • resampler ({‘cic’, ‘tsc’, ‘pcs’}, optional) – Grid assignment scheme (default is 'tsc') for catalogue interpolation.

  • interlaced (bool, optional) – If True (default), use interlacing for aliasing mitigation.

Returns

Power spectrum multipoles binned in wavenumbers.

Return type

CartesianArray

See also

CartesianMap

compute_window(orders, kmin=1e-06, **kwargs)[source]

Compute the Fourier-space window function multipoles from the synthetic catalogue.

Parameters
  • orders (sequence of int) – Orders of the power spectrum multipoles.

  • kmin (float, optional) – Minimum wavenumber (in \(h\)/Mpc) (default is 1.e-6).

  • **kwargs – Any other parameters to be passed to compute_power().

Returns

Window function multipoles binned in logarithmic wavenumbers.

Return type

CartesianArray

See also

compute_power()

compute_mask(orders, *args, **kwargs)[source]

Compute the configuration-space mask function multipoles from the synthetic catalogue by Hankel transform of the window function multipoles.

Parameters
  • orders (sequence of int) – Orders of the mask function multipoles.

  • *args, **kwargs – Any other position and keyword parameters to be passed to compute_window(). Only required when window_multipoles is not available.

Returns

Mask function multipoles binned in logarithmic separations. This is a NumPy structured array similar to array of CartesianArray.

Return type

numpy.ndarray

See also

CartesianArray

class harmonia.surveyor.synthesis.CovarianceEstimator(realisations, reference=None)[source]

Covariance matrix estimator for power spectrum multipoles at fiducial values.

Parameters
  • realisations (sequence of CartesianArray) – Independent realisations of power spectrum multipoles from which the covariance matrix is estimated. For each realisation, orders and wavenumbers of the multipoles must be sorted.

  • reference (CartesianArray or None, optional) – Underlying power spectrum multipoles binned in wavenumber that are realised for covariance estimation. Orders and wavenumbers of the multipoles must be sorted. If None (default), this is determined by the average of realisations.

Variables
  • realisations (sequence of CartesianArray) – Independent realisations of power spectrum multipoles from which the covariance matrix is estimated.

  • reference (CartesianArray or None) – Underlying power spectrum multipoles binned in wavenumber that are realised for covariance estimation.

  • wavenumbers (float ndarray) – Wavenumbers at which the covariance matrix is estimated. This is set by reference if available; otherwise it is set by the first of realisations.

See also

CartesianArray

save(output_file)[source]

Save the estimator with its attributes as a .npz file.

Parameters

output_file (str or pathlib.Path) – Output file path.

classmethod load(input_file)[source]

Load the estimator with its attributes from a .npz file.

Parameters

input_file (str or pathlib.Path) – Input file path.

get_fiducial_vector(pivot)[source]

Return the fiducial data vector for which the covariance matrix is estimated.

Parameters

pivot ({‘order’, ‘wavenumber’}) – Pivot order used for data vectorisation.

Returns

Fiducial vector.

Return type

float numpy.ndarray

get_fiducial_covariance(pivot)[source]

Return the fiducial covariance matrix estimated from data realisations.

Parameters

pivot ({‘order’, ‘wavenumber’}) – Pivot order used for data vectorisation.

Returns

Estimate fiducialcovariance.

Return type

float numpy.ndarray

harmonia.surveyor.synthesis.generate_compression_matrix(fiducial_model_kwargs, extremal_model_kwargs=None, sensitivity_threshold=0.01, discard=None)[source]

Generate a compression matrix for spherical modes.

Notes

Compression is achieved by discarding non-positive eigenvalue modes that are at least \(10^{-8}\) times smaller than the largest and in addition any of the following means:

  • discard is passed to discard a number of low-eigenvalue modes;

  • extremal_model_kwargs is passed and eigenvalues of the resulting model covariance are compared with those from fiducial_covariance. Modes corresponding to low, insensitive (i.e. relative difference less than sensitivity_threshold) are discarded.

  • A combination of the above if the appropriate parameters are passed.

Parameters
  • fiducial_model_kwargs (dict) – Fiducial model parameters to be passed to spherical_covariance().

  • extremal_model_kwargs (dict or None, optional) – Extremal model parameters to be passed to spherical_covariance().

  • sensitivity_threshold (float, optional) – Sensitivity threshold for modes deemed discardable (default is 0.01).

  • discard (int or None, optional) – Number of low-eigenvalue modes to discard from all modes (default is None).

Returns

compression_matrix – Compression matrix.

Return type

numpy.ndarray

Utilities (utils)

Provide utilities for processing and logging as well as common mathematical function and algorithms.

Processing and monitoring

Progress(task_length[, num_checkpts, …])

Progress status of tasks.

setup_logger()

Return the root logger formatted with elapsed time and piped to stdout.

clean_warning_format(message, category, …)

Clean warning message format.

restore_warnings(captured_warnings)

Emit captured warnings.

mpi_compute(data_array, mapping[, comm, …])

Multiprocess mapping of data with optional progress bar.

Mathematics

const_function(const)

Return a constant function with arbitrary arguments.

normalise_vector(vector_array)

Normalise vector(s).

binary_search(func, a, b[, maxnum, precision])

Binary seach for all zeros of a function in a real interval.

covar_to_corr(covar)

Convert a covariance matrix to its correlation matrix.

mat_logdet(matrix)

Calculate logarithm of the determinant of a matrix.

PositiveDefinitenessWarning

Emit a warning when a matrix is not positive definite.

is_positive_definite(matrix)

Check the positive definiteness of a square matrix by attempting a Cholesky decomposition.


class harmonia.utils.Progress(task_length, num_checkpts=4, process_name=None, logger=None, comm=None, root=0)[source]

Progress status of tasks.

This is an alternative to tqdm for cases where progress is not uniform and only needs to be infrequently reported to a logger. If multiple parallel processes exist, progress status is only reported for the first and last of them.

Parameters
  • task_length (int) – Total number of tasks.

  • num_checkpts (int, optional) – Number of checkpoints for reporting progress (default is 4).

  • process_name (str or None, optional) – If not None (default), this is the process name to be logged.

  • logger (logging.Logger or None, optional) – Logger. If None (default), a print statement is issued.

  • comm (mpi4py.MPI.Comm or None, optional) – MPI communicator (default is None).

  • root (int, optional) – Root process number (default is 0).

Variables
  • process_name (str or None, optional) – If not None (default), this is the process name to be logged.

  • task_length (int) – Total number of tasks.

  • progress_checkpts (float) – Scheduled progress check points, 0 < progress_checkpts <= 1.

  • last_checkpt (int) – Index of the last passed checkpoint, 0 <= last_checkpt <= num_checkpts.

Examples

>>> ntasks = 100
>>> p = Progress(ntasks, process_name='null test')
>>> for task_idx in range(ntasks):
...     p.report(task_idx)
Progress for the single 'null test' process: 25% computed.
Progress for the single 'null test' process: 50% computed.
Progress for the single 'null test' process: 75% computed.
Progress for the single 'null test' process: 100% computed.
report(current_position)[source]

Report the current position in the tasks.

Parameters

current_position (int) – Index of the current position in the tasks (starting from 0).

harmonia.utils.setup_logger()[source]

Return the root logger formatted with elapsed time and piped to stdout.

Returns

logger – Formatted root logger.

Return type

logging.Logger

harmonia.utils.clean_warning_format(message, category, filename, lineno, line=None)[source]

Clean warning message format.

Parameters
  • message, category, filename, lineno (str) – Warning message, warning catagory, origin filename, line number.

  • line (str or None, optional) – Source code line to be included in the warning message (default is None).

Returns

Warning message format.

Return type

str

harmonia.utils.restore_warnings(captured_warnings)[source]

Emit captured warnings.

Parameters

captured_warnings (list of warnings.WarningMessage) – List of recorded warnings as returned by warnings.catch_warnings(record=True).

harmonia.utils.mpi_compute(data_array, mapping, comm=None, root=0, process_name=None, update_rate=None)[source]

Multiprocess mapping of data with optional progress bar.

For each map to be applied, the input data array is scattered over the first axis for computation on difference process, and the computed results are gathered in the exact structure of the input data array on the root process.

Parameters
  • data_array (array_like) – Data array.

  • mapping (callable) – Mapping to be applied.

  • comm (mpi4py.MPI.Comm or None, optional) – MPI communicator. If None, no multiprocessing is performed.

  • root (int, optional) – Rank of the process taken as the root process (default is 0).

  • process_name (str or None) – If None (default), no progress bar is displayed; else this is the process name to be displayed.

  • update_rate (int or None, optional) – If not None (default), this is the progress bar update rate (in inverse seconds). Has no effect if process_name is not provided.

Returns

output_array – Output data processed from mapping. Returns None for process ranks other than root.

Return type

array_like or None

harmonia.utils.const_function(const)[source]

Return a constant function with arbitrary arguments.

Parameters

const (float) – Constant value.

Returns

Constant function.

Return type

callable

Binary seach for all zeros of a function in a real interval.

Parameters
  • func (callable) – Function whose zeros are to be found.

  • a, b (float) – Interval end points, a < b.

  • maxnum (int or None, optional) – Maximum number of zeros needed from below (default is None).

  • precision (float, optional) – Desired absolute precision of the zeros (default is 1.e-8).

Returns

roots – Possible roots.

Return type

float numpy.ndarray or None

Raises

ValueError – If the initial interval covers only one point (a == b).

harmonia.utils.covar_to_corr(covar)[source]

Convert a covariance matrix to its correlation matrix.

Parameters

covar (complex, array_like) – Covariance matrix.

Returns

corr – Correlation matrix.

Return type

numpy.ndarray

harmonia.utils.mat_logdet(matrix)[source]

Calculate logarithm of the determinant of a matrix.

Parameters

matrix (float or complex, array_like) – Matrix.

Returns

log_det – Logarithm of the matrix determinant.

Return type

float

exception harmonia.utils.PositiveDefinitenessWarning[source]

Emit a warning when a matrix is not positive definite.

harmonia.utils.is_positive_definite(matrix)[source]

Check the positive definiteness of a square matrix by attempting a Cholesky decomposition.

Parameters

matrix (float or complex, array_like) – Matrix.

Returns

Positive definiteness.

Return type

bool