The search algorithm¶
This section will outline skysearcher parallel search algorithm step by step.
For each halo file¶
The first loop introduces the halo grid file. The grid file is a
numpy array.¶for halo, grid_fh in grids:
# Get a list of satids and a table for counting
# satids per region (satid_list) (satid_table).
satid_list, satid_table = satid_setup(halo)
# Get the data grid for this halo (grid).
grid = load_grid(grid_fh)
# Adjust nstars for units and distance.
grid[:, :, 1] /= mod
# <commented out>
# Kernel function
# grid[:, :, 1] = gaussian_filter(grid[:,:,0], sigma=0.8, mode="constant", cval=0, order=0)
""" End of first loop """
For each radius assigned to this process¶
The second for loop concerns the radial selection. Each process is assigned a list of radii for which they are responsible. The second for loop iterates over this list.
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 |
""" End of first loop """
# Step through the radii (r_start, r_stop).
for job_id in work_index:
# Start time (a_tic).
a_tic = time()
# r = radius in Kpc
# r_start = starting radius (always < r)
# r_stop = ending radius (always > r)
r, r_start, r_stop = _radii[job_id]
# This is so we dont need to index the whole table every
# loop of the following for loop (local_satid_table).
local_satid_table = satid_table[np.logical_and(
satid_table["Rads"] >= r_start,
satid_table["Rads"] < r_stop)]
# Get the mu for this annulus (mu).
mu, r_idx = mu_idx(grid, r_start, r_stop)
# <commented out>
# The number of boxes in the annulus (log_n_boxes_in_ann).
# log_n_boxes_in_ann = np.log10(len(np.nonzero(r_idx)[0]))
# <commented out>
# The number of stars in this annulus (log_n_stars_in_ann)
# log_n_stars_in_ann = np.log10(grid[:, :, 1][r_idx].sum())
# Boolean value to indicate the existence of
# a immediately previous accepted segment(one_before).
one_before = False
# Load array of annuli segments and
# step value (annuli) & (annuli_step).
# i.e. annuli_step = float
# i.e. annuli[x] = (deg_0, deg_1)
annuli, annuli_step = get_annuli(r)
# Number of empty segments (n_mt_seg).
n_mt_seg = 0
# Number of increasing segments (n_seg_increase).
# Number of decreasing segments (n_seg_decrease).
n_seg_increase = 0
n_seg_decrease = 0
last_nstars = 0
|
# Step through the annulus (-pi, pi].
# _deg0 = starting point.
# _deg1 = ending point
for _deg0, _deg1 in annuli:
# Get the mu for this sub-annulus-section (mu).
mu, r_idx2 = mu_idx2(grid, r_idx, _deg0, _deg1)
# The grid index for grid spaces within
# this segment (idx).
idx = get_idx(grid, _deg0, _deg1, r_idx2)
# The number of grid boxes this segment covers
# (n_boxes_tot).
n_boxes_in_seg = len(idx[0])
# If there are none, then continue.
if not n_boxes_in_seg:
"""
STDOUT.write(
"rank "
+ str(MPI_RANK)
+ " [ EMPTY SEGMENT ] "
+ halo
+ " radius: " + str(r) + " Kpc"
+ " phi: " + str(_deg0)
+ "\n")
STDOUT.flush()
"""
n_mt_seg += 1
continue
# The array of grid spaces from this segment"s
# contrast density value (xbox).
xbox = get_xbox(grid, idx, mu)
# Number of stars here (n_stars_here).
n_stars_here = grid[:, :, 1][idx].sum()
# TODO
# What is this doing?
if n_stars_here >= last_nstars:
n_seg_increase += 1
else:
n_seg_decrease += 1
last_nstars = n_stars_here
# ------------------------------------------------------
# Does this segment qualify to be a feature?
# ------------------------------------------------------
# Is the local min above xbox_cut?
# If yes:
if (xbox.min() >= XBOX_CUT
and np.log10(n_stars_here) >= MIN_LOG_NSTARS):
# If this is a new feature...
if not one_before:
# open a list for the max values
# for this feature (xbmax).
xbmax = []
# open a list for max Log10(n_stars) (n_stars_max).
log_n_stars_max = []
# Start a dictionary for counting stars
# per satid (sat_stars).
sat_stars = new_sat_stars(satid_list)
# Set grid space count (n_stars).
n_boxes = 0
# Set star count (n_stars).
n_stars = 0
# Set the run length (n_segments) to 0.
n_segments = 0
# Remember the starting angular
# value (starting_deg).
starting_deg = _deg0
# Set allowed segment skips (n_skips).
n_skips = N_SKIPS
# Do this for every accepted segment:
# -----------------------------------
# Put all region info into a list (r_info).
r_info = [r_start, r_stop, _deg0, _deg1]
# Add this features min, mean and max
# to the existing lists.
xbmax.append(xbox.max())
# Count stars per sat.
sat_stars = count_strs(
sat_stars,
r_info,
local_satid_table)
# Add boxes to total boxes.
n_boxes += n_boxes_in_seg
# Add stars to total stars for feature.
n_stars += n_stars_here
# Add n_stars to list for later.
log_n_stars_max.append(np.log10(n_stars_here))
# Increase run length (run_length) by 1.
n_segments += 1
# Remember that there is a feature
# currently being processed (one_before).
one_before = True
# If no:
else:
# Use allowed segment skips:
if one_before and n_skips:
n_skips -= 1
n_segments += 1
one_before = True
else:
# If this is the end of a segment:
if (one_before and n_segments >= MIN_N_SEGMENTS):
# The feature"s angular extent
# (angular_extent).
angular_extent = _deg0 - starting_deg
# The dominate satellite number (domsat_id).
# domsats"s % of all stars (domsat_purity).
_va = dom_satid(sat_stars)
domsat_id, domsat_purity, standout, nsats = _va
# An integer value for each halo (halo_num).
halo_num = int(halo[-2:])
# Mass of parent satellite (mass).
mass = m_book[domsat_id]
# Accretion time of parent satellite (atime).
atime = t_book[domsat_id]
# Circle factor (jcirc).
jcirc = j_book[domsat_id]
# A new row for the r_table (row).
# Each feature is a row in the table.
row = [
# Halo.
halo_num,
# Annulus location values.
r,
r_start,
r_stop,
annuli_step,
# n_mt_seg,
# Annulus content values.
# log_n_boxes_in_ann,
# log_n_stars_in_ann,
np.log10(mu),
# Feature content values.
max(xbmax),
max(log_n_stars_max),
domsat_purity,
domsat_id,
standout,
nsats,
mass,
atime,
jcirc,
# Feature location values.
starting_deg,
# _deg0,
angular_extent,
# n_segments,
# n_boxes,
# MPI values.
# MPI_RANK
]
# Add the row to the table.
r_table.add_row(row)
# Save point if added a new row.
if len(r_table) % SAVE_INTERVAL == 0:
save_record_table(_table=r_table)