Merge pull request #48 from marcobonici/develop

Develop
marcobonici · Mar 14, 2021 · c642eac · c642eac
2 parents c9166db + 1fea7b7
commit c642eac
Show file tree

Hide file tree

Showing 10 changed files with 103 additions and 21 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "CosmoCentral"
 uuid = "ad8f630a-9b1f-4174-92d0-7638535e3d6b"
 authors = ["marcobonici <[email protected]>"]
-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/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/AngularCoefficients.jl b/src/AngularCoefficients.jl
@@ -63,20 +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
-    @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]
+    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
 end
diff --git 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
@@ -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
@@ -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/src/MathUtils.jl b/src/MathUtils.jl
@@ -50,3 +50,53 @@ 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
+        weight_array ./= 2
+    end
+    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
@@ -127,14 +127,15 @@ 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
+    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 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]) *
+                Weight_Matrix[idx_ZArray, idx_ZArrayInt]
             end
         end
     end

diff --git 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"

diff --git 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(