Dev hm fe

soliket/halo_model/halo_model_fe/HODS.py

@@ -0,0 +1,81 @@
+import numpy as np
+from scipy import special
+from numpy import trapz
+
+
+class hod_ngal:
+    def __init__(self, mh, redshift, clust_param, instance_200):
+        self.mh = mh
+        self.redshift = redshift
+        self.clust_param = clust_param
+        self.instance_200 = instance_200
+
+        self.HODS_EP()
+        self.mean_gal_EP()
+        self.HODS_LP()
+        self.mean_gal_LP()
+
+    def HODS_EP(self):
+        Ncent = np.zeros([len(self.mh)])
+        Nsat = np.zeros([len(self.mh)])
+        Nbra = np.zeros([len(self.mh)])
+
+        Mmin = 10 ** self.clust_param['LogMmin_EP']
+        Msat = self.clust_param['scale_EP'] * Mmin
+        Ncent = 0.5 * (1 + special.erf((np.log10(self.mh) - 
+                                        np.log10(Mmin)) / self.clust_param['sigma_EP']))
+        Nsat = (
+            0.5
+            * (1 + special.erf((np.log10(self.mh) - 
+                                np.log10(2 * Mmin)) / self.clust_param['sigma_EP']))
+            * ((self.mh) / Msat) ** self.clust_param['alpha_EP']
+        )
+        Nbra = Ncent + Nsat
+
+        self.Nbra_EP = Nbra
+        self.Ncent_EP = Ncent
+        self.Nsat_EP = Nsat
+
+        return Ncent, Nsat, Nbra
+
+    def mean_gal_EP(self):
+
+        Nbra = self.HODS_EP()[2]
+        ngal_200c = trapz(self.instance_200.dndM[:, :] * Nbra[np.newaxis, :], self.mh[:])
+
+        self.ngal_EP_200c = ngal_200c
+
+        return
+
+    def HODS_LP(self):
+        Ncent = np.zeros([len(self.mh)])
+        Nsat = np.zeros([len(self.mh)])
+        Nbra = np.zeros([len(self.mh)])
+
+        Mmin = 10 ** self.clust_param['LogMmin_LP']
+        Msat = self.clust_param['scale_LP'] * Mmin
+        Ncent = 0.5 * (1 + special.erf((np.log10(self.mh) - 
+                                        np.log10(Mmin)) / self.clust_param['sigma_LP']))
+        Nsat = (
+            0.5
+            * (1 + special.erf((np.log10(self.mh) - 
+                                np.log10(2 * Mmin)) / self.clust_param['sigma_LP']))
+            * ((self.mh) / Msat) ** self.clust_param['alpha_LP']
+        )
+        Nbra = Ncent + Nsat
+
+        self.Nbra_LP = Nbra
+        self.Ncent_LP = Ncent
+        self.Nsat_LP = Nsat
+
+        return Ncent, Nsat, Nbra
+
+    def mean_gal_LP(self):
+
+        Nbra = self.HODS_LP()[2]
+
+        ngal_200c = trapz(self.instance_200.dndM[:, :] * Nbra[np.newaxis, :], self.mh[:])
+
+        self.ngal_LP_200c = ngal_200c
+
+        return

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

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

soliket/halo_model/halo_model_fe/halo_model_fe.py

@@ -0,0 +1,101 @@
+import numpy as np
+from cobaya.theory import Theory
+from .utils import u_p_nfw_hmf_bias
+from .HODS import hod_ngal
+from .power_spectrum import mm_gg_mg_spectra
+
+
+class HaloModel(Theory):
+    _logz = np.linspace(-3, np.log10(1100), 150)
+    _default_z_sampling = 10**_logz
+    _default_z_sampling[0] = 0
+
+    def initialize(self):
+        self._var_pairs = set()
+        self._required_results = {}
+
+    def _get_Pk_mm_lin(self):
+        for pair in self._var_pairs:
+            self.k, self.z, Pk_mm = \
+                self.provider.get_Pk_grid(var_pair=pair, nonlinear=False)
+
+        return Pk_mm
+
+    def get_Pk_mm_grid(self):
+
+        return self.current_state["Pk_mm_grid"]
+
+    def get_Pk_gg_grid(self):
+
+        return self.current_state["Pk_gg_grid"]
+
+    def get_Pk_gm_grid(self):
+
+        return self.current_state["Pk_gm_grid"]
+
+
+class HaloModel_fe(HaloModel):
+    def initialize(self):
+        super().initialize()
+        self.logmass = np.linspace(self.Mmin, self.Mmax, self.nm)
+        self.clust_param = {'sigma_EP': self.sigma_EP,
+                            'sigma_LP': self.sigma_LP,
+                            'scale_EP': self.scale_EP,
+                            'scale_LP': self.scale_LP,
+                            'alpha_EP': self.alpha_EP,
+                            'alpha_LP': self.alpha_LP,
+                            'LogMmin_EP': self.LogMmin_EP,
+                            'LogMmin_LP': self.LogMmin_LP,
+                            'Dc': self.Dc,
+                            }
+
+    def must_provide(self, **requirements):
+
+        options = requirements.get("halo_model") or {}
+        self._var_pairs.update(
+            set((x, y) for x, y in
+                options.get("vars_pairs", [("delta_tot", "delta_tot")])))
+
+        self.kmax = max(self.kmax, options.get("kmax", self.kmax))
+        self.z = np.unique(np.concatenate(
+                            (np.atleast_1d(options.get("z", self._default_z_sampling)),
+                            np.atleast_1d(self.z))))
+
+        needs = {}
+
+        needs["Pk_grid"] = {
+                "vars_pairs": self._var_pairs,
+                "nonlinear": (False, False),
+                "z": self.z,
+                "k_max": self.kmax
+            }
+
+        needs["sigma_R"] = {"vars_pairs": self._var_pairs,
+                           "z": self.z,
+                           "k_max": self.kmax,
+                           "R": np.linspace(0.14, 66, 256) # list of radii required
+                           }
+        return needs
+
+    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,
+                    Pk_mm_lin,
+                    self.logmass,
+                    self.z,
+                    self.instance_HOD,
+                    self.instance_200,
+                    self.gal_mod
+                )
+
+        Pgal = spectra.halo_terms_galaxy()[0]
+
+        state['Pk_gg_grid'] = Pgal