bayesian inference notebook ready

sambit-giri · Aug 26, 2024 · 6060e8e · 6060e8e
1 parent 6a477d0
commit 6060e8e
Show file tree

Hide file tree

Showing 10 changed files with 528 additions and 554 deletions.
diff --git a/.gitmodules b/.gitmodules
diff --git a/docs/examples/tutorials.ipynb b/docs/examples/tutorials.ipynb
diff --git a/notebooks/.ipynb_checkpoints/bayesian_inference-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/bayesian_inference-checkpoint.ipynb
diff --git a/notebooks/.ipynb_checkpoints/rotation_curves-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/rotation_curves-checkpoint.ipynb
diff --git a/notebooks/bayesian_inference.ipynb b/notebooks/bayesian_inference.ipynb
diff --git a/notebooks/hubble_diagram.ipynb b/notebooks/hubble_diagram.ipynb
diff --git a/notebooks/rotation_curves.ipynb b/notebooks/rotation_curves.ipynb
diff --git a/src/AstronomyCalc/cosmo_calculator.py b/src/AstronomyCalc/cosmo_calculator.py
@@ -20,4 +20,161 @@ def age_estimator(param, z):
     a = z_to_a(z)
     I = lambda a: 1/a/Feq.H(a=a)
     t = lambda a: quad(I, 0, a)[0]*const.Mpc_to_km/const.Gyr_to_s
-    return np.vectorize(t)(a) 
+    return np.vectorize(t)(a) 
+
+def Hubble(param, z=None, a=None):
+    """
+    Calculate the Hubble parameter at a given redshift or scale factor.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+        z (float, optional): The redshift at which to evaluate the Hubble parameter.
+        a (float, optional): The scale factor at which to evaluate the Hubble parameter.
+                             Only one of `z` or `a` should be provided.
+
+    Returns:
+        float: The Hubble parameter (H) at the specified redshift or scale factor.
+    """
+    cosmo = CosmoDistances(param)
+    return cosmo.H(z=z, a=a)
+
+def cosmic_age(param, z=None, a=None):
+    """
+    Calculate the cosmic age at a given redshift or scale factor.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+        z (float, optional): The redshift at which to calculate the cosmic age.
+        a (float, optional): The scale factor at which to calculate the cosmic age.
+                             Only one of `z` or `a` should be provided.
+
+    Returns:
+        float: The age of the universe in Gyr at the specified redshift or scale factor.
+    """
+    cosmo = CosmoDistances(param)
+    return cosmo.age(z=z, a=a)
+
+def Hubble_distance(param):
+    """
+    Calculate the Hubble distance.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+
+    Returns:
+        float: The Hubble distance (c / H0) in Mpc.
+    """
+    cosmo = CosmoDistances(param)
+    return cosmo.Hubble_dist()
+
+def comoving_distance(param, z=None, a=None):
+    """
+    Calculate the comoving distance at a given redshift or scale factor.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+        z (float, optional): The redshift at which to calculate the comoving distance.
+        a (float, optional): The scale factor at which to calculate the comoving distance.
+                             Only one of `z` or `a` should be provided.
+
+    Returns:
+        float: The comoving distance in Mpc at the specified redshift or scale factor.
+    """
+    cosmo = CosmoDistances(param)
+    return cosmo.comoving_dist(z=z, a=a)
+
+def proper_distance(param, z=None, a=None):
+    """
+    Calculate the proper distance at a given redshift or scale factor.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+        z (float, optional): The redshift at which to calculate the proper distance.
+        a (float, optional): The scale factor at which to calculate the proper distance.
+                             Only one of `z` or `a` should be provided.
+
+    Returns:
+        float: The proper distance in Mpc at the specified redshift or scale factor.
+    """
+    cosmo = CosmoDistances(param)
+    return cosmo.proper_dist(z=z, a=a)
+
+def light_travel_distance(param, z=None, a=None):
+    """
+    Calculate the light travel distance (also known as lookback distance) at a given redshift or scale factor.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+        z (float, optional): The redshift at which to calculate the light travel distance.
+        a (float, optional): The scale factor at which to calculate the light travel distance.
+                             Only one of `z` or `a` should be provided.
+
+    Returns:
+        float: The light travel distance in Mpc at the specified redshift or scale factor.
+    """
+    cosmo = CosmoDistances(param)
+    return cosmo.light_travel_dist(z=z, a=a)
+
+def angular_distance(param, z=None, a=None):
+    """
+    Calculate the angular diameter distance at a given redshift or scale factor.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+        z (float, optional): The redshift at which to calculate the angular distance.
+        a (float, optional): The scale factor at which to calculate the angular distance.
+                             Only one of `z` or `a` should be provided.
+
+    Returns:
+        float: The angular diameter distance in Mpc at the specified redshift or scale factor.
+    """
+    cosmo = CosmoDistances(param)
+    return cosmo.angular_dist(z=z, a=a)
+
+def luminosity_distance(param, z=None, a=None):
+    """
+    Calculate the luminosity distance at a given redshift or scale factor.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+        z (float, optional): The redshift at which to calculate the luminosity distance.
+        a (float, optional): The scale factor at which to calculate the luminosity distance.
+                             Only one of `z` or `a` should be provided.
+
+    Returns:
+        float: The luminosity distance in Mpc at the specified redshift or scale factor.
+    """
+    cosmo = CosmoDistances(param)
+    return cosmo.luminosity_dist(z=z, a=a)
+
+def horizon_distance(param):
+    """
+    Calculate the horizon distance, which is the maximum distance from which light has traveled to the observer
+    in the age of the universe.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+
+    Returns:
+        float: The horizon distance in Mpc.
+    """
+    cosmo = CosmoDistances(param)
+    return cosmo.horizon_dist()
+
+def distance_modulus(param, z=None, a=None):
+    """
+    Calculate the distance modulus at a given redshift or scale factor.
+
+    Parameters:
+        param (dict): A dictionary containing cosmological parameters.
+        z (float, optional): The redshift at which to calculate the luminosity distance.
+        a (float, optional): The scale factor at which to calculate the luminosity distance.
+                             Only one of `z` or `a` should be provided.
+
+    Returns:
+        float: The distance modulus in Mpc at the specified redshift or scale factor.
+    """
+    if z is None: z = a_to_z(a)
+    lumdist = lambda z: luminosity_distance(param, z)
+    distmod = lambda z: 5*np.log10(lumdist(z))+25
+    return distmod(z)
diff --git a/src/AstronomyCalc/cosmo_equations.py b/src/AstronomyCalc/cosmo_equations.py
@@ -94,143 +94,3 @@ def luminosity_dist(self, z=None, a=None):
 
     def horizon_dist(self):
         return self.comoving_dist(self, a=1e-8)
