diff --git a/bin/all_sky_search/pycbc_fit_sngls_split_binned b/bin/all_sky_search/pycbc_fit_sngls_split_binned index daa992f6c13..961ec6f5892 100644 --- a/bin/all_sky_search/pycbc_fit_sngls_split_binned +++ b/bin/all_sky_search/pycbc_fit_sngls_split_binned @@ -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="", @@ -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', @@ -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, @@ -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() @@ -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) @@ -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 @@ -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) @@ -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: @@ -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) * \ @@ -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 @@ -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"