From 1e3fb2ce70d6a2f497c51422b0c09247428c4a5a Mon Sep 17 00:00:00 2001 From: lesgourg Date: Fri, 4 Oct 2024 11:50:35 +0200 Subject: [PATCH 1/4] 186 fix real space interface (#187) made RealSpaceInterface work again, credits Timothy Morton and Mowen Zhao --- explanatory.ini | 3 ++ .../Calc2D/CalculationClass.py | 4 +- .../RealSpaceInterface/Calc2D/Database.py | 18 +++---- .../Calc2D/TransferFunction.py | 1 + external/RealSpaceInterface/README | 9 +++- external/RealSpaceInterface/tornadoserver.py | 9 ++-- include/perturbations.h | 2 + source/input.c | 9 +++- source/perturbations.c | 50 ++++++++++--------- 9 files changed, 62 insertions(+), 43 deletions(-) diff --git a/explanatory.ini b/explanatory.ini index b483cd901..c74f297bd 100644 --- a/explanatory.ini +++ b/explanatory.ini @@ -151,6 +151,9 @@ gauge = synchronous # (default: set to 'no') #nbody_gauge_transfer_functions = yes +# 4.c) Do you want the source functions for total non-relativistic matter, delta_m and theta_m, and baryon+cdm, delta_cb and theta_cb, to be outputed in the current gauge (the one selected in 4.a), instead of being automatically expressed as a gauge-invariant (GI) quantity, likeL: delta_m^GI=delta_m + 3 a H theta_m/k2, theta_m^GI=theta_m+alpha*k2 (default: no, that is, convert to GI) +#matter_source_in_current_gauge= no + # 5) Hubble parameter : either 'H0' in km/s/Mpc or 'h' (or 'theta_s_100'), where # the latter is the peak scale parameter defined exactly as 100(ds_dec/da_dec) # with a decoupling time given by maximum of visibility function (quite different diff --git a/external/RealSpaceInterface/Calc2D/CalculationClass.py b/external/RealSpaceInterface/Calc2D/CalculationClass.py index 0a6b5bb40..c63af4204 100644 --- a/external/RealSpaceInterface/Calc2D/CalculationClass.py +++ b/external/RealSpaceInterface/Calc2D/CalculationClass.py @@ -71,7 +71,7 @@ def getData(self, redshiftindex): Valuenew = dict() FValue_abs = np.abs(self.FValue) _min, _max = FValue_abs.min(), FValue_abs.max() - dimensions = (self.endshape[0] / 2, self.endshape[1]) + dimensions = (self.endshape[0] // 2, self.endshape[1]) for quantity, FT in FValuenew.items(): FT_abs = np.abs(FT) FT_normalized = cv2.resize(FT_abs, dimensions).ravel() @@ -98,7 +98,7 @@ def getInitialData(self): return Value.ravel(), cv2.resize( (np.abs(self.FValue) - np.abs(self.FValue).min()) / (np.abs(self.FValue).max() - np.abs(self.FValue).min()), - (self.endshape[0] / 2, self.endshape[1])).ravel(), (minimum, + (self.endshape[0] // 2, self.endshape[1])).ravel(), (minimum, maximum) def getTransferData(self, redshiftindex): diff --git a/external/RealSpaceInterface/Calc2D/Database.py b/external/RealSpaceInterface/Calc2D/Database.py index bc8f006c2..af5e4176c 100644 --- a/external/RealSpaceInterface/Calc2D/Database.py +++ b/external/RealSpaceInterface/Calc2D/Database.py @@ -14,22 +14,22 @@ def __init__(self, directory, db_file="database.dat"): self.db_path = os.path.join(directory, db_file) if not os.path.exists(self.db_path): logging.info("No database found; Creating one at {}.".format(self.db_path)) - with open(self.db_path, "w") as f: + with open(self.db_path, "wb") as f: pickle.dump(dict(), f) self.db = self.__read_database() def __read_database(self): - with open(self.db_path) as f: - return pickle.load(f) + with open(self.db_path, "rb") as f: + return pickle.load(f, encoding='latin1') def __write_database(self): - with open(self.db_path, "w") as f: - pickle.dump(self.db, f) + with open(self.db_path, "wb") as f: + pickle.dump(self.db, f, protocol=pickle.HIGHEST_PROTOCOL) def __create_file(self, data): filename = str(uuid.uuid4()) - with open(os.path.join(self.directory, filename), "w") as f: + with open(os.path.join(self.directory, filename), "wb") as f: pickle.dump(data, f) return filename @@ -40,8 +40,8 @@ def __getitem__(self, key): frozen_key = self.__get_frozen_key(key) if frozen_key in self.db: filename = self.db[frozen_key] - with open(os.path.join(self.directory, filename)) as f: - return pickle.load(f) + with open(os.path.join(self.directory, filename), "rb") as f: + return pickle.load(f, encoding='latin1') else: raise KeyError("No data for key: {}".format(key)) @@ -55,4 +55,4 @@ def __contains__(self, key): Return whether `self` contains a record for the given `key`. """ - return self.__get_frozen_key(key) in self.db \ No newline at end of file + return self.__get_frozen_key(key) in self.db diff --git a/external/RealSpaceInterface/Calc2D/TransferFunction.py b/external/RealSpaceInterface/Calc2D/TransferFunction.py index 8747f007f..135bc195f 100644 --- a/external/RealSpaceInterface/Calc2D/TransferFunction.py +++ b/external/RealSpaceInterface/Calc2D/TransferFunction.py @@ -42,6 +42,7 @@ def ComputeTransferFunctionList(cosmologicalParameters, redshift, kperdecade=200 class_settings.update({ "output": "mTk", "gauge": "newtonian", + "matter_source_in_current_gauge": "yes", "evolver": "1", "P_k_max_h/Mpc": P_k_max, "k_per_decade_for_pk": kperdecade, diff --git a/external/RealSpaceInterface/README b/external/RealSpaceInterface/README index 4fc84d3e3..975f58989 100644 --- a/external/RealSpaceInterface/README +++ b/external/RealSpaceInterface/README @@ -1,8 +1,14 @@ +CLASS real space interface +07.2015: started by Max Beutelspacher +09.2018: improved by Georgios Samaras and released in class v2.7.0 +10.2025: for v3.2.5: made compatible with python 3, credits Thimothy Morton and Mowen Zhao + + For installation of python packages, run pip install -r requirements.txt -Launch the application with +To see the real space interface, go to this directory (external/RealSpaceInterface/) and launch the application with python tornadoserver.py @@ -15,4 +21,3 @@ Then in any browser open the URL Cache files are located in cache/, so to clear the cache, run rm cache/* - diff --git a/external/RealSpaceInterface/tornadoserver.py b/external/RealSpaceInterface/tornadoserver.py index 808119e25..299369ba4 100644 --- a/external/RealSpaceInterface/tornadoserver.py +++ b/external/RealSpaceInterface/tornadoserver.py @@ -143,8 +143,8 @@ def send_frame(self, redindex): self.write_message(json.dumps({'type': 'extrema', 'extrema': extrema})) progress = float(redindex) / len(self.calc.redshift) - real = {quantity: base64.b64encode(data.astype(np.float32)) for quantity, data in Valuenew.iteritems()} - transfer = {quantity: base64.b64encode(data.astype(np.float32)) for quantity, data in TransferData.iteritems()} + real = {quantity: base64.b64encode(data.astype(np.float32)).decode() for quantity, data in Valuenew.items()} + transfer = {quantity: base64.b64encode(data.astype(np.float32)).decode() for quantity, data in TransferData.items()} self.write_message( json.dumps({ 'type': 'data', @@ -166,9 +166,9 @@ def send_initial_state(self): extremastring = json.dumps({'type': 'extrema', 'extrema': extrema}) datastring = json.dumps({ 'type': 'data', - 'real': base64.b64encode(Value.astype(np.float32)), + 'real': base64.b64encode(Value.astype(np.float32)).decode(), 'fourier': [], - 'transfer': base64.b64encode(TransferData.astype(np.float32)), + 'transfer': base64.b64encode(TransferData.astype(np.float32)).decode(), 'k': krange.tolist() }) self.write_message(extremastring) @@ -212,6 +212,7 @@ def set_cosmological_parameters(self, cosmologicalParameters): z_of_decoupling = self.calc.z_dec frame_of_decoupling = np.argmin(np.abs(z_of_decoupling - self.calc.redshift)) + frame_of_decoupling = int(frame_of_decoupling) # this is np.int64 object, which json doesn't like if self.calc.redshift[frame_of_decoupling] > z_of_decoupling: frame_of_decoupling -= 1 messages.append({ diff --git a/include/perturbations.h b/include/perturbations.h index bf4c4e976..be28af61e 100644 --- a/include/perturbations.h +++ b/include/perturbations.h @@ -198,6 +198,8 @@ struct perturbations enum possible_gauges gauge; /**< gauge in which to perform this calculation */ + short has_matter_source_in_current_gauge; /**< whether to keep matter and baryon+CDM sources in current gauge, instead of automatic conversion to gauge-invariant variables */ + //@} /** @name - indices running on modes (scalar, vector, tensor) */ diff --git a/source/input.c b/source/input.c index 581946223..4eef79942 100644 --- a/source/input.c +++ b/source/input.c @@ -757,7 +757,7 @@ int input_shooting(struct file_content * pfc, /* Create file content structure with additional entries */ class_call(parser_extend(pfc, 1, errmsg), errmsg,errmsg); - + class_call(parser_init_from_pfc(pfc, &(fzw.fc), errmsg), errmsg,errmsg); @@ -824,7 +824,7 @@ int input_shooting(struct file_content * pfc, fzw.fc.name[pfc->size - 1], fzw.fc.value[pfc->size - 1]); } - + parser_copy(&(fzw.fc), pfc, pfc->size - 1, pfc->size); /** Free arrays allocated */ @@ -2077,6 +2077,9 @@ int input_read_parameters_general(struct file_content * pfc, } + /** 4.c) Do we want matter and baryon+CDM sources in current gauge, instead of automatic conversion to gauge-invariant variables? */ + class_read_flag("matter_source_in_current_gauge",ppt->has_matter_source_in_current_gauge); + /** 5) h in [-] and H_0/c in [1/Mpc = h/2997.9 = h*10^5/c] */ /* Read */ class_call(parser_read_double(pfc,"H0",¶m1,&flag1,errmsg), @@ -5658,6 +5661,8 @@ int input_default_params(struct background *pba, ppt->gauge=synchronous; /** 4.b) N-body gauge */ ppt->has_Nbody_gauge_transfers = _FALSE_; + /** 4.c) keep delta_m, theta_m, delta_cb, theta_cb in current gauge */ + ppt->has_matter_source_in_current_gauge = _FALSE_; /** 5) Hubble parameter */ pba->h = 0.67810; diff --git a/source/perturbations.c b/source/perturbations.c index d88662b21..9eaf301ce 100644 --- a/source/perturbations.c +++ b/source/perturbations.c @@ -866,9 +866,9 @@ int perturbations_init( ppt->z_max_pk, pth->z_rec); - class_test(ppt->has_source_delta_m == _TRUE_, + class_test((ppt->has_source_delta_m == _TRUE_) && (ppt->has_matter_source_in_current_gauge == _FALSE_), ppt->error_message, - "You requested a very high z_pk=%e, higher than z_rec=%e. This works very well when you ask only transfer functions, e.g. with 'output=mTk' or 'output=mTk,vTk'. But if you need the total matter (e.g. with 'mPk', 'dCl', etc.) there is an issue with the calculation of delta_m at very early times. By default, delta_m is a gauge-invariant variable (the density fluctuation in comoving gauge) and this quantity is hard to get accurately at very early times. The solution is to define delta_m as the density fluctuation in the current gauge, synchronous or newtonian. For the moment this must be done manually by commenting the line 'ppw->delta_m += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_m/k2;' in perturbations_sources(). In the future there will be an option for doing it in an easier way.", + "You requested a very high z_pk=%e, higher than z_rec=%e. This works very well when you ask only transfer functions, e.g. with 'output=mTk' or 'output=mTk,vTk'. But if you need the total matter (e.g. with 'mPk', 'dCl', etc.) there is an issue with the calculation of delta_m (or delta_cb) at very early times. By default, delta_m and delta_cb are expressed by the code as gauge-invariant variables (the density fluctuation in comoving gauge). This quantity is hard to get accurately at very early times. The solution is to define delta_m and delta_cb as density fluctuations in the current gauge, e.g. synchronous or newtonian. This is done by setting the input flag 'matter_source_in_current_gauge' to 'yes'.", ppt->z_max_pk, pth->z_rec); @@ -6610,31 +6610,33 @@ int perturbations_einstein( gauge-independent variables (you could comment this out if you really want gauge-dependent results) */ - if (ppt->has_source_delta_m == _TRUE_) { - ppw->delta_m += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_m/k2; - // note: until 2.4.3 there was a typo, the factor was (-2 H'/H) instead - // of (3 aH). There is the same typo in the CLASSgal paper - // 1307.1459v1,v2,v3. It came from a confusion between (1+w_total) - // and (1+w_matter)=1 [the latter is the relevant one here]. - // - // note2: at this point this gauge-invariant variable is only - // valid if all matter components are pressureless and - // stable. This relation will be generalized soon to the case - // of decaying dark matter. - } + if (ppt->has_matter_source_in_current_gauge == _FALSE_) { + if (ppt->has_source_delta_m == _TRUE_) { + ppw->delta_m += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_m/k2; + // note: until 2.4.3 there was a typo, the factor was (-2 H'/H) instead + // of (3 aH). There is the same typo in the CLASSgal paper + // 1307.1459v1,v2,v3. It came from a confusion between (1+w_total) + // and (1+w_matter)=1 [the latter is the relevant one here]. + // + // note2: at this point this gauge-invariant variable is only + // valid if all matter components are pressureless and + // stable. This relation will be generalized soon to the case + // of decaying dark matter. + } - if (ppt->has_source_delta_cb == _TRUE_) { - ppw->delta_cb += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_cb/k2;//check gauge transformation - } + if (ppt->has_source_delta_cb == _TRUE_) { + ppw->delta_cb += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_cb/k2;//check gauge transformation + } - if (ppt->has_source_theta_m == _TRUE_) { - if (ppt->gauge == synchronous) { - ppw->theta_m += ppw->pvecmetric[ppw->index_mt_alpha]*k2; + if (ppt->has_source_theta_m == _TRUE_) { + if (ppt->gauge == synchronous) { + ppw->theta_m += ppw->pvecmetric[ppw->index_mt_alpha]*k2; + } } - } - if (ppt->has_source_theta_cb == _TRUE_){ - if (ppt->gauge == synchronous) { - ppw->theta_cb += ppw->pvecmetric[ppw->index_mt_alpha]*k2; //check gauge transformation + if (ppt->has_source_theta_cb == _TRUE_){ + if (ppt->gauge == synchronous) { + ppw->theta_cb += ppw->pvecmetric[ppw->index_mt_alpha]*k2; //check gauge transformation + } } } } From 2c7e63b1eef363e6f3a7d969b240105de8175860 Mon Sep 17 00:00:00 2001 From: Nils <32747859+schoeneberg@users.noreply.github.com> Date: Fri, 4 Oct 2024 17:36:31 +0200 Subject: [PATCH 2/4] 188 shooting early stopping (#189) * First part of bugfix: changing abs to fabs, since in C only fabs handles floating point numbers, and abs only integers * Only allow early stopping if the derivative is large enough * Fixing spelling mistake that was somehow accidentially commited * Fixed bug where shooting failing was overwritten at later stages by re-initializing the shooting_failed parameter to its default value * More sensible shooting outputs * More explicit about Neff < N_eff_min inputs * Removing early stopping again, to make code more elegant and avoid potential issues (and the slowdown should be very small anyways) --------- Co-authored-by: schoeneberg --- source/input.c | 64 ++++++++++++++++++++++++++++++-------------------- 1 file changed, 38 insertions(+), 26 deletions(-) diff --git a/source/input.c b/source/input.c index 4eef79942..954311d6f 100644 --- a/source/input.c +++ b/source/input.c @@ -412,6 +412,9 @@ int input_read_from_file(struct file_content * pfc, class_read_int("input_verbose",input_verbose); if (input_verbose >0) printf("Reading input parameters\n"); + /** -- Special setting of parameter, before anything else: did shooting fail? */ + pba->shooting_failed = _FALSE_; + /** Find out if shooting necessary and, eventually, shoot and initialize read parameters */ class_call(input_shooting(pfc,ppr,pba,pth,ppt,ptr,ppm,phr,pfo,ple,psd,pop, @@ -641,10 +644,9 @@ int input_shooting(struct file_content * pfc, /* We can do 1 dimensional root finding */ if (input_verbose > 0) { - fprintf(stdout, - "Computing unknown input parameter '%s' using input parameter '%s'\n", - fzw.fc.name[fzw.unknown_parameters_index[0]], - target_namestrings[fzw.target_name[0]]); + printf("Computing unknown input parameter '%s' using input parameter '%s'\n", + fzw.fc.name[fzw.unknown_parameters_index[0]], + target_namestrings[fzw.target_name[0]]); } /* If shooting fails, postpone error to background module to play nice with MontePython. */ @@ -662,9 +664,14 @@ int input_shooting(struct file_content * pfc, // precision of around 1e-16, so 1e-20 should be good enough for the shooting class_sprintf(fzw.fc.value[fzw.unknown_parameters_index[0]],"%.20e",xzero); if (input_verbose > 0) { - fprintf(stdout," -> found '%s = %s'\n", - fzw.fc.name[fzw.unknown_parameters_index[0]], - fzw.fc.value[fzw.unknown_parameters_index[0]]); + if (shooting_failed == _FALSE_){ + printf(" -> found '%s = %s'\n", + fzw.fc.name[fzw.unknown_parameters_index[0]], + fzw.fc.value[fzw.unknown_parameters_index[0]]); + } + else{ + printf("Shooting failed! Aborting...\n"); + } } } @@ -673,7 +680,7 @@ int input_shooting(struct file_content * pfc, /* We need to do multidimensional root finding */ if (input_verbose > 0) { - fprintf(stdout,"Computing unknown input parameters\n"); + printf("Computing unknown input parameters\n"); } /* Allocate local variables */ @@ -710,9 +717,14 @@ int input_shooting(struct file_content * pfc, class_sprintf(fzw.fc.value[fzw.unknown_parameters_index[counter]], "%.20e",x_inout[counter]); if (input_verbose > 0) { - fprintf(stdout," -> found '%s = %s'\n", - fzw.fc.name[fzw.unknown_parameters_index[counter]], - fzw.fc.value[fzw.unknown_parameters_index[counter]]); + if (shooting_failed == _FALSE_){ + printf(" -> found '%s = %s'\n", + fzw.fc.name[fzw.unknown_parameters_index[counter]], + fzw.fc.value[fzw.unknown_parameters_index[counter]]); + } + else{ + printf("Shooting failed! Aborting...\n"); + } } } @@ -721,8 +733,8 @@ int input_shooting(struct file_content * pfc, free(dxdF); } - if (input_verbose > 1) { - fprintf(stdout,"Shooting completed using %d function evaluations\n",fevals); + if (input_verbose > 1 && shooting_failed == _FALSE_) { + printf("Shooting completed using %d function evaluations\n",fevals); } /** Set status of shooting */ @@ -790,10 +802,9 @@ int input_shooting(struct file_content * pfc, /* Print to the user */ if (input_verbose > 0) { - fprintf(stdout, - "Computing unknown input parameter '%s' using input parameter '%s'\n", - (flag1 ==_TRUE_?"sigma8":"S8"), - "A_s"); + printf("Computing unknown input parameter '%s' using input parameter '%s'\n", + (flag1 ==_TRUE_?"sigma8":"S8"), + "A_s"); } /* Set a guess for A_s from LCDM (doesn't need to be super accurate) */ @@ -820,9 +831,9 @@ int input_shooting(struct file_content * pfc, /* Store the derived value with high enough accuracy */ class_sprintf(fzw.fc.value[pfc->size - 1],"%.20e",A_s); if (input_verbose > 0) { - fprintf(stdout," -> found '%s = %s'\n", - fzw.fc.name[pfc->size - 1], - fzw.fc.value[pfc->size - 1]); + printf(" -> found '%s = %s'\n", + fzw.fc.name[pfc->size - 1], + fzw.fc.value[pfc->size - 1]); } parser_copy(&(fzw.fc), pfc, pfc->size - 1, pfc->size); @@ -868,7 +879,6 @@ int input_needs_shooting_for_target(struct file_content * pfc, case Omega_scf: case Omega_ini_dcdm: case omega_ini_dcdm: - case Neff: /* Check that Omega's or omega's are nonzero: */ if (target_value == 0.) *needs_shooting = _FALSE_; @@ -926,8 +936,6 @@ int input_find_root(double *xzero, /* Try fifteen times to go above and below the root (i.e. where shooting succeeds) */ for (iter=1; iter<=15; iter++){ x2 = x1 - dx; - // Early stopping if the two points are infinitessimaly close together, i.e. the answer is already reached - if(abs(x2/x1-1) < tol_x_rel){*xzero = x1; return _SUCCESS_;} /* Try three times to get a 'reasonable' value, i.e. no CLASS error */ for (iter2=1; iter2 <= 3; iter2++) { return_function = input_fzerofun_1d(x2, pfzw, &f2, errmsg); @@ -1154,6 +1162,9 @@ int input_get_guess(double *xguess, /* Cheat to read only known parameters: */ pfzw->fc.size -= pfzw->target_size; + /* Assume for now shooting did not fail */ + ba.shooting_failed = _FALSE_; + class_call(input_read_precisions(&(pfzw->fc),&pr,&ba,&th,&pt,&tr,&pm,&hr,&fo,&le,&sd,&op, errmsg), errmsg, @@ -1319,6 +1330,9 @@ int input_try_unknown_parameters(double * unknown_parameter, int param; short compute_sigma8 = _FALSE_; + /* Assume for now shooting did not fail */ + ba.shooting_failed = _FALSE_; + pfzw = (struct fzerofun_workspace *) voidpfzw; /** Read input parameters */ // This needs to be done with enough accuracy. A standard double has a relative @@ -2433,7 +2447,7 @@ int input_read_parameters_species(struct file_content * pfc, pba->Omega0_ur = param3/pba->h/pba->h; } } - class_test(pba->Omega0_ur<0,errmsg,"You cannot set the density of ultra-relativistic relics (dark radiation/neutrinos) to negative values."); + class_test(pba->Omega0_ur<0,errmsg,"You cannot set the density of ultra-relativistic relics (dark radiation/neutrinos) to negative values. You might have input a total Neff smaller than what your massive neutrinos require minimally (around 1.02 * N_ncdm * deg_ncdm)."); /** 3.a) Case of non-standard properties */ /* Read */ @@ -5831,8 +5845,6 @@ int input_default_params(struct background *pba, pba->phi_prime_ini_scf = 1; // factors of the radiation attractor values /** 9.b.3) Tuning parameter */ pba->scf_tuning_index = 0; - /** 9.b.4) Shooting parameter */ - pba->shooting_failed = _FALSE_; /** * Deafult to input_read_parameters_heating From b15e761e6559a67c261e1559d6010df29e9054a9 Mon Sep 17 00:00:00 2001 From: Thomas Tram Date: Fri, 4 Oct 2024 18:29:52 +0200 Subject: [PATCH 3/4] [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_; } From 2febfbefaf996d81fd139c9243863fbc13a96fa2 Mon Sep 17 00:00:00 2001 From: lesgourg Date: Mon, 7 Oct 2024 12:57:28 +0200 Subject: [PATCH 4/4] v3.2.5, negative m_ncdm forbidden --- external/RealSpaceInterface/README | 2 +- include/common.h | 2 +- source/input.c | 14 ++++++++++++-- 3 files changed, 14 insertions(+), 4 deletions(-) diff --git a/external/RealSpaceInterface/README b/external/RealSpaceInterface/README index 975f58989..83b3ca585 100644 --- a/external/RealSpaceInterface/README +++ b/external/RealSpaceInterface/README @@ -1,7 +1,7 @@ CLASS real space interface 07.2015: started by Max Beutelspacher 09.2018: improved by Georgios Samaras and released in class v2.7.0 -10.2025: for v3.2.5: made compatible with python 3, credits Thimothy Morton and Mowen Zhao +10.2025: for v3.2.5: made compatible with python 3, credits Timothy Morton and Mowen Zhao For installation of python packages, run diff --git a/include/common.h b/include/common.h index 2100f1cb6..699fcd18c 100644 --- a/include/common.h +++ b/include/common.h @@ -15,7 +15,7 @@ #ifndef __COMMON__ #define __COMMON__ -#define _VERSION_ "v3.2.4" +#define _VERSION_ "v3.2.5" /* @cond INCLUDE_WITH_DOXYGEN */ diff --git a/source/input.c b/source/input.c index 954311d6f..f35aa82cc 100644 --- a/source/input.c +++ b/source/input.c @@ -2560,18 +2560,28 @@ int input_read_parameters_species(struct file_content * pfc, errmsg, errmsg); - /** 5.d) Mass or Omega of each ncdm species */ + /** 5.d) Mass and/or Omega of each ncdm species */ /* Read */ class_read_list_of_doubles_or_default("m_ncdm",pba->m_ncdm_in_eV,0.0,N_ncdm); + for (n=0; nm_ncdm_in_eV[n]<0, + errmsg, + "You entered a negative non-CDM mass m_ncdm[%d], which makes no sense. This error was not caught in previous CLASS versions because the mass is always squared in the code, so CLASS returned the exact same results form +m_ncdm and -m_ncdm. If you want to define an 'effective negative neutrino mass' in the sense of e.g. 2405.00836 or 2407.10965, you can implement it in a python script following e.g. eq.(3) of 2407.10965",n); + } + class_read_list_of_doubles_or_default("Omega_ncdm",pba->Omega0_ncdm,0.0,N_ncdm); + // the name M_ncdm is borrowed temporarily to store omega_ncdm class_read_list_of_doubles_or_default("omega_ncdm",pba->M_ncdm,0.0,N_ncdm); for (n=0; nM_ncdm[n]!=0.0){ /* Test */ class_test(pba->Omega0_ncdm[n]!=0,errmsg, "You can only enter one of 'Omega_ncdm' or 'omega_ncdm' for ncdm species %d.",n); - /* Complete set of parameters */ + /* Complete set of parameters: if the user passed either + Omega_ncdm or omega_ncdm, now it's stored anyway as + Omega_0_ncdm */ pba->Omega0_ncdm[n] = pba->M_ncdm[n]/pba->h/pba->h; + // the name M_ncdm is now available for its real destination } /* Set default value this is the right place for passing the default value of the mass