Skip to content

Commit

Permalink
Merge branch 'issue/591/triaxiality' of https://github.com/LSSTDESC/CLMM
Browse files Browse the repository at this point in the history 
 into issue/591/triaxiality
  • Loading branch information
@tae-h-shin
tae-h-shin committed Jul 18, 2024
2 parents ce51076 + 3e90071 commit 243bf0f
Showing 1 changed file with 148 additions and 113 deletions.
261 changes: 148 additions & 113 deletions clmm/theory/parent_class.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,6 @@
compute_profile_mass_in_radius,
convert_profile_mass_concentration,
)
from .func_layer import (
compute_surface_density,
)
from ..utils import (
validate_argument,
compute_beta_s_func,
Expand Down Expand Up @@ -133,27 +130,21 @@ def cdelta(self, cdelta):
def massdef(self, massdef):
"""Set definition for the mass of cluster"""
self.set_halo_density_profile(
halo_profile_model=self.halo_profile_model,
massdef=massdef,
delta_mdef=self.delta_mdef,
halo_profile_model=self.halo_profile_model, massdef=massdef, delta_mdef=self.delta_mdef,
)

@delta_mdef.setter
def delta_mdef(self, delta_mdef):
"""Set number of deltas in mass definition of cluster"""
self.set_halo_density_profile(
halo_profile_model=self.halo_profile_model,
massdef=self.massdef,
delta_mdef=delta_mdef,
halo_profile_model=self.halo_profile_model, massdef=self.massdef, delta_mdef=delta_mdef,
)

@halo_profile_model.setter
def halo_profile_model(self, halo_profile_model):
"""Set halo profile model"""
self.set_halo_density_profile(
halo_profile_model=halo_profile_model,
massdef=self.massdef,
delta_mdef=self.delta_mdef,
halo_profile_model=halo_profile_model, massdef=self.massdef, delta_mdef=self.delta_mdef,
)

# 2. Functions to be implemented by children classes
Expand Down Expand Up @@ -236,7 +227,7 @@ def __integrand__(l_value, theta):

l_values = np.logspace(loglbounds[0], loglbounds[1], lsteps)
kernel = np.array([simpson(__integrand__(l_values, t), x=l_values) for t in theta])
return halobias * kernel * rho_m / (2 * np.pi * (1 + z_cl) ** 3 * da**2)
return halobias * kernel * rho_m / (2 * np.pi * (1 + z_cl) ** 3 * da ** 2)

def _eval_surface_density_2h(
self,
Expand Down Expand Up @@ -268,78 +259,50 @@ def _eval_excess_surface_density_2h(
2, r_proj, z_cl, halobias, logkbounds, ksteps, loglbounds, lsteps
)

def _eval_excess_surface_density_triaxial4theta(
self, r_proj, z_cl, ell
):

sigma = compute_surface_density(
r_proj,
self.mdelta,
self.cdelta,
z_cl,
self.cosmo,
self.delta_mdef,
self.halo_profile_model,
self.massdef,
)
def _eval_excess_surface_density_triaxial4theta(self, r_proj, z_cl, ell):
"""eval triaxial-4theta term of excess surface density"""
sigma = self.eval_surface_density(r_proj, z_cl)

eta = r_proj * np.gradient(np.log(sigma), r_proj)
f = InterpolatedUnivariateSpline(r_proj, r_proj**3 * sigma * eta, k=5)
f = InterpolatedUnivariateSpline(r_proj, r_proj ** 3 * sigma * eta, k=5)
integral_vec = np.vectorize(f.integral)

I = 3 / r_proj**4 * integral_vec(0, r_proj)
I = 3 / r_proj ** 4 * integral_vec(0, r_proj)

return 0.5 * ell * (2 * I - sigma * eta) / 1e12
return 0.5 * ell * (2 * I - sigma * eta)

def _eval_excess_surface_density_triaxialConst(
self, r_proj, z_cl, ell
):

sigma = compute_surface_density(
r_proj,
self.mdelta,
self.cdelta,
z_cl,
self.cosmo,
self.delta_mdef,
self.halo_profile_model,
self.massdef,
)
def _eval_excess_surface_density_triaxialConst(self, r_proj, z_cl, ell):
"""eval triaxial-const term of excess surface density"""
sigma = self.eval_surface_density(r_proj, z_cl)

eta = r_proj * np.gradient(np.log(sigma), r_proj)
f = InterpolatedUnivariateSpline(r_proj, sigma * eta / r_proj, k=5)
integral_vec = np.vectorize(f.integral)

I = integral_vec(r_proj, np.inf)

return 0.5 * ell * (2 * I - sigma * eta) / 1e12
return 0.5 * ell * (2 * I - sigma * eta)

def _eval_triaxial_corrected_surface_density(
self, r_proj, z_cl, ell
):
def _eval_triaxial_corrected_surface_density(self, r_proj, z_cl, ell):
"""eval surface density with triaxial corrections"""
sigma = self.eval_surface_density(r_proj, z_cl)

sigma = compute_surface_density(
r_proj,
self.mdelta,
self.cdelta,
z_cl,
self.cosmo,
self.delta_mdef,
self.halo_profile_model,
self.massdef,
)
eta = r_proj * np.gradient(np.log(sigma), r_proj)
deta_dlnr = r_proj * np.gradient(eta, r_proj)

return sigma * (1 + ell**2 * (eta + eta**2/2 + deta_dlnr/2) / 2)
return sigma * (1 + ell ** 2 * (eta + eta ** 2 / 2 + deta_dlnr / 2) / 2)

def _eval_triaxial_corrected_excess_surface_density_monopole(
self, r_proj, z_cl, ell
):
def _eval_triaxial_corrected_excess_surface_density_monopole(self, r_proj, z_cl, ell):
"""eval triaxial-monopole term of excess surface density"""
sigma = self._eval_triaxial_corrected_surface_density(r_proj, z_cl, ell)
integral = np.array([simpson(r_proj[r_proj <= ri] * sigma[r_proj<=ri], r_proj[r_proj<=ri]) for ri in r_proj])
integral = np.array(
[
simpson(r_proj[r_proj <= ri] * sigma[r_proj <= ri], r_proj[r_proj <= ri])
for ri in r_proj
]
)

return ((2/r_proj**2) * integral - sigma) / sigma_crit
return (2 / r_proj ** 2) * integral - sigma

def _eval_rdelta(self, z_cl):
return compute_rdelta(self.mdelta, z_cl, self.cosmo, self.massdef, self.delta_mdef)
Expand Down Expand Up @@ -783,54 +746,133 @@ def eval_surface_density_2h(
r_proj, z_cl, halobias, logkbounds, ksteps, loglbounds, lsteps
)

def eval_excess_surface_density_triaxial4theta(
self, r_proj, z_cl, ell
):
r"""DOCUMENTATION RANCE
def eval_excess_surface_density_triaxial4theta(self, r_proj, z_cl, ell):
r"""Compute the "4-theta"-quadrupole component
Parameters
----------
r_proj: array
Projected radial position from the cluster center in :math:`M\!pc`.
ell: float
ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
q=b/a (Ratio of major axis to the minor axis lengths)
z_cl: float
Redshift of lens cluster
Returns
-------
numpy.ndarry, float
Triaxial-4theta component of the excess surface density in units of :math:`M_\odot\ Mpc^{-2}`.
"""
## VALIDATE ARGUMENTS
## VALIDATE BACKEND

return self._eval_excess_surface_density_triaxial4theta(
r_proj, z_cl, ell
)
if self.validate_input:
validate_argument(locals(), "r_proj", "float_array", argmin=0)
validate_argument(locals(), "z_cl", float, argmin=0)
validate_argument(locals(), "ell", float, argmin=0, argmax=1)

def eval_excess_surface_density_triaxialConst(
self, r_proj, z_cl, ell
):
r"""DOCUMENTATION RANCE
if self.backend not in ("ccl", "nc"):
raise NotImplementedError(
f"Triaxial-4theta term not currently supported with the {self.backend} backend. "
"Use the CCL or NumCosmo backend instead"
)

return self._eval_excess_surface_density_triaxial4theta(r_proj, z_cl, ell)

def eval_excess_surface_density_triaxialConst(self, r_proj, z_cl, ell):
r"""Compute the "constant"-quadrupole component
Parameters
----------
r_proj: array
Projected radial position from the cluster center in :math:`M\!pc`.
ell: float
ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
q=b/a (Ratio of major axis to the minor axis lengths)
z_cl: float
Redshift of lens cluster
Returns
-------
numpy.ndarry, float
Triaxial-constant component of the excess surface density in units of :math:`M_\odot\ Mpc^{-2}`.
"""
## VALIDATE ARGUMENTS
## VALIDATE BACKEND

return self._eval_excess_surface_density_triaxialConst(
r_proj, z_cl, ell
)
if self.validate_input:
validate_argument(locals(), "r_proj", "float_array", argmin=0)
validate_argument(locals(), "z_cl", float, argmin=0)
validate_argument(locals(), "ell", float, argmin=0, argmax=1)

def eval_triaxial_corrected_surface_density(
self, r_proj, z_cl, ell
):
r"""DOCUMENTATION RANCE
if self.backend not in ("ccl", "nc"):
raise NotImplementedError(
f"Triaxial-constant term not currently supported with the {self.backend} backend. "
"Use the CCL or NumCosmo backend instead"
)

return self._eval_excess_surface_density_triaxialConst(r_proj, z_cl, ell)

def eval_triaxial_corrected_surface_density(self, r_proj, z_cl, ell):
r"""Compute surface density with triaxial corrections
Parameters
----------
r_proj: array
Projected radial position from the cluster center in :math:`M\!pc`.
ell: float
ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
q=b/a (Ratio of major axis to the minor axis lengths)
z_cl: float
Redshift of lens cluster
Returns
-------
numpy.ndarry, float
Triaxial corrected surface density in units of :math:`M_\odot\ Mpc^{-2}`.
"""
## VALIDATE ARGUMENTS
## VALIDATE BACKEND

return self._eval_triaxial_corrected_surface_density(
r_proj, z_cl, ell
)
if self.validate_input:
validate_argument(locals(), "r_proj", "float_array", argmin=0)
validate_argument(locals(), "z_cl", float, argmin=0)
validate_argument(locals(), "ell", float, argmin=0, argmax=1)

def eval_triaxial_corrected_excess_surface_density_monopole(
self, r_proj, z_cl, ell
):
r"""DOCUMENTATION RANCE
if self.backend not in ("ccl", "nc"):
raise NotImplementedError(
f"Triaxial correction to surface density not currently supported with the {self.backend} backend. "
"Use the CCL or NumCosmo backend instead"
)

return self._eval_triaxial_corrected_surface_density(r_proj, z_cl, ell)

def eval_triaxial_corrected_excess_surface_density_monopole(self, r_proj, z_cl, ell):
r"""Compute monopole component
Parameters
----------
r_proj: array
Projected radial position from the cluster center in :math:`M\!pc`.
ell: float
ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
q=b/a (Ratio of major axis to the minor axis lengths)
z_cl: float
Redshift of lens cluster
Returns
-------
numpy.ndarry, float
Monopole component of excess surface density in units of :math:`M_\odot\ Mpc^{-2}`.
"""
## VALIDATE ARGUMENTS
## VALIDATE BACKEND

return self._eval_triaxial_corrected_excess_surface_density_monopole(
r_proj, z_cl, ell
)
if self.validate_input:
validate_argument(locals(), "r_proj", "float_array", argmin=0)
validate_argument(locals(), "z_cl", float, argmin=0)
validate_argument(locals(), "ell", float, argmin=0, argmax=1)

if self.backend not in ("ccl", "nc"):
raise NotImplementedError(
f"Triaxial correction to surface density not currently supported with the {self.backend} backend. "
"Use the CCL or NumCosmo backend instead"
)

return self._eval_triaxial_corrected_excess_surface_density_monopole(r_proj, z_cl, ell)

def eval_tangential_shear(self, r_proj, z_cl, z_src, z_src_info="discrete", verbose=False):
r"""Computes the tangential shear
Expand Down Expand Up @@ -981,14 +1023,7 @@ def eval_convergence(self, r_proj, z_cl, z_src, z_src_info="discrete", verbose=F
+ "\nConvergence = 0 for those galaxies."
)
kappa = compute_for_good_redshifts(
self._eval_convergence_core,
z_cl,
z_src,
0.0,
warning_msg,
"z_cl",
"z_src",
r_proj,
self._eval_convergence_core, z_cl, z_src, 0.0, warning_msg, "z_cl", "z_src", r_proj,
)
elif z_src_info == "beta":
beta_s_mean = z_src[0]
Expand Down Expand Up @@ -1351,7 +1386,7 @@ def eval_magnification(
if approx is None:
if z_src_info == "distribution":
mu = self._pdz_weighted_avg(
lambda gammat, kappa: 1 / ((1 - kappa) ** 2 - gammat**2),
lambda gammat, kappa: 1 / ((1 - kappa) ** 2 - gammat ** 2),
z_src,
r_proj,
z_cl,
Expand Down Expand Up @@ -1383,7 +1418,7 @@ def eval_magnification(
if approx == "order2":
beta_s_square_mean = z_src[1]
# Taylor expansion with up to second-order terms
mu += 3 * beta_s_square_mean * kappa_inf**2 + beta_s_square_mean * gammat_inf**2
mu += 3 * beta_s_square_mean * kappa_inf ** 2 + beta_s_square_mean * gammat_inf ** 2

return mu

Expand Down Expand Up @@ -1519,7 +1554,7 @@ def eval_magnification_bias(
# z_src (float or array) is redshift
if z_src_info == "distribution":
mu_bias = self._pdz_weighted_avg(
lambda gammat, kappa: 1 / ((1 - kappa) ** 2 - gammat**2) ** (alpha - 1),
lambda gammat, kappa: 1 / ((1 - kappa) ** 2 - gammat ** 2) ** (alpha - 1),
z_src,
r_proj,
z_cl,
Expand Down Expand Up @@ -1553,9 +1588,9 @@ def eval_magnification_bias(
if approx == "order2":
beta_s_square_mean = z_src[1]
# Taylor expansion with up to second-order terms
mu_bias += (alpha - 1) * (beta_s_square_mean * gammat_inf**2) + (
mu_bias += (alpha - 1) * (beta_s_square_mean * gammat_inf ** 2) + (
2 * alpha - 1
) * (alpha - 1) * beta_s_square_mean * kappa_inf**2
) * (alpha - 1) * beta_s_square_mean * kappa_inf ** 2

return mu_bias

Expand Down

0 comments on commit 243bf0f

Please sign in to comment.