gh-124: use cosmology.api over the cosmology package

glass/__init__.py

@@ -19,7 +19,6 @@
     transform_cls,
 )
 from glass.galaxies import (
-    _kappa_ia_nla,
     galaxy_shear,
     gaussian_phz,
     redshifts,

glass/galaxies.py

@@ -29,8 +29,6 @@
 from glass.core.array import broadcast_leading_axes, cumtrapz
 
 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]

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
 
 
@@ -286,7 +286,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.
 
@@ -333,7 +336,7 @@ def add_window(self, delta: npt.NDArray[np.float64], w: RadialWindow) -> None:
                 zsrc,
                 w.za,
                 w.wa,
-            )
+            ),
         )
 
         self.add_plane(delta, zsrc, lens_weight)
@@ -376,15 +379,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
@@ -426,7 +433,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.
@@ -454,7 +461,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.

glass/shells.py

@@ -52,19 +52,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 = typing.Union[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.
@@ -81,12 +80,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.
@@ -103,12 +102,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.
@@ -125,7 +127,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):
@@ -793,7 +801,7 @@ def redshift_grid(
 
 
 def distance_grid(
-    cosmo: Cosmology,
+    cosmo: StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]],
     zmin: float,
     zmax: float,
     *,
@@ -826,7 +834,7 @@ def distance_grid(
         If both ``dx`` and ``num`` are given.
 
     """
-    xmin, xmax = cosmo.dc(zmin), cosmo.dc(zmax)
+    xmin, xmax = cosmo.comoving_distance(zmin), cosmo.comoving_distance(zmax)
     if dx is not None and num is None:
         x = np.arange(xmin, np.nextafter(xmax + dx, xmax), dx)
     elif dx is None and num is not None:

pyproject.toml

@@ -21,7 +21,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",
@@ -51,6 +51,7 @@ docs = [
 ]
 examples = [
     "camb",
+    "cosmology",
     "glass.ext.camb",
     "jupyter",
     "matplotlib",
@@ -157,6 +158,7 @@ lint.isort = {known-first-party = [
 ], sections = {"cosmo" = [
     "camb",
     "cosmology",
+    "cosmology.api",
 ]}}
 lint.per-file-ignores = {"__init__.py" = [
     "F401", # unused-import

tests/test_lensing.py

@@ -16,7 +16,7 @@
 from glass.shells import RadialWindow
 
 if typing.TYPE_CHECKING:
-    from cosmology import Cosmology
+    from cosmology.api import StandardCosmology
 
 
 @pytest.fixture
@@ -31,20 +31,21 @@ def shells() -> list[RadialWindow]:
 
 
 @pytest.fixture
-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
             return 0.3
 
-        def ef(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
-            return (self.omega_m * (1 + z) ** 3 + 1 - self.omega_m) ** 0.5
+        def H_over_H0(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:  # noqa: N802
+            return (self.Omega_m0 * (1 + z) ** 3 + 1 - self.Omega_m0) ** 0.5
 
         def xm(
             self,
             z: npt.NDArray[np.float64],
             z2: npt.NDArray[np.float64] | None = None,
         ) -> npt.NDArray[np.float64]:
+            """Dimensionless transverse comoving distance."""
             if z2 is None:
                 return np.array(z) * 1000
             return (np.array(z2) - np.array(z)) * 1000
@@ -97,7 +98,7 @@ def test_deflect_many(rng: np.random.Generator) -> 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)
@@ -119,7 +120,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))