From 484cfe89922183b7e372e1cc47c3f8566176cb89 Mon Sep 17 00:00:00 2001 From: Matteo Foglieni <79974922+cosmofico97@users.noreply.github.com> Date: Wed, 21 Feb 2024 14:52:50 +0100 Subject: [PATCH] test new brand --- ipynbs/gapse-computations.jl | 338 +++++++++++++++++++++++++++++++++++ 1 file changed, 338 insertions(+) create mode 100644 ipynbs/gapse-computations.jl diff --git a/ipynbs/gapse-computations.jl b/ipynbs/gapse-computations.jl new file mode 100644 index 00000000..0a9de5d1 --- /dev/null +++ b/ipynbs/gapse-computations.jl @@ -0,0 +1,338 @@ +# -*- encoding: utf-8 -*- +# +# This file is part of GaPSE +# Copyright (C) 2022 Matteo Foglieni +# +# GaPSE is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# GaPSE is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GaPSE. If not, see . +# + +using Pkg + +PATH_TO_GAPSE ="../" +#PATH_TO_GAPSE = "/dss/dsshome1/08/di75tom/julia_with_salvo/GaPSE.jl/" +Pkg.activate(normpath(PATH_TO_GAPSE)) # or @__DIR__ +using GaPSE + + +const global N_LOB = 1000 +const global N_TRAP = 1000 +const global ALG = :lobatto +const global N_CHIS = 100 +const global N_CHIS_2 = 50 +const global SUIT_SAMPLING = true +const global LENGTH = 500 +const global OM = true +COMPUTE_XIS_GNC = true # set this to true if you want to compute the TPCFs of the GNC! +# This is the directory name where to put the files computed in the +# following steps; makes sure that it's name ends with "/" ! +#DIR = "/dss/dsshome1/08/di75tom/julia_with_salvo/" +DIR="./" +@assert isdir(DIR) "ERROR: DIR=$DIR DOESN'T EXIST!!!" + +using DelimitedFiles, Dierckx, Printf + +FILE_NAME = split(PROGRAM_FILE, "/")[end] + +#main(x::Union{String, Float64, Int64}...) = main([string(var) for var in [x...]]) +function main() + + # run in the command line: + # + # $ julia gapse-computations.jl + + + ########## This is the basic configuration for a Cosmology ########### + + ps_kwargs(alg::Symbol=:fftlog) = alg == :twofast ? + Dict( + :alg => :twofast, :epl => true, :pr => false, + :N_left => 12, :N_right => 12, + :p0_left => [-2.0, 1.0], :p0_right => [-2.0, 1.0], + :int_s_min => 1e0, :int_s_max => 1200.0, + :cut_first_n => 0, :cut_last_n => 0 + ) : alg == :fftlog ? + Dict( + :alg => :fftlog, :pr => true, :ν => 1.5, + :n_extrap_low => 300, :n_extrap_high => 300, + :n_pad => 500, :cut_first_n => 0, :cut_last_n => 0, + ) : throw(AssertionError("alg = :fftlog (recommended) or alg = :twofast !")) + tf = :fftlog + + FILE_F_MAP = PATH_TO_GAPSE * "data/NEW_F_pi2.txt" + #FILE_F_MAP = PATH_TO_GAPSE * "data/F_REFERENCE_pi2.txt"; + #= + kwargs_map_F_hcub = Dict( + :θ_max => π / 2.0, :tolerance => 1e-10, + :rtol => 1e-2, :atol => 1e-3, :pr => true, + ); + + kwargs_map_F_trap = Dict( + :θ_max => π / 2.0, :tolerance => 1e-10, + :N => 1000, :pr => true, + ); + + xs = [x for x in 0:0.02:5] + μs = union( + [μ for μ in range(-1.0, -0.98, length = 50)], + [μ for μ in range(-0.98, 0.98, length = 102)], + [μ for μ in range(0.98, 1.0, length = 50)]); + GaPSE.print_map_F(FILE_F_MAP, xs, μs; + alg = :trap, Fmap_opts = kwargs_map_F_trap # we recommend to use :trap + #alg = :hcub, Fmap_opts = kwargs_map_F_hcub # but you can use also :hcub if you prefer + ) + =# + + z_min, z_max, θ_max = 1.0, 1.5, π / 2.0 + + FILE_BACKGROUND = PATH_TO_GAPSE * "data/WideA_ZA_background.dat" + + params = GaPSE.CosmoParams(z_min, z_max, θ_max; + Ω_b=0.0489, Ω_cdm=0.251020, + h_0 = 0.70, s_lim = 1e-2, z_spline_lim = 1000.0, + b1=1.5, s_b1=0.0, 𝑓_evo1=0.0, + b2 = nothing, s_b2 = nothing, 𝑓_evo2 = nothing, + IPS_opts=Dict( + :fit_left_min => 1e-6, :fit_left_max => 3e-6, + :fit_right_min => 1e1, :fit_right_max => 2e1), + IPSTools_opts=Dict( + :N => 1024, :fit_min => 0.05, :fit_max => 0.5, + :con => true, :k_min => 1e-8, :k_max => 10.0) + ) + + #FILE_F_MAP = PATH_TO_GAPSE * "data/NEW_F_pi2.txt"; + FILE_IF_MAP = PATH_TO_GAPSE * "data/NEW_IntegrF_pi2_z115.txt" + + #FILE_F_MAP = PATH_TO_GAPSE*"data/F_REFERENCE_pi2.txt"; + #FILE_IF_MAP = PATH_TO_GAPSE*"data/IntegrF_REFERENCE_pi2_z115.txt"; + + #= + calc_μs = union( + [μ for μ in range(-1.0, -0.98, length = 50)], + [μ for μ in range(-0.98, 0.98, length = 102)], + [μ for μ in range(0.98, 1.0, length = 50)]); + + GaPSE.print_map_IntegratedF( + z_min, z_max, calc_μs, + FILE_F_MAP, FILE_IF_MAP, + FILE_BACKGROUND; + alg = :trap, N_ss = 200, m = 2.1, + Dict( + :llim => nothing, :rlim => nothing, + :rtol => 1e-2, :atol => 0.0, + :N => 1000, :pr => true, + )... + ) + =# + + + FILE_PS = PATH_TO_GAPSE * "test/datatest/file_pk.txt" + cosmo = GaPSE.Cosmology(params, FILE_BACKGROUND, FILE_PS, FILE_F_MAP, FILE_IF_MAP) + + GaPSE.parameters_used(stdout, cosmo) + + + + + + + + ######## Computation of the GNC TPCFs and Newtonian TPCF (and their PS) ######### + + name_xis_GNC_L0_noF_noobsvel_file = DIR * "xis_GNC_L0_noF_noobsvel_GenWin.txt" + name_xis_GNC_L1_noF_noobsvel_file = DIR * "xis_GNC_L1_noF_noobsvel_GenWin.txt" + name_xis_GNC_L2_noF_noobsvel_file = DIR * "xis_GNC_L2_noF_noobsvel_GenWin.txt" + name_xis_GNC_L3_noF_noobsvel_file = DIR * "xis_GNC_L3_noF_noobsvel_GenWin.txt" + name_xis_GNC_L4_noF_noobsvel_file = DIR * "xis_GNC_L4_noF_noobsvel_GenWin.txt" + + + if (COMPUTE_XIS_GNC == true) + GaPSE.print_map_sum_ξ_GNC_multipole( + cosmo, name_xis_GNC_L0_noF_noobsvel_file, + 10 .^ range(0, log10(2 * cosmo.s_max), length=LENGTH); + use_windows=false, L=0, alg=ALG, obs=:noobsvel, + single=true, enhancer=1e8, + N_trap=N_TRAP, N_lob=N_LOB, atol_quad=0.0, rtol_quad=1e-2, + N_χs=N_CHIS, N_χs_2=N_CHIS_2, suit_sampling=SUIT_SAMPLING + ) + end + + if (COMPUTE_XIS_GNC == true && OM == false) + GaPSE.print_map_sum_ξ_GNC_multipole( + cosmo, name_xis_GNC_L1_noF_noobsvel_file, + 10 .^ range(0, log10(2 * cosmo.s_max), length=LENGTH); + use_windows=false, L=1, alg=ALG, obs=:noobsvel, + single=true, enhancer=1e8, + N_trap=N_TRAP, N_lob=N_LOB, atol_quad=0.0, rtol_quad=1e-2, + N_χs=N_CHIS, N_χs_2=N_CHIS_2, suit_sampling=SUIT_SAMPLING + ) + end + + if (COMPUTE_XIS_GNC == true && OM == false) + GaPSE.print_map_sum_ξ_GNC_multipole( + cosmo, name_xis_GNC_L2_noF_noobsvel_file, + 10 .^ range(0, log10(2 * cosmo.s_max), length=LENGTH); + use_windows=false, L=2, alg=ALG, obs=:noobsvel, + single=true, enhancer=1e8, + N_trap=N_TRAP, N_lob=N_LOB, atol_quad=0.0, rtol_quad=1e-2, + N_χs=N_CHIS, N_χs_2=N_CHIS_2, suit_sampling=SUIT_SAMPLING + ) + end + + if (COMPUTE_XIS_GNC == true && OM == false) + GaPSE.print_map_sum_ξ_GNC_multipole( + cosmo, name_xis_GNC_L3_noF_noobsvel_file, + 10 .^ range(0, log10(2 * cosmo.s_max), length=LENGTH); + use_windows=false, L=3, alg=ALG, obs=:noobsvel, + single=true, enhancer=1e8, + N_trap=N_TRAP, N_lob=N_LOB, atol_quad=0.0, rtol_quad=1e-2, + N_χs=N_CHIS, N_χs_2=N_CHIS_2, suit_sampling=SUIT_SAMPLING + ) + end + + if (COMPUTE_XIS_GNC == true && OM == false) + GaPSE.print_map_sum_ξ_GNC_multipole( + cosmo, name_xis_GNC_L4_noF_noobsvel_file, + 10 .^ range(0, log10(2 * cosmo.s_max), length=LENGTH); + use_windows=false, L=4, alg=ALG, obs=:noobsvel, + single=true, enhancer=1e8, + N_trap=N_TRAP, N_lob=N_LOB, atol_quad=0.0, rtol_quad=1e-2, + N_χs=N_CHIS, N_χs_2=N_CHIS_2, suit_sampling=SUIT_SAMPLING + ) + end + + + ######## create_file_for_XiMultipoles function ######### + + + + if OM == false + name_file_ximultipoles = DIR * "xis_GNC_L01234_noF_noobsvel.txt" + + vector_xis_names = [ + name_xis_GNC_L0_noF_noobsvel_file, + name_xis_GNC_L1_noF_noobsvel_file, + name_xis_GNC_L2_noF_noobsvel_file, + name_xis_GNC_L3_noF_noobsvel_file, + name_xis_GNC_L4_noF_noobsvel_file, + ] + + GaPSE.create_file_for_XiMultipoles(name_file_ximultipoles, vector_xis_names, 2, "GNC") + + #run(`head -n 18 $(name_file_ximultipoles)`) + + ps_kwargs_2(alg::Symbol=:fftlog) = alg == :twofast ? + Dict( + :alg => :twofast, :epl => true, :pr => false, + :N_left => 12, :N_right => 12, + :p0_left => [-2.0, 1.0], :p0_right => [-2.0, 1.0], + :int_s_min => 1e0, :int_s_max => 1200.0, + :cut_first_n => 0, :cut_last_n => 0 + ) : alg == :fftlog ? + Dict( + :alg => :fftlog, :pr => true, :ν => 1.5, + :n_extrap_low => 0, :n_extrap_high => 0, + :n_pad => 500, :cut_first_n => 0, :cut_last_n => 0, + ) : throw(AssertionError("alg = :fftlog (recommended) or alg = :twofast !")) + tf = :fftlog + + ks, pks = GaPSE.PS_multipole_GenWin(name_file_ximultipoles, name_file_Qmultipoles; L=0, ps_kwargs_2(:fftlog)...) + + out = DIR * "results.txt" + isfile(out) && run(`rm $out`) + open(out, "w") do io + println(io, BRAND) + + println(io, "# Power Spectrum GenWin monopole computation.") + println(io, "#\n# For this PS_multipole computation we set: ") + println(io, "# \t algorithm chosen for the computation: :$ALG") + println(io, "# \t N_LOB = $N_LOB \t N_TRAP = $N_TRAP \t LENGTH = $LENGTH") + println(io, "# \t N_LOB = $N_LOB \t N_TRAP = $N_TRAP \t LENGTH = $LENGTH") + #println(io, "# computational time needed (in s) : $(@sprintf("%.4f", time_2-time_1))") + print(io, "# kwards passed to \"print_PS_multipole\": ") + + if isempty(ps_kwargs_2) + println(io, "none") + else + print(io, "\n") + for key in keys(ps_kwargs_2) + println(io, "# \t\t$(key) = $(ps_kwargs_2[key])") + end + end + println(io, "# ") + println(io, "# k [h_0/Mpc] \t \t P [(Mpc/h_0)^3]") + for (k, pk) in zip(ks, pks) + println(io, "$k \t " * GaPSE.number_to_string(pk)) + end + end + end + + #= + GaPSE.print_all_PS_multipole(DIR * name_xis_GNC_L0_noF_noobsvel_file, + DIR * name_ps_GNC_L0_noF_noobsvel_file, "GNC"; L=0, ps_kwargs(tf)...) + ks_GNC_L0_noF_noobsvel, pks_sum_GNC_L0_noF_noobsvel, pks_all_GNC_L0_noF_noobsvel = + GaPSE.readxyall(DIR * name_ps_GNC_L0_noF_noobsvel_file) + spline_GNCsum_pks_L0_noF_noobsvel = Spline1D(ks_GNC_L0_noF_noobsvel, pks_sum_GNC_L0_noF_noobsvel; bc="error") + + GaPSE.print_all_PS_multipole(DIR * name_xis_GNC_L0_noF_noobs_file, + DIR * name_ps_GNC_L0_noF_noobs_file, "GNC"; L=0, ps_kwargs(tf)...) + ks_GNC_L0_noF_noobs, pks_sum_GNC_L0_noF_noobs, pks_all_GNC_L0_noF_noobs = + GaPSE.readxyall(DIR * name_ps_GNC_L0_noF_noobs_file) + spline_GNCsum_pks_L0_noF_noobs = Spline1D(ks_GNC_L0_noF_noobs, pks_sum_GNC_L0_noF_noobs; bc="error") + + GaPSE.print_all_PS_multipole(DIR * name_xis_GNC_L0_withF_noobsvel_file, + DIR * name_ps_GNC_L0_withF_noobsvel_file, "GNC"; L=0, ps_kwargs(tf)...) + ks_GNC_L0_withF_noobsvel, pks_sum_GNC_L0_withF_noobsvel, pks_all_GNC_L0_withF_noobsvel = + GaPSE.readxyall(DIR * name_ps_GNC_L0_withF_noobsvel_file) + spline_GNCsum_pks_L0_withF_noobsvel = Spline1D(ks_GNC_L0_withF_noobsvel, pks_sum_GNC_L0_withF_noobsvel; bc="error") + + GaPSE.print_all_PS_multipole(DIR * name_xis_GNC_L0_withF_noobs_file, + DIR * name_ps_GNC_L0_withF_noobs_file, "GNC"; L=0, ps_kwargs(tf)...) + ks_GNC_L0_withF_noobs, pks_sum_GNC_L0_withF_noobs, pks_all_GNC_L0_withF_noobs = + GaPSE.readxyall(DIR * name_ps_GNC_L0_withF_noobs_file) + spline_GNCsum_pks_L0_withF_noobs = Spline1D(ks_GNC_L0_withF_noobs, pks_sum_GNC_L0_withF_noobs; bc="error") + + + ks_Newtonian_L0_noF_noobsvel, pks_Newtonian_L0_noF_noobsvel = + ks_GNC_L0_noF_noobsvel, pks_all_GNC_L0_noF_noobsvel[GaPSE.INDEX_GR_EFFECT_GNC["auto_newton"]] + spline_Newt_pks_L0_noF_noobsvel = Spline1D(ks_Newtonian_L0_noF_noobsvel, pks_Newtonian_L0_noF_noobsvel; bc="error") + ks_Newtonian_L0_noF_noobs, pks_Newtonian_L0_noF_noobs = + ks_GNC_L0_noF_noobs, pks_all_GNC_L0_noF_noobs[GaPSE.INDEX_GR_EFFECT_GNC["auto_newton"]] + spline_Newt_pks_L0_noF_noobs = Spline1D(ks_Newtonian_L0_noF_noobs, pks_Newtonian_L0_noF_noobs; bc="error") + + ks_Newtonian_L0_withF_noobsvel, pks_Newtonian_L0_withF_noobsvel = + ks_GNC_L0_withF_noobsvel, pks_all_GNC_L0_withF_noobsvel[GaPSE.INDEX_GR_EFFECT_GNC["auto_newton"]] + spline_Newt_pks_L0_withF_noobsvel = Spline1D(ks_Newtonian_L0_withF_noobsvel, pks_Newtonian_L0_withF_noobsvel; bc="error") + =# + + + + + + +end + + + + + + + +if (ARGS == String[]) + #println("\nwithout input arguments/commands, show this help message and exit\n") + #main(["--help"]) + main() +else + #main(ARGS) + return 0 +end