import sys,os
import re
if sys.version_info.major>=3:
import _pickle as pkl
from io import StringIO
else:
import cPickle as pkl
from StringIO import StringIO
import numpy as np
from scipy import interpolate
import astropy.units as u
from astropy.cosmology import FLRW
#Option parsing method
from .settings import select_parser,LTSettings
#Parse CAMB output log
[docs]def parseLog(fname):
"""
Parse CAMB output log
:param fname: file name or file descriptor
:type fname: str. or file.
:returns: parsed log
:rtype: dict.
"""
#Get the filehandle
if type(fname)==file:
fp = fname
else:
fp = open(fname,"r")
#Dictionary with parsed log
parsed = dict()
parsed["sigma8"] = dict()
#Cycle over the lines in the log
for line in fp.readlines():
#w0/wa
match = re.match(r"\(w0, wa\) = \(([-\.0-9]+),[\s]+([-\.0-9]+)\)",line)
if match:
parsed["w0"],parsed["wa"] = [ float(v) for v in match.groups() ]
continue
#Parameters
match = re.match(r"Reion redshift[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["z_ion"] = float(match.groups()[0])
match = re.match(r"Om_b h\^2[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["Obh2"] = float(match.groups()[0])
match = re.match(r"Om_c h\^2[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["Omch2"] = float(match.groups()[0])
match = re.match(r"Om_nu h\^2[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["Onuh2"] = float(match.groups()[0])
match = re.match(r"Om_Lambda[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["Ode"] = float(match.groups()[0])
match = re.match(r"Om_K[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["Omk"] = float(match.groups()[0])
match = re.match(r"Om_m \(1-Om_K-Om_L\)[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["Om"] = float(match.groups()[0])
match = re.match(r"100 theta \(CosmoMC\)[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["100thetaMC"] = float(match.groups()[0])
match = re.match(r"Reion opt depth[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["tau_ion"] = float(match.groups()[0])
match = re.match(r"Age of universe\/GYr[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["Age"] = float(match.groups()[0]) * u.Gyr
match = re.match(r"zstar[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["zstar"] = float(match.groups()[0])
match = re.match(r"r_s\(zstar\)/Mpc[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["rs"] = float(match.groups()[0]) * u.Mpc
match = re.match(r"zdrag[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["zdrag"] = float(match.groups()[0])
match = re.match(r"r_s\(zdrag\)/Mpc[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["rs(zdrag)"] = float(match.groups()[0]) * u.Mpc
match = re.match(r"k_D\(zstar\) Mpc[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["kD(zstar)"] = float(match.groups()[0]) / u.Mpc
match = re.match(r"100\*theta_D[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["100thetaD"] = float(match.groups()[0])
match = re.match(r"z_EQ \(if v_nu=1\)[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["zEQ"] = float(match.groups()[0])
match = re.match(r"100\*theta_EQ[\s]+=[\s]+([0-9\.]+)",line)
if match:
parsed["100thetaEQ"] = float(match.groups()[0])
match = re.match(r"tau_recomb/Mpc[\s]+=[\s]+([0-9\.]+)[\s]+tau_now/Mpc =[\s]+([0-9\.]+)",line)
if match:
parsed["tau_rec"],parsed["tau_now"] = [float(v)*u.Mpc for v in match.groups()]
match = re.match(r"[\s]+at z =[\s]+([0-9E\-\+\.]+)[\s]+sigma8 \(all matter\)=[\s]+([0-9\.]+)",line)
if match:
z,sigma8 = [ float(v) for v in match.groups() ]
parsed["sigma8"][z] = sigma8
#Return
if type(fname)!=file:
fp.close()
return parsed
##################################################################################################
[docs]class CAMBSettings(LTSettings):
def __init__(self,**kwargs):
self.get_scalar_cls = True
self.get_vector_cls = False
self.get_tensor_cls = False
self.get_transfer = True
self.do_lensing = False
self.do_nonlinear = 0
self.l_max_scalar = 8000
self.k_eta_max_scalar = 16000
self.l_max_tensor = 1500
self.k_eta_max_tensor = 3000
self.use_physical = True
#####################################
self.cs2_lam = 1
#####################################
self.helium_fraction = 0.24
self.nu_mass_eigenstates = 0
self.nu_mass_degeneracies = 0
self.share_delta_neff = True
self.scalar_amplitude = 2.41e-9
self.pivot_scalar = 0.002 * u.Mpc**-1
self.pivot_tensor = 0.002 * u.Mpc**-1
#################################################
self.reionization = True
self.re_use_optical_depth = True
self.re_optical_depth = 0.087
self.re_redshift = 11
self.re_delta_redshift = 0.5
self.re_ionization_frac = -1
self.recfast_fudge = 1.14
self.recfast_fudge_he = 0.86
self.recfast_heswitch = 6
self.initial_condition = 1
###########################################################
self.initial_vector = np.array([-1,0,0,0,0])
###########################################################
self.vector_mode = 0
self.cobe_normalize = False
self.cmb_outputscale = 7.4311e12
self.transfer_high_precision = True
self.transfer_kmax = 1000
self.transfer_k_per_logint = 100
################################################################
self.transfer_interp_matterpower = True
#############################################################
self.transfer_power_var = 7
self.scalar_output_file = "scalCls.dat"
self.vector_output_file = "vecCls.dat"
self.tensor_output_file = "tensCls.dat"
self.total_output_file = "totCls.dat"
self.lensed_output_file = "lensedCls.dat"
self.lensed_total_output_file = "lensedtotCls.dat"
self.fits_filename = "scalCls.fits"
###############################################################
self.feedback_level = 1
self.lensing_method = 1
self.accurate_bb = True
self.massive_nu_approx = 3
self.accurate_polarization = True
self.accurate_reionization = True
self.do_tensor_neutrinos = False
self.do_late_rad_truncation = False
self.number_of_threads = 0
self.accuracy_boost = 3
self.l_accuracy_boost = 3
self.l_sample_boost = 3
###############################################################
#Allow for kwargs override
for key in kwargs.keys():
setattr(self,key,kwargs[key])
[docs] def write(self,output_root,cosmology,redshifts):
"""
Writes a CAMB parameter file
:param output_root: output_root for the files that CAMB will produce in output
:type output_root: str.
:param cosmology: cosmological model to generate the parameter file for
:type cosmology: FLRW
:param redshifts: redshifts on which to compute the matter power spectrum and transfer function
:type redshifts: array.
:returns: string object
:rtype: StringIO
"""
#Safety type check
assert isinstance(cosmology,FLRW)
#Sort the redshifts in chronological order
z = -1.0*np.sort(-1.0*redshifts)
#Instantiate StringIO object
s = StringIO()
s.write("output_root = {0}\n".format(output_root))
s.write("\n\n#####################################\n\n")
s.write('get_scalar_cls = {0}\n'.format(self.get_scalar_cls.__str__()[0]))
s.write('get_vector_cls = {0}\n'.format(self.get_vector_cls.__str__()[0]))
s.write('get_tensor_cls = {0}\n'.format(self.get_tensor_cls.__str__()[0]))
s.write('get_transfer = {0}\n'.format(self.get_transfer.__str__()[0]))
s.write('do_lensing = {0}\n'.format(self.do_lensing.__str__()[0]))
s.write('do_nonlinear = {0}\n'.format(self.do_nonlinear))
s.write('l_max_scalar = {0}\n'.format(self.l_max_scalar))
s.write('k_eta_max_scalar = {0}\n'.format(self.k_eta_max_scalar))
s.write('l_max_tensor = {0}\n'.format(self.l_max_tensor))
s.write('k_eta_max_tensor = {0}\n'.format(self.k_eta_max_tensor))
s.write("\n\n#####################################\n\n")
s.write('use_physical = {0}\n'.format(self.use_physical.__str__()[0]))
#############################################
#######Cosmological parameters###############
#############################################
#Baryon and dark matter densities
s.write("ombh2 = {0:.6f}\n".format(cosmology.Ob0*(cosmology.h**2)))
s.write("omch2 = {0:.6f}\n".format((cosmology.Om0 - cosmology.Ob0)*(cosmology.h**2)))
#Neutrino density
if cosmology._nmassivenu==0:
omnuh2 = 0.0
else:
omnuh2 = cosmology.Onu0 * (cosmology.h**2)
s.write("omnuh2 = {0:.6f}\n".format(omnuh2))
#Curvature parameter (enforce Om+Ol=1-Ok)
s.write("omk = {0:.6f}\n".format(1-cosmology.Om0-cosmology.Ode0))
#Hubble constant
s.write("hubble = {0:.6f}\n".format(cosmology.h * 100))
#Dark energy parameters
if hasattr(cosmology,"w0"):
w0 = cosmology.w0
else:
w0 = -1.
if hasattr(cosmology,"wa"):
wa = cosmology.wa
else:
wa = 0.
s.write("w = {0:.6f}\n".format(w0))
s.write("wa = {0:.6f}\n".format(wa))
s.write("\n\n#####################################\n\n")
s.write('cs2_lam = {0}\n'.format(self.cs2_lam))
s.write('temp_cmb = {0:.3f}\n'.format(cosmology.Tcmb0.to(u.K).value))
s.write("\n\n#####################################\n\n")
s.write('helium_fraction = {0}\n'.format(self.helium_fraction))
s.write('massless_neutrinos = {0}\n'.format(cosmology.Neff-cosmology._nmassivenu))
s.write('massive_neutrinos = {0}\n'.format(cosmology._nmassivenu))
s.write('nu_mass_eigenstates = {0}\n'.format(self.nu_mass_eigenstates))
s.write('nu_mass_degeneracies = {0}\n'.format(self.nu_mass_degeneracies))
#Compute the mass fractions of the massive species
if cosmology._nmassivenu:
fractions = (cosmology.m_nu / cosmology.m_nu.sum()).decompose().value
else:
fractions = 1
s.write('nu_mass_fractions = {0}\n'.format(fractions.__str__().strip("[").strip("]")))
s.write('share_delta_neff = {0}\n'.format(self.share_delta_neff.__str__()[0]))
s.write("\n\n#####################################\n\n")
#############################################
#######Spectral index tilt ##################
#############################################
s.write('pivot_scalar = {0:.3f}\n'.format(self.pivot_scalar.to(u.Mpc**-1).value))
s.write('pivot_tensor = {0:.3f}\n'.format(self.pivot_tensor.to(u.Mpc**-1).value))
s.write('initial_power_num = {0}\n'.format(1))
s.write('scalar_amp(1) = {0:.6e}\n'.format(self.scalar_amplitude))
if hasattr(cosmology,"ns"):
ns = cosmology.ns
else:
ns = 1.0
s.write('scalar_spectral_index(1) = {0:.6f}\n'.format(ns))
s.write('scalar_nrun(1) = {0}\n'.format(0))
s.write('tensor_spectral_index(1) = {0}\n'.format(0))
s.write('initial_ratio(1) = {0}\n'.format(0))
s.write("\n\n#####################################\n\n")
s.write('reionization = {0}\n'.format(self.reionization.__str__()[0]))
s.write('re_use_optical_depth = {0}\n'.format(self.re_use_optical_depth.__str__()[0]))
s.write('re_optical_depth = {0:.3f}\n'.format(self.re_optical_depth))
s.write('re_redshift = {0}\n'.format(self.re_redshift))
s.write('re_delta_redshift = {0:.2f}\n'.format(self.re_delta_redshift))
s.write('re_ionization_frac = {0}\n'.format(self.re_ionization_frac))
s.write("\n\n#####################################\n\n")
s.write('RECFAST_fudge = {0}\n'.format(self.recfast_fudge))
s.write('RECFAST_fudge_he = {0}\n'.format(self.recfast_fudge_he))
s.write('RECFAST_heswitch = {0}\n'.format(self.recfast_heswitch))
s.write("\n\n#####################################\n\n")
s.write('initial_condition = {0}\n'.format(self.initial_condition))
s.write('initial_vector = {0}\n'.format(self.initial_vector.__str__().strip("[").strip("]")))
s.write("\n\n#####################################\n\n")
s.write('vector_mode = {0}\n'.format(self.vector_mode))
s.write('cobe_normalize = {0}\n'.format(self.cobe_normalize.__str__()[0]))
s.write("\n\n#####################################\n\n")
s.write('cmb_outputscale = {0:4e}\n'.format(self.cmb_outputscale))
s.write("\n\n#####################################\n\n")
s.write('transfer_high_precision = {0}\n'.format(self.transfer_high_precision.__str__()[0]))
s.write('transfer_kmax = {0}\n'.format(self.transfer_kmax))
s.write('transfer_k_per_logint = {0}\n'.format(self.transfer_k_per_logint))
s.write('transfer_interp_matterpower = {0}\n'.format(self.transfer_interp_matterpower.__str__()[0]))
#############################################
#######Transfer function ####################
#############################################
s.write("transfer_num_redshifts = {0}\n\n".format(len(z)))
for n in range(len(z)):
s.write("transfer_redshift({0}) = {1:.6f}\n".format(n+1,z[n]))
s.write("transfer_filename({0}) = transferfunc_z{1:.6f}.dat\n".format(n+1,z[n]))
s.write("transfer_matterpower({0}) = matterpower_z{1:.6f}.dat\n".format(n+1,z[n]))
s.write("\n\n#####################################\n\n")
s.write('transfer_power_var = {0}\n'.format(self.transfer_power_var))
s.write("\n\n#####################################\n\n")
s.write('scalar_output_file = {0}\n'.format(self.scalar_output_file))
s.write('vector_output_file = {0}\n'.format(self.vector_output_file))
s.write('tensor_output_file = {0}\n'.format(self.tensor_output_file))
s.write('total_output_file = {0}\n'.format(self.total_output_file))
s.write('lensed_output_file = {0}\n'.format(self.lensed_output_file))
s.write('lensed_total_output_file = {0}\n'.format(self.lensed_total_output_file))
s.write('fits_filename = {0}\n'.format(self.fits_filename))
s.write("\n\n#####################################\n\n")
s.write('feedback_level = {0}\n'.format(self.feedback_level))
s.write('lensing_method = {0}\n'.format(self.lensing_method))
s.write('accurate_bb = {0}\n'.format(self.accurate_bb.__str__()[0]))
s.write("\n\n#####################################\n\n")
s.write('massive_nu_approx = {0}\n'.format(self.massive_nu_approx))
s.write("\n\n#####################################\n\n")
s.write('accurate_polarization = {0}\n'.format(self.accurate_polarization.__str__()[0]))
s.write('accurate_reionization = {0}\n'.format(self.accurate_reionization.__str__()[0]))
s.write('do_tensor_neutrinos = {0}\n'.format(self.do_tensor_neutrinos.__str__()[0]))
s.write('do_late_rad_truncation = {0}\n'.format(self.do_late_rad_truncation.__str__()[0]))
s.write("\n\n#####################################\n\n")
s.write('number_of_threads = {0}\n'.format(self.number_of_threads))
s.write('accuracy_boost = {0}\n'.format(self.accuracy_boost))
s.write('l_accuracy_boost = {0}\n'.format(self.l_accuracy_boost))
s.write('l_sample_boost = {0}\n'.format(self.l_sample_boost))
s.seek(0)
return s.read()
#########################################################################################################################################
#CAMB transfer function
class TransferFunction(object):
def __init__(self,k):
"""
:param k: wavenumbers at which the transfer function is computed at
:type k: quantity
"""
assert k.unit.physical_type=="wavenumber"
self._k = k.to((u.Mpc)**-1)
self._transfer = dict()
self._interpolated = dict()
def add(self,z,T):
"""
Add transfer function information at redshift z
:param z: redshift
:type z: float.
:param T: CDM transfer function from CAMB output
:type T: array
"""
if hasattr(self,"_sorted_z"):
del(self._sorted_z)
assert T.shape==self._k.shape,"There should be exactly one transfer function value for each wavenumber! len(T)={0} len(k)={1}".format(len(T),len(self._k))
self._transfer[z] = T
def __getitem__(self,z):
"""
Returns the tabulated transfer function at z. If z is not in the table, returns the tabulated transfer function at the closest z available
:param z: redshift at which to output the tabulated transfer function
:type z: float.
:returns: (tabulated z,k,tabulated transfer function)
:rtype: tuple.
"""
#If the transfer function is not tabulated with z, use the closest z in the table
if not hasattr(self,"_sorted_z"):
self._sorted_z = np.sort(np.array(list(self._transfer.keys())))
if z in self._transfer:
zt = z
else:
zt = self._sorted_z[np.abs(self._sorted_z - z).argmin()]
#Return
return zt,self._k,self._transfer[zt]
def __call__(self,z,k):
"""
Compute the transfer function at redshift z by linear interpolation
:param z: redshift
:type z: float.
:param k: wavenumbers at which to compute the transfer function (linearly interpolated with scipy.interp1d)
:type k: quantity
:returns: transfer function at k
:rtype: array
"""
assert k.unit.physical_type=="wavenumber"
#If the transfer function is not tabulated with z, use the closest z in the table
if not hasattr(self,"_sorted_z"):
self._sorted_z = np.sort(np.array(list(self._transfer.keys())))
if z in self._transfer:
zt = z
else:
zt = self._sorted_z[np.abs(self._sorted_z - z).argmin()]
#If interpolator has not been built yet for the current redshift, build it
if zt not in self._interpolated:
self._interpolated[zt] = interpolate.interp1d(self._k.value,self._transfer[zt],fill_value=1,bounds_error=False)
#Use interpolator to compute the transfer function
return self._interpolated[zt](k.to((u.Mpc)**-1).value)
#I/O
def save(self,filename):
"""
Pickle the TransferFunction instance
:param filename: name of the file to save the instance to
:type filename: str.
"""
with open(filename,"wb") as fp:
pkl.dump(self,fp,protocol=2)
@classmethod
def read(cls,filename):
"""
Load a previously pickled TransferFunction instance
:param filename: name of the file from which the instance should be read
:type filename: str.
:rtype: :py:class:`TransferFunction`
"""
with open(filename,"rb") as fp:
tfr = pkl.load(fp)
if isinstance(tfr,cls):
return tfr
else:
raise TypeError("Pickled instance is not of type {0}".format(cls.__name__))
##############################################################################################################################
[docs]class CAMBTransferFunction(TransferFunction):
pass
[docs]class CAMBTransferFromPower(TransferFunction):
[docs] def add(self,z,T):
"""
Add transfer function information at redshift z
:param z: redshift
:type z: float.
:param T: CDM transfer function from CAMB output
:type T: array
"""
if hasattr(self,"_sorted_z"):
del(self._sorted_z)
assert T.shape==self._k.shape,"There should be exactly one transfer function value for each wavenumber! len(T)={0} len(k)={1}".format(len(T),len(self._k))
self._transfer[z] = np.sqrt(T)
##############################################################################################################################
#k independent transfer function for testing D(z,k) = 1/1+z
class TestTransferFunction(TransferFunction):
def __call__(self,z,k):
return np.ones(k.shape)/(1+z)