Merge pull request #31 from marcobonici/develop

Develop
marcobonici · Jan 18, 2021 · 98b74c0 · 98b74c0
2 parents 1445422 + afdee7e
commit 98b74c0
Show file tree

Hide file tree

Showing 11 changed files with 381 additions and 35 deletions.
diff --git a/Project.toml b/Project.toml
@@ -8,7 +8,9 @@ Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PlotThemes = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
@@ -26,8 +28,8 @@ DocumenterTools = "0.1"
 Interpolations = "0.13"
 NumericalIntegration = "0.2"
 Parameters = "0.12"
-PlotlyBase = "0.4"
 PlotThemes = "2.0"
+PlotlyBase = "0.4"
 PyCall = "1.92"
 QuadGK = "2.4"
 SciPy = "0.1"
diff --git a/src/AngularCoefficients.jl b/src/AngularCoefficients.jl
@@ -1,18 +1,24 @@
+abstract type IntegrationMethod end
+
+struct NumericalIntegrationSimpson <: IntegrationMethod end
+struct CustomTrapz <: IntegrationMethod end
+
+
 """
     ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients,
     WeightFunctionA::GCWeightFunction, WeightFunctionB::GCWeightFunction,
     BackgroundQuantities::BackgroundQuantities,
     w0waCDMCosmology::w0waCDMCosmology, CosmologicalGrid::CosmologicalGrid,
-    PowerSpectrum::PowerSpectrum)
+    PowerSpectrum::PowerSpectrum, ::NumericalIntegrationSimpson)
 
 This function evaluates the Angular Coefficients for all tomographic bins and
 multipole values.
 """
 function  ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients,
-    WeightFunctionA::GCWeightFunction, WeightFunctionB::GCWeightFunction,
+    WeightFunctionA::WeightFunction, WeightFunctionB::WeightFunction,
     BackgroundQuantities::BackgroundQuantities,
     w0waCDMCosmology::w0waCDMCosmology, CosmologicalGrid::CosmologicalGrid,
-    PowerSpectrum::PowerSpectrum)
+    PowerSpectrum::PowerSpectrum, ::NumericalIntegrationSimpson)
     c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                        #physical constants involved in calculations
     for idx_a in 1:size(AngularCoefficients.AngularCoefficientsArray, 2)
@@ -27,7 +33,36 @@ function  ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients,
                 AngularCoefficients.AngularCoefficientsArray[idx_l, idx_a,
                 idx_b] = NumericalIntegration.integrate(CosmologicalGrid.ZArray,
                 Integrand, SimpsonEvenFast())
+                AngularCoefficients.AngularCoefficientsArray[idx_l, idx_b,
+                idx_a] = AngularCoefficients.AngularCoefficientsArray[idx_l,
+                idx_a, idx_b]
             end
         end
     end
 end
+
+function  ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients,
+    WeightFunctionA::WeightFunction, WeightFunctionB::WeightFunction,
+    BackgroundQuantities::BackgroundQuantities,
+    w0waCDMCosmology::w0waCDMCosmology, CosmologicalGrid::CosmologicalGrid,
+    PowerSpectrum::PowerSpectrum, ::CustomTrapz)
+    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
+                       #physical constants involved in calculations
+    Integrand = zeros(size(AngularCoefficients.AngularCoefficientsArray,1),
+    size(AngularCoefficients.AngularCoefficientsArray,2),
+    size(AngularCoefficients.AngularCoefficientsArray,3))
+    @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]
+        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
@@ -4,8 +4,11 @@ using QuadGK
 using Base: @kwdef
 using Parameters
 using NumericalIntegration
+using Interpolations
 using PyCall
 using Dierckx
+using HDF5
+using LoopVectorization
 numpy = pyimport("numpy")
 classy = pyimport("classy")
 
@@ -18,5 +21,6 @@ include("WeightFunctions.jl")
 include("BoltzmannSolver.jl")
 include("PowerSpectrum.jl")
 include("AngularCoefficients.jl")
+include("FileIO.jl")
 
 end # module
diff --git a/src/CosmologicalStructures.jl b/src/CosmologicalStructures.jl
@@ -13,6 +13,7 @@ abstract type ConvolvedDensity <: AsbtractDensity end
 abstract type InstrumentResponse end
 abstract type WeightFunction end
 abstract type GCWeightFunction <: WeightFunction end
+abstract type WLWeightFunction <: WeightFunction end
 abstract type PowerSpectrum end
 abstract type AngularCoefficients end
 
@@ -143,7 +144,7 @@ The parameters contained in this struct are
 @kwdef mutable struct AnalitycalDensityStruct <: AnalitycalDensity
     Z0::Float64 = 0.9/sqrt(2.)
     ZMin::Float64 = 0.001
