Python Class.compute方法代码示例

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [as 别名]
def test_class_setup():
    cosmology = astropy.cosmology.Planck13
    assert cosmology.Om0 == cosmology.Odm0 + cosmology.Ob0
    assert 1 == (cosmology.Om0 + cosmology.Ode0 + cosmology.Ok0 +
                 cosmology.Ogamma0 + cosmology.Onu0)
    class_parameters = get_class_parameters(cosmology)
    try:
        from classy import Class
        cosmo = Class()
        cosmo.set(class_parameters)
        cosmo.compute()
        assert cosmo.h() == cosmology.h
        assert cosmo.T_cmb() == cosmology.Tcmb0.value
        assert cosmo.Omega_b() == cosmology.Ob0
        # Calculate Omega(CDM)_0 two ways:
        assert abs((cosmo.Omega_m() - cosmo.Omega_b()) -
                   (cosmology.Odm0 - cosmology.Onu0)) < 1e-8
        assert abs(cosmo.Omega_m() - (cosmology.Om0 - cosmology.Onu0)) < 1e-8
        # CLASS calculates Omega_Lambda itself so this is a non-trivial test.
        calculated_Ode0 = cosmo.get_current_derived_parameters(
            ['Omega_Lambda'])['Omega_Lambda']
        assert abs(calculated_Ode0 - (cosmology.Ode0 + cosmology.Onu0)) < 1e-5
        cosmo.struct_cleanup()
        cosmo.empty()
    except ImportError:
        pass

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [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 compute [as 别名]
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()
        

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [as 别名]
    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

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [as 别名]
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

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [as 别名]
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']}

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [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 compute [as 别名]
                   '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
# returning some characteristic scales
# (of Hubble crossing, sound horizon crossing, etc.) at different time

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

# 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
print "Omega_Lambda =",derived['Omega0_lambda']

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

# import classy module
from classy import Class


# In[ ]:

# create instance of the class "Class"
LambdaCDM = Class()
# pass input parameters
LambdaCDM.set({'omega_b':0.022032,'omega_cdm':0.12038,'h':0.67556,'A_s':2.215e-9,'n_s':0.9619,'tau_reio':0.0925})
LambdaCDM.set({'output':'tCl,pCl,lCl,mPk','lensing':'yes','P_k_max_1/Mpc':3.0})
# run class
LambdaCDM.compute()


# In[ ]:

# get all C_l output
cls = LambdaCDM.lensed_cl(2500)
# To check the format of cls
cls.viewkeys()


# In[ ]:

ll = cls['ell'][2:]
clTT = cls['tt'][2:]
clEE = cls['ee'][2:]

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [as 别名]
class classy(SlikPlugin):
    """
    Plugin for CLASS.

    Credit: Brent Follin, Teresa Hamill
    """

    #{cosmoslik name : class name}
    name_mapping = {'As':'A_s',
                    'ns':'n_s',
                    'r':'r',
                    'nt':'n_t',
                    'ombh2':'omega_b',
                    'omch2':'omega_cdm',
                    'omnuh2':'omega_ncdm',
                    'tau':'tau_reio',
                    'H0':'H0',
                    'massive_neutrinos':'N_ncdm',
                    'massless_neutrinos':'N_ur',
                    'Yp':'YHe',
                    'pivot_scalar':'k_pivot'}


    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,
                 ombh2,
                 omch2,
                 H0,
                 As,
                 ns,
                 tau,
                 omnuh2, #0.006
                 w=None,
                 r=None,
                 nrun=None,
                 omk=0,
                 Yp=None,
                 Tcmb=2.7255,
                 massive_neutrinos=1,
                 massless_neutrinos=2.046,
                 l_max_scalar=3000,
                 l_max_tensor=3000,
                 pivot_scalar=0.002,
                 outputs=[],
                 **kwargs):


        
        self.model.set(output='tCl, lCl, pCl',
                       lensing='yes',
                       l_max_scalars=l_max_scalar,
                       **{self.name_mapping[k]:v for k,v in locals().items() 
                          if k in self.name_mapping and v is not None})
        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']}

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [as 别名]
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)
	params = {
	    'output': 'tCl lCl',

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [as 别名]
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..)

    """

    @classmethod
    def setUpClass(self):
        self.faulty_figs_path = os.path.join(
            os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-1]), "faulty_figs"
        )

        if os.path.isdir(self.faulty_figs_path):
            shutil.rmtree(self.faulty_figs_path)

        os.mkdir(self.faulty_figs_path)

    @classmethod
    def tearDownClass(self):
        pass

    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 tearDown(self):
        self.cosmo.struct_cleanup()
        self.cosmo.empty()
        self.cosmo_newt.struct_cleanup()
        self.cosmo_newt.empty()
        del self.scenario

    def poormansname(self, somedict):
        string = "_".join([k + "=" + str(v) for k, v in somedict.iteritems()])
        string = string.replace("/", "%")
        string = string.replace(",", "")
        string = string.replace(" ", "")
        return string

    @parameterized.expand(TUPLE_ARRAY)
    def test_0wrapper_implementation(self, inputdict):
        """Create a few instances based on different cosmologies"""
        self.scenario.update(inputdict)

        self.name = self.poormansname(inputdict)

        sys.stderr.write("\n\n---------------------------------\n")
        sys.stderr.write("| Test case %s |\n" % self.name)
        sys.stderr.write("---------------------------------\n")
        for key, value in self.scenario.iteritems():
            sys.stderr.write("%s = %s\n" % (key, value))
            sys.stdout.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_dict = {"tCl": ["tt"], "lCl": ["pp"], "pCl": ["ee", "bb"]}
        density_cl_list = ["nCl", "sCl"]

        # '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
#.........这里部分代码省略.........

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [as 别名]
class classy(SlikPlugin):
    """
    Plugin for CLASS.
    Credit: Brent Follin, Teresa Hamill, Andy Scacco
    """

    #{cosmoslik name : class name} - This needs to be done even for variables with the same name (because of for loop in self.model.set)!
    name_mapping = {'As':'A_s',
                    'ns':'n_s',
                    'r':'r',
                    'k_c':'k_c',
                    'alpha_exp':'alpha_exp',
                    'nt':'n_t',
                    'ombh2':'omega_b',
                    'omch2':'omega_cdm',
                    'omnuh2':'omega_ncdm',
                    'tau':'tau_reio',
                    'H0':'H0',
                    'massive_neutrinos':'N_ncdm',
                    'massless_neutrinos':'N_ur',
                    'Yp':'YHe',
                    'pivot_scalar':'k_pivot',
                    #'Tcmb':'T_cmb',
                    #'P_k_max_hinvMpc':'P_k_max_h/Mpc'
                    #'w':'w0_fld',
                    #'nrun':'alpha_s',
                    #'omk':'Omega_k',
                    #'l_max_scalar':'l_max_scalars',
                    #'l_max_tensor':'l_max_tensors'
                    }


    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,
                 ombh2,
                 omch2,
                 H0,
                 As,
                 ns,
                 k_c,
                 alpha_exp,
                 tau,
                 #omnuh2=0, #0.006  #None means that Class will take the default for this, maybe?
                 w=None,
                 r=None,
                 nrun=None,
                 omk=0,
                 Yp=None,
                 Tcmb=2.7255,
                 #massive_neutrinos=0,
                 massless_neutrinos=3.046,
                 l_max_scalar=3000,
                 l_max_tensor=3000,
                 pivot_scalar=0.05,
                 outputs=[],
                 **kwargs):


        
        self.model.set(output='tCl, lCl, pCl',
                       lensing='yes',
                       l_max_scalars=l_max_scalar,
                       **{self.name_mapping[k]:v for k,v in locals().items() 
                          if k in self.name_mapping and v is not None})
        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']}

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import compute [as 别名]
                  'background_verbose':1
                 }

# array of k values in 1/Mpc
kvec = np.logspace(-4,np.log10(3),100)
# array for storing legend
legarray = []

# loop over total mass values
for sum_masses in [0.1, 0.115, 0.13]:
    # 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()