From efa4154a36db2ae1f2a1ced3dd4dad1250be5ec9 Mon Sep 17 00:00:00 2001 From: Hsin Fan <57552401+hsinfan1996@users.noreply.github.com> Date: Sat, 20 Jul 2024 06:02:01 +0800 Subject: [PATCH] Fix mean_surface_density; Tweak tests --- clmm/theory/parent_class.py | 27 ++++++++-------- tests/test_theory.py | 64 +++++++++++++++++++++++++------------ 2 files changed, 56 insertions(+), 35 deletions(-) diff --git a/clmm/theory/parent_class.py b/clmm/theory/parent_class.py index 433c56214..7ee3ee4b4 100644 --- a/clmm/theory/parent_class.py +++ b/clmm/theory/parent_class.py @@ -339,8 +339,7 @@ def _eval_mean_surface_density_miscentered(self, r_proj, z_cl, r_mis, backend): if backend: integrand = self._integrand_mean_surface_density_mis res[i] = dblquad(integrand, r_lower, r, 0, np.pi, - args=(r_mis, z_cl) - )[0] + args=(r_mis, z_cl))[0] else: c = self.cdelta @@ -350,16 +349,17 @@ def _eval_mean_surface_density_miscentered(self, r_proj, z_cl, r_mis, backend): if self.halo_profile_model == "nfw": rho_s = (self.delta_mdef / 3.0 * c**3.0 * rho_def - / (np.log(1.0 + c) - c / (1.0 + c))) + / (np.log(1.0 + c) - c / (1.0 + c)) + ) integrand = self._integrand_mean_surface_density_mis_nfw - res[i] = (dblquad(integrand, r_lower, r, 0, np.pi, + res[i] = ( + dblquad(integrand, r_lower, r, 0, np.pi, args=(r_mis, r_s), epsrel=1e-6)[0] * 2 * r_s * rho_s - / np.pi ) # Einasto @@ -381,12 +381,11 @@ def _eval_mean_surface_density_miscentered(self, r_proj, z_cl, r_mis, backend): integrand = self._integrand_mean_surface_density_mis_einasto - res[i] = (tplquad(integrand, r_lower, r, 0, np.pi, 0, np.inf, - args=(r_mis, r_s, alpha_ein), epsrel=1e-6 - )[0] + res[i] = ( + tplquad(integrand, r_lower, r, 0, np.pi, 0, np.inf, + args=(r_mis, r_s, alpha_ein), epsrel=1e-6)[0] * 2 * rho_s - / np.pi ) # Hernquist @@ -394,18 +393,18 @@ def _eval_mean_surface_density_miscentered(self, r_proj, z_cl, r_mis, backend): rho_s = self.delta_mdef / 3.0 * c**3.0 * rho_def / ((c / (1.0 + c)) ** 2.0) * 2 integrand = self._integrand_mean_surface_density_mis_hernquist - res[i] = (dblquad(integrand, r_lower, r, 0, np.pi, - args=(r_mis, r_s), epsrel=1e-6)[0] + res[i] = ( + dblquad(integrand, r_lower, r, 0, np.pi, + args=(r_mis, r_s), epsrel=1e-6)[0] * r_s * rho_s - / np.pi ) - res = np.cumsum(res) * 2 / r_proj**2 + res = np.cumsum(res) * 2 / np.pi / r_proj**2 return res - def _integrand_mean_surface_density_mis(self, theta, r, z_cl, r_mis): + def _integrand_mean_surface_density_mis(self, theta, r, r_mis, z_cl): # pylint: disable=invalid-name return r * self._integrand_surface_density_mis(theta, r, r_mis, z_cl) diff --git a/tests/test_theory.py b/tests/test_theory.py index 78c630168..a44d00b69 100644 --- a/tests/test_theory.py +++ b/tests/test_theory.py @@ -323,13 +323,6 @@ def test_profiles(modeling_data, profile_init): cfg["numcosmo_profiles"]["Sigma"], reltol, ) - assert_allclose( - mod.eval_surface_density( - cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, verbose=True - )[-40:], - cfg["numcosmo_profiles"]["Sigma"][-40:], - 2.5e-2, - ) assert_allclose( mod.eval_mean_surface_density( cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], verbose=True @@ -337,13 +330,6 @@ def test_profiles(modeling_data, profile_init): cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"], reltol, ) - assert_allclose( - mod.eval_mean_surface_density( - cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, verbose=True - )[-40:], - (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:], - 8.5e-3, - ) assert_allclose( mod.eval_excess_surface_density( cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], verbose=True @@ -351,13 +337,49 @@ def test_profiles(modeling_data, profile_init): cfg["numcosmo_profiles"]["DeltaSigma"], reltol, ) - assert_allclose( - mod.eval_excess_surface_density( - cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, verbose=True - )[-40:], - cfg["numcosmo_profiles"]["DeltaSigma"][-40:], - 3e-2 - ) + if mod.backend == "nc": + assert_allclose( + mod.eval_surface_density( + cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, + verbose=True)[-40:], + cfg["numcosmo_profiles"]["Sigma"][-40:], + 2.5e-2, + ) + assert_allclose( + mod.eval_mean_surface_density( + cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, + verbose=True)[-40:], + (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:], + 8.5e-3, + ) + assert_allclose( + mod.eval_excess_surface_density( + cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, + verbose=True)[-40:], + cfg["numcosmo_profiles"]["DeltaSigma"][-40:], + 3e-2 + ) + assert_allclose( + mod.eval_surface_density( + cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, + verbose=True, backend=True)[-40:], + cfg["numcosmo_profiles"]["Sigma"][-40:], + 2.5e-2, + ) + assert_allclose( + mod.eval_mean_surface_density( + cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, + verbose=True, backend=True)[-40:], + (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:], + 8.5e-3, + ) + assert_allclose( + mod.eval_excess_surface_density( + cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, + verbose=True, backend=False)[-40:], + cfg["numcosmo_profiles"]["DeltaSigma"][-40:], + 2.5e-2 + ) if mod.backend == "ct": assert_raises( ValueError, mod.eval_excess_surface_density, 1e-12, cfg["SIGMA_PARAMS"]["z_cl"]