Add z dep

flip/covariance/carreres23/coefficients.py

@@ -5,6 +5,8 @@ def get_coefficients(
     model_type,
     parameter_values_dict,
     variant=None,
+    redshift_dict=None,
+    power_spectrum_amplitude_function=None,
 ):
     coefficients_dict = {}
     coefficients_dict["vv"] = [parameter_values_dict["fs8"] ** 2]

flip/covariance/generator.py

@@ -12,6 +12,7 @@
 from flip.covariance.adamsblake20 import flip_terms as flip_terms_adamsblake20
 from flip.covariance.carreres23 import flip_terms as flip_terms_carreres23
 from flip.covariance.lai22 import flip_terms as flip_terms_lai22
+from flip.covariance.rcrk24 import flip_terms as flip_terms_rcrk24
 from flip.covariance.ravouxcarreres import flip_terms as flip_terms_ravouxcarreres
 from flip.utils import create_log
 
@@ -22,6 +23,7 @@
     "lai22",
     "carreres23",
     "ravouxcarreres",
+    "rcrk24",
 ]
 
 

flip/covariance/rcrk24/coefficients.py

@@ -14,7 +14,7 @@ def get_coefficients(
     cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"])
 
     coefficient_vector = (
-        cosmo.Om(redshift_velocities).value ** parameter_values_dict["gamma"]
+        np.array(cosmo.Om(redshift_velocities)) ** parameter_values_dict["gamma"]
         * cosmo.H(redshift_velocities).value
         / cosmo.H0
         * power_spectrum_amplitude_function(redshift_velocities)

flip/covariance/rcrk24/fisher_terms.py

@@ -1,6 +1,23 @@
 import numpy as np
 from astropy.cosmology import FlatLambdaCDM
 
+# The flip convention is to split the power spectrum into several terms
+# where linearity assumptions are made
+# P_ab = A * B * P0_xy
+#
+# A is the coefficients where
+# P_ab = A P_xy
+# B is the cofficient where
+# P_xy = B P0_xy
+# and
+# P0_xy is the power spectrum for a fiducial cosmology at z=0
+
+# The derivative cofficients we need are
+# dA/dp B + A dB/dp
+
+# for vv
+# A = (aHfs8)_1 * (aHfs8_1)
+# B = s8_1 * s8_2
 
 def get_partial_derivative_coefficients(
     model_type,
@@ -10,35 +27,112 @@ def get_partial_derivative_coefficients(
     power_spectrum_amplitude_function=None,
 ):
 
-    cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"])
+    # s80 is considered to be fixed by the CMB and is hence not a fit parameter
+
+    s80 = 0.832
     redshift_velocities = redshift_dict["v"]
+    a = 1/(1+redshift_velocities)
+
+    cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"])
+
+
+    # The Om0-gamma model f=Omega(Om0)^gamma
+
+    # pre-calculate a few things that are used repeatedly
+    cosmoOm=np.array(cosmo.Om(redshift_velocities))
+    f = cosmoOm**parameter_values_dict["gamma"]
+    f0 = parameter_values_dict["Om0"]**parameter_values_dict["gamma"]
+
+
+
+    # Calculation of s8 and its derivatives requires an integral.  It is useful to
+    # expand Omega in terms of (1-a), which allows analytic solutions
+
+    # approximation
+    def lnD(a):
+        return np.log(a)*(f0+f0*3*parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])) \
+            + (1-a)*f0*3*parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])
+
+    def dlnDdOm0(a):
+        return parameter_values_dict["gamma"] *parameter_values_dict["Om0"]**(parameter_values_dict["gamma"]-1) \
+            * ( 3 * (a-1)*(parameter_values_dict["gamma"] * (parameter_values_dict["Om0"]-1) + parameter_values_dict["Om0"]) \
+                + np.log(a)*(-3*parameter_values_dict["gamma"] * (parameter_values_dict["Om0"]-1) - 3* parameter_values_dict["Om0"] +1 ))
+
+    def dlnDdgamma(a):
+        return 3*(1-a)*(1-parameter_values_dict["Om0"])*f0 \
+            + 3*(1-a)*parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])*f0 *np.log(parameter_values_dict["Om0"]) \
+            + np.log(a) \
+                * ( \
+                    3*(1-parameter_values_dict["Om0"])*f0 \
+                    + 3 * parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])*f0 *np.log(parameter_values_dict["Om0"]) \
+                    + f0*np.log(parameter_values_dict["Om0"]) \
+                )
+
+    def s8(a):
+        return s80 * np.exp(lnD(a))
+
+    def aHfs8(a):
+        return f * s8(a) / (1+redshift_velocities) * cosmo.H(redshift_velocities)/cosmo.H0 
+
+    # now for the partials
+
+    def dOmdOm0(a):
+        numerator = parameter_values_dict["Om0"] * a**(-3) 
+        denominator = (numerator + 1 - parameter_values_dict["Om0"])
+        return a**(-3)/denominator  - numerator/denominator**2*(a**(-3)-1)
+
+    def dfdOm0(a):
+        return parameter_values_dict["gamma"] * f / cosmoOm * dOmdOm0(a)
+
+    def dfdgamma(a):
+        return np.log(cosmoOm)  * f
+
+    def ds8dOm0(a):
+        return s8(a)*dlnDdOm0(a)
+
+    def ds8dgamma(a):
+        return s8(a)*dlnDdgamma(a)
+
+    def dAdOm0(a):
+        return a*cosmo.H(a)/cosmo.H0*(dfdOm0(a)*s8(a) + f*ds8dOm0(a))
+
+    def dAdgamma(a):
+        return a*cosmo.H(a)/cosmo.H0*(dfdgamma(a)*s8(a) + f*ds8dgamma(a))
+
+
+    aHfs8s8 = aHfs8(a)*s8(a)
 
     Omega_m_partial_derivative_coefficients = (
-        2
-        * cosmo.Om(redshift_velocities).value ** (2 * parameter_values_dict["gamma"])
-        * parameter_values_dict["gamma"]
-        * power_spectrum_amplitude_function(redshift_velocities) ** 2
-        / cosmo.Om(redshift_velocities).value
+        dAdOm0(a) * s8(a) + aHfs8(a) * ds8dOm0(a)
     )
 
     gamma_partial_derivative_coefficients = (
-        2
-        * cosmo.Om(redshift_velocities).value ** (2 * parameter_values_dict["gamma"])
-        * power_spectrum_amplitude_function(redshift_velocities) ** 2
-        * np.log(cosmo.Om(redshift_velocities).value)
+        dAdgamma(a) * s8(a) + aHfs8(a) * ds8dgamma(a)
     )
 
