Unification branch

.github/workflows/test.yml

@@ -18,8 +18,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9'
           - '1.10'
+          - '1.11'
         os:
           - ubuntu-latest
         arch:

Project.toml

@@ -5,22 +5,22 @@ version = "0.2.0"
 
 [deps]
 AbstractCosmologicalEmulators = "c83c1981-e5c4-4837-9eb8-c9b1572acfc6"
-Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 Integrals = "de52edbc-65ea-441a-8357-d3a637375a31"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 Memoization = "6fafb56a-5788-4b4e-91ca-c0cea6611c73"
+NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
-AbstractCosmologicalEmulators = "0.3.3, 0.4, 0.5, 0.6"
-Adapt = "3"
+AbstractCosmologicalEmulators = "0.6"

README.md

@@ -6,6 +6,11 @@
 
 Repository containing the EFfective Field theORy surrogaTe.
 
+The code has been used in the following publications:
+
+- H. Zhang, M. Bonici, G. D'Amico, S. Paradiso, W. J. Percival, [HOD-informed prior for EFT-based full-shape analyses of LSS](https://arxiv.org/abs/2409.12937)
+- S. Paradiso, **M. Bonici**, M. Chen, W. J. Percival, G. D'Amico, H. Zhang, G. McGee, [Reducing nuisance prior sensitivity via non-linear reparameterization, with application to EFT analyses of large-scale structure](https://arxiv.org/abs/2412.03503)
+
 ### Authors
 
 - Marco Bonici, PostDoctoral Researcher at Waterloo Centre for Astrophysics

src/Effort.jl

@@ -16,13 +16,15 @@ using Integrals
 using LinearAlgebra
 using SparseArrays
 using Tullio
+using NPZ
 using Zygote
+import JSON.parsefile
 
 const c_0 = 2.99792458e5
 
 function __init__()
-    min_y = get_y(0,0) #obvious, I knowadd OrdinaryDiffEqTsit5
-    max_y = get_y(1,10)
+    min_y = _get_y(0,0) #obvious, I knowadd OrdinaryDiffEqTsit5
+    max_y = _get_y(1,10)
     y_grid = vcat(LinRange(min_y, 100, 100), LinRange(100.1, max_y, 1000))
     F_grid = [_F(y) for y in y_grid]
     global F_interpolant = AkimaInterpolation(F_grid, y_grid)

src/background.jl

@@ -1,6 +1,6 @@
 #TODO: this part should be moved to a dedicate package. While necessary to a full Effort
 #functionality, this could be factorized to a new module, specifically taylored to this goal
-# and maybe used in other packages
+# and maybe used in other packages, maybe in AbstractCosmologicalEmulators?
 
 function _F(y)
     f(x, y) = x^2 * √(x^2 + y^2) / (1 + exp(x))
@@ -10,7 +10,7 @@
     return sol
 end
 
-function get_y(mν, a; kB=8.617342e-5, Tν=0.71611 * 2.7255)
+function _get_y(mν, a; kB=8.617342e-5, Tν=0.71611 * 2.7255)
     return mν * a / (kB * Tν)
 end
 
@@ -23,109 +23,112 @@
 end
 
 function _ΩνE2(a, Ωγ0, mν; kB=8.617342e-5, Tν=0.71611 * 2.7255, Neff=3.044)
-    Γν = (4 / 11)^(1 / 3) * (Neff / 3)^(1 / 4)#0.71649^4
-    return 15 / π^4 * Γν^4 * Ωγ0 / a^4 * F_interpolant(get_y(mν, a))
+    Γν = (4 / 11)^(1 / 3) * (Neff / 3)^(1 / 4)
+    return 15 / π^4 * Γν^4 * Ωγ0 / a^4 * F_interpolant(_get_y(mν, a))
 end
 
 function _ΩνE2(a, Ωγ0, mν::Array; kB=8.617342e-5, Tν=0.71611 * 2.7255, Neff=3.044)
     Γν = (4 / 11)^(1 / 3) * (Neff / 3)^(1 / 4)#0.71649^4
     sum_interpolant = 0.0
     for mymν in mν
-        sum_interpolant += F_interpolant(get_y(mymν, a))
+        sum_interpolant += F_interpolant(_get_y(mymν, a))
     end
     return 15 / π^4 * Γν^4 * Ωγ0 / a^4 * sum_interpolant
 end
 
 function _dΩνE2da(a, Ωγ0, mν; kB=8.617342e-5, Tν=0.71611 * 2.7255, Neff=3.044)
     Γν = (4 / 11)^(1 / 3) * (Neff / 3)^(1 / 4)#0.71649^4
