aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic/clusters/psformalism.py
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim/extragalactic/clusters/psformalism.py')
-rw-r--r--fg21sim/extragalactic/clusters/psformalism.py205
1 files changed, 136 insertions, 69 deletions
diff --git a/fg21sim/extragalactic/clusters/psformalism.py b/fg21sim/extragalactic/clusters/psformalism.py
index 0ed19b0..899d94e 100644
--- a/fg21sim/extragalactic/clusters/psformalism.py
+++ b/fg21sim/extragalactic/clusters/psformalism.py
@@ -6,7 +6,7 @@ Press-Schechter (PS) formalism
First determine the number of clusters within a sky patch (i.e., sky
coverage) according to the cluster distribution predicted by the PS
-formalism; then sampling from the PS mass function to derive the mass
+formalism; then sampling from the halo mass function to derive the mass
and redshift for each cluster.
"""
@@ -15,10 +15,12 @@ import random
import numpy as np
import pandas as pd
+import hmf
from ...share import CONFIGS, COSMO
from ...utils.interpolate import bilinear
from ...utils.units import UnitConversions as AUC
+from ...utils.io import write_dndlnm
logger = logging.getLogger(__name__)
@@ -28,106 +30,171 @@ class PSFormalism:
"""
Press-Schechter (PS) formalism
- Simulate the clusters number and their distribution (mass and z)
- within a sky patch of certain coverage.
+ Calculate the halo mass distribution with respect to mass and redshift,
+ determine the clusters number counts and generate their distribution
+ (mass and z) within a sky patch of certain coverage.
"""
def __init__(self, configs=CONFIGS):
self.configs = configs
self._set_configs()
- self._load_data()
def _set_configs(self):
"""
Load the required configurations and set them.
"""
+ comp = "extragalactic/psformalism"
+ self.model = self.configs.getn(comp+"/model")
+ self.M_min = self.configs.getn(comp+"/M_min")
+ self.M_max = self.configs.getn(comp+"/M_max")
+ self.M_step = self.configs.getn(comp+"/M_step")
+ self.z_min = self.configs.getn(comp+"/z_min")
+ self.z_max = self.configs.getn(comp+"/z_max")
+ self.z_step = self.configs.getn(comp+"/z_step")
+ self.dndlnm_outfile = self.configs.get_path(comp+"/dndlnm_outfile")
+
comp = "extragalactic/clusters"
- self.datafile = self.configs.get_path(comp+"/ps_data")
self.Mmin_cluster = self.configs.getn(comp+"/mass_min") # [Msun]
self.boost = self.configs.getn(comp+"/boost")
+ self.clobber = self.configs.getn("output/clobber")
+
+ @property
+ def hmf_model(self):
+ return {"PS": "PS",
+ "SMT": "SMT",
+ "JENKINS": "Jenkins"}[self.model.upper()]
+
+ def hmf_massfunc(self, z=0.0):
+ """
+ Halo mass function as a `~hmf.MassFunction` instance.
+ """
+ if not hasattr(self, "_hmf_massfunc"):
+ h = COSMO.h
+ cosmo = COSMO._cosmo
+ self._hmf_massfunc = hmf.MassFunction(
+ Mmin=np.log10(self.M_min*h),
+ Mmax=np.log10(self.M_max*h),
+ dlog10m=self.M_step,
+ hmf_model=self.hmf_model,
+ cosmo_model=cosmo,
+ n=COSMO.ns,
+ sigma_8=COSMO.sigma8)
+ logger.info("Initialized '%s' halo mass function." %
+ self.hmf_model)
+
+ massfunc = self._hmf_massfunc
+ massfunc.update(z=z)
+ return massfunc
+
+ @property
+ def z(self):
+ """
+ The redshift points where to calculate the dndlnm data.
+ """
+ return np.arange(self.z_min, self.z_max+self.z_step/2, self.z_step)
+
+ @property
+ def mass(self):
+ """
+ The mass points where to calculate the dndlnm data.
+
+ NOTE:
+ The maximum mass end is exclusive, to be consistent with hmf's
+ mass function!
+ """
+ return 10 ** np.arange(np.log10(self.M_min),
+ np.log10(self.M_max),
+ self.M_step)
+
+ @property
+ def dndlnm(self):
+ """
+ The calculated halo mass distributions data.
+ """
+ if not hasattr(self, "_dndlnm"):
+ self._dndlnm = self.calc_dndlnm()
+ return self._dndlnm
+
+ def calc_dndlnm(self):
+ """
+ Calculate the halo mass distributions expressed in ``dndlnm``,
+ the differential mass distribution in terms of natural log of
+ masses.
+ Unit: [Mpc^-3] (the little "h" is folded into the values)
+
+ NOTE
+ ----
+ dndlnm = d n(M,z) / d ln(M); [Mpc^-3]
+ describes the number of halos per comoving volume (Mpc^3) at
+ redshift z per unit logarithmic mass interval at mass M.
+ """
+ logger.info("Calculating dndlnm data ...")
+ dndlnm = []
+ h = COSMO.h
+ for z_ in self.z:
+ massfunc = self.hmf_massfunc(z_)
+ dndlnm.append(massfunc.dndlnm * h**3)
+ self._dndlnm = np.array(dndlnm)
+ logger.info("Calculated dndlnm within redshift: %.1f - %.1f" %
+ (self.z_min, self.z_max))
+ return self._dndlnm
+
+ def write(self, outfile=None):
+ """
+ Write the calculate dndlnm data into file as NumPy ".npz" format.
+ """
+ if outfile is None:
+ outfile = self.dndlnm_outfile
+ write_dndlnm(outfile, dndlnm=self.dndlnm, z=self.z, mass=self.mass,
+ clobber=self.clobber)
+ logger.info("Wrote dndlnm data into file: %s" % outfile)
+
@property
def Mmin_halo(self):
return self.Mmin_cluster * COSMO.darkmatter_fraction
- def _load_data(self, filepath=None):
- """
- Load dndM data and reformat into a 2D density grid together with
- redshifts and masses vectors.
-
- Data File Description
- ---------------------
- z1 mass1 density1
- z1 mass2 density2
- z1 .. density3
- z2 mass1 density4
- z2 mass2 density5
- z2 .. density6
- ...
-
- where,
- * Redshifts: 0.0 -> 3.02, even-spacing, step 0.02
- * Mass: unit 1e12 -> 9.12e15 [Msun], log-even (dark matter)
- * Density: [number]/dVc/dM
- with,
- - dVc: differential comvoing volume, [Mpc^3]/[sr]/[unit redshift]
- """
- if filepath is None:
- filepath = self.datafile
- data = np.loadtxt(filepath)
- redshifts = data[:, 0]
- masses = data[:, 1]
- densities = data[:, 2]
-
- redshifts = np.array(list(set(redshifts)))
- redshifts.sort()
- masses = np.array(list(set(masses)))
- masses.sort()
- densities = densities.reshape((len(redshifts), len(masses)))
-
- logger.info("Loaded PS data from file: %s" % filepath)
- logger.info("Number of redshift bins: %d" % len(redshifts))
- logger.info("Number of mass bins: %d" % len(masses))
- self.redshifts = redshifts
- self.masses = masses
- self.densities = densities
-
@staticmethod
def delta(x, logeven=False):
"""
Calculate the delta values for each element of a vector,
assuming they are evenly or log-evenly distributed,
- with extrapolating.
+ by extrapolating.
"""
x = np.asarray(x)
if logeven:
- x = np.log(x)
- step = x[1] - x[0]
- x1 = np.concatenate([[x[0]-step], x[:-1]])
- x2 = np.concatenate([x[1:], [x[-1]+step]])
- dx = (x2 - x1) * 0.5
- if logeven:
- dx = np.exp(dx)
+ ratio = x[1] / x[0]
+ x_left = x[0] / ratio
+ x_right = x[-1] * ratio
+ else:
+ step = x[1] - x[0]
+ x_left = x[0] - step
+ x_right = x[-1] + step
+ x2 = np.concatenate([[x_left], x, [x_right]])
+ dx = (x2[2:] - x2[:-2]) / 2
return dx
@property
def number_grid(self):
"""
- Calculate the number distribution w.r.t. redshift, mass, and
- unit coverage [sr] from the density distribution.
+ The halo number per unit solid angle [sr] distribution w.r.t.
+ mass and redshift.
+ Unit: [/sr]
"""
- dz = self.delta(self.redshifts)
- dM = self.delta(self.masses)
- dMgrip, dzgrip = np.meshgrid(dM, dz)
- Mgrip, zgrip = np.meshgrid(self.masses, self.redshifts)
- dVcgrip = COSMO.dVc(zgrip) # [Mpc^3/sr]
- numgrid = self.densities * dVcgrip * dzgrip * dMgrip
- return numgrid
+ if not hasattr(self, "_number_grid"):
+ dz = self.delta(self.z)
+ dM = self.delta(self.mass, logeven=True)
+ dlnM = dM / self.mass
+ dlnMgrid, dzgrid = np.meshgrid(dlnM, dz)
+ __, zgrid = np.meshgrid(self.mass, self.z)
+ dVcgrid = COSMO.dVc(zgrid) # [Mpc^3/sr]
+ self._number_grid = self.dndlnm * dlnMgrid * (dVcgrid*dzgrid)
+
+ return self._number_grid
def calc_cluster_counts(self, coverage):
"""
Calculate the total number of clusters (>= minimum mass) within
- the FoV coverage according to the number density distribution
- (e.g., predicted by the Press-Schechter mass function)
+ the FoV coverage according to the halo number density distribution.
Parameters
----------
@@ -145,10 +212,10 @@ class PSFormalism:
----------
counts
"""
- logger.info("Determine the total number of clusters within "
+ logger.info("Calculating the total number of clusters within "
"sky patch of coverage %.1f [deg^2]" % coverage)
coverage *= AUC.deg2rad**2 # [deg^2] -> [rad^2] = [sr]
- midx = (self.masses >= self.Mmin_halo)
+ midx = (self.mass >= self.Mmin_halo)
numgrid = self.number_grid
counts = np.sum(numgrid[:, midx]) * coverage
counts *= self.boost # number boost factor
@@ -159,7 +226,7 @@ class PSFormalism:
def sample_z_m(self, counts=None):
"""
Randomly generate the requested number of pairs of (z, M) following
- the specified number distribution.
+ the halo number distribution.
Parameters
----------