transverse distance fixed

Signed-off-by: Anto Idicherian Lonappan <[email protected]>
antolonappan · Sep 2, 2023 · 3c00469 · 3c00469
1 parent 28403a8
commit 3c00469
Show file tree

Hide file tree

Showing 6 changed files with 244 additions and 194 deletions.
diff --git a/Notebooks/Validation_cosmo.ipynb b/Notebooks/Validation_cosmo.ipynb
@@ -24,19 +24,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 3,
    "id": "a53c96d1-5fb0-42ed-b5f0-4f0e1d88f856",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "The autoreload extension is already loaded. To reload it, use:\n",
-      "  %reload_ext autoreload\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "%load_ext autoreload\n",
     "%autoreload 2\n",
@@ -58,7 +49,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 6,
    "id": "b53aa1cb-b55c-4ef2-b814-ce68fd25098f",
    "metadata": {},
    "outputs": [],
@@ -355,10 +346,123 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 67,
    "id": "14d103a7",
    "metadata": {},
    "outputs": [],
+   "source": [
+    "lcdm = Cosmology('LCDM',['H0','Omega_m'])\n",
+    "wcdm = Cosmology('wCDM',['H0','Omega_m','w0'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 73,
+   "id": "44681cac",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([   0.        , 3563.00252552, 5480.86846219, 6652.2955333 ,\n",
+       "       7453.46122345, 8043.88164985, 8501.73372259, 8870.0167495 ,\n",
+       "       9174.50291195, 9431.67380384])"
+      ]
+     },
+     "execution_count": 73,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "lcdm._transverse_distance_([70, 0.3], zlist)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 87,
+   "id": "53990820",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([   0.        , 3756.12021805, 5743.63627523, 6933.36832763,\n",
+       "       7740.79000998, 8333.81846036, 8792.9217005 , 9161.87050455,\n",
+       "       9466.74024632, 9724.14483018])"
+      ]
+     },
+     "execution_count": 87,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "wcdm._transverse_distance_( [70, 0.3, -1.3], zlist)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 78,
+   "id": "b4612f32",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([   0.        , 3756.12021805, 5743.63627523, 6933.36832763,\n",
+       "       7740.79000998, 8333.81846036, 8792.9217005 , 9161.87050455,\n",
+       "       9466.74024632, 9724.14483018])"
+      ]
+     },
+     "execution_count": 78,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "wcdm._transverse_distance_( [70, 0.3, -1.3], zlist)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "id": "b0ea5a83",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "klcdm = Cosmology('kLCDM',['H0','Omega_m','Omega_k'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 36,
+   "id": "407eb913",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([4.28157123e+00, 3.49909018e+03, 5.42973331e+03, 6.66479343e+03,\n",
+       "       7.54168971e+03, 8.20704907e+03, 8.73514182e+03, 9.16806810e+03,\n",
+       "       9.53172140e+03, 9.84303874e+03])"
+      ]
+     },
+     "execution_count": 36,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "klcdm.transverse_distance( [70, 0.3, 0.1], zlist)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b4d2effa",
+   "metadata": {},
+   "outputs": [],
    "source": []
   }
  ],

diff --git a/marcia/__init__.py b/marcia/__init__.py
@@ -6,4 +6,4 @@
 from .likelihood import Likelihood
 from .sampler import Sampler
 from .kernel import *
-from .likelihood_GP import Likelihood_GP
+from .likelihood import Likelihood_GP
diff --git a/marcia/backend/cosmology.py b/marcia/backend/cosmology.py
@@ -5,6 +5,7 @@
 @njit([f8[:](f8,f8,f8,f8,f8[:],f8[:])])
 def hubble_rate(H0,Omega_m,Omega_b,Omega_k,de,z):
     Omega_r = 4.18343*10**-5./(H0/100.)**2.
+    #print('backend_H_r',(1. - Omega_m - Omega_k - Omega_b - Omega_r)*de)
     E2 = Omega_r*((1. + z)**4) + Omega_m*((1. + z)**3) + Omega_b*((1. + z)**3) + Omega_k*((1. + z)**2) + (1. - Omega_m - Omega_k - Omega_b - Omega_r)*de
     Hofz = H0*np.sqrt(E2)
     return Hofz
@@ -44,14 +45,16 @@ def z_inp(z):
 def interpolate(z_inp,z,func):
     return np.interp(z,z_inp,func)
 
-@njit
-def Ly(y,t,H0,Omega_m,Omega_b,Omega_k,de):
-    return inv_hubble_rate(H0,Omega_m,Omega_b,Omega_k,de,np.array([t]))[0]
 
 
-def transverse_distance(H0,Omega_m,Omega_b,Omega_k,de,z,clight):
-    y=odeint(Ly,0.0,z,args=(H0,Omega_m,Omega_b,Omega_k,de))
-    return clight* y[:,0]
+@njit
+def transverse_distance(H0,clight,Omega_k,y):
+    if Omega_k > 0.0:
+        return clight/(np.sqrt(Omega_k)*H0)*np.sinh(np.sqrt(Omega_k) * H0*y)
+    elif Omega_k < 0.0:
+        return clight/(np.sqrt(np.abs(Omega_k))*H0)*np.sin((np.sqrt(np.abs(Omega_k)))* H0* y)
+    elif Omega_k == 0.0:
+        return clight* y
 
 @njit
 def distance_modulus(Mb,z2,d):

diff --git a/marcia/cosmology.py b/marcia/cosmology.py
@@ -193,6 +193,7 @@ def __check_mandatory_parameters__(self):
 
     def dark_energy_f(self, parameters, z):
         p = self.param(parameters)
+        #print('dark_energy_f',p.w0)
         return cbackend.dark_energy_f_wCDM(p.w0, p.wa, z)
 
     def dark_energy_w(self,parameters, z):
@@ -201,19 +202,14 @@ def dark_energy_w(self,parameters, z):
         return p.w0 + p.wa*(1-a)
 
     def _transverse_distance_(self,parameters, z):
-        # Use this print, in case needed for Error_handling.
-        # print ''
+        p = self.param(parameters)
 
-        # def Ly(y,t):
-        #     return self.inv_hubble_rate(parameters,np.array([t]))
+        def Ly(y,t):
+            return self.inv_hubble_rate(parameters,np.array([t]))
 
-        # y=it.odeint(Ly,0.0,z)
-        # tint = np.array(y[:,0])
-        # return self.clight* tint
+        y=it.odeint(Ly,0.0,z)
+        return cbackend.transverse_distance(p.H0,self.clight,p.Omega_k,y[:,0])
 
-        p = self.param(parameters)
-        de = self.dark_energy_f(parameters, z)
-        return cbackend.transverse_distance(p.H0,p.Omega_m,p.Omega_b,p.Omega_k,de, z,self.clight)
 
 class LCDM(wCDM):
     def __init__(self, parameters,prior_file=None):
@@ -326,22 +322,6 @@ def __check_mandatory_parameters__(self):
         assert len(self.param.parameters) == n, 'kwCDM: parameters are not correct'
 
 
-    def transverse_distance(self, parameters, z):
-        p = self.param(parameters)
-        def Ly(y,t):
-            return 1./self.hubble_rate(parameters,t)
-
-        y=it.odeint(Ly,0.0,z)
-        tint = np.array(y[:,0])
-        if p.Omega_k > 0.0:
-            return np.nan_to_num( self.clight/(np.sqrt(p.Omega_k)*p.H0)*np.sinh(np.sqrt(p.Omega_k) * p.H0*tint))
-
-        elif p.Omega_k < 0.0:
-            return np.nan_to_num( self.clight/(np.sqrt(abs(p.Omega_k))*p.H0)*np.sin((np.sqrt(abs(p.Omega_k)))* p.H0* tint))
-
-        elif p.Omega_k == 0.0:
-            return np.nan_to_num( self.clight* tint)
-
 
 class kLCDM(kwCDM):
     def __init__(self, parameters,prior_file=None):
@@ -394,21 +374,6 @@ def __check_mandatory_parameters__(self):
             n+=1
         assert len(self.param.parameters) == n, 'kCPL3: parameters are not correct'
 
-    def transverse_distance(self, parameters, z):
-        p = self.param(parameters)
-        def Ly(y,t):
-            return 1./self.hubble_rate(parameters,t)
-
-        y=it.odeint(Ly,0.0,z)
-        tint = np.array(y[:,0])
-        if p.Omega_k > 0.0:
-            return np.nan_to_num( self.clight/(np.sqrt(p.Omega_k)*p.H0)*np.sinh(np.sqrt(p.Omega_k) * p.H0*tint))
-
-        elif p.Omega_k < 0.0:
-            return np.nan_to_num( self.clight/(np.sqrt(abs(p.Omega_k))*p.H0)*np.sin((np.sqrt(abs(p.Omega_k)))* p.H0* tint))
-
-        elif p.Omega_k == 0.0:
-            return np.nan_to_num( self.clight* tint)
 
 
 class kXCDM(XCDM):
@@ -428,21 +393,5 @@ def __check_mandatory_parameters__(self):
         if self.Mbsample:
             n+=1
         assert len(self.param.parameters) == n, 'kXCDM: parameters are not correct'
-
-    def transverse_distance(self, parameters, z):
-        """transverse distance in Mpc/h"""
-        p = self.param(parameters)
-        def Ly(y,t):
-            return 1./self.hubble_rate(parameters,t)
 
-        y=it.odeint(Ly,0.0,z)
-        tint = np.array(y[:,0])
-        if p.Omega_k > 0.0:
-            return np.nan_to_num( self.clight/(np.sqrt(p.Omega_k)*p.H0)*np.sinh(np.sqrt(p.Omega_k) * p.H0*tint))
-
-        elif p.Omega_k < 0.0:
-            return np.nan_to_num( self.clight/(np.sqrt(abs(p.Omega_k))*p.H0)*np.sin((np.sqrt(abs(p.Omega_k)))* p.H0* tint))
-
-        elif p.Omega_k == 0.0:
-            return np.nan_to_num( self.clight* tint)