-    return 15 / π^4 * Γν^4 * Ωγ0 * (-4 * F_interpolant(get_y(mν, a)) / a^5 +
-                                    dFdy_interpolant(get_y(mν, a)) / a^4 * (mν / kB / Tν))
+    return 15 / π^4 * Γν^4 * Ωγ0 * (-4 * F_interpolant(_get_y(mν, a)) / a^5 +
+                                    dFdy_interpolant(_get_y(mν, a)) / a^4 * (mν / kB / Tν))
 end
 
 function _dΩνE2da(a, Ωγ0, mν::Array; kB=8.617342e-5, Tν=0.71611 * 2.7255, Neff=3.044)
     Γν = (4 / 11)^(1 / 3) * (Neff / 3)^(1 / 4)#0.71649^4
     sum_interpolant = 0.0
     for mymν in mν
-        sum_interpolant += -4 * F_interpolant(get_y(mymν, a)) / a^5 +
-                           dFdy_interpolant(get_y(mymν, a)) / a^4 * (mymν / kB / Tν)
+        sum_interpolant += -4 * F_interpolant(_get_y(mymν, a)) / a^5 +
+                           dFdy_interpolant(_get_y(mymν, a)) / a^4 * (mymν / kB / Tν)
     end
     return 15 / π^4 * Γν^4 * Ωγ0 * sum_interpolant
 end
 
