diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 6eee27c..61f428e 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -15,14 +15,14 @@ jobs:
build:
runs-on: ubuntu-latest
steps:
- - uses: actions/checkout@v2
+ - uses: actions/checkout@v4
- uses: julia-actions/setup-julia@latest
with:
- version: '1.8'
+ version: '1.10'
- name: Install dependencies
run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
- name: Build and deploy
env:
GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
- run: julia --project=docs/ docs/make.jl
\ No newline at end of file
+ run: julia --project=docs/ docs/make.jl
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
new file mode 100644
index 0000000..28528da
--- /dev/null
+++ b/.github/workflows/test.yml
@@ -0,0 +1,52 @@
+name: CI
+on:
+ schedule:
+ - cron: 0 0 * * *
+ pull_request:
+ branches:
+ - main
+ - develop
+ - abstractcosmoemu
+ push:
+ branches:
+ - '**'
+ tags: '*'
+jobs:
+ test:
+ name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+ runs-on: ${{ matrix.os }}
+ strategy:
+ fail-fast: false
+ matrix:
+ version:
+ - '1.7'
+ - '1.8'
+ - '1.9'
+ - '1.10'
+ os:
+ - ubuntu-latest
+ arch:
+ - x64
+ steps:
+ - uses: actions/checkout@v4
+ - uses: julia-actions/setup-julia@latest
+ with:
+ version: ${{ matrix.version }}
+ arch: ${{ matrix.arch }}
+ - uses: actions/cache@v4
+ env:
+ cache-name: cache-artifacts
+ with:
+ path: ~/.julia/artifacts
+ key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+ restore-keys: |
+ ${{ runner.os }}-test-${{ env.cache-name }}-
+ ${{ runner.os }}-test-
+ ${{ runner.os }}-
+ - uses: julia-actions/julia-runtest@v1
+ - uses: julia-actions/julia-processcoverage@v1
+ - uses: codecov/codecov-action@v4
+ env:
+ CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
+ with:
+ file: lcov.info
diff --git a/Project.toml b/Project.toml
index 7fb6e9f..64e9b05 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,9 +1,11 @@
name = "Effort"
uuid = "ace19e24-02c2-4bce-bc41-b57e89c0fc53"
authors = ["marcobonici <bonici.marco@gmail.com>"]
-version = "0.1.1"
+version = "0.2.0"
[deps]
+AbstractCosmologicalEmulators = "c83c1981-e5c4-4837-9eb8-c9b1572acfc6"
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
@@ -11,4 +13,7 @@ LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
Memoization = "6fafb56a-5788-4b4e-91ca-c0cea6611c73"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
-SimpleChains = "de6bee2f-e2f4-4ec7-b6ed-219cc6f6e9e5"
+
+[compat]
+AbstractCosmologicalEmulators = "0.3.3"
+Adapt = "3"
diff --git a/README.md b/README.md
index ccfeca0..68a2e65 100644
--- a/README.md
+++ b/README.md
@@ -1,13 +1,12 @@
# Effort.jl
-[](https://cosmologicalemulators.github.io/Effort.jl/stable)
-[](https://cosmologicalemulators.github.io/Effort.jl/dev)
-
-[](https://github.com/invenia/BlueStyle)
-[](https://github.com/SciML/ColPrac)
+
+| **Documentation** | **Build Status** | **Code style** |
+|:--------:|:----------------:|:----------------:|
+| [](https://cosmologicalemulators.github.io/Effort.jl/dev) [](https://cosmologicalemulators.github.io/Effort.jl/stable) | [](https://github.com/CosmologicalEmulators/Effort.jl/actions) [](https://codecov.io/gh/CosmologicalEmulators/Effort.jl) | [](https://github.com/invenia/BlueStyle) [](https://github.com/SciML/ColPrac) |
Repository containing the EFfective Field theORy surrogaTe.
### Authors
-- Marco Bonici, INAF - Institute of Space Astrophysics and Cosmic Physics (IASF), Milano
-- Guido D'Amico, Università Degli Studi di Parma
+- Marco Bonici, PostDoctoral Researcher at Waterloo Centre for Astrophysics
+- Guido D'Amico, Associate Professor at Università Degli Studi di Parma
diff --git a/src/Effort.jl b/src/Effort.jl
index 8febeb5..918821e 100644
--- a/src/Effort.jl
+++ b/src/Effort.jl
@@ -1,6 +1,8 @@
module Effort
using Base: @kwdef
+using AbstractCosmologicalEmulators
+import AbstractCosmologicalEmulators.get_emulator_description
using DataInterpolations
using FastGaussQuadrature
using LegendrePolynomials
@@ -8,7 +10,19 @@ using LoopVectorization
using Memoization
using OrdinaryDiffEq
using QuadGK
-using SimpleChains
+
+const c_0 = 2.99792458e5
+
+function __init__()
+ min_y = get_y(0,0) #obvious, I know
+ 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)
+ y_grid = vcat(LinRange(min_y, 10., 10000), LinRange(10.1, max_y, 10000))
+ dFdy_grid = [_dFdy(y) for y in y_grid]
+ global dFdy_interpolant = AkimaInterpolation(dFdy_grid, y_grid)
+end
include("background.jl")
include("neural_networks.jl")
diff --git a/src/background.jl b/src/background.jl
index c15338e..b3f9d86 100644
--- a/src/background.jl
+++ b/src/background.jl
@@ -1,28 +1,122 @@
-#TODO: this part shoul be moved to a dedicate package. While necessary to a full Effort
+#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
-function _E_z(z, ΩM, w0, wa)
- return sqrt(ΩM*(1+z)^3+(1-ΩM)*_ρDE_z(z, w0, wa))
+function _F(y)
+ result, _ = quadgk(x -> x^2*√(x^2+y^2)/(1+exp(x)), 0, Inf, rtol=1e-12)
+ return result
end
-function _H_z(z, H0, ΩM, w0, wa)
- return H0*_E_z(z, ΩM, w0, wa)
+function get_y(mν, a; kB = 8.617342e-5, Tν = 0.71611*2.7255)
+ return mν * a / (kB * Tν)
end
-function _r̃_z(z, ΩM, w0, wa)
- integral, _ = quadgk(x -> 1 / _E_z(x, ΩM, w0, wa), 0, z, rtol=1e-10)
+function _dFdy(y)
+ result, _ = quadgk(x -> x^2/((1+exp(x))*√(x^2+y^2)), 0, Inf, rtol=1e-12)
+ return result*y
+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))
+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.
+ for mymν in mν
+ 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ν))
+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.
+ 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ν)
+ end
+ return 15/π^4*Γν^4*Ωγ0*sum_interpolant
+end
+
+function _E_a(a, Ωm0, h; mν =0., w0=-1., wa=0.)
+ Ωγ0 = 2.469e-5/h^2
+ Ων0 = _ΩνE2(1., Ωγ0, mν)
+ ΩΛ0 = 1. - (Ωγ0 + Ωm0 + Ων0)
+ return sqrt(Ωγ0*a^-4 + Ωm0*a^-3 + ΩΛ0 * _ρDE_a(a, w0, wa)+ _ΩνE2(a, Ωγ0, mν))
+end
+
+function _E_z(z, Ωm0, h; mν =0., w0=-1., wa=0.)
+ a = _a_z.(z)
+ return _E_a(a, Ωm0, 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)
+
+function _χ_z(z, Ωm0, h; mν =0., w0=-1., wa=0.)
+ integral, _ = quadgk(x -> 1 /
+ _E_a(_a_z(x), Ωm0, h; mν =mν, w0=w0, wa=wa), 0, z, rtol=1e-6)
+ return integral*c_0/(100*h)
+end
+
+function _dEda(a, Ωm0, h; mν =0., w0=-1., wa=0.)
+ Ωγ0 = 2.469e-5/h^2
+ Ων0 = _ΩνE2(1., Ωγ0, mν)
+ ΩΛ0 = 1. - (Ωγ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ν))
+end
+
+function _dlogEdloga(a, Ωm0, h; mν =0., w0=-1., wa=0.)
+ Ωγ0 = 2.469e-5/h^2
+ Ων0 = _ΩνE2(1., Ωγ0, mν)
+ ΩΛ0 = 1. - (Ωγ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^-
+ 5+ΩΛ0*_dρDEda(a, w0, wa)+_dΩνE2da(a, Ωγ0, mν))
+end
+
+function _ΩMa(a, Ωm0, h; mν =0., w0=-1., wa=0.)
+ Ωγ0 = 2.469e-5/h^2
+ Ων0 = _ΩνE2(1., Ωγ0, mν)
+ return (Ωm0 + Ων0 )*a^-3 / (_E_a(a, Ωm0, h; mν =mν, w0=w0, wa=wa))^2
+end
+
+function _r̃_z(z, ΩM, h; mν =0., w0=-1., wa=0.)
+ integral, _ = quadgk(x -> 1 / _E_z(x, ΩM, h; mν =mν, w0=w0, wa=wa), 0, z, rtol=1e-10)
return integral
end
-function _d̃A_z(z, ΩM, w0, wa)
- return (1+z) * _r̃_z(z, ΩM, w0, wa)
+function _r_z(z, ΩM, h; mν =0., w0=-1., wa=0.)
+ return c_0 * _r̃_z(z, ΩM, h; mν =mν, w0=w0, wa=wa) / (100*h)
+end
+
+function _d̃A_z(z, ΩM, h; mν =0., w0=-1., wa=0.)
+ return _r̃_z(z, ΩM, h; mν =mν, w0=w0, wa=wa) / (1+z)
+end
+
+function _dA_z(z, ΩM, h; mν =0., w0=-1., wa=0.)
+ return _r_z(z, ΩM, h; mν =mν, w0=w0, wa=wa) / (1+z)
end
function _ρDE_z(z, w0, wa)
return (1+z)^(3.0 * (1.0 + w0 + wa)) * exp(-3.0 * wa * z /(1+z))
end
+function _ρDE_a(a, w0, wa)
+ return a^(-3.0 * (1.0 + w0 + wa)) * exp(3.0 * wa * (a-1))
+end
+
+function _dρDEda(a, w0, wa)
+ 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))
end
@@ -31,68 +125,107 @@ function _w_z(z, w0, wa)
return w0+wa*z/(1+z)
end
-function _growth!(du,u,p,a)
- ΩM = p[1]
- w0 = p[2]
- wa = p[3]
- z = 1.0 / a - 1.0
- G = u[1]
- dG = u[2]
- du[1] = dG
- du[2] = -(3.5-1.5*_w_z(z, w0, wa)/(1+_X_z(z, ΩM, w0, wa)))*dG/a-1.5*(1-_w_z(z, w0,wa))/(1+_X_z(z, ΩM, w0, wa))*G/(a^2)
+function _growth!(du,u,p,loga)
+ #Ωγ0 = p[1]
+ Ωm0 = p[1]
+ mν = p[2]
+ h = p[3]
+ w0 = p[4]
+ wa = p[5]
+ a = exp(loga)
+ 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
end
function _a_z(z)
return 1/(1+z)
end
-function growth_solver(ΩM, w0, wa)
- u₀ = [1.0,0.0]
+function growth_solver(Ωm0, h; mν =0., w0=-1., wa=0.)
+ amin = 1/139
+ u₀ = [amin, amin]
+
+ logaspan = (log(amin), log(1.01))
+ Ωγ0 = 2.469e-5/h^2
+
+ p = (Ωm0, 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., w0=-1., wa=0.)
+ amin = 1/139
+ loga = vcat(log.(_a_z.(z)), 0.)
+ u₀ = [amin, amin]
- aspan = (0.99e-3, 1.01)
+ logaspan = (log(amin), log(1.01))
+ Ωγ0 = 2.469e-5/h^2
- p = [ΩM, w0, wa]
+ p = [Ωm0, mν, h, w0, wa]
- prob = ODEProblem(_growth!, u₀, aspan, p)
+ prob = ODEProblem(_growth!, u₀, logaspan, p)
- sol = solve(prob, Tsit5(), abstol=1e-6, reltol=1e-6;verbose=false)
+ sol = solve(prob, OrdinaryDiffEq.Tsit5(), abstol=1e-6, reltol=1e-6; saveat=loga)[1:2,:]
return sol
end
-function _D_z(z::Array, sol::SciMLBase.ODESolution)
- [u for (u,t) in sol.(_a_z.(z))] .* _a_z.(z) ./ (sol(_a_z(0.))[1,:])
+function _D_z_old(z::Array, sol::SciMLBase.ODESolution)
+ [u for (u,t) in sol.(log.(_a_z.(z)))] ./ (sol(log(_a_z(0.)))[1,:])
+end
+
+function _D_z_old(z, sol::SciMLBase.ODESolution)
+ return (sol(log(_a_z(z)))[1,:]/sol(log(_a_z(0.)))[1,:])[1,1]
+end
+
+function _D_z(z, Ωm0, h; mν =0., w0=-1., wa=0.)
+ sol = growth_solver(z,Ωm0, 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., w0=-1., wa=0.)
+ sol = growth_solver(Ωm0, h; mν =mν, w0=w0, wa=wa)
+ return _D_z_old(z, sol)
end
-function _D_z(z, sol::SciMLBase.ODESolution)
- return (Effort._a_z(z) .* sol(Effort._a_z(z))[1,:]/sol(Effort._a_z(0.))[1,:])[1,1]
+function _D_z_unnorm(z, Ωm0, h; mν =0., w0=-1., wa=0.)
+ sol = growth_solver(Ωm0, h; mν =mν, w0=w0, wa=wa)
+ return _D_z_unnorm(z, sol)
end
-function _D_z(z, ΩM, w0, wa)
- sol = growth_solver(ΩM, w0, wa)
- return _D_z(z, sol)
+function _D_z_unnorm(z::Array, sol::SciMLBase.ODESolution)
+ [u for (u,t) in sol.(log.(_a_z.(z)))]
end
-function _f_a(a, sol::SciMLBase.ODESolution)
- G, G_prime = sol(a)
- D = G * a
- D_prime = G_prime * a + G
- return a / D * D_prime
+function _D_z_unnorm(z, sol::SciMLBase.ODESolution)
+ return (sol(log(_a_z(z)))[1,:])[1,1]
end
-function _f_a(a::Array, sol::SciMLBase.ODESolution)
- G = [u for (u,t) in sol.(a)]
- G_prime = [t for (u,t) in sol.(a)]
- D = G .* a
- D_prime = G_prime .* a .+ G
- return a ./ D .* D_prime
+function _f_a_old(a, sol::SciMLBase.ODESolution)
+ D, D_prime = sol.(log.(a))
+ return @. 1 / D * D_prime
+end
+
+function _f_z_old(z, sol::SciMLBase.ODESolution)
+ a = _a_z.(z)
+ return _f_a_old(a, sol)
end
-function _f_z(z, sol::SciMLBase.ODESolution)
+function _f_z_old(z, Ωm0, h; mν =0 , w0=-1., wa=0.)
a = _a_z.(z)
- return _f_a(a, sol)
+ sol = growth_solver(Ωm0, h; mν =mν, w0=w0, wa=wa)
+ return _f_a_old(a, sol)
end
-function _f_z(z, ΩM, w0, wa)
- sol = growth_solver(ΩM, w0, wa)
- return _f_z(z, sol)
+function _f_z(z, Ωm0, h; mν =0 , w0=-1., wa=0.)
+ sol = growth_solver(z, Ωm0, 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
end
diff --git a/src/eft_commands.jl b/src/eft_commands.jl
index bf07dfe..c87da31 100644
--- a/src/eft_commands.jl
+++ b/src/eft_commands.jl
@@ -49,14 +49,14 @@ function bias_multiplication!(input_array, bias_array, Pk_input)
end
end
-function get_stoch(cϵ0, cϵ1, cϵ2, n_bar, k_grid::Array, k_nl=0.7)
+function get_stoch(cϵ0, cϵ1, cϵ2, n_bar, k_grid::Array; k_nl=0.7)
P_stoch = zeros(3, length(k_grid))
P_stoch[1,:] = @. 1/n_bar*(cϵ0 + cϵ1*(k_grid/k_nl)^2)
P_stoch[2,:] = @. 1/n_bar*(cϵ2 * (k_grid/k_nl)^2)
return P_stoch
end
-function get_stoch_terms(cϵ0, cϵ1, cϵ2, n_bar, k_grid::Array, k_nl=0.7)
+function get_stoch_terms(cϵ0, cϵ1, cϵ2, n_bar, k_grid::Array; k_nl=0.7)
P_stoch_0 = @. 1/n_bar*(cϵ0 + cϵ1*(k_grid/k_nl)^2)
P_stoch_2 = @. 1/n_bar*(cϵ2 * (k_grid/k_nl)^2)
return P_stoch_0, P_stoch_2
@@ -64,7 +64,7 @@ end
function get_stoch_terms_binned_efficient(cϵ0, cϵ1, cϵ2, n_bar, k_edges::Array; k_nl=0.7)
n_bins = length(k_edges)-1
- P_stoch_0_c, _ = get_stoch_terms(cϵ0, cϵ1, cϵ2, n_bar, k_edges, k_nl)
+ P_stoch_0_c, _ = get_stoch_terms(cϵ0, cϵ1, cϵ2, n_bar, k_edges; k_nl = k_nl)
mytype = typeof(P_stoch_0_c[1])
P_stoch_0 = zeros(mytype, n_bins)
P_stoch_2 = zeros(mytype, n_bins)
@@ -78,7 +78,7 @@ end
function get_stoch_terms_binned(cϵ0, cϵ1, cϵ2, n_bar, k_edges::Array; k_nl=0.7, nk = 30)
n_bins = length(k_edges)-1
- P_stoch_0_c, _ = get_stoch_terms(cϵ0, cϵ1, cϵ2, n_bar, k_edges, k_nl)
+ P_stoch_0_c, _ = get_stoch_terms(cϵ0, cϵ1, cϵ2, n_bar, k_edges; k_nl = k_nl)
mytype = typeof(P_stoch_0_c[1])
P_stoch_0 = zeros(mytype, n_bins)
P_stoch_2 = zeros(mytype, n_bins)
diff --git a/src/neural_networks.jl b/src/neural_networks.jl
index f172451..cbbb3d9 100644
--- a/src/neural_networks.jl
+++ b/src/neural_networks.jl
@@ -1,26 +1,3 @@
-function _maximin_input!(x, in_MinMax)
- for i in eachindex(x)
- x[i] -= in_MinMax[i,1]
- x[i] /= (in_MinMax[i,2]-in_MinMax[i,1])
- end
-end
-
-function _inv_maximin_output!(x, out_MinMax)
- for i in eachindex(x)
- x[i] *= (out_MinMax[i,2]-out_MinMax[i,1])
- x[i] += out_MinMax[i,1]
- end
-end
-
-abstract type AbstractTrainedEmulators end
-#TODO: right now we support only SimpleChains emulators, but in the near future we may
-#support as well Flux and Lux. This is the reason behind the struct SimpleChainsEmulator
-
-@kwdef mutable struct SimpleChainsEmulator <: AbstractTrainedEmulators
- Architecture
- Weights
-end
-
abstract type AbstractComponentEmulators end
@kwdef mutable struct P11Emulator <: AbstractComponentEmulators
@@ -46,9 +23,9 @@ end
function get_component(input_params, comp_emu::AbstractComponentEmulators)
input = deepcopy(input_params)
- _maximin_input!(input, comp_emu.InMinMax)
- output = Array(_run_emulator(input, comp_emu.TrainedEmulator))
- _inv_maximin_output!(output, comp_emu.OutMinMax)
+ maximin_input!(input, comp_emu.InMinMax)
+ output = Array(run_emulator(input, comp_emu.TrainedEmulator))
+ inv_maximin_output!(output, comp_emu.OutMinMax)
As = exp(input_params[1])*1e-10
output .*= As
return reshape(output, Int(length(output)/length(comp_emu.kgrid)), :)
@@ -56,18 +33,14 @@ end
function get_component(input_params, comp_emu::PloopEmulator)
input = deepcopy(input_params)
- _maximin_input!(input, comp_emu.InMinMax)
- output = Array(_run_emulator(input, comp_emu.TrainedEmulator))
- _inv_maximin_output!(output, comp_emu.OutMinMax)
+ maximin_input!(input, comp_emu.InMinMax)
+ output = Array(run_emulator(input, comp_emu.TrainedEmulator))
+ inv_maximin_output!(output, comp_emu.OutMinMax)
As = exp(input_params[1])*1e-10
output .*= As^2
return reshape(output, Int(length(output)/length(comp_emu.kgrid)), :)
end
-function _run_emulator(input, trained_emulator::SimpleChainsEmulator)
- return trained_emulator.Architecture(input, trained_emulator.Weights)
-end
-
abstract type AbstractPℓEmulators end
@kwdef mutable struct PℓEmulator <: AbstractPℓEmulators
diff --git a/src/projection.jl b/src/projection.jl
index 843b4ea..a6ceca3 100644
--- a/src/projection.jl
+++ b/src/projection.jl
@@ -3,7 +3,7 @@ function _Pkμ(k, μ, Int_Mono, Int_Quad, Int_Hexa)
end
function _k_true(k_o, μ_o, q_perp, F)
- return k_o/q_perp*sqrt(1+μ_o^2*(1/F^2-1))
+ return @. k_o/q_perp*sqrt(1+μ_o^2*(1/F^2-1))
end
function _μ_true(μ_o, F)
@@ -18,9 +18,11 @@ function _P_obs(k_o, μ_o, q_par, q_perp, Int_Mono, Int_Quad, Int_Hexa)
end
function interp_Pℓs(Mono_array, Quad_array, Hexa_array, k_grid)
- Int_Mono = CubicSpline(Mono_array, k_grid)
- Int_Quad = CubicSpline(Quad_array, k_grid)
- Int_Hexa = CubicSpline(Hexa_array, k_grid)
+ #extrapolation might introduce some errors ar high k, when q << 1.
+ #maybe we should implement a log extrapolation?
+ Int_Mono = QuadraticSpline(Mono_array, k_grid; extrapolate = true)
+ Int_Quad = QuadraticSpline(Quad_array, k_grid; extrapolate = true)
+ Int_Hexa = QuadraticSpline(Hexa_array, k_grid; extrapolate = true)
return Int_Mono, Int_Quad, Int_Hexa
end
@@ -29,8 +31,8 @@ function k_projection(k_projection, Mono_array, Quad_array, Hexa_array, k_grid)
return int_Mono.(k_projection), int_Quad.(k_projection), int_Hexa.(k_projection)
end
-function apply_AP_check(k_grid, int_Mono::CubicSpline, int_Quad::CubicSpline,
- int_Hexa::CubicSpline, q_par, q_perp)
+function apply_AP_check(k_grid, int_Mono::QuadraticSpline, int_Quad::QuadraticSpline,
+ int_Hexa::QuadraticSpline, q_par, q_perp)
nk = length(k_grid)
result = zeros(3, nk)
ℓ_array = [0,2,4]
@@ -58,13 +60,21 @@ end
function _mygemmavx(A, B, C)
Dm = zero(eltype(C))
- @turbo for n ∈ axes(A,1)
+ @tturbo for n ∈ axes(A,1)
Dm += A[n] * B[n] * C[n]
end
return Dm
end
-function apply_AP(k_grid, int_Mono::CubicSpline, int_Quad::CubicSpline, int_Hexa::CubicSpline,
+function _mygemm(A, B, C)
+ Dm = zero(eltype(C))
+ for n ∈ axes(A,1)
+ Dm += A[n] * B[n] * C[n]
+ end
+ return Dm
+end
+
+function apply_AP(k_grid, int_Mono::QuadraticSpline, int_Quad::QuadraticSpline, int_Hexa::QuadraticSpline,
q_par, q_perp)
nk = length(k_grid)
n_GL_points = 5
@@ -95,22 +105,82 @@ function apply_AP(k_grid, int_Mono::CubicSpline, int_Quad::CubicSpline, int_Hexa
return result
end
-function q_perp_par(z, ΩM_ref, w0_ref, wa_ref, ΩM_true, w0_true, wa_true)
+function _stoch_obs(k_o, μ_o, q_par, q_perp, n_bar, cϵ0, cϵ1, cϵ2, k_nl)
+ F = q_par/q_perp
+ k_t = _k_true(k_o, μ_o, q_perp, F)
+ μ_t = _μ_true(μ_o, F)
+ return _stoch_kμ(k_t, μ_t, n_bar, cϵ0, cϵ1, cϵ2, k_nl)/(q_par*q_perp^2)
+end
+
+function _k_grid_over_nl(k_grid, k_nl)
+ return @. (k_grid/k_nl)^2
+ end
+
+ function _stoch_kμ(k_grid, μ, n_bar, cϵ0, cϵ1, cϵ2, k_nl)
+ return (cϵ0 * Pl(μ, 0) .+ Effort._k_grid_over_nl(k_grid, k_nl) .* (cϵ1 * Pl(μ, 0) +
+ cϵ2 * Pl(μ, 2)) ) ./ n_bar
+end
+
+function get_stochs_AP(k_grid, q_par, q_perp, n_bar, cϵ0, cϵ1, cϵ2; k_nl = 0.7, n_GL_points = 4)
+ nk = length(k_grid)
+ #TODO: check that the extrapolation does not create problems. Maybe logextrap?
+ nodes, weights = @memoize gausslobatto(n_GL_points*2)
+ #since the integrand is symmetric, we are gonna use only half of the points
+ μ_nodes = nodes[1:n_GL_points]
+ μ_weights = weights[1:n_GL_points]
+ result = zeros(2, nk)
+
+ Pl_array = zeros(2, n_GL_points)
+ Pl_array[1,:] .= Pl.(μ_nodes, 0)
+ Pl_array[2,:] .= Pl.(μ_nodes, 2)
+
+ temp = zeros(n_GL_points, nk)
+
+ for i in 1:n_GL_points
+ temp[i,:] = _stoch_obs(k_grid, μ_nodes[i], q_par, q_perp, n_bar, cϵ0, cϵ1, cϵ2,
+ k_nl)
+ end
+
+ multipole_weight = [2*0+1, 2*2+1]
+
+ @turbo for i in 1:2
+ for j in 1:nk
+ for k in 1:n_GL_points
+ result[i, j] += μ_weights[k] * temp[k,j] * Pl_array[i, k] *
+ multipole_weight[i]
+ end
+ end
+ end
+
+ return result
+end
+
+function q_par_perp(z, ΩM_ref, w0_ref, wa_ref, ΩM_true, w0_true, wa_true)
E_ref = _E_z(z, ΩM_ref, w0_ref, wa_ref)
E_true = _E_z(z, ΩM_true, w0_true, wa_true)
d̃A_ref = _d̃A_z(z, ΩM_ref, w0_ref, wa_ref)
d̃A_true = _d̃A_z(z, ΩM_true, w0_true, wa_true)
- q_perp = E_ref/E_true
- q_par = d̃A_true/d̃A_ref
+ q_perp = E_true/E_ref
+ q_par = d̃A_ref/d̃A_true
+ return q_par, q_perp
+end
+
+function q_par_perp(z, ΩM_ref, w0_ref, wa_ref, E_true, d̃A_true)
+ E_ref = _E_z(z, ΩM_ref, w0_ref, wa_ref)
+
+ d̃A_ref = _d̃A_z(z, ΩM_ref, w0_ref, wa_ref)
+
+ q_perp = E_true/E_ref
+ q_par = d̃A_ref/d̃A_true
return q_perp, q_par
end
function apply_AP(k_grid, Mono_array::Array, Quad_array::Array, Hexa_array::Array, z, ΩM_ref,
w0_ref, wa_ref, ΩM_true, w0_true, wa_true)
- q_perp, q_par = q_perp_par(z, ΩM_ref, w0_ref, wa_ref, ΩM_true, w0_true, wa_true)
+ q_par, q_perp = q_par_perp(z, ΩM_ref, w0_ref, wa_ref, ΩM_true, w0_true, wa_true)
int_Mono, int_Quad, int_Hexa = interp_Pℓs(Mono_array, Quad_array, Hexa_array, k_grid)
return apply_AP(k_grid, int_Mono, int_Quad, int_Hexa, q_par, q_perp)
end
@@ -118,7 +188,7 @@ end
function apply_AP(k_grid_AP, k_grid, Mono_array::Array, Quad_array::Array, Hexa_array::Array, z,
ΩM_ref, w0_ref, wa_ref, ΩM_true, w0_true, wa_true)
- q_perp, q_par = q_perp_par(z, ΩM_ref, w0_ref, wa_ref, ΩM_true, w0_true, wa_true)
+ q_par, q_perp = q_par_perp(z, ΩM_ref, w0_ref, wa_ref, ΩM_true, w0_true, wa_true)
int_Mono, int_Quad, int_Hexa = interp_Pℓs(Mono_array, Quad_array, Hexa_array, k_grid)
return apply_AP(k_grid_AP, int_Mono, int_Quad, int_Hexa, q_par, q_perp)
@@ -141,3 +211,22 @@ function apply_AP(k_grid_AP, k_interp, Mono_array::Array, Quad_array::Array, Hex
int_Mono, int_Quad, int_Hexa = interp_Pℓs(Mono_array, Quad_array, Hexa_array, k_interp)
return apply_AP(k_grid_AP, int_Mono, int_Quad, int_Hexa, q_par, q_perp)
end
+
+function window_convolution(W,v)
+ a,b,c,d = size(W)
+ result = zeros(a,c)
+ window_convolution!(result,W,v)
+ return result
+end
+
+function window_convolution!(result, W, v)
+ @turbo for i ∈ axes(W,1), k ∈ axes(W,3)
+ c = zero(eltype(v))
+ for l ∈ axes(W,4)
+ for j ∈ axes(W,2)
+ c += W[i,j,k,l] * v[j,l]
+ end
+ end
+ result[i,k] = c
+ end
+end
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..ea13341
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,10 @@
+[deps]
+AbstractCosmologicalEmulators = "c83c1981-e5c4-4837-9eb8-c9b1572acfc6"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
+SimpleChains = "de6bee2f-e2f4-4ec7-b6ed-219cc6f6e9e5"
+Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
diff --git a/test/runtests.jl b/test/runtests.jl
new file mode 100644
index 0000000..9ca8a72
--- /dev/null
+++ b/test/runtests.jl
@@ -0,0 +1,54 @@
+using Test
+using NPZ
+using SimpleChains
+using Static
+using Effort
+using ForwardDiff
+using Zygote
+using FiniteDiff
+using SciMLSensitivity
+
+mlpd = SimpleChain(
+ static(6),
+ TurboDense(tanh, 64),
+ TurboDense(tanh, 64),
+ TurboDense(tanh, 64),
+ TurboDense(tanh, 64),
+ TurboDense(tanh, 64),
+ TurboDense(identity, 40)
+)
+
+k_test = Array(LinRange(0,200, 40))
+weights = SimpleChains.init_params(mlpd)
+inminmax = rand(6,2)
+outminmax = rand(40,2)
+a, Ωγ0, Ωm0, mν, h, w0, wa = [1., 1e-5, 0.3, 0.06, 0.67, -1.1, 0.2]
+z = [0.4, 0.5]
+
+emu = Effort.SimpleChainsEmulator(Architecture = mlpd, Weights = weights)
+
+effort_emu = Effort.P11Emulator(TrainedEmulator = emu, kgrid=k_test, InMinMax = inminmax,
+ OutMinMax = outminmax)
+
+x = [Ωm0, h, mν, w0, wa]
+
+function D_z_x(z, x)
+ Ωm0, h, mν, w0, wa = x
+ sum(Effort._D_z(z, Ωm0, h; mν =mν, w0=w0, wa=wa))
+end
+
+function f_z_x(z, x)
+ Ωm0, h, mν, w0, wa = x
+ sum(Effort._f_z(z, Ωm0, h; mν =mν, w0=w0, wa=wa))
+end
+
+@testset "Effort tests" begin
+ @test isapprox(Effort._H_a(a, Ωγ0, Ωm0, mν, h, w0, wa), h*100)
+ @test isapprox(Effort._E_a(a, Ωm0, h), 1.)
+ @test isapprox(Effort._D_z_old(z, Ωm0, h), Effort._D_z(z, Ωm0, h), rtol=1e-9)
+ @test isapprox(Effort._f_z_old(0.4, Ωm0, h), Effort._f_z(0.4, Ωm0, h)[1], rtol=1e-9)
+ @test isapprox(Zygote.gradient(x->D_z_x(z, x), x)[1], ForwardDiff.gradient(x->D_z_x(z, x), x), rtol=1e-5)
+ @test isapprox(FiniteDiff.finite_difference_gradient(x->D_z_x(z, x), x), ForwardDiff.gradient(x->D_z_x(z, x), x), rtol=1e-5)
+ @test isapprox(Zygote.gradient(x->f_z_x(z, x), x)[1], ForwardDiff.gradient(x->f_z_x(z, x), x), rtol=1e-5)
+ @test isapprox(FiniteDiff.finite_difference_gradient(x->f_z_x(z, x), x), ForwardDiff.gradient(x->f_z_x(z, x), x), rtol=1e-4)
+end