bug fixing halo model theory for cobaya

soliket/halo_model/halo_model_fe/HaloModel_fe.yaml

@@ -0,0 +1,32 @@
+read_matterPS: False #if True, reads pre-computed linear matter PS
+                        #if False, computes the linear matter PS using CAMB
+gal_mod: False
+matter_matter: 'mm'
+galaxy_galaxy: #'gg'
+matter_galaxy: "gm"
+redshift: './tabulated/redshift.txt'
+
+z : [0,1,2]
+
+kmax: 10.0
+
+Mmin: 11
+Mmax: 15
+nm: 200 
+
+Dc: 200
+
+T_CMB: 2.725
+tau: 0.0544
+ns: 0.9649
+As: 1.97448e-9 #2.101e-9
+pivot_scalar: 0.05
+matter_PS:
+sigma_EP: 0.1
+sigma_LP: 0.1
+scale_EP: 10.0
+scale_LP: 10.0
+alpha_EP: 1.0
+alpha_LP: 1.0
+LogMmin_EP: 12.0
+LogMmin_LP: 10.8

soliket/halo_model/halo_model_fe/__init__.py

@@ -1,5 +1,5 @@
-from .halo_model_fe import paramfile_halomod
+from .halo_model_fe import HaloModel_fe
 
 __all__ = {
-    "paramfile_halomod"
+    "HaloModel_fe"
 }

soliket/halo_model/halo_model_fe/halo_model_fe.py

@@ -1,10 +1,10 @@
 import numpy as np
 from cobaya.theory import Theory
-from utils import *
-from HODS import *
-from power_spectrum import *
-from lin_matterPS import * 
-from cosmology import *
+from .utils import *
+from .HODS import *
+from .power_spectrum import *
+from .lin_matterPS import * 
+from .cosmology import *
 
 class HaloModel(Theory):
     _logz = np.linspace(-3, np.log10(1100), 150)
@@ -37,7 +37,8 @@ def get_Pk_gm_grid(self):
 class HaloModel_fe(HaloModel):
     def initialize(self):
         super().initialize()
-        self.logmass = np.linspace(self.Mmin, self.mmax, self.nm)
+        self.logmass = np.linspace(self.Mmin, self.Mmax, self.nm)
+        #self.k = np.linspace(0, self.kmax, 1)
         self.clust_param = {'sigma_EP': self.sigma_EP,
                             'sigma_LP': self.sigma_LP,
                             'scale_EP': self.scale_EP,
@@ -48,8 +49,8 @@ def initialize(self):
                             'LogMmin_LP': self.LogMmin_LP,
                             'Dc': self.Dc,
                             }
-        self.isinstance_200 = u_p_nfw_hmf_bias(self.k, self._get_Pk_mm_lin(), self.logmass, self.z, self.Dc)
-        self.instance_HOD = hod_ngal(self.logmass, self.z, self.clust_param, self.instance_200)
+        #self.instance_200 = u_p_nfw_hmf_bias(self.k, self._get_Pk_mm_lin(), self.logmass, self.z, self.Dc)
+        #self.instance_HOD = hod_ngal(self.logmass, self.z, self.clust_param, self.instance_200)
 
     def must_provide(self, **requirements):
 
@@ -79,10 +80,16 @@ def must_provide(self, **requirements):
                            }
         return needs
 
-    def calculate(self, state: dict):
+    def calculate(self, state: dict, want_derived: bool = True,
+                  **params_values_dict):
+        Pk_mm_lin = self._get_Pk_mm_lin()
+
+        self.instance_200 = u_p_nfw_hmf_bias(self.k, Pk_mm_lin, self.logmass, self.z, self.Dc)
+        self.instance_HOD = hod_ngal(self.logmass, self.z, self.clust_param, self.instance_200)
+
         spectra = mm_gg_mg_spectra(
                     self.k,
-                    self._get_Pk_mm_lin(),
+                    Pk_mm_lin,
                     self.logmass,
                     self.z,
                     self.instance_HOD,

soliket/halo_model/halo_model_fe/utils.py

@@ -54,7 +54,8 @@ def sigma(self, rad, red, zeta):
         lnk = np.log(k)
         uW = self.W(rk)
         integ = rest * uW ** 2
-        sigm = (0.5 / np.pi ** 2) * scipy.integrate.simps(integ, x=lnk, axis=-1)
+        sigm = (0.5 / np.pi ** 2) * trapz(integ, x=lnk, axis=-1) ####################################### gz
+        #sigm = (0.5 / np.pi ** 2) * scipy.integrate.simps(integ, x=lnk, axis=-1)
         return np.sqrt(sigm)
 
     # sigma depends on z from the linear power spectrum
@@ -128,7 +129,10 @@ def dlns2_dlnr(self, rad, red, zeta):
         inte = w * dw * rest
         lnk = np.log(k)
         s = self.sigma(rad, red, zeta)
-        return scipy.integrate.simps(inte, x=lnk, axis=-1, even="avg") / (
+        #return scipy.integrate.simps(inte, x=lnk, axis=-1, even="avg") / (
+           # np.pi ** 2 * s ** 2
+        #)
+        return trapz(inte, x=lnk, axis=-1) / (
             np.pi ** 2 * s ** 2
         )