Python Class.pk方法代码示例

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import pk [as 别名]
def m_Pk(k=np.logspace(-3, 0., 100), z=0.53, nl_model='trg'):
    print k
    cosmo = Class()

    CLASS_INPUT = {}

    CLASS_INPUT['Mnu'] = ([{'N_eff': 0.0, 'N_ncdm': 1, 'm_ncdm': 0.06, 'deg_ncdm': 3.0}], 'normal')
    CLASS_INPUT['Output_spectra'] = ([{'output': 'mPk', 'P_k_max_1/Mpc': 1, 'z_pk': z}], 'power')

    CLASS_INPUT['Nonlinear'] = ([{'non linear': nl_model}], 'power')
            
    verbose = {}
    #    'input_verbose': 1,
    #    'background_verbose': 1,
    #    'thermodynamics_verbose': 1,
    #    'perturbations_verbose': 1,
    #    'transfer_verbose': 1,
    #    'primordial_verbose': 1,
    #    'spectra_verbose': 1,
    #    'nonlinear_verbose': 1,
    #    'lensing_verbose': 1,
    #    'output_verbose': 1
    #    }

    cosmo.struct_cleanup()
    cosmo.empty()


    INPUTPOWER = []
    INPUTNORMAL = [{}]
    for key, value in CLASS_INPUT.iteritems():
        models, state = value
        if state == 'power':
            INPUTPOWER.append([{}]+models)
        else:
            INPUTNORMAL.extend(models)

        PRODPOWER = list(itertools.product(*INPUTPOWER))

        DICTARRAY = []
        for normelem in INPUTNORMAL:
            for powelem in PRODPOWER:  # itertools.product(*modpower):
                temp_dict = normelem.copy()
                for elem in powelem:
                    temp_dict.update(elem)
                DICTARRAY.append(temp_dict)

    scenario = {}
    for dic in DICTARRAY:
        scenario.update(dic)
    setting = cosmo.set(dict(verbose.items()+scenario.items()))
    cosmo.compute()
    pk_out = [] 
    for k_i in k: 
        pk_out.append(cosmo.pk(k_i,z))
    return pk_out 

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import pk [as 别名]
def calculate_power(cosmology, k_min, k_max, z=0, num_k=500, scaled_by_h=True,
                    n_s=0.9619, logA=3.0980):
    """
    Calculate the power spectrum P(k,z) over the range k_min <= k <= k_max.
    """
    try:
        from classy import Class
        cosmo = Class()
    except ImportError:
        raise RuntimeError('power.calculate_power requires classy.')

    class_parameters = get_class_parameters(cosmology)
    class_parameters['output'] = 'mPk'
    if scaled_by_h:
        class_parameters['P_k_max_h/Mpc'] = k_max
    else:
        class_parameters['P_k_max_1/Mpc'] = k_max
    class_parameters['n_s'] = n_s
    class_parameters['ln10^{10}A_s'] = logA
    cosmo.set(class_parameters)
    cosmo.compute()

    if scaled_by_h:
        k_scale = cosmo.h()
        Pk_scale = cosmo.h()**3
    else:
        k_scale = 1.
        Pk_scale = 1.

    result = np.empty((num_k,), dtype=[('k', float), ('Pk', float)])
    result['k'][:] = np.logspace(np.log10(k_min), np.log10(k_max), num_k)
    for i, k in enumerate(result['k']):
        result['Pk'][i] = cosmo.pk(k * k_scale, z) * Pk_scale

    cosmo.struct_cleanup()
    cosmo.empty()

    return result

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import pk [as 别名]

