diff --git a/Project.toml b/Project.toml index c9024f6..16e65e9 100755 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LimberJack" uuid = "6b86205d-155a-4b14-b82d-b6a149ea78f2" authors = ["Jaime Ruiz-Zapatero ", "David Alonso", "Andrina Nicola", "Arrykrishna Mootoovaloo", "James Sullivan", "Marco Bonic"] -version = "0.1.4" +version = "0.1.5" [deps] Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" diff --git a/src/core.jl b/src/core.jl index 2df0f60..2a38c33 100755 --- a/src/core.jl +++ b/src/core.jl @@ -6,28 +6,23 @@ Constructor of settings structure constructor. Kwargs: - `nz::Int=300` : number of nodes in the general redshift array. -- `nz_chi::Int=1000` : number of nodes in the redshift array used to compute matter power spectrum grid. -- `nz_t::Int=350` : number of nodes in the general redshift array. - `nk::Int=500`: number of nodes in the k-scale array used to compute matter power spectrum grid. - `nℓ::Int=300`: number of nodes in the multipoles array. - `using_As::Bool=false`: `True` if using the `As` parameter. - `cosmo_type::Type=Float64` : type of cosmological parameters. - `tk_mode::String=:EisHu` : choice of transfer function. - `Dz_mode::String=:RK2` : choice of method to compute the linear growth factor. -- `Pk_mode::String=:linear` : choice of method to apply non-linear corrections to the matter power spectrum. +- `pk_mode::String=:linear` : choice of method to apply non-linear corrections to the matter power spectrum. Returns: ``` mutable struct Settings nz::Int - nz_chi::Int - nz_t::Int nk::Int nℓ::Int xs zs - zs_t ks ℓs logk @@ -38,21 +33,17 @@ mutable struct Settings cosmo_type::DataType tk_mode::Symbol Dz_mode::Symbol - Pk_mode::Symbol + pk_mode::Symbol end ``` """ mutable struct Settings nz::Int - nz_chi::Int - nz_t::Int nk::Int nℓ::Int xs zs - zs_chi - zs_t ks ℓs logk @@ -63,20 +54,17 @@ mutable struct Settings cosmo_type::DataType tk_mode::Symbol Dz_mode::Symbol - Pk_mode::Symbol + pk_mode::Symbol end Settings(;kwargs...) = begin nz = get(kwargs, :nz, 300) - nz_chi = get(kwargs, :nz_chi, 1000) - nz_t = get(kwargs, :nz_t, 350) nk = get(kwargs, :nk, 500) nℓ = get(kwargs, :nℓ, 300) - xs = LinRange(0, log(1+1100), nz) + z_max = get(kwargs, :z_max, 3.0) + xs = LinRange(0, log(1+z_max), nz) zs = @.(exp(xs) - 1) - zs_chi = 10 .^ Vector(LinRange(-3, log10(1100), nz_chi)) - zs_t = range(0.00001, stop=3.0, length=nz_t) logk = range(log(0.0001), stop=log(100.0), length=nk) ks = exp.(logk) dlogk = log(ks[2]/ks[1]) @@ -87,11 +75,11 @@ Settings(;kwargs...) = begin cosmo_type = get(kwargs, :cosmo_type, Float64) tk_mode = get(kwargs, :tk_mode, :EisHu) Dz_mode = get(kwargs, :Dz_mode, :RK2) - Pk_mode = get(kwargs, :Pk_mode, :linear) - Settings(nz, nz_chi, nz_t, nk, nℓ, - xs, zs, zs_chi, zs_t, ks, ℓs, logk, dlogk, + pk_mode = get(kwargs, :pk_mode, :linear) + Settings(nz, nk, nℓ, + xs, zs, ks, ℓs, logk, dlogk, using_As, - cosmo_type, tk_mode, Dz_mode, Pk_mode) + cosmo_type, tk_mode, Dz_mode, pk_mode) end """ @@ -215,8 +203,8 @@ the primordial power spectrum is calculated using: Depending on the choice of power spectrum mode in the settings, \ the matter power spectrum is either: -- `Pk_mode = :linear` : the linear matter power spectrum. -- `Pk_mode = :halofit` : the Halofit non-linear matter power spectrum (arXiv:astro-ph/0207664). +- `pk_mode = :linear` : the linear matter power spectrum. +- `pk_mode = :halofit` : the Halofit non-linear matter power spectrum (arXiv:astro-ph/0207664). Arguments: - `Settings::MutableStructure` : cosmology constructure settings. @@ -243,41 +231,41 @@ end Cosmology(cpar::CosmoPar, settings::Settings; kwargs...) = begin # Load settings cosmo_type = settings.cosmo_type - zs_chi, nz_chi = settings.zs_chi, settings.nz_chi zs, nz = settings.zs, settings.nz logk, nk = settings.logk, settings.nk ks = settings.ks dlogk = settings.dlogk pk0, pki = lin_Pk0(cpar, settings; kwargs...) # Compute redshift-distance relation - chis = zeros(cosmo_type, nz_chi) - ts = zeros(cosmo_type, nz_chi) - for i in 1:nz_chi - zz = zs_chi[i] + chis = zeros(cosmo_type, nz) + ts = zeros(cosmo_type, nz) + for i in 1:nz + zz = zs[i] chis[i] = quadgk(z -> 1.0/Ez(cpar, z), 0.0, zz, rtol=1E-5)[1] chis[i] *= CLIGHT_HMPC / cpar.h ts[i] = quadgk(z -> 1.0/((1+z)*Ez(cpar, z)), 0.0, zz, rtol=1E-5)[1] ts[i] *= cpar.h*(18.356164383561644*10^9) end # Distance to LSS - chi_LSS = chis[end] + chi_LSS = quadgk(z -> 1.0/Ez(cpar, z), 0.0, 1100., rtol=1E-5)[1] + chi_LSS *= CLIGHT_HMPC / cpar.h # OPT: tolerances, interpolation method - chii = linear_interpolation(zs_chi, Vector(chis), extrapolation_bc=Line()) - ti = linear_interpolation(zs_chi, Vector(ts), extrapolation_bc=Line()) - zi = linear_interpolation(chis, Vector(zs_chi), extrapolation_bc=Line()) + chii = linear_interpolation(zs, Vector(chis), extrapolation_bc=Line()) + ti = linear_interpolation(zs, Vector(ts), extrapolation_bc=Line()) + zi = linear_interpolation(chis, Vector(zs), extrapolation_bc=Line()) Dzs, Dzi, fs8zi = get_growth(cpar, settings; kwargs...) - if settings.Pk_mode == :linear + if settings.pk_mode == :linear Pks = pk0 * (Dzs.^2)' Pki = linear_interpolation((logk, zs), log.(Pks); extrapolation_bc=Line()) - elseif settings.Pk_mode == :Halofit + elseif settings.pk_mode == :Halofit Pki = get_PKnonlin(cpar, zs, ks, pk0, Dzs; cosmo_type=cosmo_type) else @error("Pk mode not implemented") end - Cosmology(settings, cpar, chii, zi, ti, chis[end], + Cosmology(settings, cpar, chii, zi, ti, chi_LSS, chi_LSS, Dzi, fs8zi, pki, Pki) end diff --git a/src/data_utils.jl b/src/data_utils.jl index 9dce28d..8b9b69f 100644 --- a/src/data_utils.jl +++ b/src/data_utils.jl @@ -49,6 +49,33 @@ function _apply_scale_cuts(s, yaml_file) return s end + +""" + make_data(sacc_file, yaml_file) + +Process `sacc` and `yaml` files into a `Meta` structure \ +containing the instructions of how to compose the theory \ +vector and a `npz` file with the neccesary redshift distributions \ +of the tracers involved. + +Arguments: +- `sacc_file` : sacc file +- `yaml_file` : yaml_file + +Returns: +``` +struct Instructions + names : names of tracers. + pairs : pairs of tracers to compute angular spectra for. + types : types of the tracers. + idx : positions of cls in theory vector. + data : data vector. + cov : covariance of the data. + inv_cov : inverse covariance of the data. +end +``` +-files: npz file +""" function make_data(sacc_file, yaml_file; kwargs...) kwargs=Dict(kwargs) @@ -88,7 +115,7 @@ function make_data(sacc_file, yaml_file; kwargs...) # build struct instructions = Instructions(names, pairs, types, idx, - cls, cov, inv_cov) + cls, cov, inv_cov) # Initialize files = Dict{String}{Vector}() @@ -119,4 +146,4 @@ function make_data(sacc_file, yaml_file; kwargs...) end return instructions, files -end \ No newline at end of file +end \ No newline at end of file diff --git a/src/tracers.jl b/src/tracers.jl index f17a41e..fffdbb8 100755 --- a/src/tracers.jl +++ b/src/tracers.jl @@ -27,9 +27,8 @@ Kwargs: Returns: - `NumberCountsTracer::NumberCountsTracer` : Number counts tracer structure. """ -NumberCountsTracer(cosmo::Cosmology, z_n, nz; b=1.0) = begin +NumberCountsTracer(cosmo::Cosmology, z_n, nz; b=1.0, res=1000) = begin nz_int = linear_interpolation(z_n, nz, extrapolation_bc=0) - res = cosmo.settings.nz_t z_w = range(0.00001, stop=z_n[end], length=res) nz_w = nz_int(z_w) nz_norm = integrate(z_w, nz_w, SimpsonEven()) @@ -62,10 +61,10 @@ struct WeakLensingTracer <: Tracer F::Function end -WeakLensingTracer(cosmo::Cosmology, z_n, nz; IA_params = [0.0, 0.0], m=0.0, kwargs...) = begin +WeakLensingTracer(cosmo::Cosmology, z_n, nz; + IA_params = [0.0, 0.0], m=0.0, res=350, kwargs...) = begin nz_int = linear_interpolation(z_n, nz, extrapolation_bc=0) cosmo_type = cosmo.settings.cosmo_type - res =cosmo.settings.nz_t z_w = range(0.00001, stop=z_n[end], length=res) dz_w = (z_w[end]-z_w[1])/res nz_w = nz_int(z_w) @@ -127,9 +126,9 @@ Arguments: Returns: - `CMBLensingTracer::CMBLensingTracer` : CMB lensing tracer structure. """ -CMBLensingTracer(cosmo::Cosmology) = begin +CMBLensingTracer(cosmo::Cosmology; res=350) = begin # chi array - chis = range(0.0, stop=cosmo.chi_max, length=cosmo.settings.nz_t) + chis = range(0.0, stop=cosmo.chi_max, length=res) zs = cosmo.z_of_chi(chis) # Prefactor H0 = cosmo.cpar.h/CLIGHT_HMPC diff --git a/test/runtests.jl b/test/runtests.jl index 7b734f8..e74e4da 100755 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -26,14 +26,15 @@ end test_output = Dict{String}{Vector}() if test_main - cosmo_EisHu = Cosmology(nz=1000, nz_t=1000, nz_pk=1000, nk=1000, tk_mode=:EisHu) + cosmo_EisHu = Cosmology(z_max=1100, nz=1000, nz_t=1000, nz_pk=1000, nk=1000, tk_mode=:EisHu) cosmo_emul = Cosmology(Ωm=(0.12+0.022)/0.75^2, Ωb=0.022/0.75^2, h=0.75, ns=1.0, σ8=0.81, - nz=1000, nz_t=1000, nz_pk=1000, nk=1000, tk_mode=:EmuPk) + z_max=1100, nz=1000, nz_t=1000, nz_pk=1000, nk=1000, tk_mode=:EmuPk) cosmo_emul_As = Cosmology(Ωm=0.27, Ωb=0.046, h=0.7, ns=1.0, As=2.097e-9, - nz=1000, nz_t=1000, nz_pk=1000, nk=1000, tk_mode=:EmuPk) - cosmo_EisHu_nonlin = Cosmology(nk=1000, nz=1000, nz_t=1000, nz_pk=1000, tk_mode=:EisHu, Pk_mode=:Halofit) + z_max=1100, nz=1000, nz_t=1000, nz_pk=1000, nk=1000, tk_mode=:EmuPk) + cosmo_EisHu_nonlin = Cosmology(nk=1000, nz=1000, nz_t=1000, nz_pk=1000, + z_max=1100, tk_mode=:EisHu, pk_mode=:Halofit) cosmo_emul_nonlin = Cosmology(Ωm=(0.12+0.022)/0.75^2, Ωb=0.022/0.75^2, h=0.75, ns=1.0, σ8=0.81, - nk=1000, nz=1000, nz_t=1000, nz_pk=1000, tk_mode=:EmuPk, Pk_mode=:Halofit) + z_max=1100, nk=1000, nz=1000, nz_t=1000, nz_pk=1000, tk_mode=:EmuPk, pk_mode=:Halofit) @testset "Main tests" begin @testset "CreateCosmo" begin @@ -616,7 +617,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, - nz=700, nz_t=700, nz_pk=700) + z_max=1100.0, nz=700, nz_t=700, nz_pk=700) z = range(0., stop=2., length=256) nz = @. exp(-0.5*((z-0.5)/0.05)^2) tk = CMBLensingTracer(cosmo) @@ -706,7 +707,7 @@ if test_main cosmo.settings.cosmo_type = typeof(p) z = Vector(range(0., stop=2., length=1000)) .- p nz = Vector(@. exp(-0.5*((z-0.5)/0.05)^2)) - tg = NumberCountsTracer(cosmo, z, nz; b=1) + tg = NumberCountsTracer(cosmo, z, nz; b=1, res=350) ℓs = [10.0, 30.0, 100.0, 300.0] Cℓ_gg = angularCℓs(cosmo, tg, tg, ℓs) return Cℓ_gg diff --git a/test/test_output.npz b/test/test_output.npz index f81b4e7..552a811 100644 Binary files a/test/test_output.npz and b/test/test_output.npz differ