-    ZMax::Float64 = 2.5
+    ZMax::Float64 = 4.0
     SurfaceDensity::Float64 = 30.
     Normalization::Float64 = 1.
 end
@@ -211,6 +212,18 @@ for all tomographic bins and redshift values in the [`CosmologicalGridStruct`](@
     WeightFunctionArray::AbstractArray{Float64, 2}
 end
 
+"""
+    WLWeightFunctionStruct()
+
+This struct contains the array with the Weak Lensing Weight function values
+for all tomographic bins and redshift values in the [`CosmologicalGridStruct`](@ref).
+"""
+@kwdef mutable struct WLWeightFunctionStruct <: WLWeightFunction
+    WeightFunctionArray::AbstractArray{Float64, 2}
+    LensingEfficiencyArray::AbstractArray{Float64, 2}
+end
+
+
 """
     PowerSpectrumStruct()
 
@@ -235,3 +248,9 @@ This struct contains the array with the Angular Coefficients.
 @kwdef mutable struct AngularCoefficientsStruct <: AngularCoefficients
     AngularCoefficientsArray::AbstractArray{Float64, 3} = zeros(2991, 10, 10)
 end
+
+abstract type InterpolationMethod end
+
+struct RectBivSplineDierckx <: InterpolationMethod end
+struct GriddedLinear <: InterpolationMethod end
+struct BSplineCubic <: InterpolationMethod end
diff --git a/src/FileIO.jl b/src/FileIO.jl
@@ -0,0 +1,90 @@
+function WriteAngularCoefficients(AngularCoefficients::AngularCoefficients,
+    CosmologicalGrid::CosmologicalGrid, WeightFunction::WeightFunction,
+    Bias::Bias, ConvolvedDensity::ConvolvedDensity, Filename::String)
+    h5write(Filename*".h5", "cls/PhotometricGalaxy_PhotometricGalaxy/c_lij",
+    AngularCoefficients.AngularCoefficientsArray)
+    h5write(Filename*".h5",
+    "weight_functions/PhotometricGalaxy/bias/b_i_z",
+    Bias.BiasArray)
+    h5write(Filename*".h5",
+    "weight_functions/PhotometricGalaxy/density/norm_density_iz",
+    ConvolvedDensity.DensityNormalizationArray)
+end
+
+function ReadAngularCoefficients(Filename::String)
+    Filename *= ".h5"
+    file = HDF5.h5open(Filename, "r")
+    c_lij =
+    HDF5.read(file["cls"]["PhotometricGalaxy_PhotometricGalaxy"]["c_lij"])
+    AngularCoefficients = AngularCoefficientsStruct(AngularCoefficientsArray =
+    c_lij)
+    return AngularCoefficients
+end
+
+function WritePowerSpectrumBackground(PowerSpectrum::PowerSpectrum,
+    BackgroundQuantities::BackgroundQuantities,
+    CosmologicalGrid::CosmologicalGrid, Filename::String)
+    h5write(Filename*".h5",
+    "cosmology/comoving_distance_array",
+    BackgroundQuantities.rZArray)
+    h5write(Filename*".h5",
+    "cosmology/hubble_array",
+    BackgroundQuantities.HZArray)
+    h5write(Filename*".h5",
+    "cosmology/z_grid",
+    CosmologicalGrid.ZArray)
+    h5write(Filename*".h5",
+    "power_spectrum/z_grid",
+    CosmologicalGrid.ZArray)
+    h5write(Filename*".h5",
+    "power_spectrum/k_grid",
+    CosmologicalGrid.KArray)
+    h5write(Filename*".h5",
+    "power_spectrum/lin_p_mm_k_z",
+    PowerSpectrum.PowerSpectrumLinArray)
+    h5write(Filename*".h5",
+    "power_spectrum/nonlin_p_mm_k_z",
+    PowerSpectrum.PowerSpectrumNonlinArray)
+end
+
+function ReadPowerSpectrumBackground(Filename::String,
+    MultipolesArray::Vector{Float64})
+    Filename *= ".h5"
+    file = HDF5.h5open(Filename, "r")
+    nonlin_p_mm_k_z = HDF5.read(file["power_spectrum"]["nonlin_p_mm_k_z"])
+    lin_p_mm_k_z    = HDF5.read(file["power_spectrum"]["lin_p_mm_k_z"])
+    z_grid = HDF5.read(file["power_spectrum"]["z_grid"])
+    k_grid = HDF5.read(file["power_spectrum"]["k_grid"])
+    comoving_distance_array = HDF5.read(file["cosmology"]["comoving_distance_array"])
+    hubble_array = HDF5.read(file["cosmology"]["hubble_array"])
+    BackgroundQuantities = BackgroundQuantitiesStruct(HZArray = hubble_array,
+    rZArray = comoving_distance_array)
+    CosmologicalGrid = CosmologicalGridStruct(ZArray = z_grid, KArray = k_grid,
+    MultipolesArray = MultipolesArray)
+    PowerSpectrum = PowerSpectrumStruct(PowerSpectrumLinArray = lin_p_mm_k_z,
+    PowerSpectrumNonlinArray = nonlin_p_mm_k_z,
+    InterpolatedPowerSpectrum = zeros(length(CosmologicalGrid.MultipolesArray), length(CosmologicalGrid.ZArray)))
+    return PowerSpectrum, BackgroundQuantities, CosmologicalGrid
+end
+
+function ReadPowerSpectrumBackgroundSeyfert(Filename::String,
+    MultipolesArray::Vector{Float64})
+    c_0 = 2.99792458e5
+    Filename *= ".h5"
+    file = HDF5.h5open(Filename, "r")
+    nonlin_p_mm_k_z = HDF5.read(file["power_spectrum"]["nonlin_p_mm_z_k"])
+    lin_p_mm_k_z    = HDF5.read(file["power_spectrum"]["lin_p_mm_z_k"])
+    z_grid = HDF5.read(file["power_spectrum"]["z_grid"])
+    k_grid = HDF5.read(file["power_spectrum"]["k_grid"])
+    dimensionless_comoving_distance_array = HDF5.read(file["cosmology"]["dimensionless_comoving_distance_array"])
+    dimensionless_hubble_array = HDF5.read(file["cosmology"]["dimensionless_hubble_array"])
+    BackgroundQuantities = BackgroundQuantitiesStruct(HZArray =
+    dimensionless_hubble_array*67,
+    rZArray = dimensionless_comoving_distance_array*c_0/67)
+    CosmologicalGrid = CosmologicalGridStruct(ZArray = z_grid, KArray = k_grid,
+    MultipolesArray = MultipolesArray)
+    PowerSpectrum = PowerSpectrumStruct(PowerSpectrumLinArray = lin_p_mm_k_z,
+    PowerSpectrumNonlinArray = nonlin_p_mm_k_z,
+    InterpolatedPowerSpectrum = zeros(length(CosmologicalGrid.MultipolesArray), length(CosmologicalGrid.ZArray)))
+    return PowerSpectrum, BackgroundQuantities, CosmologicalGrid
+end
diff --git a/src/PowerSpectrum.jl b/src/PowerSpectrum.jl
@@ -24,7 +24,8 @@ This function interpolates the Power Spectrum on the ``k-z`` grid and evaluates
 it on the Limber grid.
 """
 function InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid::CosmologicalGrid,
-    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum)
+    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
+    ::RectBivSplineDierckx)
     InterpPmm = Dierckx.Spline2D(log10.(CosmologicalGrid.KArray),
     CosmologicalGrid.ZArray, log10.(PowerSpectrum.PowerSpectrumNonlinArray);
     kx=5, ky=5, s=0.0)
