[Draft] Increasing the test coverage (#40)

* Added hmcode for testing

* Test also sigma8

* Added valgrind mode to the nightly, where it is less expensive to test

* Increased level of testing for the actual PR

* Added mode which is very time-saving for rare/comparatively unimportant features: comparing only withe full output already activated

* Updating test script to not give false error messages

* Making sure to catch the HMcode case for P_k_ini not analytic, which is a very special failure case (!)

* Corrected Pk_ini_type notation

* Corrected full settings + check shootings

* Adding theta_s to naming

* Adding S8 and tau_reio

---------

Co-authored-by: schoeneberg <[email protected]>

.github/workflows/test_nightly.yml

@@ -93,7 +93,7 @@ jobs:
       run: |
         source venv/bin/activate
         cd main_class/python
-        TEST_LEVEL=3 CLASS_VERBOSE=1 python "$(which nosetests)" -a dump_ini_files test_class.py
+        TEST_LEVEL=3 VALGRIND_MODE=1 CLASS_VERBOSE=1 python "$(which nosetests)" -a dump_ini_files test_class.py
         deactivate
     - name: Valgrind 🤖
       run: |

.github/workflows/test_on_pull_request.yml

@@ -38,7 +38,7 @@ jobs:
       run: |
         source venv/bin/activate
         cd main_class/python
-        OMP_NUM_THREADS=16 COMPARE_OUTPUT_REF=1 TEST_LEVEL=2 nosetests -v -a test_scenario test_class.py --nologcapture --nocapture
+        OMP_NUM_THREADS=16 COMPARE_OUTPUT_REF=1 TEST_LEVEL=3 nosetests -v -a test_scenario test_class.py --nologcapture --nocapture
         deactivate
     - name: Upload plots 📤
       if: success() || failure()

include/fourier.h

@@ -356,6 +356,7 @@ extern "C" {
 
   int fourier_get_k_list(
                          struct precision *ppr,
+                         struct primordial *ppm,
                          struct perturbations * ppt,
                          struct fourier * pfo
                          );

python/test_class.py

@@ -93,6 +93,8 @@
 COMPARE_OUTPUT_GAUGE = bool(int(os.getenv('COMPARE_OUTPUT_GAUGE', '0'))) # Compare synchronous and Newtonian gauge outputs?
 COMPARE_OUTPUT_REF = bool(int(os.getenv('COMPARE_OUTPUT_REF', '0'))) # Compare classy with classyref?
 POWER_ALL = bool(int(os.getenv('POWER_ALL', '0'))) # Combine every extension with each other? (Very slow!)
+VALGRIND_MODE = bool(int(os.getenv('VALGRIND_MODE', '0')))
+
 TEST_LEVEL = int(os.getenv('TEST_LEVEL', '0')) # 0 <= TEST_LEVEL <= 3
 
 if COMPARE_OUTPUT_REF:
@@ -108,6 +110,7 @@
 COMPARE_PK_RELATIVE_ERROR = 1e-2
 COMPARE_PK_RELATIVE_ERROR_GAUGE = 5*1e-2
 COMPARE_PK_ABSOLUTE_ERROR = 1e-20
+COMPARE_SHOOTING_RELATIVE_ERROR = 1e-5
 
 # Dictionary of models to test the wrapper against. Each of these scenario will
 # be run against all the possible output choices (nothing, tCl, mPk, etc...),
@@ -116,6 +119,10 @@
 # against. Indeed, when not specifying a field, CLASS takes the default input.
 CLASS_INPUT = {}
 
+# power = compare ALL settings with each of the entries
+# normal = compare only the power settings with this entry
+# onlyfull = compare only for 'full' output settings (see below)
+
 CLASS_INPUT['Output_spectra'] = (
     [{'output': 'mPk', 'P_k_max_1/Mpc': 2},
      {'output': 'tCl'},
@@ -126,7 +133,8 @@
     'power')
 
 CLASS_INPUT['Nonlinear'] = (
-    [{'non linear': 'halofit'}],
+    [{'non_linear': 'halofit'},
+    {'non_linear': 'hmcode'}],
     'power')
 
 CLASS_INPUT['Lensing'] = (
@@ -148,44 +156,65 @@
     CLASS_INPUT['modes'] = (
         [{'modes': 't'},
         {'modes': 's, t'}],
-        'power')
+        'normal')
 
     CLASS_INPUT['Tensor_method'] = (
         [{'tensor method': 'exact'},
         {'tensor method': 'photons'}],
-        'power')
+        'onlyfull')
 
 if TEST_LEVEL > 2:
     CLASS_INPUT['Isocurvature_modes'] = (
         [{'ic': 'ad,nid,cdi', 'c_ad_cdi': -0.5}],
-        'normal')
+        'onlyfull')
 
     CLASS_INPUT['Scalar_field'] = (
         [{'Omega_scf': 0.1, 'attractor_ic_scf': 'yes',
         'scf_parameters': '10, 0, 0, 0'}],
-        'normal')
+        'onlyfull')
 
     CLASS_INPUT['Inflation'] = (
-        [{'P_k_ini type': 'inflation_V'},
-        {'P_k_ini type': 'inflation_H'}],
-        'normal')
-    # CLASS_INPUT['Inflation'] = (
-    #     [{'P_k_ini type': 'inflation_V'},
-    #     {'P_k_ini type': 'inflation_H'},
-    #     {'P_k_ini type': 'inflation_V_end'}],
-    #     'normal')
+        [{'Pk_ini_type': 'inflation_V'},
+        {'Pk_ini_type': 'inflation_H'}],
+        'onlyfull')
+
+    CLASS_INPUT['Shooting'] = (
+        [{'sigma8': 0.8},
+        {'S8': 0.8},
+        {'theta_s_100': 1.04},
+        {'Neff': 4.0},
+        {'Neff': 2.0},
+        {'Neff': 2.0,'N_ncdm':1,'m_ncdm':0.06},
+        {'Neff': 3.0,'N_ncdm':1,'m_ncdm':0.06},
+        {'Neff': 4.0,'N_ncdm':1,'m_ncdm':0.06},
+        {'Neff': 4.0,'N_ncdm':1,'m_ncdm':0.06, 'theta_s_100':1.04,'sigma8':0.8},
+        {'Neff': 4.0,'N_ncdm':1,'m_ncdm':0.06, 'theta_s_100':1.04,'S8':0.8},
+        {'tau_reio':0.05}
+        ],
+        'onlyfull')
+# TODO :: possibly another one with nonlinear=hmcode in the future, but for now this is more stable
+FULL_SETTINGS = {'output':'mPk mTk vTk tCl pCl lCl sCl nCl','P_k_max_1/Mpc': 2,'lensing':'yes','modes':'s,t'}
 
 if POWER_ALL:
     for k, v in iteritems(CLASS_INPUT):
         models, state = v
         CLASS_INPUT[k] = (models, 'power')
