Skip to content

Commit

Permalink
Allow backend-dependent calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
@hsinfan1996
hsinfan1996 committed Jul 17, 2024
1 parent e55db86 commit 00b0069
Showing 1 changed file with 60 additions and 33 deletions.
93 changes: 60 additions & 33 deletions clmm/theory/parent_class.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,43 +205,51 @@ def _set_cosmo(self, cosmo):
r"""Sets the cosmology to the internal cosmology object"""
self.cosmo = cosmo if cosmo is not None else self.cosmo_class()

def _eval_surface_density_miscentered(self, r_proj, z_cl, r_mis):
def _eval_surface_density_miscentered(self, r_proj, z_cl, r_mis, backend):
# pylint: disable=invalid-name, possibly-used-before-assignment
c = self.cdelta
rho_def = self.cosmo.get_rho_m(z_cl)
r_s = self.eval_rdelta(z_cl) / c

if self.halo_profile_model=='nfw':
rho_s = self.delta_mdef / 3. * c**3. * rho_def / (np.log(1.+c) - c/(1.+c))
integrand = self._integrand_surface_density_mis_NFW
if backend:
integrand = self._integrand_surface_density_mis
res = quad_vec(integrand, 0., np.pi, args=(r_proj, r_mis, z_cl), epsrel=1e-6)[0]/np.pi

res = (quad_vec(integrand, 0., np.pi, args=(r_proj, r_mis, r_s), epsrel=1e-6)[0]
* 2 * r_s * rho_s / np.pi)
else:
if self.halo_profile_model=='nfw':
rho_s = self.delta_mdef / 3. * c**3. * rho_def / (np.log(1.+c) - c/(1.+c))
integrand = self._integrand_surface_density_mis_NFW

res = (quad_vec(integrand, 0., np.pi, args=(r_proj, r_mis, r_s), epsrel=1e-6)[0]
* 2 * r_s * rho_s / np.pi)

# Einasto
elif self.halo_profile_model=='einasto':
alpha_ein = self._get_einasto_alpha(z_cl)
rho_s = self.delta_mdef/3.*c**3.*rho_def/(
2.**(-3./alpha_ein) * alpha_ein**(-1.+3./alpha_ein)
* np.exp(2./alpha_ein) * gamma(3./alpha_ein)
* gammainc(3./alpha_ein, 2./alpha_ein*c**alpha_ein))
# Einasto
elif self.halo_profile_model=='einasto':
alpha_ein = self._get_einasto_alpha(z_cl)
rho_s = self.delta_mdef/3.*c**3.*rho_def/(
2.**(-3./alpha_ein) * alpha_ein**(-1.+3./alpha_ein)
* np.exp(2./alpha_ein) * gamma(3./alpha_ein)
* gammainc(3./alpha_ein, 2./alpha_ein*c**alpha_ein))

integrand = self._integrand_surface_density_mis_Einasto
integrand = self._integrand_surface_density_mis_Einasto

res = (quad_vec(integrand, 0., np.pi,
args=(r_proj, r_mis, r_s, alpha_ein), epsrel=1e-6)[0]
* 2 * rho_s / np.pi)
res = (quad_vec(integrand, 0., np.pi,
args=(r_proj, r_mis, r_s, alpha_ein), epsrel=1e-6)[0]
* 2 * rho_s / np.pi)

# Hernquist
elif self.halo_profile_model=='hernquist':
rho_s = self.delta_mdef/3.*c**3.*rho_def/((c/(1. + c))**2.)*2
integrand = self._integrand_surface_density_mis_Hernquist
# Hernquist
elif self.halo_profile_model=='hernquist':
rho_s = self.delta_mdef/3.*c**3.*rho_def/((c/(1. + c))**2.)*2
integrand = self._integrand_surface_density_mis_Hernquist

res = (quad_vec(integrand, 0., np.pi, args=(r_proj, r_mis, r_s), epsrel=1e-6)[0]
* r_s * rho_s / np.pi)
res = (quad_vec(integrand, 0., np.pi, args=(r_proj, r_mis, r_s), epsrel=1e-6)[0]
* r_s * rho_s / np.pi)

return res

def _integrand_surface_density_mis(self, theta, R, Roff, z_cl):
return self.eval_surface_density(np.sqrt(R*R + Roff*Roff - 2*R*Roff*np.cos(theta)), z_cl)

