From 84f0b697c03b971e6d4e292c31dfcaba470c2bf4 Mon Sep 17 00:00:00 2001 From: "Patrick J. Roddy" Date: Wed, 13 Nov 2024 17:08:19 +0000 Subject: [PATCH] gh-188: add docstrings to all functions and tidy docs (#381) Co-authored-by: Saransh Chopra Co-authored-by: Nicolas Tessore --- docs/conf.py | 13 +- glass/core/algorithm.py | 30 ++++- glass/core/array.py | 89 +++++++++++++- glass/ext/__init__.py | 15 +++ glass/fields.py | 209 ++++++++++++++++++++++++++++---- glass/galaxies.py | 115 ++++++++++-------- glass/lensing.py | 130 +++++++++++++++----- glass/observations.py | 97 +++++++++------ glass/points.py | 87 ++++++++++---- glass/shapes.py | 76 +++++++----- glass/shells.py | 259 +++++++++++++++++++++++++++++++++------- glass/user.py | 59 ++++++++- pyproject.toml | 1 + 13 files changed, 926 insertions(+), 254 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 874ec590..6cfa10c8 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -25,9 +25,11 @@ extensions = [ "matplotlib.sphinxext.plot_directive", "nbsphinx", + "sphinx.ext.autodoc", "sphinx.ext.intersphinx", "sphinx.ext.napoleon", "sphinx.ext.viewcode", + "sphinx_autodoc_typehints", "sphinxcontrib.katex", ] @@ -62,8 +64,8 @@ # This config value contains the locations and names of other projects that # should be linked to in this documentation. intersphinx_mapping = { + "numpy": ("https://numpy.org/doc/stable", None), "python": ("https://docs.python.org/3", None), - "numpy": ("https://numpy.org/doc/stable/", None), } @@ -74,9 +76,18 @@ # -- napoleon ---------------------------------------------------------------- +napoleon_custom_sections = [ + ("Returns", "params_style"), + ("Yields", "params_style"), +] napoleon_google_docstring = False +# -- sphinx-autodoc-typehints ------------------------------------------------ + +always_use_bars_union = True + + # -- plot_directive ---------------------------------------------------------- # Whether to show a link to the source in HTML (default: True). diff --git a/glass/core/algorithm.py b/glass/core/algorithm.py index 4cd6b953..09fd7c38 100644 --- a/glass/core/algorithm.py +++ b/glass/core/algorithm.py @@ -18,13 +18,37 @@ def nnls( Implementation of the algorithm due to [1] as described in [2]. + Parameters + ---------- + a + The matrix. + b + The vector. + tol + The tolerance for convergence. + maxiter + The maximum number of iterations. + + Returns + ------- + The non-negative least squares solution. + + Raises + ------ + ValueError + If ``a`` is not a matrix. + ValueError + If ``b`` is not a vector. + ValueError + If the shapes of ``a`` and ``b`` do not match. + References ---------- * [1] Lawson, C. L. and Hanson, R. J. (1995), Solving Least Squares - Problems. doi: 10.1137/1.9781611971217 + Problems. doi: 10.1137/1.9781611971217 * [2] Bro, R. and De Jong, S. (1997), A fast - non-negativity-constrained least squares algorithm. J. - Chemometrics, 11, 393-401. + non-negativity-constrained least squares algorithm. J. + Chemometrics, 11, 393-401. """ a = np.asanyarray(a) diff --git a/glass/core/array.py b/glass/core/array.py index c963c35b..0d2eef41 100644 --- a/glass/core/array.py +++ b/glass/core/array.py @@ -15,7 +15,19 @@ def broadcast_first( *arrays: npt.NDArray[np.float64], ) -> tuple[npt.NDArray[np.float64], ...]: - """Broadcast arrays, treating the first axis as common.""" + """ + Broadcast arrays, treating the first axis as common. + + Parameters + ---------- + arrays + The arrays to broadcast. + + Returns + ------- + The broadcasted arrays. + + """ arrays = tuple(np.moveaxis(a, 0, -1) if np.ndim(a) else a for a in arrays) arrays = np.broadcast_arrays(*arrays) return tuple(np.moveaxis(a, -1, 0) if np.ndim(a) else a for a in arrays) @@ -33,8 +45,15 @@ def broadcast_leading_axes( """ Broadcast all but the last N axes. - Returns the shape of the broadcast dimensions, and all input arrays - with leading axes matching that shape. + Parameters + ---------- + args + The arrays and the number of axes to keep. + + Returns + ------- + The shape of the broadcast dimensions, and all input arrays + with leading axes matching that shape. Examples -------- @@ -75,7 +94,31 @@ def ndinterp( # noqa: PLR0913 right: float | None = None, period: float | None = None, ) -> npt.NDArray[np.float64]: - """Interpolate multi-dimensional array over axis.""" + """ + Interpolate multi-dimensional array over axis. + + Parameters + ---------- + x + The x-coordinates. + xp + The x-coordinates of the data points. + fp + The function values corresponding to the x-coordinates in *xp*. + axis + The axis to interpolate over. + left + The value to return for x < xp[0]. + right + The value to return for x > xp[-1]. + period + The period of the function, used for interpolating periodic data. + + Returns + ------- + The interpolated array. + + """ return np.apply_along_axis( partial(np.interp, x, xp), axis, @@ -94,7 +137,23 @@ def trapz_product( ], axis: int = -1, ) -> npt.NDArray[np.float64]: - """Trapezoidal rule for a product of functions.""" + """ + Trapezoidal rule for a product of functions. + + Parameters + ---------- + f + The first function. + ff + The other functions. + axis + The axis along which to integrate. + + Returns + ------- + The integral of the product of the functions. + + """ x: npt.NDArray[np.float64] x, _ = f for x_, _ in ff: @@ -118,7 +177,25 @@ def cumtrapz( dtype: npt.DTypeLike | None = None, out: npt.NDArray[np.float64] | None = None, ) -> npt.NDArray[np.float64]: - """Cumulative trapezoidal rule along last axis.""" + """ + Cumulative trapezoidal rule along last axis. + + Parameters + ---------- + f + The function values. + x + The x-coordinates. + dtype + The output data type. + out + The output array. + + Returns + ------- + The cumulative integral of the function. + + """ if out is None: out = np.empty_like(f, dtype=dtype) diff --git a/glass/ext/__init__.py b/glass/ext/__init__.py index 2b55ab0f..bb7a8501 100644 --- a/glass/ext/__init__.py +++ b/glass/ext/__init__.py @@ -8,6 +8,21 @@ def _extend_path(path: list[str], name: str) -> list[str]: + """ + Extend the path to include the "ext" submodules of packages. + + Parameters + ---------- + path + The path to extend. + name + The name of the package. + + Returns + ------- + The extended path. + + """ import os.path from pkgutil import extend_path diff --git a/glass/fields.py b/glass/fields.py index 613a939e..0d02f4de 100644 --- a/glass/fields.py +++ b/glass/fields.py @@ -47,7 +47,30 @@ def iternorm( ) -> collections.abc.Generator[ tuple[int | None, npt.NDArray[np.float64], npt.NDArray[np.float64]] ]: - """Return the vector a and variance sigma^2 for iterative normal sampling.""" + """ + Return the vector a and variance sigma^2 for iterative normal sampling. + + Parameters + ---------- + k + The number of fields. + cov + The covariance matrix for the fields. + size + The size of the covariance matrix. + + Yields + ------ + The index, vector, and standard deviation for iterative sampling. + + Raises + ------ + TypeError + If the covariance matrix is not broadcastable to the given size. + ValueError + If the covariance matrix is not positive definite. + + """ n = (size,) if isinstance(size, int) else size m = np.zeros((*n, k, k)) @@ -104,7 +127,30 @@ def cls2cov( nf: int, nc: int, ) -> collections.abc.Generator[npt.NDArray[np.float64]]: - """Return array of cls as a covariance matrix for iterative sampling.""" + """ + Return array of Cls as a covariance matrix for iterative sampling. + + Parameters + ---------- + cls + Angular matter power spectra in *GLASS* ordering. + nl + The number of modes. + nf + The number of fields. + nc + The number of correlated fields. + + Yields + ------ + The covariance matrix for iterative sampling. + + Raises + ------ + ValueError + If negative values are found in the Cls. + + """ cov = np.zeros((nl, nc + 1)) end = 0 for j in range(nf): @@ -126,7 +172,23 @@ def multalm( *, inplace: bool = False, ) -> npt.NDArray[np.complex128]: - """Multiply alm by bl.""" + """ + Multiply alm by bl. + + Parameters + ---------- + alm + The alm to multiply. + bl + The bl to multiply. + inplace + Whether to perform the operation in place. + + Returns + ------- + The product of alm and bl. + + """ n = len(bl) out = np.asanyarray(alm) if inplace else np.copy(alm) for m in range(n): @@ -139,7 +201,23 @@ def transform_cls( tfm: str | typing.Callable[[npt.NDArray[np.float64]], npt.NDArray[np.float64]], pars: tuple[typing.Any, ...] = (), ) -> Cls: - """Transform Cls to Gaussian Cls.""" + """ + Transform Cls to Gaussian Cls. + + Parameters + ---------- + cls + Angular matter power spectra in *GLASS* ordering. + tfm + The transformation to apply. + pars + The parameters for the transformation. + + Returns + ------- + The transformed angular power spectra. + + """ gls = [] for cl in cls: if len(cl) > 0: @@ -171,6 +249,26 @@ def discretized_cls( applies the HEALPix pixel window function of the given ``nside``. If no arguments are given, no action is performed. + Parameters + ---------- + cls + Angular matter power spectra in *GLASS* ordering. + lmax + The maximum mode number to truncate the spectra. + ncorr + The number of correlated fields to keep. + nside + The resolution parameter for the HEALPix maps. + + Returns + ------- + The discretised angular power spectra. + + Raises + ------ + ValueError + If the length of the Cls array is not a triangle number. + """ if ncorr is not None: n = int((2 * len(cls)) ** 0.5) @@ -202,7 +300,21 @@ def lognormal_gls( cls: Cls, shift: float = 1.0, ) -> Cls: - """Compute Gaussian Cls for a lognormal random field.""" + """ + Compute Gaussian Cls for a lognormal random field. + + Parameters + ---------- + cls + Angular matter power spectra in *GLASS* ordering. + shift + The shift parameter for the lognormal transformation. + + Returns + ------- + The Gaussian angular power spectra for a lognormal random field. + + """ return transform_cls(cls, "lognormal", (shift,)) @@ -235,6 +347,26 @@ def generate_gaussian( Missing entries can be set to ``None``. + Parameters + ---------- + gls + The Gaussian angular power spectra for a random field. + nside + The resolution parameter for the HEALPix maps. + ncorr + The number of correlated fields. If not given, all fields are correlated. + rng + Random number generator. If not given, a default RNG is used. + + Yields + ------ + The Gaussian random fields. + + Raises + ------ + ValueError + If all gls are empty. + """ # get the default RNG if not given if rng is None: @@ -299,7 +431,27 @@ def generate_lognormal( ncorr: int | None = None, rng: np.random.Generator | None = None, ) -> collections.abc.Generator[npt.NDArray[np.float64]]: - """Sample lognormal random fields from Gaussian Cls iteratively.""" + """ + Sample lognormal random fields from Gaussian Cls iteratively. + + Parameters + ---------- + gls + The Gaussian angular power spectra for a lognormal random field. + nside + The resolution parameter for the HEALPix maps. + shift + The shift parameter for the lognormal transformation. + ncorr + The number of correlated fields. If not given, all fields are correlated. + rng + Random number generator. If not given, a default RNG is used. + + Yields + ------ + The lognormal random fields. + + """ for i, m in enumerate(generate_gaussian(gls, nside, ncorr=ncorr, rng=rng)): # compute the variance of the auto-correlation gl = gls[i * (i + 1) // 2] @@ -331,20 +483,22 @@ def getcl( """ Return a specific angular power spectrum from an array. - Return the angular power spectrum for indices *i* and *j* from an - array in *GLASS* ordering. - Parameters ---------- - cls: - List of angular power spectra in *GLASS* ordering. - i: + cls + Angular matter power spectra in *GLASS* ordering. + i Indices to return. - j: + j Indices to return. - lmax: + lmax Truncate the returned spectrum at this mode number. + Returns + ------- + The angular power spectrum for indices *i* and *j* from an + array in *GLASS* ordering. + """ if j > i: i, j = j, i @@ -366,29 +520,38 @@ def effective_cls( *, lmax: int | None = None, ) -> npt.NDArray[np.float64]: - r""" + """ Compute effective angular power spectra from weights. Computes a linear combination of the angular power spectra *cls* using the factors provided by *weights1* and *weights2*. Additional axes in *weights1* and *weights2* produce arrays of spectra. - Returns a dictionary of effective angular power spectra, where keys - correspond to the leading axes of *weights1* and *weights2*. - Parameters ---------- - cls: - Angular matter power spectra to combine, in *GLASS* ordering. - weights1: + cls + Angular matter power spectra in *GLASS* ordering. + weights1 Weight factors for spectra. The first axis must be equal to the number of fields. - weights2: + weights2 Second set of weights. If not given, *weights1* is used. - lmax: + lmax Truncate the angular power spectra at this mode number. If not given, the longest input in *cls* will be used. + Returns + ------- + A dictionary of effective angular power spectra, where keys + correspond to the leading axes of *weights1* and *weights2*. + + Raises + ------ + ValueError + If the length of *cls* is not a triangle number. + ValueError + If the shapes of *weights1* and *weights2* are incompatible. + """ from itertools import combinations_with_replacement, product diff --git a/glass/galaxies.py b/glass/galaxies.py index b5125e6d..4d14af49 100644 --- a/glass/galaxies.py +++ b/glass/galaxies.py @@ -46,18 +46,20 @@ def redshifts( This function samples *n* redshifts from a distribution that follows the given radial window function *w*. - Returns random redshifts following the radial window function. - Parameters ---------- - n: + n Number of redshifts to sample. If an array is given, the results are concatenated. - w: + w Radial window function. - rng: + rng Random number generator. If not given, a default RNG is used. + Returns + ------- + Random redshifts following the radial window function. + """ return redshifts_from_nz(n, w.za, w.wa, rng=rng, warn=False) @@ -80,26 +82,28 @@ def redshifts_from_nz( and redshifts are sampled independently for each extra dimension. The results are concatenated into a flat array. - Returns redshifts sampled from the given source distribution. For - inputs with extra dimensions, returns a flattened 1-D array of - samples from all populations. - Parameters ---------- - count: + count Number of redshifts to sample. If an array is given, its shape is broadcast against the leading axes of *z* and *nz*. - z: + z Source distribution. Leading axes are broadcast against the shape of *count*. - nz: + nz Source distribution. Leading axes are broadcast against the shape of *count*. - rng: + rng Random number generator. If not given, a default RNG is used. - warn: + warn Throw relevant warnings. + Returns + ------- + Redshifts sampled from the given source distribution. For + inputs with extra dimensions, returns a flattened 1-D array of + samples from all populations. + """ if warn: warnings.warn( @@ -157,27 +161,29 @@ def galaxy_shear( # noqa: PLR0913 Takes lensing maps for convergence and shear and produces a lensed ellipticity (shear) for each intrinsic galaxy ellipticity. - Returns an array of complex-valued observed galaxy shears - (lensed ellipticities). - Parameters ---------- - lon: + lon Array for galaxy longitudes. - lat: + lat Array for galaxy latitudes. - eps: + eps Array of galaxy :term:`ellipticity`. - kappa: + kappa HEALPix map for convergence. - gamma1: + gamma1 HEALPix maps for a component of shear. - gamma2: + gamma2 HEALPix maps for a component of shear. - reduced_shear: + reduced_shear If ``False``, galaxy shears are not reduced by the convergence. Default is ``True``. + Returns + ------- + An array of complex-valued observed galaxy shears + (lensed ellipticities). + """ nside = healpix.npix2nside(np.broadcast(kappa, gamma1, gamma2).shape[-1]) @@ -223,22 +229,29 @@ def gaussian_phz( Gaussian error with redshift-dependent standard deviation :math:`\sigma(z) = (1 + z) \sigma_0` [1]. - Returns photometric redshifts assuming Gaussian errors, of the same - shape as *z*. - Parameters ---------- - z: + z True redshifts. - sigma_0: + sigma_0 Redshift error in the tomographic binning at zero redshift. - lower: + lower Bounds for the returned photometric redshifts. - upper: + upper Bounds for the returned photometric redshifts. - rng: + rng Random number generator. If not given, a default RNG is used. + Returns + ------- + Photometric redshifts assuming Gaussian errors, of the same + shape as *z*. + + Raises + ------ + ValueError + If the bounds are not consistent. + Warnings -------- The *lower* and *upper* bounds are implemented using plain rejection @@ -253,7 +266,7 @@ def gaussian_phz( References ---------- * [1] Amara A., Réfrégier A., 2007, MNRAS, 381, 1018. - doi:10.1111/j.1365-2966.2007.12271.x + doi:10.1111/j.1365-2966.2007.12271.x Examples -------- @@ -307,34 +320,36 @@ def _kappa_ia_nla( # noqa: PLR0913 r""" Effective convergence from intrinsic alignments using the NLA model. - Returns the effective convergence due to intrinsic alignments. - Parameters ---------- - delta: + delta Matter density contrast. - zeff: + zeff Effective redshift of the matter field. - a_ia: + a_ia Intrinsic alignments amplitude. - cosmo: + cosmo Cosmology instance. - z0: + z0 Reference redshift for the redshift dependence. - eta: + eta Power of the redshift dependence. - lbar: + lbar Mean luminosity of the galaxy sample. - l0: + l0 Reference luminosity for the luminosity dependence. - beta: + beta Power of the luminosity dependence. + Returns + ------- + The effective convergence due to intrinsic alignments. + Notes ----- The Non-linear Alignments Model (NLA) describes an effective convergence :math:`\kappa_{\rm IA}` that models the effect of - intrinsic alignments. It is computed from the matter density + intrinsic alignments. It is computed from the matter density contrast :math:`\delta` as [1] [3] .. math:: @@ -366,15 +381,15 @@ def _kappa_ia_nla( # noqa: PLR0913 References ---------- * [1] Catelan P., Kamionkowski M., Blandford R. D., 2001, MNRAS, - 320, L7. doi:10.1046/j.1365-8711.2001.04105.x + 320, L7. doi:10.1046/j.1365-8711.2001.04105.x * [2] Hirata C. M., Seljak U., 2004, PhRvD, 70, 063526. - doi:10.1103/PhysRevD.70.063526 + doi:10.1103/PhysRevD.70.063526 * [3] Bridle S., King L., 2007, NJPh, 9, 444. - doi:10.1088/1367-2630/9/12/444 + doi:10.1088/1367-2630/9/12/444 * [4] Johnston, H., Georgiou, C., Joachimi, B., et al., 2019, - A&A, 624, A30. doi:10.1051/0004-6361/201834714 + A&A, 624, A30. doi:10.1051/0004-6361/201834714 * [5] Tessore, N., Loureiro, A., Joachimi, B., et al., 2023, - OJAp, 6, 11. doi:10.21105/astro.2302.01942 + OJAp, 6, 11. doi:10.21105/astro.2302.01942 """ c1 = 5e-14 / cosmo.h**2 # Solar masses per cubic Mpc diff --git a/glass/lensing.py b/glass/lensing.py index abcd3db3..a9b2e6a9 100644 --- a/glass/lensing.py +++ b/glass/lensing.py @@ -61,27 +61,31 @@ def from_convergence( # noqa: PLR0913 deflection potential, deflection, and shear maps. The maps are computed via spherical harmonic transforms. - Returns the maps of: - - * deflection potential if ``potential`` is true. - * potential (complex) if ``deflection`` is true. - * shear (complex) if ``shear`` is true. - Parameters ---------- - kappa: + kappa HEALPix map of the convergence field. - lmax: + lmax Maximum angular mode number to use in the transform. - potential: + potential Which lensing maps to return. - deflection: + deflection Which lensing maps to return. - shear: + shear Which lensing maps to return. - discretized: + discretized Correct the pixel window function in output maps. + Returns + ------- + psi + Map of the lensing (or deflection) potential. Only returned if + ``potential`` is true. + alpha + Map of the deflection (complex). Only returned if ``deflection`` if true. + gamma + Map of the shear (complex). Only returned if ``shear`` is true. + Notes ----- The weak lensing fields are computed from the convergence or @@ -152,7 +156,7 @@ def from_convergence( # noqa: PLR0913 References ---------- * [1] Tessore N., et al., OJAp, 6, 11 (2023). - doi:10.21105/astro.2302.01942 + doi:10.21105/astro.2302.01942 """ # no output means no computation, return empty tuple @@ -228,14 +232,27 @@ def shear_from_convergence( *, discretized: bool = True, ) -> npt.NDArray[np.float64]: - r""" + """ Weak lensing shear from convergence. + Computes the shear from the convergence using a spherical harmonic + transform. + .. deprecated:: 2023.6 Use the more general :func:`from_convergence` function instead. - Computes the shear from the convergence using a spherical harmonic - transform. + Parameters + ---------- + kappa + The convergence map. + lmax + The maximum angular mode number to use in the transform. + discretized + Whether to correct the pixel window function in the output map. + + Returns + ------- + The shear map. """ nside = hp.get_nside(kappa) @@ -270,7 +287,15 @@ class MultiPlaneConvergence: """Compute convergence fields iteratively from multiple matter planes.""" def __init__(self, cosmo: Cosmology) -> None: - """Create a new instance to iteratively compute the convergence.""" + """ + Create a new instance to iteratively compute the convergence. + + Parameters + ---------- + cosmo + Cosmology instance. + + """ self.cosmo = cosmo # set up initial values of variables @@ -290,6 +315,13 @@ def add_window(self, delta: npt.NDArray[np.float64], w: RadialWindow) -> None: The lensing weight is computed from the window function, and the source plane redshift is the effective redshift of the window. + Parameters + ---------- + delta + The mass plane. + w + The window function. + """ zsrc = w.zeff lens_weight = float( @@ -312,7 +344,24 @@ def add_plane( zsrc: float, wlens: float = 1.0, ) -> None: - """Add a mass plane at redshift ``zsrc`` to the convergence.""" + """ + Add a mass plane at redshift ``zsrc`` to the convergence. + + Parameters + ---------- + delta + The mass plane. + zsrc + The redshift of the source plane. + wlens + The weight of the mass plane. + + Raises + ------ + ValueError + If the source redshift is not increasing. + + """ if zsrc <= self.z3: msg = "source redshift must be increasing" raise ValueError(msg) @@ -379,7 +428,21 @@ def multi_plane_matrix( shells: collections.abc.Sequence[RadialWindow], cosmo: Cosmology, ) -> npt.NDArray[np.float64]: - """Compute the matrix of lensing contributions from each shell.""" + """ + Compute the matrix of lensing contributions from each shell. + + Parameters + ---------- + shells + The shells of the mass distribution. + cosmo + Cosmology instance. + + Returns + ------- + The matrix of lensing contributions. + + """ mpc = MultiPlaneConvergence(cosmo) wmat = np.eye(len(shells)) for i, w in enumerate(shells): @@ -403,18 +466,25 @@ def multi_plane_weights( redshift distribution :math:`n(z)` into the lensing efficiency sometimes denoted :math:`g(z)` or :math:`q(z)`. - Returns the relative lensing weight of each shell. - Parameters ---------- - weights: + weights Relative weight of each shell. The first axis must broadcast against the number of shells, and is normalised internally. - shells: + shells Window functions of the shells. - cosmo: + cosmo Cosmology instance. + Returns + ------- + The relative lensing weight of each shell. + + Raises + ------ + ValueError + If the shape of *weights* does not match the number of shells. + """ # ensure shape of weights ends with the number of shells shape = np.shape(weights) @@ -442,18 +512,20 @@ def deflect( Takes an array of :term:`deflection` values and applies them to the given positions. - Returns the longitudes and latitudes after deflection. - Parameters ---------- - lon: + lon Longitudes to be deflected. - lat: + lat Latitudes to be deflected. - alpha: + alpha Deflection values. Must be complex-valued or have a leading axis of size 2 for the real and imaginary component. + Returns + ------- + The longitudes and latitudes after deflection. + Notes ----- Deflections on the sphere are :term:`defined ` as diff --git a/glass/observations.py b/glass/observations.py index d4fdac07..6581ca98 100644 --- a/glass/observations.py +++ b/glass/observations.py @@ -48,21 +48,25 @@ def vmap_galactic_ecliptic( the galactic and ecliptic planes. The location of the stripes is set with optional parameters. - Returns a HEALPix :term:`visibility map`. - Parameters ---------- - nside: + nside The NSIDE parameter of the resulting HEALPix map. - galactic: + galactic The location of the galactic plane in the respective coordinate system. - ecliptic: + ecliptic The location of the ecliptic plane in the respective coordinate system. + Returns + ------- + A HEALPix :term:`visibility map`. + Raises ------ TypeError - If the ``galactic`` or ``ecliptic`` arguments are not pairs of numbers. + If the ``galactic`` argument is not a pair of numbers. + TypeError + If the ``ecliptic`` argument is not a pair of numbers. """ if np.ndim(galactic) != 1 or len(galactic) != 2: @@ -86,7 +90,7 @@ def gaussian_nz( *, norm: npt.NDArray[np.float64] | None = None, ) -> npt.NDArray[np.float64]: - r""" + """ Gaussian redshift distribution. The redshift follows a Gaussian distribution with the given mean and @@ -95,19 +99,21 @@ def gaussian_nz( If ``mean`` or ``sigma`` are array_like, their axes will be the leading axes of the redshift distribution. - Returns the redshift distribution at the given ``z`` values. - Parameters ---------- - z: + z Redshift values of the distribution. - mean: + mean Mean(s) of the redshift distribution. - sigma: + sigma Standard deviation(s) of the redshift distribution. - norm: + norm If given, the normalisation of the distribution. + Returns + ------- + The redshift distribution at the given ``z`` values. + """ mean = np.reshape(mean, np.shape(mean) + (1,) * np.ndim(z)) sigma = np.reshape(sigma, np.shape(sigma) + (1,) * np.ndim(z)) @@ -138,21 +144,23 @@ def smail_nz( The redshift follows the Smail et al. [1] redshift distribution. - Returns the redshift distribution at the given ``z`` values. - Parameters ---------- - z: + z Redshift values of the distribution. - z_mode: + z_mode Mode of the redshift distribution, must be positive. - alpha: + alpha Power law exponent (z/z0)^\alpha, must be positive. - beta: + beta Log-power law exponent exp[-(z/z0)^\beta], must be positive. - norm: + norm If given, the normalisation of the distribution. + Returns + ------- + The redshift distribution at the given ``z`` values. + Notes ----- The probability distribution function :math:`p(z)` for redshift :math:`z` @@ -201,19 +209,26 @@ def fixed_zbins( This function creates contiguous tomographic redshift bins of fixed size. It takes either the number or size of the bins. - Returns a list of redshift bin edges. - Parameters ---------- - zmin: + zmin Extent of the redshift binning. - zmax: + zmax Extent of the redshift binning. - nbins: + nbins Number of redshift bins. Only one of ``nbins`` and ``dz`` can be given. - dz: + dz Size of redshift bin. Only one of ``nbins`` and ``dz`` can be given. + Returns + ------- + A list of redshift bin edges. + + Raises + ------ + ValueError + If both ``nbins`` and ``dz`` are given. + """ if nbins is not None and dz is None: zbinedges = np.linspace(zmin, zmax, nbins + 1) @@ -237,17 +252,19 @@ def equal_dens_zbins( This function subdivides a source redshift distribution into ``nbins`` tomographic redshift bins with equal density. - Returns a list of redshift bin edges. - Parameters ---------- - z: + z The source redshift distribution. Must be one-dimensional. - nz: + nz The source redshift distribution. Must be one-dimensional. - nbins: + nbins Number of redshift bins. + Returns + ------- + A list of redshift bin edges. + """ # compute the normalised cumulative distribution function # first compute the cumulative integral (by trapezoidal rule) @@ -276,20 +293,22 @@ def tomo_nz_gausserr( standard deviation of the Gaussian depends on redshift and is given by ``sigma(z) = sigma_0*(1 + z)``. - Returns the tomographic redshift bins convolved with a gaussian error. - Array has a shape (nbins, len(z)) - Parameters ---------- - z: + z The true source redshift distribution. Must be one-dimensional. - nz: + nz The true source redshift distribution. Must be one-dimensional. - sigma_0: + sigma_0 Redshift error in the tomographic binning at zero redshift. - zbins: + zbins List of redshift bin edges. + Returns + ------- + The tomographic redshift bins convolved with a gaussian error. + Array has a shape (nbins, len(z)) + See Also -------- equal_dens_zbins: @@ -300,7 +319,7 @@ def tomo_nz_gausserr( References ---------- * [1] Amara A., Réfrégier A., 2007, MNRAS, 381, 1018. - doi:10.1111/j.1365-2966.2007.12271.x + doi:10.1111/j.1365-2966.2007.12271.x """ # converting zbins into an array: diff --git a/glass/points.py b/glass/points.py index 7c0a5c2b..831c083a 100644 --- a/glass/points.py +++ b/glass/points.py @@ -60,17 +60,19 @@ def effective_bias( and computes an effective bias parameter :math:`\bar{b}` for a given window function :math:`w(z)`. - Returns the effective bias parameter for the window. - Parameters ---------- - z: + z Redshifts and values of the bias function :math:`b(z)`. - bz: + bz Redshifts and values of the bias function :math:`b(z)`. - w: + w The radial window function :math:`w(z)`. + Returns + ------- + The effective bias parameter for the window. + Notes ----- The effective bias parameter :math:`\bar{b}` is computed using the @@ -93,7 +95,21 @@ def linear_bias( delta: npt.NDArray[np.float64], b: float | npt.NDArray[np.float64], ) -> npt.NDArray[np.float64]: - r"""Linear bias model :math:`\delta_g = b \, \delta`.""" + r""" + Linear bias model :math:`\delta_g = b \, \delta`. + + Parameters + ---------- + delta + The input density contrast. + b + The bias parameter. + + Returns + ------- + The density contrast after biasing. + + """ return b * delta @@ -101,7 +117,21 @@ def loglinear_bias( delta: npt.NDArray[np.float64], b: float | npt.NDArray[np.float64], ) -> npt.NDArray[np.float64]: - r"""log-linear bias model :math:`\ln(1 + \delta_g) = b \ln(1 + \delta)`.""" + r""" + Log-linear bias model :math:`\ln(1 + \delta_g) = b \ln(1 + \delta)`. + + Parameters + ---------- + delta + The input density contrast. + b + The bias parameter. + + Returns + ------- + The density contrast after biasing. + + """ delta_g = np.log1p(delta) delta_g *= b np.expm1(delta_g, out=delta_g) @@ -148,38 +178,43 @@ def positions_from_delta( # noqa: PLR0912, PLR0913, PLR0915 Parameters ---------- - ngal: + ngal Number density, expected number of points per arcmin2. - delta: + delta Map of the input density contrast. This is fed into the bias model to produce the density contrast for sampling. - bias: + bias Bias parameter, is passed as an argument to the bias model. - vis: + vis Visibility map for the observed points. This is multiplied with the full sky number count map, and must hence be of compatible shape. - bias_model: + bias_model The bias model to apply. If a string, refers to a function in the :mod:`~glass.points` module, e.g. ``'linear'`` for :func:`linear_bias()` or ``'loglinear'`` for :func:`loglinear_bias`. - remove_monopole: + remove_monopole If true, the monopole of the density contrast after biasing is fixed to zero. - batch: + batch Maximum number of positions to yield in one batch. - rng: + rng Random number generator. If not given, a default RNG is used. Yields ------ - lon: + lon Columns of longitudes for the sampled points. - lat: + lat Columns of latitudes for the sampled points. - count: + count The number of sampled points If multiple populations are sampled, an array of counts in the shape of the extra dimensions is returned. + Raises + ------ + TypeError + If the bias model is not a string or callable. + """ # get default RNG if not given if rng is None: @@ -300,10 +335,10 @@ def uniform_positions( Parameters ---------- - ngal: + ngal Number density, expected number of positions per arcmin2. - rng: - Random number generator. If not given, a default RNG will be used. + rng + Random number generator. If not given, a default RNG is used. Yields ------ @@ -361,16 +396,18 @@ def position_weights( redshift distribution and bias factor :math:`n(z) \, b(z)` for the discretised shells. - Returns the relative weight of each shell for angular clustering. - Parameters ---------- - densities: + densities Density of points in each shell. The first axis must broadcast against the number of shells, and is normalised internally. - bias: + bias Value or values of the linear bias parameter for each shell. + Returns + ------- + The relative weight of each shell for angular clustering. + """ # bring densities and bias into the same shape if bias is not None: diff --git a/glass/shapes.py b/glass/shapes.py index 33cd36dc..8bc0a7f0 100644 --- a/glass/shapes.py +++ b/glass/shapes.py @@ -35,25 +35,26 @@ def triaxial_axis_ratio( *, rng: np.random.Generator | None = None, ) -> npt.NDArray[np.float64]: - r""" + """ Axis ratio of a randomly projected triaxial ellipsoid. Given the two axis ratios `1 >= zeta >= xi` of a randomly oriented triaxial ellipsoid, computes the axis ratio `q` of the projection. - Returns the axis ratio of the randomly projected ellipsoid. - Parameters ---------- - zeta: + zeta Axis ratio of intermediate and major axis. - xi: + xi Axis ratio of minor and major axis. - size: - Size of the random draw. If `None` is given, size is inferred from - other inputs. - rng: - Random number generator. If not given, a default RNG will be used. + size + Size of the random draw. + rng + Random number generator. If not given, a default RNG is used. + + Returns + ------- + The axis ratio of the randomly projected ellipsoid. Notes ----- @@ -115,22 +116,24 @@ def ellipticity_ryden04( # noqa: PLR0913 standard deviation :math:`\sigma_\gamma` [2]. Both distributions are truncated to produce ratios in the range 0 to 1 using rejection sampling. - Returns an array of :term:`ellipticity` from projected axis ratios. - Parameters ---------- - mu: + mu Mean of the truncated normal for :math:`\log(1 - B/A)`. - sigma: + sigma Standard deviation for :math:`\log(1 - B/A)`. - gamma: + gamma Mean of the truncated normal for :math:`1 - C/B`. - sigma_gamma: + sigma_gamma Standard deviation for :math:`1 - C/B`. - size: - Sample size. If ``None``, the size is inferred from the parameters. - rng: - Random number generator. If not given, a default RNG will be used. + size + Sample size. + rng + Random number generator. If not given, a default RNG is used. + + Returns + ------- + An array of :term:`ellipticity` from projected axis ratios. References ---------- @@ -183,7 +186,7 @@ def ellipticity_gaussian( *, rng: np.random.Generator | None = None, ) -> npt.NDArray[np.complex128]: - r""" + """ Sample Gaussian galaxy ellipticities. The ellipticities are sampled from a normal distribution with @@ -192,17 +195,19 @@ def ellipticity_gaussian( this function with too large values of ``sigma``, or the sampling will become inefficient. - Returns an array of galaxy :term:`ellipticity`. - Parameters ---------- - count: + count Number of ellipticities to be sampled. - sigma: + sigma Standard deviation in each component. - rng: + rng Random number generator. If not given, a default RNG is used. + Returns + ------- + An array of galaxy :term:`ellipticity`. + """ # default RNG if not provided if rng is None: @@ -237,23 +242,30 @@ def ellipticity_intnorm( *, rng: np.random.Generator | None = None, ) -> npt.NDArray[np.complex128]: - r""" + """ Sample galaxy ellipticities with intrinsic normal distribution. The ellipticities are sampled from an intrinsic normal distribution with standard deviation ``sigma`` for each component. - Returns an array of galaxy :term:`ellipticity`. - Parameters ---------- - count: + count Number of ellipticities to sample. - sigma: + sigma Standard deviation of the ellipticity in each component. - rng: + rng Random number generator. If not given, a default RNG is used. + Returns + ------- + An array of galaxy :term:`ellipticity`. + + Raises + ------ + ValueError + If the standard deviation is not in the range [0, sqrt(0.5)]. + """ # default RNG if not provided if rng is None: diff --git a/glass/shells.py b/glass/shells.py index 9414e746..34f5a2ca 100644 --- a/glass/shells.py +++ b/glass/shells.py @@ -65,7 +65,21 @@ def distance_weight( z: npt.NDArray[np.float64], cosmo: Cosmology, ) -> npt.NDArray[np.float64]: - """Uniform weight in comoving distance.""" + """ + Uniform weight in comoving distance. + + Parameters + ---------- + z + The redshifts at which to evaluate the weight. + cosmo + Cosmology instance. + + Returns + ------- + The weight function evaluated at redshifts *z*. + + """ return 1 / cosmo.ef(z) # type: ignore[no-any-return] @@ -73,7 +87,21 @@ def volume_weight( z: npt.NDArray[np.float64], cosmo: Cosmology, ) -> npt.NDArray[np.float64]: - """Uniform weight in comoving volume.""" + """ + Uniform weight in comoving volume. + + Parameters + ---------- + z + The redshifts at which to evaluate the weight. + cosmo + Cosmology instance. + + Returns + ------- + The weight function evaluated at redshifts *z*. + + """ return cosmo.xm(z) ** 2 / cosmo.ef(z) # type: ignore[no-any-return] @@ -81,7 +109,21 @@ def density_weight( z: npt.NDArray[np.float64], cosmo: Cosmology, ) -> npt.NDArray[np.float64]: - """Uniform weight in matter density.""" + """ + Uniform weight in matter density. + + Parameters + ---------- + z + The redshifts at which to evaluate the weight. + cosmo + Cosmology instance. + + Returns + ------- + The weight function evaluated at redshifts *z*. + + """ return cosmo.rho_m_z(z) * cosmo.xm(z) ** 2 / cosmo.ef(z) # type: ignore[no-any-return] @@ -116,16 +158,17 @@ class RadialWindow(typing.NamedTuple): Attributes ---------- - za: + za Redshift array; the abscissae of the window function. - wa: + wa Weight array; the values (ordinates) of the window function. - zeff: + zeff Effective redshift of the window. Methods ------- _replace + Create a new instance with changed attribute values. """ @@ -153,17 +196,24 @@ def tophat_windows( Their effective redshifts are the mean redshifts of the (weighted) tophat bins. - Returns a list of window functions. - Parameters ---------- - zbins: + zbins Redshift bin edges for the tophat window functions. - dz: + dz Approximate spacing of the redshift grid. - weight: + weight If given, a weight function to be applied to the window functions. + Returns + ------- + A list of window functions. + + Raises + ------ + ValueError + If the number of redshift bins is less than 2. + See Also -------- :ref:`user-window-functions` @@ -215,17 +265,24 @@ def linear_windows( The resulting windows functions are :class:`RadialWindow` instances. Their effective redshifts correspond to the given nodes. - Returns a list of window functions. - Parameters ---------- - zgrid: + zgrid Redshift grid for the triangular window functions. - dz: + dz Approximate spacing of the redshift grid. - weight: + weight If given, a weight function to be applied to the window functions. + Returns + ------- + A list of window functions. + + Raises + ------ + ValueError + If the number of nodes is less than 3. + See Also -------- :ref:`user-window-functions` @@ -273,17 +330,24 @@ def cubic_windows( The resulting windows functions are :class:`RadialWindow` instances. Their effective redshifts correspond to the given nodes. - Returns a list of window functions. - Parameters ---------- - zgrid: + zgrid Redshift grid for the cubic spline window functions. - dz: + dz Approximate spacing of the redshift grid. - weight: + weight If given, a weight function to be applied to the window functions. + Returns + ------- + A list of window functions. + + Raises + ------ + ValueError + If the number of nodes is less than 3. + See Also -------- :ref:`user-window-functions` @@ -333,17 +397,19 @@ def restrict( the function and window over the support of the window. Intermediate function values are found by linear interpolation - Returns the restricted function - Parameters ---------- - z: + z The function to be restricted. - f: + f The function to be restricted. - w: + w The window function for the restriction. + Returns + ------- + The restricted function + """ z_ = np.compress(np.greater(z, w.za[0]) & np.less(z, w.za[-1]), z) zr = np.union1d(w.za, z_) @@ -373,23 +439,30 @@ def partition( The window functions are given by the sequence *shells* of :class:`RadialWindow` or compatible entries. - Returns the weights of the partition, where the leading axis corresponds to - *shells*. - Parameters ---------- - z: + z The function to be partitioned. If *f* is multi-dimensional, its last axis must agree with *z*. - fz: + fz The function to be partitioned. If *f* is multi-dimensional, its last axis must agree with *z*. - shells: + shells Ordered sequence of window functions for the partition. - method: + method Method for the partition. See notes for description. The options are "lstsq", "nnls", "restrict". + Returns + ------- + The weights of the partition, where the leading axis corresponds to + *shells*. + + Raises + ------ + ValueError + If the method is not recognised. + Notes ----- Formally, if :math:`w_i` are the normalised window functions, @@ -398,6 +471,7 @@ def partition( approximate solution of .. math:: + \begin{pmatrix} w_1(z_1) \Delta z_1 & w_2(z_1) \, \Delta z_1 & \cdots \\ w_1(z_2) \Delta z_2 & w_2(z_2) \, \Delta z_2 & \cdots \\ @@ -418,6 +492,7 @@ def partition( equals the integral of the target function, .. math:: + \begin{pmatrix} w_1(z_1) \Delta z_1 & w_2(z_1) \, \Delta z_1 & \cdots \\ w_1(z_2) \Delta z_2 & w_2(z_2) \, \Delta z_2 & \cdots \\ @@ -469,7 +544,25 @@ def partition_lstsq( *, sumtol: float = 0.01, ) -> npt.NDArray[np.float64]: - """Least-squares partition.""" + """ + Least-squares partition. + + Parameters + ---------- + z + The function to be partitioned. + fz + The function to be partitioned. + shells + Ordered sequence of window functions. + sumtol + Tolerance for the sum of the partition. + + Returns + ------- + The partition. + + """ # make sure nothing breaks sumtol = max(sumtol, 1e-4) @@ -535,8 +628,20 @@ def partition_nnls( """ Non-negative least-squares partition. - Uses the ``nnls()`` algorithm from ``scipy.optimize`` and thus - requires SciPy. + Parameters + ---------- + z + The function to be partitioned. + fz + The function to be partitioned. + shells + Ordered sequence of window functions. + sumtol + Tolerance for the sum of the partition. + + Returns + ------- + The partition. """ from glass.core.algorithm import nnls @@ -619,7 +724,23 @@ def partition_restrict( fz: npt.NDArray[np.float64], shells: collections.abc.Sequence[RadialWindow], ) -> npt.NDArray[np.float64]: - """Partition by restriction and integration.""" + """ + Partition by restriction and integration. + + Parameters + ---------- + z + The function to be partitioned. + fz + The function to be partitioned. + shells + Ordered sequence of window functions. + + Returns + ------- + The partition. + + """ part = np.empty((len(shells),) + np.shape(fz)[:-1]) for i, w in enumerate(shells): zr, fr = restrict(z, fz, w) @@ -638,7 +759,30 @@ def redshift_grid( dz: float | None = None, num: int | None = None, ) -> npt.NDArray[np.float64]: - """Redshift grid with uniform spacing in redshift.""" + """ + Redshift grid with uniform spacing in redshift. + + Parameters + ---------- + zmin + The minimum redshift. + zmax + The maximum redshift. + dz + The redshift spacing. + num + The number redshift samples. + + Returns + ------- + The redshift grid. + + Raises + ------ + ValueError + If both ``dz`` and ``num`` are given. + + """ if dz is not None and num is None: z = np.arange(zmin, np.nextafter(zmax + dz, zmax), dz) elif dz is None and num is not None: @@ -657,7 +801,32 @@ def distance_grid( dx: float | None = None, num: int | None = None, ) -> npt.NDArray[np.float64]: - """Redshift grid with uniform spacing in comoving distance.""" + """ + Redshift grid with uniform spacing in comoving distance. + + Parameters + ---------- + cosmo + Cosmology instance. + zmin + The minimum redshift. + zmax + The maximum redshift. + dx + The comoving distance spacing. + num + The number of samples. + + Returns + ------- + The redshift grid. + + Raises + ------ + ValueError + If both ``dx`` and ``num`` are given. + + """ xmin, xmax = cosmo.dc(zmin), cosmo.dc(zmax) if dx is not None and num is None: x = np.arange(xmin, np.nextafter(xmax + dx, xmax), dx) @@ -685,18 +854,20 @@ def combine( The window functions are given by the sequence *shells* of :class:`RadialWindow` or compatible entries. - Returns a linear combination of window functions, evaluated in *z*. - Parameters ---------- - z: + z Redshifts *z* in which to evaluate the combined function. - weights: + weights Weights of the linear combination, where the leading axis corresponds to *shells*. - shells: + shells Ordered sequence of window functions to be combined. + Returns + ------- + A linear combination of window functions, evaluated in *z*. + See Also -------- partition: diff --git a/glass/user.py b/glass/user.py index 11edd66c..58f1860e 100644 --- a/glass/user.py +++ b/glass/user.py @@ -46,6 +46,13 @@ def save_cls( Uses :func:`numpy.savez` internally. The filename should therefore have a ``.npz`` suffix, or it will be given one. + Parameters + ---------- + filename + The name of the file to save to. + cls + Angular matter power spectra in *GLASS* ordering. + """ split = np.cumsum([len(cl) for cl in cls[:-1]]) values = np.concatenate(cls) @@ -60,6 +67,15 @@ def load_cls( Uses :func:`numpy.load` internally. + Parameters + ---------- + filename + The name of the file to load from. + + Returns + ------- + The list of Cls. + """ with np.load(filename) as npz: values = npz["values"] @@ -75,7 +91,17 @@ class _FitsWriter: """ def __init__(self, fits: fitsio.FITS, ext: str | None = None) -> None: - """Create a new, uninitialised writer.""" + """ + Create a new, uninitialised writer. + + Parameters + ---------- + fits + The fits object. + ext + The file extension. + + """ self.fits = fits self.ext = ext @@ -84,7 +110,17 @@ def _append( data: npt.NDArray[np.float64] | list[npt.NDArray[np.float64]], names: list[str] | None = None, ) -> None: - """Write the FITS file.""" + """ + Write the FITS file. + + Parameters + ---------- + data + The data to write. + names + The names of the columns. + + """ if self.ext is None or self.ext not in self.fits: self.fits.write_table(data, names=names, extname=self.ext) if self.ext is None: @@ -105,6 +141,14 @@ def write( Pass either a positional variable (data) or multiple named arguments (**columns) + + Parameters + ---------- + data + The data to write. + columns + The columns to write. + """ # if data is given, write it as it is if data is not None: @@ -137,6 +181,17 @@ def write_catalog( .. note:: Requires the ``fitsio`` package. + Parameters + ---------- + filename + The name of the file to write to. + ext + The file extension. + + Yields + ------ + The writer object. + """ import fitsio diff --git a/pyproject.toml b/pyproject.toml index b3ea3561..6e510f68 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,6 +46,7 @@ docs = [ "matplotlib", "nbsphinx", "sphinx", + "sphinx_autodoc_typehints", "sphinxcontrib-katex", ] examples = [