Skip to content

Commit

Permalink
updated tests to newest ccl
Browse files Browse the repository at this point in the history
  • Loading branch information
@JaimeRZP
JaimeRZP committed Nov 18, 2024
1 parent bce68a6 commit f6bcf20
Show file tree
Hide file tree
Showing 6 changed files with 54 additions and 110 deletions.
16 changes: 4 additions & 12 deletions src/theory.jl
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,7 @@
function Theory(cosmology::Cosmology,
names, types, pairs,
idx, files;
Nuisances=Dict(),
res_wl=350, res_gc=1000,
int_wl="linear", int_gc="linear")
Nuisances=Dict())

nui_type = eltype(valtype(Nuisances))
if !(nui_type <: Float64) & (nui_type != Any)
Expand All @@ -22,22 +20,16 @@ function Theory(cosmology::Cosmology,
b = get(Nuisances, string(name, "_", "b"), 1.0)
nz = get(Nuisances, string(name, "_", "nz"), nz_mean)
zs = get(Nuisances, string(name, "_", "zs"), zs_mean)
dz = get(Nuisances, string(name, "_", "dz"), 0.0)
tracer = NumberCountsTracer(cosmology, zs .+ dz, nz;
b=b, res=res_gc,
nz_interpolation=int_gc)
tracer = NumberCountsTracer(cosmology, zs, nz; b=b)
elseif t_type == "galaxy_shear"
zs_mean, nz_mean = files[string("nz_", name)]
m = get(Nuisances, string(name, "_", "m"), 0.0)
IA_params = [get(Nuisances, "A_IA", 0.0),
get(Nuisances, "alpha_IA", 0.0)]
nz = get(Nuisances, string(name, "_", "nz"), nz_mean)
zs = get(Nuisances, string(name, "_", "zs"), zs_mean)
dz = get(Nuisances, string(name, "_", "dz"), 0.0)
tracer = WeakLensingTracer(cosmology, zs .+ dz, nz;
m=m, IA_params=IA_params,
res=res_wl,
nz_interpolation=int_wl)
tracer = WeakLensingTracer(cosmology, zs, nz;
m=m, IA_params=IA_params)

elseif t_type == "cmb_convergence"
tracer = CMBLensingTracer(cosmology)
Expand Down
74 changes: 19 additions & 55 deletions src/tracers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,16 +27,11 @@ Kwargs:
Returns:
- `NumberCountsTracer::NumberCountsTracer` : Number counts tracer structure.
"""
NumberCountsTracer(cosmo::Cosmology, z_n, nz;
b=1.0, res=1000, nz_interpolation="linear", smooth=0) = begin
z_w, nz_w = nz_interpolate(z_n, nz, res;
mode=nz_interpolation)
nz_norm = integrate(z_w, nz_w, SimpsonEven())

chi = cosmo.chi(z_w)
hz = Hmpc(cosmo, z_w)

w_arr = @. (nz_w*hz/nz_norm)
NumberCountsTracer(cosmo::Cosmology, z, nz; b=1.0) = begin
nz_norm = integrate(z, nz, SimpsonEven())
chi = cosmo.chi(z)
hz = Hmpc(cosmo, z)
w_arr = @. (nz*hz/nz_norm)
wint = linear_interpolation(chi, b .* w_arr, extrapolation_bc=0.0)
F::Function = ℓ -> 1
NumberCountsTracer(wint, F)
Expand All @@ -61,40 +56,29 @@ struct WeakLensingTracer <: Tracer
F::Function
end

WeakLensingTracer(cosmo::Cosmology, z_n, nz;
IA_params = [0.0, 0.0], m=0.0,
res=350, nz_interpolation="linear") = begin
WeakLensingTracer(cosmo::Cosmology, z, nz;
IA_params = [0.0, 0.0], m=0.0) = begin
cosmo_type = cosmo.settings.cosmo_type
z_w, nz_w = nz_interpolate(z_n, nz, res;
mode=nz_interpolation)
res = length(z_w)
nz_norm = integrate(z_w, nz_w, SimpsonEven())
chi = cosmo.chi(z_w)

# Calculate chis at which to precalculate the lensing kernel
# OPT: perhaps we don't need to sample the lensing kernel
# at all zs.
# Calculate integral at each chi
w_itg(chii) = @.(nz_w*(1-chii/chi))
res = length(z)
nz_norm = integrate(z, nz, SimpsonEven())
chi = cosmo.chi(z)
# Compute kernels
w_itg(chii) = @.(nz*(1-chii/chi))
w_arr = zeros(cosmo_type, res)
@inbounds for i in 1:res-3
w_arr[i] = integrate(z_w[i:res], w_itg(chi[i])[i:res], SimpsonEven())
w_arr[i] = integrate(z[i:res], w_itg(chi[i])[i:res], SimpsonEven())
end

# Normalize
H0 = cosmo.cpar.h/CLIGHT_HMPC
lens_prefac = 1.5*cosmo.cpar.Ωm*H0^2
# Fix first element
chi[1] = 0.0
w_arr = @. w_arr * chi * lens_prefac * (1+z_w) / nz_norm

w_arr = @. w_arr * chi * lens_prefac * (1+z) / nz_norm
if IA_params != [0.0, 0.0]
hz = Hmpc(cosmo, z_w)
As = get_IA(cosmo, z_w, IA_params)
corr = @. As * (nz_w * hz / nz_norm)
hz = Hmpc(cosmo, z)
As = get_IA(cosmo, z, IA_params)
corr = @. As * (nz * hz / nz_norm)
w_arr = @. w_arr - corr
end

# Interpolate
b = m+1.0
wint = linear_interpolation(chi, b.*w_arr, extrapolation_bc=0.0)
Expand Down Expand Up @@ -127,13 +111,12 @@ Returns:
CMBLensingTracer(cosmo::Cosmology; res=350) = begin
# chi array
chis = range(0.0, stop=cosmo.chi_max, length=res)
zs = cosmo.z_of_chi(chis)
z = cosmo.z_of_chi(chis)
# Prefactor
H0 = cosmo.cpar.h/CLIGHT_HMPC
lens_prefac = 1.5*cosmo.cpar.Ωm*H0^2
# Kernel
w_arr = @. lens_prefac*chis*(1-chis/cosmo.chi_LSS)*(1+zs)

w_arr = @. lens_prefac*chis*(1-chis/cosmo.chi_LSS)*(1+z)
# Interpolate
wint = linear_interpolation(chis, w_arr, extrapolation_bc=0.0)
F::Function = ℓ -> @.(((ℓ+1)*ℓ)/(ℓ+0.5)^2)
Expand All @@ -156,22 +139,3 @@ function get_IA(cosmo::Cosmology, zs, IA_params)
return @. A_IA*((1 + zs)/1.62)^alpha_IA * (0.0134 * cosmo.cpar.Ωm / cosmo.Dz(zs))
end

function nz_interpolate(z, nz, res; mode="linear")
if mode!="none"
if mode=="linear"
nz_int = linear_interpolation(z, nz;
extrapolation_bc=0.0)
end
if mode=="cubic"
dz = mean(z[2:end] - z[1:end-1])
z_range = z[1]:dz:z[end]
nz_int = cubic_spline_interpolation(z_range, nz;
extrapolation_bc=0.0)
end
zz_range = range(0.00001, stop=z[end], length=res)
nzz = nz_int(zz_range)
return zz_range, nzz
else
return z, nz
end
end
42 changes: 21 additions & 21 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ if test_main
end

@testset "Traces" begin
z = Vector(range(0., stop=2., length=1000))
z = Vector(range(0.01, stop=2., length=1000))
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
clustering_t = NumberCountsTracer(cosmo_EisHu, z, nz)
lensing_t = WeakLensingTracer(cosmo_EisHu, z, nz)
Expand Down Expand Up @@ -178,7 +178,7 @@ if test_main
else
ℓs = [10.0, 30.0, 100.0, 300.0, 1000.0]
end
z = Vector(range(0., stop=2., length=1000))
z = Vector(range(0.01, stop=2., length=1000))
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
tg = NumberCountsTracer(cosmo_EisHu, z, nz; b=1.0)
ts = WeakLensingTracer(cosmo_EisHu, z, nz;
Expand Down Expand Up @@ -220,7 +220,7 @@ if test_main
else
ℓs = [10.0, 30.0, 100.0, 300.0, 1000.0]
end
z = Vector(range(0., stop=2., length=1000))
z = Vector(range(0.01, stop=2., length=1000))
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
tg = NumberCountsTracer(cosmo_EisHu, z, nz; b=1.0)
ts = WeakLensingTracer(cosmo_EisHu, z, nz;
Expand Down Expand Up @@ -262,7 +262,7 @@ if test_main
else
ℓs = [10.0, 30.0, 100.0, 300.0, 1000.0]
end
z = Vector(range(0., stop=2., length=1000))
z = Vector(range(0.01, stop=2., length=1000))
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
tg = NumberCountsTracer(cosmo_emul, z, nz; b=1.0)
ts = WeakLensingTracer(cosmo_emul, z, nz;
Expand Down Expand Up @@ -303,7 +303,7 @@ if test_main
else
ℓs = [10.0, 30.0, 100.0, 300.0, 1000.0]
end
z = Vector(range(0.0, stop=2., length=256))
z = Vector(range(0.01, stop=2., length=1000))
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
tg = NumberCountsTracer(cosmo_EisHu_nonlin, z, nz; b=1.0)
ts = WeakLensingTracer(cosmo_EisHu_nonlin, z, nz;
Expand Down Expand Up @@ -345,7 +345,7 @@ if test_main
else
ℓs = [10.0, 30.0, 100.0, 300.0, 1000.0]
end
z = Vector(range(0., stop=2., length=1000))
z = Vector(range(0.01, stop=2., length=1000))
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
tg = NumberCountsTracer(cosmo_emul_nonlin, z, nz; b=1.0)
ts = WeakLensingTracer(cosmo_emul_nonlin, z, nz;
Expand Down Expand Up @@ -559,7 +559,7 @@ if test_main
function Cl_gg(p::T)::Array{T,1} where T<:Real
cosmo = LimberJack.Cosmology(Ωm=p, tk_mode=:EisHu, Pk_mode=:Halofit,
nz=700, nz_t=700, nz_pk=700)
z = Vector(range(0., stop=2., length=1000))
z = Vector(range(0.01, stop=2., length=1000))
nz = Vector(@. exp(-0.5*((z-0.5)/0.05)^2))
tg = NumberCountsTracer(cosmo, z, nz; b=1.0)
Cℓ_gg = angularCℓs(cosmo, tg, tg, ℓs)
Expand All @@ -569,7 +569,7 @@ if test_main
function Cl_gs(p::T)::Array{T,1} where T<:Real
cosmo = LimberJack.Cosmology(Ωm=p, tk_mode=:EisHu, Pk_mode=:Halofit,
nz=700, nz_t=700, nz_pk=700)
z = Vector(range(0., stop=2., length=1000))
z = Vector(range(0.01, stop=2., length=1000))
nz = Vector(@. exp(-0.5*((z-0.5)/0.05)^2))
tg = NumberCountsTracer(cosmo, z, nz; b=1.0)
ts = WeakLensingTracer(cosmo, z, nz;
Expand All @@ -582,7 +582,7 @@ if test_main
function Cl_ss(p::T)::Array{T,1} where T<:Real
cosmo = LimberJack.Cosmology(Ωm=p, tk_mode=:EisHu, Pk_mode=:Halofit,
nz=700, nz_t=700, nz_pk=700)
z = Vector(range(0., stop=2., length=1000))
z = Vector(range(0.01, stop=2., length=1000))
nz = Vector(@. exp(-0.5*((z-0.5)/0.05)^2))
ts = WeakLensingTracer(cosmo, z, nz;
m=0.0,
Expand All @@ -594,7 +594,7 @@ if test_main
function Cl_sk(p::T)::Array{T,1} where T<:Real
cosmo = LimberJack.Cosmology(Ωm=p, tk_mode=:EisHu, Pk_mode=:Halofit,
nz=500, nz_t=500, nz_pk=700)
z = range(0., stop=2., length=256)
z = range(0.01, stop=2., length=2000)
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
ts = WeakLensingTracer(cosmo, z, nz;
m=0.0,
Expand All @@ -607,7 +607,7 @@ if test_main
function Cl_gk(p::T)::Array{T,1} where T<:Real
cosmo = LimberJack.Cosmology(Ωm=p, tk_mode=:EisHu, Pk_mode=:Halofit,
nz=700, nz_t=700, nz_pk=700)
z = range(0., stop=2., length=256)
z = range(0., stop=2., length=1000)
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
tg = NumberCountsTracer(cosmo, z, nz; b=1.0)
tk = CMBLensingTracer(cosmo)
Expand All @@ -618,7 +618,7 @@ if test_main
function Cl_kk(p::T)::Array{T,1} where T<:Real
cosmo = LimberJack.Cosmology(Ωm=p, tk_mode=:EisHu, Pk_mode=:Halofit,
z_max=1100.0, nz=700, nz_t=700, nz_pk=700)
z = range(0., stop=2., length=256)
z = range(0., stop=2., length=1000)
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
tk = CMBLensingTracer(cosmo)
Cℓ_kk = angularCℓs(cosmo, tk, tk, ℓs)
Expand Down Expand Up @@ -694,7 +694,7 @@ if test_main
function bias(p::T)::Array{T,1} where T<:Real
cosmo = Cosmology(tk_mode=:EisHu, Pk_mode=:Halofit)
cosmo.settings.cosmo_type = typeof(p)
z = Vector(range(0., stop=2., length=1000))
z = Vector(range(0.01, stop=2., length=2000))
nz = Vector(@. exp(-0.5*((z-0.5)/0.05)^2))
tg = NumberCountsTracer(cosmo, z, nz; b=p)
ℓs = [10.0, 30.0, 100.0, 300.0]
Expand All @@ -705,18 +705,18 @@ if test_main
function dz(p::T)::Array{T,1} where T<:Real
cosmo = Cosmology(tk_mode=:EisHu, Pk_mode=:Halofit, nz=300)
cosmo.settings.cosmo_type = typeof(p)
z = Vector(range(0., stop=2., length=1000)) .- p
z = Vector(range(0.01, stop=2., length=2000)) .- p
nz = Vector(@. exp(-0.5*((z-0.5)/0.05)^2))
tg = NumberCountsTracer(cosmo, z, nz; b=1, res=350)
ts = WeakLensingTracer(cosmo, z, nz; m=p, IA_params=[0.0, 0.0])
ℓs = [10.0, 30.0, 100.0, 300.0]
Cℓ_gg = angularCℓs(cosmo, tg, tg, ℓs)
return Cℓ_gg
Cℓ_ss = angularCℓs(cosmo, ts, ts, ℓs)
return Cℓ_ss
end

function mbias(p::T)::Array{T,1} where T<:Real
cosmo = Cosmology(tk_mode=:EisHu, Pk_mode=:Halofit)
cosmo.settings.cosmo_type = typeof(p)
z = range(0., stop=2., length=256)
z = range(0.01, stop=2., length=2000)
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
ts = WeakLensingTracer(cosmo, z, nz; m=p, IA_params=[0.0, 0.0])
ℓs = [10.0, 30.0, 100.0, 300.0]
Expand All @@ -727,7 +727,7 @@ if test_main
function IA_A(p::T)::Array{T,1} where T<:Real
cosmo = Cosmology(tk_mode=:EisHu, Pk_mode=:Halofit, )
cosmo.settings.cosmo_type = typeof(p)
z = range(0., stop=2., length=256)
z = range(0.01, stop=2., length=2000)
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
ts = WeakLensingTracer(cosmo, z, nz; m=2, IA_params=[p, 0.1])
ℓs = [10.0, 30.0, 100.0, 300.0]
Expand All @@ -738,15 +738,15 @@ if test_main
function IA_alpha(p::T)::Array{T,1} where T<:Real
cosmo = Cosmology(tk_mode=:EisHu, Pk_mode=:Halofit)
cosmo.settings.cosmo_type = typeof(p)
z = range(0., stop=2., length=256)
z = range(0.01, stop=2., length=2000)
nz = @. exp(-0.5*((z-0.5)/0.05)^2)
ts = WeakLensingTracer(cosmo, z, nz; m=2, IA_params=[0.3, p])
ℓs = [10.0, 30.0, 100.0, 300.0]
Cℓ_ss = angularCℓs(cosmo, ts, ts, ℓs)
return Cℓ_ss
end

d = 0.00005
d = 0.0005
b_autodiff = ForwardDiff.derivative(bias, 2.0)
b_anal = (bias(2.0+d)-bias(2.0-d))/2d
dz_autodiff = ForwardDiff.derivative(dz, -0.1)
Expand Down
Loading

0 comments on commit f6bcf20

Please sign in to comment.