ENH: Subhalo abundance (SHAM) code

skypy/halos/__init__.py

@@ -14,9 +14,18 @@
 
 __all__ = [
     'colossus_mf',
+    'quenching_funct',
+    'find_min',
+    'run_file',
+    'gen_sub_cat',
+    'galaxy_cat',
+    'assignment',
+    'sham_arrays',
+    'run_sham',
 ]
 
 from . import abundance_matching  # noqa F401,F403
 from . import mass  # noqa F401,F403
 from . import quenching  # noqa F401,F403
 from ._colossus import colossus_mf  # noqa F401,F403
+from . import sham # noqa F401,F403

skypy/halos/bc_file.yaml

@@ -0,0 +1,20 @@
+m_star: !numpy.power [10, 10.65]
+phi_star: !numpy.power [10, -2.98]
+alpha_val: -1.48
+m_min: !numpy.power [10, 7.0]
+m_max: !numpy.power [10, 14.0]
+sky_area: 100.0 deg2
+z_range: !numpy.linspace [0.010004503125093451, 0.09999991945507442, 100]
+cosmology: !astropy.cosmology.FlatLambdaCDM
+  H0: 70.0
+  Om0: 0.3
+tables:
+  galaxy:
+    z, sm: !skypy.galaxies.schechter_smf
+      redshift: $z_range
+      m_star: $m_star
+      phi_star: $phi_star
+      alpha: $alpha_val
+      m_min: $m_min
+      m_max: $m_max
+      sky_area: $sky_area

skypy/halos/bs_file.yaml

@@ -0,0 +1,20 @@
+m_star: !numpy.power [10, 10.55]
+phi_star: !numpy.power [10, -3.96]
+alpha_val: -1.53
+m_min: !numpy.power [10, 7.0]
+m_max: !numpy.power [10, 14.0]
+sky_area: 100.0 deg2
+z_range: !numpy.linspace [0.010004503125093451, 0.09999991945507442, 100]
+cosmology: !astropy.cosmology.FlatLambdaCDM
+  H0: 70.0
+  Om0: 0.3
+tables:
+  galaxy:
+    z, sm: !skypy.galaxies.schechter_smf
+      redshift: $z_range
+      m_star: $m_star
+      phi_star: $phi_star
+      alpha: $alpha_val
+      m_min: $m_min
+      m_max: $m_max
+      sky_area: $sky_area

skypy/halos/halo.yaml

@@ -0,0 +1,26 @@
+parameters:
+  model: 'sheth99'
+  mdef: 'fof'
+  m_min: 1.e+9
+  m_max: 1.e+15
+  sky_area: 600. deg2 # Same as galaxies
+  sigma_8: 0.8
+  ns: 0.96
+cosmology: !astropy.cosmology.FlatLambdaCDM
+  H0: 70
+  Om0: 0.3
+  name: 'FlatLambdaCDM' # Requires user defined name
+  Ob0: 0.045
+  Tcmb0: 2.7
+z_range: !numpy.linspace [0.01, 0.1, 100] # Same as galaxies 
+tables:
+  halo:
+    z, mass: !skypy.halos.colossus_mf
+      redshift: $z_range
+      model: $model
+      mdef: $mdef
+      m_min: $m_min
+      m_max: $m_max
+      sky_area: $sky_area
+      sigma8: $sigma_8
+      ns: $ns

skypy/halos/mass.py

@@ -348,7 +348,7 @@ def number_subhalos(halo_mass, alpha, beta, gamma_M, x, m_min, noise=True):
         subhalos.
     m_min : array_like
         Original mass of the least massive subhalo, in units of solar mass.
-        Current stripped mass is given by :math:`x m_{min}`.
+        Current stripped mass is given by :math:`m_{min}/x`.
     noise : bool, optional
         Poisson-sample the number of subhalos. Default is `True`.
 
@@ -383,6 +383,7 @@ def number_subhalos(halo_mass, alpha, beta, gamma_M, x, m_min, noise=True):
 
     # Random number of subhalos following a Poisson distribution
     # with mean n_subhalos
+    # TODO Poisson rounds the value, otherwise doesn't. mass_sampler needs rounded numbers
     return np.random.poisson(n_subhalos) if noise else n_subhalos
 
 
@@ -440,9 +441,12 @@ def subhalo_mass_sampler(halo_mass, nsubhalos, alpha, beta,
     halo_mass = np.atleast_1d(halo_mass)
     nsubhalos = np.atleast_1d(nsubhalos)
     subhalo_list = []
+    # TODO fix for not using Poisson in above function
+    # nsubhalos = [int(x) for x in nsubhalos]
+
     for M, n in zip(halo_mass, nsubhalos):
         x_min = m_min / (x * beta * M)
         x_max = 0.5 / (x * beta)
-        subhalo_mass = schechter(alpha, x_min, x_max, resolution, size=n, scale=x * beta * M)
+        subhalo_mass = schechter(-alpha, x_min, x_max, resolution, size=n, scale=x * beta * M)
         subhalo_list.append(subhalo_mass)
     return np.concatenate(subhalo_list)

skypy/halos/rc_file.yaml

@@ -0,0 +1,20 @@
+m_star: !numpy.power [10, 10.58]
+phi_star: !numpy.power [10, -2.77]
+alpha_val: -0.33
+m_min: !numpy.power [10, 7.0]
+m_max: !numpy.power [10, 14.0]
+sky_area: 100.0 deg2
+z_range: !numpy.linspace [0.010004503125093451, 0.09999991945507442, 100]
+cosmology: !astropy.cosmology.FlatLambdaCDM
+  H0: 70.0
+  Om0: 0.3
+tables:
+  galaxy:
+    z, sm: !skypy.galaxies.schechter_smf
+      redshift: $z_range
+      m_star: $m_star
+      phi_star: $phi_star
+      alpha: $alpha_val
+      m_min: $m_min
+      m_max: $m_max
+      sky_area: $sky_area

skypy/halos/red_central.yaml

@@ -0,0 +1,20 @@
+m_star: !numpy.power [10, 10.75]
+phi_star: !numpy.power [10, -2.37]
+alpha_val: -0.18
+m_min: !numpy.power [10, 7.0]
+m_max: !numpy.power [10, 14.0]
+sky_area: 600.0 deg2
+z_range: !numpy.linspace [0.01, 0.1, 100]
+cosmology: !astropy.cosmology.FlatLambdaCDM
+  H0: 70.0
+  Om0: 0.3
+tables:
+  galaxy:
+    z, sm: !skypy.galaxies.schechter_smf
+      redshift: $z_range
+      m_star: $m_star
+      phi_star: $phi_star
+      alpha: $alpha_val
+      m_min: $m_min
+      m_max: $m_max
+      sky_area: $sky_area

skypy/halos/rs_file.yaml

@@ -0,0 +1,20 @@
+m_star: !numpy.power [10, 10.64]
+phi_star: !numpy.power [10, -4.24]
+alpha_val: -1.54
+m_min: !numpy.power [10, 7.0]
+m_max: !numpy.power [10, 14.0]
+sky_area: 100.0 deg2
+z_range: !numpy.linspace [0.010004503125093451, 0.09999991945507442, 100]
+cosmology: !astropy.cosmology.FlatLambdaCDM
+  H0: 70.0
+  Om0: 0.3
+tables:
+  galaxy:
+    z, sm: !skypy.galaxies.schechter_smf
+      redshift: $z_range
+      m_star: $m_star
+      phi_star: $phi_star
+      alpha: $alpha_val
+      m_min: $m_min
+      m_max: $m_max
+      sky_area: $sky_area