-
-
-def Hubble(param, z=None, a=None):
-    """
-    Calculate the Hubble parameter at a given redshift or scale factor.
-
-    Parameters:
-        param (dict): A dictionary containing cosmological parameters.
-        z (float, optional): The redshift at which to evaluate the Hubble parameter.
-        a (float, optional): The scale factor at which to evaluate the Hubble parameter.
-                             Only one of `z` or `a` should be provided.
-
-    Returns:
-        float: The Hubble parameter (H) at the specified redshift or scale factor.
-    """
-    cosmo = CosmoDistances(param)
-    return cosmo.H(z=z, a=a)
-
-def cosmic_age(param, z=None, a=None):
-    """
-    Calculate the cosmic age at a given redshift or scale factor.
-
-    Parameters:
-        param (dict): A dictionary containing cosmological parameters.
-        z (float, optional): The redshift at which to calculate the cosmic age.
-        a (float, optional): The scale factor at which to calculate the cosmic age.
-                             Only one of `z` or `a` should be provided.
-
-    Returns:
-        float: The age of the universe in Gyr at the specified redshift or scale factor.
-    """
-    cosmo = CosmoDistances(param)
-    return cosmo.age(z=z, a=a)
-
-def Hubble_distance(param):
-    """
-    Calculate the Hubble distance.
-
-    Parameters:
-        param (dict): A dictionary containing cosmological parameters.
-
-    Returns:
-        float: The Hubble distance (c / H0) in Mpc.
-    """
-    cosmo = CosmoDistances(param)
-    return cosmo.Hubble_dist()
-
-def comoving_distance(param, z=None, a=None):
-    """
-    Calculate the comoving distance at a given redshift or scale factor.
-
-    Parameters:
-        param (dict): A dictionary containing cosmological parameters.
-        z (float, optional): The redshift at which to calculate the comoving distance.
-        a (float, optional): The scale factor at which to calculate the comoving distance.
-                             Only one of `z` or `a` should be provided.
-
-    Returns:
-        float: The comoving distance in Mpc at the specified redshift or scale factor.
-    """
-    cosmo = CosmoDistances(param)
-    return cosmo.comoving_dist(z=z, a=a)
-
-def proper_distance(param, z=None, a=None):
-    """
-    Calculate the proper distance at a given redshift or scale factor.
-
-    Parameters:
-        param (dict): A dictionary containing cosmological parameters.
-        z (float, optional): The redshift at which to calculate the proper distance.
-        a (float, optional): The scale factor at which to calculate the proper distance.
-                             Only one of `z` or `a` should be provided.
-
-    Returns:
-        float: The proper distance in Mpc at the specified redshift or scale factor.
-    """
-    cosmo = CosmoDistances(param)
-    return cosmo.proper_dist(z=z, a=a)
-
-def light_travel_distance(param, z=None, a=None):
-    """
-    Calculate the light travel distance (also known as lookback distance) at a given redshift or scale factor.
-
-    Parameters:
-        param (dict): A dictionary containing cosmological parameters.
-        z (float, optional): The redshift at which to calculate the light travel distance.
-        a (float, optional): The scale factor at which to calculate the light travel distance.
-                             Only one of `z` or `a` should be provided.
-
-    Returns:
-        float: The light travel distance in Mpc at the specified redshift or scale factor.
-    """
-    cosmo = CosmoDistances(param)
-    return cosmo.light_travel_dist(z=z, a=a)
-
-def angular_distance(param, z=None, a=None):
-    """
-    Calculate the angular diameter distance at a given redshift or scale factor.
-
-    Parameters:
-        param (dict): A dictionary containing cosmological parameters.
-        z (float, optional): The redshift at which to calculate the angular distance.
-        a (float, optional): The scale factor at which to calculate the angular distance.
-                             Only one of `z` or `a` should be provided.
-
-    Returns:
-        float: The angular diameter distance in Mpc at the specified redshift or scale factor.
-    """
-    cosmo = CosmoDistances(param)
-    return cosmo.angular_dist(z=z, a=a)
-
-def luminosity_distance(param, z=None, a=None):
-    """
-    Calculate the luminosity distance at a given redshift or scale factor.
-
-    Parameters:
-        param (dict): A dictionary containing cosmological parameters.
-        z (float, optional): The redshift at which to calculate the luminosity distance.
-        a (float, optional): The scale factor at which to calculate the luminosity distance.
-                             Only one of `z` or `a` should be provided.
-
-    Returns:
-        float: The luminosity distance in Mpc at the specified redshift or scale factor.
-    """
-    cosmo = CosmoDistances(param)
-    return cosmo.luminosity_dist(z=z, a=a)
-
-def horizon_distance(param):
-    """
-    Calculate the horizon distance, which is the maximum distance from which light has traveled to the observer
-    in the age of the universe.
-
-    Parameters:
-        param (dict): A dictionary containing cosmological parameters.
-
-    Returns:
-        float: The horizon distance in Mpc.
-    """
-    cosmo = CosmoDistances(param)
-    return cosmo.horizon_dist()
diff --git a/src/AstronomyCalc/create_data.py b/src/AstronomyCalc/create_data.py
@@ -83,14 +83,64 @@ def Hubble1929_data(data_link=None):
     # print(df.head())
     return np.array(df['distance']), np.array(df['velocity'])
 
