fixed todos

gravity_smg/gravity_models_smg.c

@@ -43,13 +43,8 @@ int gravity_models_gravity_properties_smg(
   *  -> real models: "parameters_smg" to pba->parameters_smg
   *  -> parameterizations: "parameters_smg" to pba->parameters_2_smg
   *                        "expansion_smg" to pba->parameters_smg
-  * NOTE: can change class_read_list_of_doubles_or_default <-> class_read_list_of_doubles
-  * to make it mandatory or allow for default values
   */
 
-  // TODO_EB: if parameters_smg is not specified it gives segfault. There are two cases:
-  // - gravity_model is not specified either. In this case we redirect to propto_omega. We should have default values also for the parameters
-  // - gravity_model is specified. In this case it should return a documented error rather than segfault.
   if (strcmp(string1,"propto_omega") == 0) {
      pba->gravity_model_smg = propto_omega;
      pba->field_evolution_smg = _FALSE_;
@@ -525,7 +520,7 @@ int gravity_models_expansion_properties_smg(
     flag2=_TRUE_;
     pba->parameters_size_smg = 1;
     pba->rho_evolution_smg=_FALSE_;
-    class_read_list_of_doubles_or_default("expansion_smg",pba->parameters_smg,0.0,pba->parameters_size_smg);
+    class_read_list_of_doubles("expansion_smg",pba->parameters_smg,pba->parameters_size_smg);
   }
   //accept different names
 
@@ -534,23 +529,24 @@ int gravity_models_expansion_properties_smg(
     flag2=_TRUE_;
     pba->parameters_size_smg = 3;
     pba->rho_evolution_smg=_FALSE_;
-    class_read_list_of_doubles_or_default("expansion_smg",pba->parameters_smg,0.0,pba->parameters_size_smg);
+    class_read_list_of_doubles("expansion_smg",pba->parameters_smg,pba->parameters_size_smg);
   }
 
   if (strcmp(string1,"wowa_w") == 0 || strcmp(string1,"w0wa_w") == 0 || strcmp(string1,"cpl_w") == 0 ) {
     pba->expansion_model_smg = wowa_w;
     flag2=_TRUE_;
     pba->parameters_size_smg = 3;
     pba->rho_evolution_smg=_TRUE_;
-    class_read_list_of_doubles_or_default("expansion_smg",pba->parameters_smg,0.0,pba->parameters_size_smg);
+    class_read_list_of_doubles("expansion_smg",pba->parameters_smg,pba->parameters_size_smg);
   }
 
   if (strcmp(string1,"wede") == 0) {
     //ILSWEDE
     pba->expansion_model_smg = wede;
     flag2=_TRUE_;
     pba->parameters_size_smg = 3;
-    class_read_list_of_doubles_or_default("expansion_smg",pba->parameters_smg,0.0,pba->parameters_size_smg);
+    class_read_list_of_doubles("expansion_smg",pba->parameters_smg,pba->parameters_size_smg);
+    class_read_double("wede_Omega_e_regularizer_smg",pba->wede_Omega_e_regularizer_smg);
     // 	//optimize the guessing BUG: eventually leads to problem in the MCMC, perhaps the guess is too good?
     // 	if(pba->tuning_index_smg == 0){
     // 	  pba->parameters_smg[0] = pba->Omega0_smg;
@@ -787,9 +783,7 @@ int gravity_models_get_back_par_smg(
 
     double Om0 = pba->parameters_smg[0];
     double w0 = pba->parameters_smg[1];
-    // TODO_EB: add a precision parameter for 1e-10
-    //NOTE: I've regularized the expression adding a tiny Omega_e
-    double Ome = pba->parameters_smg[2] + 1e-10;
+    double Ome = pba->parameters_smg[2] + pba->wede_Omega_e_regularizer_smg;
 
     double Om = ((Om0 - Ome*(1.-pow(a,-3.*w0)))/(Om0 + (1.-Om0)*pow(a,3*w0)) + Ome*(1.-pow(a,-3*w0)));
     double dOm_da = (3*pow(a,-1 - 3*w0)*(-1 + Om0)*(-2*pow(a,3*w0)*(-1 + Om0)*Ome + Om0*Ome + pow(a,6*w0)*(Om0 - 2*Ome + Om0*Ome))*w0)/pow(-(pow(a,3*w0)*(-1 + Om0)) + Om0,2); //from Mathematica

gravity_smg/input_smg.c

@@ -165,16 +165,14 @@ int input_read_parameters_smg(
   class_call(parser_read_string(pfc,"gravity_model",&string1,&flag1,errmsg),
              errmsg,
              errmsg);
-
-  if (flag1 == _FALSE_) {
-    printf(" gravity_model not read, default will be used \n");
-  }
-  else {
-    /** Fix flags for each gravity model */
-     class_call(gravity_models_gravity_properties_smg(pfc, pba, string1, has_tuning_index_smg, has_dxdy_guess_smg, errmsg),
-                errmsg,
-                errmsg);
-  }
+  class_test(flag1 == _FALSE_,
+              errmsg,
+              "gravity_model not read, you should specify one!\n");
+
+  /** Fix flags for each gravity model */
+    class_call(gravity_models_gravity_properties_smg(pfc, pba, string1, has_tuning_index_smg, has_dxdy_guess_smg, errmsg),
+              errmsg,
+              errmsg);
   // end of loop over models
 
   if(pba->field_evolution_smg == _TRUE_){
@@ -192,8 +190,9 @@ int input_read_parameters_smg(
     class_call(parser_read_string(pfc,"expansion_model",&string1,&flag1,errmsg),
                errmsg,
                errmsg);
-    if (flag1 == _FALSE_)
-      printf("No expansion model specified, will take default one \n");
+    class_test(flag1 == _FALSE_,
+                errmsg,
+                "expansion_model not read, you should specify one!\n");
 
     /** Fix flags for each expansion model */
     class_call(gravity_models_expansion_properties_smg(pfc, pba, string1, errmsg),
@@ -337,6 +336,7 @@ int input_default_params_smg(
   pba->hubble_evolution = _TRUE_; /** dynamical evolution of Friedmann eq. */
   pba->hubble_friction = 3.; /** friction coefficient in H' equation: H' = ... + H_friction*(H^2 - rho_crit) [NOT ONLY IN SMG!] */
   pba->rho_evolution_smg = _FALSE_; /*does the model require to evolve the background energy density? (only for parameterizations)*/
+  pba->wede_Omega_e_regularizer_smg = 1.e-10; /** regularize adding a tiny Omega_e */
 
   pba->kineticity_safe_smg = 0; /* value added to the kineticity, useful to cure perturbations at early time in some models */
   pba->cs2_safe_smg = 0; /* threshold to consider the sound speed of scalars negative in the stability check */

hi_class.ini

@@ -292,7 +292,7 @@ n_max_qs_smg = 1e4
 
 start_small_k_at_tau_c_over_tau_h = 1e-4
 start_large_k_at_tau_h_over_tau_k = 1e-4
-perturb_sampling_stepsize = 0.05
+perturbations_sampling_stepsize = 0.05
 l_logstep = 1.045
 l_linstep = 50
 
@@ -310,20 +310,26 @@ output = tCl,pCl,lCl,mPk
 
 8b) file name root 'root' for all output files
 
-root = output/test_
+root = output/test
+overwrite_root = yes
 
 8c) Do you want to write a file with the input parameters you used? Do you want to write a table
     of background quantitites in a file?
 
-write parameters = yeap
-write background = yeah
+write_parameters = yeap
+write_background = yeah
 k_output_values = 0.1, 0.01, 0.0001
 
 
 8d) Verbose parameters
 input_verbose = 1
 background_verbose = 1
-output_verbose = 1
 thermodynamics_verbose = 1
 perturbations_verbose = 1
-spectra_verbose = 1
+transfer_verbose = 1
+primordial_verbose = 1
+harmonic_verbose = 1
+fourier_verbose = 1
+lensing_verbose = 1
+distortions_verbose = 1
+output_verbose = 1

include/background.h

@@ -149,6 +149,7 @@ struct background
 
   double hubble_friction; /** friction coefficient in H' equation: H' = ... + H_friction*(H^2 - rho_crit) [NOT ONLY IN SMG!] */
   int hubble_evolution; /** whether to evolve H' from the equation */
+  double wede_Omega_e_regularizer_smg; /** regularize adding a tiny Omega_e */
 
   enum gravity_model gravity_model_smg; /** Horndeski model */
   //   enum gravity_model_subclass gravity_submodel_smg; /** Horndeski model */

source/background.c

@@ -624,8 +624,7 @@ int background_functions(
 
   /** - compute critical density */
   rho_crit = rho_tot-pba->K/a/a;
-  // TODO_EB: here the condition in hi_class was (rho_crit <= 0.) && (pba->has_smg == _FALSE_). Why?
-  class_test(rho_crit <= 0.,
+  class_test((rho_crit <= 0.) && (pba->initial_conditions_set_smg == _TRUE_),
              pba->error_message,
              "rho_crit = %e instead of strictly positive",rho_crit);
 
@@ -1975,9 +1974,10 @@ int background_solve(
   class_alloc(pvecback_integration,pba->bi_size*sizeof(double),pba->error_message);
 
   /** - impose initial conditions with background_initial_conditions() */
-  class_call(background_initial_conditions(ppr,pba,pvecback,pvecback_integration,&(loga_ini)),
+  class_call_except(background_initial_conditions(ppr,pba,pvecback,pvecback_integration,&(loga_ini)),
              pba->error_message,
-             pba->error_message);
+             pba->error_message,
+             free(pvecback);free(pvecback_integration););
 
   /** - Determine output vector */
   loga_final = 0.; // with our conventions, loga is in fact log(a/a_0); we integrate until today, when log(a/a_0) = 0
@@ -2403,7 +2403,6 @@ int background_initial_conditions(
 
   /* Just checking that our initial time indeed is deep enough in the radiation
      dominated regime (_smg) */
-  // TODO_EB: There was a class_test_except with free(pvecback);free(pvecback_integration);background_free(pba)
   class_test(fabs(pvecback[pba->index_bg_Omega_r]+pvecback[pba->index_bg_Omega_de]-1.) > ppr->tol_initial_Omega_r,
              pba->error_message,
              "Omega_r = %e, Omega_de = %e, not close enough to 1. Decrease a_ini_over_a_today_default in order to start from radiation domination.",