Skip to content

Commit

Permalink
Revert "split/bin parameters in pycbc_fit_sngls_split_binned (gwastro…
Browse files Browse the repository at this point in the history 
…#4577)"

This reverts commit d7e72a7.
  • Loading branch information
@maxtrevor
maxtrevor committed Dec 11, 2023
1 parent f124032 commit 88281a2
Showing 1 changed file with 58 additions and 50 deletions.
108 changes: 58 additions & 50 deletions bin/all_sky_search/pycbc_fit_sngls_split_binned
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@ from pycbc.events import triggers as trigs
from pycbc.events import trigger_fits as trstats
from pycbc.events import stat as pystat
from pycbc.types.optparse import MultiDetOptionAction
from pycbc.tmpltbank import bank_conversions
import pycbc.version

parser = argparse.ArgumentParser(usage="",
Expand All @@ -41,7 +40,6 @@ parser.add_argument('--output-file', required=True,
help="Output image file. Required")
parser.add_argument('--verbose', action='store_true')
parser.add_argument('--bin-param', default='template_duration',
choices=bank_conversions.conversion_options,
help="Parameter for binning within plots, default "
"'template_duration'")
parser.add_argument('--bin-spacing', choices=['linear', 'log'], default='log',
Expand All @@ -52,11 +50,9 @@ parser.add_argument('--num-bins', type=int, default=6,
parser.add_argument('--max-bin-param', type=float, default=None,
help="Maximum allowed value of bin-param")
parser.add_argument('--split-param-one', default='eta',
choices=bank_conversions.conversion_options,
help="Parameter for splitting plot grid in y-direction, "
"default 'eta'")
parser.add_argument('--split-param-two', default='chi_eff',
choices=bank_conversions.conversion_options,
help="Parameter for splitting plot grid in x-direction, "
"default 'chi_eff'")
parser.add_argument('--split-one-nbins', default=3, type=int,
Expand Down Expand Up @@ -114,46 +110,60 @@ trigf = h5py.File(args.trigger_file, 'r')

logging.info('Opening template file: %s' % args.bank_file)
bank = h5py.File(args.bank_file, 'r')
logging.info('Getting template bank parameters')
logging.info('setting up template bank parameters')
params = {
"mass1" : bank['mass1'][:],
"mass2" : bank['mass2'][:],
"spin1z": bank['spin1z'][:],
"spin2z": bank['spin2z'][:]
}
bank.close()

# Calculate params needed
# Only calculate params if needed
usedparams = [args.split_param_two, args.split_param_one, args.bin_param]
extparams = [par for par in ['mchirp', 'eta', 'chi_eff', 'template_duration']
if par in usedparams]

params = {}
for par in usedparams:
if par in ['template_duration', 'duration']:
for ex_p in extparams:
if ex_p == 'template_duration':
logging.info('Reading duration from trigger file')
# Loop over template regions, return the first duration value from each region.
dur = []
for ref in trigf[args.ifo + '/template_duration_template'][:]:
dur.append(trigf[args.ifo + '/template_duration'][ref][0])
params['template_duration'] = np.array(dur)
# List comprehension loops over templates; if a template has no triggers, accessing
# the 0th entry of its region reference will return zero due to a quirk of h5py.
params[ex_p] = np.array(
[
trigf[args.ifo + '/template_duration'][ref][0]
for ref in trigf[args.ifo + '/template_duration_template'][:]
]
)
else:
logging.info("Calculating %s from template parameters", par)
params[par] = bank_conversions.get_bank_property(
par, bank, np.arange(bank['mass1'].size))

bank.close()

logging.info('setting up %s bins', ', '.join(usedparams))
logging.info("Calculating " + ex_p + " from template parameters")
params[ex_p] = trigs.get_param(ex_p, args, params['mass1'],
params['mass2'], params['spin1z'],
params['spin2z'])

# string formats for labels, logging etc.
formats = {
"mchirp": '{:.2f}',
"eta": '{:.3f}',
"chi_eff": '{:.3f}',
'template_duration': '{:.0f}'
}

logging.info('setting up {}'.format(usedparams))
sp_one_bin_input = (params[args.split_param_one].min(),
params[args.split_param_one].max(), args.split_one_nbins)
sp_two_bin_input = (params[args.split_param_two].min(),
params[args.split_param_two].max(), args.split_two_nbins)
if args.max_bin_param:
logging.info(
'setting maximum %s value: %.3f',
args.bin_param,
args.max_bin_param
)
logging.info(('setting maximum {} value: '
+ formats[args.bin_param]).format(args.bin_param,
args.max_bin_param))
pbin_upper_lim = float(args.max_bin_param)
else:
pbin_upper_lim = params[args.bin_param].max()

# For templates with no triggers, the duration will be read as zero
if args.bin_param == 'template_duration' and params[args.bin_param].min() == 0:
# Accessing the 0th entry of an empty region reference will return
# zero due to a quirk of h5py.
logging.warn('WARNING: Some templates do not contain triggers')
# Use the lowest nonzero template duration as lower limit for bins
pbin_lower_lim = params[args.bin_param][params[args.bin_param] > 0].min()
Expand All @@ -162,7 +172,7 @@ else:

bb_input = (pbin_lower_lim, pbin_upper_lim, args.num_bins)

logging.info('splitting %s into bins', args.bin_param)
logging.info('splitting {} into bins'.format(args.bin_param))
if args.bin_spacing == 'log':
assert pbin_lower_lim > 0
pbins = bin_utils.LogarithmicBins(*bb_input)
Expand Down Expand Up @@ -353,11 +363,8 @@ for x in range(args.split_one_nbins):
time_inbin = np.array(time_inbin)

count_pruned = 0
logging.info(
'Pruning in split %s-%i %s-%i',
args.split_param_one, x,
args.split_param_two, y
)
logging.info('Pruning in split {}-{} {}-{}'.format(
args.split_param_one, x, args.split_param_two, y))
logging.info('Currently have %d triggers', len(vals_inbin))
while count_pruned < args.prune_number:
# Getting loudest statistic value in split
Expand Down Expand Up @@ -407,7 +414,7 @@ fig, axes = plt.subplots(args.split_one_nbins, args.split_two_nbins,
lines = []
labels = []
for i, lower, upper in zip(range(args.num_bins), pbins.lower(), pbins.upper()):
binlabel = f"{lower:#.3g} - {upper:#.3g}"
binlabel = r"%.3g - %.3g" % (lower, upper)
line, = axes[0,0].plot([0,0], [0,0], linewidth=2,
color=histcolors[i], alpha=0.6)
lines.append(line)
Expand All @@ -431,16 +438,13 @@ for x in range(args.split_one_nbins):
for y in range(args.split_two_nbins):
id_bin2 = id_in_bin2[y]
id_in_both = np.intersect1d(id_bin1, id_bin2)
logging.info(
'Split %s-%i %s-%i',
args.split_param_one, x,
args.split_param_two, y,
)
logging.info('split {}: {}, {}: {}'.format(args.split_param_one, x,
args.split_param_two, y))
ax = axes[x,y]
for i, lower, upper in zip(range(args.num_bins), pbins.lower(),
pbins.upper()):
indices_all_conditions = np.intersect1d(pidx[i], id_in_both)
logging.info('%s split %#.3g-%#.3g', args.bin_param, lower, upper)
logging.info('{} split {}-{}'.format(args.bin_param, lower, upper))
if len(indices_all_conditions) == 0: continue
vals_inbin = []
for idx in indices_all_conditions:
Expand All @@ -452,8 +456,8 @@ for x in range(args.split_one_nbins):
logging.info('No triggers above threshold')
continue
else:
logging.info('%d triggers out of %d above threshold',
len(vals_above_thresh), len(vals_inbin))
logging.info('{} triggers out of {} above threshold'.format(
len(vals_above_thresh), len(vals_inbin)))
alpha, sig_alpha = trstats.fit_above_thresh(args.fit_function,
vals_above_thresh, args.stat_fit_threshold)
fitted_cum_counts = len(vals_above_thresh) * \
Expand All @@ -476,7 +480,7 @@ for x in range(args.split_one_nbins):
ax.semilogy(edges[:-1], cum_counts, linewidth=2,
color=histcolors[i], alpha=0.6)
ax.semilogy(fitrange, fitted_cum_counts, "--", color=histcolors[i],
label=f"$\\alpha = ${alpha:#.3f}" % alpha)
label=r"$\alpha = $%.2f" % alpha)
lgd_sub = ax.legend(fontsize='small', framealpha=0.5)

del vals_inbin, vals_above_thresh
Expand All @@ -495,21 +499,25 @@ for j in range(args.split_two_nbins):
axes[args.split_one_nbins - 1, j].set_xlabel(args.sngl_ranking, size="large")
if args.split_two_nbins == 1:
break
xrange_string = f'{sp_two_bounds.lower()[j]:#.3g} to {sp_two_bounds.upper()[j]:#.3g}'
axes[0, j].set_xlabel(f'{args.split_param_two}: {xrange_string}', size="large")
axes[0, j].set_xlabel(args.split_param_two + ': ' +
(formats[args.split_param_two] + ' to ' +
formats[args.split_param_two]).format(sp_two_bounds.lower()[j],
sp_two_bounds.upper()[j]), size="large")
axes[0, j].xaxis.set_label_position("top")

for i in range(args.split_one_nbins):
if args.split_one_nbins == 1:
axes[0, 0].set_ylabel('Cumulative number', size='large')
break
yrange_string = f'{sp_one_bounds.lower()[i]:#.3g} to {sp_one_bounds.upper()[i]:#.3g}'
axes[i, 0].set_ylabel(f'{args.split_param_one}: {yrange_string}\ncumulative number',
size="large")
axes[i, 0].set_ylabel(args.split_param_one + ': ' +
(formats[args.split_param_one] + ' to ' +
formats[args.split_param_one]).format(sp_one_bounds.lower()[i],
sp_one_bounds.upper()[i]) + '\ncumulative number',
size="large")

fig.tight_layout(rect=(1./(args.split_two_nbins+1), 0, 1, 1))

logging.info('Saving to file %s', args.output_file)
logging.info('Saving to file ' + args.output_file)
results.save_fig_with_metadata(
fig, args.output_file,
title="{}: {} histogram of single detector triggers split by"
Expand Down

0 comments on commit 88281a2

Please sign in to comment.