-def PantheonPlus_distance_modulus():
+def PantheonPlus_distance_modulus(zmin=0.023, zmax=6, zval='hd'):
+    #string values for the redshift parameter
+    zval_keys = ['hd', 'helio', 'cmb']
+    assert zval in zval_keys
+
     package_folder = str(pkg_resources.files('AstronomyCalc').joinpath('input_data'))
     data_folder = os.path.join(package_folder, "PantheonPlus")
     mean_file = os.path.join(data_folder, "Pantheon+SH0ES.dat")
-    cov_file = os.path.join(data_folder, "Pantheon+SH0ES_STAT+SYS.cov")
-    mean = np.loadtxt(mean_file)
-    cov  = np.loadtxt(cov_file)
-    return {'mean': mean, 'cov': cov}
+    cmat_file = os.path.join(data_folder, "Pantheon+SH0ES_STAT+SYS.cov")
+
+    ## Data vector ## 
+    panthplus = np.loadtxt(mean_file, dtype=str)[1:]
+    panthcols = np.loadtxt(mean_file, dtype=str)[0]
+    hd_indexx = np.where(panthcols == 'zHD')[0][0]
+    cmb_indexx = np.where(panthcols == 'zCMB')[0][0]
+    hel_indexx = np.where(panthcols == 'zHEL')[0][0]
+    if zval == 'hd':
+        zindexx = hd_indexx
+    elif zval == 'helio':
+        zindexx = hel_indexx
+    elif zval == 'cmb':
+        zindexx = cmb_indexx
+
+    nocalib_cond = (panthplus[:, np.where(panthcols == 'IS_CALIBRATOR')[0][0]] == '0')  & (panthplus[:, zindexx].astype('float32') >= zmin) & (panthplus[:, np.where(panthcols == 'zHD')[0][0]].astype('float32') <= zmax)
+    newcal_cond =  (panthplus[:, np.where(panthcols == 'IS_CALIBRATOR')[0][0]] == '1') | ((panthplus[:, zindexx].astype('float32') >= zmin) & (panthplus[:, np.where(panthcols == 'zHD')[0][0]].astype('float32') <= zmax))     
+    hf_cond =  (panthplus[:, zindexx].astype('float32') >= zmin) & (panthplus[:, np.where(panthcols == 'zHD')[0][0]].astype('float32') <= zmax)
+
+    panth_noshoes = panthplus[nocalib_cond]
+    panth_newcal = panthplus[newcal_cond]
+    panth_nocal = panthplus[hf_cond]
+
+    if zval == 'hd':
+        zarr = panth_noshoes[:,hd_indexx].astype('float32')
+        znew = panth_newcal[:,hd_indexx].astype('float32')
+    elif zval == 'helio':
+        zarr = panth_noshoes[:,hel_indexx].astype('float32')
+        znew = panth_newcal[:,hel_indexx].astype('float32')
+    elif zval == 'cmb':
+        zarr = panth_noshoes[:,cmb_indexx].astype('float32')
+        znew = panth_newcal[:,cmb_indexx].astype('float32')
+
+    zarr2 = panth_nocal[:,hd_indexx].astype('float32')
+    mu_hd_newcal = panth_newcal[:,np.where(panthcols == 'm_b_corr')[0][0]].astype('float32')
+
+    calibcol_index = np.where(panthcols == 'IS_CALIBRATOR')[0][0]
+    mu_shoes = panth_newcal[panth_newcal[:, calibcol_index] == '1'][:,np.where(panthcols == 'MU_SH0ES')[0][0]].astype('float32')
+
+    ## Covariance matrix ##
+    cov_mat_sys = np.loadtxt(cmat_file, skiprows=1)
+    cov_mat_sys_shape = cov_mat_sys.reshape(len(panthplus), len(panthplus))
+    cond_indices = np.where(nocalib_cond)[0]
+
+    cond_hf_indices = np.where(hf_cond)[0]
+    cond_newcal_indices = np.where(newcal_cond)[0]
+    cov_mat_newcal = cov_mat_sys_shape[cond_newcal_indices][:, cond_newcal_indices]
+    c_inv_newcal = np.linalg.inv(cov_mat_newcal)
+
+    return {'z': znew, 'data': mu_hd_newcal, 'cov': cov_mat_newcal, 'cov_inv': c_inv_newcal}
 
 def SPARC_galaxy_rotation_curves_data(filename=None, name=None):
     """