diff --git a/glass/galaxies.py b/glass/galaxies.py index 327fc1de..df873a33 100644 --- a/glass/galaxies.py +++ b/glass/galaxies.py @@ -29,8 +29,6 @@ from glass.core.array import broadcast_leading_axes, cumulative_trapezoid if typing.TYPE_CHECKING: - from cosmology import Cosmology - from glass.shells import RadialWindow @@ -303,105 +301,3 @@ def gaussian_phz( trunc = trunc[(znew < lower) | (znew > upper)] return zphot - - -def _kappa_ia_nla( # noqa: PLR0913 - delta: npt.NDArray[np.float64], - zeff: float, - a_ia: float, - cosmo: Cosmology, - *, - z0: float = 0.0, - eta: float = 0.0, - lbar: float = 0.0, - l0: float = 1e-9, - beta: float = 0.0, -) -> npt.NDArray[np.float64]: - r""" - Effective convergence from intrinsic alignments using the NLA model. - - Parameters - ---------- - delta - Matter density contrast. - zeff - Effective redshift of the matter field. - a_ia - Intrinsic alignments amplitude. - cosmo - Cosmology instance. - z0 - Reference redshift for the redshift dependence. - eta - Power of the redshift dependence. - lbar - Mean luminosity of the galaxy sample. - l0 - Reference luminosity for the luminosity dependence. - 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 - contrast :math:`\delta` as [1] [3] - - .. math:: - - \kappa_{\rm IA} = f_{\rm NLA} \, \delta \;, - - where the NLA factor :math:`f_{\rm NLA}` is defined as [4] [5] - - .. math:: - - f_{\rm{NLA}} - = -A_{\rm IA} \, \frac{C_1 \, \bar{\rho}(z)}{D(z)} \, - \biggl(\frac{1+z}{1+z_0}\biggr)^\eta \, - \biggl(\frac{\bar{L}}{L_0}\biggr)^\beta \;, - - with - - * :math:`A_{\rm IA}` the intrinsic alignments amplitude, - * :math:`C_1` a normalisation constant [2], - * :math:`z` the effective redshift of the model, - * :math:`\bar{\rho}` the mean matter density, - * :math:`D` the growth factor, - * :math:`\eta` the power that describes the redshift-dependence with - respect to :math:`z_0`, - * :math:`\bar{L}` the mean luminosity of the galaxy sample, and - * :math:`\beta` the power that describes the luminosity-dependence - :math:`\bar{L}` with respect to :math:`L_0`. - - References - ---------- - * [1] Catelan P., Kamionkowski M., Blandford R. D., 2001, MNRAS, - 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 - * [3] Bridle S., King L., 2007, NJPh, 9, 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 - * [5] Tessore, N., Loureiro, A., Joachimi, B., et al., 2023, - OJAp, 6, 11. doi:10.21105/astro.2302.01942 - - """ - c1 = 5e-14 / cosmo.h**2 # Solar masses per cubic Mpc - rho_c1 = c1 * cosmo.rho_c0 - - prefactor = -a_ia * rho_c1 * cosmo.Om - inverse_linear_growth = 1.0 / cosmo.gf(zeff) - redshift_dependence = ((1 + zeff) / (1 + z0)) ** eta - luminosity_dependence = (lbar / l0) ** beta - - f_nla = ( - prefactor * inverse_linear_growth * redshift_dependence * luminosity_dependence - ) - - return delta * f_nla # type: ignore[no-any-return] diff --git a/glass/lensing.py b/glass/lensing.py index 981c025f..b46ff389 100644 --- a/glass/lensing.py +++ b/glass/lensing.py @@ -37,11 +37,11 @@ import numpy as np import numpy.typing as npt +from cosmology.api import CosmologyConstantsNamespace, StandardCosmology + if typing.TYPE_CHECKING: import collections.abc - from cosmology import Cosmology - from glass.shells import RadialWindow @@ -397,7 +397,10 @@ def shear_from_convergence( class MultiPlaneConvergence: """Compute convergence fields iteratively from multiple matter planes.""" - def __init__(self, cosmo: Cosmology) -> None: + def __init__( + self, + cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]], + ) -> None: """ Create a new instance to iteratively compute the convergence. @@ -477,15 +480,19 @@ def add_plane( w2, self.w3 = self.w3, wlens # extrapolation law - x2, self.x3 = self.x3, self.cosmo.xm(self.z3) + hubble_length = CosmologyConstantsNamespace.c / StandardCosmology.H0 + x2 = self.x3 + self.x3 = self.cosmo.transverse_comoving_distance(self.z3) / hubble_length r12 = self.r23 - r13, self.r23 = self.cosmo.xm([z1, self.z2], self.z3) / self.x3 + r13, self.r23 = self.cosmo.transverse_comoving_distance( + [z1, self.z2], self.z3 + ) / (hubble_length * self.x3) t = r13 / r12 # lensing weight of mass plane to be added - f = 3 * self.cosmo.omega_m / 2 + f = 3 * self.cosmo.Omega_m0 / 2 f *= x2 * self.r23 - f *= (1 + self.z2) / self.cosmo.ef(self.z2) + f *= (1 + self.z2) / self.cosmo.H_over_H0(self.z2) f *= w2 # create kappa planes on first iteration @@ -527,7 +534,7 @@ def wlens(self) -> float: def multi_plane_matrix( shells: collections.abc.Sequence[RadialWindow], - cosmo: Cosmology, + cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]], ) -> npt.NDArray[np.float64]: """ Compute the matrix of lensing contributions from each shell. @@ -555,7 +562,7 @@ def multi_plane_matrix( def multi_plane_weights( weights: npt.NDArray[np.float64], shells: collections.abc.Sequence[RadialWindow], - cosmo: Cosmology, + cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]], ) -> npt.NDArray[np.float64]: """ Compute effective weights for multi-plane convergence. diff --git a/glass/shells.py b/glass/shells.py index 4ee414d5..1e0c2e32 100644 --- a/glass/shells.py +++ b/glass/shells.py @@ -53,19 +53,18 @@ import numpy as np import numpy.typing as npt +from cosmology.api import CosmologyConstantsNamespace, StandardCosmology + from glass.core.algorithm import nnls from glass.core.array import ndinterp -if typing.TYPE_CHECKING: - from cosmology import Cosmology - ArrayLike1D = collections.abc.Sequence[float] | npt.NDArray[np.float64] WeightFunc = typing.Callable[[ArrayLike1D], npt.NDArray[np.float64]] def distance_weight( z: npt.NDArray[np.float64], - cosmo: Cosmology, + cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]], ) -> npt.NDArray[np.float64]: """ Uniform weight in comoving distance. @@ -82,12 +81,12 @@ def distance_weight( The weight function evaluated at redshifts *z*. """ - return 1 / cosmo.ef(z) # type: ignore[no-any-return] + return 1 / cosmo.H_over_H0(z) # type: ignore[no-any-return] def volume_weight( z: npt.NDArray[np.float64], - cosmo: Cosmology, + cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]], ) -> npt.NDArray[np.float64]: """ Uniform weight in comoving volume. @@ -104,12 +103,15 @@ def volume_weight( The weight function evaluated at redshifts *z*. """ - return cosmo.xm(z) ** 2 / cosmo.ef(z) # type: ignore[no-any-return] + hubble_length = CosmologyConstantsNamespace.c / StandardCosmology.H0 + return ( # type: ignore[no-any-return] + cosmo.transverse_comoving_distance(z) / hubble_length + ) ** 2 / cosmo.H_over_H0(z) def density_weight( z: npt.NDArray[np.float64], - cosmo: Cosmology, + cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]], ) -> npt.NDArray[np.float64]: """ Uniform weight in matter density. @@ -126,7 +128,13 @@ def density_weight( 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] + hubble_length = CosmologyConstantsNamespace.c / StandardCosmology.H0 + return ( # type: ignore[no-any-return] + cosmo.critical_density0 + * cosmo.Omega_m(z) + * (cosmo.transverse_comoving_distance(z) / hubble_length) ** 2 + / cosmo.H_over_H0(z) + ) class RadialWindow(typing.NamedTuple): @@ -777,7 +785,7 @@ def redshift_grid( def distance_grid( - cosmo: Cosmology, + cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]], zmin: float, zmax: float, *, @@ -805,7 +813,7 @@ def distance_grid( The redshift grid. """ - xmin, xmax = cosmo.dc(zmin), cosmo.dc(zmax) + xmin, xmax = cosmo.comoving_distance(zmin), cosmo.comoving_distance(zmax) x = _uniform_grid(xmin, xmax, step=dx, num=num) return cosmo.dc_inv(x) # type: ignore[no-any-return] diff --git a/pyproject.toml b/pyproject.toml index 3c292882..be67c35d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,7 +20,7 @@ classifiers = [ "Typing :: Typed", ] dependencies = [ - "cosmology>=2022.10.9", + "cosmology.api>=0.1.0", "gaussiancl>=2022.10.21", "healpix>=2022.11.1", "healpy>=1.15.0", @@ -50,6 +50,7 @@ docs = [ ] examples = [ "camb", + "cosmology", "glass.ext.camb", "jupyter", "matplotlib", @@ -154,6 +155,7 @@ lint.isort = {known-first-party = [ ], sections = {"cosmo" = [ "camb", "cosmology", + "cosmology.api", ]}} lint.per-file-ignores = {"__init__.py" = [ "F401", # unused-import diff --git a/tests/conftest.py b/tests/conftest.py index 844b432e..7e9640fa 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -8,7 +8,7 @@ import packaging.version import pytest -from cosmology import Cosmology +from cosmology import StandardCosmology from glass import RadialWindow @@ -98,13 +98,12 @@ def _import_and_add_jax(xp_available_backends: dict[str, types.ModuleType]) -> N # use this as a decorator for tests involving array API compatible functions array_api_compatible = pytest.mark.parametrize("xp", xp_available_backends.values()) - # Pytest fixtures @pytest.fixture(scope="session") -def cosmo() -> Cosmology: +def cosmo() -> StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]]: class MockCosmology: @property - def omega_m(self) -> float: + def Omega_m0(self) -> float: # noqa: N802 """Matter density parameter at redshift 0.""" return 0.3 @@ -113,9 +112,9 @@ def rho_c(self) -> float: """Critical density at redshift 0 in Msol Mpc-3.""" return 3e4 - def ef(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + def H_over_H0(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: # noqa: N802 """Standardised Hubble function :math:`E(z) = H(z)/H_0`.""" - return (self.omega_m * (1 + z) ** 3 + 1 - self.omega_m) ** 0.5 + return (self.Omega_m0 * (1 + z) ** 3 + 1 - self.Omega_m0) ** 0.5 def xm( self, @@ -133,7 +132,7 @@ def xm( def rho_m_z(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: """Redshift-dependent matter density in Msol Mpc-3.""" - return self.rho_c * self.omega_m * (1 + z) ** 3 + return self.rho_c * self.Omega_m0 * (1 + z) ** 3 def dc( self, diff --git a/tests/test_lensing.py b/tests/test_lensing.py index 54487aae..00928c70 100644 --- a/tests/test_lensing.py +++ b/tests/test_lensing.py @@ -4,6 +4,7 @@ import healpix import numpy as np +import numpy.typing as npt import pytest from glass import ( @@ -17,7 +18,7 @@ ) if typing.TYPE_CHECKING: - from cosmology import Cosmology + from cosmology.api import StandardCosmology def test_from_convergence(rng: np.random.Generator) -> None: @@ -63,7 +64,7 @@ def test_shear_from_convergence() -> None: def test_multi_plane_matrix( shells: list[RadialWindow], - cosmo: Cosmology, + cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]], rng: np.random.Generator, ) -> None: mat = multi_plane_matrix(shells, cosmo) @@ -85,7 +86,7 @@ def test_multi_plane_matrix( def test_multi_plane_weights( shells: list[RadialWindow], - cosmo: Cosmology, + cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]], rng: np.random.Generator, ) -> None: w_in = np.eye(len(shells))