Merge pull request #10 from JaimeRZP/inference_bug

Inference bug

Project.toml

@@ -1,7 +1,7 @@
 name = "LimberJack"
 uuid = "6b86205d-155a-4b14-b82d-b6a149ea78f2"
 authors = ["Jaime Ruiz-Zapatero <[email protected]>", "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"

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
 

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
+end 

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

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