From 8ae6d9aa52cd4935e86a150bb6992aa105fc3e9d Mon Sep 17 00:00:00 2001 From: SBouchard01 Date: Thu, 4 Jun 2026 06:24:26 -0700 Subject: [PATCH 1/9] menv: replace print statement with logger for mass environment calculation --- abacusnbody/hod/menv.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/abacusnbody/hod/menv.py b/abacusnbody/hod/menv.py index f4099604..84116d23 100644 --- a/abacusnbody/hod/menv.py +++ b/abacusnbody/hod/menv.py @@ -9,12 +9,16 @@ import numpy as np from scipy.spatial import KDTree +import logging +from .utils import setup_logging + from ..util import cumsum __all__ = ['do_Menv_from_tree'] DEFAULT_BATCH_SIZE = 10**5 +logger = logging.getLogger(__file__) def do_Menv_from_tree( pos, @@ -53,7 +57,7 @@ def do_Menv_from_tree( if r_outer.ndim > 0: r_outer = r_outer[mmask] - print('Building and querying trees for mass env calculation') + logger.info('Building and querying trees for mass env calculation') tree = KDTree(pos, boxsize=treebox) # we're taking potentially large differences, use float64 From dd58fe72b0f0aef166acc217f44b610a7a403611 Mon Sep 17 00:00:00 2001 From: SBouchard01 Date: Thu, 23 Oct 2025 10:53:47 -0700 Subject: [PATCH 2/9] Change halo & particle filter for BGS --- abacusnbody/hod/prepare_sim.py | 103 ++++++--------------------------- 1 file changed, 17 insertions(+), 86 deletions(-) diff --git a/abacusnbody/hod/prepare_sim.py b/abacusnbody/hod/prepare_sim.py index b579ded5..f107a0fe 100644 --- a/abacusnbody/hod/prepare_sim.py +++ b/abacusnbody/hod/prepare_sim.py @@ -79,98 +79,29 @@ def unwrap_x_for_slab(x, i, numslabs, Lbox): # https://arxiv.org/pdf/2001.06018.pdf Figure 13 shows redshift evolution of LRG HOD -# standard power law satellites + +# NOTE : Changed subsample_halos & submask_particles to match BGS requirements in density, change proposed by Sandy Yuan. def subsample_halos(m, MT): x = np.log10(m) downfactors = np.zeros(len(x)) - if MT: - # for elgs - mask1 = x < 11.4 - mask2 = x < 11.6 - downfactors[mask1] = 0.2 / (1.0 + 10 * np.exp(-(x[mask1] - 11.2) * 25)) - downfactors[mask2 & (~mask1)] = 0.4 / ( - 1.0 + 10 * np.exp(-(x[mask2 & (~mask1)] - 11.3) * 25) - ) - downfactors[~mask2] = 1.0 / (1.0 + 0.1 * np.exp(-(x[~mask2] - 11.7) * 10)) - - # # for bgs - # mask1 = x < 11.0 - # mask2 = x < 11.2 - # downfactors[mask2&(~mask1)] = 0.1 # 0.4/(1.0 + 10*np.exp(-(x[mask2&(~mask1)] - 10.9)*25)) - # downfactors[~mask2] = 1 # 1.0/(1.0 + 0.1*np.exp(-(x[~mask2] - 11.3)*10)) - return downfactors - else: - downfactors = 1.0 / ( - 1.0 + 0.1 * np.exp(-(x - 11.8) * 10) - ) # LRG only, default 12.3, set to 12.0 for z = 1.1 - downfactors[x > 13.0] = 1 - return downfactors - - -# # new version for negative alpha ELG satellites -# def subsample_halos(m, MT): -# x = np.log10(m) -# downfactors = np.zeros(len(x)) -# if MT: -# mask1 = x < 11.35 -# mask2 = x < 11.52 -# downfactors[mask1] = 0.5/(1.0 + 10*np.exp(-(x[mask1] - 11.1)*25)) -# downfactors[mask2&(~mask1)] = 0.8/(1.0 + 10*np.exp(-(x[mask2&(~mask1)] - 11.25)*25)) -# downfactors[~mask2] = 1.0/(1.0 + 0.1*np.exp(-(x[~mask2] - 11.6)*10)) -# # save all halos that could host a satellite -# # downfactors[x>11.2] = 1 -# return downfactors -# else: -# downfactors = 1.0/(1.0 + 0.1*np.exp(-(x - 12.0)*10)) # LRG only, might be able to step back to 12.5 depending on bestfit -# downfactors[x > 13.0] = 1 -# return downfactors - -# # standard satellites -# def submask_particles(m_in, n_in, MT): -# x = np.log10(m_in) - -# if MT: -# if m_in < 1e11: -# return np.zeros(n_in) -# else: -# # a target number of particles -# ntarget = np.minimum(n_in, int(1 + 1.5*10**(x-13))) -# submask = np.zeros(n_in).astype(int) -# submask[np.random.choice(n_in, ntarget, replace = False)] = 1 -# return submask -# else: -# if 10**x < 1e12: -# return np.zeros(n_in) # essentially removing particles in halos below Mmin -# else: -# ntarget = np.minimum(n_in, int(1 + 1.5*10**(x-13))) -# submask = np.zeros(n_in).astype(int) -# submask[np.random.choice(n_in, ntarget, replace = False)] = 1 -# return submask - - -# conformity fix + # for bgs + mask1 = x < 11.0 + mask2 = x < 11.2 + downfactors[mask2&(~mask1)] = 0.1 # 0.4/(1.0 + 10*np.exp(-(x[mask2&(~mask1)] - 10.9)*25)) + downfactors[~mask2] = 1 # 1.0/(1.0 + 0.1*np.exp(-(x[~mask2] - 11.3)*10)) + return downfactors + def submask_particles(m_in, n_in, MT): x = np.log10(m_in) - - if MT: - if m_in < 1e11: - return np.zeros(n_in) - else: - # a target number of particles - # ntarget = np.minimum(n_in, int(1 + 1.5*10**(x-11.8))) - ntarget = np.minimum(n_in, int(1 + 1.5 * 10 ** (x - 12.5))) - ntarget = np.minimum(ntarget, 100) - submask = np.zeros(n_in).astype(int) - submask[np.random.choice(n_in, ntarget, replace=False)] = 1 - return submask + if m_in < 1e11: + return np.zeros(n_in) else: - if 10**x < 1e12: - return np.zeros(n_in) # essentially removing particles in halos below Mmin - else: - ntarget = np.minimum(n_in, int(1 + 1.5 * 10 ** (x - 13))) - submask = np.zeros(n_in).astype(int) - submask[np.random.choice(n_in, ntarget, replace=False)] = 1 - return submask + # a target number of particles + ntarget = np.minimum(n_in, int(1 + 1.5*10**(x-11.8))) + ntarget = np.minimum(ntarget, 100) + submask = np.zeros(n_in).astype(int) + submask[np.random.choice(n_in, ntarget, replace = False)] = 1 + return submask def get_vertices_cube(units=0.5, N=3): From ff33ff3e0e78e1444b5c4323425df9dc3df19cfe Mon Sep 17 00:00:00 2001 From: SBouchard01 Date: Thu, 23 Oct 2025 11:11:12 -0700 Subject: [PATCH 3/9] Refactor logging in prepare_sim.py to use logger instead of print statements --- abacusnbody/hod/prepare_sim.py | 73 +++++++++++++++++----------------- 1 file changed, 37 insertions(+), 36 deletions(-) diff --git a/abacusnbody/hod/prepare_sim.py b/abacusnbody/hod/prepare_sim.py index f107a0fe..0cb3ca40 100644 --- a/abacusnbody/hod/prepare_sim.py +++ b/abacusnbody/hod/prepare_sim.py @@ -30,9 +30,13 @@ from ..analysis.tsc import tsc_parallel from .menv import do_Menv_from_tree +import logging +from .utils import setup_logging + DEFAULTS = {} DEFAULTS['path2config'] = 'config/abacus_hod.yaml' +logger = logging.getLogger('prepare_sim.py') # ------------------------------------------------------------ # local env fix @@ -200,7 +204,7 @@ def gen_rand(N, chi_min, chi_max, fac, Lbox, offset, origins, rng): mask = mask0 | mask1 | mask2 else: mask = mask0 - print('masked randoms = ', np.sum(mask) * 100.0 / len(mask)) + logger.info('%f masked randoms', np.sum(mask) * 100.0 / len(mask)) rands_pos = np.vstack((x_cart[mask], y_cart[mask], z_cart[mask])).T rands_chis = rands_chis[mask] @@ -264,7 +268,7 @@ def prepare_slab( outfilename_env = ( savedir + '/env_xcom_' + str(i) + '_abacushod_localenv' + '_new.h5' ) - print('processing slab ', i) + logger.info('Processing slab %d', i) if MT: outfilename_halos += '_MT' outfilename_particles += '_MT' @@ -285,11 +289,11 @@ def prepare_slab( and (os.path.exists(outfilename_particles)) and ((not need_env_file) or os.path.exists(outfilename_env)) ): - print('files exists, skipping ', i) + logger.info('Files exist, skipping slab %d', i) return 0 # load the halo catalog slab - print('loading halo catalog ') + logger.info('Loading halo catalog ') if halo_lc: slabname = ( simdir @@ -377,7 +381,7 @@ def prepare_slab( # creating a mask of which halos to keep, which halos to drop p_halos = subsample_halos(halos['N'] * Mpart, MT) mask_halos = np.random.random(len(halos)) < p_halos - print('total number of halos, ', len(halos), 'keeping ', np.sum(mask_halos)) + logger.info('Total number of halos: %d, keeping %d', len(halos), np.sum(mask_halos)) halos['mask_subsample'] = mask_halos halos['multi_halos'] = 1.0 / p_halos @@ -686,7 +690,7 @@ def prepare_slab( # fenv will be computed later globally inside abacushod.py halos['fenv_rank'] = np.zeros(len(halos)) - print('computing c rank') + logger.info('Computing c rank (delta concentration)') halos_c = halos['r98_L2com'] / halos['r25_L2com'] deltac_rank = np.zeros(len(halos)) for ibin in range(nbins): @@ -723,7 +727,7 @@ def prepare_slab( new_shear_rank = halo_shears.argsort().argsort() shear_rank[mmask] = new_shear_rank / np.max(new_shear_rank) - 0.5 halos['shear_rank'] = shear_rank - print('finished shear compute') + logger.info('Finished shear compute') else: halos['shear_rank'] = np.zeros(len(halos)) @@ -753,11 +757,11 @@ def prepare_slab( fenvh_parts = np.full(len_old, -1.0) shearh_parts = np.full(len_old, -1.0) - print('compiling particle subsamples') + logger.info('Compiling particle subsamples') start_tracker = 0 for j in range(len(halos)): if j % 10000 == 0: - print('halo id', j, end='\r') + print('halo id', j, end='\r') # NOTE: left to print to monitor progress and avoid many lines of logging if mask_halos[j] and halos['npoutA'][j] > 0: # subsample_factor = subsample_particles(halos['N'][j] * Mpart, halos['npoutA'][j], MT) # submask = np.random.binomial(n = 1, p = subsample_factor, size = halos_pnum[j]) @@ -927,7 +931,7 @@ def prepare_slab( ) # attaching random numbers # output halo file - print('outputting new halo file ') + logger.info('Outputting new halo file ') # output_dir = savedir+'/halos_xcom_'+str(i)+'_seed'+str(newseed)+'_abacushodMT_new.h5' if os.path.exists(outfilename_halos): os.remove(outfilename_halos) @@ -937,15 +941,11 @@ def prepare_slab( # output the new particle file if z_type == 'primary' or z_type == 'lightcone': - print('adding rank fields to particle data ') + logger.info('Adding rank fields to particle data ') mask_parts = mask_parts.astype(bool) parts = parts[mask_parts] - print( - 'pre process particle number ', - len_old, - ' post process particle number ', - len(parts), - ) + logger.info('Pre-process particle number: %d, ', len_old) + logger.info('Post-process particle number: %d', len(parts)) if want_ranks: parts['ranks'] = ranks_parts[mask_parts] parts['ranksv'] = ranksv_parts[mask_parts] @@ -962,12 +962,12 @@ def prepare_slab( parts['halo_fenv'] = fenvh_parts[mask_parts] parts['halo_shear'] = shearh_parts[mask_parts] - print( - 'are there any negative particle values? ', + logger.info( + 'There are %d negative particles in downsample_halo and %d negative particles in halo_mass', np.sum(parts['downsample_halo'] < 0), np.sum(parts['halo_mass'] < 0), ) - print('outputting new particle file ') + logger.info('Outputting new particle file ') # output_dir = savedir+'/particles_xcom_'+str(i)+'_seed'+str(newseed)+'_abacushodMT_new.h5' if os.path.exists(outfilename_particles): os.remove(outfilename_particles) @@ -975,12 +975,8 @@ def prepare_slab( newfile.create_dataset('particles', data=parts) newfile.close() - print( - 'pre process particle number ', - len_old, - ' post process particle number ', - len(parts), - ) + logger.info('Pre-process particle number: %d', len_old) + logger.info('Post-process particle number: %d', len(parts)) def calc_shearmark(simdir, simname, z_mock, N_dim, R, fn, partdown=100): @@ -1028,7 +1024,7 @@ def calc_shearmark(simdir, simname, z_mock, N_dim, R, fn, partdown=100): ] pos_parts = np.concatenate(partpos) - print('compiled all particles', len(pos_parts), 'took time', time.time() - start) + logger.info('Compiled all particles: %d, took time: %f', len(pos_parts), time.time() - start) start = time.time() cat = CompaSOHaloCatalog( @@ -1039,19 +1035,19 @@ def calc_shearmark(simdir, simname, z_mock, N_dim, R, fn, partdown=100): header = cat.header Lbox = header['BoxSizeHMpc'] Lbox / N_dim - print('compiled all halos', 'took time', time.time() - start) + logger.info('Compiled all halos, took time: %f', time.time() - start) start = time.time() # dens = np.zeros((N_dim, N_dim, N_dim)) # numba_tsc_3D(pos_parts, dens, Lbox) dens = tsc_parallel(pos_parts, N_dim, Lbox) - print('finished TSC, took time', time.time() - start) + logger.info('Finished TSC, took time: %f', time.time() - start) start = time.time() dens_smooth = smooth_density(dens, R, N_dim, Lbox) - print('finished smoothing, took time', time.time() - start) + logger.info('Finished smoothing, took time: %f', time.time() - start) start = time.time() shearmark = get_shear(dens_smooth, N_dim, Lbox) - print('finished shear mark, took time', time.time() - start) + logger.info('Finished shear mark, took time: %f', time.time() - start) # output file np.save(fn + '.npy', shearmark) @@ -1067,7 +1063,7 @@ def main( halo_lc=False, overwrite=1, ): - print('compiling compaso halo catalogs into subsampled catalogs') + logger.info('Compiling compaso halo catalogs into subsampled catalogs') config = yaml.safe_load(open(path2config)) # update params if needed @@ -1162,7 +1158,7 @@ def main( if os.path.exists(shear_fn + '.npy'): shearmark = np.load(shear_fn + '.npy') else: - print('computing shear field') + logger.info('Computing shear field') shearmark = calc_shearmark( simdir, simname, z_mock, Ndim, Rsm, shear_fn, partdown ) @@ -1174,7 +1170,7 @@ def main( nthread = ( len(os.sched_getaffinity(0)) // config['prepare_sim']['Nparallel_load'] ) - print(f'prepare_sim inferred Nthread_per_load = {nthread}') + logger.info(f'prepare_sim inferred Nthread_per_load = {nthread}') else: nthread = int(nthread) @@ -1215,7 +1211,7 @@ def main( raise RuntimeError( 'A subprocess died in prepare_sim. Did prepare_slab() run out of memory?' ) from bpp - # print("done, took time ", time.time() - start) + # logger.info("Done, took time: %f", time.time() - start) class ArgParseFormatter( @@ -1232,6 +1228,9 @@ class ArgParseFormatter( parser.add_argument( '--path2config', help='Path to the config file', default=DEFAULTS['path2config'] ) + parser.add_argument( + '--path2log', help='Path to the log file', default=None + ) parser.add_argument( '--alt_simname', help='alternative simname to process, like "AbacusSummit_base_c000_ph003"', @@ -1251,6 +1250,8 @@ class ArgParseFormatter( '--overwrite', help='overwrite existing subsamples', default=1, type=int ) args = vars(parser.parse_args()) + setup_logging(filename=args.pop('path2log')) + main(**args) - print('done') + logger.info('Done !') \ No newline at end of file From 2452cf8f5f219c0788fb41b897bf5f4437c2f161 Mon Sep 17 00:00:00 2001 From: SBouchard01 Date: Thu, 23 Oct 2025 11:24:45 -0700 Subject: [PATCH 4/9] Add option to mask slab progress display (to avoid verbose) --- abacusnbody/hod/prepare_sim.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/abacusnbody/hod/prepare_sim.py b/abacusnbody/hod/prepare_sim.py index 0cb3ca40..7ecad5ba 100644 --- a/abacusnbody/hod/prepare_sim.py +++ b/abacusnbody/hod/prepare_sim.py @@ -248,6 +248,7 @@ def prepare_slab( mcut=1e11, rad_outer=10, numslabs=None, + show_slab_progress=False, ): outfilename_halos = ( savedir @@ -760,7 +761,7 @@ def prepare_slab( logger.info('Compiling particle subsamples') start_tracker = 0 for j in range(len(halos)): - if j % 10000 == 0: + if j % 10000 == 0 and show_slab_progress: print('halo id', j, end='\r') # NOTE: left to print to monitor progress and avoid many lines of logging if mask_halos[j] and halos['npoutA'][j] > 0: # subsample_factor = subsample_particles(halos['N'][j] * Mpart, halos['npoutA'][j], MT) @@ -1062,6 +1063,7 @@ def main( newseed=600, halo_lc=False, overwrite=1, + show_slab_progress=False, ): logger.info('Compiling compaso halo catalogs into subsampled catalogs') @@ -1199,6 +1201,7 @@ def main( nthread=nthread, overwrite=overwrite, numslabs=numslabs, + show_slab_progress=show_slab_progress, ) for i in range(numslabs) ] @@ -1249,6 +1252,12 @@ class ArgParseFormatter( parser.add_argument( '--overwrite', help='overwrite existing subsamples', default=1, type=int ) + parser.add_argument( + '--show_slab_progress', + help='Show progress of slab preparation (can be verbose)', + default=False, + action='store_true', + ) args = vars(parser.parse_args()) setup_logging(filename=args.pop('path2log')) From 3e98fa6dc8e7c71dfbed0c597156718df872851b Mon Sep 17 00:00:00 2001 From: SBouchard01 Date: Thu, 4 Jun 2026 06:57:39 -0700 Subject: [PATCH 5/9] prepare_sim: add back original subsample & submask methods w/ bgs filter override option when no MT --- abacusnbody/hod/prepare_sim.py | 115 +++++++++++++++++++++++++++------ 1 file changed, 95 insertions(+), 20 deletions(-) diff --git a/abacusnbody/hod/prepare_sim.py b/abacusnbody/hod/prepare_sim.py index 7ecad5ba..08dee169 100644 --- a/abacusnbody/hod/prepare_sim.py +++ b/abacusnbody/hod/prepare_sim.py @@ -36,7 +36,7 @@ DEFAULTS = {} DEFAULTS['path2config'] = 'config/abacus_hod.yaml' -logger = logging.getLogger('prepare_sim.py') +logger = logging.getLogger(__file__) # ------------------------------------------------------------ # local env fix @@ -83,29 +83,95 @@ def unwrap_x_for_slab(x, i, numslabs, Lbox): # https://arxiv.org/pdf/2001.06018.pdf Figure 13 shows redshift evolution of LRG HOD - -# NOTE : Changed subsample_halos & submask_particles to match BGS requirements in density, change proposed by Sandy Yuan. -def subsample_halos(m, MT): +# standard power law satellites +def subsample_halos(m, MT: bool, is_bgs: bool = False): x = np.log10(m) downfactors = np.zeros(len(x)) - # for bgs - mask1 = x < 11.0 - mask2 = x < 11.2 - downfactors[mask2&(~mask1)] = 0.1 # 0.4/(1.0 + 10*np.exp(-(x[mask2&(~mask1)] - 10.9)*25)) - downfactors[~mask2] = 1 # 1.0/(1.0 + 0.1*np.exp(-(x[~mask2] - 11.3)*10)) + if MT: # for elgs + mask1 = x < 11.4 + mask2 = x < 11.6 + downfactors[mask1] = 0.2 / (1.0 + 10 * np.exp(-(x[mask1] - 11.2) * 25)) + downfactors[mask2 & (~mask1)] = 0.4 / ( + 1.0 + 10 * np.exp(-(x[mask2 & (~mask1)] - 11.3) * 25) + ) + downfactors[~mask2] = 1.0 / (1.0 + 0.1 * np.exp(-(x[~mask2] - 11.7) * 10)) + elif is_bgs: # for bgs + mask1 = x < 11.0 + mask2 = x < 11.2 + downfactors[mask2&(~mask1)] = 0.1 # 0.4/(1.0 + 10*np.exp(-(x[mask2&(~mask1)] - 10.9)*25)) + downfactors[~mask2] = 1 # 1.0/(1.0 + 0.1*np.exp(-(x[~mask2] - 11.3)*10)) + else: + downfactors = 1.0 / ( + 1.0 + 0.1 * np.exp(-(x - 11.8) * 10) + ) # LRG only, default 12.3, set to 12.0 for z = 1.1 + downfactors[x > 13.0] = 1 return downfactors -def submask_particles(m_in, n_in, MT): + +# # new version for negative alpha ELG satellites +# def subsample_halos(m, MT): +# x = np.log10(m) +# downfactors = np.zeros(len(x)) +# if MT: +# mask1 = x < 11.35 +# mask2 = x < 11.52 +# downfactors[mask1] = 0.5/(1.0 + 10*np.exp(-(x[mask1] - 11.1)*25)) +# downfactors[mask2&(~mask1)] = 0.8/(1.0 + 10*np.exp(-(x[mask2&(~mask1)] - 11.25)*25)) +# downfactors[~mask2] = 1.0/(1.0 + 0.1*np.exp(-(x[~mask2] - 11.6)*10)) +# # save all halos that could host a satellite +# # downfactors[x>11.2] = 1 +# return downfactors +# else: +# downfactors = 1.0/(1.0 + 0.1*np.exp(-(x - 12.0)*10)) # LRG only, might be able to step back to 12.5 depending on bestfit +# downfactors[x > 13.0] = 1 +# return downfactors + +# # standard satellites +# def submask_particles(m_in, n_in, MT): +# x = np.log10(m_in) + +# if MT: +# if m_in < 1e11: +# return np.zeros(n_in) +# else: +# # a target number of particles +# ntarget = np.minimum(n_in, int(1 + 1.5*10**(x-13))) +# submask = np.zeros(n_in).astype(int) +# submask[np.random.choice(n_in, ntarget, replace = False)] = 1 +# return submask +# else: +# if 10**x < 1e12: +# return np.zeros(n_in) # essentially removing particles in halos below Mmin +# else: +# ntarget = np.minimum(n_in, int(1 + 1.5*10**(x-13))) +# submask = np.zeros(n_in).astype(int) +# submask[np.random.choice(n_in, ntarget, replace = False)] = 1 +# return submask + + +# conformity fix +def submask_particles(m_in, n_in, MT: bool, is_bgs: bool = False): x = np.log10(m_in) - if m_in < 1e11: - return np.zeros(n_in) - else: - # a target number of particles - ntarget = np.minimum(n_in, int(1 + 1.5*10**(x-11.8))) + + M_min = 1e11 if (MT or is_bgs) else 1e12 + + # Edge cases + if m_in < M_min: + return np.zeros(n_in) # essentially removing particles in halos below Mmin + + # a target number of particles + if MT: + ntarget = np.minimum(n_in, int(1 + 1.5 * 10 ** (x - 12.5))) + ntarget = np.minimum(ntarget, 100) + elif is_bgs: + ntarget = np.minimum(n_in, int(1 + 1.5 * 10 ** (x - 11.8))) ntarget = np.minimum(ntarget, 100) - submask = np.zeros(n_in).astype(int) - submask[np.random.choice(n_in, ntarget, replace = False)] = 1 - return submask + else: + ntarget = np.minimum(n_in, int(1 + 1.5 * 10 ** (x - 13))) + + submask = np.zeros(n_in).astype(int) + submask[np.random.choice(n_in, ntarget, replace=False)] = 1 + return submask def get_vertices_cube(units=0.5, N=3): @@ -248,6 +314,7 @@ def prepare_slab( mcut=1e11, rad_outer=10, numslabs=None, + is_bgs=False, show_slab_progress=False, ): outfilename_halos = ( @@ -380,7 +447,7 @@ def prepare_slab( # # form a halo table of the columns i care about # creating a mask of which halos to keep, which halos to drop - p_halos = subsample_halos(halos['N'] * Mpart, MT) + p_halos = subsample_halos(halos['N'] * Mpart, MT, is_bgs) mask_halos = np.random.random(len(halos)) < p_halos logger.info('Total number of halos: %d, keeping %d', len(halos), np.sum(mask_halos)) @@ -767,7 +834,7 @@ def prepare_slab( # subsample_factor = subsample_particles(halos['N'][j] * Mpart, halos['npoutA'][j], MT) # submask = np.random.binomial(n = 1, p = subsample_factor, size = halos_pnum[j]) submask = submask_particles( - halos['N'][j] * Mpart, halos['npoutA'][j], MT + halos['N'][j] * Mpart, halos['npoutA'][j], MT, is_bgs ) # updating the particles' masks, downsample factors, halo mass @@ -1063,6 +1130,7 @@ def main( newseed=600, halo_lc=False, overwrite=1, + is_bgs=False, show_slab_progress=False, ): logger.info('Compiling compaso halo catalogs into subsampled catalogs') @@ -1201,6 +1269,7 @@ def main( nthread=nthread, overwrite=overwrite, numslabs=numslabs, + is_bgs=is_bgs, show_slab_progress=show_slab_progress, ) for i in range(numslabs) @@ -1252,6 +1321,12 @@ class ArgParseFormatter( parser.add_argument( '--overwrite', help='overwrite existing subsamples', default=1, type=int ) + parser.add_argument( + '--is_bgs', + help='Use the BGS filter instead of LRG if Multi-Tracer option is not enabled.', + default=False, + action='store_true', + ) parser.add_argument( '--show_slab_progress', help='Show progress of slab preparation (can be verbose)', From 8b1cf2112a7e41318b828489d8458ea54e383fa6 Mon Sep 17 00:00:00 2001 From: SBouchard01 Date: Thu, 4 Jun 2026 06:58:03 -0700 Subject: [PATCH 6/9] add __init__.py file to abacusnbody.analysis module to allow import when used outside of the folder --- abacusnbody/analysis/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 abacusnbody/analysis/__init__.py diff --git a/abacusnbody/analysis/__init__.py b/abacusnbody/analysis/__init__.py new file mode 100644 index 00000000..e69de29b From 4f2c379dd9e76b726b5034160955620c7f15fcbd Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 4 Jun 2026 14:20:22 +0000 Subject: [PATCH 7/9] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- abacusnbody/hod/menv.py | 2 +- abacusnbody/hod/prepare_sim.py | 37 +++++++++++++++++++--------------- 2 files changed, 22 insertions(+), 17 deletions(-) diff --git a/abacusnbody/hod/menv.py b/abacusnbody/hod/menv.py index 84116d23..a0259014 100644 --- a/abacusnbody/hod/menv.py +++ b/abacusnbody/hod/menv.py @@ -10,7 +10,6 @@ from scipy.spatial import KDTree import logging -from .utils import setup_logging from ..util import cumsum @@ -20,6 +19,7 @@ logger = logging.getLogger(__file__) + def do_Menv_from_tree( pos, mass, diff --git a/abacusnbody/hod/prepare_sim.py b/abacusnbody/hod/prepare_sim.py index 08dee169..fb237162 100644 --- a/abacusnbody/hod/prepare_sim.py +++ b/abacusnbody/hod/prepare_sim.py @@ -38,6 +38,7 @@ logger = logging.getLogger(__file__) + # ------------------------------------------------------------ # local env fix # ------------------------------------------------------------ @@ -87,7 +88,7 @@ def unwrap_x_for_slab(x, i, numslabs, Lbox): def subsample_halos(m, MT: bool, is_bgs: bool = False): x = np.log10(m) downfactors = np.zeros(len(x)) - if MT: # for elgs + if MT: # for elgs mask1 = x < 11.4 mask2 = x < 11.6 downfactors[mask1] = 0.2 / (1.0 + 10 * np.exp(-(x[mask1] - 11.2) * 25)) @@ -95,11 +96,13 @@ def subsample_halos(m, MT: bool, is_bgs: bool = False): 1.0 + 10 * np.exp(-(x[mask2 & (~mask1)] - 11.3) * 25) ) downfactors[~mask2] = 1.0 / (1.0 + 0.1 * np.exp(-(x[~mask2] - 11.7) * 10)) - elif is_bgs: # for bgs + elif is_bgs: # for bgs mask1 = x < 11.0 mask2 = x < 11.2 - downfactors[mask2&(~mask1)] = 0.1 # 0.4/(1.0 + 10*np.exp(-(x[mask2&(~mask1)] - 10.9)*25)) - downfactors[~mask2] = 1 # 1.0/(1.0 + 0.1*np.exp(-(x[~mask2] - 11.3)*10)) + downfactors[mask2 & (~mask1)] = ( + 0.1 # 0.4/(1.0 + 10*np.exp(-(x[mask2&(~mask1)] - 10.9)*25)) + ) + downfactors[~mask2] = 1 # 1.0/(1.0 + 0.1*np.exp(-(x[~mask2] - 11.3)*10)) else: downfactors = 1.0 / ( 1.0 + 0.1 * np.exp(-(x - 11.8) * 10) @@ -152,12 +155,12 @@ def subsample_halos(m, MT: bool, is_bgs: bool = False): # conformity fix def submask_particles(m_in, n_in, MT: bool, is_bgs: bool = False): x = np.log10(m_in) - + M_min = 1e11 if (MT or is_bgs) else 1e12 - + # Edge cases if m_in < M_min: - return np.zeros(n_in) # essentially removing particles in halos below Mmin + return np.zeros(n_in) # essentially removing particles in halos below Mmin # a target number of particles if MT: @@ -168,7 +171,7 @@ def submask_particles(m_in, n_in, MT: bool, is_bgs: bool = False): ntarget = np.minimum(ntarget, 100) else: ntarget = np.minimum(n_in, int(1 + 1.5 * 10 ** (x - 13))) - + submask = np.zeros(n_in).astype(int) submask[np.random.choice(n_in, ntarget, replace=False)] = 1 return submask @@ -829,7 +832,9 @@ def prepare_slab( start_tracker = 0 for j in range(len(halos)): if j % 10000 == 0 and show_slab_progress: - print('halo id', j, end='\r') # NOTE: left to print to monitor progress and avoid many lines of logging + print( + 'halo id', j, end='\r' + ) # NOTE: left to print to monitor progress and avoid many lines of logging if mask_halos[j] and halos['npoutA'][j] > 0: # subsample_factor = subsample_particles(halos['N'][j] * Mpart, halos['npoutA'][j], MT) # submask = np.random.binomial(n = 1, p = subsample_factor, size = halos_pnum[j]) @@ -1031,7 +1036,7 @@ def prepare_slab( parts['halo_shear'] = shearh_parts[mask_parts] logger.info( - 'There are %d negative particles in downsample_halo and %d negative particles in halo_mass', + 'There are %d negative particles in downsample_halo and %d negative particles in halo_mass', np.sum(parts['downsample_halo'] < 0), np.sum(parts['halo_mass'] < 0), ) @@ -1092,7 +1097,9 @@ def calc_shearmark(simdir, simname, z_mock, N_dim, R, fn, partdown=100): ] pos_parts = np.concatenate(partpos) - logger.info('Compiled all particles: %d, took time: %f', len(pos_parts), time.time() - start) + logger.info( + 'Compiled all particles: %d, took time: %f', len(pos_parts), time.time() - start + ) start = time.time() cat = CompaSOHaloCatalog( @@ -1300,9 +1307,7 @@ class ArgParseFormatter( parser.add_argument( '--path2config', help='Path to the config file', default=DEFAULTS['path2config'] ) - parser.add_argument( - '--path2log', help='Path to the log file', default=None - ) + parser.add_argument('--path2log', help='Path to the log file', default=None) parser.add_argument( '--alt_simname', help='alternative simname to process, like "AbacusSummit_base_c000_ph003"', @@ -1335,7 +1340,7 @@ class ArgParseFormatter( ) args = vars(parser.parse_args()) setup_logging(filename=args.pop('path2log')) - + main(**args) - logger.info('Done !') \ No newline at end of file + logger.info('Done !') From 9f7ac1760e18a1970000f0f485d63e8a99519222 Mon Sep 17 00:00:00 2001 From: SBouchard01 Date: Thu, 4 Jun 2026 07:38:04 -0700 Subject: [PATCH 8/9] forgotten print to logger statements from last update --- abacusnbody/hod/prepare_sim.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/abacusnbody/hod/prepare_sim.py b/abacusnbody/hod/prepare_sim.py index fb237162..57079f2d 100644 --- a/abacusnbody/hod/prepare_sim.py +++ b/abacusnbody/hod/prepare_sim.py @@ -723,7 +723,7 @@ def prepare_slab( env_id = env_id[uniq_idx] nbr_count = len(env_mass) - Ncentral - print( + logger.info( f'[slab {i}] env centers = {Ncentral:,}, ' f'neighbor halos = {nbr_count:,}, ' f'total env halos = {len(env_mass):,}, ' @@ -745,7 +745,7 @@ def prepare_slab( Menv_central = Menv_all[:Ncentral] n_nonzero = np.count_nonzero(Menv_central) - print( + logger.info( f'[slab {i}] computed padded Menv for {Ncentral:,} central halos; ' f'nonzero Menv count = {n_nonzero:,}' ) From cbe0dc14c27d6a84115fda4505e3534e8ecd4d62 Mon Sep 17 00:00:00 2001 From: SBouchard01 Date: Fri, 19 Jun 2026 02:32:03 -0700 Subject: [PATCH 9/9] fix: replace logger names to use __name__ --- abacusnbody/hod/menv.py | 2 +- abacusnbody/hod/prepare_sim.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/abacusnbody/hod/menv.py b/abacusnbody/hod/menv.py index a0259014..d141b36a 100644 --- a/abacusnbody/hod/menv.py +++ b/abacusnbody/hod/menv.py @@ -17,7 +17,7 @@ DEFAULT_BATCH_SIZE = 10**5 -logger = logging.getLogger(__file__) +logger = logging.getLogger(__name__) def do_Menv_from_tree( diff --git a/abacusnbody/hod/prepare_sim.py b/abacusnbody/hod/prepare_sim.py index 57079f2d..0378b010 100644 --- a/abacusnbody/hod/prepare_sim.py +++ b/abacusnbody/hod/prepare_sim.py @@ -36,7 +36,7 @@ DEFAULTS = {} DEFAULTS['path2config'] = 'config/abacus_hod.yaml' -logger = logging.getLogger(__file__) +logger = logging.getLogger(__name__) # ------------------------------------------------------------