BUG fix and add test for high-z halofit (#932)

* TST add test for high-z halofit

* figure out why it passes

* TST redo test as ratio at high-z

* TST try again to make it fail

* try larger range for roots

* see what is going on

* see if it fails

* ENH try again

* try again

* grr

* do not fail for a bit

* try this

* try this

* try this

* try this

* try this

* try this

* try this

* try this

* clean it up

* full ci now

* full ci now

* TST test a bigger range of values

* Update pyccl/tests/test_halofit_highz.py

* DOC added changelog entry

* TST add test at very low sigma9

* BUG expand range for root

.github/workflows/ci.yml

@@ -37,7 +37,7 @@ jobs:
           channel-priority: strict
           show-channel-urls: true
           miniforge-version: latest
-          miniforge-variant: Mambaforge          
+          miniforge-variant: Mambaforge
 
       - name: remove homebrew
         if: matrix.os == 'macos-10.15'

CHANGELOG.md

@@ -1,5 +1,8 @@
 # Unreleased
 
+## Python library
+- Fixed high-z halofit behavior (#932)
+
 # v2.3.0 Changes
 
 ## Python library

pyccl/tests/test_halofit_highz.py

@@ -0,0 +1,31 @@
+import numpy as np
+import pyccl as ccl
+import pytest
+
+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='halofit')
+
+COSMO_LOWS8 = ccl.Cosmology(
+    Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.1, n_s=0.96,
+    transfer_function='bbks', matter_power_spectrum='halofit')
+
+
+@pytest.mark.parametrize("cosmo", [COSMO, COSMO_LOWS8])
+def test_halofit_highz(cosmo):
+    vals = [(25, 75)] + list(zip(range(0, 98), range(1, 99)))
+    for zl, zh in vals:
+        al = 1.0/(1 + zl)
+        ah = 1.0/(1 + zh)
+
+        k = np.logspace(0, 2, 10)
+        pkratl = (
+            ccl.nonlin_matter_power(cosmo, k, al)
+            / ccl.linear_matter_power(cosmo, k, al)
+        )
+        pkrath = (
+            ccl.nonlin_matter_power(cosmo, k, ah)
+            / ccl.linear_matter_power(cosmo, k, ah)
+        )
+
+        assert np.all(pkratl >= pkrath), (zl, zh, pkratl, pkrath)

src/ccl_halofit.c

@@ -62,8 +62,8 @@ static ccl_cosmology *create_w0eff_cosmo(double w0eff, ccl_cosmology *cosmo, int
     cosmo->params.Neff, mnu, cosmo->params.N_nu_mass,
     w0eff, 0, cosmo->params.h, norm_pk,
     cosmo->params.n_s, cosmo->params.bcm_log10Mc, cosmo->params.bcm_etab,
-    cosmo->params.bcm_ks, cosmo->params.mu_0, cosmo->params.sigma_0, 
-    cosmo->params.c1_mg, cosmo->params.c2_mg, cosmo->params.lambda_mg, 
+    cosmo->params.bcm_ks, cosmo->params.mu_0, cosmo->params.sigma_0,
+    cosmo->params.c1_mg, cosmo->params.c2_mg, cosmo->params.lambda_mg,
     cosmo->params.nz_mgrowth,
     cosmo->params.z_mgrowth, cosmo->params.df_mgrowth, status);
 
@@ -236,7 +236,7 @@ static double rsigma_func(double rsigma, void *p) {
   int gsl_status;
 
   lnkmin = hfd->plin->lkmin;
-  lnkmax = hfd->plin->lkmax;
+  lnkmax = fmax(hfd->plin->lkmax, log(30/rsigma));
   hfd->r = rsigma;
   hfd->r2 = rsigma * rsigma;
   F.function = &gauss_norm_int_func;
@@ -257,22 +257,26 @@ static double rsigma_func(double rsigma, void *p) {
   return result - 1.0;
 }
 
+static double lnrsigma_func(double lnrsigma, void *p) {
+  return rsigma_func(exp(lnrsigma), p);
+}
+
 static double get_rsigma(double a, struct hf_int_data data) {
-  double rsigma, rlow = 1e-2, rhigh = 1e2;
+  double rsigma, rlow = log(1e-64), rhigh = log(1e16);
   double flow, fhigh;
   int itr, max_itr = 1000, gsl_status;
   const gsl_root_fsolver_type *T;
   gsl_root_fsolver *s;
   gsl_function F;
 
   data.a = a;
-  F.function = &rsigma_func;
+  F.function = &lnrsigma_func;
   F.params = &data;
 
   // we have to bound the root, otherwise return -1
   // we will fiil in any -1's in the calling routine
-  flow = rsigma_func(rlow, &data);
-  fhigh = rsigma_func(rhigh, &data);
+  flow = lnrsigma_func(rlow, &data);
+  fhigh = lnrsigma_func(rhigh, &data);
   if (flow * fhigh > 0) {
     return -1;
   }
@@ -310,6 +314,7 @@ static double get_rsigma(double a, struct hf_int_data data) {
       *(data.status) |= gsl_status;
     }
   }
+  rsigma = exp(rsigma);
 
   return rsigma;
 }
@@ -564,7 +569,7 @@ halofit_struct* ccl_halofit_struct_new(ccl_cosmology *cosmo,
 
     for (i=0; i<n_a; ++i) {
       vals[i] = get_rsigma(a_vec[i], data);
-      if (*status != 0) {
+      if ((*status != 0) || (vals[i] <= 0)) {
         *status = CCL_ERROR_ROOT;
         ccl_cosmology_set_status_message(
           cosmo,
@@ -573,51 +578,6 @@ halofit_struct* ccl_halofit_struct_new(ccl_cosmology *cosmo,
         break;
       }
     }
-
-    // now go backwards and fill any -1's
-    // one must work, so set an error if not
-    if (vals[n_a-1] == -1) {
-      *status = CCL_ERROR_ROOT;
-      ccl_cosmology_set_status_message(
-        cosmo,
-        "ccl_halofit.c: ccl_halofit_struct_new(): "
-        "could not solve for non-linear scale for halofit at least once\n");
-    }
-    if (*status == 0) {
-      // linearly set rsigma to a very small number at high-z
-      double min_a = -1, max_a = -1;
-      double max_val = -1;
-      double w;
-
-      // first non -1 value and scale factor of last -1 value
-      for (i=1; i<n_a; ++i) {
-        if (vals[i] != -1) {
-          max_a = a_vec[i];
-          max_val = vals[i];
-          break;
-        }
-      }
-      // scale factor of first -1 value
-      for (i=0; i<n_a; ++i) {
-        if (vals[i] == -1) {
-          min_a = a_vec[i];
-          break;
-        }
-      }
-
-      if (min_a != -1) {
-        // at least one value is -1 so set the zeroth value
-        vals[0] = 1e-6;
-
-        // interp any values that remain -1
-        for(i=1; i<n_a-1; ++i) {
-          if (vals[i] == -1) {
-            w = (a_vec[i] - min_a) / (max_a - min_a);
-            vals[i] = w * max_val + (1.0 - w) * vals[0];
-          }
-        }
-      }
-    }
   }
 
   // spline the non-linear scales
@@ -697,7 +657,8 @@ halofit_struct* ccl_halofit_struct_new(ccl_cosmology *cosmo,
 
       gsl_status = gsl_integration_cquad(
         &F,
-        lnkmin, lnkmax,
+        lnkmin,
+        fmax(lnkmax, log(30/rsigma)),
         0.0, cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL,
         workspace, &result, NULL, NULL);
 
@@ -761,7 +722,8 @@ halofit_struct* ccl_halofit_struct_new(ccl_cosmology *cosmo,
 
       gsl_status = gsl_integration_cquad(
         &F,
-        lnkmin, lnkmax,
+        lnkmin,
+        fmax(lnkmax, log(30/rsigma)),
         0.0, cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL,
         workspace, &result, NULL, NULL);