From 3f623443558adedd2faf53823bd1270ba21f2d50 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 20 Sep 2024 15:36:12 +0000 Subject: [PATCH] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- abacusnbody/hod/GRAND_HOD.py | 887 +++++++++++++++++++++------------- abacusnbody/hod/abacus_hod.py | 137 ++++-- 2 files changed, 647 insertions(+), 377 deletions(-) diff --git a/abacusnbody/hod/GRAND_HOD.py b/abacusnbody/hod/GRAND_HOD.py index 0723944..9c9afbf 100644 --- a/abacusnbody/hod/GRAND_HOD.py +++ b/abacusnbody/hod/GRAND_HOD.py @@ -19,119 +19,167 @@ G = 4.302e-6 # in kpc/Msol (km.s)^2 - @njit(fastmath=True) def n_cen_CSMF(M_h, Mstar_low, Mstar_up, M_1, M_0, gamma1, gamma2, sigma_c): """ Standard Cacciato et al. (2008) centrals HOD parametrization for CSMF """ M_c_value = M_c(M_h, M_1, M_0, gamma1, gamma2) - x_low = np.log10(Mstar_low/M_c_value)/(1.41421356*sigma_c) - x_up = np.log10(Mstar_up/M_c_value)/(1.41421356*sigma_c) - - return 0.5*(math.erf(x_up) - math.erf(x_low)) + x_low = np.log10(Mstar_low / M_c_value) / (1.41421356 * sigma_c) + x_up = np.log10(Mstar_up / M_c_value) / (1.41421356 * sigma_c) + + return 0.5 * (math.erf(x_up) - math.erf(x_low)) + @njit(fastmath=True) def CSMF_centrals(M_h, Mstar, M_1, M_0, gamma1, gamma2, sigma_c): """ Eq. (34) from Cacciato et al. (2008) """ - + M_c_value = M_c(M_h, M_1, M_0, gamma1, gamma2) - - return 1/(1.41421356*np.sqrt(np.pi)*np.log(10)*sigma_c*Mstar)*np.exp(-(np.log10(Mstar)-np.log10(M_c_value))**2/(2*sigma_c**2)) + + return ( + 1 + / (1.41421356 * np.sqrt(np.pi) * np.log(10) * sigma_c * Mstar) + * np.exp(-((np.log10(Mstar) - np.log10(M_c_value)) ** 2) / (2 * sigma_c**2)) + ) + @njit(fastmath=True) def M_c(M_h, M_1, M_0, gamma1, gamma2): """ Eq. (37) from Cacciato et al. (2008) """ - return M_0 * (M_h/M_1)**gamma1/(1+M_h/M_1)**(gamma1-gamma2) - + return M_0 * (M_h / M_1) ** gamma1 / (1 + M_h / M_1) ** (gamma1 - gamma2) @njit(fastmath=True) -def get_random_cen_stellarmass(M_h, Mstar_low, Mstar_up, M_1, M_0, gamma1, gamma2, sigma_c): - +def get_random_cen_stellarmass( + M_h, Mstar_low, Mstar_up, M_1, M_0, gamma1, gamma2, sigma_c +): nbins = 1000 - stellarmass = np.logspace(np.log10(Mstar_low),np.log10(Mstar_up),nbins) - stellarmass_centers = stellarmass[1:]/2+stellarmass[:-1]/2 - delta_stellar_masses = stellarmass[1:]-stellarmass[:-1] + stellarmass = np.logspace(np.log10(Mstar_low), np.log10(Mstar_up), nbins) + stellarmass_centers = stellarmass[1:] / 2 + stellarmass[:-1] / 2 + delta_stellar_masses = stellarmass[1:] - stellarmass[:-1] + + CSMF_cen = CSMF_centrals( + M_h, stellarmass_centers, M_1, M_0, gamma1, gamma2, sigma_c + ) - CSMF_cen = CSMF_centrals(M_h, stellarmass_centers, M_1, M_0, gamma1, gamma2, sigma_c) - - cdf = [] + cdf = [] cdf_current_value = 0 for k in range(len(delta_stellar_masses)): - cdf_current_value = cdf_current_value + CSMF_cen[k]*delta_stellar_masses[k] + cdf_current_value = cdf_current_value + CSMF_cen[k] * delta_stellar_masses[k] cdf.append(cdf_current_value) - cdf=np.array(cdf) - - cdf=cdf/cdf[-1] - - random_rv = np.random.uniform(cdf.min(),cdf.max()) + cdf = np.array(cdf) + + cdf = cdf / cdf[-1] + + random_rv = np.random.uniform(cdf.min(), cdf.max()) bin_clostest = (np.abs(cdf - random_rv)).argmin() - - return np.random.uniform(stellarmass[bin_clostest],stellarmass[bin_clostest+1]) - + + return np.random.uniform(stellarmass[bin_clostest], stellarmass[bin_clostest + 1]) - @njit(fastmath=True) -def get_random_cen_stellarmass_linearinterpolation(M_h, Mstar_low, Mstar_up, M_1, M_0, gamma1, gamma2, sigma_c): - +def get_random_cen_stellarmass_linearinterpolation( + M_h, Mstar_low, Mstar_up, M_1, M_0, gamma1, gamma2, sigma_c +): nbins = 1000 - stellarmass = np.logspace(np.log10(Mstar_low),np.log10(Mstar_up),nbins) + stellarmass = np.logspace(np.log10(Mstar_low), np.log10(Mstar_up), nbins) CSMF_cen = CSMF_centrals(M_h, stellarmass, M_1, M_0, gamma1, gamma2, sigma_c) - - delta=stellarmass[1:]-stellarmass[:-1] - - cdf = [] + + delta = stellarmass[1:] - stellarmass[:-1] + + cdf = [] cdf_current_value = 0 for i in range(len(delta)): - cdf_current_value = cdf_current_value + CSMF_cen[i]*delta[i] + cdf_current_value = cdf_current_value + CSMF_cen[i] * delta[i] cdf.append(cdf_current_value) - cdf=np.array(cdf) - - cdf=cdf/cdf[-1] - - random_rv = np.random.uniform(cdf.min(),cdf.max()) - bin = np.where(cdf>random_rv)[0][0] - - m = (stellarmass[bin]-stellarmass[bin-1])/(cdf[bin]-cdf[bin-1]) - return m * (random_rv - cdf[bin-1]) + stellarmass[bin-1] + cdf = np.array(cdf) + + cdf = cdf / cdf[-1] + + random_rv = np.random.uniform(cdf.min(), cdf.max()) + bin = np.where(cdf > random_rv)[0][0] + + m = (stellarmass[bin] - stellarmass[bin - 1]) / (cdf[bin] - cdf[bin - 1]) + return m * (random_rv - cdf[bin - 1]) + stellarmass[bin - 1] @njit(fastmath=True) -def n_sat_CSMF(M_h, Mstar_low, Mstar_up, M_1, M_0, gamma1, gamma2, sigma_c, a1, a2, M2, b0, b1, b2, delta1, delta2): +def n_sat_CSMF( + M_h, + Mstar_low, + Mstar_up, + M_1, + M_0, + gamma1, + gamma2, + sigma_c, + a1, + a2, + M2, + b0, + b1, + b2, + delta1, + delta2, +): """ Standard Cacciato et al. (2008) satellite HOD parametrization for CSMF """ nbins = 1000 - stellarmass = np.logspace(np.log10(Mstar_low),np.log10(Mstar_up),nbins) - - CSMF_sat = CSMF_satelites(M_h, stellarmass, M_1, M_0, gamma1, gamma2, a1, a2, M2, b0, b1, b2, delta1, delta2) + stellarmass = np.logspace(np.log10(Mstar_low), np.log10(Mstar_up), nbins) + + CSMF_sat = CSMF_satelites( + M_h, + stellarmass, + M_1, + M_0, + gamma1, + gamma2, + a1, + a2, + M2, + b0, + b1, + b2, + delta1, + delta2, + ) nsat = 0 - for i in range(nbins-1): - nsat += (CSMF_sat[i+1]-CSMF_sat[i])*(stellarmass[i+1]-stellarmass[i])/2 + (stellarmass[i+1]-stellarmass[i])*CSMF_sat[i] - - return nsat#*ncen + for i in range(nbins - 1): + nsat += (CSMF_sat[i + 1] - CSMF_sat[i]) * ( + stellarmass[i + 1] - stellarmass[i] + ) / 2 + (stellarmass[i + 1] - stellarmass[i]) * CSMF_sat[i] + + return nsat # *ncen @njit(fastmath=True) -def CSMF_satelites(M_h, Mstar, M_1, M_0, gamma1, gamma2, a1, a2, M2, b0, b1, b2, delta1, delta2): +def CSMF_satelites( + M_h, Mstar, M_1, M_0, gamma1, gamma2, a1, a2, M2, b0, b1, b2, delta1, delta2 +): """ Eq. (36) from Cacciato et al. (2008) """ M_s_value = M_s(M_h, M_1, M_0, gamma1, gamma2) alpha_s_value = alpha_s(M_h, a1, a2, M2) phi_s_value = phi_s(M_h, b0, b1, b2) - + delta = 10 ** (delta1 + delta2 * (np.log10(M_h) - 12)) - - return phi_s_value/M_s_value * (Mstar/M_s_value)**alpha_s_value * np.exp(-delta*(Mstar/M_s_value)**2) + + return ( + phi_s_value + / M_s_value + * (Mstar / M_s_value) ** alpha_s_value + * np.exp(-delta * (Mstar / M_s_value) ** 2) + ) + @njit(fastmath=True) def M_s(M_h, M_1, M_0, gamma1, gamma2): @@ -140,73 +188,134 @@ def M_s(M_h, M_1, M_0, gamma1, gamma2): """ return 0.562 * M_c(M_h, M_1, M_0, gamma1, gamma2) + @njit(fastmath=True) def alpha_s(M_h, a1, a2, M2): """ Eq. (39) from Cacciato et al. (2008) """ - return -2.0 + a1 * (1-2/np.pi*np.arctan(a2*np.log10(M_h/M2))) + return -2.0 + a1 * (1 - 2 / np.pi * np.arctan(a2 * np.log10(M_h / M2))) + @njit(fastmath=True) def phi_s(M_h, b0, b1, b2): """ Eq. (40) from Cacciato et al. (2008) """ - M12 = M_h/1e12 - log_phi_s = b0 + b1 * np.log10(M12) + b2 * np.log10(M12)**2 + M12 = M_h / 1e12 + log_phi_s = b0 + b1 * np.log10(M12) + b2 * np.log10(M12) ** 2 return 10**log_phi_s @njit(fastmath=True) -def get_random_sat_stellarmass(M_h, Mstar_low, Mstar_up, M_1, M_0, gamma1, gamma2, a1, a2, M2, b0, b1, b2, delta1, delta2): - +def get_random_sat_stellarmass( + M_h, + Mstar_low, + Mstar_up, + M_1, + M_0, + gamma1, + gamma2, + a1, + a2, + M2, + b0, + b1, + b2, + delta1, + delta2, +): nbins = 1000 - stellarmass = np.logspace(np.log10(Mstar_low),np.log10(Mstar_up),nbins) - stellarmass_centers = stellarmass[1:]/2+stellarmass[:-1]/2 - delta_stellar_masses = stellarmass[1:]-stellarmass[:-1] - - CSMF_sat = CSMF_satelites(M_h, stellarmass_centers, M_1, M_0, gamma1, gamma2, a1, a2, M2, b0, b1, b2, delta1, delta2) - - - cdf = [] + stellarmass = np.logspace(np.log10(Mstar_low), np.log10(Mstar_up), nbins) + stellarmass_centers = stellarmass[1:] / 2 + stellarmass[:-1] / 2 + delta_stellar_masses = stellarmass[1:] - stellarmass[:-1] + + CSMF_sat = CSMF_satelites( + M_h, + stellarmass_centers, + M_1, + M_0, + gamma1, + gamma2, + a1, + a2, + M2, + b0, + b1, + b2, + delta1, + delta2, + ) + + cdf = [] cdf_current_value = 0 for k in range(len(delta_stellar_masses)): - cdf_current_value = cdf_current_value + CSMF_sat[k]*delta_stellar_masses[k] + cdf_current_value = cdf_current_value + CSMF_sat[k] * delta_stellar_masses[k] cdf.append(cdf_current_value) - cdf=np.array(cdf) - - cdf=cdf/cdf[-1] - - random_rv = np.random.uniform(cdf.min(),cdf.max()) + cdf = np.array(cdf) + + cdf = cdf / cdf[-1] + + random_rv = np.random.uniform(cdf.min(), cdf.max()) bin_clostest = (np.abs(cdf - random_rv)).argmin() - - return np.random.uniform(stellarmass[bin_clostest],stellarmass[bin_clostest+1]) - + + return np.random.uniform(stellarmass[bin_clostest], stellarmass[bin_clostest + 1]) + @njit(fastmath=True) -def get_random_sat_stellarmass_linearinterpolation(M_h, Mstar_low, Mstar_up, M_1, M_0, gamma1, gamma2, a1, a2, M2, b0, b1, b2, delta1, delta2): - +def get_random_sat_stellarmass_linearinterpolation( + M_h, + Mstar_low, + Mstar_up, + M_1, + M_0, + gamma1, + gamma2, + a1, + a2, + M2, + b0, + b1, + b2, + delta1, + delta2, +): nbins = 100 - stellarmass = np.logspace(np.log10(Mstar_low),np.log10(Mstar_up),nbins) - - CSMF_sat = CSMF_satelites(M_h, stellarmass, M_1, M_0, gamma1, gamma2, a1, a2, M2, b0, b1, b2, delta1, delta2) - - delta=stellarmass[1:]-stellarmass[:-1] - - cdf = [] + stellarmass = np.logspace(np.log10(Mstar_low), np.log10(Mstar_up), nbins) + + CSMF_sat = CSMF_satelites( + M_h, + stellarmass, + M_1, + M_0, + gamma1, + gamma2, + a1, + a2, + M2, + b0, + b1, + b2, + delta1, + delta2, + ) + + delta = stellarmass[1:] - stellarmass[:-1] + + cdf = [] cdf_current_value = 0 for i in range(len(delta)): - cdf_current_value = cdf_current_value + CSMF_sat[i]*delta[i] + cdf_current_value = cdf_current_value + CSMF_sat[i] * delta[i] cdf.append(cdf_current_value) - cdf=np.array(cdf) - - cdf=cdf/cdf[-1] - - random_rv = np.random.uniform(cdf.min(),cdf.max()) - bin = np.where(cdf>random_rv)[0][0] - - m = (stellarmass[bin]-stellarmass[bin-1])/(cdf[bin]-cdf[bin-1]) - return m * (random_rv - cdf[bin-1]) + stellarmass[bin-1] + cdf = np.array(cdf) + + cdf = cdf / cdf[-1] + + random_rv = np.random.uniform(cdf.min(), cdf.max()) + bin = np.where(cdf > random_rv)[0][0] + + m = (stellarmass[bin] - stellarmass[bin - 1]) / (cdf[bin] - cdf[bin - 1]) + return m * (random_rv - cdf[bin - 1]) + stellarmass[bin - 1] @njit(fastmath=True) @@ -391,21 +500,20 @@ def gen_cent( ) if want_CSMF: Mstar_low_C, Mstar_up_C, M_1_C, M_0_C, gamma1_C, gamma2_C, sigma_c_C = ( - CSMF_hod_dict['Mstar_low'], - CSMF_hod_dict['Mstar_up'], - CSMF_hod_dict['M_1'], - CSMF_hod_dict['M_0'], - CSMF_hod_dict['gamma_1'], - CSMF_hod_dict['gamma_2'], - CSMF_hod_dict['sigma_c'] + CSMF_hod_dict['Mstar_low'], + CSMF_hod_dict['Mstar_up'], + CSMF_hod_dict['M_1'], + CSMF_hod_dict['M_0'], + CSMF_hod_dict['gamma_1'], + CSMF_hod_dict['gamma_2'], + CSMF_hod_dict['sigma_c'], ) ic_C, alpha_c_C, Ac_C, Bc_C = ( - CSMF_hod_dict['ic'], - CSMF_hod_dict['alpha_c'], - CSMF_hod_dict['Acent'], - CSMF_hod_dict['Bcent'] + CSMF_hod_dict['ic'], + CSMF_hod_dict['alpha_c'], + CSMF_hod_dict['Acent'], + CSMF_hod_dict['Bcent'], ) - H = len(mass) @@ -446,13 +554,22 @@ def gen_cent( QSO_marker += ( N_cen_QSO(mass[i], logM_cut_Q_temp, sigma_Q) * ic_Q * multis[i] ) - + CSMF_marker = QSO_marker if want_CSMF: - M_1_C_temp = 10**(np.log10(M_1_C) + Ac_C * deltac[i] + Bc_C * fenv[i]) - - ncen = n_cen_CSMF(mass[i], Mstar_low_C, Mstar_up_C, M_1_C_temp, M_0_C, gamma1_C, gamma2_C, sigma_c_C) - + M_1_C_temp = 10 ** (np.log10(M_1_C) + Ac_C * deltac[i] + Bc_C * fenv[i]) + + ncen = n_cen_CSMF( + mass[i], + Mstar_low_C, + Mstar_up_C, + M_1_C_temp, + M_0_C, + gamma1_C, + gamma2_C, + sigma_c_C, + ) + CSMF_marker += ncen * ic_C * multis[i] if randoms[i] <= LRG_marker: @@ -465,7 +582,7 @@ def gen_cent( Nout[tid, 2, 0] += 1 # counting keep[i] = 3 elif randoms[i] <= CSMF_marker: - Nout[tid, 3, 0] += 1 # counting + Nout[tid, 3, 0] += 1 # counting keep[i] = 4 else: keep[i] = 0 @@ -510,18 +627,18 @@ def gen_cent( qso_vz = np.empty(N_qso, dtype=mass.dtype) qso_mass = np.empty(N_qso, dtype=mass.dtype) qso_id = np.empty(N_qso, dtype=ids.dtype) - + # galaxy arrays N_CSMF = gstart[-1, 3] - CSMF_x = np.empty(N_CSMF, dtype = mass.dtype) - CSMF_y = np.empty(N_CSMF, dtype = mass.dtype) - CSMF_z = np.empty(N_CSMF, dtype = mass.dtype) - CSMF_vx = np.empty(N_CSMF, dtype = mass.dtype) - CSMF_vy = np.empty(N_CSMF, dtype = mass.dtype) - CSMF_vz = np.empty(N_CSMF, dtype = mass.dtype) - CSMF_mass = np.empty(N_CSMF, dtype = mass.dtype) - CSMF_stellarmass = np.empty(N_CSMF, dtype = mass.dtype) - CSMF_id = np.empty(N_CSMF, dtype = ids.dtype) + CSMF_x = np.empty(N_CSMF, dtype=mass.dtype) + CSMF_y = np.empty(N_CSMF, dtype=mass.dtype) + CSMF_z = np.empty(N_CSMF, dtype=mass.dtype) + CSMF_vx = np.empty(N_CSMF, dtype=mass.dtype) + CSMF_vy = np.empty(N_CSMF, dtype=mass.dtype) + CSMF_vz = np.empty(N_CSMF, dtype=mass.dtype) + CSMF_mass = np.empty(N_CSMF, dtype=mass.dtype) + CSMF_stellarmass = np.empty(N_CSMF, dtype=mass.dtype) + CSMF_id = np.empty(N_CSMF, dtype=ids.dtype) # fill in the galaxy arrays for tid in numba.prange(Nthread): @@ -613,34 +730,43 @@ def gen_cent( j3 += 1 elif keep[i] == 4: # loop thru three directions to assign galaxy velocities and positions - CSMF_x[j4] = pos[i,0] - CSMF_vx[j4] = vel[i,0] + alpha_c_C * vdev[i,0] # velocity bias - CSMF_y[j4] = pos[i,1] - CSMF_vy[j4] = vel[i,1] + alpha_c_C * vdev[i,1] # velocity bias - CSMF_z[j4] = pos[i,2] - CSMF_vz[j4] = vel[i,2] + alpha_c_C * vdev[i,2] # velocity bias - + CSMF_x[j4] = pos[i, 0] + CSMF_vx[j4] = vel[i, 0] + alpha_c_C * vdev[i, 0] # velocity bias + CSMF_y[j4] = pos[i, 1] + CSMF_vy[j4] = vel[i, 1] + alpha_c_C * vdev[i, 1] # velocity bias + CSMF_z[j4] = pos[i, 2] + CSMF_vz[j4] = vel[i, 2] + alpha_c_C * vdev[i, 2] # velocity bias + # rsd only applies to the z direction if rsd and origin is not None: nx = CSMF_x[j4] - origin[0] ny = CSMF_y[j4] - origin[1] nz = CSMF_z[j4] - origin[2] - inv_norm = 1./np.sqrt(nx*nx + ny*ny + nz*nz) + inv_norm = 1.0 / np.sqrt(nx * nx + ny * ny + nz * nz) nx *= inv_norm ny *= inv_norm nz *= inv_norm - proj = inv_velz2kms*(CSMF_vx[j4]*nx+CSMF_vy[j4]*ny+CSMF_vz[j4]*nz) - CSMF_x[j4] = CSMF_x[j4]+proj*nx - CSMF_y[j4] = CSMF_y[j4]+proj*ny - CSMF_z[j4] = CSMF_z[j4]+proj*nz + proj = inv_velz2kms * ( + CSMF_vx[j4] * nx + CSMF_vy[j4] * ny + CSMF_vz[j4] * nz + ) + CSMF_x[j4] = CSMF_x[j4] + proj * nx + CSMF_y[j4] = CSMF_y[j4] + proj * ny + CSMF_z[j4] = CSMF_z[j4] + proj * nz elif rsd: - CSMF_z[j4] = wrap(pos[i,2] + CSMF_vz[j4] * inv_velz2kms, lbox) - + CSMF_z[j4] = wrap(pos[i, 2] + CSMF_vz[j4] * inv_velz2kms, lbox) + CSMF_mass[j4] = mass[i] - M_1_C_temp = 10**(np.log10(M_1_C) + Ac_C * deltac[i] + Bc_C * fenv[i]) + M_1_C_temp = 10 ** (np.log10(M_1_C) + Ac_C * deltac[i] + Bc_C * fenv[i]) CSMF_stellarmass[j4] = get_random_cen_stellarmass_linearinterpolation( - mass[i], Mstar_low_C, Mstar_up_C, - M_1_C_temp, M_0_C, gamma1_C, gamma2_C, sigma_c_C) + mass[i], + Mstar_low_C, + Mstar_up_C, + M_1_C_temp, + M_0_C, + gamma1_C, + gamma2_C, + sigma_c_C, + ) CSMF_id[j4] = ids[i] j4 += 1 # assert j == gstart[tid + 1] @@ -648,7 +774,7 @@ def gen_cent( LRG_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) ELG_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) QSO_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) - CSMF_dict = Dict.empty(key_type = types.unicode_type, value_type = float_array) + CSMF_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) ID_dict = Dict.empty(key_type=types.unicode_type, value_type=int_array) LRG_dict['x'] = lrg_x LRG_dict['y'] = lrg_y @@ -676,7 +802,7 @@ def gen_cent( QSO_dict['vz'] = qso_vz QSO_dict['mass'] = qso_mass ID_dict['QSO'] = qso_id - + CSMF_dict['x'] = CSMF_x CSMF_dict['y'] = CSMF_y CSMF_dict['z'] = CSMF_z @@ -897,31 +1023,47 @@ def gen_sats_nfw( QSO_hod_dict['ic'], ) f_sigv_Q = QSO_hod_dict['f_sigv'] - + if want_CSMF: - Mstar_low_C, Mstar_up_C, M_1_C, M_0_C, gamma1_C, gamma2_C, sigma_c_C, a1_C, a2_C, M2_C, b0_C, b1_C, b2_C, delta1_C, delta2_C = ( - CSMF_hod_dict['Mstar_low'], - CSMF_hod_dict['Mstar_up'], - CSMF_hod_dict['M_1'], - CSMF_hod_dict['M_0'], - CSMF_hod_dict['gamma_1'], - CSMF_hod_dict['gamma_2'], - CSMF_hod_dict['sigma_c'], - CSMF_hod_dict['a_1'], - CSMF_hod_dict['a_2'], - CSMF_hod_dict['M_2'], - CSMF_hod_dict['b_0'], - CSMF_hod_dict['b_1'], - CSMF_hod_dict['b_2'], - CSMF_hod_dict['delta_1'], - CSMF_hod_dict['delta_2'] + ( + Mstar_low_C, + Mstar_up_C, + M_1_C, + M_0_C, + gamma1_C, + gamma2_C, + sigma_c_C, + a1_C, + a2_C, + M2_C, + b0_C, + b1_C, + b2_C, + delta1_C, + delta2_C, + ) = ( + CSMF_hod_dict['Mstar_low'], + CSMF_hod_dict['Mstar_up'], + CSMF_hod_dict['M_1'], + CSMF_hod_dict['M_0'], + CSMF_hod_dict['gamma_1'], + CSMF_hod_dict['gamma_2'], + CSMF_hod_dict['sigma_c'], + CSMF_hod_dict['a_1'], + CSMF_hod_dict['a_2'], + CSMF_hod_dict['M_2'], + CSMF_hod_dict['b_0'], + CSMF_hod_dict['b_1'], + CSMF_hod_dict['b_2'], + CSMF_hod_dict['delta_1'], + CSMF_hod_dict['delta_2'], ) Ac_C, As_C, Bc_C, Bs_C, ic_C = ( - CSMF_hod_dict['Acent'], - CSMF_hod_dict['Asat'], + CSMF_hod_dict['Acent'], + CSMF_hod_dict['Asat'], CSMF_hod_dict['Bcent'], - CSMF_hod_dict['Bsat'], - CSMF_hod_dict['ic'] + CSMF_hod_dict['Bsat'], + CSMF_hod_dict['ic'], ) f_sigv_C = CSMF_hod_dict['f_sigv'] @@ -932,8 +1074,8 @@ def gen_sats_nfw( num_sats_L = np.zeros(len(hid), dtype=np.int64) num_sats_E = np.zeros(len(hid), dtype=np.int64) num_sats_Q = np.zeros(len(hid), dtype=np.int64) - num_sats_C = np.zeros(len(hid), dtype = np.int64) - stellarmass_C = np.zeros(len(hid), dtype = np.int64) + num_sats_C = np.zeros(len(hid), dtype=np.int64) + stellarmass_C = np.zeros(len(hid), dtype=np.int64) hstart = np.rint(np.linspace(0, len(hid), Nthread + 1)).astype( np.int64 ) # starting index of each thread @@ -1009,32 +1151,35 @@ def gen_sats_nfw( * ic_Q ) num_sats_Q[i] = np.random.poisson(base_p_Q) - + if want_CSMF: - M_1_C_temp = 10**(np.log10(M_1_C) + Ac_C * hdeltac[i] + Bc_C * hfenv[i]) + M_1_C_temp = 10 ** ( + np.log10(M_1_C) + Ac_C * hdeltac[i] + Bc_C * hfenv[i] + ) a1_C_temp = a1_C + As_C * hdeltac[i] + Bs_C * hfenv[i] base_p_C = ( n_sat_CSMF( - hmass[i], - Mstar_low_C, - Mstar_up_C, - M_1_C_temp, - M_0_C, - gamma1_C, - gamma2_C, - sigma_c_C, - a1_C_temp, - a2_C, - M2_C, - b0_C, - b1_C, - b2_C, - delta1_C, - delta2_C) + hmass[i], + Mstar_low_C, + Mstar_up_C, + M_1_C_temp, + M_0_C, + gamma1_C, + gamma2_C, + sigma_c_C, + a1_C_temp, + a2_C, + M2_C, + b0_C, + b1_C, + b2_C, + delta1_C, + delta2_C, + ) * ic_C ) - - num_sats_C[i] = np.random.poisson(base_p_C) + + num_sats_C[i] = np.random.poisson(base_p_C) # generate rdpos rd_pos_L = getPointsOnSphere(np.sum(num_sats_L), Nthread) @@ -1115,33 +1260,32 @@ def gen_sats_nfw( nfw_rescale, ) ) - + h_id_C, x_sat_C, y_sat_C, z_sat_C, vx_sat_C, vy_sat_C, vz_sat_C, M_C = ( compute_fast_NFW( - NFW_draw, - hid, - hpos[:, 0], - hpos[:, 1], - hpos[:, 2], - hvel[:, 0], - hvel[:, 1], + NFW_draw, + hid, + hpos[:, 0], + hpos[:, 1], + hpos[:, 2], + hvel[:, 0], + hvel[:, 1], hvel[:, 2], - hvrms, - hc, - hmass, - hrvir, - rd_pos_C, - num_sats_C, - f_sigv_C, - vel_sat, + hvrms, + hc, + hmass, + hrvir, + rd_pos_C, + num_sats_C, + f_sigv_C, + vel_sat, Nthread, - exp_frac, - exp_scale, - nfw_rescale + exp_frac, + exp_scale, + nfw_rescale, ) ) - - + # do rsd if rsd: z_sat_L = (z_sat_L + vz_sat_L * inv_velz2kms) % lbox @@ -1152,7 +1296,7 @@ def gen_sats_nfw( LRG_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) ELG_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) QSO_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) - CSMF_dict = Dict.empty(key_type = types.unicode_type, value_type = float_array) + CSMF_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) ID_dict = Dict.empty(key_type=types.unicode_type, value_type=int_array) LRG_dict['x'] = x_sat_L LRG_dict['y'] = y_sat_L @@ -1180,7 +1324,7 @@ def gen_sats_nfw( QSO_dict['vz'] = vz_sat_Q QSO_dict['mass'] = M_Q ID_dict['QSO'] = h_id_Q - + CSMF_dict['x'] = x_sat_C CSMF_dict['y'] = y_sat_C CSMF_dict['z'] = z_sat_C @@ -1194,12 +1338,29 @@ def gen_sats_nfw( hstart = np.rint(np.linspace(0, num_sats_C.sum(), Nthread + 1)) for tid in numba.prange(Nthread): for i in range(int(hstart[tid]), int(hstart[tid + 1])): - M_1_C_temp = 10**(np.log10(M_1_C) + Ac_C * hdeltac[i] + Bc_C * hfenv[i]) + M_1_C_temp = 10 ** ( + np.log10(M_1_C) + Ac_C * hdeltac[i] + Bc_C * hfenv[i] + ) a1_C_temp = a1_C + As_C * hdeltac[i] + Bs_C * hfenv[i] stellarmass_C[i] = get_random_sat_stellarmass_linearinterpolation( - M_C[i], Mstar_low_C, Mstar_up_C, M_1_C_temp, M_0_C, gamma1_C, gamma2_C, - a1_C_temp, s, a2_C, M2_C, b0_C, b1_C, b2_C, delta1_C, delta2_C) - + M_C[i], + Mstar_low_C, + Mstar_up_C, + M_1_C_temp, + M_0_C, + gamma1_C, + gamma2_C, + a1_C_temp, + s, + a2_C, + M2_C, + b0_C, + b1_C, + b2_C, + delta1_C, + delta2_C, + ) + CSMF_dict['stellarmass'] = stellarmass_C ID_dict['CSMF'] = h_id_C @@ -1328,37 +1489,53 @@ def gen_sats( QSO_hod_dict['Bsat'], QSO_hod_dict['ic'], ) - + if want_CSMF: - Mstar_low_C, Mstar_up_C, M_1_C, M_0_C, gamma1_C, gamma2_C, sigma_c_C, a1_C, a2_C, M2_C, b0_C, b1_C, b2_C, delta1_C, delta2_C = ( - CSMF_hod_dict['Mstar_low'], - CSMF_hod_dict['Mstar_up'], - CSMF_hod_dict['M_1'], - CSMF_hod_dict['M_0'], + ( + Mstar_low_C, + Mstar_up_C, + M_1_C, + M_0_C, + gamma1_C, + gamma2_C, + sigma_c_C, + a1_C, + a2_C, + M2_C, + b0_C, + b1_C, + b2_C, + delta1_C, + delta2_C, + ) = ( + CSMF_hod_dict['Mstar_low'], + CSMF_hod_dict['Mstar_up'], + CSMF_hod_dict['M_1'], + CSMF_hod_dict['M_0'], CSMF_hod_dict['gamma_1'], - CSMF_hod_dict['gamma_2'], - CSMF_hod_dict['sigma_c'], - CSMF_hod_dict['a_1'], - CSMF_hod_dict['a_2'], + CSMF_hod_dict['gamma_2'], + CSMF_hod_dict['sigma_c'], + CSMF_hod_dict['a_1'], + CSMF_hod_dict['a_2'], CSMF_hod_dict['M_2'], - CSMF_hod_dict['b_0'], - CSMF_hod_dict['b_1'], - CSMF_hod_dict['b_2'], - CSMF_hod_dict['delta_1'], - CSMF_hod_dict['delta_2'] + CSMF_hod_dict['b_0'], + CSMF_hod_dict['b_1'], + CSMF_hod_dict['b_2'], + CSMF_hod_dict['delta_1'], + CSMF_hod_dict['delta_2'], ) - + alpha_s_C, s_C, s_v_C, s_p_C, s_r_C, Ac_C, As_C, Bc_C, Bs_C, ic_C = ( - CSMF_hod_dict['alpha_s'], + CSMF_hod_dict['alpha_s'], CSMF_hod_dict['s'], - CSMF_hod_dict['s_v'], - CSMF_hod_dict['s_p'], - CSMF_hod_dict['s_r'], + CSMF_hod_dict['s_v'], + CSMF_hod_dict['s_p'], + CSMF_hod_dict['s_r'], CSMF_hod_dict['Acent'], - CSMF_hod_dict['Asat'], - CSMF_hod_dict['Bcent'], - CSMF_hod_dict['Bsat'], - CSMF_hod_dict['ic'] + CSMF_hod_dict['Asat'], + CSMF_hod_dict['Bcent'], + CSMF_hod_dict['Bsat'], + CSMF_hod_dict['ic'], ) H = len(hmass) # num of particles @@ -1492,35 +1669,43 @@ def gen_sats( else: exp_sat = base_p_Q QSO_marker += exp_sat - - + CSMF_marker = QSO_marker if want_CSMF: - M_1_C_temp = 10**(np.log10(M_1_C) + Ac_C * hdeltac[i] + Bc_C * hfenv[i]) + M_1_C_temp = 10 ** ( + np.log10(M_1_C) + Ac_C * hdeltac[i] + Bc_C * hfenv[i] + ) a1_C_temp = a1_C + As_C * hdeltac[i] + Bs_C * hfenv[i] base_p_C = ( n_sat_CSMF( - hmass[i], - Mstar_low_C, - Mstar_up_C, - M_1_C_temp, - M_0_C, - gamma1_C, - gamma2_C, - sigma_c_C, - a1_C_temp, - a2_C, - M2_C, - b0_C, - b1_C, - b2_C, - delta1_C, - delta2_C - ) - * weights[i] - * ic_C) + hmass[i], + Mstar_low_C, + Mstar_up_C, + M_1_C_temp, + M_0_C, + gamma1_C, + gamma2_C, + sigma_c_C, + a1_C_temp, + a2_C, + M2_C, + b0_C, + b1_C, + b2_C, + delta1_C, + delta2_C, + ) + * weights[i] + * ic_C + ) if enable_ranks: - decorator_C = 1 + s_C * ranks[i] + s_v_C * ranksv[i] + s_p_C * ranksp[i] + s_r_C * ranksr[i] + decorator_C = ( + 1 + + s_C * ranks[i] + + s_v_C * ranksv[i] + + s_p_C * ranksp[i] + + s_r_C * ranksr[i] + ) exp_sat = base_p_C * decorator_C else: exp_sat = base_p_C @@ -1536,7 +1721,7 @@ def gen_sats( Nout[tid, 2, 0] += 1 # counting keep[i] = 3 elif randoms[i] <= CSMF_marker: - Nout[tid, 3, 0] += 1 # counting + Nout[tid, 3, 0] += 1 # counting keep[i] = 4 else: keep[i] = 0 @@ -1581,18 +1766,18 @@ def gen_sats( qso_vz = np.empty(N_qso, dtype=hmass.dtype) qso_mass = np.empty(N_qso, dtype=hmass.dtype) qso_id = np.empty(N_qso, dtype=hid.dtype) - + # galaxy arrays N_CSMF = gstart[-1, 3] - CSMF_x = np.empty(N_CSMF, dtype = hmass.dtype) - CSMF_y = np.empty(N_CSMF, dtype = hmass.dtype) - CSMF_z = np.empty(N_CSMF, dtype = hmass.dtype) - CSMF_vx = np.empty(N_CSMF, dtype = hmass.dtype) - CSMF_vy = np.empty(N_CSMF, dtype = hmass.dtype) - CSMF_vz = np.empty(N_CSMF, dtype = hmass.dtype) - CSMF_mass = np.empty(N_CSMF, dtype = hmass.dtype) - CSMF_stellarmass = np.empty(N_CSMF, dtype = hmass.dtype) - CSMF_id = np.empty(N_CSMF, dtype = hid.dtype) + CSMF_x = np.empty(N_CSMF, dtype=hmass.dtype) + CSMF_y = np.empty(N_CSMF, dtype=hmass.dtype) + CSMF_z = np.empty(N_CSMF, dtype=hmass.dtype) + CSMF_vx = np.empty(N_CSMF, dtype=hmass.dtype) + CSMF_vy = np.empty(N_CSMF, dtype=hmass.dtype) + CSMF_vz = np.empty(N_CSMF, dtype=hmass.dtype) + CSMF_mass = np.empty(N_CSMF, dtype=hmass.dtype) + CSMF_stellarmass = np.empty(N_CSMF, dtype=hmass.dtype) + CSMF_id = np.empty(N_CSMF, dtype=hid.dtype) # fill in the galaxy arrays for tid in numba.prange(Nthread): @@ -1696,46 +1881,55 @@ def gen_sats( j3 += 1 elif keep[i] == 4: CSMF_x[j4] = ppos[i, 0] - CSMF_vx[j4] = hvel[i, 0] + alpha_s_C * (pvel[i, 0] - hvel[i, 0]) # velocity bias + CSMF_vx[j4] = hvel[i, 0] + alpha_s_C * ( + pvel[i, 0] - hvel[i, 0] + ) # velocity bias CSMF_y[j4] = ppos[i, 1] - CSMF_vy[j4] = hvel[i, 1] + alpha_s_C * (pvel[i, 1] - hvel[i, 1]) # velocity bias + CSMF_vy[j4] = hvel[i, 1] + alpha_s_C * ( + pvel[i, 1] - hvel[i, 1] + ) # velocity bias CSMF_z[j4] = ppos[i, 2] - CSMF_vz[j4] = hvel[i, 2] + alpha_s_C * (pvel[i, 2] - hvel[i, 2]) # velocity bias + CSMF_vz[j4] = hvel[i, 2] + alpha_s_C * ( + pvel[i, 2] - hvel[i, 2] + ) # velocity bias if rsd and origin is not None: nx = CSMF_x[j4] - origin[0] ny = CSMF_y[j4] - origin[1] nz = CSMF_z[j4] - origin[2] - inv_norm = 1./np.sqrt(nx*nx + ny*ny + nz*nz) + inv_norm = 1.0 / np.sqrt(nx * nx + ny * ny + nz * nz) nx *= inv_norm ny *= inv_norm nz *= inv_norm - proj = inv_velz2kms*(CSMF_vx[j4]*nx+CSMF_vy[j4]*ny+CSMF_vz[j4]*nz) - CSMF_x[j4] = CSMF_x[j4]+proj*nx - CSMF_y[j4] = CSMF_y[j4]+proj*ny - CSMF_z[j4] = CSMF_z[j4]+proj*nz + proj = inv_velz2kms * ( + CSMF_vx[j4] * nx + CSMF_vy[j4] * ny + CSMF_vz[j4] * nz + ) + CSMF_x[j4] = CSMF_x[j4] + proj * nx + CSMF_y[j4] = CSMF_y[j4] + proj * ny + CSMF_z[j4] = CSMF_z[j4] + proj * nz elif rsd: - CSMF_z[j4] = wrap(CSMF_z[j4] + CSMF_vz[j4] * inv_velz2kms, lbox) - - - M_1_C_temp = 10**(np.log10(M_1_C) + Ac_C * hdeltac[i] + Bc_C * hfenv[i]) + CSMF_z[j4] = wrap(CSMF_z[j4] + CSMF_vz[j4] * inv_velz2kms, lbox) + + M_1_C_temp = 10 ** ( + np.log10(M_1_C) + Ac_C * hdeltac[i] + Bc_C * hfenv[i] + ) a1_C_temp = a1_C + As_C * hdeltac[i] + Bs_C * hfenv[i] CSMF_stellarmass[j4] = get_random_sat_stellarmass( - hmass[i], - Mstar_low_C, - Mstar_up_C, - M_1_C_temp, - M_0_C, - gamma1_C, - gamma2_C, - a1_C_temp, - a2_C, - M2_C, - b0_C, - b1_C, - b2_C, - delta1_C, - delta2_C - ) + hmass[i], + Mstar_low_C, + Mstar_up_C, + M_1_C_temp, + M_0_C, + gamma1_C, + gamma2_C, + a1_C_temp, + a2_C, + M2_C, + b0_C, + b1_C, + b2_C, + delta1_C, + delta2_C, + ) CSMF_mass[j4] = hmass[i] CSMF_id[j4] = hid[i] j4 += 1 @@ -1744,7 +1938,7 @@ def gen_sats( LRG_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) ELG_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) QSO_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) - CSMF_dict = Dict.empty(key_type = types.unicode_type, value_type = float_array) + CSMF_dict = Dict.empty(key_type=types.unicode_type, value_type=float_array) ID_dict = Dict.empty(key_type=types.unicode_type, value_type=int_array) LRG_dict['x'] = lrg_x LRG_dict['y'] = lrg_y @@ -1772,7 +1966,7 @@ def gen_sats( QSO_dict['vz'] = qso_vz QSO_dict['mass'] = qso_mass ID_dict['QSO'] = qso_id - + CSMF_dict['x'] = CSMF_x CSMF_dict['y'] = CSMF_y CSMF_dict['z'] = CSMF_z @@ -1782,7 +1976,7 @@ def gen_sats( CSMF_dict['mass'] = CSMF_mass CSMF_dict['stellarmass'] = CSMF_stellarmass ID_dict['CSMF'] = CSMF_id - + return LRG_dict, ELG_dict, QSO_dict, CSMF_dict, ID_dict @@ -1992,11 +2186,13 @@ def gen_gals( QSO_hod_dict = nb.typed.Dict.empty( key_type=nb.types.unicode_type, value_type=nb.types.float64 ) - + if 'CSMF' in tracers.keys(): want_CSMF = True - CSMF_hod_dict = nb.typed.Dict.empty(key_type=nb.types.unicode_type, value_type= nb.types.float64) + CSMF_hod_dict = nb.typed.Dict.empty( + key_type=nb.types.unicode_type, value_type=nb.types.float64 + ) for key, value in CSMF_HOD.items(): CSMF_hod_dict[key] = value @@ -2009,7 +2205,9 @@ def gen_gals( else: want_CSMF = False - CSMF_hod_dict = nb.typed.Dict.empty(key_type=nb.types.unicode_type, value_type= nb.types.float64) + CSMF_hod_dict = nb.typed.Dict.empty( + key_type=nb.types.unicode_type, value_type=nb.types.float64 + ) start = time.time() @@ -2018,7 +2216,14 @@ def gen_gals( lbox = params['Lbox'] origin = params['origin'] - LRG_dict_cent, ELG_dict_cent, QSO_dict_cent, CSMF_dict_cent, ID_dict_cent, keep_cent = gen_cent( + ( + LRG_dict_cent, + ELG_dict_cent, + QSO_dict_cent, + CSMF_dict_cent, + ID_dict_cent, + keep_cent, + ) = gen_cent( halos_array['hpos'], halos_array['hvel'], halos_array['hmass'], @@ -2051,31 +2256,33 @@ def gen_gals( warnings.warn( 'NFW profile is unoptimized. It has different velocity bias. It does not support lightcone.' ) - LRG_dict_sat, ELG_dict_sat, QSO_dict_sat, CSMF_dict_sat, ID_dict_sat = gen_sats_nfw( - NFW_draw, - halos_array['hpos'], - halos_array['hvel'], - halos_array['hmass'], - halos_array['hid'], - halos_array.get('hdeltac', np.zeros(len(halos_array['hmass']))), - halos_array.get('hfenv', np.zeros(len(halos_array['hmass']))), - halos_array.get('hshear', np.zeros(len(halos_array['hmass']))), - halos_array['hsigma3d'], - halos_array['hc'], - halos_array['hrvir'], - LRG_hod_dict, - ELG_hod_dict, - QSO_hod_dict, - CSMF_hod_dict, - want_LRG, - want_ELG, - want_QSO, - want_CSMF, - rsd, - inv_velz2kms, - lbox, - keep_cent, - Nthread=Nthread, + LRG_dict_sat, ELG_dict_sat, QSO_dict_sat, CSMF_dict_sat, ID_dict_sat = ( + gen_sats_nfw( + NFW_draw, + halos_array['hpos'], + halos_array['hvel'], + halos_array['hmass'], + halos_array['hid'], + halos_array.get('hdeltac', np.zeros(len(halos_array['hmass']))), + halos_array.get('hfenv', np.zeros(len(halos_array['hmass']))), + halos_array.get('hshear', np.zeros(len(halos_array['hmass']))), + halos_array['hsigma3d'], + halos_array['hc'], + halos_array['hrvir'], + LRG_hod_dict, + ELG_hod_dict, + QSO_hod_dict, + CSMF_hod_dict, + want_LRG, + want_ELG, + want_QSO, + want_CSMF, + rsd, + inv_velz2kms, + lbox, + keep_cent, + Nthread=Nthread, + ) ) else: LRG_dict_sat, ELG_dict_sat, QSO_dict_sat, CSMF_dict_sat, ID_dict_sat = gen_sats( @@ -2115,8 +2322,18 @@ def gen_gals( print('generating satellites took ', time.time() - start) # B.H. TODO: need a for loop above so we don't need to do this by hand - HOD_dict_sat = {'LRG': LRG_dict_sat, 'ELG': ELG_dict_sat, 'QSO': QSO_dict_sat, 'CSMF': CSMF_dict_sat} - HOD_dict_cent = {'LRG': LRG_dict_cent, 'ELG': ELG_dict_cent, 'QSO': QSO_dict_cent, 'CSMF': CSMF_dict_cent} + HOD_dict_sat = { + 'LRG': LRG_dict_sat, + 'ELG': ELG_dict_sat, + 'QSO': QSO_dict_sat, + 'CSMF': CSMF_dict_sat, + } + HOD_dict_cent = { + 'LRG': LRG_dict_cent, + 'ELG': ELG_dict_cent, + 'QSO': QSO_dict_cent, + 'CSMF': CSMF_dict_cent, + } # do a concatenate in numba parallel start = time.time() diff --git a/abacusnbody/hod/abacus_hod.py b/abacusnbody/hod/abacus_hod.py index cdd9438..9d1962b 100644 --- a/abacusnbody/hod/abacus_hod.py +++ b/abacusnbody/hod/abacus_hod.py @@ -34,7 +34,7 @@ N_cen_QSO, N_sat_generic, n_cen_CSMF, - n_sat_CSMF + n_sat_CSMF, ) # TODO B.H.: staging can be shorter and prettier; perhaps asdf for h5 and ecsv? @@ -312,7 +312,7 @@ def staging(self): if ( ('ELG' not in self.tracers.keys()) and ('QSO' not in self.tracers.keys()) - and (not self.force_mt) + and (not self.force_mt) and ('CSMF' not in self.tracers.keys()) ): halofilename = subsample_dir / ( @@ -896,35 +896,34 @@ def compute_ngal(self, tracers=None, Nthread=16): ) ngal_dict[etracer] = newngal[0] + newngal[1] fsat_dict[etracer] = newngal[1] / (newngal[0] + newngal[1]) - - + elif etracer == 'CSMF': newngal = AbacusHOD._compute_ngal_CSMF( - self.logMbins, - self.deltacbins, - self.fenvbins, + self.logMbins, + self.deltacbins, + self.fenvbins, self.halo_mass_func, - tracer_hod['Mstar_low'], - tracer_hod['Mstar_up'], - tracer_hod['M_1'], - tracer_hod['M_0'], - tracer_hod['gamma_1'], - tracer_hod['gamma_2'], - tracer_hod['sigma_c'], - tracer_hod['a_1'], - tracer_hod['a_2'], - tracer_hod['M_2'], - tracer_hod['b_0'], - tracer_hod['b_1'], - tracer_hod['b_2'], - tracer_hod.get('delta_1', 0), - tracer_hod.get('delta_2', 0), + tracer_hod['Mstar_low'], + tracer_hod['Mstar_up'], + tracer_hod['M_1'], + tracer_hod['M_0'], + tracer_hod['gamma_1'], + tracer_hod['gamma_2'], + tracer_hod['sigma_c'], + tracer_hod['a_1'], + tracer_hod['a_2'], + tracer_hod['M_2'], + tracer_hod['b_0'], + tracer_hod['b_1'], + tracer_hod['b_2'], + tracer_hod.get('delta_1', 0), + tracer_hod.get('delta_2', 0), tracer_hod.get('Acent', 0), - tracer_hod.get('Asat', 0), - tracer_hod.get('Bcent', 0), - tracer_hod.get('Bsat', 0), - tracer_hod.get('ic', 1), - Nthread + tracer_hod.get('Asat', 0), + tracer_hod.get('Bcent', 0), + tracer_hod.get('Bsat', 0), + tracer_hod.get('ic', 1), + Nthread, ) ngal_dict[etracer] = newngal[0] + newngal[1] fsat_dict[etracer] = newngal[1] / (newngal[0] + newngal[1]) @@ -1130,35 +1129,89 @@ def _compute_ngal_qso( ngal_cent += halo_mass_func[i, j, k] * ncent_temp * ic ngal_sat += halo_mass_func[i, j, k] * nsat_temp * ic return ngal_cent, ngal_sat - - + @staticmethod - @njit(fastmath = True, parallel = True) - def _compute_ngal_CSMF(logMbins, deltacbins, fenvbins, halo_mass_func, - Mstar_low, Mstar_up, M_1, M_0, gamma1, gamma2, sigma_c, a1, a2, M2, b0, b1, b2, delta1, delta2, Acent, Asat, Bcent, Bsat, ic, Nthread): + @njit(fastmath=True, parallel=True) + def _compute_ngal_CSMF( + logMbins, + deltacbins, + fenvbins, + halo_mass_func, + Mstar_low, + Mstar_up, + M_1, + M_0, + gamma1, + gamma2, + sigma_c, + a1, + a2, + M2, + b0, + b1, + b2, + delta1, + delta2, + Acent, + Asat, + Bcent, + Bsat, + ic, + Nthread, + ): """ internal helper to compute number of CSMFs """ numba.set_num_threads(Nthread) - logMs = 0.5*(logMbins[1:] + logMbins[:-1]) - deltacs = 0.5*(deltacbins[1:] + deltacbins[:-1]) - fenvs = 0.5*(fenvbins[1:] + fenvbins[:-1]) + logMs = 0.5 * (logMbins[1:] + logMbins[:-1]) + deltacs = 0.5 * (deltacbins[1:] + deltacbins[:-1]) + fenvs = 0.5 * (fenvbins[1:] + fenvbins[:-1]) ngal_cent = 0 ngal_sat = 0 for i in numba.prange(len(logMbins) - 1): for j in range(len(deltacbins) - 1): for k in range(len(fenvbins) - 1): - Mh_temp = 10**logMs[i] - M_1_temp = 10**(np.log10(M_1)) #+ Acent * deltacs[j] + Bcent * fenvs[k]) - M2_temp = 10**(np.log10(M2)) #+ Asat * deltacs[j] + Bsat * fenvs[k]) - - ncent_temp = n_cen_CSMF(Mh_temp, Mstar_low, Mstar_up, M_1_temp, M_0, gamma1, gamma2, sigma_c) - nsat_temp = n_sat_CSMF(Mh_temp, Mstar_low, Mstar_up, M_1_temp, M_0, gamma1, gamma2, sigma_c, a1, a2, M2_temp, b0, b1, b2, delta1, delta2) + Mh_temp = 10 ** logMs[i] + M_1_temp = 10 ** ( + np.log10(M_1) + ) # + Acent * deltacs[j] + Bcent * fenvs[k]) + M2_temp = 10 ** ( + np.log10(M2) + ) # + Asat * deltacs[j] + Bsat * fenvs[k]) + + ncent_temp = n_cen_CSMF( + Mh_temp, + Mstar_low, + Mstar_up, + M_1_temp, + M_0, + gamma1, + gamma2, + sigma_c, + ) + nsat_temp = n_sat_CSMF( + Mh_temp, + Mstar_low, + Mstar_up, + M_1_temp, + M_0, + gamma1, + gamma2, + sigma_c, + a1, + a2, + M2_temp, + b0, + b1, + b2, + delta1, + delta2, + ) ngal_cent += halo_mass_func[i, j, k] * ncent_temp * ic ngal_sat += halo_mass_func[i, j, k] * nsat_temp * ic - + return ngal_cent, ngal_sat def compute_clustering(self, mock_dict, *args, **kwargs):