From b15e761e6559a67c261e1559d6010df29e9054a9 Mon Sep 17 00:00:00 2001 From: Thomas Tram Date: Fri, 4 Oct 2024 18:29:52 +0200 Subject: [PATCH] [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 --- .github/workflows/test_nightly.yml | 2 +- .github/workflows/test_on_pull_request.yml | 2 +- include/fourier.h | 1 + python/test_class.py | 99 ++++++++++++++++++---- source/fourier.c | 11 ++- 5 files changed, 95 insertions(+), 20 deletions(-) diff --git a/.github/workflows/test_nightly.yml b/.github/workflows/test_nightly.yml index f32773617..87bf7d711 100644 --- a/.github/workflows/test_nightly.yml +++ b/.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: | diff --git a/.github/workflows/test_on_pull_request.yml b/.github/workflows/test_on_pull_request.yml index ad0e6cb5d..f6cc3145b 100644 --- a/.github/workflows/test_on_pull_request.yml +++ b/.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() diff --git a/include/fourier.h b/include/fourier.h index af8356d4d..e63136cfc 100644 --- a/include/fourier.h +++ b/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 ); diff --git a/python/test_class.py b/python/test_class.py index 4e0b6baf6..a28b0f066 100644 --- a/python/test_class.py +++ b/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): diff --git a/source/fourier.c b/source/fourier.c index 392b4983c..1d2186ee1 100644 --- a/source/fourier.c +++ b/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_; }