+if VALGRIND_MODE:
+    for k, v in iteritems(CLASS_INPUT):
+        if k == 'Output_spectra':
+            continue
+        models, state = v
+        CLASS_INPUT[k] = (models, 'normal')
 
 INPUTPOWER = []
 INPUTNORMAL = [{}]
+INPUTONLYFULL = [{}]
 for key, value in list(CLASS_INPUT.items()):
     models, state = value
     if state == 'power':
         INPUTPOWER.append([{}]+models)
+    elif state == 'onlyfull':
+        INPUTONLYFULL.extend(models)
     else:
         INPUTNORMAL.extend(models)
 
@@ -199,6 +228,10 @@
                 temp_dict.update(elem)
             DICTARRAY.append(temp_dict)
 
+    for onlyfull in INPUTONLYFULL:
+        temp_dict = onlyfull.copy()
+        temp_dict.update(FULL_SETTINGS)
+        DICTARRAY.append(temp_dict)
 
 TUPLE_ARRAY = []
 for e in DICTARRAY:
@@ -212,7 +245,7 @@ def powerset(iterable):
         itertools.combinations(xs, n) for n in range(1, len(xs)+1))
 
 def custom_name_func(testcase_func, param_num, param):
-    special_keys = ['N_ncdm']
+    special_keys = ['N_ncdm','Pk_ini_type','ic','Neff','sigma8','S8','theta_s_100','tau_reio']
     somekeys = []
     for key in param.args[0]:
         if key in special_keys:
@@ -435,9 +468,9 @@ def has_incompatible_input(self):
             if not has_tensor(self.scenario):
                 should_fail = True
 
-        # If we have specified non linear, we must have some form of
+        # If we have specified non_linear, we must have some form of
         # perturbations output.
-        if 'non linear' in self.scenario:
+        if 'non_linear' in self.scenario:
             if 'output' not in self.scenario:
                 should_fail = True
 