#.........这里部分代码省略.........
        self.assertTrue(setting, "Class failed to initialize with input dict")

        cl_list = ["tCl", "lCl", "pCl", "nCl", "sCl"]

        # Depending on the cases, the compute should fail or not
        should_fail = True
        output = self.scenario["output"].split()
        for elem in output:
            if elem in ["tCl", "pCl"]:
                for elem2 in output:
                    if elem2 == "lCl":
                        should_fail = False
                        break

        if not should_fail:
            self.cosmo.compute()
        else:
            self.assertRaises(CosmoSevereError, self.cosmo.compute)
            return

        self.assertTrue(self.cosmo.state, "Class failed to go through all __init__ methods")
        if self.cosmo.state:
            print "--> Class is ready"
        # Depending
        if "output" in self.scenario.keys():
            # Positive tests
            output = self.scenario["output"]
            for elem in output.split():
                if elem in cl_list:
                    print "--> testing raw_cl function"
                    cl = self.cosmo.raw_cl(100)
                    self.assertIsNotNone(cl, "raw_cl returned nothing")
                    self.assertEqual(np.shape(cl["tt"])[0], 101, "raw_cl returned wrong size")
                if elem == "mPk":
                    print "--> testing pk function"
                    pk = self.cosmo.pk(0.1, 0)
                    self.assertIsNotNone(pk, "pk returned nothing")
            # Negative tests of output functions
            if not any([elem in cl_list for elem in output.split()]):
                print "--> testing absence of any Cl"
                self.assertRaises(CosmoSevereError, self.cosmo.raw_cl, 100)
            if "mPk" not in self.scenario["output"].split():
                print "--> testing absence of mPk"
                # args = (0.1, 0)
                self.assertRaises(CosmoSevereError, self.cosmo.pk, 0.1, 0)

    @parameterized.expand(
        itertools.product(("massless", "massive", "both"), ("photons", "massless", "exact"), ("t", "s, t"))
    )
    def test_tensors(self, scenario, method, modes):
        """Test the new tensor mode implementation"""
        self.scenario = {}
        if scenario == "massless":
            self.scenario.update({"N_eff": 3.046, "N_ncdm": 0})
        elif scenario == "massiv":
            self.scenario.update({"N_eff": 0, "N_ncdm": 2, "m_ncdm": "0.03, 0.04", "deg_ncdm": "2, 1"})
        elif scenario == "both":
            self.scenario.update({"N_eff": 1.5, "N_ncdm": 2, "m_ncdm": "0.03, 0.04", "deg_ncdm": "1, 0.5"})

        sys.stderr.write("\n\n---------------------------------\n")
        sys.stderr.write("| Test case: %s %s %s |\n" % (scenario, method, modes))
        sys.stderr.write("---------------------------------\n")
        self.scenario.update({"tensor method": method, "modes": modes, "output": "tCl, pCl"})
        for key, value in self.scenario.iteritems():
            sys.stderr.write("%s = %s\n" % (key, value))
        sys.stderr.write("\n")
        self.cosmo.set(dict(self.verbose.items() + self.scenario.items()))
        self.cosmo.compute()

    @parameterized.expand(itertools.izip(powerset(["100*theta_s", "Omega_dcdmdr"]), powerset([1.04, 0.20])))
    def test_shooting_method(self, variables, values):
        Omega_cdm = 0.25

        scenario = {"Omega_b": 0.05}

        for variable, value in zip(variables, values):
            scenario.update({variable: value})

        if "Omega_dcdmdr" in variables:
            scenario.update({"Gamma_dcdm": 100, "Omega_cdm": Omega_cdm - scenario["Omega_dcdmdr"]})
        else:
            scenario.update({"Omega_cdm": Omega_cdm})

        sys.stderr.write("\n\n---------------------------------\n")
        sys.stderr.write("| Test shooting: %s |\n" % (", ".join(variables)))
        sys.stderr.write("---------------------------------\n")
        for key, value in scenario.iteritems():
            sys.stderr.write("%s = %s\n" % (key, value))
        sys.stderr.write("\n")

        scenario.update(self.verbose)
        self.assertTrue(self.cosmo.set(scenario), "Class failed to initialise with this input")
        self.assertRaises
        self.cosmo.compute()

        # Now, check that the values are properly extracted
        for variable, value in zip(variables, values):
            if variable == "100*theta_s":
                computed_value = self.cosmo.get_current_derived_parameters([variable])[variable]
                self.assertAlmostEqual(value, computed_value, places=5)

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import pk [as 别名]

# In[ ]:

