Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Draft
wants to merge 23 commits into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from 11 commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion glass/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,6 @@
from glass.galaxies import (
galaxy_shear,
gaussian_phz,
kappa_ia_nla,
redshifts,
redshifts_from_nz,
)
Expand Down
106 changes: 0 additions & 106 deletions glass/galaxies.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,6 @@
.. autofunction:: galaxy_shear
.. autofunction:: gaussian_phz

Intrinsic alignments
--------------------
.. autofunction:: kappa_ia_nla
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have reverted #21 as I believe it is broken/needs testing. Don't want to waste time porting.


""" # noqa: D205, D400

from __future__ import annotations
Expand All @@ -33,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


Expand Down Expand Up @@ -294,103 +288,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.

Returns the effective convergence due to intrinsic alignments.

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.

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]
23 changes: 16 additions & 7 deletions glass/lensing.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@
if typing.TYPE_CHECKING:
import collections.abc

from cosmology import Cosmology
import cosmology.api

from glass.shells import RadialWindow

Expand Down Expand Up @@ -269,7 +269,12 @@ def shear_from_convergence(
class MultiPlaneConvergence:
"""Compute convergence fields iteratively from multiple matter planes."""

def __init__(self, cosmo: Cosmology) -> None:
def __init__(
self,
cosmo: cosmology.api.StandardCosmology[
npt.NDArray[np.float64], npt.NDArray[np.float64]
],
) -> None:
"""Create a new instance to iteratively compute the convergence."""
self.cosmo = cosmo

Expand Down Expand Up @@ -301,7 +306,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)
Expand Down Expand Up @@ -333,9 +338,9 @@ def add_plane(
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
Expand Down Expand Up @@ -377,7 +382,9 @@ def wlens(self) -> float:

def multi_plane_matrix(
shells: collections.abc.Sequence[RadialWindow],
cosmo: Cosmology,
cosmo: cosmology.api.StandardCosmology[
npt.NDArray[np.float64], npt.NDArray[np.float64]
],
) -> npt.NDArray[np.float64]:
"""Compute the matrix of lensing contributions from each shell."""
mpc = MultiPlaneConvergence(cosmo)
Expand All @@ -391,7 +398,9 @@ def multi_plane_matrix(
def multi_plane_weights(
weights: npt.NDArray[np.float64],
shells: collections.abc.Sequence[RadialWindow],
cosmo: Cosmology,
cosmo: cosmology.api.StandardCosmology[
npt.NDArray[np.float64], npt.NDArray[np.float64]
],
) -> npt.NDArray[np.float64]:
"""
Compute effective weights for multi-plane convergence.
Expand Down
31 changes: 22 additions & 9 deletions glass/shells.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,34 +55,45 @@
from glass.core.array import ndinterp

if typing.TYPE_CHECKING:
from cosmology import Cosmology
import cosmology.api

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: cosmology.api.StandardCosmology[
npt.NDArray[np.float64], npt.NDArray[np.float64]
],
) -> npt.NDArray[np.float64]:
"""Uniform weight in comoving distance."""
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: cosmology.api.StandardCosmology[
npt.NDArray[np.float64], npt.NDArray[np.float64]
],
) -> npt.NDArray[np.float64]:
"""Uniform weight in comoving volume."""
return cosmo.xm(z) ** 2 / cosmo.ef(z) # type: ignore[no-any-return]
return cosmo.xm(z) ** 2 / cosmo.H_over_H0(z) # type: ignore[no-any-return]


def density_weight(
z: npt.NDArray[np.float64],
cosmo: Cosmology,
cosmo: cosmology.api.StandardCosmology[
npt.NDArray[np.float64], npt.NDArray[np.float64]
],
) -> npt.NDArray[np.float64]:
"""Uniform weight in matter density."""
return cosmo.rho_m_z(z) * cosmo.xm(z) ** 2 / cosmo.ef(z) # type: ignore[no-any-return]
return ( # type: ignore[no-any-return]
cosmo.critical_density0
* cosmo.Omega_m(z)
* cosmo.xm(z) ** 2
/ cosmo.H_over_H0(z)
)


class RadialWindow(typing.NamedTuple):
Expand Down Expand Up @@ -650,15 +661,17 @@ def redshift_grid(


def distance_grid(
cosmo: Cosmology,
cosmo: cosmology.api.StandardCosmology[
npt.NDArray[np.float64], npt.NDArray[np.float64]
],
zmin: float,
zmax: float,
*,
dx: float | None = None,
num: int | None = None,
) -> npt.NDArray[np.float64]:
"""Redshift grid with uniform spacing in comoving distance."""
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:
Expand Down
4 changes: 3 additions & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -50,6 +50,7 @@ docs = [
]
examples = [
"camb",
"cosmology",
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cosmology is still used in the notebooks currently

"glass.ext.camb",
"jupyter",
"matplotlib",
Expand Down Expand Up @@ -156,6 +157,7 @@ lint.isort = {known-first-party = [
], sections = {"cosmo" = [
"camb",
"cosmology",
"cosmology.api",
]}}
lint.per-file-ignores = {"__init__.py" = [
"F401", # unused-import
Expand Down
20 changes: 13 additions & 7 deletions tests/test_lensing.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
from glass.shells import RadialWindow

if typing.TYPE_CHECKING:
from cosmology import Cosmology
import cosmology.api


@pytest.fixture
Expand All @@ -31,14 +31,16 @@ def shells() -> list[RadialWindow]:


@pytest.fixture
def cosmo() -> Cosmology:
def cosmo() -> (
cosmology.api.StandardCosmology[npt.NDArray[np.float64], npt.NDArray[np.float64]]
):
class MockCosmology:
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've changed the mock to match the names from cosmology.api, hence the noqa

@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(
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is no equivalently named function in cosmology.api. In practice, this would be

cosmology.api.Standard.Cosmology.transverse_comoving_distance(z) / (cosmology.api.CosmologyConstantsNamespace.c / cosmology.api.Standard.Cosmology.H0)

self,
Expand Down Expand Up @@ -97,7 +99,9 @@ def test_deflect_many(rng: np.random.Generator) -> None:

def test_multi_plane_matrix(
shells: list[RadialWindow],
cosmo: Cosmology,
cosmo: cosmology.api.StandardCosmology[
npt.NDArray[np.float64], npt.NDArray[np.float64]
],
rng: np.random.Generator,
) -> None:
mat = multi_plane_matrix(shells, cosmo)
Expand All @@ -119,7 +123,9 @@ def test_multi_plane_matrix(

def test_multi_plane_weights(
shells: list[RadialWindow],
cosmo: Cosmology,
cosmo: cosmology.api.StandardCosmology[
npt.NDArray[np.float64], npt.NDArray[np.float64]
],
rng: np.random.Generator,
) -> None:
w_in = np.eye(len(shells))
Expand Down