aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorAaron LI <aly@aaronly.me>2017-10-04 23:36:10 +0800
committerAaron LI <aly@aaronly.me>2017-10-04 23:36:10 +0800
commit86b3eefa1efe237aacbfb43d354f95c1171f5dde (patch)
tree7152e2be834ba14b0af2ef4a30dfb306d3d6a371
parent0fe535c2ccbe408c37a6c54629c507c2788485b9 (diff)
downloadfg21sim-86b3eefa1efe237aacbfb43d354f95c1171f5dde.tar.bz2
clusters/psformalism.py: Rewrite sample_z_m() method
Also update clusters/main.py to calculate the halo mass distributions before sampling the mass and redshifts for clusters.
-rw-r--r--fg21sim/extragalactic/clusters/main.py2
-rw-r--r--fg21sim/extragalactic/clusters/psformalism.py90
2 files changed, 54 insertions, 38 deletions
diff --git a/fg21sim/extragalactic/clusters/main.py b/fg21sim/extragalactic/clusters/main.py
index 3d71593..a5c2497 100644
--- a/fg21sim/extragalactic/clusters/main.py
+++ b/fg21sim/extragalactic/clusters/main.py
@@ -104,6 +104,8 @@ class GalaxyClusters:
"""
logger.info("Simulating the clusters (z, mass) catalog ...")
psform = PSFormalism(configs=self.configs)
+ psform.calc_dndlnm()
+ psform.write()
counts = psform.calc_cluster_counts(coverage=self.sky.area)
self.catalog, self.catalog_comment = psform.sample_z_m(counts)
logger.info("Simulated cluster catalog of counts %d." % counts)
diff --git a/fg21sim/extragalactic/clusters/psformalism.py b/fg21sim/extragalactic/clusters/psformalism.py
index 899d94e..da2796d 100644
--- a/fg21sim/extragalactic/clusters/psformalism.py
+++ b/fg21sim/extragalactic/clusters/psformalism.py
@@ -18,7 +18,6 @@ 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
@@ -228,6 +227,19 @@ class PSFormalism:
Randomly generate the requested number of pairs of (z, M) following
the halo number distribution.
+ NOTE
+ ----
+ First derive the cluster (M>=Mmin) number distribution w.r.t.
+ redshifts, from which the redshift for each cluster is sampled
+ using the acceptance-rejection algorithm. Then for each cluster
+ at redshift z, the corresponding halo mass distribution is used
+ to generate the cluster mass using the same algorithm.
+
+ NOTE
+ ----
+ Sampling masses in logarithmic scale improve the speed very
+ significantly (~30x)!
+
Parameters
----------
counts : int, optional
@@ -250,52 +262,54 @@ class PSFormalism:
counts = self.counts
logger.info("Sampling (z, mass) pairs for %d clusters ..." % counts)
- redshifts = self.redshifts
- masses = self.masses
- zmin = redshifts.min()
- zmax = redshifts.max()
- Mmax = masses.max()
- midx = (masses >= self.Mmin_halo)
- numgrid = self.number_grid
- numgrid2 = numgrid[:, midx]
- NM = numgrid2.max()
+ z = self.z
+ zmin = z.min()
+ zmax = z.max()
+ log10mass = np.log10(self.mass)
+ log10Mmin = np.log10(self.Mmin_halo)
+ log10Mmax = log10mass.max()
+ midx = (log10mass >= log10Mmin)
+ log10mass = log10mass[midx]
+ Ngrid = self.number_grid[:, midx]
+
+ logger.info("Sampling redshifts ...")
z_list = []
- M_list = []
+ zi_list = []
+ Nzdist = Ngrid.sum(axis=1)
+ NMax = Nzdist.max()
i = 0
while i < counts:
- z = random.uniform(zmin, zmax)
- M = random.uniform(self.Mmin_halo, Mmax)
+ zc = random.uniform(zmin, zmax)
+ zi = (z < zc).sum()
+ Nzc = Nzdist[zi]
r = random.random()
- zi1 = (self.redshifts < z).sum()
- zi2 = zi1 - 1
- if zi2 < 0:
- zi2 += 1
- zi1 += 1
- Mi1 = (self.masses < M).sum()
- Mi2 = Mi1 - 1
- if Mi2 < 0:
- Mi2 += 1
- Mi1 += 1
- N = bilinear(
- z, np.log(M),
- p11=(redshifts[zi1], np.log(masses[Mi1]), numgrid[zi1, Mi1]),
- p12=(redshifts[zi1], np.log(masses[Mi2]), numgrid[zi1, Mi2]),
- p21=(redshifts[zi2], np.log(masses[Mi1]), numgrid[zi2, Mi1]),
- p22=(redshifts[zi2], np.log(masses[Mi2]), numgrid[zi2, Mi2]))
- if r < N/NM:
- z_list.append(z)
- M_list.append(M)
+ if r < Nzc/NMax:
+ z_list.append(zc)
+ zi_list.append(zi)
i += 1
- if i % 100 == 0:
- logger.info("[%d/%d] %.1f%% ..." %
- (i, counts, 100*i/counts))
- logger.info("Sampled %d pairs of (z, mass) for each cluster" % counts)
- df = pd.DataFrame(np.column_stack([z_list, M_list]),
+ logger.info("Sampling masses ...")
+ mass_list = []
+ NMax_list = Ngrid.max(axis=1)
+ i = 0
+ while i < counts:
+ zi = zi_list[i]
+ NMax = NMax_list[zi]
+ Nmassdist = Ngrid[zi, :]
+ log10Mc = random.uniform(log10Mmin, log10Mmax)
+ Mi = (log10mass < log10Mc).sum()
+ NMc = Nmassdist[Mi]
+ r = random.random()
+ if r < NMc/NMax:
+ mass_list.append(10**log10Mc)
+ i += 1
+
+ logger.info("Sampled %d pairs of (z, mass) for each cluster" % counts)
+ df = pd.DataFrame(np.column_stack([z_list, mass_list]),
columns=["z", "mass"])
df["mass"] /= COSMO.darkmatter_fraction
comment = [
- "dndM data: %s" % self.datafile,
+ "halo mass function model: %s" % self.hmf_model,
"cluster minimum mass: %.2e [Msun]" % self.Mmin_cluster,
"dark matter fraction: %.2f" % COSMO.darkmatter_fraction,
"cluster counts: %d" % counts,