Dev hm fe

soliket/halo_model_fe/HODS.py

@@ -0,0 +1,97 @@
+import numpy as np
+from numpy import trapz
+from scipy import special
+
+
+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_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_fe/__init__.py

@@ -0,0 +1 @@
+from .halo_model_fe import HaloModel, HaloModel_fe

soliket/halo_model_fe/halo_model_fe.py

@@ -0,0 +1,115 @@
+import numpy as np
+from cobaya.theory import Theory
+
+from .HODS import hod_ngal
+from .power_spectrum import mm_gg_mg_spectra
+from .utils import u_p_nfw_hmf_bias
+
+
+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):
+        if not self._var_pairs:
+            raise ValueError(
+                "Error: `_var_pairs` is empty. "
+                "Check `must_provide()` for missing dependencies."
+            )
+        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(
+            {(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()
+        Pmat = spectra.halo_terms_matter()
+
+        state["Pk_gg_grid"] = Pgal.T
+        state["Pk_mm_grid"] = Pmat.T