Fix mean_surface_density; Tweak tests

LSSTDESC · Jul 19, 2024 · efa4154 · efa4154
1 parent c807e67
commit efa4154
Show file tree

Hide file tree

Showing 2 changed files with 56 additions and 35 deletions.
diff --git 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,31 +381,30 @@ 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
                 elif self.halo_profile_model == "hernquist":
                     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
@@ -323,41 +323,63 @@ 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
             ),
             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
             ),
             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"]