def _integrand_surface_density_mis_NFW(self, theta, R, Roff, r_s):
# pylint: disable=invalid-name
x = np.sqrt(R**2. + Roff**2. - 2.*R*Roff*np.cos(theta)) / r_s
Expand Down Expand Up @@ -288,19 +296,19 @@ def f2(x):

return np.piecewise(x, [x<1, x>1], [f1, f2, 4./15.])

def _eval_mean_surface_density_miscentered(self, r_proj, z_cl, r_mis):
def _eval_mean_surface_density_miscentered(self, r_proj, z_cl, r_mis, backend):
res = np.zeros_like(r_proj)
# pylint: disable=invalid-name
for i, R in enumerate(r_proj):
R_lower = 0 if i==0 else r_proj[i-1]
res[i] = quad(self._integrand_mean_surface_density_mis, R_lower, R,
args=(z_cl, r_mis))[0]
args=(z_cl, r_mis, backend))[0]
res = np.cumsum(res)*2/r_proj**2
return res

def _integrand_mean_surface_density_mis(self, R, z_cl, r_mis):
def _integrand_mean_surface_density_mis(self, R, z_cl, r_mis, backend):
# pylint: disable=invalid-name
return R * self._eval_surface_density_miscentered(R, z_cl, r_mis)
return R * self._eval_surface_density_miscentered(R, z_cl, r_mis, backend)

def _eval_excess_surface_density_miscentered(self, r_proj, z_cl, r_mis):
return (self._eval_mean_surface_density_miscentered(r_proj, z_cl, r_mis)
Expand Down Expand Up @@ -630,7 +638,7 @@ def inv_sigmac(redshift):

return 1.0 / _integ_pzfuncs(pzpdf, pzbins, kernel=inv_sigmac)

def eval_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False):
def eval_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False, backend=False):
r"""Computes the surface mass density
Parameters
Expand All @@ -641,6 +649,12 @@ def eval_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False):
Redshift of the cluster
r_mis : float, optional
Projected miscenter distance in :math:`M\!pc`
verbose : bool, optional
If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and
CCL backends. (Default: False)
backend : bool, optional
If True, use the projected surface density from the backend for miscentering
calculations. (Default: False)
Returns
-------
Expand All @@ -658,10 +672,11 @@ def eval_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False):
print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")

if r_mis is not None:
return self._eval_surface_density_miscentered(r_proj=r_proj, z_cl=z_cl, r_mis=r_mis)
return self._eval_surface_density_miscentered(r_proj=r_proj, z_cl=z_cl, r_mis=r_mis,
backend=backend)
return self._eval_surface_density(r_proj=r_proj, z_cl=z_cl)

def eval_mean_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False):
def eval_mean_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False, backend=False):
r"""Computes the mean value of surface density inside radius `r_proj`
Parameters
Expand All @@ -672,6 +687,12 @@ def eval_mean_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False):
Redshift of the cluster
r_mis : float, optional
Projected miscenter distance in :math:`M\!pc`.
verbose : bool, optional
If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and
CCL backends. (Default: False)
backend : bool, optional
If True, use the projected surface density from the backend for miscentering
calculations. (Default: False)
Returns
-------
Expand All @@ -689,11 +710,11 @@ def eval_mean_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False):

if r_mis is not None:
return self._eval_mean_surface_density_miscentered(r_proj=r_proj, z_cl=z_cl,
r_mis=r_mis)
r_mis=r_mis, backend=backend)

return self._eval_mean_surface_density(r_proj=r_proj, z_cl=z_cl)

def eval_excess_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False):
def eval_excess_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False, backend=False):
r"""Computes the excess surface density
Parameters
Expand All @@ -704,6 +725,12 @@ def eval_excess_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False):
Redshift of the cluster
r_mis : float, optional
Projected miscenter distance in :math:`M\!pc`.
verbose : bool, optional
If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and
CCL backends. (Default: False)
backend : bool, optional
If True, use the projected surface density from the backend for miscentering
calculations. (Default: False)
Returns
-------
Expand All @@ -721,7 +748,7 @@ def eval_excess_surface_density(self, r_proj, z_cl, r_mis=None, verbose=False):

if r_mis is not None:
return self._eval_excess_surface_density_miscentered(r_proj=r_proj, z_cl=z_cl,
r_mis=r_mis)
r_mis=r_mis, backend=backend)
return self._eval_excess_surface_density(r_proj=r_proj, z_cl=z_cl)

def eval_excess_surface_density_2h(
Expand Down

0 comments on commit 00b0069

Please sign in to comment.