Bug in magnification bias (#963)

* fixed bug in use of magnification bias

* added test

* fixed pure lensing case

* additional test

pyccl/tests/test_tracers.py

@@ -41,6 +41,38 @@ def get_tracer(tracer_type, cosmo=None, **tracer_kwargs):
     return tr, ntr
 
 
+def test_tracer_mag_0p4():
+    z = np.linspace(0., 1., 2000)
+    n = dndz(z)
+    # Small bias so the magnification contribution
+    # is significant (if there at all)
+    b = np.ones_like(z)*0.1
+    s_no = np.ones_like(z)*0.4
+    s_yes = np.zeros_like(z)
+    # Tracer with no magnification by construction
+    t1 = ccl.NumberCountsTracer(COSMO, True,
+                                dndz=(z, n),
+                                bias=(z, b))
+    # Tracer with s=0.4
+    t2 = ccl.NumberCountsTracer(COSMO, True,
+                                dndz=(z, n),
+                                bias=(z, b),
+                                mag_bias=(z, s_no))
+    # Tracer with magnification
+    t3 = ccl.NumberCountsTracer(COSMO, True,
+                                dndz=(z, n),
+                                bias=(z, b),
+                                mag_bias=(z, s_yes))
+    ls = np.array([2, 200, 2000])
+    cl1 = ccl.angular_cl(COSMO, t1, t1, ls)
+    cl2 = ccl.angular_cl(COSMO, t2, t2, ls)
+    cl3 = ccl.angular_cl(COSMO, t3, t3, ls)
+    # Check cl1 == cl2
+    assert np.all(np.fabs(cl2/cl1-1) < 1E-5)
+    # Check cl1 != cl3
+    assert np.all(np.fabs(cl3/cl1-1) > 1E-2)
+
+
 @pytest.mark.parametrize('tracer_type', ['nc', 'wl'])
 def test_tracer_dndz_smoke(tracer_type):
     tr, _ = get_tracer(tracer_type)
@@ -220,6 +252,52 @@ def test_tracer_lensing_kernel_spline_vs_gsl_intergation(z_min, z_max,
         assert np.allclose(w_wl_spline[0], w_wl_gsl[0], atol=5e-9, rtol=1e-5)
 
 
+@pytest.mark.parametrize('z_min, z_max, n_z_samples', [(0.0, 1.0, 2000),
+                                                       (0.0, 1.0, 1000),
+                                                       (0.0, 1.0, 500),
+                                                       (0.0, 1.0, 100),
+                                                       (0.3, 1.0, 1000)])
+def test_tracer_magnification_kernel_spline_vs_gsl_intergation(z_min, z_max,
+                                                               n_z_samples):
+    # Create a new Cosmology object so that we're not messing with the other
+    # tests
+    cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67,
+                          sigma8=0.8, n_s=0.96,
+                          transfer_function='bbks',
+                          matter_power_spectrum='linear')
+    z = np.linspace(z_min, z_max, n_z_samples)
+    n = dndz(z)
+    b = np.zeros_like(z)
+
+    # Make sure case where z[0] > 0 and n[0] > 0 is tested for
+    if z_min > 0:
+        assert n[0] > 0
+
+    cosmo.cosmo.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION = True
+    tr_mg = ccl.NumberCountsTracer(cosmo, False, dndz=(z, n),
+                                   bias=(z, b), mag_bias=(z, b))
+    w_mg_spline, _ = tr_mg.get_kernel(chi=None)
+
+    cosmo.cosmo.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION = False
+    tr_mg = ccl.NumberCountsTracer(cosmo, False, dndz=(z, n),
+                                   bias=(z, b), mag_bias=(z, b))
+    w_mg_gsl, chi = tr_mg.get_kernel(chi=None)
+    tr_wl = ccl.WeakLensingTracer(cosmo, dndz=(z, n))
+    w_wl_gsl, _ = tr_wl.get_kernel(chi=None)
+
+    # Peak of kernel is ~1e-5
+    if n_z_samples >= 1000:
+        assert np.allclose(w_mg_spline[1], w_mg_gsl[1],
+                           atol=2e-10, rtol=1e-9)
+        assert np.allclose(w_mg_spline[1], -2*w_wl_gsl[0],
+                           atol=2e-10, rtol=1e-9)
+    else:
+        assert np.allclose(w_mg_spline[1], w_mg_gsl[1],
+                           atol=1e-8, rtol=1e-5)
+        assert np.allclose(w_mg_spline[1], -2*w_wl_gsl[0],
+                           atol=1e-8, rtol=1e-5)
+
+
 def test_tracer_delta_function_nz():
     z = np.linspace(0., 1., 2000)
     z_s_idx = int(z.size*0.8)

src/ccl_tracers.c

@@ -350,31 +350,31 @@ static void integrate_lensing_kernel_spline(ccl_cosmology *cosmo,
                                               int nchi, double* chi_arr, double* wL_arr,
                                               int* status) {
   double* chi_of_z_array = malloc(nz*sizeof(double));
-  double* sz_array = malloc(nz*sizeof(double));
-  if(chi_of_z_array == NULL || sz_array == NULL) {
+  double* qz_array = malloc(nz*sizeof(double));
+  if(chi_of_z_array == NULL || qz_array == NULL) {
     *status = CCL_ERROR_MEMORY;
     ccl_cosmology_set_status_message(
       cosmo,
       "ccl_tracers.c: lensing_kernel_integrate_spline(): error allocating memory\n");
   }
 
   if(*status == 0) {
-    // Fill chi and sz arrays
+    // Fill chi and qz arrays
     for(int i=0; i<nz; i++) {
       double a = 1./(1+z_arr[i]);
       chi_of_z_array[i] = ccl_comoving_radial_distance(cosmo, a, status);
       if(sz_f != NULL) {
-        sz_array[i] = ccl_f1d_t_eval(sz_f, z_arr[i]);
+        qz_array[i] = 1-2.5*ccl_f1d_t_eval(sz_f, z_arr[i]);
       } else {
-        sz_array[i] = 1.0;
+        qz_array[i] = 1.0;
       }
     }
   }
 
   if(*status == 0) {
     #pragma omp parallel default(none) \
-                        shared(cosmo, nz, z_arr, nz_arr, nz_norm, sz_f, \
-                               chi_of_z_array, sz_array, \
+                        shared(cosmo, nz, z_arr, nz_arr, nz_norm, \
+                               chi_of_z_array, qz_array, \
                                nchi, chi_arr, wL_arr, status, gsl_interp_akima)
     {
       int local_status = *status;
@@ -403,7 +403,11 @@ static void integrate_lensing_kernel_spline(ccl_cosmology *cosmo,
               integrand_array[i] = 0.0;
               i_chi_end = i+1;
             } else {
-              integrand_array[i] = lensing_kernel_integrand(cosmo, chi_of_z_array[i], chi_end, nz_arr[i], sz_array[i], &local_status);
+	      integrand_array[i] = lensing_kernel_integrand(cosmo,
+							    chi_of_z_array[i],
+							    chi_end, nz_arr[i],
+							    qz_array[i],
+							    &local_status);
             }
           }
           if(local_status) {
@@ -446,7 +450,7 @@ static void integrate_lensing_kernel_spline(ccl_cosmology *cosmo,
   }
 
   free(chi_of_z_array);
-  free(sz_array);
+  free(qz_array);
 }
 
 //Returns number of divisions on which