Lot of things

work_in_class/compute_models.py

@@ -94,7 +94,7 @@ def add_derived(self, varied_name, keys, value):
 
     def compute_models(self, params, varied_name, index_variable, values,
                        back=[], thermo=[], prim=[], pert=[], trans=[],
-                       pk=[0.0001, 0.1, 100], extra=[], update=True,
+                       pk=[-4, -1, 100], extra=[], update=True,
                        cosmo_msg=False, texname=""):
         """
         Fill dic with the hi_class output structures for the model with given
@@ -157,7 +157,7 @@ def compute_models(self, params, varied_name, index_variable, values,
             except Exception, e:
                 print "Error: skipping {}={}".format(key, val)
                 if cosmo_msg:
-                    print e
+                    print str(e) + '\n'
 
                 continue
 
@@ -222,9 +222,8 @@ def compute_models(self, params, varied_name, index_variable, values,
             except CosmoSevereError:
                 pass
 
-
             if ("output" in self.cosmo.pars) and ('mPk' in self.cosmo.pars['output']):
-                k_array = np.linspace(*pk)
+                k_array = np.logspace(*pk)
                 pk_array = np.array([self.cosmo.pk(k, 0) for k in k_array])
 
                 d['pk'] = {'k': k_array, 'pk': pk_array}

work_in_class/cosmo_extra.py

@@ -0,0 +1,50 @@
+#!/usr/bin/python
+import numpy as np
+from classy import CosmoSevereError
+import scipy.integrate as integrate
+
+#####
+# NOTE: Emilio's functions from classy.pyx. self should not be used like this, but let's do an exception.
+#####
+
+def growthrate_at_k_and_z(cosmo, k, z):
+    dz = 0.005
+    f=-0.5*(1+z)*(-3.*np.log(cosmo.pk(k,z))+4.*np.log(cosmo.pk(k,z+dz))-np.log(cosmo.pk(k,z+2.*dz)))/(2.*dz)
+    if (f==0.):
+        raise CosmoSevereError(
+            "The step in redshift is too small to compute growth rate")
+    return f
+
+def growthrate_at_z(cosmo, z):
+    k_fid = 0.01
+    return growthrate_at_k_and_z(cosmo, k_fid, z)
+#########
+
+def fprime(OmegaDEw, OmegaM):
+    """
+    Return function df/dz(z, f).
+
+    OmegaDEw and OmegaM must be functions of z
+    """
+    def wrap(z, f):
+        try:
+            #output = 1.5 * OmegaM(z) - 0.5 * (1 - 3*w(z)*OmegaDE(z)) * f - f**2
+            output = ((0.5 * (1 - 3*OmegaDEw(z))) * f + f**2 - 1.5 * OmegaM(z)) / (1+z)
+        except Exception as e:
+            raise e
+        return output
+    return wrap
+def angular_distance(z, H, z_points=None):
+    """
+    Return the angular distance.
+    z, H must be CLASS-like arrays
+    """
+    DA_fit = []
+    if z_points is None:
+        z_points = z
+
+    for i in z_points:
+        #DA_fit.append(1/(1+i)*integrate.trapz(1/H_fit[zTMP<=i][::-1], zTMP[zTMP<=i][::-1]))
+        DA_fit.append(1/(1+i)*integrate.simps(1/H[z <= i][::-1], z[z <= i][::-1], even='last'))
+
+    return np.array(DA_fit)

work_in_class/data_mp_likelihoods.py

@@ -194,7 +194,7 @@ def remove_experiments(self, experiments):
 
         for exp in self.__name_str_to_list(experiments):
             if exp in self.experiments:
-                self.experiments.remove(experiment)
+                self.experiments.remove(exp)
 
     def get_mcmc_parameters(self, table_of_strings):
         """
@@ -298,8 +298,9 @@ def clear(self, params=''):
             self.mcmc_parameters.clear()
 
     def compute_lkl(self, cosmo, params, experiments=[], nuisance=[], overwrite=False, cosmo_struct_cleanup=False):
-        """Compute the likelihood for the parms given. Params will update the
-        cosmo_arguments and nuisance_parameters unless overwrite is True.
+        """Compute the likelihood for the parms and experiments given. Params
+        will update the cosmo_arguments and nuisance_parameters unless overwrite
+        is True. If no experment is passed, it will used the initialised ones.
 
         For the time being, params must be given in classy form.
 
@@ -315,24 +316,27 @@ def compute_lkl(self, cosmo, params, experiments=[], nuisance=[], overwrite=Fals
                     them could have been set by the experiments so that \
                     their removal will likely cause an error."
             self.clear()
+        else:
+            self.clear("cosmo")
 
         if not experiments:
             experiments = self.experiments
-
-        lkl = 0
-        for experiment in self.__name_str_to_list(experiments):
-            if experiment not in self.experiments:
-                print "Adding {} to list of initialised experiments".format(experiment)
-                self.add_experiments([experiment], nuisance)
+        else:
+            for experiment in self.__name_str_to_list(experiments):
+                if experiment not in self.experiments:
+                    print "Adding {} to list of initialised experiments".format(experiment)
+                    self.add_experiments([experiment], nuisance)
 
         # set and compute must be called after initializing the lkl's. Some of them add
         # some mandatory parameters (e.g. {'output': 'mPk'}).
 
         self.update_params(params)
+        cosmo.empty()
         cosmo.set(self.__params_to_cosmo_classy()) # TODO: Modify to accept params as in MP input
         cosmo.compute()
 
-        for experiment in self.experiments:
+        lkl = 0
+        for experiment in experiments:
             lkl += self.lkl[experiment].loglkl(cosmo, self)
 
         if cosmo_struct_cleanup: