Python classy.Class類代碼示例

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 

 def setup(self):
     """
     Create an instance of Class and attach it to self.
     """
     self.cosmo = Class()
     self.cosmo.set(self.constants)
     self.cosmo.compute()
     self.cosmo.struct_cleanup()
     

    def __init__(self):
        super(classy,self).__init__()

        try:
            from classy import Class
        except ImportError:
            raise Exception("Failed to import CLASS python wrapper 'Classy'.")

        self.model = Class()

    def setUp(self):
        """
        set up data used in the tests.
        setUp is called before each test function execution.
        """
        self.cosmo = Class()
        self.cosmo_newt = Class()

        self.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}
        self.scenario = {}

    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 setUp(self):
        """
        set up data used in the tests.
        setUp is called before each test function execution.
        """
        self.cosmo = Class()

        self.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,
        }
        self.scenario = {"lensing": "yes"}

class ClassCoreModule(object):
    
    def __init__(self, mapping=DEFAULT_PARAM_MAPPING, constants=CLASS_DEFAULT_PARAMS):
        """
        Core Module for the delegation of the computation of the cmb power
        spectrum to the Class wrapper classy.
        The defaults are for the 6 LambdaCDM cosmological parameters.
        
        :param mapping: (optional) dict mapping name of the parameter to the index
        :param constants: (optional) dict with constants overwriting CLASS defaults
        """
        self.mapping = mapping
        if constants is None:
            constants = {}
        self.constants = constants
        
    def __call__(self, ctx):
        p1 = ctx.getParams()
        
        params = self.constants.copy()
        for k,v in self.mapping.items():
            params[k] = p1[v]
        self.cosmo.set(params)
        self.cosmo.compute()
        if self.constants['lensing'] == 'yes':
            cls = self.cosmo.lensed_cl()
        else:
            cls = self.cosmo.raw_cl()
        Tcmb = self.cosmo.T_cmb()*1e6
        frac = Tcmb**2 * cls['ell'][2:] * (cls['ell'][2:] + 1) / 2. / pi
        ctx.add(CL_TT_KEY, frac*cls['tt'][2:])
        ctx.add(CL_TE_KEY, frac*cls['te'][2:])
        ctx.add(CL_EE_KEY, frac*cls['ee'][2:])
        ctx.add(CL_BB_KEY, frac*cls['bb'][2:])
        self.cosmo.struct_cleanup()

    def setup(self):
        """
        Create an instance of Class and attach it to self.
        """
        self.cosmo = Class()
        self.cosmo.set(self.constants)
        self.cosmo.compute()
        self.cosmo.struct_cleanup()
        

    def loglkl(self, params):
        cosmo = Class()
        cosmo.set(params)
        cosmo.compute()

        chi2 = 0.

        # for each point, compute angular distance da, radial distance dr,
        # volume distance dv, sound horizon at baryon drag rs_d,
        # theoretical prediction and chi2 contribution
        for i in range(self.num_points):

            da = cosmo.angular_distance(self.z[i])
            dr = self.z[i] / cosmo.Hubble(self.z[i])
            dv = pow(da * da * (1 + self.z[i]) * (1 + self.z[i]) * dr, 1. / 3.)
            rs = cosmo.rs_drag()

            if self.type[i] == 3:
                theo = dv / rs

            elif self.type[i] == 4:
                theo = dv

            elif self.type[i] == 5:
                theo = da / rs

            elif self.type[i] == 6:
                theo = 1. / cosmo.Hubble(self.z[i]) / rs

            elif self.type[i] == 7:
                theo = rs / dv

            chi2 += ((theo - self.data[i]) / self.error[i]) ** 2

        # return ln(L)
        # lkl = - 0.5 * chi2
        # return -2ln(L)
        lkl = chi2

        return lkl

def ComputeTransferData(settings, redshift):
    database_key = settings.copy()
    database_key.update({'redshift': tuple(redshift)})

    database = Database.Database(config.DATABASE_DIR)
    if database_key in database:
        return database[database_key], redshift
    else:
        cosmo = Class()
        cosmo.set(settings)
        cosmo.compute()

        outputData = [cosmo.get_transfer(z) for z in redshift]
        # Calculate d_g/4+psi
        for transfer_function_dict in outputData:
            transfer_function_dict["d_g/4 + psi"] = transfer_function_dict["d_g"]/4 + transfer_function_dict["psi"]
        # Now filter the relevant fields
        fields = TRANSFER_QUANTITIES + ["k (h/Mpc)"]
        outputData = [{field: outputData[i][field] for field in fields} for i in range(len(redshift))]

        database[database_key] = outputData
        return outputData, redshift

params = {
	    'output': 'tCl lCl',
	    'l_max_scalars': 2508,
	    'lensing': 'yes',
	    'P_k_ini type': 'external_Pk',
	    'command': 'python /home/andrew/Research/tools/class_public-2.4.3/external_Pk/generate_Pk_cosines.py',
	    'custom1': 0,
	    'custom2': 0,
	    'custom3': 0,
	    'custom4': 0,
	    'custom5': 0}

#Get the unperturbed cls for comparison
cosmo = Class()
cosmo.set(params)
cosmo.compute()
clso=cosmo.lensed_cl(2508)['tt'][30:]
ell = cosmo.lensed_cl(2508)['ell'][30:]

for i in range(len(clso)):
	clso[i]=ell[i]*(ell[i]+1)/(4*np.pi)*((2.726e6)**2)*clso[i]
a=np.zeros(5)
cosmo.struct_cleanup()
cosmo.empty()
dcls=np.zeros([clso.shape[0],5])
h=1e-6
for m in range(5):
	a[m]=h
	# Define your cosmology (what is not specified will be set to CLASS default parameters)

                   'recfast_z_initial':z_max_pk,
                   #'k_step_sub':'0.01',
                   'k_per_decade_for_pk':k_per_decade,
                   'k_per_decade_for_bao':k_per_decade,
                   'k_min_tau0':k_min_tau0, # this value controls the minimum k value in the figure
                   'perturb_sampling_stepsize':'0.05',
                   'P_k_max_1/Mpc':P_k_max_inv_Mpc,
                   'compute damping scale':'yes', # needed to output and plot Silk damping scale
                   'gauge':'newtonian'}

###############
#
# call CLASS
#
###############
M = Class()
M.set(common_settings)
M.compute()
#
# define conformal time sampling array
#
times = M.get_current_derived_parameters(['tau_rec','conformal_age'])
tau_rec=times['tau_rec']
tau_0 = times['conformal_age']
tau1 = np.logspace(math.log10(tau_ini),math.log10(tau_rec),tau_num_early)
tau2 = np.logspace(math.log10(tau_rec),math.log10(tau_0),tau_num_late)[1:]
tau2[-1] *= 0.999 # this tiny shift avoids interpolation errors
tau = np.concatenate((tau1,tau2))
tau_num = len(tau)
#
# use table of background and thermodynamics quantitites to define some functions

class TestClass(unittest.TestCase):
    """
    Testing Class and its wrapper classy on different cosmologies

    To run it, do
    ~] nosetest test_class.py

    It will run many times Class, on different cosmological scenarios, and
    everytime testing for different output possibilities (none asked, only mPk,
    etc..)

    """

    def setUp(self):
        """
        set up data used in the tests.
        setUp is called before each test function execution.
        """
        self.cosmo = Class()

        self.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,
        }
        self.scenario = {"lensing": "yes"}

    def tearDown(self):
        self.cosmo.struct_cleanup()
        self.cosmo.empty()
        del self.scenario

    @parameterized.expand(
        itertools.product(
            ("LCDM", "Mnu", "Positive_Omega_k", "Negative_Omega_k", "Isocurvature_modes"),
            (
                {"output": ""},
                {"output": "mPk"},
                {"output": "tCl"},
                {"output": "tCl pCl lCl"},
                {"output": "mPk tCl lCl", "P_k_max_h/Mpc": 10},
                {"output": "nCl sCl"},
                {"output": "tCl pCl lCl nCl sCl"},
            ),
            ({"gauge": "newtonian"}, {"gauge": "sync"}),
            ({}, {"non linear": "halofit"}),
        )
    )
    def test_wrapper_implementation(self, name, scenario, gauge, nonlinear):
        """Create a few instances based on different cosmologies"""
        if name == "Mnu":
            self.scenario.update({"N_ncdm": 1, "m_ncdm": 0.06})
        elif name == "Positive_Omega_k":
            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")
#.........這裏部分代碼省略.........

import numpy as np
from classy import Class


# In[ ]:

font = {'size'   : 20, 'family':'STIXGeneral'}
axislabelfontsize='large'
matplotlib.rc('font', **font)
matplotlib.mathtext.rcParams['legend.fontsize']='medium'


# In[ ]:

#Lambda CDM
LCDM = Class()
LCDM.set({'Omega_cdm':0.25,'Omega_b':0.05})
LCDM.compute()


# In[ ]:

#Einstein-de Sitter
CDM = Class()
CDM.set({'Omega_cdm':0.95,'Omega_b':0.05})
CDM.compute()

# Just to cross-check that Omega_Lambda is negligible
# (but not exactly zero because we neglected radiation)
derived = CDM.get_current_derived_parameters(['Omega0_lambda'])
print derived

class classy(SlikPlugin):
    """
    Plugin for CLASS.
    Credit: Brent Follin, Teresa Hamill, Andy Scacco
    """
    
    def __init__(self):
        super(classy,self).__init__()

        try:
            from classy import Class
        except ImportError:
            raise Exception("Failed to import CLASS python wrapper 'Classy'.")

        self.model = Class()


    def __call__(self,
        **kwargs):


        self.model.set(**kwargs)
        self.model.compute()

        ell = arange(l_max_scalar+1)
        self.cmb_result = {'cl_%s'%x:(self.model.lensed_cl(l_max_scalar)[x.lower()])*Tcmb**2*1e12*ell*(ell+1)/2/pi
                           for x in ['TT','TE','EE','BB','PP','TP']}

        self.model.struct_cleanup()
        self.model.empty()
        
        return self.cmb_result

    def get_bao_observables(self, z):
        return {'H':self.model.Hubble(z),
                'D_A':self.model.angular_distance(z),
                'c':1.0,
                'r_d':(self.model.get_current_derived_parameters(['rs_rec']))['rs_rec']}

                   # LambdaCDM parameters
                   'h':0.67556,
                   'omega_b':0.022032,
                   'omega_cdm':0.12038,
                   'A_s':2.215e-9,
                   'n_s':0.9619,
                   'tau_reio':0.0925,
                   # Take fixed value for primordial Helium (instead of automatic BBN adjustment)
                   'YHe':0.246,
                   # other output and precision parameters
                   'l_max_scalars':5000}
###############
#
# call CLASS
#
M = Class()
M.set(common_settings)
M.compute()
cl_tot = M.raw_cl(3000)
cl_lensed = M.lensed_cl(3000)
M.struct_cleanup()  # clean output
M.empty()           # clean input
#
M.set(common_settings) # new input
M.set({'temperature contributions':'tsw'})
M.compute()
cl_tsw = M.raw_cl(3000)
M.struct_cleanup()
M.empty()
#
M.set(common_settings)