-function _E_a(a, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
+function _E_a(a, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
     Ωγ0 = 2.469e-5 / h^2
     Ων0 = _ΩνE2(1.0, Ωγ0, mν)
-    ΩΛ0 = 1.0 - (Ωγ0 + Ωm0 + Ων0)
-    return sqrt(Ωγ0 * a^-4 + Ωm0 * a^-3 + ΩΛ0 * _ρDE_a(a, w0, wa) + _ΩνE2(a, Ωγ0, mν))
+    ΩΛ0 = 1.0 - (Ωγ0 + Ωcb0 + Ων0)
+    return @. sqrt(Ωγ0 * a^-4 + Ωcb0 * a^-3 + ΩΛ0 * _ρDE_a(a, w0, wa) + _ΩνE2(a, Ωγ0, mν))
 end
 
-function _E_z(z, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
+function _E_z(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
     a = _a_z.(z)
-    return _E_a(a, Ωm0, h; mν=mν, w0=w0, wa=wa)
+    return _E_a(a, Ωcb0, h; mν=mν, w0=w0, wa=wa)
 end
 
-_H_a(a, Ωγ0, Ωm0, mν, h, w0, wa) = 100 * h * _E_a(a, Ωm0, h; mν=mν, w0=w0, wa=wa)
+#TODO check whether to cancel this one
+_H_a(a, Ωcb0, mν, h, w0, wa) = 100 * h * _E_a(a, Ωcb0, h; mν=mν, w0=w0, wa=wa)
 
-function _χ_z(z, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
-    p = [Ωm0, h, mν, w0, wa]
+#TODO check whether to cancel this one
+function _χ_z(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
+    p = [Ωcb0, h, mν, w0, wa]
     f(x, p) = 1 / _E_a(_a_z(x), p[1], p[2]; mν=p[3], w0=p[4], wa=p[5])
     domain = (zero(eltype(z)), z) # (lb, ub)
     prob = IntegralProblem(f, domain, y; preltol=1e-12)
     sol = solve(prob, QuadGKJL())[1]
     return sol * c_0 / (100 * h)
 end
 
-function _dEda(a, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
+#TODO check whether to cancel this one
+function _dEda(a, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
     Ωγ0 = 2.469e-5 / h^2
     Ων0 = _ΩνE2(1.0, Ωγ0, mν)
-    ΩΛ0 = 1.0 - (Ωγ0 + Ωm0 + Ων0)
-    return 0.5 / _E_a(a, Ωm0, h; mν=mν, w0=w0, wa=wa) * (-3(Ωm0)a^-4 - 4Ωγ0 * a^-5 +
-                                                         ΩΛ0 * _dρDEda(a, w0, wa) + _dΩνE2da(a, Ωγ0, mν))
+    ΩΛ0 = 1.0 - (Ωγ0 + Ωcb0 + Ων0)
+    return 0.5 / _E_a(a, Ωcb0, h; mν=mν, w0=w0, wa=wa) *
+    (-3(Ωcb0)a^-4 - 4Ωγ0 * a^-5 + ΩΛ0 * _dρDEda(a, w0, wa) + _dΩνE2da(a, Ωγ0, mν))
 end
 
-function _dlogEdloga(a, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
+function _dlogEdloga(a, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
     Ωγ0 = 2.469e-5 / h^2
     Ων0 = _ΩνE2(1.0, Ωγ0, mν)
-    ΩΛ0 = 1.0 - (Ωγ0 + Ωm0 + Ων0)
-    return a * 0.5 / (_E_a(a, Ωm0, h; mν=mν, w0=w0, wa=wa)^2) * (-3(Ωm0)a^-4 - 4Ωγ0 * a^-
+    ΩΛ0 = 1.0 - (Ωγ0 + Ωcb0 + Ων0)
+    return a * 0.5 / (_E_a(a, Ωcb0, h; mν=mν, w0=w0, wa=wa)^2) * (-3(Ωcb0)a^-4 - 4Ωγ0 * a^-
                                                                  5 + ΩΛ0 * _dρDEda(a, w0, wa) + _dΩνE2da(a, Ωγ0, mν))
 end
 
-function _ΩMa(a, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
+function _Ωcba(a, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
     Ωγ0 = 2.469e-5 / h^2
     Ων0 = _ΩνE2(1.0, Ωγ0, mν)
-    return (Ωm0 + Ων0) * a^-3 / (_E_a(a, Ωm0, h; mν=mν, w0=w0, wa=wa))^2
+    return (Ωcb0 + Ων0) * a^-3 / (_E_a(a, Ωcb0, h; mν=mν, w0=w0, wa=wa))^2
 end
 
-function _r̃_z_check(z, ΩM, h; mν=0.0, w0=-1.0, wa=0.0)
-    p = [ΩM, h, mν, w0, wa]
+function _r̃_z_check(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
+    p = [Ωcb0, h, mν, w0, wa]
     f(x, p) = 1 / _E_a(_a_z(x), p[1], p[2]; mν=p[3], w0=p[4], wa=p[5])
     domain = (zero(eltype(z)), z) # (lb, ub)
     prob = IntegralProblem(f, domain, p; reltol=1e-12)
     sol = solve(prob, QuadGKJL())[1]
     return sol
 end
 
-function _r̃_z(z, ΩM, h; mν=0.0, w0=-1.0, wa=0.0)
+function _r̃_z(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
     z_array, weigths_array = _transformed_weights(FastGaussQuadrature.gausslegendre, 9, 0, z)
-    integrand_array = @. 1.0 / _E_a(_a_z(z_array), ΩM, h; mν=mν, w0=w0, wa=wa)
+    integrand_array = 1.0 ./ _E_a(_a_z(z_array), Ωcb0, h; mν=mν, w0=w0, wa=wa)
     return dot(weigths_array, integrand_array)
 end
 
-function _r_z_check(z, ΩM, h; mν=0.0, w0=-1.0, wa=0.0)
-    return c_0 * _r̃_z_check(z, ΩM, h; mν=mν, w0=w0, wa=wa) / (100 * h)
+function _r_z_check(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
+    return c_0 * _r̃_z_check(z, Ωcb0, h; mν=mν, w0=w0, wa=wa) / (100 * h)
 end
 
-function _r_z(z, ΩM, h; mν=0.0, w0=-1.0, wa=0.0)
-    return c_0 * _r̃_z(z, ΩM, h; mν=mν, w0=w0, wa=wa) / (100 * h)
+function _r_z(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
+    return c_0 * _r̃_z(z, Ωcb0, h; mν=mν, w0=w0, wa=wa) / (100 * h)
 end
 
-function _d̃A_z(z, ΩM, h; mν=0.0, w0=-1.0, wa=0.0)
-    return _r̃_z(z, ΩM, h; mν=mν, w0=w0, wa=wa) / (1 + z)
+function _d̃A_z(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
+    return _r̃_z(z, Ωcb0, h; mν=mν, w0=w0, wa=wa) / (1 + z)
 end
 
-function _dA_z(z, ΩM, h; mν=0.0, w0=-1.0, wa=0.0)
-    return _r_z(z, ΩM, h; mν=mν, w0=w0, wa=wa) / (1 + z)
+function _dA_z(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
+    return _r_z(z, Ωcb0, h; mν=mν, w0=w0, wa=wa) / (1 + z)
 end
 
 function _ρDE_z(z, w0, wa)
@@ -140,16 +143,17 @@
     return 3 * (-(1 + w0 + wa) / a + wa) * _ρDE_a(a, w0, wa)
 end
 
-function _X_z(z, ΩM, w0, wa)
-    return ΩM * ((1 + z)^3) / ((1 - ΩM) * _ρDE_z(z, w0, wa))
+#TODO check whether we can remove this function
+function _X_z(z, Ωcb0, w0, wa)
+    return Ωcb0 * ((1 + z)^3) / ((1 - Ωcb0) * _ρDE_z(z, w0, wa))
 end
 
 function _w_z(z, w0, wa)
     return w0 + wa * z / (1 + z)
 end
 
 function _growth!(du, u, p, loga)
-    Ωm0 = p[1]
+    Ωcb0 = p[1]
     mν = p[2]
     h = p[3]
     w0 = p[4]
@@ -158,38 +162,38 @@
     D = u[1]
     dD = u[2]
     du[1] = dD
-    du[2] = -(2 + _dlogEdloga(a, Ωm0, h; mν=mν, w0=w0, wa=wa)) * dD +
-            1.5 * _ΩMa(a, Ωm0, h; mν=mν, w0=w0, wa=wa) * D
+    du[2] = -(2 + _dlogEdloga(a, Ωcb0, h; mν=mν, w0=w0, wa=wa)) * dD +
+            1.5 * _Ωcba(a, Ωcb0, h; mν=mν, w0=w0, wa=wa) * D
 end
 
 function _a_z(z)
-    return 1 / (1 + z)
+    return @. 1 / (1 + z)
 end
 
-function growth_solver(Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
+function _growth_solver(Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
     amin = 1 / 139
     u₀ = [amin, amin]
 
-    logaspan = (log(amin), log(1.01))
-    Ωγ0 = 2.469e-5 / h^2
+    logaspan = (log(amin), log(1.01))#to ensure we cover the relevant range
+    #Ωγ0 = 2.469e-5 / h^2
 
-    p = (Ωm0, mν, h, w0, wa)
+    p = (Ωcb0, mν, h, w0, wa)
 
     prob = ODEProblem(_growth!, u₀, logaspan, p)
 
     sol = solve(prob, OrdinaryDiffEq.Tsit5(), abstol=1e-6, reltol=1e-6; verbose=false)
     return sol
 end
 
-function growth_solver(z, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
+function _growth_solver(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
     amin = 1 / 139
     loga = vcat(log.(_a_z.(z)), 0.0)
     u₀ = [amin, amin]
 
     logaspan = (log(amin), log(1.01))
-    Ωγ0 = 2.469e-5 / h^2
+    #Ωγ0 = 2.469e-5 / h^2
 
-    p = [Ωm0, mν, h, w0, wa]
+    p = [Ωcb0, mν, h, w0, wa]
 
     prob = ODEProblem(_growth!, u₀, logaspan, p)
 
@@ -205,19 +209,19 @@
     return (sol(log(_a_z(z)))[1, :]/sol(log(_a_z(0.0)))[1, :])[1, 1]
 end
 
-function _D_z(z, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
-    sol = growth_solver(z, Ωm0, h; mν=mν, w0=w0, wa=wa)
+function _D_z(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
+    sol = _growth_solver(z, Ωcb0, h; mν=mν, w0=w0, wa=wa)
     D_z = reverse(sol[1, 1:end-1]) ./ sol[1, end]
     return D_z
 end
 
-function _D_z_old(z, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
-    sol = growth_solver(Ωm0, h; mν=mν, w0=w0, wa=wa)
+function _D_z_old(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
+    sol = _growth_solver(Ωcb0, h; mν=mν, w0=w0, wa=wa)
     return _D_z_old(z, sol)
 end
 
-function _D_z_unnorm(z, Ωm0, h; mν=0.0, w0=-1.0, wa=0.0)
-    sol = growth_solver(Ωm0, h; mν=mν, w0=w0, wa=wa)
+function _D_z_unnorm(z, Ωcb0, h; mν=0.0, w0=-1.0, wa=0.0)
+    sol = _growth_solver(Ωcb0, h; mν=mν, w0=w0, wa=wa)
     return _D_z_unnorm(z, sol)
 end
 
@@ -239,14 +243,14 @@
     return _f_a_old(a, sol)
 end
 
-function _f_z_old(z, Ωm0, h; mν=0, w0=-1.0, wa=0.0)
+function _f_z_old(z, Ωcb0, h; mν=0, w0=-1.0, wa=0.0)
     a = _a_z.(z)
-    sol = growth_solver(Ωm0, h; mν=mν, w0=w0, wa=wa)
+    sol = _growth_solver(Ωcb0, h; mν=mν, w0=w0, wa=wa)
     return _f_a_old(a, sol)
 end
 
-function _f_z(z, Ωm0, h; mν=0, w0=-1.0, wa=0.0)
-    sol = growth_solver(z, Ωm0, h; mν=mν, w0=w0, wa=wa)
+function _f_z(z, Ωcb0, h; mν=0, w0=-1.0, wa=0.0)
+    sol = _growth_solver(z, Ωcb0, h; mν=mν, w0=w0, wa=wa)
     D = sol[1, 1:end-1]
     D_prime = sol[2, 1:end-1]
     return @. 1 / D * D_prime