ENH make sure to call numcosmo properly for benchmarks (#789)

* ENH make sure to call numcosmo propetly for benchmarks

* ENH add a high accuracy benchmark

* REF moved tests for number counts

benchmarks/data/codes/numcosmo_halo_number_counts.py

@@ -66,14 +66,15 @@
 # now compute the mass function
 psf = Ncm.PowspecFilter.new(psml, Ncm.PowspecFilterType.TOPHAT)
 psf.set_best_lnr0()
+psf.set_zf(2.5)
 mulf = Nc.MultiplicityFunc.new_from_name("NcMultiplicityFuncTinkerMean")
 mulf.set_Delta(200)
 mf = Nc.HaloMassFunction.new(dist, psf, mulf)
 
 psf.prepare(cosmo)
 mf.set_area_sd(200.0)
 mf.set_prec(1.0e-9)
-mf.set_eval_limits(cosmo, math.log(1e14), math.log(1e16), 0.0, 2.0)
+mf.set_eval_limits(cosmo, math.log(1e14), math.log(1e16), 0.0, 2.5)
 mf.prepare(cosmo)
 
 gf = psml.peek_gf()

benchmarks/data/numcosmo_cluster_counts.txt

@@ -1,4 +1,4 @@
 # counts  mmin  mmax  zmin  zmax
-        3578.9272152                1e+14                2e+14                    0                    2
-       959.715661575                2e+14                1e+15                    0                    2
-       7.19881988858                1e+15                1e+16                    0                    2
+       3612.50288061                1e+14                2e+14                    0                    2
+       962.829291777                2e+14                1e+15                    0                    2
+       7.19914009689                1e+15                1e+16                    0                    2

benchmarks/test_halomod_numbercounts.py

@@ -57,6 +57,73 @@ def sel(m, a):
         assert np.isfinite(nc)
         assert not np.allclose(nc, 0)
 
-        tol = max(0.01, np.sqrt(bench[0]) / bench[0] / 10)
+        tol = max(0.013, np.sqrt(bench[0]) / bench[0] / 10)
+        print(nc, bench[0], nc/bench[0]-1, tol)
+        assert np.allclose(nc, bench[0], atol=0, rtol=tol)
+
+
+def test_hmcalculator_number_counts_numcosmo_highacc():
+    cosmo = ccl.Cosmology(
+        Omega_c=0.25,
+        Omega_b=0.05,
+        h=0.7,
+        n_s=0.96,
+        sigma8=0.8,
+        Omega_k=0.0,
+        Omega_g=0,
+        Neff=0.0,
+        w0=-1.0,
+        wa=0.0,
+        T_CMB=2.7245,
+        mu_0=0.0,
+        transfer_function='eisenstein_hu',
+        matter_power_spectrum='linear'
+    )
+    mdef = ccl.halos.MassDef(200, 'matter')
+    hmf = ccl.halos.MassFuncTinker08(cosmo, mdef,
+                                     mass_def_strict=False)
+    hbf = ccl.halos.HaloBiasTinker10(cosmo, mass_def=mdef,
+                                     mass_def_strict=False)
+
+    benches = np.loadtxt("./benchmarks/data/numcosmo_cluster_counts.txt")
+
+    for i in range(benches.shape[0]):
+        bench = benches[i, :]
+        hmc = ccl.halos.HMCalculator(
+            cosmo, hmf, hbf, mdef,
+            log10M_min=np.log10(bench[1]),
+            log10M_max=np.log10(bench[2]),
+            integration_method_M='spline',
+            nlog10M=4096,
+        )
+
+        a_2 = 1.0 / (1.0 + bench[4])
+        a_1 = 1.0 / (1.0 + bench[3])
+
+        def sel(m, a):
+            m = np.atleast_1d(m)
+            a = np.atleast_1d(a)
+            val = np.zeros_like(m.reshape(-1, 1) * a.reshape(1, -1))
+            msk_a = (a > a_2) & (a < a_1)
+            msk_m = (m > bench[1]) & (m < bench[2])
+            val[msk_m, :] += 1
+            val[:, msk_a] += 1
+            msk = val == 2
+            val[~msk] = 0
+            val[msk] = 1.0
+            return val
+
+        area = 200 * (np.pi / 180)**2
+
+        nc = hmc.number_counts(
+            cosmo, sel,
+            amin=a_2,
+            amax=a_1,
+            na=4096,
+        ) * area
+        assert np.isfinite(nc)
+        assert not np.allclose(nc, 0)
+
+        tol = 1e-3
         print(nc, bench[0], nc/bench[0]-1, tol)
         assert np.allclose(nc, bench[0], atol=0, rtol=tol)