allow bbn_predictor to be the string name of a predictor class or int…

…erpolation table

camb/bbn.py

@@ -28,6 +28,8 @@
 n_photon = (kB * TCMB / hbar / c) ** 3 * zeta3 * 2 / np.pi ** 2
 omegafac = (1e5 / Mpc) ** 2 / (8 * np.pi * G) * 3
 
+default_interpolation_table = 'PArthENoPE_880.2_standard.dat'
+
 
 def yhe_to_ypBBN(Yp):
     return -4 * m_H * Yp / (Yp * m_He - 4 * Yp * m_H - m_He)
@@ -76,7 +78,7 @@ class BBN_table_interpolator(BBNPredictor):
 
     """
 
-    def __init__(self, interpolation_table='PArthENoPE_880.2_standard.dat', function_of=['ombh2', 'DeltaN']):
+    def __init__(self, interpolation_table=default_interpolation_table, function_of=['ombh2', 'DeltaN']):
 
         if os.sep not in interpolation_table and '/' not in interpolation_table:
             interpolation_table = os.path.normpath(os.path.join(os.path.dirname(__file__), interpolation_table))
@@ -194,17 +196,23 @@ def DH(self, ombh2, delta_neff, tau_neutron=None):
                ) * pow((tau_neutron or self.taun) / 880.3, 0.418) * 1e-5
 
 
-_default_predictor = None
+_predictors = {}
 
 
-def get_default_predictor():
+def get_predictor(predictor_name=None):
     """
     Get instance of default BBNPredictor class. Currently numerical table interpolation as Planck 2018 analysis.
     """
-    global _default_predictor
-    if _default_predictor is None:
-        _default_predictor = BBN_table_interpolator()
-    return _default_predictor
+    global _predictors
+    predictor_name = predictor_name or default_interpolation_table
+    predictor = _predictors.get(predictor_name, None)
+    if predictor is None:
+        if predictor_name == 'BBN_fitting_parthenope':
+            predictor = BBN_fitting_parthenope()
+        else:
+            predictor = BBN_table_interpolator(interpolation_table=predictor_name)
+        _predictors[predictor_name] = predictor
+    return predictor
 
 
 if __name__ == "__main__":

camb/model.py

@@ -401,9 +401,10 @@ def f(H0):
     def set_cosmology(self, H0=None, ombh2=0.022, omch2=0.12, omk=0.0,
                       cosmomc_theta=None, thetastar=None,
                       neutrino_hierarchy='degenerate', num_massive_neutrinos=1,
-                      mnu=0.06, nnu=constants.default_nnu, YHe=None, meffsterile=0.0, standard_neutrino_neff=constants.default_nnu,
+                      mnu=0.06, nnu=constants.default_nnu, YHe=None, meffsterile=0.0,
+                      standard_neutrino_neff=constants.default_nnu,
                       TCMB=constants.COBE_CMBTemp, tau=None, deltazrei=None, Alens=1.0,
-                      bbn_predictor=None, theta_H0_range=[10, 100]):
+                      bbn_predictor=None, bbn_interpolation_table=None, theta_H0_range=[10, 100]):
         r"""
         Sets cosmological parameters in terms of physical densities and parameters (e.g. as used in Planck analyses).
         Default settings give a single distinct neutrino mass eigenstate, by default one neutrino with mnu = 0.06eV.
@@ -439,14 +440,18 @@ def set_cosmology(self, H0=None, ombh2=0.022, omch2=0.12, omk=0.0,
         :param tau: optical depth; if None, current Reion settings are not changed
         :param deltazrei: redshift width of reionization; if None, uses default
         :param Alens: (non-physical) scaling of the lensing potential compared to prediction
-        :param bbn_predictor: :class:`.bbn.BBNPredictor` instance used to get YHe from BBN consistency if YHe is None
+        :param bbn_predictor: :class:`.bbn.BBNPredictor` instance used to get YHe from BBN consistency if YHe is None,
+         or name of a BBN predictor class, or file name of an interpolation table
         :param theta_H0_range: if thetastar or cosmomc_theta is specified, the min, max interval of H0 values to map to;
           if H0 is outside this range it will raise an exception.
         """
 
         if YHe is None:
             # use BBN prediction
-            self.bbn_predictor = bbn_predictor or bbn.get_default_predictor()
+            if isinstance(bbn_predictor, six.string_types):
+                self.bbn_predictor = bbn.get_predictor(bbn_predictor)
+            else:
+                self.bbn_predictor = bbn_predictor or bbn.get_predictor()
             YHe = self.bbn_predictor.Y_He(ombh2 * (constants.COBE_CMBTemp / TCMB) ** 3, nnu - standard_neutrino_neff)
         self.YHe = YHe
         self.TCMB = TCMB
@@ -465,7 +470,7 @@ def set_cosmology(self, H0=None, ombh2=0.022, omch2=0.12, omk=0.0,
             omnuh2 = mnu / neutrino_mass_fac * (nnu / 3.0) ** 0.75
         omnuh2_sterile = meffsterile / neutrino_mass_fac
         if omnuh2_sterile > 0 and nnu < standard_neutrino_neff:
-            raise CAMBError('sterile neutrino mass required Neff> %.3g'%(constants.default_nnu))
+            raise CAMBError('sterile neutrino mass required Neff> %.3g' % (constants.default_nnu))
         if omnuh2 and not num_massive_neutrinos:
             raise CAMBError('non-zero mnu with zero num_massive_neutrinos')
 
@@ -474,7 +479,7 @@ def set_cosmology(self, H0=None, ombh2=0.022, omch2=0.12, omk=0.0,
         self.omk = omk
         if omnuh2_sterile > 0:
             if nnu < standard_neutrino_neff:
-                raise CAMBError('nnu < %.3g with massive sterile'%(constants.default_nnu))
+                raise CAMBError('nnu < %.3g with massive sterile' % (constants.default_nnu))
         assert num_massive_neutrinos == int(num_massive_neutrinos)
         self.f_SetNeutrinoHierarchy(byref(c_double(omnuh2)), byref(c_double(omnuh2_sterile)),
                                     byref(c_double(nnu)),