SSC covariance using halo model response (#855)

* SSC, tk3d Co-authored-by: Tilman Troester <tilman.troester@gmail.com>
LSSTDESC · Jul 9, 2021 · e287373 · e287373
1 parent 14cb88a
commit e287373
Show file tree

Hide file tree

Showing 20 changed files with 1,609 additions and 33 deletions.
diff --git a/benchmarks/ccl_test_f2d.c b/benchmarks/ccl_test_f2d.c
@@ -171,7 +171,7 @@ CTEST2(f2d,a_overflow) {
 CTEST2(f2d,sanity) {
   int status=0;
   ccl_f2d_t *psp;
-  double fka;
+  double fka,dfka;
   double lktest=-2.,atest=0.5;
   double alo=0.02;
 
@@ -195,23 +195,32 @@ CTEST2(f2d,sanity) {
   fka=ccl_f2d_t_eval(psp,lktest,atest,NULL,&status);
   ASSERT_TRUE(status==0);
   ASSERT_DBL_NEAR(1,fka/fka_model_analytical(exp(lktest),atest));
+  dfka=ccl_f2d_t_dlogf_dlk_eval(psp, lktest, atest,NULL,&status);
+  ASSERT_TRUE(status==0);
+  ASSERT_DBL_NEAR(-1.,dfka);
 
   //Evaluate at very high z and see if it checks out
   fka=ccl_f2d_t_eval(psp,lktest,alo,NULL,&status);
   ASSERT_TRUE(status==0);
   ASSERT_DBL_NEAR(1,fka/fka_model_analytical(exp(lktest),alo));
+  dfka = ccl_f2d_t_dlogf_dlk_eval(psp, lktest, alo,NULL,&status);
+  ASSERT_DBL_NEAR(-1.,dfka);
 
   //Evaluate at very high k and see if it checks out
   double lkhi=data->lk_arr[data->n_k-1]*1.1;
   fka=ccl_f2d_t_eval(psp,lkhi,atest,NULL,&status);
   ASSERT_TRUE(status==0);
   ASSERT_DBL_NEAR(1,fka/fka_model_analytical(exp(lkhi),atest));
+  dfka = ccl_f2d_t_dlogf_dlk_eval(psp, lkhi, atest,NULL,&status);
+  ASSERT_DBL_NEAR(-1.,dfka);
 
   //Evaluate at very low k and see if it checks out
   double lklo=data->lk_arr[0]/1.1;
   fka=ccl_f2d_t_eval(psp,lklo,atest,NULL,&status);
   ASSERT_TRUE(status==0);
   ASSERT_DBL_NEAR(1,fka/fka_model_analytical(exp(lklo),atest));
+  dfka = ccl_f2d_t_dlogf_dlk_eval(psp, lklo, atest,NULL,&status);
+  ASSERT_DBL_NEAR(-1.,dfka);
 
   ccl_f2d_t_free(psp);
 }

diff --git a/benchmarks/data/covariances/ssc_WL_cov_matrix.txt b/benchmarks/data/covariances/ssc_WL_cov_matrix.txt
@@ -0,0 +1,10 @@
+3.5827589908e-18	2.1933032306e-18	1.0195609562e-18	3.9498025255e-19	1.5006849349e-19	5.8863418831e-20	1.9539477848e-20	4.6939126814e-21	8.3192018555e-22	1.1336282393e-22	
+2.1933032306e-18	1.4768783335e-18	7.3887186615e-19	2.9525366302e-19	1.0947282832e-19	4.3019749519e-20	1.5489453850e-20	4.1033885830e-21	7.7985249459e-22	1.1248585689e-22	
+1.0195609562e-18	7.3887186615e-19	3.9390130562e-19	1.6323203150e-19	5.9907638875e-20	2.3173464828e-20	8.7427311235e-21	2.5079467781e-21	5.0782251995e-22	7.6922542286e-23	
+3.9498025255e-19	2.9525366302e-19	1.6323203150e-19	6.9629693452e-20	2.5603803143e-20	9.6882447681e-21	3.7135387390e-21	1.1204304046e-21	2.3665509805e-22	3.6995745170e-23	
+1.5006849349e-19	1.0947282832e-19	5.9907638875e-20	2.5603803143e-20	9.5379183119e-21	3.5923698595e-21	1.3497466983e-21	4.0830265040e-22	8.7423838992e-23	1.3831025304e-23	
+5.8863418831e-20	4.3019749519e-20	2.3173464828e-20	9.6882447681e-21	3.5923698595e-21	1.3879675434e-21	5.2023575248e-22	1.5198271589e-22	3.1713716352e-23	4.9275462171e-24	
+1.9539477848e-20	1.5489453850e-20	8.7427311235e-21	3.7135387390e-21	1.3497466983e-21	5.2023575248e-22	2.0461796895e-22	6.1623620544e-23	1.2919437394e-23	2.0071723283e-24	
+4.6939126814e-21	4.1033885830e-21	2.5079467781e-21	1.1204304046e-21	4.0830265040e-22	1.5198271589e-22	6.1623620544e-23	2.0168158521e-23	4.4820186218e-24	7.2448013820e-25	
+8.3192018555e-22	7.7985249459e-22	5.0782251995e-22	2.3665509805e-22	8.7423838992e-23	3.1713716352e-23	1.2919437394e-23	4.4820186218e-24	1.0537038760e-24	1.7683998813e-25	
+1.1336282393e-22	1.1248585689e-22	7.6922542286e-23	3.6995745170e-23	1.3831025304e-23	4.9275462171e-24	2.0071723283e-24	7.2448013820e-25	1.7683998813e-25	3.0411198865e-26	
diff --git a/benchmarks/data/covariances/ssc_WL_ell.txt b/benchmarks/data/covariances/ssc_WL_ell.txt
@@ -0,0 +1,10 @@
+2.984599475985548089e+00
+6.646579196058059580e+00
+1.480165608984571257e+01
+3.296267396196881094e+01
+7.340650722647497162e+01
+1.634732458115384190e+02
+3.640481355925554681e+02
+8.107200928842213443e+02
+1.805440008465850724e+03
+4.020639987560632562e+03
diff --git a/benchmarks/data/covariances/ssc_WL_nofz.txt b/benchmarks/data/covariances/ssc_WL_nofz.txt
diff --git a/benchmarks/test_covariances.py b/benchmarks/test_covariances.py
@@ -0,0 +1,73 @@
+import os
+
+import numpy as np
+
+import pyccl as ccl
+
+
+def test_ssc_WL():
+    # Compare against Benjamin Joachimi's code. An overview of the methodology
+    # is given in appendix E.2 of 2007.01844.
+    data_dir = os.path.join(os.path.dirname(__file__), "data/covariances/")
+
+    h = 0.7
+    cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=h, n_s=0.97,
+                          sigma8=0.8, m_nu=0.0)
+
+    mass_def = ccl.halos.MassDef200m()
+    hmf = ccl.halos.MassFuncTinker10(cosmo,
+                                     mass_def=mass_def)
+    hbf = ccl.halos.HaloBiasTinker10(cosmo,
+                                     mass_def=mass_def)
+    nfw = ccl.halos.HaloProfileNFW(ccl.halos.ConcentrationDuffy08(mass_def),
+                                   fourier_analytic=True)
+    hmc = ccl.halos.HMCalculator(cosmo, hmf, hbf, mass_def)
+
+    n_z = 100
+
+    n_k = 200
+    k_min = 1e-4
+    k_max = 1e2
+
+    a = np.linspace(1/(1+6), 1, n_z)
+    k = np.geomspace(k_min, k_max, n_k)
+
+    tk3D = ccl.halos.halomod_Tk3D_SSC(cosmo=cosmo, hmc=hmc,
+                                      prof1=nfw,
+                                      prof2=nfw,
+                                      prof12_2pt=None,
+                                      normprof1=True, normprof2=True,
+                                      lk_arr=np.log(k), a_arr=a,
+                                      use_log=True)
+
+    z, nofz = np.loadtxt(os.path.join(data_dir, "ssc_WL_nofz.txt"),
+                         unpack=True)
+    WL_tracer = ccl.WeakLensingTracer(cosmo, (z, nofz))
+
+    ell = np.loadtxt(os.path.join(data_dir, "ssc_WL_ell.txt"))
+
+    fsky = 0.05
+
+    sigma2_B = ccl.sigma2_B_disc(cosmo, a=a, fsky=fsky)
+    cov_ssc = ccl.covariances.angular_cl_cov_SSC(cosmo,
+                                                 cltracer1=WL_tracer,
+                                                 cltracer2=WL_tracer,
+                                                 ell=ell, tkka=tk3D,
+                                                 sigma2_B=(a, sigma2_B),
+                                                 fsky=None)
+    var_ssc_ccl = np.diag(cov_ssc)
+    off_diag_1_ccl = np.diag(cov_ssc, k=1)
+
+    cov_ssc_bj = np.loadtxt(os.path.join(data_dir, "ssc_WL_cov_matrix.txt"))
+
+    # At large scales, CCL uses a different convention for the Limber
+    # approximation. This factor accounts for this difference
+    ccl_limber_shear_fac = np.sqrt((ell-1)*ell*(ell+1)*(ell+2))/(ell+1/2)**2
+    cov_ssc_bj_corrected = cov_ssc_bj * np.outer(ccl_limber_shear_fac**2,
+                                                 ccl_limber_shear_fac**2)
+    var_bj = np.diag(cov_ssc_bj_corrected)
+    off_diag_1_bj = np.diag(cov_ssc_bj_corrected, k=1)
+
+    assert np.all(np.fabs(var_ssc_ccl/var_bj - 1) < 3e-2)
+    assert np.all(np.fabs(off_diag_1_ccl/off_diag_1_bj - 1) < 3e-2)
+    assert np.all(np.fabs(cov_ssc/cov_ssc_bj_corrected - 1) < 3e-2)
diff --git a/include/ccl_f2d.h b/include/ccl_f2d.h
@@ -88,6 +88,17 @@ ccl_f2d_t *ccl_f2d_t_new(int na,double *a_arr,
 double ccl_f2d_t_eval(ccl_f2d_t *fka,double lk,double a,void *cosmo,
 		      int *status);
 
+/**
+ * Evaluate logarithmic derivative of 2D function of k and a defined by ccl_f2d_t structure wrt k.
+ * @param fka ccl_f2d_t structure defining f(k,a).
+ * @param lk Natural logarithm of the wavenumber.
+ * @param a Scale factor.
+ * @param cosmo ccl_cosmology structure, only needed if evaluating f(k,a) at small scale factors outside the interpolation range, and if fka was initialized with extrap_linear_growth = ccl_f2d_cclgrowth.
+ * @param status Status flag. 0 if there are no errors, nonzero otherwise.
+ */
+double ccl_f2d_t_dlogf_dlk_eval(ccl_f2d_t *f2d,double lk,double a,void *cosmo, int *status);
+
+
 /**
  * F2D structure destructor.
  * Frees up all memory associated with a f2d structure.

diff --git a/include/ccl_power.h b/include/ccl_power.h
@@ -17,6 +17,35 @@ void ccl_rescale_linpower(ccl_cosmology* cosmo, ccl_f2d_t *psp,
                           int rescale_mg, int rescale_norm,
                           int *status);
 
+/**
+ * Variance of the projected matter density field with 2D (top-hat) smoothing scale R [Mpc].
+ * Returns sigma2(R) for specified cosmology at a = 1.
+ * @param cosmo Cosmology parameters and configurations
+ * @param R Smoothing scale, in [Mpc] units
+ * @param a scale factor
+ * @param psp input power spectrum.
+ * @param status Status flag. 0 if there are no errors, nonzero otherwise.
+ * For specific cases see documentation for ccl_error.c
+ * @return sigma(R).
+ */
+double ccl_sigma2B(ccl_cosmology *cosmo,double R,double a,
+                   ccl_f2d_t *psp, int *status);
+
+/**
+ * As `ccl_sigma2B`, calculated for an array of scale factors and smoothing scales.
+ * @param cosmo Cosmology parameters and configurations
+ * @param na number of scale factor values
+ * @param a scale factor values
+ * @param R Smoothing scale values, in [Mpc] units
+ * @param sigma2B_out output values of the variance calculated for the input a and R.
+ * @param psp input power spectrum.
+ * @param status Status flag. 0 if there are no errors, nonzero otherwise.
+ * For specific cases see documentation for ccl_error.c
+ * @return sigma(R,a).
+ */
+void ccl_sigma2Bs(ccl_cosmology *cosmo,int na, double *a, double *R,
+                  double *sigma2B_out, ccl_f2d_t *psp, int *status);
+
 /**
  * Variance of the matter density field with (top-hat) smoothing scale R [Mpc].
  * Returns sigma(R) for specified cosmology at a = 1.
@@ -26,7 +55,6 @@ void ccl_rescale_linpower(ccl_cosmology* cosmo, ccl_f2d_t *psp,
  * @param psp input power spectrum.
  * @param status Status flag. 0 if there are no errors, nonzero otherwise.
  * For specific cases see documentation for ccl_error.c
- * @return sigma(R).
  */
 double ccl_sigmaR(ccl_cosmology *cosmo, double R, double a,
                   ccl_f2d_t *psp, int * status);

diff --git a/pyccl/__init__.py b/pyccl/__init__.py
@@ -47,7 +47,7 @@
 from .tracers import Tracer, NumberCountsTracer, WeakLensingTracer, CMBLensingTracer, \
     tSZTracer, get_density_kernel, get_kappa_kernel, get_lensing_kernel
 from .cls import angular_cl
-from .covariances import angular_cl_cov_cNG
+from .covariances import angular_cl_cov_cNG, angular_cl_cov_SSC, sigma2_B_disc, sigma2_B_from_mask
 
 # Useful constants and unit conversions
 physical_constants = lib.cvar.constants

diff --git a/pyccl/ccl_covs.i b/pyccl/ccl_covs.i
@@ -9,9 +9,31 @@
 // Enable vectorised arguments for arrays
 %apply (double* IN_ARRAY1, int DIM1) {(double* ell1, int nell1)};
 %apply (double* IN_ARRAY1, int DIM1) {(double* ell2, int nell2)};
+%apply (double* IN_ARRAY1, int DIM1) {(double* s2b, int ns2b)};
+%apply (double* IN_ARRAY1, int DIM1) {(double* a, int na)};
+%apply (double* IN_ARRAY1, int DIM1) {(double* R, int nR)};
 %apply (int DIM1, double* ARGOUT_ARRAY1) {(int nout, double* output)};
 
-%feature("pythonprepend") angular_cl_vec %{
+%feature("pythonprepend") sigma2b_vec %{
+    if len(R) != nout:
+        raise CCLError("Input shape for `R` must match `(nout,)`!")
+%}
+
+%inline %{
+
+void sigma2b_vec(ccl_cosmology * cosmo,
+                 double *a, int na,
+                 double *R, int nR,
+                 ccl_f2d_t *psp,
+                 int nout, double* output,
+                 int *status)
+{
+  ccl_sigma2Bs(cosmo, na, a, R, output, psp, status);
+}
+
+%}
+
+%feature("pythonprepend") angular_cov_vec %{
     if len(ell1)*len(ell2) != nout:
         raise CCLError("Input shape for `ell1` and `ell2` must match `(nout,)`!")
 %}
@@ -38,3 +60,37 @@ void angular_cov_vec(ccl_cosmology * cosmo,
 }
 
 %}
+
+%feature("pythonprepend") angular_cov_ssc_vec %{
+    if len(ell1)*len(ell2) != nout:
+        raise CCLError("Input shape for `ell1` and `ell2` must match `(nout,)`!")
+%}
+
+%inline %{
+
+void angular_cov_ssc_vec(ccl_cosmology * cosmo,
+                         ccl_cl_tracer_collection_t *clt1,
+                         ccl_cl_tracer_collection_t *clt2,
+                         ccl_cl_tracer_collection_t *clt3,
+                         ccl_cl_tracer_collection_t *clt4,
+                         ccl_f3d_t *tspec,
+                         double *a, int na, 
+                         double *s2b, int ns2b, 
+                         double* ell1, int nell1,
+                         double* ell2, int nell2,
+                         int integration_type,
+                         int chi_exponent, double prefac,
+                         int nout, double* output,
+                         int *status)
+{
+  ccl_f1d_t *s2b_f=ccl_f1d_t_new(na, a, s2b,
+                                 s2b[0], s2b[na-1],
+                                 0, 0, status);
+  ccl_angular_cl_covariance(cosmo, clt1, clt2, clt3, clt4, tspec,
+                            nell1,  ell1, nell2, ell2, output,
+                            integration_type, chi_exponent, s2b_f,
+                            prefac, status);
+  ccl_f1d_t_free(s2b_f);
+}
+
+%}
diff --git a/pyccl/ccl_pk2d.i b/pyccl/ccl_pk2d.i
@@ -49,4 +49,17 @@ void pk2d_eval_multi(ccl_f2d_t *psp,double* lkarr,int nk,
   for(int ii=0;ii<ndout;ii++)
     doutput[ii]=ccl_f2d_t_eval(psp,lkarr[ii],a,cosmo,status);
 }
+
+double pk2d_der_eval_single(ccl_f2d_t *psp,double lk,double a,ccl_cosmology *cosmo,int *status)
+{
+  return ccl_f2d_t_dlogf_dlk_eval(psp,lk,a,cosmo,status);
+}
+
+ void pk2d_der_eval_multi(ccl_f2d_t *psp,double* lkarr,int nk,
+                          double a,ccl_cosmology *cosmo,
+                          int ndout,double *doutput,int *status)
+{
+  for(int ii=0;ii<ndout;ii++)
+    doutput[ii]=ccl_f2d_t_dlogf_dlk_eval(psp,lkarr[ii],a,cosmo,status);
+}
 %}