From 6a9936afbefee312b1cad66c38f9012d283bdd72 Mon Sep 17 00:00:00 2001 From: Marco Bonici Date: Mon, 8 Mar 2021 12:23:05 +0100 Subject: [PATCH 1/9] Added @avx to the lensing efficiency evaluation (25%speedup) --- input_files/Angular.json | 2 +- src/WeightFunctions.jl | 15 +++++++++++---- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/input_files/Angular.json b/input_files/Angular.json index 7cb85510..a005ae11 100644 --- a/input_files/Angular.json +++ b/input_files/Angular.json @@ -1,6 +1,6 @@ { "Lensing": { - "present": false + "present": true }, "PhotometricGalaxy":{ "present": true, diff --git a/src/WeightFunctions.jl b/src/WeightFunctions.jl index afaf2a34..0bbcb70f 100644 --- a/src/WeightFunctions.jl +++ b/src/WeightFunctions.jl @@ -127,14 +127,21 @@ function ComputeLensingEfficiencyOverGridCustom( BackgroundQuantities::BackgroundQuantities, w0waCDMCosmology::w0waCDMCosmologyStruct) WLWeightFunction.LensingEfficiencyArray .*= 0 - #TODO add avx here to improve performance - for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1 + elimination_matrix = zeros(length(CosmologicalGrid.ZArray), + length(CosmologicalGrid.ZArray)) + for i in 1:length(CosmologicalGrid.ZArray) + for j in i:length(CosmologicalGrid.ZArray) + elimination_matrix[i,j] = 1 + end + end + @avx for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1 for idx_ZArray in 1:length(CosmologicalGrid.ZArray) - for idx_ZArrayInt in idx_ZArray:length(CosmologicalGrid.ZArray) + for idx_ZArrayInt in 1:length(CosmologicalGrid.ZArray) WLWeightFunction.LensingEfficiencyArray[idx_ZBinArray, idx_ZArray] += ConvolvedDensity.DensityGridArray[idx_ZBinArray, idx_ZArrayInt] * (1 - BackgroundQuantities.rZArray[idx_ZArray] / - BackgroundQuantities.rZArray[idx_ZArrayInt]) + BackgroundQuantities.rZArray[idx_ZArrayInt]) * + elimination_matrix[idx_ZArray, idx_ZArrayInt] end end end From f750c93c7380db260d72935d80176994dbf0bb3a Mon Sep 17 00:00:00 2001 From: Marco Bonici Date: Mon, 8 Mar 2021 17:59:06 +0100 Subject: [PATCH 2/9] Added Simpson integration to angular coefficients --- src/AngularCoefficients.jl | 4 +++- src/MathUtils.jl | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 37 insertions(+), 1 deletion(-) diff --git a/src/AngularCoefficients.jl b/src/AngularCoefficients.jl index bac6a860..97f6384c 100644 --- a/src/AngularCoefficients.jl +++ b/src/AngularCoefficients.jl @@ -64,6 +64,7 @@ function ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients, c_0 = 2.99792458e5 #TODO: find a package containing the exact value of #physical constants involved in calculations Integrand = similar(AngularCoefficients.AngularCoefficientsArray) .*0 + Simpson_weights = SimpsonWeightArray(length(CosmologicalGrid.ZArray)) @avx for i ∈ axes(AngularCoefficients.AngularCoefficientsArray,2), j ∈ axes(AngularCoefficients.AngularCoefficientsArray,3), l ∈ axes(AngularCoefficients.AngularCoefficientsArray,1) @@ -73,7 +74,8 @@ function ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients, WeightFunctionB.WeightFunctionArray[j, z] / (BackgroundQuantities.HZArray[z] * BackgroundQuantities.rZArray[z]^2) * - PowerSpectrum.InterpolatedPowerSpectrum[l,z] + PowerSpectrum.InterpolatedPowerSpectrum[l,z] * + Simpson_weights[z] end end Integrand .*= (last(CosmologicalGrid.ZArray)- diff --git a/src/MathUtils.jl b/src/MathUtils.jl index 68009f65..2de6ee63 100644 --- a/src/MathUtils.jl +++ b/src/MathUtils.jl @@ -50,3 +50,37 @@ function CustomRegression(x::Vector{Float64}, y::Vector{Float64}) c = y_mean - m* x_mean return c, m end + +""" + SimpsonWeightArray(n::Int64) + +This function evaluates an array with the Simpson weight, for Simpson +Integration +""" +function SimpsonWeightArray(n::Int64) + number_intervals = floor((n-1)/2) + weight_array = zeros(n) + if n == number_intervals*2+1 + for i in 1:number_intervals + weight_array[Int((i-1)*2+1)] += 1/3 + weight_array[Int((i-1)*2+2)] += 4/3 + weight_array[Int((i-1)*2+3)] += 1/3 + end + else + weight_array[1] +=0.5 + weight_array[2] +=0.5 + for i in 1:number_intervals + weight_array[Int((i-1)*2+1)+1] += 1/3 + weight_array[Int((i-1)*2+2)+1] += 4/3 + weight_array[Int((i-1)*2+3)+1] += 1/3 + end + weight_array[length(weight_array)] +=0.5 + weight_array[length(weight_array)-1] +=0.5 + for i in 1:number_intervals + weight_array[Int((i-1)*2+1)] += 1/3 + weight_array[Int((i-1)*2+2)] += 4/3 + weight_array[Int((i-1)*2+3)] += 1/3 + end + end + return weight_array ./2 +end From 54565a27b354850a7a71151c0cd9a026268d0a13 Mon Sep 17 00:00:00 2001 From: Marco Bonici Date: Mon, 8 Mar 2021 21:28:49 +0100 Subject: [PATCH 3/9] Added simpson integration in the lensing efficiency --- src/MathUtils.jl | 26 +++++++++++++++++++++----- src/WeightFunctions.jl | 10 ++-------- 2 files changed, 23 insertions(+), 13 deletions(-) diff --git a/src/MathUtils.jl b/src/MathUtils.jl index 2de6ee63..35316a75 100644 --- a/src/MathUtils.jl +++ b/src/MathUtils.jl @@ -67,20 +67,36 @@ function SimpsonWeightArray(n::Int64) weight_array[Int((i-1)*2+3)] += 1/3 end else - weight_array[1] +=0.5 - weight_array[2] +=0.5 + weight_array[1] += 0.5 + weight_array[2] += 0.5 for i in 1:number_intervals weight_array[Int((i-1)*2+1)+1] += 1/3 weight_array[Int((i-1)*2+2)+1] += 4/3 weight_array[Int((i-1)*2+3)+1] += 1/3 end - weight_array[length(weight_array)] +=0.5 - weight_array[length(weight_array)-1] +=0.5 + weight_array[length(weight_array)] += 0.5 + weight_array[length(weight_array)-1] += 0.5 for i in 1:number_intervals weight_array[Int((i-1)*2+1)] += 1/3 weight_array[Int((i-1)*2+2)] += 4/3 weight_array[Int((i-1)*2+3)] += 1/3 end + weight_array ./= 2 end - return weight_array ./2 + return weight_array +end + +""" + SimpsonWeightMatrix(n::Int64) + +This function evaluates an array with the Simpson weight, for Simpson +Integration of the Lensing Efficiency +""" +function SimpsonWeightMatrix(n::Int64) + number_intervals = floor((n-1)/2) + weight_matrix = zeros(n, n) + for i in 1:n + weight_matrix[i,i:n] = SimpsonWeightArray(n-i+1) + end + return weight_matrix end diff --git a/src/WeightFunctions.jl b/src/WeightFunctions.jl index 0bbcb70f..4ed552db 100644 --- a/src/WeightFunctions.jl +++ b/src/WeightFunctions.jl @@ -127,13 +127,7 @@ function ComputeLensingEfficiencyOverGridCustom( BackgroundQuantities::BackgroundQuantities, w0waCDMCosmology::w0waCDMCosmologyStruct) WLWeightFunction.LensingEfficiencyArray .*= 0 - elimination_matrix = zeros(length(CosmologicalGrid.ZArray), - length(CosmologicalGrid.ZArray)) - for i in 1:length(CosmologicalGrid.ZArray) - for j in i:length(CosmologicalGrid.ZArray) - elimination_matrix[i,j] = 1 - end - end + Weight_Matrix = SimpsonWeightMatrix(length(CosmologicalGrid.ZArray)) @avx for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1 for idx_ZArray in 1:length(CosmologicalGrid.ZArray) for idx_ZArrayInt in 1:length(CosmologicalGrid.ZArray) @@ -141,7 +135,7 @@ function ComputeLensingEfficiencyOverGridCustom( idx_ZArray] += ConvolvedDensity.DensityGridArray[idx_ZBinArray, idx_ZArrayInt] * (1 - BackgroundQuantities.rZArray[idx_ZArray] / BackgroundQuantities.rZArray[idx_ZArrayInt]) * - elimination_matrix[idx_ZArray, idx_ZArrayInt] + Weight_Matrix[idx_ZArray, idx_ZArrayInt] end end end From 8ca7b0b1dbe556c6f890559ba73d7be1c3591a44 Mon Sep 17 00:00:00 2001 From: Marco Bonici Date: Mon, 8 Mar 2021 21:31:51 +0100 Subject: [PATCH 4/9] Changed tollerance angular coefficients test --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 2a20ad6c..a783e73d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -251,5 +251,5 @@ end AngularCoefficientsReloaded = CosmoCentral.ReadAngularCoefficients( "new_cl") @test isapprox(AngularCoefficientsReloaded.AngularCoefficientsArray, - AngularCoefficients.AngularCoefficientsArray, rtol=1e-9) + AngularCoefficients.AngularCoefficientsArray, rtol=1e-5) end From 73e6bbc1043be7667f075e9427494adfe3dd70cc Mon Sep 17 00:00:00 2001 From: Marco Bonici Date: Mon, 8 Mar 2021 21:42:53 +0100 Subject: [PATCH 5/9] Changed tollerance test --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index a783e73d..5f4e9d7f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -245,7 +245,7 @@ end GCWeightFunction, GCWeightFunction, BackgroundQuantities, w0waCDMCosmology, CosmologicalGrid, PowerSpectrum, CosmoCentral.CustomTrapz()) @test isapprox(AngularCoefficientsLoaded.AngularCoefficientsArray, - AngularCoefficients.AngularCoefficientsArray, rtol=1e-9) + AngularCoefficients.AngularCoefficientsArray, rtol=1e-5) CosmoCentral.WriteAngularCoefficients("PhotometricGalaxy_PhotometricGalaxy", AngularCoefficients, "new_cl") AngularCoefficientsReloaded = CosmoCentral.ReadAngularCoefficients( From 0c8c412131bcfd8c935e0a6183aba9957964ca6c Mon Sep 17 00:00:00 2001 From: Marco Bonici Date: Sat, 13 Mar 2021 19:17:30 +0100 Subject: [PATCH 6/9] Added cosmology output json --- Project.toml | 3 ++- src/AngularCoefficients.jl | 6 ++++++ src/CosmoCentral.jl | 1 + src/FileIO.jl | 16 ++++++++++++++++ src/Forecaster.jl | 1 + test/Project.toml | 1 + 6 files changed, 27 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1d9d97c0..88a418c8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "CosmoCentral" uuid = "ad8f630a-9b1f-4174-92d0-7638535e3d6b" authors = ["marcobonici "] -version = "0.0.1" +version = "0.0.2" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" @@ -12,6 +12,7 @@ DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" diff --git a/src/AngularCoefficients.jl b/src/AngularCoefficients.jl index 97f6384c..283e4acc 100644 --- a/src/AngularCoefficients.jl +++ b/src/AngularCoefficients.jl @@ -81,4 +81,10 @@ function ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients, Integrand .*= (last(CosmologicalGrid.ZArray)- first(CosmologicalGrid.ZArray))/(length(CosmologicalGrid.ZArray)-1) AngularCoefficients.AngularCoefficientsArray = Integrand + if any(isnan,Integrand) + println(w0waCDMCosmology) + throw("There is a NaN in the array you evaluated!") + else + println("Bella!") + end end diff --git a/src/CosmoCentral.jl b/src/CosmoCentral.jl index 606e81ac..21079283 100644 --- a/src/CosmoCentral.jl +++ b/src/CosmoCentral.jl @@ -11,6 +11,7 @@ using HDF5 using LoopVectorization using Statistics using JSON +using JSON3 numpy = pyimport("numpy") classy = pyimport("classy") diff --git a/src/FileIO.jl b/src/FileIO.jl index 22b03e40..ec9336e0 100644 --- a/src/FileIO.jl +++ b/src/FileIO.jl @@ -4,6 +4,22 @@ function WriteAngularCoefficients(Probes::String, AngularCoefficients.AngularCoefficientsArray) end +function WriteCosmology(Cosmology::w0waCDMCosmologyStruct, Filename::String) + CosmoDict = Dict{String,Float64}() + CosmoDict["w0"] = Cosmology.w0 + CosmoDict["wa"] = Cosmology.wa + CosmoDict["Mν"] = Cosmology.Mν + CosmoDict["H0"] = Cosmology.H0 + CosmoDict["ΩM"] = Cosmology.ΩM + CosmoDict["ΩB"] = Cosmology.ΩB + CosmoDict["ΩDE"] = Cosmology.ΩDE + CosmoDict["Ωk"] = Cosmology.Ωk + CosmoDict["Ωr"] = Cosmology.Ωr + CosmoDict["ns"] = Cosmology.ns + CosmoDict["σ8"] = Cosmology.σ8 + JSON3.write(Filename*"/cosmology.json", CosmoDict) +end + function ReadAngularCoefficients(Filename::String) Filename *= ".h5" file = HDF5.h5open(Filename, "r") diff --git a/src/Forecaster.jl b/src/Forecaster.jl index b11a9d2a..fa11b20e 100644 --- a/src/Forecaster.jl +++ b/src/Forecaster.jl @@ -299,6 +299,7 @@ function InitializeComputeAngularCoefficients(ProbesDict::Dict, CosmoCentral.CustomTrapz()) WriteAngularCoefficients(key_A*"_"*key_B, AngularCoefficients, PathOutput*key*"/cl") + WriteCosmology(w0waCDMCosmology, PathOutput*key) end end end diff --git a/test/Project.toml b/test/Project.toml index 2ec9d656..4ebef31b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,6 +5,7 @@ Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" From d6ab63d2f4087682d2aae7537701f532a2d8f152 Mon Sep 17 00:00:00 2001 From: Marco Bonici Date: Sat, 13 Mar 2021 20:06:57 +0100 Subject: [PATCH 7/9] Added safety check in Cl computation --- src/AngularCoefficients.jl | 47 ++++++++++++++++++++------------------ 1 file changed, 25 insertions(+), 22 deletions(-) diff --git a/src/AngularCoefficients.jl b/src/AngularCoefficients.jl index 283e4acc..37d4b4fd 100644 --- a/src/AngularCoefficients.jl +++ b/src/AngularCoefficients.jl @@ -63,28 +63,31 @@ function ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients, PowerSpectrum::PowerSpectrum, ::CustomTrapz) c_0 = 2.99792458e5 #TODO: find a package containing the exact value of #physical constants involved in calculations - Integrand = similar(AngularCoefficients.AngularCoefficientsArray) .*0 - Simpson_weights = SimpsonWeightArray(length(CosmologicalGrid.ZArray)) - @avx for i ∈ axes(AngularCoefficients.AngularCoefficientsArray,2), - j ∈ axes(AngularCoefficients.AngularCoefficientsArray,3), - l ∈ axes(AngularCoefficients.AngularCoefficientsArray,1) - for z ∈ axes(CosmologicalGrid.ZArray,1) - Integrand[l,i,j] += c_0 * - WeightFunctionA.WeightFunctionArray[i, z] * - WeightFunctionB.WeightFunctionArray[j, z] / - (BackgroundQuantities.HZArray[z] * - BackgroundQuantities.rZArray[z]^2) * - PowerSpectrum.InterpolatedPowerSpectrum[l,z] * - Simpson_weights[z] + check = true + while check == true + Integrand = similar(AngularCoefficients.AngularCoefficientsArray) .*0 + Simpson_weights = SimpsonWeightArray(length(CosmologicalGrid.ZArray)) + @avx for i ∈ axes(AngularCoefficients.AngularCoefficientsArray,2), + j ∈ axes(AngularCoefficients.AngularCoefficientsArray,3), + l ∈ axes(AngularCoefficients.AngularCoefficientsArray,1) + for z ∈ axes(CosmologicalGrid.ZArray,1) + Integrand[l,i,j] += c_0 * + WeightFunctionA.WeightFunctionArray[i, z] * + WeightFunctionB.WeightFunctionArray[j, z] / + (BackgroundQuantities.HZArray[z] * + BackgroundQuantities.rZArray[z]^2) * + PowerSpectrum.InterpolatedPowerSpectrum[l,z] * + Simpson_weights[z] + end + end + Integrand .*= (last(CosmologicalGrid.ZArray)- + first(CosmologicalGrid.ZArray))/(length(CosmologicalGrid.ZArray)-1) + AngularCoefficients.AngularCoefficientsArray = Integrand + if any(isnan,Integrand) + println("There is a problem, we need to evaluate the coefficients + for this cosmology again!!",w0waCDMCosmology) + else + check = false end - end - Integrand .*= (last(CosmologicalGrid.ZArray)- - first(CosmologicalGrid.ZArray))/(length(CosmologicalGrid.ZArray)-1) - AngularCoefficients.AngularCoefficientsArray = Integrand - if any(isnan,Integrand) - println(w0waCDMCosmology) - throw("There is a NaN in the array you evaluated!") - else - println("Bella!") end end From 2cd543cfd1e681fd7fa89e3a0aa1a8df9b487f34 Mon Sep 17 00:00:00 2001 From: Marco Bonici Date: Sun, 14 Mar 2021 09:56:44 +0100 Subject: [PATCH 8/9] Modified runtests tolerance --- test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 5f4e9d7f..2a20ad6c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -245,11 +245,11 @@ end GCWeightFunction, GCWeightFunction, BackgroundQuantities, w0waCDMCosmology, CosmologicalGrid, PowerSpectrum, CosmoCentral.CustomTrapz()) @test isapprox(AngularCoefficientsLoaded.AngularCoefficientsArray, - AngularCoefficients.AngularCoefficientsArray, rtol=1e-5) + AngularCoefficients.AngularCoefficientsArray, rtol=1e-9) CosmoCentral.WriteAngularCoefficients("PhotometricGalaxy_PhotometricGalaxy", AngularCoefficients, "new_cl") AngularCoefficientsReloaded = CosmoCentral.ReadAngularCoefficients( "new_cl") @test isapprox(AngularCoefficientsReloaded.AngularCoefficientsArray, - AngularCoefficients.AngularCoefficientsArray, rtol=1e-5) + AngularCoefficients.AngularCoefficientsArray, rtol=1e-9) end From 1fea7b7132c53344550a5a3083a62d433f384d0e Mon Sep 17 00:00:00 2001 From: Marco Bonici Date: Sun, 14 Mar 2021 10:05:56 +0100 Subject: [PATCH 9/9] Fixing test tolerance --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 2a20ad6c..4587e553 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -245,7 +245,7 @@ end GCWeightFunction, GCWeightFunction, BackgroundQuantities, w0waCDMCosmology, CosmologicalGrid, PowerSpectrum, CosmoCentral.CustomTrapz()) @test isapprox(AngularCoefficientsLoaded.AngularCoefficientsArray, - AngularCoefficients.AngularCoefficientsArray, rtol=1e-9) + AngularCoefficients.AngularCoefficientsArray, rtol=1e-6) CosmoCentral.WriteAngularCoefficients("PhotometricGalaxy_PhotometricGalaxy", AngularCoefficients, "new_cl") AngularCoefficientsReloaded = CosmoCentral.ReadAngularCoefficients(