@@ -34,3 +35,36 @@ function InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid::CosmologicalGrid,
         CosmologicalGrid.ZArray))
     end
 end
+
+function InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
+    ::GriddedLinear)
+    InterpPmm = Interpolations.interpolate((log10.(CosmologicalGrid.KArray),
+    CosmologicalGrid.ZArray), log10.(PowerSpectrum.PowerSpectrumNonlinArray),
+    Gridded(Linear()))
+    InterpPmm = Interpolations.extrapolate(InterpPmm, Line())
+    for (idx_l, myl) in enumerate(CosmologicalGrid.MultipolesArray)
+        PowerSpectrum.InterpolatedPowerSpectrum[idx_l, :] =
+        10 .^(InterpPmm.(log10.(CosmologicalGrid.KLimberArray[idx_l, :]),
+        CosmologicalGrid.ZArray))
+    end
+end
+
+function InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
+    ::BSplineCubic)
+    x = LinRange(log10(first(CosmologicalGrid.KArray)),
+    log10(last(CosmologicalGrid.KArray)), length(CosmologicalGrid.KArray))
+    y = LinRange(first(CosmologicalGrid.ZArray), last(CosmologicalGrid.ZArray),
+    length(CosmologicalGrid.ZArray))
+    InterpPmm = Interpolations.interpolate(
+    log10.(PowerSpectrum.PowerSpectrumNonlinArray),
+    BSpline(Cubic(Line(OnGrid()))))
+    InterpPmm = scale(InterpPmm, x, y)
+    InterpPmm = Interpolations.extrapolate(InterpPmm, Line())
+    for (idx_l, myl) in enumerate(CosmologicalGrid.MultipolesArray)
+        PowerSpectrum.InterpolatedPowerSpectrum[idx_l, :] =
+        10 .^(InterpPmm.(log10.(CosmologicalGrid.KLimberArray[idx_l, :]),
+        CosmologicalGrid.ZArray))
+    end
+end
diff --git a/src/WeightFunctions.jl b/src/WeightFunctions.jl
@@ -47,3 +47,83 @@ function ComputeWeightFunctionOverGrid(
         end
     end
 end
+
+"""
+    ComputeWeightFunction(z::Float64, i::Int64,
+        ConvolvedDensity::ConvolvedDensity,
+        w0waCDMCosmology::w0waCDMCosmology)
+
+This function returns the Galaxy Clustering Weight function for a given redshift
+``z`` and tomographic bin ``i``.
+"""
+function ComputeLensingEfficiency(z::Float64, i::Int64,
+    ConvolvedDensity::ConvolvedDensity, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, w0waCDMCosmology::w0waCDMCosmology,
+    CosmologicalGrid::CosmologicalGrid, WLWeightFunction::WLWeightFunction)
+    int, err = QuadGK.quadgk(x -> CosmoCentral.ComputeConvolvedDensityFunction(x, i,
+    ConvolvedDensity, AnalitycalDensity, InstrumentResponse)*
+    (1. - CosmoCentral.ComputeComovingDistance(z, w0waCDMCosmology)/
+    CosmoCentral.ComputeComovingDistance(x, w0waCDMCosmology)), z,
+    last(CosmologicalGrid.ZArray) , rtol=1e-12)
+    return int
+end
+
+"""
+    ComputeWeightFunctionOverGrid(
+    GCWeightFunction::GCWeightFunction, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, ConvolvedDensity::ConvolvedDensity,
+    Bias::Bias, CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities,
+    w0waCDMCosmology::w0waCDMCosmology)
+
+This function evaluates the Galaxy Clustering Weight Function over the ``z``
+grid and for all tomographic bins ``i``.
+"""
+function ComputeLensingEfficiencyOverGrid(
+    WLWeightFunction::WLWeightFunction, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, ConvolvedDensity::ConvolvedDensity,
+    CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities,
+    w0waCDMCosmology::w0waCDMCosmology)
+    for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
+        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
+            WLWeightFunction.LensingEfficiencyArray[idx_ZBinArray, idx_ZArray] =
+            ComputeLensingEfficiency(CosmologicalGrid.ZArray[idx_ZArray],
+            idx_ZBinArray, ConvolvedDensity, AnalitycalDensity,
+            InstrumentResponse, w0waCDMCosmology, CosmologicalGrid,
+            WLWeightFunction)
+        end
+    end
+end
+
+function ComputeWeightFunction(z::Float64, i::Int64,
+    ConvolvedDensity::ConvolvedDensity, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, w0waCDMCosmology::w0waCDMCosmology,
+    CosmologicalGrid::CosmologicalGrid, WLWeightFunction::WLWeightFunction)
+    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
+                       #physical constants involved in calculations
+    return 1.5 * ComputeLensingEfficiency(z, i, ConvolvedDensity,
+    AnalitycalDensity,InstrumentResponse, w0waCDMCosmology, CosmologicalGrid,
+    WLWeightFunction) * (w0waCDMCosmology.H0/c_0)^2 * w0waCDMCosmology.ΩM *
+    (1. + z) * ComputeComovingDistance(z, w0waCDMCosmology)
+end
+
+
+function ComputeWeightFunctionOverGrid(
+    WLWeightFunction::WLWeightFunction, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, ConvolvedDensity::ConvolvedDensity,
+    CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities,
+    w0waCDMCosmology::w0waCDMCosmology)
+    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
+                       #physical constants involved in calculations
+    for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
+        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
+            WLWeightFunction.WeightFunctionArray[idx_ZBinArray, idx_ZArray] =
+            1.5 * (w0waCDMCosmology.H0/c_0)^2 * w0waCDMCosmology.ΩM *
+            (1. + CosmologicalGrid.ZArray[idx_ZArray]) *
+            BackgroundQuantities.rZArray[idx_ZArray] *
+            WLWeightFunction.LensingEfficiencyArray[idx_ZBinArray, idx_ZArray]
+        end
+    end
+end
diff --git a/test/Project.toml b/test/Project.toml
@@ -1,6 +1,9 @@
 [deps]
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"

diff --git a/test/cl.h5 b/test/cl.h5
diff --git a/test/p_mm.h5 b/test/p_mm.h5