plt.savefig('warmup_cltt.pdf')


# In[ ]:

# get P(k) at redhsift z=0
import numpy as np
kk = np.logspace(-4,np.log10(3),1000) # k in h/Mpc
Pk = [] # P(k) in (Mpc/h)**3
h = LambdaCDM.h() # get reduced Hubble for conversions to 1/Mpc
for k in kk:
    Pk.append(LambdaCDM.pk(k*h,0.)*h**3) # function .pk(k,z)


# In[ ]:

# plot P(k)
plt.figure(2)
plt.xscale('log');plt.yscale('log');plt.xlim(kk[0],kk[-1])
plt.xlabel(r'$k \,\,\,\, [h/\mathrm{Mpc}]$')
plt.ylabel(r'$P(k) \,\,\,\, [\mathrm{Mpc}/h]^3$')
plt.plot(kk,Pk,'b-')


# In[ ]:

plt.savefig('warmup_pk.pdf')

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import pk [as 别名]

#.........这里部分代码省略.........
            self.scenario.update({'Omega_k': 0.01})
        elif name == 'Negative_Omega_k':
            self.scenario.update({'Omega_k': -0.01})
        elif name == 'Isocurvature_modes':
            self.scenario.update({'ic': 'ad,nid,cdi', 'c_ad_cdi': -0.5})

        self.scenario.update(scenario)
        if scenario != {}:
            self.scenario.update(gauge)
        self.scenario.update(nonlinear)

        sys.stderr.write('\n\n---------------------------------\n')
        sys.stderr.write('| Test case %s |\n' % name)
        sys.stderr.write('---------------------------------\n')
        for key, value in self.scenario.iteritems():
            sys.stderr.write("%s = %s\n" % (key, value))
        sys.stderr.write("\n")

        setting = self.cosmo.set(
            dict(self.verbose.items()+self.scenario.items()))
        self.assertTrue(setting, "Class failed to initialize with input dict")

        cl_list = ['tCl', 'lCl', 'pCl', 'nCl', 'sCl']

        # Depending on the cases, the compute should fail or not
        should_fail = True
        output = self.scenario['output'].split()
        for elem in output:
            if elem in ['tCl', 'pCl']:
                for elem2 in output:
                    if elem2 == 'lCl':
                        should_fail = False
                        break

        if not should_fail:
            self.cosmo.compute()
        else:
            self.assertRaises(CosmoSevereError, self.cosmo.compute)
            return

        self.assertTrue(
            self.cosmo.state,
            "Class failed to go through all __init__ methods")
        if self.cosmo.state:
            print '--> Class is ready'
        # Depending
        if 'output' in self.scenario.keys():
            # Positive tests
            output = self.scenario['output']
            for elem in output.split():
                if elem in cl_list:
                    print '--> testing raw_cl function'
                    cl = self.cosmo.raw_cl(100)
                    self.assertIsNotNone(cl, "raw_cl returned nothing")
                    self.assertEqual(
                        np.shape(cl['tt'])[0], 101,
                        "raw_cl returned wrong size")
                if elem == 'mPk':
                    print '--> testing pk function'
                    pk = self.cosmo.pk(0.1, 0)
                    self.assertIsNotNone(pk, "pk returned nothing")
            # Negative tests of output functions
            if not any([elem in cl_list for elem in output.split()]):
                print '--> testing absence of any Cl'
                self.assertRaises(CosmoSevereError, self.cosmo.raw_cl, 100)
            if 'mPk' not in self.scenario['output'].split():
                print '--> testing absence of mPk'
                #args = (0.1, 0)
                self.assertRaises(CosmoSevereError, self.cosmo.pk, 0.1, 0)

        print '~~~~~~~~ passed ? '

    @parameterized.expand(
        itertools.product(
            ('massless', 'massive', 'both'),
            ('photons', 'massless', 'exact'),
            ('t', 's, t')))
    def test_tensors(self, scenario, method, modes):
        """Test the new tensor mode implementation"""
        self.scenario = {}
        if scenario == 'massless':
            self.scenario.update({'N_eff': 3.046, 'N_ncdm':0})
        elif scenario == 'massiv':
            self.scenario.update(
                {'N_eff': 0, 'N_ncdm': 2, 'm_ncdm': '0.03, 0.04',
                 'deg_ncdm': '2, 1'})
        elif scenario == 'both':
            self.scenario.update(
                {'N_eff': 1.5, 'N_ncdm': 2, 'm_ncdm': '0.03, 0.04',
                 'deg_ncdm': '1, 0.5'})

        self.scenario.update({
            'tensor method': method, 'modes': modes,
            'output': 'tCl, pCl'})
        for key, value in self.scenario.iteritems():
            sys.stderr.write("%s = %s\n" % (key, value))
        sys.stderr.write("\n")
        self.cosmo.set(
            dict(self.verbose.items()+self.scenario.items()))
        self.cosmo.compute()

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import pk [as 别名]