-    s8_partial_derivative_coefficients = (
-        2
-        * cosmo.Om(redshift_velocities).value ** (2 * parameter_values_dict["gamma"])
-        * power_spectrum_amplitude_function(redshift_velocities)
+    # in the fs8 case
+    def s8_fs8(a):
+        return  s80 + parameter_values_dict["fs8"] * np.log(a)
+
+    def ds8dfs8(a):
+        return np.log(a)
+
+    fs8_partial_derivative_coefficients = (
+        a*cosmo.H(redshift_velocities)/cosmo.H0 *
+        (s8_fs8(a) + parameter_values_dict["fs8"]*ds8dfs8(a))
     )
 
+    aHfs8s8_fs8 = a*cosmo.H(redshift_velocities)/cosmo.H0*parameter_values_dict["fs8"]*s8_fs8(a)
+
     partial_coefficients_dict = {
         "Omegam": {
             "vv": [
                 np.outer(
                     Omega_m_partial_derivative_coefficients,
+                    aHfs8s8,
+                ) +
+                np.outer(
+                    aHfs8s8,
                     Omega_m_partial_derivative_coefficients,
                 ),
             ],
@@ -47,15 +141,23 @@ def get_partial_derivative_coefficients(
             "vv": [
                 np.outer(
                     gamma_partial_derivative_coefficients,
+                    aHfs8s8,
+                ) +
+                np.outer(
+                    aHfs8s8,
                     gamma_partial_derivative_coefficients,
                 ),
             ],
         },
-        "s8": {
+        "fs8": {
             "vv": [
                 np.outer(
-                    s8_partial_derivative_coefficients,
-                    s8_partial_derivative_coefficients,
+                    fs8_partial_derivative_coefficients,
+                    aHfs8s8_fs8,
+                ) +
+                np.outer(
+                    aHfs8s8_fs8,
+                    fs8_partial_derivative_coefficients,
                 ),
             ],
         },

flip/covariance/rcrk24/flip_terms.py

@@ -23,6 +23,8 @@ def func(k):
 def N_vv_0_2_0(theta, phi):
     return (9 / 2) * np.cos(2 * phi) + (3 / 2) * np.cos(theta)
 
+def power_spectrum_amplitude_function(r):
+    return 1.
 
 dictionary_terms = {"vv": ["0"]}
 dictionary_lmax = {"vv": [2]}

scripts/develop.py

@@ -0,0 +1,132 @@
+import os
+
+import numpy as np
+import pandas as pd
+from pkg_resources import resource_filename
+
+from flip import fisher, utils
+from flip.covariance import covariance
+
+def main():
+    flip_base = resource_filename("flip", ".")
+    data_path = os.path.join(flip_base, "data")
+
+    ### Load data
+    sn_data = pd.read_parquet(os.path.join(data_path, "velocity_data.parquet"))
+
+    sn_data = sn_data[np.array(sn_data["status"]) != False]
+    sn_data = sn_data[np.array(sn_data["status"]) != None]
+    coordinates_velocity = np.array([sn_data["ra"], sn_data["dec"], sn_data["rcom_zobs"], sn_data["zobs"]])
+
+    data_velocity = sn_data.to_dict("list")
+    for key in data_velocity.keys():
+        data_velocity[key] = np.array(data_velocity[key])
+    data_velocity["velocity"] = data_velocity.pop("vpec")
+    data_velocity["velocity_error"] = np.zeros_like(data_velocity["velocity"])
+
+
+    ktt, ptt = np.loadtxt(os.path.join(data_path, "power_spectrum_tt.txt"))
+    kmt, pmt = np.loadtxt(os.path.join(data_path, "power_spectrum_mt.txt"))
+    kmm, pmm = np.loadtxt(os.path.join(data_path, "power_spectrum_mm.txt"))
+
+    sigmau_fiducial = 15
+
+    power_spectrum_dict = {"vv": [[ktt, ptt * utils.Du(ktt, sigmau_fiducial) ** 2]]}
+
+    ### Compute covariance
+    size_batch = 10_000
+    number_worker = 16
+
+
+    from flip.covariance.rcrk24.flip_terms import power_spectrum_amplitude_function
+    covariance_fit = covariance.CovMatrix.init_from_flip(
+        "rcrk24",
+        # "agk24"
+        # 'carreres23',
+        "velocity",
+        power_spectrum_dict,
+        coordinates_velocity=coordinates_velocity,
+        size_batch=size_batch,
+        number_worker=number_worker,
+        power_spectrum_amplitude_function=power_spectrum_amplitude_function,
+    )
+
+    ### Load fitter
+
+    fisher_properties = {
+        "inversion_method": "inverse",
+        "velocity_type": "scatter",
+    }
+
+    variant = None  # can be replaced by growth_index
+
+    parameter_dict = {
+        "fs8": 0.4,
+        "Om0": 0.3,
+        "gamma": 0.55,        
+        "sigv": 200,
+        "sigma_M": 0.12,
+    }
+
+    # parameter_dict = {
+    #     "Om0": 0.3,
+    #     "gamma": 0.55,
+    #     "sigv": 200,        
+    #     "sigma_M": 0.12,
+    # }
+
+    Fisher = fisher.FisherMatrix.init_from_covariance(
+        covariance_fit,
+        data_velocity,
+        parameter_dict,
+        fisher_properties=fisher_properties,
+    )
+
+    parameter_name_list, fisher_matrix = Fisher.compute_fisher_matrix(
+        parameter_dict, variant=variant
+    )
+    return parameter_name_list, fisher_matrix
+
+def dlnDdOm0(parameter_values_dict):
+    a=(1/(1+0.05))
+    lna=np.log(a)
+    return (
+        parameter_values_dict["gamma"] * parameter_values_dict["Om0"]**(parameter_values_dict["gamma"]-1) *
+            (
+                3 * parameter_values_dict["gamma"] * (parameter_values_dict["Om0"]-1) * (a - lna -1) +
+                3 * (a-1) * parameter_values_dict["Om0"] -
+                3 * np.log(a) * parameter_values_dict["Om0"] +  lna
+            )
+        )
+
+def dlnDdgamma(parameter_values_dict):
+    a=(1/(1+0.05))
+    lna=np.log(a)
+    f0 = parameter_values_dict["Om0"]**parameter_values_dict["gamma"]
+    return (
+            f0 *
+                (
+                    np.log(parameter_values_dict["Om0"]) *
+                        (
+                            3 * parameter_values_dict["gamma"] * (parameter_values_dict["Om0"]-1) * (a - lna -1) + lna
+                        ) +
+                    3 * (parameter_values_dict["Om0"]-1) * (a - lna -1)
+                )
+        )
+
+if __name__ == "__main__":
+    parameter_dict = {
+        "Om0": 0.3,
+        "gamma": 0.55,
+        "sigv": 200,        
+        "sigma_M": 0.12,
+    }
+    parameter_name_list, fisher_matrix = main()
+    cov = np.linalg.inv(fisher_matrix[0:2,0:2])
+
+    s80 = 0.832
+    partials = s80*np.array([parameter_dict['gamma']*parameter_dict['Om0']**(parameter_dict['gamma']-1),np.log(parameter_dict['Om0'])*parameter_dict['Om0']**parameter_dict['gamma']])
+    partials = partials + parameter_dict['Om0']**parameter_dict['gamma'] *s80 * np.array([dlnDdOm0(parameter_dict), dlnDdgamma(parameter_dict)])
+
+    print(np.sqrt(partials.T @ cov[0:2,0:2] @ partials))
+    print(1/np.sqrt(fisher_matrix[2,2]))