diff --git a/py/BaseCosmology.py b/Cosmo/BaseCosmology.py similarity index 80% rename from py/BaseCosmology.py rename to Cosmo/BaseCosmology.py index bceb2a0..c8ae41b 100644 --- a/py/BaseCosmology.py +++ b/Cosmo/BaseCosmology.py @@ -2,31 +2,29 @@ # parameterization of the equation of state or densities or anything. # However, it does know about Hubble's constant at z=0 OR the prefactor # c/(H0*rd) which should be fit for in the case of "rd agnostic" fits. -# That is why you should let it declare those parameterd based on its settings +# That is why you should let it declare those parameters based on its settings ## # However, to get the angular diameter distance you need to pass it -# its Curvature parameter (Omega_k basically), so you need to update -# it +# its Curvature parameter (Omega_k basically), so you need to update it +# -from scipy import * -from RadiationAndNeutrinos import * -from Parameter import * -from ParamDefs import * + +import scipy as sp +from scipy import constants import scipy.integrate as intg -from scipy.interpolate import griddata -from scipy.integrate import odeint -import CosmoApprox as CA + +from ParamDefs import h_par, Pr_par class BaseCosmology: - # speed of light - c_ = 299792.45 + # speed of light in km s^-1 + c_ = constants.c/1000. def __init__(self, h=h_par.value): - self.Curv = 0 - self.rd = 149.50 - self.h = h + self.Curv = 0 + self.rd = 149.50 + self.h = h self.prefact = Pr_par.value self.varyPrefactor = False BaseCosmology.updateParams(self, []) @@ -34,21 +32,26 @@ def __init__(self, h=h_par.value): def setCurvature(self, R): self.Curv = R + def setrd(self, rd): self.rd = rd + def setVaryPrefactor(self, T=True): self.varyPrefactor = T + def setPrefactor(self, p): self.prefact = p + def prefactor(self): if self.varyPrefactor: return self.prefact else: return self.c_/(self.rd*self.h*100) + def freeParameters(self): if (self.varyPrefactor): Pr_par.setValue(self.prefact) @@ -58,14 +61,17 @@ def freeParameters(self): l = [h_par] return l + def printFreeParameters(self): print("Free parameters:") self.printParameters(self.freeParameters()) + def printParameters(self, params): for p in params: print(p.name, '=', p.value, '+/-', p.error) + def updateParams(self, pars): for p in pars: if p.name == "h": @@ -74,26 +80,32 @@ def updateParams(self, pars): self.setPrefactor(p.value) # h shouldn't matter here # we do not want it to enter secondarily through - # say neutrinos, so let's keep it - # sane + # say neutrinos, so let's keep it sane + # # self.h=p.value*self.rd*100/self.c_ return True + def prior_loglike(self): return 0 + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): print("You should not instatiate BaseCosmology") - error("BAD") + print("BAD") + return 0 + def Hinv_z(self, z): - return 1/sqrt(self.RHSquared_a(1.0/(1+z))) + return 1./sp.sqrt(self.RHSquared_a(1.0/(1+z))) + # @autojit def DistIntegrand_a(self, a): - return 1/sqrt(self.RHSquared_a(a))/a**2 + return 1./sp.sqrt(self.RHSquared_a(a))/a**2 + # @autojit def Da_z(self, z): @@ -105,47 +117,51 @@ def Da_z(self, z): if self.Curv == 0: return r elif (self.Curv > 0): - q = sqrt(self.Curv) + q = sp.sqrt(self.Curv) # someone check this eq # Pure ADD has a 1+z fact, but have # comoving one - return sinh(r*q)/(q) + return sp.sinh(r*q)/(q) else: - q = sqrt(-self.Curv) - return sin(r*q)/(q) + q = sp.sqrt(-self.Curv) + return sp.sin(r*q)/(q) + # D_a / rd def DaOverrd(self, z): return self.prefactor()*self.Da_z(z) - # H^{-1} / rd + + # H^{-1} / rd def HIOverrd(self, z): return self.prefactor()*self.Hinv_z(z) - # Dv / rd + + # Dv / rd def DVOverrd(self, z): return self.prefactor()*(self.Da_z(z)**(2./3.)*(z*self.Hinv_z(z))**(1./3.)) + # distance modulus def distance_modulus(self, z): - # I think this should also work with varyPrefactor as long as BAO is there too # assert(not self.varyPrefactor) # note that our Da_z is comoving, so we're only # multilpyting with a single (1+z) factor - return 5*log10(self.Da_z(z)*(1+z)) + return 5*sp.log10(self.Da_z(z)*(1+z)) - # returns the growth factor as a function of redshift + # returns the growth factor as a function of redshift def GrowthIntegrand_a(self, a): - return 1/(self.RHSquared_a(a)*a*a)**(1.5) + return 1./(self.RHSquared_a(a)*a*a)**(1.5) + def growth(self, z): # Equation 7.80 from Doddie af = 1/(1.+z) r = intg.quad(self.GrowthIntegrand_a, 1e-7, af) - gr = sqrt(self.RHSquared_a(af))*r[0] # assume precision is ok + gr = sp.sqrt(self.RHSquared_a(af))*r[0] # assume precision is ok # If we have Omega_m, let's normalize that way if hasattr(self, "Om"): gr *= 5/2.*self.Om diff --git a/py/CosmoApprox.py b/Cosmo/CosmoApprox.py similarity index 97% rename from py/CosmoApprox.py rename to Cosmo/CosmoApprox.py index 8eea9f8..b23309b 100644 --- a/py/CosmoApprox.py +++ b/Cosmo/CosmoApprox.py @@ -2,6 +2,7 @@ # and rd calculators in the second section ## +import sys import math as N Tcmb = 2.7255 @@ -72,6 +73,7 @@ def rd_anderson_approx(obh2, ocbh2, onuh2, Nnu): if (abs(Nnu-3) > 0.1): print("ERROR, cannot use anderson approx with Nnu") print("Nnu=", Nnu) + sys.exit(1) return 55.234 / (ocbh2**0.2538 * obh2**0.1278 * (1+onuh2)**0.3794) @@ -79,7 +81,7 @@ def rd_cuesta_approx(obh2, ocbh2, onuh2, Nnu): if (abs(Nnu-3) > 0.1): print("ERROR, Tony Cuesta says: 'not in this ceral box.'") print("Nnu=", Nnu) - stop() + sys.exit(1) return 55.154 / (ocbh2**0.25351 * (obh2)**0.12807 * N.exp( (onuh2+0.0006)**2.0/0.1176**2)) @@ -92,5 +94,5 @@ def rd_EH_approx(obh2, ocbh2, onuh2, Nnu): if (abs(Nnu-3) > 0.1): print("ERROR, cannot use EH approx with Nnu.") print("Nnu=", Nnu) - stop() + sys.exit(1) return self.CA.soundhorizon_eh(ocbh2, obh2, Nnu) diff --git a/py/NuDensity.py b/Cosmo/NuDensity.py similarity index 71% rename from py/NuDensity.py rename to Cosmo/NuDensity.py index 959fb02..c158844 100644 --- a/py/NuDensity.py +++ b/Cosmo/NuDensity.py @@ -3,10 +3,13 @@ # of neutrino energy densities. # -from scipy import * + +import scipy as sp +from scipy import constants as ct +from scipy.special import zeta from scipy.integrate import quad from scipy.interpolate import interp1d -from CosmoApprox import * +import CosmoApprox as CA #from numba import autojit @@ -14,29 +17,30 @@ class NuIntegral: def __init__(self): "this integrates the stupid integral" print("Initalizing nu density look up table...", end=' ') - rat = 10**(arange(-4, 5, 0.1)) + rat = 10**(sp.arange(-4, 5, 0.1)) intg = [] for r in rat: # min below is to supress the stupid overlow warning - res = quad(lambda x: sqrt(x*x+r**2) / - (exp(min(x, 400))+1.0)*x**2, 0, 1000) + res = quad(lambda x: sp.sqrt(x*x+r**2) / + (sp.exp(min(x, 400))+1.0)*x**2, 0, 1000) intg.append(res[0]/(1+r)) - intg = array(intg) + intg = sp.array(intg) intg *= 7/8./intg[0] # the right normalization - self.interpolator = interp1d(log(rat), intg) + self.interpolator = interp1d(sp.log(rat), intg) # type this into maple - # evalf(45*Zeta(3)/(2*Pi^4)); - self.int_infty = 0.2776566337 + # evalf(45*Zeta(3)/(2*Pi^4)); 0.2776566337 + self.int_infty = 45*zeta(3)/(2*ct.pi**4) print("Done") # self.int_infty,intg[-1]*(1+r)/r + def SevenEights(self, mnuOT): if (mnuOT < 1e-4): return 7/8. elif (mnuOT > 1e4): - return self.int_infty*mnuOT # I don't think this every matters + return self.int_infty*mnuOT # I don't think this ever matters else: - return self.interpolator(log(mnuOT))*(1+mnuOT) + return self.interpolator(sp.log(mnuOT))*(1+mnuOT) class ZeroNuDensity: @@ -63,45 +67,49 @@ def __init__(self, TCMB, Nnu=3.046, mnu=0.06, degenerate=False): # one neutrino worth of radiation at the nominal temperature, or heated on? # See CAMB notes Eq. 4-7. - self.gfact = (3.046/3.0) + self.gfact = (3.046/3.0) self.gfact_o4 = self.gfact**(0.25) # ideal neutrino temp - self.Tnu0 = (4./11.)**(1./3.)*TCMB + self.Tnu0 = (4./11.)**(1./3.)*TCMB # actual neutrino temp - self.Tnu = self.Tnu0*self.gfact_o4 + self.Tnu = self.Tnu0*self.gfact_o4 # same for prefactors - self.prefix0 = 4.48130979e-7 * TCMB**4 * ((4./11.)**(4./3.)) - self.prefix = self.prefix0*self.gfact + self.prefix0 = 4.48130979e-7 * TCMB**4 * ((4./11.)**(4./3.)) + self.prefix = self.prefix0*self.gfact self.degenerate = degenerate self.set_mnuone_() + def set_mnuone_(self): if self.degenerate: - self.mnuone = self.mnu_/self.Nnu_ + self.mnuone = self.mnu_/self.Nnu_ else: - self.mnuone = self.mnu_*1.0 + self.mnuone = self.mnu_*1.0 self.omnuh2today = self.rho(1) + def setMnu(self, mnu): self.mnu_ = mnu self.set_mnuone_() + def setNnu(self, Nnu): self.Nnu_ = Nnu self.set_mnuone_() - # @autojit + # @autojit def rho(self, a): # This returns the density at a normalized so that # we get nuh2 at a=0 - # (1 eV) / (Boltzmann constant * 1 kelvin) = - # 11 604.5193 + # (1 eV) / (Boltzmann constant * 1 kelvin) = 11 604.5193 + if (self.mnuone == 0): return self.Nnu_*7/8.*self.prefix0/a**4 - mnuOT = self.mnuone/(self.Tnu/a)*11604.5193 + mnuOT = self.mnuone/(self.Tnu/a)*(1./ct.value(u'Boltzmann constant in eV/K')) #11604.5193 + # Here for massive we use 1*prefix (accounting for 1.015 in Tnu) # For massles we use Neff*prefix0 (so we account if self.degenerate: @@ -110,16 +118,17 @@ def rho(self, a): return ((self.I.SevenEights(mnuOT)*self.prefix+(self.Nnu_-1.015)*7/8.*self.prefix0))/a**4 + if __name__ == "__main__": - A = NuDensity(Tcmb, 3.046, 0.06) + A = NuDensity(CA.Tcmb, 3.046, 0.06) print(A.omnuh2today, '=including massless neutrinos') - B = NuDensity(Tcmb, 2.030, 0.00) + B = NuDensity(CA.Tcmb, 2.030, 0.00) print(B.omnuh2today, end=' ') print(A.omnuh2today-B.omnuh2today, '=excluding massless neutrinos') - A = NuDensity(Tcmb, 1.015, 60.) + A = NuDensity(CA.Tcmb, 1.015, 60.) print(A.omnuh2today/1000., '=assuming very cold') - A = NuDensity(Tcmb, 1.015, 0.06) + A = NuDensity(CA.Tcmb, 1.015, 0.06) print(A.omnuh2today, '=assuming real temperature') diff --git a/py/ParamDefs.py b/Cosmo/ParamDefs.py similarity index 61% rename from py/ParamDefs.py rename to Cosmo/ParamDefs.py index d1b0404..038558a 100644 --- a/py/ParamDefs.py +++ b/Cosmo/ParamDefs.py @@ -1,23 +1,23 @@ ## -# This class has parameter defintions for all -# parameter used in this code. +# This class has parameter definitions for all +# parameters used in this code. ## # Change here for bounds, or import and rewrite. ## ## -from Parameter import * +from Parameter import Parameter # Parameters are value, variation, bounds -Om_par = Parameter("Om", 0.3038, 0.1, (0.05, 1.5), "\Omega_m") +Om_par = Parameter("Om", 0.3038, 0.1, (0.05, 0.5), "\Omega_m") Obh2_par = Parameter("Obh2", 0.02234, 0.0002, (0.02, 0.025), "\Omega_{b}h^2") -h_par = Parameter("h", 0.6821, 0.05, (0.4, 1.0), "h") -mnu_par = Parameter("mnu", 0.06, 0.1, (0, 1.0), "\Sigma m_{\\nu}") -Nnu_par = Parameter("Nnu", 3.046, 0.5, (3.046, 5.046), "N_{\\rm eff}") +h_par = Parameter("h", 0.6821, 0.05, (0.4, 1.0), "h") +mnu_par = Parameter("mnu", 0.06, 0.1, (0, 1.0), "\Sigma m_{\\nu}") +Nnu_par = Parameter("Nnu", 3.046, 0.5, (3.046, 5.046),"N_{\\rm eff}") -Ok_par = Parameter("Ok", 0.0, 0.1, (-1.5, 1.5), "\Omega_k") -w_par = Parameter("w", -1.0, 0.1, (-2.0, 0.0), "w_0") +Ok_par = Parameter("Ok", 0.0, 0.1, (-0.5, 0.5), "\Omega_k") +w_par = Parameter("w", -1.0, 0.1, (-2.0, 0.0), "w_0") wa_par = Parameter("wa", 0.0, 0.1, (-2.0, 2.0), "w_a") @@ -29,16 +29,16 @@ Om2_par = Parameter("Om2", 0.0, 0.7, (-3, 3), "\Omega_2") # JordiCDM Cosmology Parameters -q_par = Parameter("q", 0.0, 0.2, (0, 1), "q") -za_par = Parameter("za", 3, 1.0, (2, 10), "z_a") -zb_par = Parameter("zb", 1, 0.5, (0, 2), "z_b") -wd_par = Parameter("wd", -1, 0.5, (-2.0, 0.0), "w_d") -Od_par = Parameter("Od", 0.0, 0.2, (0.0, 1.0), "O_d") +q_par = Parameter("q", 0.0, 0.2, (0, 1), "q") +za_par = Parameter("za", 3, 1.0, (2, 10), "z_a") +zb_par = Parameter("zb", 1, 0.5, (0, 2), "z_b") +wd_par = Parameter("wd", -1, 0.5, (-2.0, 0.0), "w_d") +Od_par = Parameter("Od", 0.0, 0.2, (0.0, 1.0), "O_d") # Weird model -mu_par = Parameter("mu", 2.13, 0.1, (0.01, 5), "\mu") -Amp_par = Parameter("Amp", -0.1839, 0.01, (-2, 2), "Amp") -sig_par = Parameter("sig", 0.1, 0.2, (0, 3.0), "\sig") +mu_par = Parameter("mu", 2.13, 0.1, (0.01, 5), "\mu") +Amp_par = Parameter("Amp", -0.1839, 0.01, (-2, 2), "Amp") +sig_par = Parameter("sig", 0.1, 0.2, (0, 3.0), "\sig") # Tired Light beta_par = Parameter("beta", 1, 0.1, (-1, 1), "\\beta") @@ -70,7 +70,7 @@ # Decaying Dark Matter lambda_par = Parameter("lambda", 1.0, 1.0, (0., 20.0), "\lambda") -xfrac_par = Parameter("xfrac", 0.1, 0.1, (0.0, 1.0), "f_x") +xfrac_par = Parameter("xfrac", 0.1, 0.1, (0.0, 1.0), "f_x") # Early Dark Energy Ode_par = Parameter("Ode", 0.05, 0.01, (0, 0.9), "\Omega^e_{\\rm de}") @@ -79,10 +79,10 @@ dw_par = Parameter("dw", 0, 0.1, (-1.0, 1.0), "\delta w_0") # QuinDE -lam_par = Parameter("lam", 1, 0.5, (0, 5), "\lambda") -V0_par = Parameter("V0", 1, 1, (0, 30), "V_0") -A_par = Parameter("A", 10, 1, (-20, 20), "A") -B_par = Parameter("B", -20, 2, (-100, 100), "B") +lam_par = Parameter("lam", 1, 0.5, (0, 5), "\lambda") +V0_par = Parameter("V0", 1, 1, (0, 30), "V_0") +A_par = Parameter("A", 10, 1, (-20, 20), "A") +B_par = Parameter("B", -20, 2, (-100, 100), "B") # wDark Matter wDM_par = Parameter("wDM", 0.0, 0.1, (-1.0, 1.0), "w_{DM}") diff --git a/py/Parameter.py b/Cosmo/Parameter.py similarity index 91% rename from py/Parameter.py rename to Cosmo/Parameter.py index ab5a0b3..d2e2879 100644 --- a/py/Parameter.py +++ b/Cosmo/Parameter.py @@ -5,7 +5,7 @@ class Parameter: - def __init__(self, name, value, err=0, bounds=None, Ltxname=None): + def __init__(self, name, value, err=0.0, bounds=None, Ltxname=None): self.name = name if Ltxname: self.Ltxname = Ltxname @@ -19,17 +19,22 @@ def __init__(self, name, value, err=0, bounds=None, Ltxname=None): else: self.bounds = bounds + def sameParam(self, param2): return self.name == param2.name + def setLatexName(self, Ltx): self.Ltxname = Ltx + def setValue(self, val): self.value = val + def setError(self, err): self.error = err + def setBounds(self, low, high): self.bounds = [low, high] diff --git a/py/RadiationAndNeutrinos.py b/Cosmo/RadiationAndNeutrinos.py similarity index 90% rename from py/RadiationAndNeutrinos.py rename to Cosmo/RadiationAndNeutrinos.py index efd03d6..3b8696b 100644 --- a/py/RadiationAndNeutrinos.py +++ b/Cosmo/RadiationAndNeutrinos.py @@ -5,8 +5,9 @@ # but it became clutterish there. ## -from NuDensity import * -from ParamDefs import * +import sys +from ParamDefs import mnu_par, Nnu_par +from NuDensity import ZeroNuDensity, NuDensity import CosmoApprox as CA @@ -22,9 +23,9 @@ def __init__(self, mnu=mnu_par.value, Nnu=Nnu_par.value, varyMnu=False, varyNnu=False, degenerate=False, disable=False): self.disabled = disable if (self.disabled): - self.Omrad = 0 - self.Omnuh2 = 0 - self.Omnu = 0 + self.Omrad = 0 + self.Omnuh2 = 0 + self.Omnu = 0 self.NuDensity = ZeroNuDensity() return self.varyMnu = False @@ -32,40 +33,47 @@ def __init__(self, mnu=mnu_par.value, Nnu=Nnu_par.value, self.NuDensity = NuDensity(CA.Tcmb, Nnu, mnu, degenerate) # print "Relic neutrino density:",self.NuDensity.omnuh2today + def setVaryMnu(self, T=True): if self.disabled: print("Cannot vary radiation parameter if disabled") - stop() + sys.exit(1) self.varyMnu = T + def setVaryNnu(self, T=True): if self.disabled: print("Cannot vary radiation parameter if disabled") - stop() + sys.exit(1) self.varyNnu = T + def setMnu(self, mnu): if self.disabled: print("Cannot vary radiation parameter if disabled") - stop() + sys.exit(1) self.NuDensity.setMnu(mnu) + def setNnu(self, mnu): if self.disabled: print("Cannot vary radiation parameter if disabled") - stop() + sys.exit(1) self.NuDensity.setNnu(mnu) + def mnu(self): if self.disabled: return 0 return self.NuDensity.mnu_ + def Nnu(self): if self.disabled: return 0 return self.NuDensity.Nnu_ + def freeParameters(self): if self.disabled: return [] @@ -78,6 +86,7 @@ def freeParameters(self): l.append(Nnu_par) return l + def updateParams(self, pars): if self.disabled: return True @@ -93,5 +102,4 @@ def updateParams(self, pars): self.Omrad = self.omrad_pref_/(self.h**2) self.Omnuh2 = self.NuDensity.omnuh2today self.Omnu = self.Omnuh2/self.h**2 - return True diff --git a/py/BAOLikelihoods.py b/Likelihoods/BAOLikelihoods.py similarity index 80% rename from py/BAOLikelihoods.py rename to Likelihoods/BAOLikelihoods.py index 113109a..732ee10 100644 --- a/py/BAOLikelihoods.py +++ b/Likelihoods/BAOLikelihoods.py @@ -2,19 +2,19 @@ # The BAO likelihoods. # -from TabulatedBAOLikelihood import * -from TabulatedBAODVLikelihood import * -from GaussBAODVLikelihood import * -from LCDMCosmology import * -from scipy import * +from TabulatedBAOLikelihood import TabulatedBAOLikelihood +from TabulatedBAODVLikelihood import TabulatedBAODVLikelihood +from GaussBAODVLikelihood import GaussBAODVLikelihood +from LCDMCosmology import LCDMCosmology + class DR11LOWZ(GaussBAODVLikelihood): def __init__(self): obh2 = 0.0224 - Om = 0.274 - h = 0.7 - mnu = 0 + Om = 0.274 + h = 0.7 + mnu = 0 fidTheory = LCDMCosmology(obh2, Om, h, mnu) GaussBAODVLikelihood.__init__( self, "DR11LOWZ", 0.32, 1264.0, 25.0, fidTheory) @@ -25,9 +25,9 @@ def __init__(self): # fiducial cosmology for LOWZ/CMASS data. # see Anderson et al, page 28 obh2 = 0.0224 - Om = 0.274 - h = 0.7 - mnu = 0 # rd=149.28 + Om = 0.274 + h = 0.7 + mnu = 0 # rd=149.28 fidTheory = LCDMCosmology(obh2, Om, h, mnu) # negative col means the cols is probability and not chi2 TabulatedBAOLikelihood.__init__(self, "DR11CMASS", 'data/sdss_DR11CMASS_consensus.dat', @@ -40,9 +40,9 @@ def __init__(self): # see e.g. Busca's email on 12/3/13 obh2 = 0.0227 - Om = 0.27 - h = 0.7 - mnu = 0.06 # ;# rd=149.77 + Om = 0.27 + h = 0.7 + mnu = 0.06 # ;# rd=149.77 fidTheory = LCDMCosmology(obh2, Om, h, mnu) # File from 5/16 from Nicolas. TabulatedBAOLikelihood.__init__(self, "DR11LyaAuto", 'data/chi2_surface_dr11_baseline_fit.txt', @@ -52,36 +52,38 @@ def __init__(self): class DR11LyaCross(TabulatedBAOLikelihood): def __init__(self): obh2 = 0.0227 - Om = 0.27 - h = 0.7 - mnu = 0 # ;# rd=149.77 + Om = 0.27 + h = 0.7 + mnu = 0 # ;# rd=149.77 fidTheory = LCDMCosmology(obh2, Om, h, mnu) TabulatedBAOLikelihood.__init__(self, "DR11LyaCross", 'data/lyabaocross.scan', 2, fidTheory, 2.36) + class DR14LyaAuto(TabulatedBAOLikelihood): def __init__(self): # fiducial cosmology for Lya data. # Taken from https://github.com/igmhub/picca/blob/master/data/deSainteAgatheetal2019/auto_alone_stdFit/auto_alone_stdFit..ap.at.scan.dat # fiducial model -- see Table 2 of Victoria's paper obh2 = 0.02222 - h = 0.6731 - Om = 0.1426/h**2 - mnu = 0.06 # rd=147.33 + h = 0.6731 + Om = 0.1426/h**2 + mnu = 0.06 # rd=147.33 fidTheory = LCDMCosmology(obh2, Om, h, mnu) TabulatedBAOLikelihood.__init__(self, "DR14LyaAuto", 'data/deSainteAgatheetal2019_ap_at_scan.dat', 2, fidTheory, 2.34, aperp_col=1, apar_col=0, skiprows=1) + class DR14LyaCross(TabulatedBAOLikelihood): def __init__(self): # fiducial cosmology for Lya data. # Taken from https://github.com/igmhub/picca/tree/master/data/Blomqvistetal2019/cross_alone_stdFit # fiducial model -- double check obh2 = 0.02222 - h = 0.6731 - Om = 0.1426/h**2 - mnu = 0.06 # rd=147.33 + h = 0.6731 + Om = 0.1426/h**2 + mnu = 0.06 # rd=147.33 fidTheory = LCDMCosmology(obh2, Om, h, mnu) # File from 5/16 from Nicolas. TabulatedBAOLikelihood.__init__(self, "DR14LyaCross", 'data/Blomqvistetal2019_ap_at_scan.dat', @@ -96,22 +98,22 @@ def __init__(self): class SixdFGS(GaussBAODVLikelihood): def __init__(self): obh2 = 0.02227 - Om = 0.27 - h = 0.7 - mnu = 0 + Om = 0.27 + h = 0.7 + mnu = 0 fidTheory = LCDMCosmology(obh2, Om, h, mnu) GaussBAODVLikelihood.__init__( self, "SixdFGS", 0.106, 456.0, 27.0, fidTheory, maxchi2=4) -# SDSS Main Galaxy Sample BAO +# SDSS Main Galaxy Sample BAO class SDSSMGS(TabulatedBAODVLikelihood): def __init__(self): obh2 = 0.021547 - Om = 0.31 - h = 0.67 - mnu = 0 + Om = 0.31 + h = 0.67 + mnu = 0 fidTheory = LCDMCosmology(obh2, Om, h, mnu) TabulatedBAODVLikelihood.__init__( self, "MGS", "data/chidavexi8stavePk5staverec.dat", fidTheory, 0.15) diff --git a/py/BaseLikelihood.py b/Likelihoods/BaseLikelihood.py similarity index 94% rename from py/BaseLikelihood.py rename to Likelihoods/BaseLikelihood.py index 8938654..ddc0f46 100644 --- a/py/BaseLikelihood.py +++ b/Likelihoods/BaseLikelihood.py @@ -26,7 +26,7 @@ def updateParams(self, params): def loglike(self): return 0.0 - def theory_loglike_prior(): + def theory_loglike_prior(self): return self.theory_.prior_loglike() def loglike_wprior(self): diff --git a/py/CompositeLikelihood.py b/Likelihoods/CompositeLikelihood.py similarity index 78% rename from py/CompositeLikelihood.py rename to Likelihoods/CompositeLikelihood.py index cb2ef8a..3f6abb8 100644 --- a/py/CompositeLikelihood.py +++ b/Likelihoods/CompositeLikelihood.py @@ -2,43 +2,51 @@ # This thing takes likelihoods and cooks one # which is a composite of many. # -from BaseLikelihood import * -from scipy import * +from BaseLikelihood import BaseLikelihood +import scipy as sp -class CompositeLikelihood (BaseLikelihood): + +class CompositeLikelihood(BaseLikelihood): def __init__(self, llist=[]): BaseLikelihood.__init__(self, "Composite") self.llist_ = llist + def setTheory(self, theory): BaseLikelihood.setTheory(self, theory) for l in self.llist_: l.setTheory(theory) assert(self.theory_ is l.theory_) + def addLikelihood(self, like): self.llist_.append(like) if hasattr(self, 'theory_'): like.setTheory(self.theory_) assert(self.theory_ is like.theory_) + def addLikelihoods(self, likes): for like in likes: self.addLikelihood(like) + def compositeNames(self): return [like.name() for like in self.llist_] + def compositeLogLikes(self): - return array([like.loglike() for like in self.llist_]) + return sp.array([like.loglike() for like in self.llist_]) + def compositeLogLikes_wprior(self): - return array([like.loglike() for like in self.llist_]+[self.theory_.prior_loglike()]) + return sp.array([like.loglike() for like in self.llist_] + [self.theory_.prior_loglike()]) + def loglike(self): likes = [like.loglike() for like in self.llist_] - return array(likes).sum() + return sp.array(likes).sum() # The base si good enough as # Everybody is holding the same reference diff --git a/py/CompressedSNLikelihood.py b/Likelihoods/CompressedSNLikelihood.py similarity index 100% rename from py/CompressedSNLikelihood.py rename to Likelihoods/CompressedSNLikelihood.py diff --git a/py/GaussBAODVLikelihood.py b/Likelihoods/GaussBAODVLikelihood.py similarity index 70% rename from py/GaussBAODVLikelihood.py rename to Likelihoods/GaussBAODVLikelihood.py index 442e5ac..1593f4d 100644 --- a/py/GaussBAODVLikelihood.py +++ b/Likelihoods/GaussBAODVLikelihood.py @@ -3,27 +3,28 @@ # at which chi2 is cut. # -from BaseLikelihood import * -from scipy import * -from scipy.interpolate import RectBivariateSpline +from BaseLikelihood import BaseLikelihood -class GaussBAODVLikelihood (BaseLikelihood): + +class GaussBAODVLikelihood(BaseLikelihood): def __init__(self, name, z, DV, DVErr, fidtheory, maxchi2=1e30): BaseLikelihood.__init__(self, name) self.z = z - rd = fidtheory.rd - DV /= rd + rd = fidtheory.rd + DV /= rd DVErr /= rd self.maxchi2 = maxchi2 print(name, "measurement in ", fidtheory.rd_approx, ":", DV, "+-", DVErr) self.setData(DV, DVErr) + def setData(self, DV, DVErr): - self.DV = DV + self.DV = DV self.DVErr2 = DVErr**2 + def loglike(self): - DVT = self.theory_.DVOverrd(self.z) + DVT = self.theory_.DVOverrd(self.z) chi2 = min(self.maxchi2, (DVT-self.DV)**2/(self.DVErr2)) return -chi2/2.0 diff --git a/py/HubbleParameterLikelihood.py b/Likelihoods/HubbleParameterLikelihood.py similarity index 100% rename from py/HubbleParameterLikelihood.py rename to Likelihoods/HubbleParameterLikelihood.py diff --git a/py/LikelihoodMultiplier.py b/Likelihoods/LikelihoodMultiplier.py similarity index 92% rename from py/LikelihoodMultiplier.py rename to Likelihoods/LikelihoodMultiplier.py index 68c1bfb..96b4fa0 100644 --- a/py/LikelihoodMultiplier.py +++ b/Likelihoods/LikelihoodMultiplier.py @@ -1,5 +1,5 @@ # This is a helper class that takes a likelihood -# and just multiplies loglike with some value to simulater +# and just multiplies loglike with some value to simulate # better or shittier data. ## # You have to give python a credit for how effortlessly @@ -16,6 +16,7 @@ def __init__(self, LikeInstance, mult): self.LikeType = LikeInstance.__class__ self.mult = mult + def loglike(self): return self.LikeType.loglike(self)*self.mult diff --git a/py/SimpleCMB.py b/Likelihoods/SimpleCMB.py similarity index 59% rename from py/SimpleCMB.py rename to Likelihoods/SimpleCMB.py index 0b48f20..5292fbc 100644 --- a/py/SimpleCMB.py +++ b/Likelihoods/SimpleCMB.py @@ -3,15 +3,14 @@ # where we compress to just obh2, obh2+och2 and thetas ## -from BaseLikelihood import * -from scipy import * +from BaseLikelihood import BaseLikelihood import scipy.linalg as la - +import scipy as sp class SimpleCMB (BaseLikelihood): def __init__(self, name, mean, cov, kill_Da=False, kill_rd=False): print("Initializing CMB likelihood:", name) - cov = array(cov) + cov = sp.array(cov) if kill_Da: cov[2, :] = 0.0 cov[:, 2] = 0.0 @@ -23,39 +22,44 @@ def __init__(self, name, mean, cov, kill_Da=False, kill_rd=False): cov[1, 1] = 1e10 name += "_noRd" BaseLikelihood.__init__(self, name) - self.mean = array(mean) - self.cov = cov + self.mean = sp.array(mean) + self.cov = cov self.icov = la.inv(cov) + def setTheory(self, theory): self.theory_ = theory self.theory_.setNoObh2prior() + def loglike(self): - delt = self.theory_.CMBSimpleVec()-self.mean - return -dot(delt, dot(self.icov, delt))/2.0 + delt = self.theory_.CMBSimpleVec() - self.mean + return -sp.dot(delt, sp.dot(self.icov, delt))/2.0 + class PlanckLikelihood(SimpleCMB): def __init__(self, kill_Da=False, kill_rd=False): - mean = array([2.24519776e-02, 1.38572404e-01, 9.43303000e+01]) - cov = array([[1.28572727e-07, -6.03323687e-07, 1.44305285e-05], - [-6.03323687e-07, 7.54205794e-06, -3.60547663e-05], - [1.44305285e-05, -3.60547663e-05, 4.26414740e-03]]) + mean = sp.array([2.24519776e-02, 1.38572404e-01, 9.43303000e+01]) + cov = sp.array([[1.28572727e-07, -6.03323687e-07, 1.44305285e-05], + [-6.03323687e-07, 7.54205794e-06, -3.60547663e-05], + [1.44305285e-05, -3.60547663e-05, 4.26414740e-03]]) name = "SPlanck" SimpleCMB.__init__(self, name, mean, cov, kill_Da, kill_rd) + class WMAP9Likelihood(SimpleCMB): def __init__(self, kill_Da=False, kill_rd=False): - mean = array([2.25946978e-02, 1.35359318e-01, 9.45118918e+01]) - cov = array([[2.86459327e-07, -4.80929954e-07, -1.11081266e-05], - [-4.80929954e-07, 1.90757225e-05, 7.49542945e-06], - [-1.11081266e-05, 7.49542945e-06, 2.54207102e-02]]) + mean = sp.array([2.25946978e-02, 1.35359318e-01, 9.45118918e+01]) + cov = sp.array([[2.86459327e-07, -4.80929954e-07, -1.11081266e-05], + [-4.80929954e-07, 1.90757225e-05, 7.49542945e-06], + [-1.11081266e-05, 7.49542945e-06, 2.54207102e-02]]) name = "SWMAP" SimpleCMB.__init__(self, name, mean, cov, kill_Da, kill_rd) + if __name__ == "__main__": D = PlanckLikelihood() for v in D.mean: diff --git a/py/TabulatedBAODVLikelihood.py b/Likelihoods/TabulatedBAODVLikelihood.py similarity index 71% rename from py/TabulatedBAODVLikelihood.py rename to Likelihoods/TabulatedBAODVLikelihood.py index 4d0b969..36ed94e 100644 --- a/py/TabulatedBAODVLikelihood.py +++ b/Likelihoods/TabulatedBAODVLikelihood.py @@ -1,32 +1,35 @@ ## # This is a DV likelihood that comes in form of a table. ## -from BaseLikelihood import * -from scipy import * + +import scipy as sp +from BaseLikelihood import BaseLikelihood from scipy.interpolate import interp1d from scipy.optimize import minimize from scipy.integrate import quad -class TabulatedBAODVLikelihood (BaseLikelihood): +class TabulatedBAODVLikelihood(BaseLikelihood): def __init__(self, name, filename, fid_theory, z): BaseLikelihood.__init__(self, name) print() print("Loading ", filename) - data = loadtxt(filename) + data = sp.loadtxt(filename) self.chi2i = interp1d(data[:, 0], data[:, 1]) self.fidDVOverrd = fid_theory.DVOverrd(z) + # find minimum alphamin = minimize(self.chi2i, 1.0, bounds=[ [0.9, 1.1]], method='SLSQP').x[0] - rms = quad(lambda x: exp(-self.chi2i(alphamin+x)/2)*x*x, -0.1, 0.1)[0] - rms /= quad(lambda x: exp(-self.chi2i(alphamin+x)/2), -0.1, 0.1)[0] - DVf = self.fidDVOverrd*fid_theory.rd + rms = quad(lambda x: sp.exp(-self.chi2i(alphamin+x)/2)*x*x, -0.1, 0.1)[0] + rms /= quad(lambda x: sp.exp(-self.chi2i(alphamin+x)/2), -0.1, 0.1)[0] + DVf = self.fidDVOverrd*fid_theory.rd print(name, "measurement in", fid_theory.rd_approx, - " : DV=", DVf*alphamin, "+-", sqrt(rms)*DVf) + " : DV=", DVf*alphamin, "+-", sp.sqrt(rms)*DVf) print("with rd=", fid_theory.rd, "DV_fid=", DVf, "alphamin=", alphamin) self.z = z + def loglike(self): alpha = self.theory_.DVOverrd(self.z)/self.fidDVOverrd try: diff --git a/py/TabulatedBAOLikelihood.py b/Likelihoods/TabulatedBAOLikelihood.py similarity index 75% rename from py/TabulatedBAOLikelihood.py rename to Likelihoods/TabulatedBAOLikelihood.py index 469565a..f4683ef 100644 --- a/py/TabulatedBAOLikelihood.py +++ b/Likelihoods/TabulatedBAOLikelihood.py @@ -3,33 +3,35 @@ # as chi2 table. See BAOLikelihoods.py # -from BaseLikelihood import * -from scipy import * + +import scipy as sp +from BaseLikelihood import BaseLikelihood from scipy.interpolate import RectBivariateSpline -class TabulatedBAOLikelihood (BaseLikelihood): + +class TabulatedBAOLikelihood(BaseLikelihood): def __init__(self, name, filename, chi2col, fid_theory, z, skiprows=0, aperp_col=0, apar_col=1): BaseLikelihood.__init__(self, name) print() print("Loading ", filename) - data = loadtxt(filename,skiprows=skiprows) + data = sp.loadtxt(filename, skiprows=skiprows) # first let's autolearn the binning aperp = set() aparl = set() for line in data: aperp.add(line[aperp_col]) aparl.add(line[apar_col]) - aperp = sorted(list(aperp)) - aparl = sorted(list(aparl)) - logltab = zeros((len(aperp), len(aparl))) + aperp = sorted(list(aperp)) + aparl = sorted(list(aparl)) + logltab = sp.zeros((len(aperp), len(aparl))) print("Aperp min,max,step,N:", aperp[0], aperp[-1], aperp[1]-aperp[0], len(aperp)) print("Aparl min,max,step,N:", aparl[0], aparl[-1], aparl[1]-aparl[0], len(aparl)) - self.aperp = array(aperp) - self.aparl = array(aparl) + self.aperp = sp.array(aperp) + self.aparl = sp.array(aparl) # now fill in the table for line in data: @@ -40,22 +42,25 @@ def __init__(self, name, filename, chi2col, fid_theory, z, skiprows=0, aperp_col logltab[ii, jj] = -chi2/2.0 else: # col is probability, add 1e-50 to avoid taking log of zero. - logltab[ii, jj] = log(line[chi2col*-1]+1e-50) + logltab[ii, jj] = sp.log(line[chi2col*-1] + 1e-50) logltab = logltab-logltab.max() self.loglint = RectBivariateSpline( self.aperp, self.aparl, logltab, kx=1, ky=1) print("Loading done") - print("rd = ",fid_theory.rd," Mpc") + print("rd = ", fid_theory.rd, " Mpc") self.fidDaOverrd = fid_theory.DaOverrd(z) self.fidHIOverrd = fid_theory.HIOverrd(z) print("Fiducials at z=", z, ":", self.fidDaOverrd, self.fidHIOverrd) self.z = z + + def loglike_aperp_apar(self, aperp, apar): return self.loglint(aperp, apar)[0][0] + def loglike(self): alphaperp = self.theory_.DaOverrd(self.z)/self.fidDaOverrd - alphapar = self.theory_.HIOverrd(self.z)/self.fidHIOverrd + alphapar = self.theory_.HIOverrd(self.z)/self.fidHIOverrd return self.loglike_aperp_apar(alphaperp, alphapar) diff --git a/py/BinnedWCosmology.py b/Models/BinnedWCosmology.py similarity index 74% rename from py/BinnedWCosmology.py rename to Models/BinnedWCosmology.py index de5838f..4169fad 100644 --- a/py/BinnedWCosmology.py +++ b/Models/BinnedWCosmology.py @@ -1,26 +1,29 @@ -# This is a CDM cosmology with w +# This is a CDM cosmology with binned w +# Still testing the file -- don't trust it -from LCDMCosmology import * -import numpy as np +from LCDMCosmology import LCDMCosmology from scipy.interpolate import interp1d from scipy.integrate import quad - +from Parameter import Parameter +import numpy as np class BinnedWCosmology(LCDMCosmology): def __init__(self, dz=0.2, zmax=1.0): - # two parameters: Om and h - self.zbins = np.arange(0, zmax, dz) - self.Nb = len(self.zbins) - self.wvals = np.ones(self.Nb)*-1.0 + # bunch of parameters: + self.zbins = np.arange(0, zmax, dz) + self.Nb = len(self.zbins) + self.wvals = np.ones(self.Nb)*-1.0 self.pnames = ["w%i" % i for i in range(self.Nb)] LCDMCosmology.__init__(self) self.integrateOmega() + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): - wpars = [Parameter(name, self.wvals[i]) + wpars = [Parameter(name, self.wvals[i], err=0.05) for i, name in enumerate(self.pnames)] - return wpars+LCDMCosmology.freeParameters(self) + return LCDMCosmology.freeParameters(self) + wpars + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) @@ -28,7 +31,10 @@ def updateParams(self, pars): return False gotone = False for p in pars: + ## Something's happening here, check it later + print ('**', p.name, self.pnames) i = self.pnames.index(p.name) + if i > 0: self.wvals[i] = p.value gotone = True @@ -36,6 +42,7 @@ def updateParams(self, pars): self.integrateOmega() return True + def integrateOmega(self): abins = np.hstack((1./(1+self.zbins), [1e-4])) w = np.hstack((self.wvals, [self.wvals[-1]])) @@ -46,9 +53,9 @@ def integrateOmega(self): print(np.exp(olnrho)) self.DEomega = interp1d(oabins, np.exp(olnrho)) + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 - def RHSquared_a(self, a): NuContrib = self.NuDensity.rho(a)/self.h**2 return (self.Ocb/a**3+self.Omrad/a**4+NuContrib+(1.0-self.Om)*self.DEomega(a)) diff --git a/py/DecayLCDMCosmology.py b/Models/DecayLCDMCosmology.py similarity index 77% rename from py/DecayLCDMCosmology.py rename to Models/DecayLCDMCosmology.py index 069a0b0..d154662 100644 --- a/py/DecayLCDMCosmology.py +++ b/Models/DecayLCDMCosmology.py @@ -2,12 +2,13 @@ # dark matter component. ## -from pylab import * +import sys import numpy as np from scipy.integrate import odeint from scipy.interpolate import interp1d from scipy.optimize import minimize -from LCDMCosmology import * +from LCDMCosmology import LCDMCosmology +from ParamDefs import xfrac_par, lambda_par class DecayLCDMCosmology(LCDMCosmology): @@ -15,7 +16,7 @@ class DecayLCDMCosmology(LCDMCosmology): # density at early a is zero. def __init__(self, varylam=True, varyxfrac=True, xfrac=xfrac_par.value): - self.varylam = varylam + self.varylam = varylam self.varyxfrac = varyxfrac self.lam = lambda_par.value @@ -23,22 +24,21 @@ def __init__(self, varylam=True, varyxfrac=True, xfrac=xfrac_par.value): LCDMCosmology.__init__(self) - self.logar = linspace(0.0, -7.1, 100) + self.logar = np.linspace(0.0, -7.1, 100) self.ilogar = self.logar[::-1] # force caching self.updateParams([]) - # my free parameters. We add lam, xfrac + # my free parameters. We add lam, xfrac def freeParameters(self): l = LCDMCosmology.freeParameters(self) - if (self.varylam): - l.append(lambda_par) - if (self.varyxfrac): - l.append(xfrac_par) + if (self.varylam): l.append(lambda_par) + if (self.varyxfrac): l.append(xfrac_par) return l + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) if not ok: @@ -55,26 +55,29 @@ def updateParams(self, pars): assert(abs(self.RHSquared_a(1.0)-1) < 1e-4) return True + def H2_rxrr_a(self, a, rx, rr): NuContrib = self.NuDensity.rho(a)/self.h**2 return self.Ocb_std/a**3+self.Omrad/a**4+NuContrib+(1.0-self.Om-self.Or)+rx/a**3+rr/a**4 + def RHS_(self, y, lna): # we are solving rx, so that rhox=rx/a**3 and rhor=rr/a**4 ## - a = exp(lna) + a = np.exp(lna) H2 = self.H2_rxrr_a(a, y[0], y[1]) - H = sqrt(abs(H2)) + H = np.sqrt(abs(H2)) factor = self.lam*y[0]/H - return array([- factor, + factor * a]) + return np.array([- factor, + factor * a]) + def FuncMin_(self, x): fractoday, self.Or = x - self.Odm_dec = self.Odm*fractoday + self.Odm_dec = self.Odm*fractoday self.Odm_ndec = self.Odm*(1-fractoday) - self.Ocb_std = self.Ob+self.Odm_ndec - yinit = array([self.Odm_dec, self.Or]) - sol = odeint(self.RHS_, yinit, self.logar) + self.Ocb_std = self.Ob+self.Odm_ndec + yinit = np.array([self.Odm_dec, self.Or]) + sol = odeint(self.RHS_, yinit, self.logar) rxe, rre = sol[-1, :] # we want Ore early be as small as possible eps = rre**2 @@ -82,9 +85,10 @@ def FuncMin_(self, x): eps += ((rxe/(rxe+self.Odm_ndec))-self.xfrac)**2 return eps + def SolveEq(self): self.Odm = self.Ocb-self.Obh2/(self.h**2) - self.Ob = self.Ocb-self.Odm + self.Ob = self.Ocb-self.Odm if (self.lam == 0): self.fractoday, self.Or = self.xfrac, 0.0 if (self.xfrac == 1.0): @@ -102,30 +106,33 @@ def SolveEq(self): # stupid interp1d doesn't take inverses self.Ocb_std = self.Ob+self.Odm*(1-self.fractoday) self.Odm_dec = self.Odm*self.fractoday - yinit = array([self.Odm_dec, self.Or]) - sol = odeint(self.RHS_, yinit, self.logar) + yinit = np.array([self.Odm_dec, self.Or]) + sol = odeint(self.RHS_, yinit, self.logar) + self.sol = sol - self.rx = interp1d(self.ilogar, sol[::-1, 0]) - self.rr = interp1d(self.ilogar, sol[::-1, 1]) + self.rx = interp1d(self.ilogar, sol[::-1, 0]) + self.rr = interp1d(self.ilogar, sol[::-1, 1]) # take early time solution self.Ocbh2_early = (self.Ocb_std+sol[-1, 0])*self.h**2 + def RHSquared_a(self, a): - lna = log(a) + lna = np.log(a) return self.H2_rxrr_a(a, self.rx(lna), self.rr(lna)) + def WangWangVec(self): print("no WW with Decay") - stop() + sys.exit(1) return None # this returns the "SimpleCMB" variables in a vec def CMBSimpleVec(self): zstar = 1090. Dastar = self.Da_z(zstar)*self.c_/(self.h*100) - return array([self.Obh2, self.Ocbh2_early, Dastar/self.rd]) + return np.array([self.Obh2, self.Ocbh2_early, Dastar/self.rd]) def Om_z(self, a): - lna = log(a) + lna = np.log(a) return self.Ocb_std/a**3 + self.rx(lna)/a**3 diff --git a/py/EarlyDECosmology.py b/Models/EarlyDECosmology.py similarity index 83% rename from py/EarlyDECosmology.py rename to Models/EarlyDECosmology.py index 4300963..0ac2bda 100644 --- a/py/EarlyDECosmology.py +++ b/Models/EarlyDECosmology.py @@ -2,8 +2,9 @@ # Early Dark Energy cosmology. # -import math as N -from LCDMCosmology import * + +from LCDMCosmology import LCDMCosmology +from ParamDefs import w_par, Ode_par class EarlyDECosmology(LCDMCosmology): @@ -13,24 +14,24 @@ def __init__(self, varyw=True, varyOde=True, userd_DE=True): self.userd_DE = userd_DE print('userd', userd_DE) - self.varyw = varyw + self.varyw = varyw self.varyOde = varyOde - self.w0 = w_par.value + self.w0 = w_par.value self.Ode = Ode_par.value self.oC = LCDMCosmology() LCDMCosmology.__init__(self) + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): l = LCDMCosmology.freeParameters(self) - if (self.varyw): - l.append(w_par) - if (self.varyOde): - l.append(Ode_par) + if (self.varyw): l.append(w_par) + if (self.varyOde): l.append(Ode_par) return l + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) self.oC.updateParams(pars) @@ -43,6 +44,7 @@ def updateParams(self, pars): self.Ode = p.value return True + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): @@ -53,6 +55,7 @@ def RHSquared_a(self, a): factor*a**(3*(1.0+self.w0))) + self.Ode*(1.0-a**(-3*self.w0)) return factor/(1.0-Omega_d) + def prefactor(self): if self.userd_DE: self.rd = self.oC.rd*((1.-self.Ode)**(0.5)) @@ -61,10 +64,10 @@ def prefactor(self): return self.c_/(self.rd*self.h*100) -# for printing purposes only +# for printing purposes only def Omega_de(self, a): NuContrib = self.NuDensity.rho(a)/self.h**2 - Omega_d0 = 1.0-self.Ocb-self.Omrad-self.NuDensity.rho(1.0)/self.h**2 - factor = self.Ocb/a**3+self.Omrad/a**4+NuContrib + Omega_d0 = 1.0-self.Ocb-self.Omrad-self.NuDensity.rho(1.0)/self.h**2 + factor = self.Ocb/a**3+self.Omrad/a**4+NuContrib return (Omega_d0-self.Ode*(1.0-a**(-3*self.w0)))/(Omega_d0+factor*a**(3*(1.0+self.w0))) + self.Ode*(1.0-a**(-3*self.w0)) diff --git a/py/JordiCDMCosmology.py b/Models/JordiCDMCosmology.py similarity index 76% rename from py/JordiCDMCosmology.py rename to Models/JordiCDMCosmology.py index 4a5b718..7996cbf 100644 --- a/py/JordiCDMCosmology.py +++ b/Models/JordiCDMCosmology.py @@ -3,47 +3,45 @@ # The world is full of mystery. # - +import sys import math as N -from LCDMCosmology import * +from LCDMCosmology import LCDMCosmology +from ParamDefs import Ok_par, q_par, za_par, zb_par, wd_par, Od_par class JordiCDMCosmology(LCDMCosmology): def __init__(self): - # two parameters: Om and h + # Adding bunch of parameters. self.Ok = Ok_par.value - self.q = q_par.value + self.q = q_par.value self.za = za_par.value self.zb = zb_par.value self.wd = wd_par.value self.Od = Od_par.value LCDMCosmology.__init__(self) - # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): - return [Ok_par, q_par, za_par, zb_par, wd_par, Od_par]+LCDMCosmology.freeParameters(self) + return LCDMCosmology.freeParameters(self) + [Ok_par, q_par, za_par, zb_par, wd_par, Od_par] + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) if not ok: return False for p in pars: - if p.name == "q": - self.q = p.value - if p.name == "za": - self.za = p.value - if p.name == "zb": - self.zb = p.value - if p.name == "wd": - self.wd = p.value - if p.name == "Od": - self.Od = p.value + if p.name == "q": self.q = p.value + if p.name == "za": self.za = p.value + if p.name == "zb": self.zb = p.value + if p.name == "wd": self.wd = p.value + if p.name == "Od": self.Od = p.value elif p.name == "Ok": self.Ok = p.value self.setCurvature(self.Ok) return True + def Az(self, a): z = 1.0/a-1.0 if((self.zb < z) & (z < self.za)): @@ -53,7 +51,9 @@ def Az(self, a): elif(z < self.zb): return 1-self.q else: - stop("BAD MODEL") + print("BAD MODEL") + sys.exit(1) + def Qdrz(self, a): z = 1.0/a-1.0 @@ -66,6 +66,7 @@ def Qdrz(self, a): else: return 0.0 + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): @@ -77,6 +78,7 @@ def RHSquared_a(self, a): Omegad = self.Od*a**(-3*(1.0+self.wd)) return (self.Omrad/a**4+self.Az(a)*self.Om/a**3+self.Ok/a**2+NuContrib+Olambda+self.Qdrz(a)+Omegad) + def prior_loglike(self): if (self.za <= self.zb): return -100 diff --git a/py/LCDMCosmology.py b/Models/LCDMCosmology.py similarity index 79% rename from py/LCDMCosmology.py rename to Models/LCDMCosmology.py index 0f88d90..0e69c31 100644 --- a/py/LCDMCosmology.py +++ b/Models/LCDMCosmology.py @@ -1,23 +1,28 @@ # This is LCDM cosmology -# I didn't indend it that way, but it is used as a base class +# I didn't intend it that way, but it is used as a base class # for most other cosmologies, mostly because it treats Neutrinos and Radiation # hassle. -from BaseCosmology import * +import sys +import scipy as sp import CosmoApprox as CA +from BaseCosmology import BaseCosmology +from RadiationAndNeutrinos import RadiationAndNeutrinos +from ParamDefs import Obh2_par, Om_par, h_par, mnu_par, Nnu_par -class LCDMCosmology(BaseCosmology, RadiationAndNeutrinos): +class LCDMCosmology(BaseCosmology, RadiationAndNeutrinos): # possible options: "Anderson", "Cuesta", "CuestaNeff", "EH" rd_approx = "Cuesta" def __init__(self, Obh2=Obh2_par.value, Om=Om_par.value, h=h_par.value, mnu=mnu_par.value, Nnu=Nnu_par.value, degenerate_nu=False, disable_radiation=False, fixOm=False): # two parameters: Om and h - self.Om = Om - self.Obh2 = Obh2 + self.Om = Om + self.Obh2 = Obh2 self.fixOm = fixOm + BaseCosmology.__init__(self, h) RadiationAndNeutrinos.__init__( self, mnu, Nnu, degenerate_nu, disable=disable_radiation) @@ -32,7 +37,8 @@ def __init__(self, Obh2=Obh2_par.value, Om=Om_par.value, h=h_par.value, mnu=mnu_ elif (self.rd_approx == "EH"): self.rd_func_ = self.rd_EH_approx else: - Error("Bad rd Approx specified") + print("Bad rd Approx specified") + sys.exit(1) # Force rd update LCDMCosmology.updateParams(self, []) @@ -40,10 +46,12 @@ def __init__(self, Obh2=Obh2_par.value, Om=Om_par.value, h=h_par.value, mnu=mnu_ # by default we keep Obh2 prior self.noObh2prior = False + def setNoObh2prior(self, val=True): print("Disabling obh2 prior.") self.noObh2prior = val + # my free parameters, note that h is defined in Base # to change parameters/priors see ParamDefs.py def freeParameters(self): @@ -59,6 +67,7 @@ def freeParameters(self): l += RadiationAndNeutrinos.freeParameters(self) return l + def updateParams(self, pars): BaseCosmology.updateParams(self, pars) RadiationAndNeutrinos.updateParams(self, pars) @@ -72,7 +81,7 @@ def updateParams(self, pars): if (not self.varyPrefactor): Nnu = self.Nnu() # if we have disabled radiation, let use just use std value - # to preven rd calculator from crashing. Not ever realistically + # to prevent rd calculator from crashing. Not ever realistically # used if (Nnu == 0): Nnu = 3.04 @@ -80,6 +89,7 @@ def updateParams(self, pars): self.Obh2, self.Ocb*self.h**2, self.Omnuh2, Nnu)) return True + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 # @autojit @@ -87,6 +97,7 @@ def RHSquared_a(self, a): NuContrib = self.NuDensity.rho(a)/self.h**2 return (self.Ocb/a**3+self.Omrad/a**4+NuContrib+(1.0-self.Om)) + # Obh2 prior def prior_loglike(self): if (self.varyPrefactor or self.noObh2prior): @@ -95,24 +106,24 @@ def prior_loglike(self): # Cooke et al, http://arxiv.org/abs/1308.3240 # 2.202 +/- 0.046 #return -(self.Obh2-0.02202)**2/(2*0.00046**2) - return 0 + # this returns the Wang+Wang variables in a vec def WangWangVec(self): - Omh2 = self.Ocb*self.h**2+self.Omnuh2 - omt = Omh2/self.h**2 - zstar = self.CA.z_lastscattering(Omh2, self.Obh2) + Omh2 = self.Ocb*self.h**2+self.Omnuh2 + omt = Omh2/self.h**2 + zstar = self.CA.z_lastscattering(Omh2, self.Obh2) Dastar = self.Da_z(zstar)*self.c_/(self.h*100) - R = sqrt(omt)*self.h*100*Dastar/self.c_ - la = pi*Dastar/self.CA.soundhorizon_star(Omh2, self.Obh2) + R = sp.sqrt(omt)*self.h*100*Dastar/self.c_ + la = sp.pi*Dastar/self.CA.soundhorizon_star(Omh2, self.Obh2) # print la, R, self.Obh2 - # stop() - return array([la, R, self.Obh2]) + return sp.array([la, R, self.Obh2]) + # this returns the "SimpleCMB" variables in a vec def CMBSimpleVec(self): - Ocbh2 = self.Ocb*self.h**2 - zstar = 1090 + Ocbh2 = self.Ocb*self.h**2 + zstar = 1090 Dastar = self.Da_z(zstar)*self.c_/(self.h*100) - return array([self.Obh2, Ocbh2, Dastar/self.rd]) + return sp.array([self.Obh2, Ocbh2, Dastar/self.rd]) diff --git a/py/PolyCDMCosmology.py b/Models/PolyCDMCosmology.py similarity index 82% rename from py/PolyCDMCosmology.py rename to Models/PolyCDMCosmology.py index 1619f68..7892aaf 100644 --- a/py/PolyCDMCosmology.py +++ b/Models/PolyCDMCosmology.py @@ -1,8 +1,8 @@ # This is LCDM cosmology with optional -# curvature which you can set up with -# setVaryOk() +# free parameters on the Hubble function -from LCDMCosmology import * +from LCDMCosmology import LCDMCosmology +from ParamDefs import Ok_par, Om1_par, Om2_par class PolyCDMCosmology(LCDMCosmology): @@ -10,24 +10,26 @@ def __init__(self, polyvary=['Om1','Om2','Ok'], Ok_prior=0.1): # Ok, LCDM has Omega_m, we also have Omega_1 and Omega_2 # and Lambda is then what remains ## - self.Ok = Ok_par.value + self.Ok = Ok_par.value self.Om1 = Om1_par.value self.Om2 = Om2_par.value - self.varyOm1 = 'Om1' in polyvary - self.varyOm2 = 'Om2' in polyvary - self.varyOk = 'Ok' in polyvary - self.Ok_prior=Ok_prior + self.varyOm1 = 'Om1' in polyvary + self.varyOm2 = 'Om2' in polyvary + self.varyOk = 'Ok' in polyvary + self.Ok_prior = Ok_prior LCDMCosmology.__init__(self, mnu=0) + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): l = LCDMCosmology.freeParameters(self) if self.varyOm1: l.append(Om1_par) if self.varyOm2: l.append(Om2_par) Ok_par.setError(0.7) - if self.varyOk: l.append(Ok_par) + if self.varyOk: l.append(Ok_par) return l + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) if not ok: @@ -44,6 +46,7 @@ def updateParams(self, pars): return False return True + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): @@ -51,6 +54,7 @@ def RHSquared_a(self, a): # as 3/23 -- note Om1 and Om2 used to be swapped return (self.Om/a**3+self.Om2/a**2+self.Ok/a**2+self.Om1/a+(1-self.Om-self.Om1-self.Om2-self.Ok)) + def prior_loglike(self): return (-self.Ok**2/(2*self.Ok_prior**2) # A 0.1 prior in Ok as discussed at the telecon + LCDMCosmology.prior_loglike(self)) diff --git a/py/QuintCosmology.py b/Models/QuintCosmology.py similarity index 100% rename from py/QuintCosmology.py rename to Models/QuintCosmology.py diff --git a/py/SlowRDECosmology.py b/Models/SlowRDECosmology.py similarity index 82% rename from py/SlowRDECosmology.py rename to Models/SlowRDECosmology.py index 187736d..62253d1 100644 --- a/py/SlowRDECosmology.py +++ b/Models/SlowRDECosmology.py @@ -3,31 +3,31 @@ # I do not. # -import math as N -from LCDMCosmology import * +from LCDMCosmology import LCDMCosmology +from ParamDefs import Ok_par, dw_par class SlowRDECosmology(LCDMCosmology): def __init__(self, varyw=True, varyOk=True): # two parameters: Om and h - self.varyw = varyw + self.varyw = varyw self.varyOk = varyOk - self.Ok = Ok_par.value + self.Ok = Ok_par.value self.dw0 = dw_par.value LCDMCosmology.__init__(self) + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): l = LCDMCosmology.freeParameters(self) - if (self.varyw): - l.append(dw_par) - if (self.varyOk): - l.append(Ok_par) + if (self.varyw): l.append(dw_par) + if (self.varyOk): l.append(Ok_par) return l + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) if not ok: @@ -42,11 +42,12 @@ def updateParams(self, pars): return False return True + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): onepz = 1.0/a NuContrib = self.NuDensity.rho(a)/self.h**2 - Ode = 1.0-self.Om-self.Ok + Ode = 1.0-self.Om-self.Ok rhow = (onepz**3.0/(self.Om*onepz**3.0+Ode))**(self.dw0/Ode) return (self.Ocb/a**3+self.Ok/a**2+self.Omrad/a**4+NuContrib+Ode*rhow) diff --git a/py/SplineLCDMCosmology.py b/Models/SplineLCDMCosmology.py similarity index 80% rename from py/SplineLCDMCosmology.py rename to Models/SplineLCDMCosmology.py index 115bb48..1c9b2f1 100644 --- a/py/SplineLCDMCosmology.py +++ b/Models/SplineLCDMCosmology.py @@ -2,9 +2,10 @@ # This si a cosmology where w(z) is defined by splines. ## -from numpy import linspace + from scipy.interpolate import InterpolatedUnivariateSpline -from LCDMCosmology import * +from LCDMCosmology import LCDMCosmology +from ParamDefs import Sp1_par, Sp2_par, Sp3_par, Sp4_par import math as N @@ -22,28 +23,27 @@ def __init__(self, varySp1=True, varySp2=True, varySp3=True, varySp4=True): self.Sp4 = Sp4_par.value # Nodes are equally-spaced in log10(z) - self.zmin = 0.1 # first-node position - self.zmax = 2.5 # last-node position + self.zmin = 0.1 # first-node position + self.zmax = 2.5 # last-node position self.Nnodes = 6 # number of nodes used - self.lzmin = log10(self.zmin) - self.lzmax = log10(self.zmax) + self.lzmin = N.log10(self.zmin) + self.lzmax = N.log10(self.zmax) LCDMCosmology.__init__(self) + + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): l = LCDMCosmology.freeParameters(self) - if (self.varySp1): - l.append(Sp1_par) - if (self.varySp2): - l.append(Sp2_par) - if (self.varySp3): - l.append(Sp3_par) - if (self.varySp4): - l.append(Sp4_par) + if (self.varySp1): l.append(Sp1_par) + if (self.varySp2): l.append(Sp2_par) + if (self.varySp3): l.append(Sp3_par) + if (self.varySp4): l.append(Sp4_par) return l + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) if not ok: @@ -59,16 +59,18 @@ def updateParams(self, pars): self.Sp4 = p.value return True + def Spline(self, a): z = 1.0/a-1.0 x = [10**(self.lzmin+(self.lzmax-self.lzmin)/(self.Nnodes-1)*i) for i in range(0, int(self.Nnodes))] - # print x + y = [1.0, self.Sp1, self.Sp2, self.Sp3, self.Sp4, 1.0] s = InterpolatedUnivariateSpline(x, y) ys = s(z) return ys + def Rho_de(self, a): z = 1.0/a-1.0 if (z > self.zmax or z < self.zmin): @@ -76,6 +78,7 @@ def Rho_de(self, a): else: return self.Spline(a) + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): diff --git a/py/StepCDMCosmology.py b/Models/StepCDMCosmology.py similarity index 81% rename from py/StepCDMCosmology.py rename to Models/StepCDMCosmology.py index fd264f3..04eb5a6 100644 --- a/py/StepCDMCosmology.py +++ b/Models/StepCDMCosmology.py @@ -3,18 +3,16 @@ # import numpy -from numpy import linspace -#from scipy.interpolate import InterpolatedUnivariateSpline -from LCDMCosmology import * -import math as N +from LCDMCosmology import LCDMCosmology +from ParamDefs import * class StepCDMCosmology(LCDMCosmology): def __init__(self): self.NZ = int(step_nz_par.value) - self.Z = numpy.zeros((self.NZ)) # redshifts of bin boundaries - self.R = numpy.zeros((self.NZ+1)) # rho_de/rho_c in bins + self.Z = numpy.zeros((self.NZ)) # redshifts of bin boundaries + self.R = numpy.zeros((self.NZ+1)) # rho_de/rho_c in bins # see ParamDefs.py for the values of parameters if(self.NZ > 0): @@ -41,23 +39,20 @@ def __init__(self): self.Ok = 1.-self.RHSquared_a(1.) self.setCurvature(self.Ok) - # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): l = LCDMCosmology.freeParameters(self) if(self.NZ > 0): l.append(step_rho0_par) l.append(step_rho1_par) - if(self.NZ > 1): - l.append(step_rho2_par) - if(self.NZ > 2): - l.append(step_rho3_par) - if(self.NZ > 3): - l.append(step_rho4_par) - if(self.NZ > 4): - l.append(step_rho5_par) + if(self.NZ > 1): l.append(step_rho2_par) + if(self.NZ > 2): l.append(step_rho3_par) + if(self.NZ > 3): l.append(step_rho4_par) + if(self.NZ > 4): l.append(step_rho5_par) return l + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) if not ok: @@ -73,6 +68,7 @@ def updateParams(self, pars): self.setCurvature(self.Ok) return True + def Rho_de(self, a): z = 1.0/a-1.0 for i in range(self.NZ): @@ -80,6 +76,7 @@ def Rho_de(self, a): return self.R[i] return self.R[self.NZ] + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): diff --git a/py/TiredLightDecorator.py b/Models/TiredLightDecorator.py similarity index 75% rename from py/TiredLightDecorator.py rename to Models/TiredLightDecorator.py index 3503cc8..70305ad 100644 --- a/py/TiredLightDecorator.py +++ b/Models/TiredLightDecorator.py @@ -2,8 +2,8 @@ # This helper class adds a "tired light" parameter # to any class ## -from ParamDefs import * -from scipy import * +from ParamDefs import beta_par +import math as N def TiredLightDecorator(LikeInstance): @@ -15,18 +15,20 @@ def __init__(self, LikeInstance): self.LikeType = LikeInstance.__class__ self.beta = beta_par.value + def freeParameters(self): - return [beta_par]+self.LikeType.freeParameters(self) + return self.LikeType.freeParameters(self) + [beta_par] + def updateParams(self, pars): for p in pars: - if p.name == "beta": - self.beta = p.value + if p.name == "beta": self.beta = p.value self.LikeType.updateParams(self, pars) + # distance modulus def distance_modulus(self, z): assert(not self.varyPrefactor) - return 5*log10(self.Da_z(z)*(1+z)**(1+self.beta)) + return 5*N.log10(self.Da_z(z)*(1+z)**(1+self.beta)) return TiredLightLikelihood(LikeInstance) diff --git a/py/WeirdCDMCosmology.py b/Models/WeirdCDMCosmology.py similarity index 73% rename from py/WeirdCDMCosmology.py rename to Models/WeirdCDMCosmology.py index a6e5902..7fc5a3d 100644 --- a/py/WeirdCDMCosmology.py +++ b/Models/WeirdCDMCosmology.py @@ -3,49 +3,46 @@ import math as N - -from LCDMCosmology import * +from LCDMCosmology import LCDMCosmology +from ParamDefs import Ok_par, mu_par, Amp_par, sig_par class WeirdCDMCosmology(LCDMCosmology): def __init__(self, varymu=True, varyAmp=True, varysig=True, varyCos=True): - # two parameters: Om and h + # several parameters: # we start with false here... varyOk = False # this is my "original cosmology" -- outside gaussian not much will change. - self.varyOk = varyOk - self.varymu = varymu + self.varyOk = varyOk + self.varymu = varymu self.varyAmp = varyAmp self.varysig = varysig self.varyCos = varyCos - self.Ok = Ok_par.value - self.mu = mu_par.value + self.Ok = Ok_par.value + self.mu = mu_par.value self.Amp = Amp_par.value self.sig = sig_par.value # auxiliary self cosmology self.oC = LCDMCosmology() - LCDMCosmology.__init__(self) + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): l = [] if (self.varyCos): l += LCDMCosmology.freeParameters(self) - if (self.varyOk): - l.append(Ok_par) - if (self.varymu): - l.append(mu_par) - if (self.varyAmp): - l.append(Amp_par) - if (self.varysig): - l.append(sig_par) + if (self.varyOk): l.append(Ok_par) + if (self.varymu): l.append(mu_par) + if (self.varyAmp): l.append(Amp_par) + if (self.varysig): l.append(sig_par) return l + def updateParams(self, pars): LCDMCosmology.updateParams(self, pars) self.oC.updateParams(pars) @@ -64,24 +61,27 @@ def updateParams(self, pars): self.sig = p.value return True + def Weird(self, a): # return W(z) and dW(z)/dz - z = 1.0/a-1 - lmu = log(self.mu+1) - lz = log(z+1) - W = self.Amp*N.exp(-(lz-lmu)**2/(2*self.sig**2)) - WP = W*(-(lz-lmu)/self.sig**2)*a # a=1/(1+z) = dlog(1+z)/dz + z = 1.0/a-1 + lmu = N.log(self.mu+1) + lz = N.log(z+1) + W = self.Amp*N.exp(-(lz-lmu)**2/(2*self.sig**2)) + WP = W*(-(lz-lmu)/self.sig**2)*a # a=1/(1+z) = dlog(1+z)/dz return W, WP + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): W, WP = self.Weird(a) H2 = LCDMCosmology.RHSquared_a(self, a) - HI = 1/sqrt(H2) + HI = 1/N.sqrt(H2) HI = HI*(1+W) + self.oC.Da_z(1/a-1)*WP return 1/HI**2 + def Da_z(self, z): W, WP = self.Weird(1/(1.0+z)) return self.oC.Da_z(z)*(1+W) diff --git a/py/oLCDMCosmology.py b/Models/oLCDMCosmology.py similarity index 83% rename from py/oLCDMCosmology.py rename to Models/oLCDMCosmology.py index c35fc19..27e7c57 100644 --- a/py/oLCDMCosmology.py +++ b/Models/oLCDMCosmology.py @@ -2,17 +2,18 @@ # curvature which you can set up with # setVaryOk() -from LCDMCosmology import * - +from LCDMCosmology import LCDMCosmology +from ParamDefs import Ok_par class oLCDMCosmology(LCDMCosmology): # zeroDE forces Ol to zero. def __init__(self, zeroDE=False, kwargs_LCDM={}): # two parameters: Om and h - self.Ok = Ok_par.value + self.Ok = Ok_par.value self.zeroDE = zeroDE LCDMCosmology.__init__(self, **kwargs_LCDM) + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): Ok_par.setValue(self.Ok) @@ -21,6 +22,7 @@ def freeParameters(self): l.append(Ok_par) return l+LCDMCosmology.freeParameters(self) + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) if not ok: @@ -28,10 +30,10 @@ def updateParams(self, pars): for p in pars: if p.name == "Ok": self.Ok = p.value - # se comments to base class + # se comments to BaseCosmology class # in short, WE hold the Ok parameter # (maybe you want to parameterize it differently), - # but the base class needs to know about Ok in order + # but the BaseCosmology class needs to know about Ok in order # to get Da and its sin/sinhs right. self.setCurvature(self.Ok) if (abs(self.Ok) > 1.0): @@ -42,7 +44,8 @@ def updateParams(self, pars): self.setCurvature(self.Ok) return True - # this is relative hsquared as a function of a + + # this is relative Hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): NuContrib = self.NuDensity.rho(a)/self.h**2 diff --git a/py/owa0CDMCosmology.py b/Models/owa0CDMCosmology.py similarity index 82% rename from py/owa0CDMCosmology.py rename to Models/owa0CDMCosmology.py index 40bbd42..85e151f 100644 --- a/py/owa0CDMCosmology.py +++ b/Models/owa0CDMCosmology.py @@ -2,14 +2,14 @@ import math as N -from LCDMCosmology import * - +from LCDMCosmology import LCDMCosmology +from ParamDefs import w_par, wa_par, Ok_par class owa0CDMCosmology(LCDMCosmology): def __init__(self, varyw=True, varywa=True, varyOk=True): - # two parameters: Om and h + # three parameters: w, wa, Ok - self.varyw = varyw + self.varyw = varyw self.varywa = varywa self.varyOk = varyOk @@ -18,17 +18,16 @@ def __init__(self, varyw=True, varywa=True, varyOk=True): self.wa = wa_par.value LCDMCosmology.__init__(self) + # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) def freeParameters(self): l = LCDMCosmology.freeParameters(self) - if (self.varyw): - l.append(w_par) - if (self.varywa): - l.append(wa_par) - if (self.varyOk): - l.append(Ok_par) + if (self.varyw): l.append(w_par) + if (self.varywa): l.append(wa_par) + if (self.varyOk): l.append(Ok_par) return l + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) if not ok: @@ -45,6 +44,7 @@ def updateParams(self, pars): return False return True + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): diff --git a/py/wLCDMCosmology.py b/Models/wCDMCosmology.py similarity index 58% rename from py/wLCDMCosmology.py rename to Models/wCDMCosmology.py index e295e8f..a5a4651 100644 --- a/py/wLCDMCosmology.py +++ b/Models/wCDMCosmology.py @@ -1,17 +1,24 @@ -# This is a CDM cosmology with w +# This is a CDM cosmology with constant eos w for DE -from LCDMCosmology import * +from LCDMCosmology import LCDMCosmology +from ParamDefs import w_par -class wLCDMCosmology(LCDMCosmology): - def __init__(self): +class wCDMCosmology(LCDMCosmology): + def __init__(self, varyw=True): # two parameters: Om and h - self.w = -1 + + self.varyw = varyw + self.w = w_par.value LCDMCosmology.__init__(self) - # my free parameters. We add Ok on top of LCDM ones (we inherit LCDM) + + # my free parameters. We add w on top of LCDM ones (we inherit LCDM) def freeParameters(self): - return [w_par]+LCDMCosmology.freeParameters(self) + l = LCDMCosmology.freeParameters(self) + if (self.varyw): l.append(w_par) + return l + def updateParams(self, pars): ok = LCDMCosmology.updateParams(self, pars) @@ -22,6 +29,7 @@ def updateParams(self, pars): self.w = p.value return True + # this is relative hsquared as a function of a ## i.e. H(z)^2/H(z=0)^2 def RHSquared_a(self, a): diff --git a/README.md b/README.md index a244b5e..8bb0339 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -april +april. ===== A simple MCMC code for cosmological parameter estimation where only @@ -113,7 +113,7 @@ L=CompositeLikelihood([ DR11LyaCross() ]) -and then connect it to a theory, like +and a then connect it to a theory, like L.setTheory(oLCDMCosmology()) This is now a nice package: diff --git a/Run/RunBase.py b/Run/RunBase.py index 96ea655..7d0604a 100644 --- a/Run/RunBase.py +++ b/Run/RunBase.py @@ -1,43 +1,51 @@ +# coding=utf-8 # Add paths, we want to be able to run in either root or Run/ -import sys -import os -# print sys.path -sys.path = ["py", "../py"]+sys.path -from CosmoMCImportanceSampler import * -from MCMCAnalyzer import * -from MaxLikeAnalyzer import * -from LikelihoodMultiplier import * -from CompositeLikelihood import * -from HubbleParameterLikelihood import * -from CompressedSNLikelihood import * -from SimpleCMB import * -from BAOLikelihoods import * -from QuintCosmology import * -from SlowRDECosmology import * -from EarlyDECosmology import * -from StepCDMCosmology import * -from DecayLCDMCosmology import * -from SplineLCDMCosmology import * -from TiredLightDecorator import * -from WeirdCDMCosmology import * -from JordiCDMCosmology import * -from owa0CDMCosmology import * -from PolyCDMCosmology import * -from wLCDMCosmology import * -from oLCDMCosmology import * -from LCDMCosmology import * - -# Cosmologies -#from wDMCosmology import * - -# Like modules +import sys +sys.path = ["py", "../py", "Models", "Cosmo", "Likelihoods"] + sys.path + +#TODO -- Include model used in several papers +#TODO -- Add Choronometers data +#TODO -- Add Planck 15 +#TODO -- Add DR12 Galaxies + +# Cosmologies already included +from LCDMCosmology import LCDMCosmology +from oLCDMCosmology import oLCDMCosmology +from wCDMCosmology import wCDMCosmology +from owa0CDMCosmology import owa0CDMCosmology +from PolyCDMCosmology import PolyCDMCosmology +from JordiCDMCosmology import JordiCDMCosmology +from WeirdCDMCosmology import WeirdCDMCosmology +from TiredLightDecorator import TiredLightDecorator +from SplineLCDMCosmology import SplineLCDMCosmology +from DecayLCDMCosmology import DecayLCDMCosmology +from StepCDMCosmology import StepCDMCosmology +from EarlyDECosmology import EarlyDECosmology +from SlowRDECosmology import SlowRDECosmology +from BinnedWCosmology import BinnedWCosmology +from QuintCosmology import QuintCosmology # Composite Likelihood +from CompositeLikelihood import CompositeLikelihood # Likelihood Multiplier +from LikelihoodMultiplier import LikelihoodMultiplier + +# Likelihood modules +from BAOLikelihoods import DR11LOWZ, DR11CMASS, DR14LyaAuto, DR14LyaCross, \ + SixdFGS, SDSSMGS, DR11LyaAuto, DR11LyaCross +from SimpleCMB import PlanckLikelihood, WMAP9Likelihood +from CompressedSNLikelihood import * +from HubbleParameterLikelihood import * # Analyzers +from MCMCAnalyzer import * +from MaxLikeAnalyzer import * + +#Importance Sampling +from CosmoMCImportanceSampler import * + # String parser Aux routines model_list = "LCDOM, LCDMasslessnu, nuLCDM, NeffLCDM, noradLCDM, nuoLCDM, nuwLCDM, oLCDM, wCDM, waCDM, owCDM,"\ @@ -46,6 +54,18 @@ def ParseModel(model): + """  + Parameters + ----------- + model: + name of the model, i.e. LCDM + + Returns + ----------- + object - info/calculations based on this model: i.e. d_L, d_A, d_H + + """ + if model == "LCDM": T = LCDMCosmology() elif model == "LCDMmasslessnu": @@ -65,9 +85,9 @@ def ParseModel(model): T = oLCDMCosmology() T.setVaryMnu() elif model == "wCDM": - T = wLCDMCosmology() + T = wCDMCosmology() elif model == "nuwCDM": - T = wLCDMCosmology() + T = wCDMCosmology() T.setVaryMnu() elif model == "waCDM": T = owa0CDMCosmology(varyOk=False) @@ -96,21 +116,21 @@ def ParseModel(model): elif model == "PolyCDM": T = PolyCDMCosmology() elif model == "fPolyCDM": - T = PolyCDMCosmology(polyvary=['Om1','Om2']) + T = PolyCDMCosmology(polyvary=['Om1', 'Om2']) elif model == "PolyOk": ## polycdm for OK T = PolyCDMCosmology(Ok_prior=10.) elif model == "PolyOkc": ## polycdm sans Om2 term to couple two - T = PolyCDMCosmology(polyvary=['Om1','Ok'],Ok_prior=10.) + T = PolyCDMCosmology(polyvary=['Om1', 'Ok'], Ok_prior=10.) elif model == "PolyOkf": ## polycdm sans Om2 term to couple two - T = PolyCDMCosmology(polyvary=['Om1','Om2']) - - + T = PolyCDMCosmology(polyvary=['Om1', 'Om2']) elif model == "EarlyDE": T = EarlyDECosmology(varyw=False, userd_DE=False) elif model == "EarlyDE_rd_DE": T = EarlyDECosmology(varyw=False) elif model == "SlowRDE": T = SlowRDECosmology(varyOk=False) + elif model == "Binned": + T = BinnedWCosmology() elif model == "Quint_last": T = QuintCosmology() elif model == 'wDM': @@ -127,6 +147,17 @@ def ParseModel(model): def ParseDataset(datasets): + """  + Parameters + ----------- + datasets: + name of datasets, i.e. BBAO + + Returns + ----------- + object - likelihood + + """ dlist = datasets.split('+') L = CompositeLikelihood([]) for name in dlist: @@ -148,9 +179,9 @@ def ParseDataset(datasets): ]) elif name == 'GBAOx10': L.addLikelihoods([ - LikelihoodMultiplier(DR11LOWZ(), 100.0), + LikelihoodMultiplier(DR11LOWZ(), 100.0), LikelihoodMultiplier(DR11CMASS(), 100.0), - LikelihoodMultiplier(SixdFGS(), 100.0) + LikelihoodMultiplier(SixdFGS(), 100.0) ]) elif name == 'GBAO_no6dF': L.addLikelihoods([ diff --git a/Run/driver.py b/Run/driver.py index b883c5b..7a38342 100755 --- a/Run/driver.py +++ b/Run/driver.py @@ -1,6 +1,8 @@ #!/usr/bin/env python import sys -from RunBase import * +from RunBase import ParseModel, ParseDataset +from MCMCAnalyzer import MCMCAnalyzer +from MaxLikeAnalyzer import MaxLikeAnalyzer # # This little code allows running something @@ -9,15 +11,15 @@ # if (len(sys.argv) < 4): - print("""Usage: + print("""Usage: -Run/driver.py pre/phy model dataset [chain #] [# to skip] [# to take] [directory] +Run/driver.py pre/phy model dataset [N chain #] [# steps to skip] [# to take] [chains directory] -where model is one of +where model is one of -%s +%s -and datasets is one of +and datasets is one of %s @@ -30,6 +32,8 @@ prefact, model, datasets = sys.argv[1:4] + + if len(sys.argv) > 4: chainno = int(sys.argv[4]) else: @@ -61,7 +65,7 @@ L.setTheory(T) if chainno > 0: - M = MCMCAnalyzer(L, chainsdir + "/" + model+"_"+prefact+"_" + + M = MCMCAnalyzer(L, chainsdir + "/" + model + "_"+prefact + "_" + datasets, skip=skip, nsamp=nsamp, temp=temp, chain_num=chainno) else: A = MaxLikeAnalyzer(L) diff --git a/chains/LCDM_phy_SN.paramnames b/chains/LCDM_phy_SN.paramnames new file mode 100644 index 0000000..3343835 --- /dev/null +++ b/chains/LCDM_phy_SN.paramnames @@ -0,0 +1,5 @@ +Om \Omega_m +Obh2 \Omega_{b}h^2 +h h +BetouleSN_like BetouleSN +theory_prior None diff --git a/chains/LCDM_phy_SN_1.maxlike b/chains/LCDM_phy_SN_1.maxlike new file mode 100644 index 0000000..ce295cf --- /dev/null +++ b/chains/LCDM_phy_SN_1.maxlike @@ -0,0 +1 @@ +0.000247475 16.6084 0.295122 0.0223056 0.731594 -16.6084 0 diff --git a/chains/LCDM_phy_SN_1.txt b/chains/LCDM_phy_SN_1.txt new file mode 100644 index 0000000..f1b63d7 --- /dev/null +++ b/chains/LCDM_phy_SN_1.txt @@ -0,0 +1,38 @@ +0.00024613 16.6193 0.291577 0.0223791 0.696686 -16.6193 0 +0.000247416 16.6089 0.294842 0.022397 0.708031 -16.6089 0 +0.000247401 16.609 0.29477 0.0223988 0.715375 -16.609 0 +0.00024002 16.6696 0.284988 0.022414 0.717949 -16.6696 0 +0.000247414 16.6089 0.294825 0.0224278 0.72003 -16.6089 0 +0.000240668 16.6642 0.285506 0.0224105 0.71777 -16.6642 0 +0.000223333 16.8137 0.275451 0.0224091 0.725563 -16.8137 0 +0.000211195 16.9255 0.270394 0.0224018 0.729729 -16.9255 0 +0.000214935 16.8903 0.271866 0.0224152 0.738237 -16.8903 0 +0.000237397 16.6916 0.283082 0.0224035 0.748507 -16.6916 0 +0.000236727 16.6972 0.282635 0.0224235 0.756308 -16.6972 0 +0.000246571 16.6157 0.292404 0.0224239 0.754242 -16.6157 0 +0.00023321 16.7271 0.280473 0.0224279 0.757173 -16.7271 0 +0.000241557 16.6568 0.286257 0.022402 0.75523 -16.6568 0 +0.000246715 16.6146 0.292714 0.022379 0.752391 -16.6146 0 +0.000242783 16.6467 0.287387 0.0223796 0.753064 -16.6467 0 +0.000237857 16.6877 0.283398 0.0223723 0.744302 -16.6877 0 +0.000212031 16.9176 0.270718 0.0223617 0.735883 -16.9176 0 +0.000223301 16.814 0.275436 0.0223432 0.729648 -16.814 0 +0.000231006 16.7461 0.279241 0.0223447 0.732616 -16.7461 0 +0.000235375 16.7087 0.281768 0.0223223 0.729448 -16.7087 0 +0.000247475 16.6084 0.295122 0.0223056 0.731594 -16.6084 0 +0.000246076 16.6197 0.29148 0.0222694 0.730857 -16.6197 0 +0.000212777 16.9105 0.271009 0.022265 0.734153 -16.9105 0 +0.000211505 16.9225 0.270513 0.0222564 0.727493 -16.9225 0 +0.000163559 17.4367 0.254673 0.0222296 0.726396 -17.4367 0 +0.000109414 18.2407 0.238431 0.0222352 0.723488 -18.2407 0 +9.50413e-05 18.5224 0.233832 0.0222494 0.723817 -18.5224 0 +0.000174107 17.3117 0.25788 0.0222406 0.720943 -17.3117 0 +0.000163774 17.434 0.254737 0.0222292 0.724456 -17.434 0 +0.000108578 18.2561 0.238168 0.0222194 0.713438 -18.2561 0 +7.19561e-05 19.0789 0.225735 0.0222478 0.720108 -19.0789 0 +0.000114559 18.1488 0.240023 0.0222661 0.715567 -18.1488 0 +0.000100891 18.4029 0.235735 0.0222807 0.731686 -18.4029 0 +0.000127746 17.9309 0.244025 0.0222865 0.74648 -17.9309 0 +0.0001315 17.873 0.245146 0.0222813 0.741444 -17.873 0 +0.000137636 17.7818 0.246969 0.0222642 0.731177 -17.7818 0 +0.000103437 18.3531 0.236548 0.0222351 0.731022 -18.3531 0