#.........这里部分代码省略.........
        # 'lensing' is always set to yes. Therefore, trying to compute 'tCl' or
        # 'pCl' will fail except if we also ask for 'lCl'. The flag
        # 'should_fail' stores this status.
        sys.stderr.write("Should")
        should_fail = self.test_incompatible_input()
        if should_fail:
            sys.stderr.write(" fail...\n")
        else:
            sys.stderr.write(" not fail...\n")

        if not should_fail:
            self.cosmo.compute()
        else:
            self.assertRaises(CosmoSevereError, self.cosmo.compute)
            return

        self.assertTrue(self.cosmo.state, "Class failed to go through all __init__ methods")
        if self.cosmo.state:
            print "--> Class is ready"
        # Depending
        if "output" in self.scenario.keys():
            # Positive tests of raw cls
            output = self.scenario["output"]
            for elem in output.split():
                if elem in cl_dict.keys():
                    for cl_type in cl_dict[elem]:
                        sys.stderr.write("--> testing raw_cl for %s\n" % cl_type)
                        cl = self.cosmo.raw_cl(100)
                        self.assertIsNotNone(cl, "raw_cl returned nothing")
                        self.assertEqual(np.shape(cl[cl_type])[0], 101, "raw_cl returned wrong size")
                    # TODO do the same for lensed if 'lCl' is there, and for
                    # density cl
                if elem == "mPk":
                    sys.stderr.write("--> testing pk function\n")
                    pk = self.cosmo.pk(0.1, 0)
                    self.assertIsNotNone(pk, "pk returned nothing")
            # Negative tests of output functions
            if not any([elem in cl_dict.keys() for elem in output.split()]):
                sys.stderr.write("--> testing absence of any Cl\n")
                self.assertRaises(CosmoSevereError, self.cosmo.raw_cl, 100)
            if "mPk" not in output.split():
                sys.stderr.write("--> testing absence of mPk\n")
                self.assertRaises(CosmoSevereError, self.cosmo.pk, 0.1, 0)

        if COMPARE_OUTPUT:
            # Now, compute with Newtonian gauge, and compare the results
            self.cosmo_newt.set(dict(self.verbose.items() + self.scenario.items()))
            self.cosmo_newt.set({"gauge": "newtonian"})
            self.cosmo_newt.compute()
            # Check that the computation worked
            self.assertTrue(self.cosmo_newt.state, "Class failed to go through all __init__ methods in Newtonian gauge")

            self.compare_output(self.cosmo, self.cosmo_newt)

    def test_incompatible_input(self):

        should_fail = False

        # If we have tensor modes, we must have one tensor observable,
        # either tCl or pCl.
        if has_tensor(self.scenario):
            if "output" not in self.scenario.keys():
                should_fail = True
            else:
                output = self.scenario["output"].split()
                if "tCl" not in output and "pCl" not in output:

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import pk [as 别名]
    # normal hierarchy
    [m1, m2, m3] = get_masses(2.45e-3,7.50e-5, sum_masses, 'NH')
    NH = Class()
    NH.set(commonsettings)
    NH.set({'m_ncdm':str(m1)+','+str(m2)+','+str(m3)})
    NH.compute()
    # inverted hierarchy
    [m1, m2, m3] = get_masses(2.45e-3,7.50e-5, sum_masses, 'IH')
    IH = Class()
    IH.set(commonsettings)
    IH.set({'m_ncdm':str(m1)+','+str(m2)+','+str(m3)})
    IH.compute()
    pkNH = []
    pkIH = []
    for k in kvec:
        pkNH.append(NH.pk(k,0.))
        pkIH.append(IH.pk(k,0.))
    NH.struct_cleanup()
    IH.struct_cleanup()
    # extract h value to convert k from 1/Mpc to h/Mpc
    h = NH.h()
    plt.semilogx(kvec/h,1-np.array(pkNH)/np.array(pkIH))
    legarray.append(r'$\Sigma m_i = '+str(sum_masses)+'$eV')