@@ -460,7 +493,7 @@ def has_incompatible_input(self):
 
         # If we use inflation module, we must have scalar modes,
         # tensor modes, no vector modes and we should only have adiabatic IC:
-        if 'P_k_ini type' in self.scenario and self.scenario['P_k_ini type'].find('inflation') != -1:
+        if 'Pk_ini_type' in self.scenario and self.scenario['Pk_ini_type'].find('inflation') != -1:
             if 'modes' not in self.scenario:
                 should_fail = True
             else:
@@ -472,6 +505,8 @@ def has_incompatible_input(self):
                     should_fail = True
             if 'ic' in self.scenario and self.scenario['ic'].find('i') != -1:
                 should_fail = True
+            if 'non_linear' in self.scenario and self.scenario['non_linear'].find('hmcode') != -1:
+                should_fail = True
 
 
         return should_fail
@@ -545,6 +580,36 @@ def compare_output(self, reference, reference_name, candidate, candidate_name, r
                     self.pk_faulty_plot(k, reference_pk, reference_name, candidate_pk, candidate_name, rtol_pk)
                     status_pass = False
 
+        if 'output' in self.scenario and self.scenario['output'].find('mPk') != -1:
+            ref_sigma8 = reference.sigma8()
+            can_sigma8 = candidate.sigma8()
+            ref_sigma = reference.sigma(8, 0, h_units=True)
+            can_sigma = candidate.sigma(8, 0, h_units=True)
+            try:
+                np.testing.assert_allclose([ref_sigma8, ref_sigma, can_sigma8],
+                                           [can_sigma8, can_sigma, can_sigma],
+                                           rtol=rtol_pk,
+                                           atol=COMPARE_PK_ABSOLUTE_ERROR)
+            except AssertionError:
+                status_pass = False
+
+        if 'sigma8' in self.scenario:
+            if not np.isclose(self.scenario['sigma8'], candidate.sigma8(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+            if not np.isclose(candidate.sigma8(), reference.sigma8(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+
+        if 'Neff' in self.scenario:
+            if not np.isclose(self.scenario['Neff'], candidate.Neff(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+            if not np.isclose(candidate.Neff(), reference.Neff(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+
+        if 'theta_s_100' in self.scenario:
+            if not np.isclose(self.scenario['theta_s_100'], candidate.theta_s_100(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+            if not np.isclose(candidate.theta_s_100(), reference.theta_s_100(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
         return status_pass
 
     def store_ini_file(self, path):

source/fourier.c

@@ -1776,7 +1776,7 @@ int fourier_indices(
 
   /** - get list of k values */
 
-  class_call(fourier_get_k_list(ppr,ppt,pfo),
+  class_call(fourier_get_k_list(ppr,ppm,ppt,pfo),
              pfo->error_message,
              pfo->error_message);
 
@@ -1857,6 +1857,7 @@ int fourier_indices(
 
 int fourier_get_k_list(
                        struct precision *ppr,
+                       struct primordial * ppm,
                        struct perturbations * ppt,
                        struct fourier * pfo
                        ) {
@@ -1881,6 +1882,7 @@ int fourier_get_k_list(
                "could not reach extrapolated value k = %.10e starting from k = %.10e with k_per_decade of %.10e in _MAX_NUM_INTERPOLATION_=%i steps",
                ppr->hmcode_max_k_extra,k_max,ppr->k_per_decade_for_pk,_MAX_NUM_EXTRAPOLATION_
                );
+
     pfo->k_size_extra = pfo->k_size+index_k;
   }
   /** - otherwise, same number of values as in perturbation module */
@@ -1906,6 +1908,13 @@ int fourier_get_k_list(
     pfo->ln_k[index_k] = log(k) + exponent*log(10.);
   }
 
+  class_test(pfo->k[pfo->k_size_extra-1]>exp(ppm->lnk[ppm->lnk_size-1]) && ppm->primordial_spec_type != analytic_Pk,
+             pfo->error_message,
+             "Setting the output to HMcode with a large 'hmcode_max_k_extra' and using the primordial spectrum to not analytic is incompatible. Either use the analytic power spectrum or set a smaller 'hmcode_max_k_extra' (k_max_hmcode=%.5e , k_max_primordial=%.5e)",
+             pfo->k[pfo->k_size_extra-1],
+             exp(ppm->lnk[ppm->lnk_size-1])
+             )
+
   return _SUCCESS_;
 }