plt.axhline(0,color='k')
plt.xlim(kvec[0]/h,kvec[-1]/h)
plt.xlabel(r'$k [h \mathrm{Mpc}^{-1}]$')
plt.ylabel(r'$1-P(k)^\mathrm{NH}/P(k)^\mathrm{IH}$')
plt.legend(legarray)

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import pk [as 别名]
class Model():
    def __init__(self, cosmo=None):
        """
        Initialize the Model class. By default Model uses its own Class
        instance.

        cosmo = external Class instance. Default is None
        """
        if cosmo:
            self.cosmo = cosmo
        else:
            self.cosmo = Class()
        self.computed = {}
        self.texnames = {}

    def __set_scale(self, axes, xscale, yscale):
        """
        Set scales for axes in axes array.

        axes = axes array (e.g. f, ax = plt.subplots(2,2))
        xscale = linear array of xscale.
        yscale = linear array of yscale.

        Scales are set once axes is flatten. Each plot is counted from left to
        right an from top to bottom.
        """
        for i, ax in enumerate(axes.flat):
            ax.set_xscale(xscale[i])
            ax.set_yscale(yscale[i])

    def __set_label(self, axes, xlabel, ylabel):
        """
        Set labels for axes in axes array.

        axes = axes array (e.g. f, ax = plt.subplots(2,2))
        xlabel = linear array of xlabels.
        ylabel = linear array of ylabels.

        Labels are set once axes is flatten. Each plot is counted from left to
        right an from top to bottom.
        """
        for i, ax in enumerate(axes.flat):
            ax.set_xlabel(xlabel[i])
            ax.set_ylabel(ylabel[i])

    def __store_cl(self, cl_dic):
        """
        Store cl's as (l*(l+1)/2pi)*cl, which is much more useful.
        """

        ell = cl_dic['ell'][2:]

        for cl, list_val in cl_dic.iteritems():
            list_val = list_val[2:]
            if (list_val == ell).all():
                cl_dic[cl] = list_val
                continue
            list_val = (ell * (ell + 1) / (2 * np.pi)) * list_val
            cl_dic[cl] = list_val  # Remove first two null items (l=0,1)

        return cl_dic

    def add_derived(self, varied_name, keys, value):
        """
        Add a derived parameter for varied_name dictionary.

        varied_name = varied variable's name.
        keys = list of keys in descending level.
        value = value to store for new dictionary key.
        """

        dic = self.computed[varied_name]

        for key in keys:
            if key not in dic:
                dic[key] = {}

            dic = dic[key]

        dic.update(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,
                       cosmo_msg=False, texname=""):
        """
        Fill dic with the hi_class output structures for the model with given
        params, modifying the varied_name value with values.

        params = parameters to be set in Class. They must be in agreement with
                what is asked for.
        varied_name = the name of the variable you are modifying. It will be
                      used as key in dic assigned to its background structures.
        index_variable = variable's index in parameters_smg array.
        values = varied variable values you want to compute the cosmology for.
        back = list of variables to store from background. If 'all', store the
              whole dictionary.
        thermo = list of variables to store from thermodynamics. If 'all',
                  store the whole dictionary.
        prim = list of variables to store from primordial. If 'all', store the
#.........这里部分代码省略.........