From 1c1a91aa3f602d79717175c04d993a502b7e0dbb Mon Sep 17 00:00:00 2001 From: Ian Harry Date: Fri, 2 Feb 2024 12:43:16 +0000 Subject: [PATCH] Additions for v2.3.4 release (#4616) * Fix error querying Virgo frame type (#4523) * Fix error querying Virgo frame type * Typo * Typo * Implement Ian's suggestion * Make it work * Use an actual DeprecationWarning * Scripts for dealing with preparation of offline gracedb upload files (#4534) * Scripts for dealing with preparation of offline gracedb upload files * Update bin/all_sky_search/pycbc_make_bayestar_skymap Co-authored-by: Tito Dal Canton * reinstate SNR timeseries being saved into HDF * TDC comments * TDC comments 2 * Remove unneeded import --------- Co-authored-by: Tito Dal Canton Co-authored-by: Tito Dal Canton * Use vetoed times followup in example (#4597) * Use vetoed times followup in example * Fix statistic files double-definition * Rationalize some calls to waveform properties (#4540) * rationalize some calls to waveform properties * Only allow explictly listed names * don't try to change function name * g g * Add the ability to choose how to sort injections in minifollowups (#4602) * Add the ability to choose how to sort injections in minifollowups * Minor cleanup * Minor refactor * Tom's comments * Update example search workflow * Further thoughts and suggestions * Tom's suggestion * Fix bug and Tom's suggestion * Standardise logging: bin/all_sky_search and pycbc/events (#4576) * add level to default logging * Decrease logging level for debug information in coinc_findtrigs * decrease logging level for some information in sngls_findtrigs * add named logger in pycbc/events/cuts.py * add named logger in pycbc/events/significance.py * REVERT ME: an example of usage * CC fix * Revert "REVERT ME: an example of usage" This reverts commit eb647e0b164c56091b4ad41a5c92df3e04ecc307. * Use init_logging throughout all_sky_search * give pycbc/events modules named loggers * missed the fact that the level was set to warn, not info, as default before * CC * missed an import * start moving verbose argparse into common facility * add common_options to all_sky_search executables * set --verbose default value to zero * pycbc not imported in some codes * CC * rename function to be more decriptive * O4b idq offline update (#4603) * Start offline dq workflow rework * Modified stat * Rewrote dq workflow * Added note for future suggested changes to dq workflow * Finished first draft of offline workflow * Fixed a comment * Removed unneeded argument to make dq workflow * Made bin templates executable * WOrking version of offline dq workflow * Reverting non-functional changes in stat.py * linting * Reverted more non-functional changes * Reverted low-latency related features * Rearrange imports * Addressed Ian's comments * Simplified masking in dq workflow * Modify dq workflow to avoid using numpy arrays * Use vetoed times followup in example (#4597) * Use vetoed times followup in example * Fix statistic files double-definition * Addressed more comments from Tom * Addressed Gareth's comments on parser arguments * Don't allow dq stat to uprank candidates * Modify dq trigger rate plots to not use dq trigger rates below 1 * Address Tom's most recent comment * Readded [Min 1] to dq plot's y-axis label * Pass bank plotting tags to dq template bin plots * Update bin/plotting/pycbc_plot_dq_flag_likelihood Co-authored-by: Thomas Dent * Update bin/workflows/pycbc_make_offline_search_workflow Co-authored-by: Thomas Dent * Use abs() correctly * Break up 2 try/ecept cases --------- Co-authored-by: Gareth S Cabourn Davies Co-authored-by: Thomas Dent * Bumping version number --------- Co-authored-by: Tito Dal Canton Co-authored-by: Gareth S Cabourn Davies Co-authored-by: Tito Dal Canton Co-authored-by: Thomas Dent Co-authored-by: maxtrevor <65971534+maxtrevor@users.noreply.github.com> --- bin/all_sky_search/pycbc_add_statmap | 2 +- bin/all_sky_search/pycbc_apply_rerank | 2 +- bin/all_sky_search/pycbc_average_psd | 2 +- bin/all_sky_search/pycbc_bin_templates | 56 ++++ bin/all_sky_search/pycbc_bin_trigger_rates_dq | 254 +++++++++-------- bin/all_sky_search/pycbc_calculate_dq | 85 ------ bin/all_sky_search/pycbc_calculate_dqflag | 73 ----- bin/all_sky_search/pycbc_calculate_psd | 2 +- bin/all_sky_search/pycbc_coinc_findtrigs | 12 +- bin/all_sky_search/pycbc_coinc_hdfinjfind | 8 +- bin/all_sky_search/pycbc_coinc_mergetrigs | 6 +- bin/all_sky_search/pycbc_coinc_statmap | 2 +- bin/all_sky_search/pycbc_coinc_statmap_inj | 9 +- .../pycbc_combine_coincident_events | 2 +- bin/all_sky_search/pycbc_combine_statmap | 2 +- .../pycbc_cut_merge_triggers_to_tmpltbank | 3 +- .../pycbc_distribute_background_bins | 2 +- bin/all_sky_search/pycbc_dtphase | 2 +- bin/all_sky_search/pycbc_exclude_zerolag | 2 +- bin/all_sky_search/pycbc_fit_sngls_binned | 4 +- .../pycbc_fit_sngls_by_template | 14 +- .../pycbc_fit_sngls_over_multiparam | 7 +- bin/all_sky_search/pycbc_fit_sngls_over_param | 10 +- .../pycbc_fit_sngls_split_binned | 2 +- bin/all_sky_search/pycbc_followup_file | 2 +- bin/all_sky_search/pycbc_foreground_censor | 2 +- bin/all_sky_search/pycbc_get_loudest_params | 4 + bin/all_sky_search/pycbc_make_bayestar_skymap | 85 ++++++ bin/all_sky_search/pycbc_merge_psds | 2 +- bin/all_sky_search/pycbc_plot_kde_vals | 6 +- .../pycbc_prepare_xml_for_gracedb | 238 ++++++++++++++++ bin/all_sky_search/pycbc_reduce_template_bank | 3 +- bin/all_sky_search/pycbc_rerank_dq | 121 --------- bin/all_sky_search/pycbc_rerank_passthrough | 3 +- bin/all_sky_search/pycbc_sngls_findtrigs | 8 +- bin/all_sky_search/pycbc_sngls_pastro | 2 +- bin/all_sky_search/pycbc_sngls_statmap | 2 +- bin/all_sky_search/pycbc_sngls_statmap_inj | 2 +- bin/all_sky_search/pycbc_strip_injections | 2 +- bin/all_sky_search/pycbc_template_kde_calc | 7 +- bin/all_sky_search/pycbc_template_kde_max | 6 +- .../pycbc_template_recovery_hist | 7 +- .../pycbc_upload_single_event_to_gracedb | 182 +++++++++++++ .../pycbc_injection_minifollowup | 193 +++++++++---- bin/plotting/pycbc_plot_dq_flag_likelihood | 87 +++--- .../pycbc_make_offline_search_workflow | 54 ++-- examples/search/plotting.ini | 18 +- pycbc/__init__.py | 23 +- pycbc/events/coinc.py | 52 ++-- pycbc/events/coinc_rate.py | 8 +- pycbc/events/cuts.py | 32 ++- pycbc/events/eventmgr.py | 32 ++- pycbc/events/significance.py | 18 +- pycbc/events/single.py | 6 +- pycbc/events/stat.py | 154 +++++++---- pycbc/events/trigger_fits.py | 6 +- pycbc/frame/frame.py | 125 +++++++-- pycbc/io/hdf.py | 4 + pycbc/io/live.py | 31 +-- pycbc/pnutils.py | 12 +- pycbc/results/__init__.py | 1 + pycbc/results/snr.py | 73 +++++ pycbc/strain/strain.py | 15 +- pycbc/workflow/dq.py | 257 +++++++----------- pycbc/workflow/plotting.py | 100 +++---- setup.py | 2 +- 66 files changed, 1516 insertions(+), 1034 deletions(-) create mode 100755 bin/all_sky_search/pycbc_bin_templates delete mode 100644 bin/all_sky_search/pycbc_calculate_dq delete mode 100644 bin/all_sky_search/pycbc_calculate_dqflag create mode 100644 bin/all_sky_search/pycbc_make_bayestar_skymap create mode 100755 bin/all_sky_search/pycbc_prepare_xml_for_gracedb delete mode 100644 bin/all_sky_search/pycbc_rerank_dq create mode 100755 bin/all_sky_search/pycbc_upload_single_event_to_gracedb create mode 100644 pycbc/results/snr.py diff --git a/bin/all_sky_search/pycbc_add_statmap b/bin/all_sky_search/pycbc_add_statmap index 5ea7a4cad17..1b8b7006e67 100755 --- a/bin/all_sky_search/pycbc_add_statmap +++ b/bin/all_sky_search/pycbc_add_statmap @@ -25,9 +25,9 @@ def get_ifo_string(fi): return istring parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action="version", version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--statmap-files', nargs='+', help="List of coinc files to be combined") parser.add_argument('--background-files', nargs='+', default=None, diff --git a/bin/all_sky_search/pycbc_apply_rerank b/bin/all_sky_search/pycbc_apply_rerank index 5a25d0bcd62..1946ea004d1 100644 --- a/bin/all_sky_search/pycbc_apply_rerank +++ b/bin/all_sky_search/pycbc_apply_rerank @@ -8,9 +8,9 @@ from pycbc.events import significance from shutil import copyfile parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--stat-files', nargs='+', help="Statistic files produced by candidate followup codes") parser.add_argument('--followup-file', diff --git a/bin/all_sky_search/pycbc_average_psd b/bin/all_sky_search/pycbc_average_psd index d85ebc29442..8070c5d2153 100644 --- a/bin/all_sky_search/pycbc_average_psd +++ b/bin/all_sky_search/pycbc_average_psd @@ -31,8 +31,8 @@ from pycbc.types import MultiDetOptionAction, FrequencySeries parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=version) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--input-files', nargs='+', required=True, metavar='PATH', help='HDF5 files from pycbc_calculate_psd (one per ' 'detector) containing the input PSDs to average.') diff --git a/bin/all_sky_search/pycbc_bin_templates b/bin/all_sky_search/pycbc_bin_templates new file mode 100755 index 00000000000..8da9c9bc6f5 --- /dev/null +++ b/bin/all_sky_search/pycbc_bin_templates @@ -0,0 +1,56 @@ +#!/usr/bin/env python +""" Bin templates by their duration +""" +import logging +import argparse +import h5py as h5 +import numpy as np + +import pycbc +import pycbc.pnutils +from pycbc.version import git_verbose_msg as version +from pycbc.events import background_bin_from_string + +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument('--version', action='version', version=version) +parser.add_argument('--verbose', action="count") +parser.add_argument("--ifo", type=str, required=True) +parser.add_argument("--f-lower", type=float, default=15., + help='Enforce a uniform low frequency cutoff to ' + 'calculate template duration over the bank') +parser.add_argument('--bank-file', help='hdf format template bank file', + required=True) +parser.add_argument('--background-bins', nargs='+', + help='Used to provide a list of ' + 'precomputed background bins') +parser.add_argument("--output-file", required=True) + +args = parser.parse_args() + +pycbc.init_logging(args.verbose) +logging.info('Starting template binning') + +with h5.File(args.bank_file, 'r') as bank: + logging.info('Sorting bank into bins') + data = { + 'mass1': bank['mass1'][:], + 'mass2': bank['mass2'][:], + 'spin1z': bank['spin1z'][:], + 'spin2z': bank['spin2z'][:], + 'f_lower': np.ones_like(bank['mass1'][:]) * args.f_lower + } + + bin_dict = background_bin_from_string(args.background_bins, data) + bin_names = [b.split(':')[0] for b in args.background_bins] + +logging.info('Writing bin template ids to file') +with h5.File(args.output_file, 'w') as f: + ifo_grp = f.create_group(args.ifo) + for bin_name in bin_names: + bin_tids = bin_dict[bin_name] + grp = ifo_grp.create_group(bin_name) + grp['tids'] = bin_tids + f.attrs['f_lower'] = args.f_lower + f.attrs['background_bins'] = args.background_bins + +logging.info('Finished') diff --git a/bin/all_sky_search/pycbc_bin_trigger_rates_dq b/bin/all_sky_search/pycbc_bin_trigger_rates_dq index cd18fd3863c..1bb84299bd0 100644 --- a/bin/all_sky_search/pycbc_bin_trigger_rates_dq +++ b/bin/all_sky_search/pycbc_bin_trigger_rates_dq @@ -3,52 +3,48 @@ """ import logging import argparse -import pycbc -import pycbc.events -from pycbc.events import stat as pystat -from pycbc.types.timeseries import load_timeseries + import numpy as np import h5py as h5 + +from ligo.segments import segmentlist + +import pycbc +from pycbc.events import stat as pystat +from pycbc.events.veto import (select_segments_by_definer, + start_end_to_segments, + segments_to_start_end) +from pycbc.types.optparse import MultiDetOptionAction from pycbc.io.hdf import HFile, SingleDetTriggers from pycbc.version import git_verbose_msg as version parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=version) -parser.add_argument('--verbose', action="store_true") -parser.add_argument("--ifo", type=str, required=True) +parser.add_argument("--template-bins-file", required=True) parser.add_argument("--trig-file", required=True) +parser.add_argument("--flag-file", required=True) +parser.add_argument("--flag-name", required=True) +parser.add_argument("--analysis-segment-file", required=True) +parser.add_argument("--analysis-segment-name", required=True) +parser.add_argument("--gating-windows", nargs='+', + action=MultiDetOptionAction, + help="Seconds to reweight before and after the central" + "time of each gate. Given as detector-values pairs, " + "e.g. H1:-1,2.5 L1:-1,2.5 V1:0,0") parser.add_argument("--stat-threshold", type=float, default=1., help="Only consider triggers with --sngl-ranking value " "above this threshold") -parser.add_argument("--f-lower", type=float, default=15., - help='Enforce a uniform low frequency cutoff to ' - 'calculate template duration over the bank') -parser.add_argument("--dq-file", required=True, nargs='+') -parser.add_argument("--dq-channel", required=False, type=str, - help='name of channel to read in from the provided ' - 'dq file. Required if file is .hdf5 format') -parser.add_argument('--bank-file', help='hdf format template bank file', - required=True) -parser.add_argument('--background-bins', nargs='+', - help='list of background bin format strings') -parser.add_argument('--n-time-bins', type=int, default=200, - help='Number of time bins to use') parser.add_argument("--output-file", required=True) -parser.add_argument("--prune-number", type=int, default=0, - help="Number of loudest events to remove from each " - "split histogram, default 0") -parser.add_argument("--prune-window", type=float, default=0.1, - help="Time (s) to remove all triggers around a trigger " - "which is loudest in each split, default 0.1s") - -pystat.insert_statistic_option_group(parser, - default_ranking_statistic='single_ranking_only') + +pystat.insert_statistic_option_group( + parser, default_ranking_statistic='single_ranking_only') args = parser.parse_args() pycbc.init_logging(args.verbose) logging.info('Start') -ifo = args.ifo +ifo, flag_name = args.flag_name.split(':') # Setup a data mask to remove any triggers with SNR below threshold # This works as a pre-filter as SNR is always greater than or equal @@ -75,6 +71,18 @@ with HFile(args.trig_file, 'r') as trig_file: data_mask = np.zeros(n_triggers_orig, dtype=bool) data_mask[idx] = True + # also get the gated times while we have the file open + gate_times = [] + if args.gating_windows: + logging.info('Getting gated times') + try: + gating_types = trig_file[f'{ifo}/gating'].keys() + for gt in gating_types: + gate_times += list(trig_file[f'{ifo}/gating/{gt}/time'][:]) + gate_times = np.unique(gate_times) + except KeyError: + logging.warning('No gating found in trigger file') + logging.info("Getting %s triggers from file with pre-cut SNR > %.3f", idx.size, args.stat_threshold) @@ -92,9 +100,6 @@ trigs = SingleDetTriggers( tmplt_ids = trigs.template_id trig_times = trigs.end_time stat = trigs.get_ranking(args.sngl_ranking) -trig_times_int = trig_times.astype(np.int64) - -del trigs n_triggers = tmplt_ids.size @@ -103,108 +108,95 @@ logging.info("Applying %s > %.3f cut", args.sngl_ranking, keep = stat >= args.stat_threshold tmplt_ids = tmplt_ids[keep] trig_times = trig_times[keep] -trig_times_int = trig_times_int[keep] logging.info("Removed %d triggers, %d remain", n_triggers - tmplt_ids.size, tmplt_ids.size) -dq_logl = np.array([]) -dq_times = np.array([]) - -for filename in args.dq_file: - logging.info('Reading DQ file %s', filename) - dq_data = load_timeseries(filename, group=args.dq_channel) - dq_logl = np.concatenate((dq_logl, dq_data[:])) - dq_times = np.concatenate((dq_times, dq_data.sample_times)) - del dq_data - -n_bins = args.n_time_bins -percentiles = np.linspace(0, 100, n_bins+1) -bin_times = np.zeros(n_bins) -dq_percentiles = np.percentile(dq_logl, percentiles)[1:] - -# seconds bin tells what bin each second ends up -seconds_bin = np.array([n_bins - np.sum(dq_percentiles >= dq_ll) - for dq_ll in dq_logl]).astype(np.int64) -del dq_percentiles - -# bin times tells how much time ends up in each bin -bin_times = np.array([np.sum(seconds_bin == i) - for i in range(n_bins)]).astype(np.float64) -full_time = float(len(seconds_bin)) -times_nz = (bin_times > 0) -del dq_logl - -# create a dict to look up dq percentile at any time -dq_percentiles_time = dict(zip(dq_times, seconds_bin/n_bins)) -del dq_times - -with h5.File(args.bank_file, 'r') as bank: - if args.background_bins: - logging.info('Sorting bank into bins') - data = { - 'mass1': bank['mass1'][:], - 'mass2': bank['mass2'][:], - 'spin1z': bank['spin1z'][:], - 'spin2z': bank['spin2z'][:], - 'f_lower': np.ones_like(bank['mass1'][:]) * args.f_lower - } - locs_dict = pycbc.events.background_bin_from_string( - args.background_bins, data) - del data - locs_names = [b.split(':')[0] for b in args.background_bins] - else: - locs_dict = {'all_bin': np.arange(0, len(bank['mass1'][:]), 1)} - locs_names = ['all_bin'] - -if args.prune_number > 0: - for bin_name in locs_names: - logging.info('Pruning bin %s', bin_name) - bin_locs = locs_dict[bin_name] - inbin = np.isin(tmplt_ids, bin_locs) - trig_times_bin = trig_times[inbin] - trig_stats_bin = stat[inbin] - - for j in range(args.prune_number): - max_stat_arg = np.argmax(trig_stats_bin) - remove = abs(trig_times_bin[max_stat_arg] - trig_times) \ - < args.prune_window - logging.info("Prune %d: pruning %d triggers", j, sum(remove)) - remove_inbin = abs(trig_times_bin[max_stat_arg] - - trig_times_bin) < args.prune_window - stat[remove] = 0 - trig_stats_bin[remove_inbin] = 0 - keep = np.flatnonzero(stat) - logging.info("%d triggers removed through pruning", keep.size) - trig_times_int = trig_times_int[keep] - tmplt_ids = tmplt_ids[keep] - del stat - del keep - -del trig_times - +# Get the template bins +bin_tids_dict = {} +with h5.File(args.template_bins_file, 'r') as f: + ifo_grp = f[ifo] + for bin_name in ifo_grp.keys(): + bin_tids_dict[bin_name] = ifo_grp[bin_name]['tids'][:] + + +# get analysis segments +analysis_segs = select_segments_by_definer( + args.analysis_segment_file, + segment_name=args.analysis_segment_name, + ifo=ifo) + +livetime = abs(analysis_segs) + +# get flag segments +flag_segs = select_segments_by_definer(args.flag_file, + segment_name=flag_name, + ifo=ifo) + +# construct gate segments +gating_segs = segmentlist([]) +if args.gating_windows: + gating_windows = args.gating_windows[ifo].split(',') + gate_before = float(gating_windows[0]) + gate_after = float(gating_windows[1]) + if gate_before > 0 or gate_after < 0: + raise ValueError("Gating window values must be negative " + "before gates and positive after gates.") + if not (gate_before == 0 and gate_after == 0): + gating_segs = start_end_to_segments( + gate_times + gate_before, + gate_times + gate_after + ).coalesce() + +# make segments into mutually exclusive dq states +gating_segs = gating_segs & analysis_segs +flag_segs = flag_segs & analysis_segs + +dq_state_segs_dict = {} +dq_state_segs_dict[2] = gating_segs +dq_state_segs_dict[1] = flag_segs - gating_segs +dq_state_segs_dict[0] = analysis_segs - flag_segs - gating_segs + + +# utility function to get the dq state at a given time +def dq_state_at_time(t): + for state, segs in dq_state_segs_dict.items(): + if t in segs: + return state + return None + + +# compute and save results with h5.File(args.output_file, 'w') as f: - for bin_name in locs_names: - bin_locs = locs_dict[bin_name] - inbin = np.isin(tmplt_ids, bin_locs) - trig_times_bin = trig_times_int[inbin] - trig_percentile = np.array([dq_percentiles_time[t] - for t in trig_times_bin]) - logging.info('Processing %d triggers in bin %s', - len(trig_percentile), bin_name) - - (counts, bins) = np.histogram(trig_percentile, bins=(percentiles)/100) - counts = counts.astype(np.float64) - rates = np.zeros_like(bin_times) - rates[times_nz] = counts[times_nz]/bin_times[times_nz] - mean_rate = len(trig_percentile) / full_time - if mean_rate > 0.: - rates = rates / mean_rate - - logging.info('Writing rates to output file %s', args.output_file) - grp = f.create_group(bin_name) - grp['rates'] = rates - grp['locs'] = locs_dict[bin_name] - - f.attrs['names'] = locs_names + ifo_grp = f.create_group(ifo) + all_bin_grp = ifo_grp.create_group('bins') + all_dq_grp = ifo_grp.create_group('dq_segments') + + # setup data for each template bin + for bin_name, bin_tids in bin_tids_dict.items(): + bin_grp = all_bin_grp.create_group(bin_name) + bin_grp['tids'] = bin_tids + + # get the dq states of the triggers in this bin + inbin = np.isin(tmplt_ids, bin_tids) + trig_times_bin = trig_times[inbin] + trig_states = np.array([dq_state_at_time(t) for t in trig_times_bin]) + + # calculate the dq rates in this bin + dq_rates = np.zeros(3, dtype=np.float64) + for state, segs in dq_state_segs_dict.items(): + frac_eff = np.mean(trig_states == state) + frac_dt = abs(segs) / livetime + dq_rates[state] = frac_eff / frac_dt + bin_grp['dq_rates'] = dq_rates + + # save dq state segments + for dq_state, segs in dq_state_segs_dict.items(): + name = f'dq_state_{dq_state}' + dq_grp = all_dq_grp.create_group(name) + starts, ends = segments_to_start_end(segs) + dq_grp['segment_starts'] = starts + dq_grp['segment_ends'] = ends + + f.attrs['stat'] = f'{ifo}-dq_stat_info' logging.info('Done!') diff --git a/bin/all_sky_search/pycbc_calculate_dq b/bin/all_sky_search/pycbc_calculate_dq deleted file mode 100644 index 92fbb2df647..00000000000 --- a/bin/all_sky_search/pycbc_calculate_dq +++ /dev/null @@ -1,85 +0,0 @@ -#!/usr/bin/env python -""" Ingest and process DQ data for use in the ranking statistic. - Any data that can be read by PyCBC as a timeseries can - be read in as a DQ timeseries. - This DQ timeseries is then processed to - lower the sample rate, either by maximizing - or calculating the BLRMS over each timestep. -""" -import logging, argparse, numpy, copy -import pycbc, pycbc.strain -from pycbc.version import git_verbose_msg as version -from pycbc.workflow import resolve_td_option -from pycbc.types.timeseries import TimeSeries -from ligo.segments import segment - -parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument('--version', action='version', version=version) -parser.add_argument('--verbose', action="store_true") -parser.add_argument("--dq-sample-rate", default=1, type=int, - help="Sample rate to calculate dq timeseries, " - "in Hz [default=1Hz]") -parser.add_argument("--dq-sample-method", required=True, - choices=['max', 'blrms'], - help="Method used to calculate dq timeseries. " - "Options are 'max' or 'blrms'") -parser.add_argument("--gps-start-time", type=int,required=True, - help="The gps start time of the data " - "(integer seconds)") -parser.add_argument("--gps-end-time", type=int,required=True, - help="The gps end time of the data " - "(integer seconds)") -parser.add_argument("--output-file", required=True, - help="name of output file") -parser.add_argument("--output-channel", required=False, - help="name of output channel") - -pycbc.strain.insert_strain_option_group(parser, gps_times=False) - -parser.set_defaults(pad_data=0) - -args = parser.parse_args() -pycbc.init_logging(args.verbose) - -seg = (args.gps_start_time,args.gps_end_time) - -logging.info('Getting dq values for %.1f-%.1f (%.1f s)', seg[0], - seg[1], abs(seg[1]-seg[0])) -argstmp = copy.deepcopy(args) -argstmp.gps_start_time = int(seg[0]) -argstmp.gps_end_time = int(seg[1]) -tmp_segment = segment([argstmp.gps_start_time, argstmp.gps_end_time]) -argstmp.channel_name = resolve_td_option(args.channel_name, tmp_segment) - -dq_data = pycbc.strain.from_cli(argstmp) - -logging.info('Normalizing dq values for %.1f-%.1f (%.1f s)', seg[0], - seg[1], abs(seg[1]-seg[0])) - -dq_sample_rate = float(args.dq_sample_rate) -dq_step_size = dq_data.sample_rate/dq_sample_rate - -dq_method = args.dq_sample_method - -if dq_method == "max": - # maximize dq_data over given interval - val_dq = numpy.array([max(dq_data.numpy() - [int(n*dq_step_size): - int((n+1)*dq_step_size)]) - for n in numpy.arange(0, len(dq_data.numpy()) / - dq_step_size)]) -elif dq_method=="blrms": - # Calculate the blrms over given interval - val_dq = numpy.array([numpy.mean(dq_data.numpy() - [int(n*dq_step_size): - int((n+1)*dq_step_size)]**2)**0.5 - for n in numpy.arange(0, len(dq_data.numpy()) / - dq_step_size)]) - -dq = TimeSeries(val_dq, epoch=dq_data.start_time, - delta_t=1./dq_sample_rate) - -dq.save(args.output_file, group=args.output_channel) - -logging.info('Done!') - diff --git a/bin/all_sky_search/pycbc_calculate_dqflag b/bin/all_sky_search/pycbc_calculate_dqflag deleted file mode 100644 index 981d190ad33..00000000000 --- a/bin/all_sky_search/pycbc_calculate_dqflag +++ /dev/null @@ -1,73 +0,0 @@ -#!/usr/bin/env python - -""" - Convert veto definer into timeseries -""" - -import logging, argparse, numpy, h5py -import pycbc -from pycbc.version import git_verbose_msg as version -from pycbc.events import veto -from pycbc.types.timeseries import TimeSeries - -parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument('--version', action='version', version=version) -parser.add_argument('--verbose', action="store_true") -parser.add_argument("--gps-start-time", type=int,required=True, - help="The gps start time of the data " - "(integer seconds)") -parser.add_argument("--gps-end-time", type=int,required=True, - help="The gps end time of the data " - "(integer seconds)") -parser.add_argument('--dq-segments', required=True, - help="segment file containing the " - "relevant data quality flag segments") -parser.add_argument('--flag', type=str, required=True, - help="name of the data quality flag") -parser.add_argument('--ifo', type=str, required=True, - help="interferometer to analyze") -parser.add_argument('--sample-frequency', required=True, type=int, - help="timeseries sampling frequency in Hz") -parser.add_argument("--output-file", required=True, - help="name of output file") -parser.add_argument("--output-channel", required=False, - help="name of output channel") - -args = parser.parse_args() -pycbc.init_logging(args.verbose) - - -def make_dq_ts(dq_segs, start, end, ifo, flag, freq): - """ Create a data quality timeseries - """ - logging.info('Creating data quality timeseries for flag %s', - flag) - # Keep just times which belong to science_segments - dur = end - start - dq_times = numpy.linspace(start, end, int(freq*dur), endpoint=False) - dq = numpy.zeros(len(dq_times)) - # Identify times within segments for the chosen flag - indexes, _ = veto.indices_within_segments(dq_times, [dq_segs], ifo=ifo, - segment_name=flag) - if indexes.size: - dq[indexes] = 1 - else: - logging.warning('Veto definer segment list is empty for flag %s-%s', - ifo, flag) - return TimeSeries(dq, epoch=start, delta_t=1.0/freq) - - -ifo = args.ifo -flag = args.flag - -flag_list = flag.split(':') -if len(flag_list)>1: - flag = flag_list[1] - -# Create data quality time series -dq = make_dq_ts(args.dq_segments, args.gps_start_time, args.gps_end_time, - ifo, flag, args.sample_frequency) - -dq.save(args.output_file, group=args.output_channel) - -logging.info('Done!') diff --git a/bin/all_sky_search/pycbc_calculate_psd b/bin/all_sky_search/pycbc_calculate_psd index 412c4163d30..29b6ab824a0 100755 --- a/bin/all_sky_search/pycbc_calculate_psd +++ b/bin/all_sky_search/pycbc_calculate_psd @@ -11,8 +11,8 @@ from ligo.segments import segmentlist, segment set_measure_level(0) parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=version) -parser.add_argument('--verbose', action="store_true") parser.add_argument("--low-frequency-cutoff", type=float, required=True, help="The low frequency cutoff to use for filtering (Hz)") parser.add_argument("--analysis-segment-file", required=True, diff --git a/bin/all_sky_search/pycbc_coinc_findtrigs b/bin/all_sky_search/pycbc_coinc_findtrigs index 234da72f53f..d76e3b9109a 100644 --- a/bin/all_sky_search/pycbc_coinc_findtrigs +++ b/bin/all_sky_search/pycbc_coinc_findtrigs @@ -4,13 +4,13 @@ import shutil, uuid, os.path, atexit from ligo.segments import infinity from pycbc.events import veto, coinc, stat, ranking, cuts import pycbc.version -from pycbc import pool +from pycbc import pool, init_logging from numpy.random import seed, shuffle from pycbc.io.hdf import ReadByTemplate from pycbc.types.optparse import MultiDetOptionAction parser = argparse.ArgumentParser() -parser.add_argument("--verbose", action="count") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) parser.add_argument("--veto-files", nargs='*', action='append', default=[], help="Optional veto file. Triggers within veto segments " @@ -78,9 +78,7 @@ args.segment_name = sum(args.segment_name, []) args.veto_files = sum(args.veto_files, []) args.trigger_files = sum(args.trigger_files, []) -if args.verbose: - logging.basicConfig(format='%(asctime)s : %(message)s', level=logging.DEBUG) - +init_logging(args.verbose) def parse_template_range(num_templates, rangestr): part = int(rangestr.split('/')[0]) @@ -238,7 +236,7 @@ def process_template(tnum): times_full = {} sds_full = {} tids_full = {} - logging.info('Obtaining trigs for template %i ..' % (tnum)) + logging.debug('Obtaining trigs for template %i ..' % (tnum)) for i, sngl in zip(trigs.ifos, trigs.singles): # Apply cuts to triggers tids_uncut = sngl.set_template(tnum) @@ -247,7 +245,7 @@ def process_template(tnum): tids_full[i] = tids_uncut[trigger_keep_ids] times_full[i] = sngl['end_time'][trigger_keep_ids] - logging.info('%s:%s', i, len(tids_uncut)) + logging.debug('%s:%s', i, len(tids_uncut)) if len(tids_full[i]) < len(tids_uncut): logging.info("%s triggers cut", len(tids_uncut) - len(tids_full[i])) diff --git a/bin/all_sky_search/pycbc_coinc_hdfinjfind b/bin/all_sky_search/pycbc_coinc_hdfinjfind index c9a6aab5cc2..0dfbe809958 100755 --- a/bin/all_sky_search/pycbc_coinc_hdfinjfind +++ b/bin/all_sky_search/pycbc_coinc_hdfinjfind @@ -6,7 +6,7 @@ files. import argparse, h5py, logging, types, numpy, os.path from ligo.lw import lsctables, utils as ligolw_utils from ligo import segments -from pycbc import events +from pycbc import events, init_logging from pycbc.events import indices_within_segments from pycbc.types import MultiDetOptionAction from pycbc.inject import CBCHDFInjectionSet @@ -55,6 +55,7 @@ def xml_to_hdf(table, hdf_file, hdf_key, columns): dtype=numpy.float32)) parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) parser.add_argument('--trigger-files', nargs='+', required=True) @@ -73,13 +74,10 @@ parser.add_argument('--optimal-snr-column', nargs='+', parser.add_argument('--redshift-column', default=None, help='Name of sim_inspiral column containing redshift. ' 'Optional') -parser.add_argument('--verbose', action='count') parser.add_argument('--output-file', required=True) args = parser.parse_args() -if args.verbose: - log_level = logging.INFO - logging.basicConfig(format='%(asctime)s : %(message)s', level=log_level) +init_logging(args.verbose) fo = h5py.File(args.output_file, 'w') diff --git a/bin/all_sky_search/pycbc_coinc_mergetrigs b/bin/all_sky_search/pycbc_coinc_mergetrigs index c5fd24c283c..c37456e725c 100755 --- a/bin/all_sky_search/pycbc_coinc_mergetrigs +++ b/bin/all_sky_search/pycbc_coinc_mergetrigs @@ -5,7 +5,7 @@ import numpy, argparse, h5py, logging import pycbc.version - +from pycbc import init_logging def changes(arr): l = numpy.where(arr[:-1] != arr[1:])[0] @@ -30,6 +30,7 @@ def region(f, key, boundaries, ids): dtype=h5py.special_dtype(ref=h5py.RegionReference)) parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) parser.add_argument('--trigger-files', nargs='+') @@ -38,10 +39,9 @@ parser.add_argument('--bank-file', required=True) parser.add_argument('--compression-level', type=int, default=6, help='Set HDF compression level in the output file ' '(default 6)') -parser.add_argument('--verbose', '-v', action='count') args = parser.parse_args() -logging.basicConfig(format='%(asctime)s : %(message)s', level=logging.INFO) +init_logging(args.verbose) f = h5py.File(args.output_file, 'w') diff --git a/bin/all_sky_search/pycbc_coinc_statmap b/bin/all_sky_search/pycbc_coinc_statmap index bfe97b866ac..d7b661606ff 100755 --- a/bin/all_sky_search/pycbc_coinc_statmap +++ b/bin/all_sky_search/pycbc_coinc_statmap @@ -34,6 +34,7 @@ class fw(object): parser = argparse.ArgumentParser() # General required options +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) parser.add_argument('--coinc-files', nargs='+', @@ -41,7 +42,6 @@ parser.add_argument('--coinc-files', nargs='+', 'FAP, FAR, etc.') parser.add_argument('--ifos', nargs='+', help='List of ifos used in these coincidence files') -parser.add_argument('--verbose', action='count') parser.add_argument('--cluster-window', type=float, default=10, help='Length of time window in seconds to cluster coinc ' 'events [default=10s]') diff --git a/bin/all_sky_search/pycbc_coinc_statmap_inj b/bin/all_sky_search/pycbc_coinc_statmap_inj index 03e4eab26c6..ed1fd819264 100644 --- a/bin/all_sky_search/pycbc_coinc_statmap_inj +++ b/bin/all_sky_search/pycbc_coinc_statmap_inj @@ -8,10 +8,11 @@ import argparse, h5py, logging, itertools, copy, pycbc.io, numpy, lal from pycbc.events import veto, coinc, significance import pycbc.version import pycbc.conversions as conv +from pycbc import init_logging parser = argparse.ArgumentParser() # General required options -parser.add_argument('--verbose', action='count') +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) parser.add_argument('--cluster-window', type=float, default=10, @@ -31,11 +32,9 @@ significance.insert_significance_option_group(parser) parser.add_argument('--output-file') args = parser.parse_args() -significance.check_significance_options(args, parser) +init_logging(args.verbose) -if args.verbose: - log_level = logging.INFO - logging.basicConfig(format='%(asctime)s : %(message)s', level=log_level) +significance.check_significance_options(args, parser) window = args.cluster_window diff --git a/bin/all_sky_search/pycbc_combine_coincident_events b/bin/all_sky_search/pycbc_combine_coincident_events index b35f748939b..4cf57fe2eb9 100644 --- a/bin/all_sky_search/pycbc_combine_coincident_events +++ b/bin/all_sky_search/pycbc_combine_coincident_events @@ -55,8 +55,8 @@ def com_with_detector_key(f, files, group): parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--statmap-files', nargs='+', help="List of coinc files to be redistributed") parser.add_argument('--output-file', help="name of output file") diff --git a/bin/all_sky_search/pycbc_combine_statmap b/bin/all_sky_search/pycbc_combine_statmap index 4e451ecfdf8..182559dab0d 100755 --- a/bin/all_sky_search/pycbc_combine_statmap +++ b/bin/all_sky_search/pycbc_combine_statmap @@ -9,8 +9,8 @@ import pycbc.version from ligo import segments parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--statmap-files', nargs='+', help="List of coinc files to be redistributed") parser.add_argument('--cluster-window', type=float) diff --git a/bin/all_sky_search/pycbc_cut_merge_triggers_to_tmpltbank b/bin/all_sky_search/pycbc_cut_merge_triggers_to_tmpltbank index 86625b4dfb3..e24d787e4ba 100644 --- a/bin/all_sky_search/pycbc_cut_merge_triggers_to_tmpltbank +++ b/bin/all_sky_search/pycbc_cut_merge_triggers_to_tmpltbank @@ -29,10 +29,9 @@ import pycbc import pycbc.version parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) -parser.add_argument("-V", "--verbose", action="store_true", - help="print extra debugging information", default=False) parser.add_argument("--input-file", required=True, help="Input merge triggers HDF file.") parser.add_argument("--output-file", required=True, diff --git a/bin/all_sky_search/pycbc_distribute_background_bins b/bin/all_sky_search/pycbc_distribute_background_bins index 1a2a4e56855..3a58dcab8b8 100644 --- a/bin/all_sky_search/pycbc_distribute_background_bins +++ b/bin/all_sky_search/pycbc_distribute_background_bins @@ -4,7 +4,7 @@ import pycbc.version parser = argparse.ArgumentParser() parser.add_argument("--version", action=pycbc.version.Version) -parser.add_argument('--verbose', action='store_true') +pycbc.add_common_pycbc_options(parser) parser.add_argument('--coinc-files', nargs='+', help="List of coinc files to be redistributed") parser.add_argument('--background-bins', nargs='+', diff --git a/bin/all_sky_search/pycbc_dtphase b/bin/all_sky_search/pycbc_dtphase index bdf74e9a7ab..e9e2f66f7da 100644 --- a/bin/all_sky_search/pycbc_dtphase +++ b/bin/all_sky_search/pycbc_dtphase @@ -20,6 +20,7 @@ from copy import deepcopy parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--ifos', nargs='+', help="The ifos to generate a histogram for") parser.add_argument('--sample-size', type=int, required=True, @@ -34,7 +35,6 @@ parser.add_argument('--relative-sensitivities', nargs='+', type=float, "expected SNR at fixed distance, one for each ifo") parser.add_argument('--seed', type=int, default=124) parser.add_argument('--output-file', required=True) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--bin-density', type=int, default=1, help="Number of bins per 1 sigma uncertainty in a " "parameter. Higher values increase the resolution of " diff --git a/bin/all_sky_search/pycbc_exclude_zerolag b/bin/all_sky_search/pycbc_exclude_zerolag index 82de048c973..c5b71571ef1 100644 --- a/bin/all_sky_search/pycbc_exclude_zerolag +++ b/bin/all_sky_search/pycbc_exclude_zerolag @@ -10,9 +10,9 @@ from pycbc.events import veto, coinc, significance import pycbc.conversions as conv parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--statmap-file', type=str, help="Coinc statmap file to be recalculated based on foreground removal") parser.add_argument('--other-statmap-files', nargs='+', diff --git a/bin/all_sky_search/pycbc_fit_sngls_binned b/bin/all_sky_search/pycbc_fit_sngls_binned index b65830df7b3..41a210acbaa 100644 --- a/bin/all_sky_search/pycbc_fit_sngls_binned +++ b/bin/all_sky_search/pycbc_fit_sngls_binned @@ -33,10 +33,8 @@ import pycbc.version parser = argparse.ArgumentParser(usage="", description="Perform maximum-likelihood fits of single inspiral trigger" " distributions to various functions") - +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action=pycbc.version.Version) -parser.add_argument("-V", "--verbose", action="store_true", - help="Print extra debugging information", default=False) parser.add_argument("--trigger-file", help="Input hdf5 file containing single triggers. " "Required") diff --git a/bin/all_sky_search/pycbc_fit_sngls_by_template b/bin/all_sky_search/pycbc_fit_sngls_by_template index 9fbedd9d1b9..3ac27771f12 100755 --- a/bin/all_sky_search/pycbc_fit_sngls_by_template +++ b/bin/all_sky_search/pycbc_fit_sngls_by_template @@ -18,7 +18,7 @@ import argparse, logging import copy, numpy as np -from pycbc import io, events +from pycbc import io, events, init_logging from pycbc.events import triggers, trigger_fits as trstats from pycbc.events import stat as statsmod from pycbc.types.optparse import MultiDetOptionAction @@ -59,10 +59,8 @@ def get_stat(statname, trigs, threshold): parser = argparse.ArgumentParser(usage="", description="Perform maximum-likelihood fits of single inspiral trigger" " distributions to various functions") - +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action=pycbc.version.Version) -parser.add_argument("-V", "--verbose", action="store_true", - help="Print extra debugging information", default=False) parser.add_argument("--trigger-file", help="Input hdf5 file containing single triggers. " "Required") @@ -124,6 +122,8 @@ statsmod.insert_statistic_option_group(parser, default_ranking_statistic='single_ranking_only') args = parser.parse_args() +init_logging(args.verbose) + args.veto_segment_name = sum(args.veto_segment_name, []) args.veto_file = sum(args.veto_file, []) @@ -135,12 +135,6 @@ if (args.prune_param or args.prune_bins or args.prune_number) and not \ raise RuntimeError("To prune, need to specify param, number of bins and " "nonzero number to prune in each bin!") -if args.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -logging.basicConfig(format='%(asctime)s : %(message)s', level=log_level) - logging.info('Fitting above threshold %f' % args.stat_threshold) logging.info('Opening trigger file: %s' % args.trigger_file) diff --git a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam index 06a7ec67220..d4b2dac0f76 100755 --- a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam +++ b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam @@ -16,7 +16,7 @@ import sys, h5py, argparse, logging, pycbc.version, numpy from scipy.stats import norm from pycbc.events import triggers - +from pycbc import init_logging def dist(i1, i2, parvals, smoothing_width): """ @@ -172,9 +172,8 @@ parser = argparse.ArgumentParser(usage="", "parameter, to suppress random noise in the resulting " "background model.") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action=pycbc.version.Version) -parser.add_argument("-V", "--verbose", action="store_true", - help="Print extra debugging information", default=False) parser.add_argument("--template-fit-file", help="hdf5 file containing fit coefficients for each" " individual template. Required") @@ -250,7 +249,7 @@ for inputstr in smooth_kwargs: assert len(args.log_param) == len(args.fit_param) == len(args.smoothing_width) -pycbc.init_logging(args.verbose) +init_logging(args.verbose) fits = h5py.File(args.template_fit_file, 'r') diff --git a/bin/all_sky_search/pycbc_fit_sngls_over_param b/bin/all_sky_search/pycbc_fit_sngls_over_param index 9dc9c3bfd30..8fe805f0221 100644 --- a/bin/all_sky_search/pycbc_fit_sngls_over_param +++ b/bin/all_sky_search/pycbc_fit_sngls_over_param @@ -18,6 +18,7 @@ import argparse, logging import numpy as np +from pycbc import init_logging from pycbc.events import triggers import pycbc.version @@ -27,9 +28,8 @@ parser = argparse.ArgumentParser(usage="", "parameter, to suppress random noise in the resulting " "background model.") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action=pycbc.version.Version) -parser.add_argument("-V", "--verbose", action="store_true", - help="Print extra debugging information", default=False) parser.add_argument("--template-fit-file", help="Input hdf5 file containing fit coefficients for each" " individual template. Required") @@ -78,11 +78,7 @@ args = parser.parse_args() if args.regression_method == 'nn' and args.num_neighbors < 1: raise RuntimeError("Need to give a positive number of nearest neighbors!") -if args.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -logging.basicConfig(format='%(asctime)s : %(message)s', level=log_level) +init_logging(args.verbose) fits = h5py.File(args.template_fit_file, 'r') diff --git a/bin/all_sky_search/pycbc_fit_sngls_split_binned b/bin/all_sky_search/pycbc_fit_sngls_split_binned index 961ec6f5892..37a384d761b 100644 --- a/bin/all_sky_search/pycbc_fit_sngls_split_binned +++ b/bin/all_sky_search/pycbc_fit_sngls_split_binned @@ -31,6 +31,7 @@ import pycbc.version parser = argparse.ArgumentParser(usage="", description="Plot histograms of triggers split over various parameters") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--trigger-file", required=True, help="Input hdf5 file containing single triggers. " "Required") @@ -38,7 +39,6 @@ parser.add_argument("--bank-file", default=None, required=True, help="hdf file containing template parameters. Required") 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', help="Parameter for binning within plots, default " "'template_duration'") diff --git a/bin/all_sky_search/pycbc_followup_file b/bin/all_sky_search/pycbc_followup_file index 559e711a582..88ab2c45f3d 100644 --- a/bin/all_sky_search/pycbc_followup_file +++ b/bin/all_sky_search/pycbc_followup_file @@ -5,9 +5,9 @@ follow-up. import h5py, numpy, argparse, logging, pycbc parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--statmap-file', help="Statmap file containing the candidates/background to follow up") parser.add_argument('--bank-file', diff --git a/bin/all_sky_search/pycbc_foreground_censor b/bin/all_sky_search/pycbc_foreground_censor index dfc7b334a19..546cd14adf2 100755 --- a/bin/all_sky_search/pycbc_foreground_censor +++ b/bin/all_sky_search/pycbc_foreground_censor @@ -7,8 +7,8 @@ import pycbc.events from pycbc.workflow import SegFile parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--foreground-triggers', help="HDF file containing the zerolag foreground triggers " "from the analysis") diff --git a/bin/all_sky_search/pycbc_get_loudest_params b/bin/all_sky_search/pycbc_get_loudest_params index 7b440cdcf5c..14fef61e945 100644 --- a/bin/all_sky_search/pycbc_get_loudest_params +++ b/bin/all_sky_search/pycbc_get_loudest_params @@ -9,12 +9,14 @@ import numpy as np import h5py import argparse import logging +from pycbc import init_logging import pycbc.events from pycbc.pnutils import mass1_mass2_to_mchirp_eta logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO) parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--single-ifo-trigs', type=str, required=True, help='HDF file containing single IFO CBC triggers') parser.add_argument('--tmpltbank-file', type=str, required=True, @@ -33,6 +35,8 @@ parser.add_argument('--print-params', action='store_true', required=False, help='Toggle printing parameters to stdout') args = parser.parse_args() +init_logging(args.verbose) + logging.info('Reading in HDF files') trigs = h5py.File(args.single_ifo_trigs,'r') template_file = h5py.File(args.tmpltbank_file,'r') diff --git a/bin/all_sky_search/pycbc_make_bayestar_skymap b/bin/all_sky_search/pycbc_make_bayestar_skymap new file mode 100644 index 00000000000..7c4d7c3f255 --- /dev/null +++ b/bin/all_sky_search/pycbc_make_bayestar_skymap @@ -0,0 +1,85 @@ +#!/usr/bin/env python + +""" +Simple wrapper around bayestar-localize-coincs to get sky localization for a +specific compact binary merger event. + +Uses an XML containing SNR timeseries and PSD and calls BAYESTAR to produce a +sky localization. + +The wrapper needs to exist as the XML format doesnt allow us to indicate +the low frequency cutoff or waveform used. We also want to be able to control +the output filename which is not possible in bayestar-localize-coincs + +Unrecognised options will be passed straight to the bayestar subprocess +""" + +import argparse +import subprocess +import shutil +import logging +import os +import glob +import sys +import random +import tempfile +import pycbc.version +from pycbc import init_logging + +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument('--version', action="version", + version=pycbc.version.git_verbose_msg) +parser.add_argument('--verbose', action='count', + help="Make the logging increasingly more verbose") +parser.add_argument('--bayestar-executable', + help="The bayestar-localize-coinc executable to be run. " + "If not given, will use whatever is available in " + "the current environment.") +parser.add_argument('--event-xml', required=True, + help="XML file containing event information, SNR " + "timeseries and PSD to pass to bayestar") +parser.add_argument('--waveform', default='TaylorF2', + help="Waveform used in the matched-filtering " + "to generate the SNR timeseries.") +parser.add_argument('--low-frequency-cutoff', type=float, default=20, + help="Low-frequency cutoff used in the matched-filtering " + "to generate the SNR timeseries") +parser.add_argument('--output-file', required=True, + help="Filename to output the fits file to.") +args, unknown = parser.parse_known_args() + +# Default logging is set higher than normal for this job +logging_level = args.verbose + 1 if args.verbose else None +init_logging(logging_level) +logging.info("Starting") + +bayestar_exe = args.bayestar_executable or 'bayestar-localize-coincs' + +tmpdir = tempfile.mkdtemp() + +# Set up the command to pass to bayestar. +# The XML file path being passed twice is a legacy requirement, not a mistake. +cmd = [bayestar_exe, + args.event_xml, + args.event_xml, + '--waveform', args.waveform, + '--f-low', str(args.low_frequency_cutoff), + '-o', tmpdir] + +# Pass any unknown options straight to the subprocess +cmd += unknown + +logging.info("Running %s", ' '.join(cmd)) +subprocess.check_output(cmd) + +# Find the fits file in the temporary directory: +# It would be nice to do this better! - maybe use the input xml +# file to find the number which is going to be used? +fits_filenames = glob.glob(os.path.join(tmpdir, '*.fits*')) +if len(fits_filenames) != 1: + raise ValueError(f'argh, got {len(fits_filenames)} FITS files after running BAYESTAR, I want one!') + +logging.info("Moving output to %s", args.output_file) +shutil.move(fits_filenames[0], args.output_file) + +logging.info("Done") diff --git a/bin/all_sky_search/pycbc_merge_psds b/bin/all_sky_search/pycbc_merge_psds index 2c6f2690e94..7e1eda05070 100755 --- a/bin/all_sky_search/pycbc_merge_psds +++ b/bin/all_sky_search/pycbc_merge_psds @@ -21,8 +21,8 @@ import logging, argparse, numpy, h5py, pycbc.types from pycbc.version import git_verbose_msg as version parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=version) -parser.add_argument('--verbose', action="store_true") parser.add_argument('--psd-files', nargs='+') parser.add_argument("--output-file", required=True) diff --git a/bin/all_sky_search/pycbc_plot_kde_vals b/bin/all_sky_search/pycbc_plot_kde_vals index 54a897041ec..ae107580c5a 100644 --- a/bin/all_sky_search/pycbc_plot_kde_vals +++ b/bin/all_sky_search/pycbc_plot_kde_vals @@ -3,8 +3,11 @@ import numpy, h5py, argparse, logging import matplotlib.pyplot as plt from matplotlib.colors import LogNorm +from pycbc import init_logging, add_common_pycbc_options +import logging parser = argparse.ArgumentParser(description=__doc__) +add_common_pycbc_options(parser) parser.add_argument('--signal-file') parser.add_argument('--template-file', required=True) parser.add_argument('--param', nargs='+', required=True, @@ -17,9 +20,10 @@ parser.add_argument('--log-axis', nargs='+', choices=['True', 'False'], required parser.add_argument('--plot-type', choices=['kde_vs_param', 'param_vs_param']) parser.add_argument('--which-kde', choices=['signal_kde', 'template_kde', 'ratio_kde']) parser.add_argument('--plot-dir', required=True) -parser.add_argument('--verbose', action='count') args = parser.parse_args() +init_logging(args.verbose) + if args.plot_type == 'kde_vs_param': if len(args.param) != 1: parser.error('For kde_vs_param, give exactly one parameter name') diff --git a/bin/all_sky_search/pycbc_prepare_xml_for_gracedb b/bin/all_sky_search/pycbc_prepare_xml_for_gracedb new file mode 100755 index 00000000000..0baf6cb19dc --- /dev/null +++ b/bin/all_sky_search/pycbc_prepare_xml_for_gracedb @@ -0,0 +1,238 @@ +#!/usr/bin/env python + +# Copyright (C) 2015-2023 Ian Harry, Gareth Cabourn Davies +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +""" +Take a coinc xml file containing multiple events and prepare one of them +for upload to gracedb. +""" + +import argparse +import logging +import numpy as np +import h5py +import matplotlib +matplotlib.use('agg') + +import pycbc +import lal +import lal.series +from ligo.lw import lsctables +from ligo.lw import utils as ligolw_utils +from ligo.segments import segment, segmentlist +from pycbc.io.ligolw import ( + LIGOLWContentHandler, + snr_series_to_xml, +) +from pycbc.psd import interpolate +from pycbc.types import FrequencySeries, TimeSeries, load_timeseries +from pycbc.types import MultiDetOptionAction +from pycbc.results import generate_asd_plot, generate_snr_plot + +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument("--verbose", action='count', + help="Increase logging level, default=info") +parser.add_argument("--psd-files", nargs='+', required=True, + help='HDF file(s) containing the PSDs to upload') +parser.add_argument("--snr-timeseries", nargs='+', required=True, + help='HDF file(s) containing the SNR timeseries to upload') +parser.add_argument('--input-file', required=True, type=str, + help='Input LIGOLW XML file of coincidences.') +parser.add_argument('--event-id', type=int, required=True, + help='ID for the event being prepared.') +parser.add_argument('--output-file', + help="Output filename for the locally stored XML. ") +parser.add_argument('--snr-timeseries-plot', + help="Output filename for the plot of SNR timeseries. " + "Must be a plot filetype as used by matplotlib.") +parser.add_argument('--psd-plot', + help="Output filename for the plot of the PSD. " + "Must be a plot filetype as used by matplotlib.") +parser.add_argument('--channel-name', action=MultiDetOptionAction, + required=True, + help="Channel name for adding to the uploaded single " + "inspiral table.") +parser.add_argument('--delta-f', type=float, default=0.25, + help="Frequency spacing of the PSD to be uploaded, Hz. " + "Default 0.25") + +args = parser.parse_args() + +if args.verbose: + args.verbose += 1 +else: + args.verbose = 1 +pycbc.init_logging(args.verbose) + +xmldoc = ligolw_utils.load_filename(args.input_file, + contenthandler=LIGOLWContentHandler) + +class psd_segment(segment): + def __new__(cls, psd, *args): + return segment.__new__(cls, *args) + def __init__(self, psd, *args): + self.psd = psd + +def read_psds(psd_files): + logging.info("Reading PSDs") + psds = {} + + for psd_file in psd_files: + (ifo, group), = h5py.File(psd_file, "r").items() + psd = [group["psds"][str(i)] for i in range(len(group["psds"].keys()))] + psds[ifo] = segmentlist(psd_segment(*segargs) for segargs in zip( + psd, group["start_time"], group["end_time"])) + return psds + +psds = read_psds(args.psd_files) + +def read_snr_timeseries(snr_timeseries_files): + """ + Get information from the single_template snr timeseries files + """ + logging.info("Reading SNR timeseries") + snr_timeseries = {} + for snr_filename in snr_timeseries_files: + with h5py.File(snr_filename, 'r') as snr_f: + ifo = snr_f.attrs['ifo'] + time = snr_f.attrs['event_time'] + snr_timeseries[(ifo, time)] = load_timeseries(snr_filename, 'snr') + return snr_timeseries + +snr_timeseries = read_snr_timeseries(args.snr_timeseries) + +coinc_table = lsctables.CoincTable.get_table(xmldoc) +coinc_inspiral_table = lsctables.CoincInspiralTable.get_table(xmldoc) +coinc_event_map_table = lsctables.CoincMapTable.get_table(xmldoc) +sngl_inspiral_table = lsctables.SnglInspiralTable.get_table(xmldoc) + +xmldoc.childNodes[-1].removeChild(sngl_inspiral_table) +xmldoc.childNodes[-1].removeChild(coinc_event_map_table) +xmldoc.childNodes[-1].removeChild(coinc_inspiral_table) +xmldoc.childNodes[-1].removeChild(coinc_table) + +# Find the coinc_table entry corresponding to the desired event ID +for event in coinc_table: + coinc_event_id = event.coinc_event_id + if args.event_id == coinc_event_id: + # This is the event we want + break +else: + raise ValueError(f"event {args.event_id} not found in coinc table") + +coinc_event_table_curr = lsctables.New(lsctables.CoincTable) +coinc_event_table_curr.append(event) +coinc_inspiral_table_curr = lsctables.New(lsctables.CoincInspiralTable) +coinc_event_map_table_curr = lsctables.New(lsctables.CoincMapTable) +sngl_inspiral_table_curr = lsctables.New(lsctables.SnglInspiralTable) + +for coinc_insp in coinc_inspiral_table: + if coinc_insp.coinc_event_id == event.coinc_event_id: + coinc_inspiral_table_curr.append(coinc_insp) + +time = coinc_inspiral_table_curr[0].end_time \ + + coinc_inspiral_table_curr[0].end_time_ns * 1e-9 + +sngl_ids = [] +for coinc_map in coinc_event_map_table: + if coinc_map.coinc_event_id == event.coinc_event_id: + coinc_event_map_table_curr.append(coinc_map) + sngl_ids.append(coinc_map.event_id) + +# Get the SNR timeseries and PSDs at the time of this event +# these may not have contributed a trigger + +# Check that the SNR timeseries correspond to the event +# - allow 0.5s either side of the time +snr_ts = {} +for ifo, t in snr_timeseries.keys(): + if not abs(t - time) < 0.5: + raise ValueError("SNR timeseries for IFO %s does not look like it " + "corresponds to this event, event time %.3f, SNR timeseries is " + "around time %.3f" % (ifo, time, t)) + # Convert the SNR dict to be keyed on IFO-only: + snr_ts[ifo] = snr_timeseries[(ifo, t)] + +# IFOs from SNR timeseries: +psds_event = {} +psddict = {} +for ifo in snr_ts.keys(): + psd = psds[ifo] + + psd = psd[psd.find(time)].psd + # resample the psd to new spacing + psd_fs = FrequencySeries(psd, delta_f=psd.attrs["delta_f"], + dtype=np.float64) + psd_fs = interpolate(psd_fs, args.delta_f) + + psds_event[ifo] = psd + psddict[ifo] = psd_fs + +lal_psddict = {} +sample_freqs = {} + +for sngl in sngl_inspiral_table: + if sngl.event_id not in sngl_ids: + continue + sngl_inspiral_table_curr.append(sngl) + sngl.channel = args.channel_name[sngl.ifo] + + # Two versions of the PSD dictionary have useful info here + psd = psds_event[sngl.ifo] + psd_fs = psddict[sngl.ifo] + + flow = psd.file.attrs['low_frequency_cutoff'] + kmin = int(flow / args.delta_f) + + fseries = lal.CreateREAL8FrequencySeries( + "psd", + lal.LIGOTimeGPS(int(psd.attrs["epoch"])), + kmin * args.delta_f, + args.delta_f, + lal.StrainUnit**2 / lal.HertzUnit, + len(psd_fs) - kmin) + fseries.data.data = psd_fs[kmin:] / np.square(pycbc.DYN_RANGE_FAC) + lal_psddict[sngl.ifo] = fseries + + snr_series_to_xml(snr_ts[sngl.ifo], xmldoc, sngl.event_id) + +xmldoc.childNodes[-1].appendChild(coinc_event_table_curr) +xmldoc.childNodes[-1].appendChild(coinc_inspiral_table_curr) +xmldoc.childNodes[-1].appendChild(coinc_event_map_table_curr) +xmldoc.childNodes[-1].appendChild(sngl_inspiral_table_curr) +lal.series.make_psd_xmldoc(lal_psddict, xmldoc.childNodes[-1]) + +ligolw_utils.write_filename(xmldoc, args.output_file) +logging.info("Saved XML file %s", args.output_file) + +if args.psd_plot: + logging.info("Saving ASD plot %s", args.psd_plot) + generate_asd_plot(psddict, args.psd_plot) + +if args.snr_timeseries_plot: + triggers = {sngl.ifo: (sngl.end_time + sngl.end_time_ns * 1e-9, sngl.snr) + for sngl in sngl_inspiral_table_curr} + base_time = int(np.floor(time)) + logging.info("Saving SNR plot %s", args.snr_timeseries_plot) + generate_snr_plot( + snr_ts, + args.snr_timeseries_plot, + triggers, + base_time + ) + +logging.info('Done!') diff --git a/bin/all_sky_search/pycbc_reduce_template_bank b/bin/all_sky_search/pycbc_reduce_template_bank index 21b07a00c13..99865b39566 100644 --- a/bin/all_sky_search/pycbc_reduce_template_bank +++ b/bin/all_sky_search/pycbc_reduce_template_bank @@ -30,10 +30,9 @@ import pycbc import pycbc.version parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) -parser.add_argument("-V", "--verbose", action="store_true", - help="print extra debugging information", default=False) parser.add_argument("--input-bank", required=True, help="Input template bank HDF file.") parser.add_argument("--output-bank", required=True, diff --git a/bin/all_sky_search/pycbc_rerank_dq b/bin/all_sky_search/pycbc_rerank_dq deleted file mode 100644 index 1eafbdc79cb..00000000000 --- a/bin/all_sky_search/pycbc_rerank_dq +++ /dev/null @@ -1,121 +0,0 @@ -#!/usr/bin/env python -""" Merge and rerank clustered DQ data -""" -import logging, argparse, numpy, h5py -import pycbc -from pycbc.types.timeseries import load_timeseries -from os import path -from pycbc.version import git_verbose_msg as version - -parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument('--version', action='version', version=version) -parser.add_argument('--verbose', action="store_true",default=False) -parser.add_argument('--no-low-dq-cutoff', action="store_true", - help="Allow dq values less than 0. " - "Otherwise values less than 0 are set to 0") -parser.add_argument("--ifo", type=str,required=True, - help="interferometer to analyze") -parser.add_argument("--max-stat", type=float,default=8, - help="maximum allowed dq value. " - "Only used if --rate-file is not used. " - "[default=8]") -parser.add_argument("--input-file", required=True,nargs='+', - help="List of files containing the dq timeseries") -parser.add_argument("--input-channel", required=False,type=str, - help='name of channel to read in from ' - 'provided input file. Required if ' - 'provided file is .hdf5 format') -rate_group = parser.add_mutually_exclusive_group() -rate_group.add_argument("--rate-file", type = str, default = None, - help="File containing the mapping from each " - "dq bin to the trigger rate") -rate_group.add_argument("--bank-file", type = str, default = None, - help="File containing the template bank") -parser.add_argument("--dq-type", type = str, default = 'dq', - help="Name of dq type to be included in output file " - "metadata [default='dq']") -parser.add_argument("--output-file", required=True, - help="name of output file") - -args = parser.parse_args() -pycbc.init_logging(args.verbose) - -ifo = args.ifo - -log_like = numpy.array([]) -times = numpy.array([]) -for filename in args.input_file: - logging.info('Reading file %s...', filename) - dq_data = load_timeseries(filename, group=args.input_channel) - log_like = numpy.concatenate((log_like,dq_data[:])) - times = numpy.concatenate((times,dq_data.sample_times)) - del dq_data - -log_like_ranking = numpy.argsort(log_like) -len_log_like = len(log_like) - -times = times[log_like_ranking] -log_like_vals = numpy.array(log_like[log_like_ranking]) - -log_like_dict = {} -percent_dict = {} - -if args.rate_file is not None: - rates = {} - locs_dict = {} - with h5py.File(args.rate_file,'r') as f: - bin_names = f.attrs['names'][:] - for bin_name in bin_names: - logging.info('Processing bin %s...', bin_name) - rates = f['%s/rates'%bin_name][:] - locs_dict[bin_name] = f['%s/locs'%bin_name][:] - log_like_dict[bin_name] = numpy.zeros(len(log_like)) - - if not args.no_low_dq_cutoff: - rates = numpy.array([max(r,1.) for r in rates]) - - log_rates = numpy.log(rates) - percent_dict[bin_name] = log_rates - - n_bins = len(log_rates) - percentiles = numpy.linspace(0,100,n_bins+1) - dq_percentiles = numpy.percentile(log_like_vals,percentiles)[1:] - - rates_bin = numpy.array([n_bins - \ - len(dq_percentiles[dq_percentiles >= dq_ll]) \ - for dq_ll in log_like_vals]) - log_like_dict[bin_name] = log_rates[rates_bin.astype('int')] - -else: - # Use analytic model from Godwin et al. arXiv:2010.15282 - bin_names = ['all_bin'] - with h5py.File(args.bank_file, 'r') as bank: - locs_dict = {'all_bin': numpy.arange(0, len(bank['mass1'][:]), 1)} - - log_like = numpy.arange(len(log_like))*100./len(log_like) - - p_min = int(len_log_like/2.) - p_max = int(len_log_like/(1.+numpy.exp(-args.max_stat))) - - log_like[p_min:p_max] = numpy.log(log_like[p_min:p_max] \ - / (100.-log_like[p_min:p_max])) - log_like[p_max:] = args.max_stat - log_like[:p_min] = 0 - log_like_dict['all_bin'] = log_like - - -with h5py.File(args.output_file, 'w') as f: - f.attrs['names'] = bin_names - for bin_name in bin_names: - f[ifo + '/dq_vals/' + bin_name] = numpy.array(log_like_dict[bin_name], - dtype=numpy.float32) - f[ifo + '/dq_percentiles/' + bin_name] = numpy.array(percent_dict[bin_name], - dtype=numpy.float32) - f[ifo + '/locs/' + bin_name] = numpy.array(locs_dict[bin_name], - dtype=numpy.float32) - f[ifo + '/times'] = numpy.array(times, dtype=numpy.uint32) - f[ifo + '/dq_input_vals'] = numpy.array(log_like_vals, dtype=numpy.float32) - f.attrs.create('stat', data=ifo+'-'+args.dq_type+'-dq_ts_reference') - -logging.info('Done!') - diff --git a/bin/all_sky_search/pycbc_rerank_passthrough b/bin/all_sky_search/pycbc_rerank_passthrough index 6c0a4b274b1..acc0340b08c 100644 --- a/bin/all_sky_search/pycbc_rerank_passthrough +++ b/bin/all_sky_search/pycbc_rerank_passthrough @@ -3,10 +3,9 @@ import h5py, numpy, argparse, logging, pycbc parser = argparse.ArgumentParser() - +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--output-file', help="File containing the newly assigned statistic values") diff --git a/bin/all_sky_search/pycbc_sngls_findtrigs b/bin/all_sky_search/pycbc_sngls_findtrigs index c779608448b..4dd447182e7 100644 --- a/bin/all_sky_search/pycbc_sngls_findtrigs +++ b/bin/all_sky_search/pycbc_sngls_findtrigs @@ -12,7 +12,7 @@ from pycbc.types.optparse import MultiDetOptionAction from pycbc import init_logging parser = argparse.ArgumentParser() -parser.add_argument("--verbose", action='count') +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action='version', version=pycbc.version.git_verbose_msg) # Basic file input options @@ -139,10 +139,10 @@ for tnum in template_ids: trigger_keep_ids = cuts.apply_trigger_cuts(trigs, trigger_cut_dict, statistic=rank_method) tids_full = tids_uncut[trigger_keep_ids] - logging.info('%s:%s', tnum, len(tids_uncut)) + logging.debug('%s:%s', tnum, len(tids_uncut)) if len(tids_full) < len(tids_uncut): - logging.info("%s triggers cut", - len(tids_uncut) - len(tids_full)) + logging.debug("%s triggers cut", + len(tids_uncut) - len(tids_full)) n_tot_trigs = tids_full.size if not n_tot_trigs: continue diff --git a/bin/all_sky_search/pycbc_sngls_pastro b/bin/all_sky_search/pycbc_sngls_pastro index 34a093d4f80..ccaa252ca67 100644 --- a/bin/all_sky_search/pycbc_sngls_pastro +++ b/bin/all_sky_search/pycbc_sngls_pastro @@ -39,7 +39,7 @@ mchirp_power = { } parser = argparse.ArgumentParser() -parser.add_argument("--verbose", action='count') +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action='version', version=pycbc.version.git_verbose_msg) parser.add_argument("--single-statmap-files", nargs='+', required=True, diff --git a/bin/all_sky_search/pycbc_sngls_statmap b/bin/all_sky_search/pycbc_sngls_statmap index e3c2f03ebdb..5ebd6b39303 100755 --- a/bin/all_sky_search/pycbc_sngls_statmap +++ b/bin/all_sky_search/pycbc_sngls_statmap @@ -34,6 +34,7 @@ class fw(object): parser = argparse.ArgumentParser() # General required options +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) parser.add_argument('--sngls-files', nargs='+', @@ -41,7 +42,6 @@ parser.add_argument('--sngls-files', nargs='+', 'information.') parser.add_argument('--ifos', nargs=1, help='List of ifos used in these coincidence files') -parser.add_argument('--verbose', action='count') parser.add_argument('--cluster-window', type=float, default=10, help='Length of time window in seconds to cluster coinc ' 'events [default=10s]') diff --git a/bin/all_sky_search/pycbc_sngls_statmap_inj b/bin/all_sky_search/pycbc_sngls_statmap_inj index c7ba970d7c4..7d164c9383d 100644 --- a/bin/all_sky_search/pycbc_sngls_statmap_inj +++ b/bin/all_sky_search/pycbc_sngls_statmap_inj @@ -34,6 +34,7 @@ class fw(object): parser = argparse.ArgumentParser() # General required options +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) parser.add_argument('--sngls-files', nargs='+', @@ -44,7 +45,6 @@ parser.add_argument('--full-data-background', required=True, 'injection coincs') parser.add_argument('--ifos', nargs=1, help='List of ifos used in these coincidence files') -parser.add_argument('--verbose', action='count') parser.add_argument('--cluster-window', type=float, default=10, help='Length of time window in seconds to cluster coinc ' 'events [default=10s]') diff --git a/bin/all_sky_search/pycbc_strip_injections b/bin/all_sky_search/pycbc_strip_injections index f77d458f5df..0ee3840ad97 100644 --- a/bin/all_sky_search/pycbc_strip_injections +++ b/bin/all_sky_search/pycbc_strip_injections @@ -12,8 +12,8 @@ def remove(l, i): l.remove(r) parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--injection-file') parser.add_argument('--veto-file', help="File containing segments used to veto injections") diff --git a/bin/all_sky_search/pycbc_template_kde_calc b/bin/all_sky_search/pycbc_template_kde_calc index 456e9b0acdb..ae4fe3b8315 100644 --- a/bin/all_sky_search/pycbc_template_kde_calc +++ b/bin/all_sky_search/pycbc_template_kde_calc @@ -13,15 +13,15 @@ # Public License for more details. import numpy, h5py, operator, argparse, logging -from pycbc import init_logging +from pycbc import init_logging, add_common_pycbc_options import pycbc.conversions as convert from pycbc import libutils from pycbc.events import triggers akde = libutils.import_optional('awkde') kf = libutils.import_optional('sklearn.model_selection') - parser = argparse.ArgumentParser(description=__doc__) +add_common_pycbc_options(parser) parser.add_argument('--signal-file', help='File with parameters of GW signals ' 'for KDE calculation') parser.add_argument('--template-file', required=True, help='Hdf5 file with ' @@ -64,9 +64,8 @@ parser.add_argument('--mchirp-downsample-power', type=float, help='Exponent value for the power law distribution') parser.add_argument('--min-ratio', type=float, help='Minimum ratio for template_kde relative to the maximum') -parser.add_argument('--verbose', action='store_true') args = parser.parse_args() -init_logging(verbose=args.verbose, format='%(asctime)s %(message)s') +init_logging(args.verbose) assert len(args.fit_param) == len(args.log_param) diff --git a/bin/all_sky_search/pycbc_template_kde_max b/bin/all_sky_search/pycbc_template_kde_max index 3474a401d00..31172746847 100644 --- a/bin/all_sky_search/pycbc_template_kde_max +++ b/bin/all_sky_search/pycbc_template_kde_max @@ -1,17 +1,17 @@ #!/usr/bin/env python import numpy, h5py, argparse, logging -from pycbc import init_logging +from pycbc import init_logging, add_common_pycbc_options parser = argparse.ArgumentParser(description=__doc__) +add_common_pycbc_options(parser) parser.add_argument('--kde-files', nargs='+', required=True, help='HDF files with KDE values') parser.add_argument('--output-file', required=True, help='Name of output HDF file') parser.add_argument('--min-ratio', type=float, help='Minimum ratio for template_kde relative to the maximum') -parser.add_argument('--verbose', action='store_true') args = parser.parse_args() -init_logging(verbose=args.verbose, format='%(asctime)s %(message)s') +init_logging(args.verbose) input_kdes = [h5py.File(kfile, 'r') for kfile in args.kde_files] diff --git a/bin/all_sky_search/pycbc_template_recovery_hist b/bin/all_sky_search/pycbc_template_recovery_hist index ed2900826a9..5ae76eb214c 100644 --- a/bin/all_sky_search/pycbc_template_recovery_hist +++ b/bin/all_sky_search/pycbc_template_recovery_hist @@ -6,9 +6,10 @@ from matplotlib import use use('Agg') from matplotlib import pyplot as plt from pycbc.events import triggers +from pycbc import init_logging, add_common_pycbc_options parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument('--verbose', action='count') +add_common_pycbc_options(parser) parser.add_argument('--output', required=True) parser.add_argument('--found-injection-files', dest='found', nargs='+', help='hdf file(s) with found injections') @@ -29,9 +30,7 @@ parser.add_argument('--min-ifar', type=float, help='only plot injections above given ifar value') args = parser.parse_args() -if args.verbose: - log_level = logging.INFO - logging.basicConfig(format='%(asctime)s : %(message)s', level=log_level) +init_logging(args.verbose) # should be same number of inj and trig files # and either 0 or 1 bank files or 1 per inj file diff --git a/bin/all_sky_search/pycbc_upload_single_event_to_gracedb b/bin/all_sky_search/pycbc_upload_single_event_to_gracedb new file mode 100755 index 00000000000..78903ebbceb --- /dev/null +++ b/bin/all_sky_search/pycbc_upload_single_event_to_gracedb @@ -0,0 +1,182 @@ +#!/usr/bin/env python + +# Copyright (C) 2015-2023 Ian Harry, Gareth Cabourn Davies +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +""" +Take a coinc xml file containing a single event and upload to gracedb. +""" + +import os +import sys +import argparse +import logging +from ligo.gracedb.rest import GraceDb +import numpy as np + +import pycbc +import lal +from pycbc.io.live import gracedb_tag_with_version + +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument('--xml-file-for-upload', required=True, type=str, + help='LIGOLW XML file containing the event.') +parser.add_argument('--log-message', type=str, metavar='MESSAGE', + help='Add a log entry to each upload with the given message') +parser.add_argument('--testing', action="store_true", default=False, + help="Upload event to the TEST group of gracedb.") +parser.add_argument('--production-server', action="store_true", default=False, + help="Upload event to production graceDB. If not given " + "events will be uploaded to playground server.") +parser.add_argument('--search-tag', default='AllSky', + help="Specify the search tag. Default: AllSky") +parser.add_argument('--snr-timeseries-plot', + help="SNR timeseries plot to be uploaded with the event.") +parser.add_argument('--asd-plot', + help="ASD plot to be uploaded with the event.") +parser.add_argument('--skymap-plot', + help="Skymap plot to be uploaded with the event.") +parser.add_argument('--skymap-fits-file', + help="Skymap .fits file to be uploaded with the event.") +parser.add_argument('--source-probabilities', + help="Source probabilities file to be uploaded with the " + "event.") +parser.add_argument('--source-probabilities-plot', + help="Source probabilities plot to be uploaded with the " + "event.") +parser.add_argument('--labels', nargs='+', + help="Labels to add to the event in GraceDB") + +args = parser.parse_args() + +pycbc.init_logging(logging.INFO) +# Make the scitokens logger a little quieter +# (it is called through GraceDB) +logging.getLogger('scitokens').setLevel(logging.root.level + 10) + +if args.production_server: + gracedb = GraceDb() +else: + gracedb = GraceDb(service_url='https://gracedb-playground.ligo.org/api/') + +labels = [l.upper() for l in args.labels] +allowed_labels = gracedb.allowed_labels + +if set(labels) - set(allowed_labels): + err_msg = "One or more supplied labels is not available on the server. " + err_msg += f"Supplied {','.join(labels)}, allowed " + err_msg += f"{','.join(allowed_labels)}." + raise RuntimeError(err_msg) + +group_tag = 'Test' if args.testing else 'CBC' +r = gracedb.create_event( + group_tag, + 'pycbc', + args.xml_file_for_upload, + filecontents=open(args.xml_file_for_upload, "rb").read(), + search=args.search_tag, + offline=True, + labels=labels +).json() + +logging.info("Uploaded event %s.", r["graceid"]) + +# add info for tracking code version +gracedb_tag_with_version(gracedb, r['graceid']) + +# document the absolute path to the input file +input_file_str = 'Candidate uploaded from ' \ + + os.path.abspath(args.xml_file_for_upload) +gracedb.write_log(r['graceid'], input_file_str) + +# document the command line used in the event log +log_str = 'Upload command: ' + ' '.join(sys.argv) + +# add the custom log message, if provided +if args.log_message is not None: + log_str += '. ' + args.log_message + +gracedb.write_log( + r['graceid'], + log_str, + tag_name=['analyst_comments'] +) + + +def upload_file(upload_filename, displayname, comment, tag): + """ + Helper function to upload files associated with the event. + """ + logging.info("Uploading %s file %s to event %s.", + displayname, upload_filename, r["graceid"]) + gracedb.write_log( + r["graceid"], + comment, + filename=upload_filename, + tag_name=[tag], + displayName=[displayname] + ) + + +if args.asd_plot: + upload_file( + args.asd_plot, + "ASDs", + "PyCBC ASD estimate from the time of event", + "psd" + ) + +if args.snr_timeseries_plot: + upload_file( + args.snr_timeseries_plot, + "SNR timeseries", + "SNR timeseries plot upload", + "background" + ) + +if args.skymap_plot: + upload_file( + args.skymap_plot, + "", + "Skymap plot upload", + "sky_loc" + ) + +if args.skymap_fits_file: + upload_file( + args.skymap_fits_file, + "", + "sky localization complete", + "sky_loc" + ) + +if args.source_probabilities: + upload_file( + args.source_probabilities, + "", + "Source probabilities JSON file upload", + "pe" + ) + +if args.source_probabilities_plot: + upload_file( + args.source_probabilities_plot, + "", + "Source probabilities plot upload", + "pe" + ) + +logging.info('Done!') diff --git a/bin/minifollowups/pycbc_injection_minifollowup b/bin/minifollowups/pycbc_injection_minifollowup index d16e7c1cd71..198ca2a34bf 100644 --- a/bin/minifollowups/pycbc_injection_minifollowup +++ b/bin/minifollowups/pycbc_injection_minifollowup @@ -14,20 +14,112 @@ # You should have received a copy of the GNU General Public License along # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -""" Followup foreground events """ -import os, argparse, numpy, logging, h5py, copy +"""Create a workflow for following up missed loud injections.""" + +import os +import argparse +import logging +import copy +import numpy +import h5py import pycbc.workflow as wf +import pycbc.workflow.minifollowups as mini +import pycbc.version from pycbc.types import MultiDetOptionAction from pycbc.events import select_segments_by_definer, coinc from pycbc.results import layout from pycbc.detector import Detector -import pycbc.workflow.minifollowups as mini from pycbc.workflow.core import resolve_url_to_file -import pycbc.version, pycbc.pnutils from pycbc.io.hdf import SingleDetTriggers, HFile -parser = argparse.ArgumentParser(description=__doc__[1:]) + +legal_distance_types = [ + 'decisive_optimal_snr', + 'comb_optimal_snr', + 'dec_chirp_distance' +] + + +def sort_injections(args, inj_group, missed): + """Return an array of indices to sort the missed injections from most to + least likely to be detected, according to a metric of choice. + + Parameters + ---------- + args : object + CLI arguments parsed by argparse. Must have a `distance_type` attribute + to specify how to sort, which must be one of the values in + `legal_distance_types`. + inj_group : h5py group object + HDF5 group object containing the injection definition. + missed : array + Array of indices of missed injections into `inj_group`. + + Returns + ------- + missed_sorted : array + Array of indices of missed injections sorted as requested. + """ + if not hasattr(args, 'distance_type'): + raise ValueError('Distance type not provided') + if args.distance_type not in legal_distance_types: + raise ValueError( + f'Invalid distance type "{args.distance_type}", ' + f'allowed types are {", ".join(legal_distance_types)}' + ) + + if 'optimal_snr' in args.distance_type: + optimal_snrs = [ + inj_group[dsn][:][missed] for dsn in inj_group.keys() + if dsn.startswith('optimal_snr_') + ] + assert optimal_snrs, 'These injections do not have optimal SNRs' + + if args.distance_type == 'decisive_optimal_snr': + # descending order of decisive (2nd largest) optimal SNR + dec_snr = numpy.array([ + sorted(snrs)[-2] for snrs in zip(*optimal_snrs) + ]) + if args.maximum_decisive_snr is not None: + # By setting to 0, these injections will not be considered + dec_snr[dec_snr > args.maximum_decisive_snr] = 0 + sorter = dec_snr.argsort()[::-1] + return missed[sorter] + + if args.distance_type == 'comb_optimal_snr': + # descending order of network optimal SNR + optimal_snrs = numpy.vstack(optimal_snrs) + net_opt_snrs_squared = (optimal_snrs ** 2).sum(axis=0) + sorter = net_opt_snrs_squared.argsort()[::-1] + return missed[sorter] + + if args.distance_type == 'dec_chirp_distance': + # ascending order of decisive (2nd smallest) chirp distance + from pycbc.conversions import mchirp_from_mass1_mass2, chirp_distance + + eff_dists = [] + for ifo in args.single_detector_triggers: + eff_dist = Detector(ifo).effective_distance( + inj_group['distance'][:][missed], + inj_group['ra'][:][missed], + inj_group['dec'][:][missed], + inj_group['polarization'][:][missed], + inj_group['tc'][:][missed], + inj_group['inclination'][:][missed] + ) + eff_dists.append(eff_dist) + dec_eff_dist = sorted(eff_dists)[-2] + mchirp = mchirp_from_mass1_mass2( + inj_group['mass1'][:][missed], + inj_group['mass2'][:][missed] + ) + dec_chirp_dist = chirp_distance(dec_dist, mchirp) + sorter = dec_chirp_dist.argsort() + return missed[sorter] + + +parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) parser.add_argument('--bank-file', help="HDF format template bank file") @@ -57,6 +149,11 @@ parser.add_argument('--nearby-triggers-window', type=float, default=0.05, help="Maximum time difference between the missed " "injection and the loudest SNR nearby trigger to " "display, seconds. Default=0.05") +parser.add_argument('--distance-type', + required=True, + choices=legal_distance_types, + help="How to sort missed injections from most to least " + "likely to be detected") wf.add_workflow_command_line_group(parser) wf.add_workflow_settings_cli(parser, include_subdax_opts=True) args = parser.parse_args() @@ -84,12 +181,16 @@ for ifo in args.single_detector_triggers: strig_file = resolve_url_to_file(os.path.abspath(fname), attrs={'ifos': ifo}) single_triggers.append(strig_file) - insp_data_seglists[ifo] = select_segments_by_definer\ - (args.inspiral_segments, segment_name=args.inspiral_data_read_name, - ifo=ifo) - insp_analysed_seglists[ifo] = select_segments_by_definer\ - (args.inspiral_segments, segment_name=args.inspiral_data_analyzed_name, - ifo=ifo) + insp_data_seglists[ifo] = select_segments_by_definer( + args.inspiral_segments, + segment_name=args.inspiral_data_read_name, + ifo=ifo + ) + insp_analysed_seglists[ifo] = select_segments_by_definer( + args.inspiral_segments, + segment_name=args.inspiral_data_analyzed_name, + ifo=ifo + ) # NOTE: make_singles_timefreq needs a coalesced set of segments. If this is # being used to determine command-line options for other codes, # please think if that code requires coalesced, or not, segments. @@ -97,56 +198,21 @@ for ifo in args.single_detector_triggers: insp_analysed_seglists[ifo].coalesce() f = h5py.File(args.injection_file, 'r') +inj_def = f['injections'] missed = f['missed/after_vetoes'][:] if args.ifar_threshold is not None: try: # injections may not have (inclusive) IFAR present ifars = f['found_after_vetoes']['ifar'][:] except KeyError: ifars = f['found_after_vetoes']['ifar_exc'][:] - logging.warn('Inclusive IFAR not found, using exclusive') + logging.warning('Inclusive IFAR not found, using exclusive') lgc_arr = ifars < args.ifar_threshold missed = numpy.append(missed, f['found_after_vetoes']['injection_index'][lgc_arr]) -num_events = int(workflow.cp.get_opt_tags('workflow-injection_minifollowups', 'num-events', '')) - -try: - optstrings = [os for os in f['injections'].keys() if \ - os.startswith('optimal_snr_')] - optimal_snr = [f['injections/%s' % os][:][missed] for os in optstrings] - # 2nd largest opt SNR - dec_snr = [sorted(snrs)[-2] for snrs in zip(*optimal_snr)] - dec_snr = numpy.array(dec_snr) - - if args.maximum_decisive_snr is not None: - # By setting to 0, these injections will not be considered - dec_snr[dec_snr > args.maximum_decisive_snr] = 0 - sorting = dec_snr.argsort() - sorting = sorting[::-1] # descending order of dec opt SNR -except KeyError: - # Fall back to effective distance if optimal SNR not available - eff_dist = {} - for trig in single_triggers: - ifo = trig.ifo - eff_dist[ifo] = Detector(ifo).effective_distance( - f['injections/distance'][:][missed], - f['injections/ra'][:][missed], - f['injections/dec'][:][missed], - f['injections/polarization'][:][missed], - f['injections/tc'][:][missed], - f['injections/inclination'][:][missed]) - - dec_dist = numpy.maximum(eff_dist[single_triggers[0].ifo], - eff_dist[single_triggers[1].ifo]) - mchirp, eta = pycbc.pnutils.mass1_mass2_to_mchirp_eta(\ - f['injections/mass1'][:][missed], - f['injections/mass2'][:][missed]) - dec_chirp_dist = pycbc.pnutils.chirp_distance(dec_dist, mchirp) - sorting = dec_chirp_dist.argsort() # ascending order of dec chirp distance - # Get the trigger SNRs and times # But only ones which are within a small window of the missed injection -missed_inj_times = numpy.sort(f['injections/tc'][:][missed]) +missed_inj_times = numpy.sort(inj_def['tc'][:][missed]) # Note: Adding Earth diameter in light seconds to the window here # to allow for different IFO's arrival times of the injection @@ -188,23 +254,32 @@ for trig in single_triggers: f'{ifo}/snr', return_indices=True) +# figure out how many injections to follow up +num_events = int(workflow.cp.get_opt_tags( + 'workflow-injection_minifollowups', + 'num-events', + '' +)) if len(missed) < num_events: num_events = len(missed) -# loop over loudest events to be followed up +# sort the injections +missed = sort_injections(args, inj_def, missed) + +# loop over sorted missed injections to be followed up found_inj_idxes = f['found_after_vetoes/injection_index'][:] for num_event in range(num_events): files = wf.FileList([]) - injection_index = missed[sorting][num_event] - time = f['injections/tc'][injection_index] - lon = f['injections/ra'][injection_index] - lat = f['injections/dec'][injection_index] + injection_index = missed[num_event] + time = inj_def['tc'][injection_index] + lon = inj_def['ra'][injection_index] + lat = inj_def['dec'][injection_index] ifo_times = '' inj_params = {} for val in ['mass1', 'mass2', 'spin1z', 'spin2z', 'tc']: - inj_params[val] = f['injections/%s' %(val,)][injection_index] + inj_params[val] = inj_def[val][injection_index] for single in single_triggers: ifo = single.ifo det = Detector(ifo) @@ -271,9 +346,10 @@ for num_event in range(num_events): tags=args.tags + [str(num_event)]) break else: - logging.info('Trigger time {} is not valid in {}, ' \ - 'skipping singles plots'.format(checkedtime, - single.ifo)) + logging.info( + 'Trigger time %s is not valid in %s, skipping singles plots', + checkedtime, single.ifo + ) files += mini.make_single_template_plots(workflow, insp_segs, args.inspiral_data_read_name, @@ -349,7 +425,6 @@ for num_event in range(num_events): params_str='loudest template in %s' % curr_ifo ) layouts += list(layout.grouper(files, 2)) - num_event += 1 workflow.save() layout.two_column_layout(args.output_dir, layouts) diff --git a/bin/plotting/pycbc_plot_dq_flag_likelihood b/bin/plotting/pycbc_plot_dq_flag_likelihood index 3da0feca016..a1d7beed5c1 100644 --- a/bin/plotting/pycbc_plot_dq_flag_likelihood +++ b/bin/plotting/pycbc_plot_dq_flag_likelihood @@ -6,15 +6,17 @@ import argparse import numpy import pycbc import h5py -from matplotlib import use; use('Agg') +from matplotlib import use as matplotlib_use from matplotlib import pyplot from pycbc.version import git_verbose_msg as version import pycbc.results +matplotlib_use('Agg') parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--version', action='version', version=version) parser.add_argument('--verbose', action="store_true") parser.add_argument("--dq-file", required=True) +parser.add_argument("--dq-label", required=True) parser.add_argument("--ifo", type=str, required=True) parser.add_argument("--output-file", required=True) args = parser.parse_args() @@ -24,45 +26,60 @@ pycbc.init_logging(args.verbose) ifo = args.ifo f = h5py.File(args.dq_file, 'r') -dq_name = f.attrs['stat'].split('-')[1] -logrates_grp = f[ifo + '/dq_percentiles'] -bin_names = logrates_grp.keys() - -# when using a flag, should have exactly 2 dq bins (off and on) -# times when the flag was on go to the last dq bin -for b in bin_names: - num_bins = len(logrates_grp[b]) - assert num_bins == 2, f'Flag should have exactly 2 dq bins, found {num_bins}' -logrates = [logrates_grp[b][-1] for b in bin_names] -x = numpy.arange(len(logrates)) - -ymax = 1.2 * max(logrates) -color = pycbc.results.ifo_color(ifo) - -fig = pyplot.figure(0) -ax = fig.add_subplot(111) +ifo_grp = f[ifo] + +dq_states = ifo_grp['dq_segments'].keys() +bin_names = ifo_grp['bins'].keys() + +num_bins = len(bin_names) +x = numpy.arange(num_bins) +width = 0.25 + +dq_states = { + 0: 'Clean', + 1: 'DQ Flag', + 2: 'Autogating', +} + +colors = { + 0: 'green', + 1: 'gold', + 2: 'red' +} + +fig, ax = pyplot.subplots(figsize=(9, 6)) ax2 = ax.twinx() -ax2.bar(x, height=logrates, label=ifo, color=color) -ax2.set_ylabel('DQ Log Likelihood Penalty') -ax2.legend(loc='upper left', markerscale=5) -ax2.set_ylim(0, ymax) -ax2.grid() +ymax = 1 +ymin = 0.9 +for n, dqstate_name in dq_states.items(): + dq_rates = numpy.array( + [ifo_grp['bins'][b]['dq_rates'][n] for b in bin_names]) + dq_rates = numpy.maximum(dq_rates, 1) + logrates = numpy.log(dq_rates) + + ymax = max(ymax, numpy.max(dq_rates)) -yticks = ax.get_yticks() + offset = width * (n - 1) + ax.bar(x + offset, dq_rates, width, label=dqstate_name, color=colors[n]) -ax_ymax = numpy.exp(ymax) -ax.set_ylim(1, ax_ymax) -new_ticks = range(0, int(numpy.ceil(numpy.log10(ax_ymax)))) -ax.set_yticks([10**t for t in new_ticks]) -ax.set_ylabel('(Rate during DQ Flag)/(Mean Rate) [Min 1]') +ymax = ymax**1.05 +ax.set_ylim(ymin, ymax) ax.set_yscale('log') -ax.set_xlabel('Template Bin Number') +ax.set_ylabel('(Trigger rate during DQ state)/(Mean Rate) [Min 1]') + ax.set_xticks(x) +ax.set_xlabel('Template Bin Number (Longer duration -> Lower number)') +ax.legend() +ax.grid() + +ax2.set_ylabel('DQ Log Likelihood Penalty') +ax2.set_ylim(numpy.log(ymin), numpy.log(ymax)) # add meta data and save figure -plot_title = f'{ifo}: {dq_name} flag log likelihood' -plot_caption = f'The log likelihood correction during during times flagged \ - by the {dq_name} flag' -pycbc.results.save_fig_with_metadata(fig, args.output_file, title=plot_title, - caption=plot_caption, cmd=' '.join(sys.argv)) +plot_title = f'{ifo} DQ Trigger Rates' +plot_caption = 'The log likelihood correction \ + during during each dq state for each template bin.' +pycbc.results.save_fig_with_metadata( + fig, args.output_file, title=plot_title, + caption=plot_caption, cmd=' '.join(sys.argv)) diff --git a/bin/workflows/pycbc_make_offline_search_workflow b/bin/workflows/pycbc_make_offline_search_workflow index e0995d33cc1..a0b16f347de 100755 --- a/bin/workflows/pycbc_make_offline_search_workflow +++ b/bin/workflows/pycbc_make_offline_search_workflow @@ -191,18 +191,14 @@ insps = wf.merge_single_detector_hdf_files(workflow, hdfbank, # 'statfiles' is list of files used in calculating statistic # 'dqfiles' is the subset of files containing data quality information statfiles = [] -dqfiles = [] -dqfile_labels = [] -dq_labels = workflow.cp.get_subsections('workflow-data_quality') -for dq_label in dq_labels: - dq_label_files = wf.setup_dq_reranking(workflow, dq_label, insps, hdfbank, - analyzable_segs, analyzable_file, - dq_segment_file, - output_dir='dq', - tags=['full_data']) - statfiles += dq_label_files - dqfiles += dq_label_files - dqfile_labels += len(dq_label_files) * [dq_label] + +dqfiles, dqfile_labels = wf.setup_dq_reranking(workflow, insps, + hdfbank, analyzable_file, + analyzable_name, + dq_segment_file, + output_dir='dq', + tags=['full_data']) +statfiles += dqfiles statfiles += wf.setup_trigger_fitting(workflow, insps, hdfbank, final_veto_file, final_veto_name, @@ -554,32 +550,18 @@ for key in final_bg_files: # DQ log likelihood plots for dqf, dql in zip(dqfiles, dqfile_labels): - dq_type = workflow.cp.get_opt_tags('workflow-data_quality', 'dq-type', tags=[dql]) - tags = [dql, 'full_data'] - if dq_type == 'timeseries': - # plot rates as function of time - outdir = rdir[f'data_quality/{dqf.ifo}_{dql}_trigger_rates'] - wf.make_dq_timeseries_trigger_rate_plot(workflow, [dqf], outdir, - tags=tags) - # plot rates as function of dq percentile - outdir = rdir[f'data_quality/{dqf.ifo}_{dql}_percentiles'] - wf.make_dq_percentile_plot(workflow, [dqf], outdir, tags=tags) - elif dq_type == 'flag': - # plot rates when flag was on - outdir = rdir[f'data_quality/{dqf.ifo}_{dql}_trigger_rates'] - wf.make_dq_flag_trigger_rate_plot(workflow, [dqf], outdir, tags=[dql]) + # plot rates when flag was on + outdir = rdir[f'data_quality/{dqf.ifo}_{dql}_trigger_rates'] + wf.make_dq_flag_trigger_rate_plot(workflow, dqf, dql, outdir, tags=[dql]) # DQ background binning plot -if workflow.cp.has_option_tags('bin_trigger_rates_dq', 'background-bins', - tags=None): - background_bins = workflow.cp.get_opt_tags('bin_trigger_rates_dq', - 'background-bins', tags=None) - bank_tags = workflow.cp.get_subsections('plot_bank') - if len(bank_tags) == 0: - bank_tags = ['bank'] - for b_tag in bank_tags: - wf.make_template_plot(workflow, hdfbank, rdir['data_quality'], - bins=background_bins, tags=['dq_bins', b_tag]) +ifo_bins = workflow.cp.get_subsections('bin_templates') +for ifo in ifo_bins: + background_bins = workflow.cp.get_opt_tags( + 'bin_templates', 'background-bins', tags=[ifo]) + wf.make_template_plot(workflow, hdfbank, rdir['data_quality'], + bins=background_bins, + tags=[ifo, 'dq_bins'] + bank_tags) ############################## Setup the injection runs ####################### diff --git a/examples/search/plotting.ini b/examples/search/plotting.ini index 4ecf6cd01c2..331adb9abbf 100644 --- a/examples/search/plotting.ini +++ b/examples/search/plotting.ini @@ -11,6 +11,9 @@ section-header = loudest_noncoinc_time [workflow-sngl_minifollowups-all] section-header = all +[workflow-sngl_minifollowups-vetoed] +section-header = loudest_within_vetoes + [workflow-injection_minifollowups] num-events=1 subsection-suffix=with_ifar_lt_1_year @@ -24,7 +27,6 @@ analysis-category = background_exc sort-variable = stat [singles_minifollowup] -ranking-statistic = single_ranking_only sngl-ranking = newsnr_sgveto_psdvar [singles_minifollowup-noncoinc] @@ -32,13 +34,23 @@ non-coinc-time-only = [singles_minifollowup-all] +[singles_minifollowup-vetoed] +vetoed-time-only = +ranking-statistic = ${sngls|ranking-statistic} +statistic-files = ${sngls|statistic-files} + [injection_minifollowup] ifar-threshold = 1 +[injection_minifollowup&plot_foundmissed-sub_mchirp_grad&plot_foundmissed-all_mchirp_grad&plot_foundmissed-summary] +distance-type = comb_optimal_snr + [page_snglinfo] -ranking-statistic = single_ranking_only sngl-ranking = newsnr_sgveto_psdvar +[page_snglinfo-vetoed] +ranking-statistic = ${sngls|ranking-statistic} + [single_template_plot] [single_template_plot-p1] @@ -139,14 +151,12 @@ log-dist = far-type = exclusive [plot_foundmissed-sub_mchirp_grad&plot_foundmissed-all_mchirp_grad&plot_foundmissed-summary] -distance-type = comb_optimal_snr axis-type = mchirp log-x = log-distance = gradient-far = [plot_foundmissed-sub_mchirp_gradm&plot_foundmissed-all_mchirp_gradm&plot_foundmissed-summarym] -distance-type = comb_optimal_snr axis-type = mchirp log-x = log-distance = diff --git a/pycbc/__init__.py b/pycbc/__init__.py index 53c35debfff..6f2d5a6c3c4 100644 --- a/pycbc/__init__.py +++ b/pycbc/__init__.py @@ -66,7 +66,28 @@ def formatTime(self, record, datefmt=None): return s -def init_logging(verbose=False, format='%(asctime)s %(message)s'): +def add_common_pycbc_options(parser): + """ + Common utility to add standard options to each PyCBC executable. + + Parameters + ---------- + parser : argparse.ArgumentParser + The argument parser to which the options will be added + """ + group = parser.add_argument_group( + title="PyCBC common options", + description="Common options for PyCBC executables.", + ) + group.add_argument('-V', '--verbose', action='count', default=0, + help='Add verbosity to logging. Adding the option ' + 'multiple times makes logging progressively ' + 'more verbose, e.g. --verbose or -V provides ' + 'logging at the info level, but -VV or ' + '--verbose --verbose provides debug logging.') + +def init_logging(verbose=False, + format='%(asctime)s %(levelname)s : %(message)s'): """Common utility for setting up logging in PyCBC. Installs a signal handler such that verbosity can be activated at diff --git a/pycbc/events/coinc.py b/pycbc/events/coinc.py index 7b8ee9840d2..8fd9420b076 100644 --- a/pycbc/events/coinc.py +++ b/pycbc/events/coinc.py @@ -35,13 +35,7 @@ from .eventmgr_cython import timecoincidence_findidxlen from .eventmgr_cython import timecluster_cython -# Mapping used in background_bin_from_string to select approximant for -# duration function, if duration-based binning is used. -_APPROXIMANT_DURATION_MAP = { - 'SEOBNRv2duration': 'SEOBNRv2', - 'SEOBNRv4duration': 'SEOBNRv4', - 'SEOBNRv5duration': 'SEOBNRv5_ROM' -} +logger = logging.getLogger('pycbc.events.coinc') def background_bin_from_string(background_bins, data): @@ -104,7 +98,7 @@ def background_bin_from_string(background_bins, data): elif bin_type == 'chi_eff': vals = pycbc.conversions.chi_eff(data['mass1'], data['mass2'], data['spin1z'], data['spin2z']) - elif bin_type in ['SEOBNRv2Peak', 'SEOBNRv4Peak']: + elif bin_type.endswith('Peak'): vals = pycbc.pnutils.get_freq( 'f' + bin_type, data['mass1'], @@ -113,14 +107,14 @@ def background_bin_from_string(background_bins, data): data['spin2z'] ) cached_values[bin_type] = vals - elif bin_type in _APPROXIMANT_DURATION_MAP: + elif bin_type.endswith('duration'): vals = pycbc.pnutils.get_imr_duration( data['mass1'], data['mass2'], data['spin1z'], data['spin2z'], data['f_lower'], - approximant=_APPROXIMANT_DURATION_MAP[bin_type] + approximant=bin_type.replace('duration', '') ) cached_values[bin_type] = vals else: @@ -308,7 +302,7 @@ def win(ifo1, ifo2): # tested against fixed and pivot are now present for testing with new # dependent ifos for ifo2 in ids: - logging.info('added ifo %s, testing against %s' % (ifo1, ifo2)) + logger.info('added ifo %s, testing against %s' % (ifo1, ifo2)) w = win(ifo1, ifo2) left = time1.searchsorted(ctimes[ifo2] - w) right = time1.searchsorted(ctimes[ifo2] + w) @@ -325,9 +319,9 @@ def win(ifo1, ifo2): # However there are rare corner cases at starts/ends of inspiral # jobs. For these, arbitrarily keep the first trigger and # discard the second (and any subsequent ones). - logging.warning('Triggers in %s are closer than coincidence ' - 'window, 1 or more coincs will be discarded. ' - 'This is a warning, not an error.' % ifo1) + logger.warning('Triggers in %s are closer than coincidence ' + 'window, 1 or more coincs will be discarded. ' + 'This is a warning, not an error.' % ifo1) # identify indices of times in ifo1 that form coincs with ifo2 dep_ids = left[nz] # slide is array of slide ids attached to pivot ifo @@ -372,7 +366,7 @@ def cluster_coincs(stat, time1, time2, timeslide_id, slide, window, **kwargs): The set of indices corresponding to the surviving coincidences. """ if len(time1) == 0 or len(time2) == 0: - logging.info('No coinc triggers in one, or both, ifos.') + logger.info('No coinc triggers in one, or both, ifos.') return numpy.array([]) if numpy.isfinite(slide): @@ -387,9 +381,9 @@ def cluster_coincs(stat, time1, time2, timeslide_id, slide, window, **kwargs): span = (time.max() - time.min()) + window * 10 time = time + span * tslide - logging.info('Clustering events over %s s window', window) + logger.info('Clustering events over %s s window', window) cidx = cluster_over_time(stat, time, window, **kwargs) - logging.info('%d triggers remaining', len(cidx)) + logger.info('%d triggers remaining', len(cidx)) return cidx @@ -418,7 +412,7 @@ def cluster_coincs_multiifo(stat, time_coincs, timeslide_id, slide, window, """ time_coinc_zip = list(zip(*time_coincs)) if len(time_coinc_zip) == 0: - logging.info('No coincident triggers.') + logger.info('No coincident triggers.') return numpy.array([]) time_avg_num = [] @@ -442,9 +436,9 @@ def cluster_coincs_multiifo(stat, time_coincs, timeslide_id, slide, window, span = (time_avg.max() - time_avg.min()) + window * 10 time_avg = time_avg + span * tslide - logging.info('Clustering events over %s s window', window) + logger.info('Clustering events over %s s window', window) cidx = cluster_over_time(stat, time_avg, window, **kwargs) - logging.info('%d triggers remaining', len(cidx)) + logger.info('%d triggers remaining', len(cidx)) return cidx @@ -506,7 +500,7 @@ def cluster_over_time(stat, time, window, method='python', right = time.searchsorted(time + window) indices = numpy.zeros(len(left), dtype=numpy.uint32) - logging.debug('%d triggers before clustering', len(time)) + logger.debug('%d triggers before clustering', len(time)) if method == 'cython': j = timecluster_cython(indices, left, right, stat, len(left)) @@ -546,7 +540,7 @@ def cluster_over_time(stat, time, window, method='python', indices = indices[:j] - logging.debug('%d triggers remaining', len(indices)) + logger.debug('%d triggers remaining', len(indices)) return time_sorting[indices] @@ -941,9 +935,9 @@ def pick_best_coinc(cls, coinc_results): # apply trials factor for the best coinc if mresult: mresult['foreground/ifar'] = mifar / float(trials) - logging.info('Found %s coinc with ifar %s', - mresult['foreground/type'], - mresult['foreground/ifar']) + logger.info('Found %s coinc with ifar %s', + mresult['foreground/type'], + mresult['foreground/ifar']) return mresult # If no coinc, just return one of the results dictionaries. They will # all contain the same results (i.e. single triggers) in this case. @@ -1090,7 +1084,7 @@ def _add_singles_to_buffer(self, results, ifos): # convert to single detector trigger values # FIXME Currently configured to use pycbc live output # where chisq is the reduced chisq and chisq_dof is the actual DOF - logging.info("adding singles to the background estimate...") + logger.info("adding singles to the background estimate...") updated_indices = {} for ifo in ifos: trigs = results[ifo] @@ -1246,7 +1240,7 @@ def _find_coincs(self, results, valid_ifos): for ifo in valid_ifos: trigger_ids[ifo] = numpy.concatenate(trigger_ids[ifo]).astype(numpy.int32) - logging.info( + logger.info( "%s: %s background and zerolag coincs", ppdets(self.ifos, "-"), len(cstat) ) @@ -1259,7 +1253,7 @@ def _find_coincs(self, results, valid_ifos): offsets = numpy.concatenate(offsets) ctime0 = numpy.concatenate(ctimes[self.ifos[0]]).astype(numpy.float64) ctime1 = numpy.concatenate(ctimes[self.ifos[1]]).astype(numpy.float64) - logging.info("Clustering %s coincs", ppdets(self.ifos, "-")) + logger.info("Clustering %s coincs", ppdets(self.ifos, "-")) cidx = cluster_coincs(cstat, ctime0, ctime1, offsets, self.timeslide_interval, self.analysis_block + 2*self.time_window, @@ -1342,7 +1336,7 @@ def add_singles(self, results): A dictionary of arrays containing the coincident results. """ # Let's see how large everything is - logging.info( + logger.info( "%s: %s coincs, %s bytes", ppdets(self.ifos, "-"), len(self.coincs), self.coincs.nbytes ) diff --git a/pycbc/events/coinc_rate.py b/pycbc/events/coinc_rate.py index 2b402b2e9b9..30ffa14e546 100644 --- a/pycbc/events/coinc_rate.py +++ b/pycbc/events/coinc_rate.py @@ -14,6 +14,8 @@ import numpy import pycbc.detector +logger = logging.getLogger('pycbc.events.coinc_rate') + def multiifo_noise_lograte(log_rates, slop): """ @@ -82,9 +84,9 @@ def combination_noise_rate(rates, slop): numpy array Expected coincidence rate in the combination, units Hz """ - logging.warning('combination_noise_rate() is liable to numerical ' - 'underflows, use combination_noise_lograte ' - 'instead') + logger.warning('combination_noise_rate() is liable to numerical ' + 'underflows, use combination_noise_lograte ' + 'instead') log_rates = {k: numpy.log(r) for (k, r) in rates.items()} # exp may underflow return numpy.exp(combination_noise_lograte(log_rates, slop)) diff --git a/pycbc/events/cuts.py b/pycbc/events/cuts.py index 750b06289d3..601fb283820 100644 --- a/pycbc/events/cuts.py +++ b/pycbc/events/cuts.py @@ -37,6 +37,8 @@ # Only used to check isinstance: from pycbc.io.hdf import ReadByTemplate +logger = logging.getLogger('pycbc.events.cuts') + # sngl_rank_keys are the allowed names of reweighted SNR functions sngl_rank_keys = ranking.sngls_ranking_function_dict.keys() @@ -97,8 +99,8 @@ def convert_inputstr(inputstr, choices): try: cut_param, cut_value_str, cut_limit = inputstr.split(':') except ValueError as value_e: - logging.warning("ERROR: Cut string format not correct, please " - "supply as PARAMETER:VALUE:LIMIT") + logger.warning("ERROR: Cut string format not correct, please " + "supply as PARAMETER:VALUE:LIMIT") raise value_e if cut_param.lower() not in choices: @@ -113,8 +115,8 @@ def convert_inputstr(inputstr, choices): try: cut_value = float(cut_value_str) except ValueError as value_e: - logging.warning("ERROR: Cut value must be convertible into a float, " - "got '%s'.", cut_value_str) + logger.warning("ERROR: Cut value must be convertible into a float, " + "got '%s'.", cut_value_str) raise value_e return {(cut_param, ineq_functions[cut_limit]): cut_value} @@ -137,9 +139,9 @@ def check_update_cuts(cut_dict, new_cut): new_cut_key = list(new_cut.keys())[0] if new_cut_key in cut_dict: # The cut has already been called - logging.warning("WARNING: Cut parameter %s and function %s have " - "already been used. Utilising the strictest cut.", - new_cut_key[0], new_cut_key[1].__name__) + logger.warning("WARNING: Cut parameter %s and function %s have " + "already been used. Utilising the strictest cut.", + new_cut_key[0], new_cut_key[1].__name__) # Extract the function and work out which is strictest cut_function = new_cut_key[1] value_new = list(new_cut.values())[0] @@ -148,17 +150,17 @@ def check_update_cuts(cut_dict, new_cut): # The new threshold would survive the cut of the # old threshold, therefore the new threshold is stricter # - update it - logging.warning("WARNING: New threshold of %.3f is " - "stricter than old threshold %.3f, " - "using cut at %.3f.", - value_new, value_old, value_new) + logger.warning("WARNING: New threshold of %.3f is " + "stricter than old threshold %.3f, " + "using cut at %.3f.", + value_new, value_old, value_new) cut_dict.update(new_cut) else: # New cut would not make a difference, ignore it - logging.warning("WARNING: New threshold of %.3f is less " - "strict than old threshold %.3f, using " - "cut at %.3f.", - value_new, value_old, value_old) + logger.warning("WARNING: New threshold of %.3f is less " + "strict than old threshold %.3f, using " + "cut at %.3f.", + value_new, value_old, value_old) else: # This is a new cut - add it cut_dict.update(new_cut) diff --git a/pycbc/events/eventmgr.py b/pycbc/events/eventmgr.py index 56ad265c1d1..ec9bd672a10 100644 --- a/pycbc/events/eventmgr.py +++ b/pycbc/events/eventmgr.py @@ -37,6 +37,8 @@ from .eventmgr_cython import findchirp_cluster_over_window_cython +logger = logging.getLogger('pycbc.events.eventmgr') + @schemed("pycbc.events.threshold_") def threshold(series, value): """Return list of values and indices values over threshold in series. @@ -196,7 +198,7 @@ def save_state(self, tnum_finished, filename): from pycbc.io.hdf import dump_state self.tnum_finished = tnum_finished - logging.info('Writing checkpoint file at template %s', tnum_finished) + logger.info('Writing checkpoint file at template %s', tnum_finished) fp = h5py.File(filename, 'w') dump_state(self, fp, protocol=pickle.HIGHEST_PROTOCOL) fp.close() @@ -214,7 +216,7 @@ def restore_state(filename): raise e fp.close() next_template = mgr.tnum_finished + 1 - logging.info('Restoring with checkpoint at template %s', next_template) + logger.info('Restoring with checkpoint at template %s', next_template) return mgr.tnum_finished + 1, mgr @classmethod @@ -351,35 +353,35 @@ def finalize_template_events(self): def consolidate_events(self, opt, gwstrain=None): self.events = numpy.concatenate(self.accumulate) - logging.info("We currently have %d triggers", len(self.events)) + logger.info("We currently have %d triggers", len(self.events)) if opt.chisq_threshold and opt.chisq_bins: - logging.info("Removing triggers with poor chisq") + logger.info("Removing triggers with poor chisq") self.chisq_threshold(opt.chisq_threshold, opt.chisq_bins, opt.chisq_delta) - logging.info("%d remaining triggers", len(self.events)) + logger.info("%d remaining triggers", len(self.events)) if opt.newsnr_threshold and opt.chisq_bins: - logging.info("Removing triggers with NewSNR below threshold") + logger.info("Removing triggers with NewSNR below threshold") self.newsnr_threshold(opt.newsnr_threshold) - logging.info("%d remaining triggers", len(self.events)) + logger.info("%d remaining triggers", len(self.events)) if opt.keep_loudest_interval: - logging.info("Removing triggers not within the top %s " - "loudest of a %s second interval by %s", - opt.keep_loudest_num, opt.keep_loudest_interval, - opt.keep_loudest_stat) + logger.info("Removing triggers not within the top %s " + "loudest of a %s second interval by %s", + opt.keep_loudest_num, opt.keep_loudest_interval, + opt.keep_loudest_stat) self.keep_loudest_in_interval\ (opt.keep_loudest_interval * opt.sample_rate, opt.keep_loudest_num, statname=opt.keep_loudest_stat, log_chirp_width=opt.keep_loudest_log_chirp_window) - logging.info("%d remaining triggers", len(self.events)) + logger.info("%d remaining triggers", len(self.events)) if opt.injection_window and hasattr(gwstrain, 'injections'): - logging.info("Keeping triggers within %s seconds of injection", - opt.injection_window) + logger.info("Keeping triggers within %s seconds of injection", + opt.injection_window) self.keep_near_injection(opt.injection_window, gwstrain.injections) - logging.info("%d remaining triggers", len(self.events)) + logger.info("%d remaining triggers", len(self.events)) self.accumulate = [self.events] diff --git a/pycbc/events/significance.py b/pycbc/events/significance.py index 9a58d2c3826..86f2a99e7ea 100644 --- a/pycbc/events/significance.py +++ b/pycbc/events/significance.py @@ -32,6 +32,8 @@ import numpy as np from pycbc.events import trigger_fits as trstats +logger = logging.getLogger('pycbc.events.significance') + def count_n_louder(bstat, fstat, dec, **kwargs): # pylint:disable=unused-argument @@ -280,8 +282,8 @@ def positive_float(inp): """ fl_in = float(inp) if fl_in < 0: - logging.warning("Value provided to positive_float is less than zero, " - "this is not allowed") + logger.warning("Value provided to positive_float is less than zero, " + "this is not allowed") raise ValueError return fl_in @@ -406,8 +408,8 @@ def digest_significance_options(combo_keys, args): # Allow options for detector combos that are not actually # used/required for a given job. Such options have # no effect, but emit a warning for (e.g.) diagnostic checks - logging.warning("Key %s not used by this code, uses %s", - combo, combo_keys) + logger.warning("Key %s not used by this code, uses %s", + combo, combo_keys) significance_dict[combo] = copy.deepcopy(_default_opt_dict) significance_dict[combo][argument_key] = conv_func(value) @@ -442,8 +444,8 @@ def apply_far_limit(far, significance_dict, combo=None): if significance_dict[combo]['far_limit'] == 0: return far_out far_limit_str = f"{significance_dict[combo]['far_limit']:.3e}" - logging.info("Applying FAR limit of %s to %s events", - far_limit_str, combo) + logger.info("Applying FAR limit of %s to %s events", + far_limit_str, combo) far_out = np.maximum(far, significance_dict[combo]['far_limit']) else: # IFO combo supplied as an array, by e.g. pycbc_add_statmap @@ -453,8 +455,8 @@ def apply_far_limit(far, significance_dict, combo=None): if significance_dict[ifo_combo]['far_limit'] == 0: continue far_limit_str = f"{significance_dict[ifo_combo]['far_limit']:.3e}" - logging.info("Applying FAR limit of %s to %s events", - far_limit_str, ifo_combo) + logger.info("Applying FAR limit of %s to %s events", + far_limit_str, ifo_combo) this_combo_idx = combo == ifo_combo.encode('utf-8') far_out[this_combo_idx] = np.maximum( far[this_combo_idx], diff --git a/pycbc/events/single.py b/pycbc/events/single.py index f41b2dfb131..626efe81807 100644 --- a/pycbc/events/single.py +++ b/pycbc/events/single.py @@ -8,6 +8,8 @@ from pycbc import conversions as conv from pycbc import bin_utils +logger = logging.getLogger('pycbc.events.single') + class LiveSingle(object): def __init__(self, ifo, @@ -198,7 +200,7 @@ def check(self, trigs, data_reader): return candidate def calculate_ifar(self, sngl_ranking, duration): - logging.info("Calculating IFAR") + logger.info("Calculating IFAR") if self.fixed_ifar and self.ifo in self.fixed_ifar: return self.fixed_ifar[self.ifo] @@ -211,7 +213,7 @@ def calculate_ifar(self, sngl_ranking, duration): rates = dist_grp['counts'][:] / live_time coeffs = dist_grp['fit_coeff'][:] except FileNotFoundError: - logging.error( + logger.error( 'Single fit file %s not found; ' 'dropping a potential single-detector candidate!', self.fit_file diff --git a/pycbc/events/stat.py b/pycbc/events/stat.py index 4151571e356..59295cc0263 100644 --- a/pycbc/events/stat.py +++ b/pycbc/events/stat.py @@ -33,6 +33,8 @@ from .eventmgr_cython import logsignalrateinternals_computepsignalbins from .eventmgr_cython import logsignalrateinternals_compute2detrate +logger = logging.getLogger('pycbc.events.stat') + class Stat(object): """Base class which should be extended to provide a coincident statistic""" @@ -64,7 +66,7 @@ def __init__(self, sngl_ranking, files=None, ifos=None, **kwargs): if stat in self.files: raise RuntimeError("We already have one file with stat attr =" " %s. Can't provide more than one!" % stat) - logging.info("Found file %s for stat %s", filename, stat) + logger.info("Found file %s for stat %s", filename, stat) self.files[stat] = filename # Provide the dtype of the single detector method's output @@ -381,7 +383,7 @@ def get_hist(self, ifos=None): if selected is None and len(ifos) > 1: raise RuntimeError("Couldn't figure out which stat file to use") - logging.info("Using signal histogram %s for ifos %s", selected, ifos) + logger.info("Using signal histogram %s for ifos %s", selected, ifos) weights = {} param = {} @@ -390,7 +392,7 @@ def get_hist(self, ifos=None): # Patch for pre-hdf5=3.0 histogram files try: - logging.info("Decoding hist ifos ..") + logger.info("Decoding hist ifos ..") self.hist_ifos = [i.decode('UTF-8') for i in self.hist_ifos] except (UnicodeDecodeError, AttributeError): pass @@ -1996,24 +1998,18 @@ def __init__(self, sngl_ranking, files=None, ifos=None, """ ExpFitFgBgNormStatistic.__init__(self, sngl_ranking, files=files, ifos=ifos, **kwargs) - self.dq_val_by_time = {} - self.dq_bin_by_id = {} - for k in self.files.keys(): - parsed_attrs = k.split('-') - if len(parsed_attrs) < 3: - continue - if parsed_attrs[2] == 'dq_ts_reference': - ifo = parsed_attrs[0] - dq_type = parsed_attrs[1] - dq_vals = self.assign_dq_val(k) - dq_bins = self.assign_bin_id(k) - if ifo not in self.dq_val_by_time: - self.dq_val_by_time[ifo] = {} - self.dq_bin_by_id[ifo] = {} - self.dq_val_by_time[ifo][dq_type] = dq_vals - self.dq_bin_by_id[ifo][dq_type] = dq_bins - - def assign_bin_id(self, key): + self.dq_rates_by_state = {} + self.dq_bin_by_tid = {} + self.dq_state_segments = {} + + for ifo in self.ifos: + key = f'{ifo}-dq_stat_info' + if key in self.files.keys(): + self.dq_rates_by_state[ifo] = self.assign_dq_rates(key) + self.dq_bin_by_tid[ifo] = self.assign_template_bins(key) + self.dq_state_segments[ifo] = self.setup_segments(key) + + def assign_template_bins(self, key): """ Assign bin ID values Assign each template id to a bin name based on a @@ -2031,18 +2027,18 @@ def assign_bin_id(self, key): """ ifo = key.split('-')[0] with h5py.File(self.files[key], 'r') as dq_file: - bin_names = dq_file.attrs['names'][:] - locs = [] - names = [] - for bin_name in bin_names: - bin_locs = dq_file[ifo + '/locs/' + bin_name][:] - locs = list(locs) + list(bin_locs.astype(int)) - names = list(names) + list([bin_name] * len(bin_locs)) - - bin_dict = dict(zip(locs, names)) + tids = [] + bin_nums = [] + bin_grp = dq_file[f'{ifo}/bins'] + for bin_name in bin_grp.keys(): + bin_tids = bin_grp[f'{bin_name}/tids'][:] + tids = list(tids) + list(bin_tids.astype(int)) + bin_nums = list(bin_nums) + list([bin_name] * len(bin_tids)) + + bin_dict = dict(zip(tids, bin_nums)) return bin_dict - def assign_dq_val(self, key): + def assign_dq_rates(self, key): """ Assign dq values to each time for every bin based on a referenced statistic file. @@ -2061,38 +2057,73 @@ def assign_dq_val(self, key): """ ifo = key.split('-')[0] with h5py.File(self.files[key], 'r') as dq_file: - times = dq_file[ifo + '/times'][:] - bin_names = dq_file.attrs['names'][:] + bin_grp = dq_file[f'{ifo}/bins'] dq_dict = {} - for bin_name in bin_names: - dq_vals = dq_file[ifo + '/dq_vals/' + bin_name][:] - dq_dict[bin_name] = dict(zip(times, dq_vals)) + for bin_name in bin_grp.keys(): + dq_dict[bin_name] = bin_grp[f'{bin_name}/dq_rates'][:] return dq_dict - def find_dq_val(self, trigs): - """Get dq values for a specific ifo and times""" - time = trigs['end_time'].astype(int) + def setup_segments(self, key): + """ + Check if segments definitions are in stat file + If they are, we are running offline and need to store them + If they aren't, we are running online + """ + ifo = key.split('-')[0] + with h5py.File(self.files[key], 'r') as dq_file: + ifo_grp = dq_file[ifo] + dq_state_segs_dict = {} + for k in ifo_grp['dq_segments'].keys(): + seg_dict = {} + seg_dict['start'] = \ + ifo_grp[f'dq_segments/{k}/segment_starts'][:] + seg_dict['end'] = \ + ifo_grp[f'dq_segments/{k}/segment_ends'][:] + dq_state_segs_dict[k] = seg_dict + + return dq_state_segs_dict + + def find_dq_noise_rate(self, trigs, dq_state): + """Get dq values for a specific ifo and dq states""" + try: tnum = trigs.template_num - ifo = trigs.ifo except AttributeError: tnum = trigs['template_id'] - assert len(self.ifos) == 1 + + try: + ifo = trigs.ifo + except AttributeError: + ifo = trigs['ifo'] + assert len(numpy.unique(ifo)) == 1 # Should be exactly one ifo provided - ifo = self.ifos[0] - dq_val = numpy.zeros(len(time)) - if ifo in self.dq_val_by_time: - for (i, t) in enumerate(time): - for k in self.dq_val_by_time[ifo].keys(): - if isinstance(tnum, numpy.ndarray): - bin_name = self.dq_bin_by_id[ifo][k][tnum[i]] - else: - bin_name = self.dq_bin_by_id[ifo][k][tnum] - val = self.dq_val_by_time[ifo][k][bin_name][int(t)] - dq_val[i] = max(dq_val[i], val) + ifo = ifo[0] + + dq_val = numpy.zeros(len(dq_state)) + + if ifo in self.dq_rates_by_state: + for (i, st) in enumerate(dq_state): + if isinstance(tnum, numpy.ndarray): + bin_name = self.dq_bin_by_tid[ifo][tnum[i]] + else: + bin_name = self.dq_bin_by_tid[ifo][tnum] + dq_val[i] = self.dq_rates_by_state[ifo][bin_name][st] return dq_val + def find_dq_state_by_time(self, ifo, times): + """Get the dq state for an ifo at times""" + dq_state = numpy.zeros(len(times), dtype=numpy.uint8) + if ifo in self.dq_state_segments: + from pycbc.events.veto import indices_within_times + for k in self.dq_state_segments[ifo]: + starts = self.dq_state_segments[ifo][k]['start'] + ends = self.dq_state_segments[ifo][k]['end'] + inds = indices_within_times(times, starts, ends) + # states are named in file as 'dq_state_N', need to extract N + dq_state[inds] = int(k[9:]) + return dq_state + def lognoiserate(self, trigs): """ Calculate the log noise rate density over single-ifo ranking @@ -2107,12 +2138,27 @@ def lognoiserate(self, trigs): Returns --------- - lognoisel: numpy.array + lognoiserate: numpy.array Array of log noise rate density for each input trigger. """ + + # make sure every trig has a dq state + + try: + ifo = trigs.ifo + except AttributeError: + ifo = trigs['ifo'] + assert len(numpy.unique(ifo)) == 1 + # Should be exactly one ifo provided + ifo = ifo[0] + + dq_state = self.find_dq_state_by_time(ifo, trigs['end_time'][:]) + dq_rate = self.find_dq_noise_rate(trigs, dq_state) + dq_rate = numpy.maximum(dq_rate, 1) + logr_n = ExpFitFgBgNormStatistic.lognoiserate( self, trigs) - logr_n += self.find_dq_val(trigs) + logr_n += numpy.log(dq_rate) return logr_n diff --git a/pycbc/events/trigger_fits.py b/pycbc/events/trigger_fits.py index 961fdafea48..e2a6002e05b 100644 --- a/pycbc/events/trigger_fits.py +++ b/pycbc/events/trigger_fits.py @@ -52,6 +52,8 @@ import numpy from scipy.stats import kstest +logger = logging.getLogger('pycbc.events.trigger_fits') + def exponential_fitalpha(vals, thresh, w): """ Maximum likelihood estimator for the fit factor for @@ -127,8 +129,8 @@ def fit_above_thresh(distr, vals, thresh=None, weights=None): above_thresh = vals >= thresh if numpy.count_nonzero(above_thresh) == 0: # Nothing is above threshold - warn and return -1 - logging.warning("No values are above the threshold, %.2f, " - "maximum is %.2f.", thresh, vals.max()) + logger.warning("No values are above the threshold, %.2f, " + "maximum is %.2f.", thresh, vals.max()) return -1., -1. vals = vals[above_thresh] diff --git a/pycbc/frame/frame.py b/pycbc/frame/frame.py index 08658e380cd..c72309d31b0 100644 --- a/pycbc/frame/frame.py +++ b/pycbc/frame/frame.py @@ -17,14 +17,19 @@ This modules contains functions for reading in data from frame files or caches """ -import lalframe, logging -import lal -import numpy +import logging +import warnings +import os.path +import glob +import time import math -import os.path, glob, time +import re +from urllib.parse import urlparse +import numpy +import lalframe +import lal from gwdatafind import find_urls as find_frame_urls import pycbc -from urllib.parse import urlparse from pycbc.types import TimeSeries, zeros @@ -256,13 +261,15 @@ def read_frame(location, channels, start_time=None, else: return _read_channel(channels, stream, start_time, duration) -def frame_paths(frame_type, start_time, end_time, server=None, url_type='file'): - """Return the paths to a span of frame files +def frame_paths( + frame_type, start_time, end_time, server=None, url_type='file', site=None +): + """Return the paths to a span of frame files. Parameters ---------- frame_type : string - The string representation of the frame type (ex. 'H1_ER_C00_L1') + The string representation of the frame type (ex. 'H1_ER_C00_L1'). start_time : int The start time that we need the frames to span. end_time : int @@ -274,6 +281,11 @@ def frame_paths(frame_type, start_time, end_time, server=None, url_type='file'): Returns only frame URLs with a particular scheme or head such as "file" or "https". Default is "file", which queries locally stored frames. Option can be disabled if set to None. + site : string, optional + One-letter string specifying which site you want data from (H, L, V, + etc). If not given, the site is assumed to be the first letter of + `frame_type`, which is usually (but not always) a safe assumption. + Returns ------- paths : list of paths @@ -283,22 +295,73 @@ def frame_paths(frame_type, start_time, end_time, server=None, url_type='file'): -------- >>> paths = frame_paths('H1_LDAS_C02_L2', 968995968, 968995968+2048) """ - site = frame_type[0] + if site is None: + # this case is tolerated for backward compatibility + site = frame_type[0] + warnings.warn( + f'Guessing site {site} from frame type {frame_type}', + DeprecationWarning + ) cache = find_frame_urls(site, frame_type, start_time, end_time, urltype=url_type, host=server) return [urlparse(entry).path for entry in cache] + +def get_site_from_type_or_channel(frame_type, channels): + """Determine the site for querying gwdatafind (H, L, V, etc) based on + substrings of the frame type and channel(s). + + The type should begin with S: or SN:, in which case S is taken as the + site. Otherwise, the same is done with the channel (with the first + channel if more than one are given). If that also fails, the site is + taken to be the first letter of the frame type, which is usually + (but not always) a correct assumption. + + Parameters + ---------- + frame_type : string + The frame type, ideally prefixed by the site indicator. + channels : string or list of strings + The channel name or names. + + Returns + ------- + site : string + The site letter. + frame_type : string + The frame type with the site prefix (if any) removed. + """ + site_re = '^([^:])[^:]?:' + m = re.match(site_re, frame_type) + if m: + return m.groups(1)[0], frame_type[m.end():] + chan = channels + if isinstance(chan, list): + chan = channels[0] + m = re.match(site_re, chan) + if m: + return m.groups(1)[0], frame_type + warnings.warn( + f'Guessing site {frame_type[0]} from frame type {frame_type}', + DeprecationWarning + ) + return frame_type[0], frame_type + + def query_and_read_frame(frame_type, channels, start_time, end_time, sieve=None, check_integrity=False): """Read time series from frame data. - Query for the locatin of physical frames matching the frame type. Return + Query for the location of physical frames matching the frame type. Return a time series containing the channel between the given start and end times. Parameters ---------- frame_type : string - The type of frame file that we are looking for. + The type of frame file that we are looking for. The string should begin + with S: or SN:, in which case S is taken as the site to query. If this + is not the case, the site will be guessed from the channel name or from + the type in a different way, which may not work. channels : string or list of strings Either a string that contains the channel name or a list of channel name strings. @@ -324,6 +387,8 @@ def query_and_read_frame(frame_type, channels, start_time, end_time, >>> ts = query_and_read_frame('H1_LDAS_C02_L2', 'H1:LDAS-STRAIN', >>> 968995968, 968995968+2048) """ + site, frame_type = get_site_from_type_or_channel(frame_type, channels) + # Allows compatibility with our standard tools # We may want to place this into a higher level frame getting tool if frame_type in ['LOSC_STRAIN', 'GWOSC_STRAIN']: @@ -337,17 +402,23 @@ def query_and_read_frame(frame_type, channels, start_time, end_time, from pycbc.frame.gwosc import read_frame_gwosc return read_frame_gwosc(channels, start_time, end_time) - logging.info('querying datafind server') - paths = frame_paths(frame_type, int(start_time), int(numpy.ceil(end_time))) - logging.info('found files: %s' % (' '.join(paths))) - return read_frame(paths, channels, - start_time=start_time, - end_time=end_time, - sieve=sieve, - check_integrity=check_integrity) + logging.info('Querying datafind server') + paths = frame_paths( + frame_type, + int(start_time), + int(numpy.ceil(end_time)), + site=site + ) + logging.info('Found frame file paths: %s', ' '.join(paths)) + return read_frame( + paths, + channels, + start_time=start_time, + end_time=end_time, + sieve=sieve, + check_integrity=check_integrity + ) -__all__ = ['read_frame', 'frame_paths', - 'query_and_read_frame'] def write_frame(location, channels, timeseries): """Write a list of time series to a single frame file. @@ -887,3 +958,15 @@ def null_advance(self, blocksize): """ self.idq.null_advance(blocksize) self.idq_state.null_advance(blocksize) + + +__all__ = [ + 'locations_to_cache', + 'read_frame', + 'query_and_read_frame', + 'frame_paths', + 'write_frame', + 'DataBuffer', + 'StatusBuffer', + 'iDQBuffer' +] diff --git a/pycbc/io/hdf.py b/pycbc/io/hdf.py index 2cd0dd1e450..754c63bf663 100644 --- a/pycbc/io/hdf.py +++ b/pycbc/io/hdf.py @@ -963,6 +963,10 @@ def to_coinc_xml_object(self, file_name): val = sngl_col_vals[name][ifo][0][idx] if name == 'end_time': sngl.end = LIGOTimeGPS(val) + elif name == 'chisq': + # Use reduced chisquared to be consistent with Live + dof = 2. * sngl_col_vals['chisq_dof'][ifo][0][idx] - 2. + sngl.chisq = val / dof else: setattr(sngl, name, val) for name in bank_col_names: diff --git a/pycbc/io/live.py b/pycbc/io/live.py index 754fcbcf563..07474396c1e 100644 --- a/pycbc/io/live.py +++ b/pycbc/io/live.py @@ -17,8 +17,7 @@ make_psd_xmldoc, snr_series_to_xml ) -from pycbc.results import generate_asd_plot -from pycbc.results import ifo_color +from pycbc.results import generate_asd_plot, generate_snr_plot from pycbc.results import source_color from pycbc.mchirp_area import calc_probabilities @@ -446,27 +445,21 @@ def upload(self, fname, gracedb_server=None, testing=True, snr_series_plot_fname = self.basename + '_snr.png' asd_series_plot_fname = self.basename + '_asd.png' - pl.figure() + triggers = { + ifo: (self.coinc_results[f'foreground/{ifo}/end_time'] + + self.time_offset, + self.coinc_results[f'foreground/{ifo}/snr']) + for ifo in self.et_ifos + } ref_time = int(self.merger_time) + generate_snr_plot(self.snr_series, snr_series_plot_fname, + triggers, ref_time) + + generate_asd_plot(self.psds, asd_series_plot_fname) + for ifo in sorted(self.snr_series): curr_snrs = self.snr_series[ifo] curr_snrs.save(snr_series_fname, group='%s/snr' % ifo) - pl.plot(curr_snrs.sample_times - ref_time, abs(curr_snrs), - c=ifo_color(ifo), label=ifo) - if ifo in self.et_ifos: - base = 'foreground/{}/'.format(ifo) - snr = self.coinc_results[base + 'snr'] - mt = (self.coinc_results[base + 'end_time'] - + self.time_offset) - pl.plot([mt - ref_time], [snr], c=ifo_color(ifo), - marker='x') - pl.legend() - pl.xlabel('GPS time from {:d} (s)'.format(ref_time)) - pl.ylabel('SNR') - pl.savefig(snr_series_plot_fname) - pl.close() - - generate_asd_plot(self.psds, asd_series_plot_fname) # Additionally save the PSDs into the snr_series file for ifo in sorted(self.psds): diff --git a/pycbc/pnutils.py b/pycbc/pnutils.py index 95900040a59..9582ce7414e 100644 --- a/pycbc/pnutils.py +++ b/pycbc/pnutils.py @@ -521,7 +521,7 @@ def frequency_cutoff_from_name(name, m1, m2, s1z, s2z): f : float or numpy.array Frequency in Hz """ - params = {"mass1":m1, "mass2":m2, "spin1z":s1z, "spin2z":s2z} + params = {"mass1": m1, "mass2": m2, "spin1z": s1z, "spin2z": s2z} return named_frequency_cutoffs[name](params) def _get_imr_duration(m1, m2, s1z, s2z, f_low, approximant="SEOBNRv4"): @@ -532,20 +532,20 @@ def _get_imr_duration(m1, m2, s1z, s2z, f_low, approximant="SEOBNRv4"): chi = lalsim.SimIMRPhenomBComputeChi(m1, m2, s1z, s2z) time_length = lalsim.SimIMRSEOBNRv2ChirpTimeSingleSpin( m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, chi, f_low) - elif approximant == 'IMRPhenomXAS': + elif approximant == "IMRPhenomXAS": time_length = lalsim.SimIMRPhenomXASDuration( m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z, f_low) elif approximant == "IMRPhenomD": time_length = lalsim.SimIMRPhenomDChirpTime( m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z, f_low) - elif approximant == "SEOBNRv4": - # NB for no clear reason this function has f_low as first argument + elif approximant in ["SEOBNRv4", "SEOBNRv4_ROM"]: + # NB the LALSim function has f_low as first argument time_length = lalsim.SimIMRSEOBNRv4ROMTimeOfFrequency( f_low, m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z) - elif approximant == 'SEOBNRv5_ROM': + elif approximant in ["SEOBNRv5", "SEOBNRv5_ROM"]: time_length = lalsim.SimIMRSEOBNRv5ROMTimeOfFrequency( f_low, m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z) - elif approximant == 'SPAtmplt' or approximant == 'TaylorF2': + elif approximant in ["SPAtmplt", "TaylorF2"]: chi = lalsim.SimInspiralTaylorF2ReducedSpinComputeChi( m1, m2, s1z, s2z ) diff --git a/pycbc/results/__init__.py b/pycbc/results/__init__.py index a12111a2395..63e0aca882b 100644 --- a/pycbc/results/__init__.py +++ b/pycbc/results/__init__.py @@ -4,6 +4,7 @@ from pycbc.results.color import * from pycbc.results.plot import * from pycbc.results.psd import * +from pycbc.results.snr import * from pycbc.results.layout import * from pycbc.results.dq import * from pycbc.results.str_utils import * diff --git a/pycbc/results/snr.py b/pycbc/results/snr.py new file mode 100644 index 00000000000..d96476d6282 --- /dev/null +++ b/pycbc/results/snr.py @@ -0,0 +1,73 @@ +# Copyright (C) 2022 +# Tito Dal Canton, Gareth Cabourn Davies, Ian Harry +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + +# +# ============================================================================= +# +# Preamble +# +# ============================================================================= +# +""" +Module to generate SNR figures +""" +import pylab as pl +from pycbc.results import ifo_color + + +def generate_snr_plot(snrdict, output_filename, triggers, ref_time): + """ + Generate an SNR timeseries plot as used for upload to GraceDB. + + Parameters + ---------- + + snrdict: dictionary + A dictionary keyed on ifo containing the SNR + TimeSeries objects + output_filename: string + The filename for the plot to be saved to + triggers : dictionary of tuples + A dictionary keyed on IFO, containing (trigger time, trigger snr) + ref_time : number, GPS seconds + Reference time which will be used as the zero point of the plot + This should be an integer value, but doesn't need to be an integer + + Returns + ------- + None + """ + pl.figure() + ref_time = int(ref_time) + for ifo in sorted(snrdict): + curr_snrs = snrdict[ifo] + + pl.plot(curr_snrs.sample_times - ref_time, abs(curr_snrs), + c=ifo_color(ifo), label=ifo) + if ifo in triggers: + pl.plot(triggers[ifo][0] - ref_time, + triggers[ifo][1], marker='x', c=ifo_color(ifo)) + + pl.legend() + pl.xlabel(f'GPS time from {ref_time:d} (s)') + pl.ylabel('SNR') + pl.savefig(output_filename) + pl.close() + + +__all__ = ["generate_snr_plot"] diff --git a/pycbc/strain/strain.py b/pycbc/strain/strain.py index 2793012d565..27aff7ccc63 100644 --- a/pycbc/strain/strain.py +++ b/pycbc/strain/strain.py @@ -504,8 +504,10 @@ def insert_strain_option_group(parser, gps_times=True): # Use datafind to get frame files data_reading_group.add_argument("--frame-type", type=str, + metavar="S:TYPE", help="(optional), replaces frame-files. Use datafind " - "to get the needed frame file(s) of this type.") + "to get the needed frame file(s) of this type " + "from site S.") # Filter frame files by URL data_reading_group.add_argument("--frame-sieve", type=str, @@ -704,7 +706,7 @@ def insert_strain_option_group_multi_ifo(parser, gps_times=True): help="Store of time series data in hdf format") # Use datafind to get frame files data_reading_group_multi.add_argument("--frame-type", type=str, nargs="+", - action=MultiDetOptionAction, + action=MultiDetOptionActionSpecial, metavar='IFO:FRAME_TYPE', help="(optional) Replaces frame-files. " "Use datafind to get the needed frame " @@ -1975,9 +1977,12 @@ def from_cli(cls, ifo, args, maxlen): idq_state_channel = ':'.join([ifo, args.idq_state_channel[ifo]]) if args.frame_type: - frame_src = pycbc.frame.frame_paths(args.frame_type[ifo], - args.start_time, - args.end_time) + frame_src = pycbc.frame.frame_paths( + args.frame_type[ifo], + args.start_time, + args.end_time, + site=ifo[0] + ) else: frame_src = [args.frame_src[ifo]] strain_channel = ':'.join([ifo, args.channel_name[ifo]]) diff --git a/pycbc/workflow/dq.py b/pycbc/workflow/dq.py index 7aa727a03e9..1068225dde0 100644 --- a/pycbc/workflow/dq.py +++ b/pycbc/workflow/dq.py @@ -22,183 +22,118 @@ # ============================================================================= # -import os import logging -from ligo import segments -from pycbc.workflow.core import (FileList, Executable, Node, - File, SegFile, make_analysis_dir) -from pycbc.workflow.datafind import setup_datafind_workflow - -class PyCBCCalculateDQExecutable(Executable): - current_retention_level = Executable.ALL_TRIGGERS - def create_node(self, segment, frames): - start = int(segment[0]) - end = int(segment[1]) - node = Node(self) - node.add_input_list_opt('--frame-files', frames) - node.add_opt('--gps-start-time', start) - node.add_opt('--gps-end-time', end) - node.new_output_file_opt(segment, '.hdf', '--output-file') - return node +from pycbc.workflow.core import (FileList, Executable, Node, make_analysis_dir) -class PyCBCRerankDQExecutable(Executable): + +class PyCBCBinTemplatesDQExecutable(Executable): current_retention_level = Executable.MERGED_TRIGGERS - def create_node(self, workflow, ifo, dq_type, dq_files, binned_rate_file): + + def create_node(self, workflow, ifo, template_bank_file): node = Node(self) - node.add_opt('--dq-type', dq_type) node.add_opt('--ifo', ifo) - node.add_input_list_opt('--input-file', dq_files) - node.add_input_opt('--rate-file', binned_rate_file) - node.new_output_file_opt(workflow.analysis_time, '.hdf', - '--output-file') + node.add_input_opt('--bank-file', template_bank_file) + node.new_output_file_opt( + workflow.analysis_time, '.hdf', '--output-file') return node + class PyCBCBinTriggerRatesDQExecutable(Executable): current_retention_level = Executable.MERGED_TRIGGERS - def create_node(self, workflow, ifo, dq_files, trig_file, bank_file): + + def create_node(self, workflow, flag_file, flag_name, + analysis_segment_file, analysis_segment_name, + trig_file, template_bins_file): node = Node(self) - node.add_opt('--ifo', ifo) - node.add_input_opt('--bank-file', bank_file) + node.add_input_opt('--template-bins-file', template_bins_file) node.add_input_opt('--trig-file', trig_file) - node.add_input_list_opt('--dq-file', dq_files) - node.new_output_file_opt(workflow.analysis_time,'.hdf', + node.add_input_opt('--flag-file', flag_file) + node.add_opt('--flag-name', flag_name) + node.add_input_opt('--analysis-segment-file', analysis_segment_file) + node.add_opt('--analysis-segment-name', analysis_segment_name) + node.new_output_file_opt(workflow.analysis_time, '.hdf', '--output-file') return node -class PyCBCCalculateDQFlagExecutable(Executable): - # current_retention_level = Executable.ALL_TRIGGERS - current_retention_level = Executable.MERGED_TRIGGERS - def create_node(self, workflow, segment, dq_file, flag): - node = Node(self) - # Executable objects are initialized with ifo information - start = int(segment[0]) - end = int(segment[1]) - node.add_opt('--ifo', self.ifo_string) - node.add_opt('--flag', flag) - node.add_opt('--gps-start-time', start) - node.add_opt('--gps-end-time', end) - node.add_input_opt('--dq-segments', dq_file) - node.new_output_file_opt(segment, '.hdf', - '--output-file') - return node - -def setup_dq_reranking(workflow, dq_label, insps, bank, - segs, analyzable_file, dq_file, - output_dir=None, tags=None): +def setup_dq_reranking(workflow, insps, bank, + analyzable_seg_file, + analyzable_name, + dq_seg_file, + output_dir=None, tags=None): + logging.info("Setting up dq reranking") make_analysis_dir(output_dir) - output = FileList() - if tags: + output_files = FileList() + output_labels = [] + if tags is None: + tags = [] + + dq_labels = workflow.cp.get_subsections('workflow-data_quality') + + dq_ifos = {} + dq_names = {} + dq_types = {} + for dql in dq_labels: + dq_ifos[dql] = workflow.cp.get_opt_tags( + 'workflow-data_quality', 'dq-ifo', [dql]) + dq_names[dql] = workflow.cp.get_opt_tags( + 'workflow-data_quality', 'dq-name', [dql]) + dq_types[dql] = workflow.cp.get_opt_tags( + 'workflow-data_quality', 'dq-type', [dql]) + + ifos = set(dq_ifos.values()) + + for ifo in ifos: + # get the dq label, type, and name for this ifo + ifo_dq_labels = [dql for dql in dq_labels if (dq_ifos[dql] == ifo)] + assert len(ifo_dq_labels) < 2, f"Received multiple dq files for {ifo}" + dq_label = ifo_dq_labels[0] + dq_name = dq_names[dq_label] + dq_type = dq_types[dq_label] + dq_tags = tags + [dq_label] - else: - dq_tags = [dq_label] - dq_type = workflow.cp.get_opt_tags("workflow-data_quality", - 'dq-type', [dq_label]) - if dq_type == 'timeseries': - if dq_label not in workflow.cp.get_subsections('workflow-datafind'): - msg = """No workflow-datafind section with dq tag. - Tags must be used in workflow-datafind sections " - if more than one source of data is used. - Strain data source must be tagged - workflow-datafind-hoft. - Consult the documentation for more info.""" - raise ValueError(msg) - dq_ifos = workflow.cp.get_opt_tags("workflow-data_quality", - 'ifos', [dq_label]) - dq_ifos = dq_ifos.split(',') - dq_segs = {} - dq_segs_for_file = {} - for ifo in dq_ifos: - dq_segs[ifo] = segs[ifo] - dq_segs_for_file[ifo+':'+dq_label] = segs[ifo] - dq_segs_file = SegFile.from_segment_list_dict(dq_label, - dq_segs_for_file, - extension='.xml', - valid_segment=workflow.analysis_time, - directory=output_dir) - datafind_files, dq_file, dq_segs, dq_name = \ - setup_datafind_workflow(workflow, - dq_segs, "datafind_dq", - seg_file=dq_segs_file, - tags=dq_tags) - for ifo in dq_ifos: - ifo_insp = [insp for insp in insps if (insp.ifo == ifo)] - assert len(ifo_insp)==1 - ifo_insp = ifo_insp[0] - - dq_files = FileList() - for seg in dq_segs[ifo]: - seg_frames = datafind_files.find_all_output_in_range(ifo, seg) - raw_exe = PyCBCCalculateDQExecutable(workflow.cp, - 'calculate_dq', ifos=ifo, - out_dir=output_dir, - tags=dq_tags) - raw_node = raw_exe.create_node(seg, seg_frames) - workflow += raw_node - dq_files += raw_node.output_files - - intermediate_exe = PyCBCBinTriggerRatesDQExecutable(workflow.cp, - 'bin_trigger_rates_dq', - ifos=ifo, - out_dir=output_dir, - tags=dq_tags) - intermediate_node = intermediate_exe.create_node(workflow, ifo, - dq_files, - ifo_insp, bank) - workflow += intermediate_node - binned_rate_file = intermediate_node.output_file - - new_exe = PyCBCRerankDQExecutable(workflow.cp, - 'rerank_dq', ifos=ifo, - out_dir=output_dir, - tags=dq_tags) - new_node = new_exe.create_node(workflow, ifo, dq_label, - dq_files, binned_rate_file) - workflow += new_node - output += new_node.output_files - elif dq_type == 'flag': - flag_str = workflow.cp.get_opt_tags("workflow-data_quality", - 'flag-name', dq_tags) - ifo = flag_str[:2] + + # get triggers for this ifo ifo_insp = [insp for insp in insps if (insp.ifo == ifo)] - assert len(ifo_insp)==1 + assert len(ifo_insp) == 1, \ + f"Received more than one inspiral file for {ifo}" ifo_insp = ifo_insp[0] - flag_name = flag_str - logging.info("Creating job for flag %s", flag_name) - dq_files = FileList() - for seg in segs[ifo]: - raw_exe = PyCBCCalculateDQFlagExecutable(workflow.cp, - 'calculate_dqflag', ifos=ifo, - out_dir=output_dir, - tags=dq_tags) - raw_node = raw_exe.create_node(workflow, seg, dq_file, - flag_name) - workflow += raw_node - dq_files += raw_node.output_files - intermediate_exe = PyCBCBinTriggerRatesDQExecutable(workflow.cp, - 'bin_trigger_rates_dq', - ifos=ifo, - out_dir=output_dir, - tags=dq_tags) - intermediate_node = intermediate_exe.create_node(workflow, ifo, - dq_files, - ifo_insp, bank) - workflow += intermediate_node - binned_rate_file = intermediate_node.output_file - - new_exe = PyCBCRerankDQExecutable(workflow.cp, - 'rerank_dq', ifos=ifo, - out_dir=output_dir, - tags=dq_tags) - new_node = new_exe.create_node(workflow, ifo, dq_label, - dq_files, binned_rate_file) - workflow += new_node - output += new_node.output_files - else: - msg = """Incorrect DQ type specified. - Only valid DQ types are 'flag' - and 'timeseries'. - Consult the documentation for more info.""" - raise ValueError(msg) - - return output + + # calculate template bins for this ifo + bin_templates_exe = PyCBCBinTemplatesDQExecutable( + workflow.cp, + 'bin_templates', + ifos=ifo, + out_dir=output_dir, + tags=tags) + bin_templates_node = bin_templates_exe.create_node(workflow, ifo, bank) + workflow += bin_templates_node + template_bins_file = bin_templates_node.output_file + + if dq_type == 'flag': + flag_file = dq_seg_file + flag_name = dq_name + else: + raise ValueError(f"{dq_type} dq support not yet implemented") + + # calculate trigger rates during dq flags + bin_triggers_exe = PyCBCBinTriggerRatesDQExecutable( + workflow.cp, + 'bin_trigger_rates_dq', + ifos=ifo, + out_dir=output_dir, + tags=dq_tags) + bin_triggers_node = bin_triggers_exe.create_node( + workflow, + flag_file, + flag_name, + analyzable_seg_file, + analyzable_name, + ifo_insp, + template_bins_file) + workflow += bin_triggers_node + output_files += bin_triggers_node.output_files + output_labels += [dq_label] + + logging.info("Finished setting up dq reranking") + return output_files, output_labels diff --git a/pycbc/workflow/plotting.py b/pycbc/workflow/plotting.py index 77d34221bc6..165786c5715 100644 --- a/pycbc/workflow/plotting.py +++ b/pycbc/workflow/plotting.py @@ -29,6 +29,7 @@ from urllib.parse import urljoin from pycbc.workflow.core import File, FileList, makedir, Executable + def excludestr(tags, substr): if substr is None: return tags @@ -38,11 +39,13 @@ def excludestr(tags, substr): substr = substr[0] return [tag for tag in tags if substr not in tag] + def requirestr(tags, substr): if substr is None: return tags return [tag for tag in tags if substr in tag] + class PlotExecutable(Executable): """ plot executable """ @@ -55,7 +58,8 @@ def create_node(self, **kwargs): node.set_priority(1000) return node -def make_template_plot(workflow, bank_file, out_dir,bins=None, + +def make_template_plot(workflow, bank_file, out_dir, bins=None, tags=None): tags = [] if tags is None else tags makedir(out_dir) @@ -72,6 +76,7 @@ def make_template_plot(workflow, bank_file, out_dir,bins=None, workflow += node return node.output_files[0] + def make_range_plot(workflow, psd_files, out_dir, exclude=None, require=None, tags=None): tags = [] if tags is None else tags @@ -89,6 +94,7 @@ def make_range_plot(workflow, psd_files, out_dir, exclude=None, require=None, files += node.output_files return files + def make_spectrum_plot(workflow, psd_files, out_dir, tags=None, hdf_group=None, precalc_psd_files=None): tags = [] if tags is None else tags @@ -106,6 +112,7 @@ def make_spectrum_plot(workflow, psd_files, out_dir, tags=None, workflow += node return node.output_files[0] + def make_segments_plot(workflow, seg_files, out_dir, tags=None): tags = [] if tags is None else tags makedir(out_dir) @@ -115,6 +122,7 @@ def make_segments_plot(workflow, seg_files, out_dir, tags=None): node.new_output_file_opt(workflow.analysis_time, '.html', '--output-file') workflow += node + def make_gating_plot(workflow, insp_files, out_dir, tags=None): tags = [] if tags is None else tags makedir(out_dir) @@ -124,6 +132,7 @@ def make_gating_plot(workflow, insp_files, out_dir, tags=None): node.new_output_file_opt(workflow.analysis_time, '.html', '--output-file') workflow += node + def make_throughput_plot(workflow, insp_files, out_dir, tags=None): tags = [] if tags is None else tags makedir(out_dir) @@ -133,6 +142,7 @@ def make_throughput_plot(workflow, insp_files, out_dir, tags=None): node.new_output_file_opt(workflow.analysis_time, '.png', '--output-file') workflow += node + def make_foreground_table(workflow, trig_file, bank_file, out_dir, singles=None, extension='.html', tags=None, hierarchical_level=None): @@ -160,6 +170,7 @@ def make_foreground_table(workflow, trig_file, bank_file, out_dir, workflow += node return node.output_files[0] + def make_sensitivity_plot(workflow, inj_file, out_dir, exclude=None, require=None, tags=None): tags = [] if tags is None else tags @@ -177,6 +188,7 @@ def make_sensitivity_plot(workflow, inj_file, out_dir, exclude=None, files += node.output_files return files + def make_coinc_snrchi_plot(workflow, inj_file, inj_trig, stat_file, trig_file, out_dir, exclude=None, require=None, tags=None): tags = [] if tags is None else tags @@ -199,6 +211,7 @@ def make_coinc_snrchi_plot(workflow, inj_file, inj_trig, stat_file, trig_file, files += node.output_files return files + def make_inj_table(workflow, inj_file, out_dir, missed=False, singles=None, tags=None): tags = [] if tags is None else tags @@ -217,6 +230,7 @@ def make_inj_table(workflow, inj_file, out_dir, missed=False, singles=None, workflow += node return node.output_files[0] + def make_seg_table(workflow, seg_files, seg_names, out_dir, tags=None, title_text=None, description=None): """ Creates a node in the workflow for writing the segment summary @@ -242,6 +256,7 @@ def make_seg_table(workflow, seg_files, seg_names, out_dir, tags=None, workflow += node return node.output_files[0] + def make_veto_table(workflow, out_dir, vetodef_file=None, tags=None): """ Creates a node in the workflow for writing the veto_definer table. Returns a File instances for the output file. @@ -268,6 +283,7 @@ def make_veto_table(workflow, out_dir, vetodef_file=None, tags=None): workflow += node return node.output_files[0] + def make_seg_plot(workflow, seg_files, out_dir, seg_names=None, tags=None): """ Creates a node in the workflow for plotting science, and veto segments. """ @@ -286,6 +302,7 @@ def make_seg_plot(workflow, seg_files, out_dir, seg_names=None, tags=None): workflow += node return node.output_files[0] + def make_ifar_plot(workflow, trigger_file, out_dir, tags=None, hierarchical_level=None, executable='page_ifar'): """ Creates a node in the workflow for plotting cumulative histogram @@ -311,6 +328,7 @@ def make_ifar_plot(workflow, trigger_file, out_dir, tags=None, workflow += node return node.output_files[0] + def make_snrchi_plot(workflow, trig_files, veto_file, veto_name, out_dir, exclude=None, require=None, tags=None): tags = [] if tags is None else tags @@ -337,6 +355,7 @@ def make_snrchi_plot(workflow, trig_files, veto_file, veto_name, files += node.output_files return files + def make_foundmissed_plot(workflow, inj_file, out_dir, exclude=None, require=None, tags=None): if tags is None: @@ -357,6 +376,7 @@ def make_foundmissed_plot(workflow, inj_file, out_dir, exclude=None, files += node.output_files return files + def make_snrratehist_plot(workflow, bg_file, out_dir, closed_box=False, tags=None, hierarchical_level=None): @@ -384,6 +404,7 @@ def make_snrratehist_plot(workflow, bg_file, out_dir, closed_box=False, workflow += node return node.output_files[0] + def make_snrifar_plot(workflow, bg_file, out_dir, closed_box=False, cumulative=True, tags=None, hierarchical_level=None): @@ -413,6 +434,7 @@ def make_snrifar_plot(workflow, bg_file, out_dir, closed_box=False, workflow += node return node.output_files[0] + def make_results_web_page(workflow, results_dir, template='orange', explicit_dependencies=None): template_path = 'templates/'+template+'.html' @@ -428,6 +450,7 @@ def make_results_web_page(workflow, results_dir, template='orange', for dep in explicit_dependencies: workflow.add_explicit_dependancy(dep, node) + def make_single_hist(workflow, trig_file, veto_file, veto_name, out_dir, bank_file=None, exclude=None, require=None, tags=None): @@ -453,6 +476,7 @@ def make_single_hist(workflow, trig_file, veto_file, veto_name, files += node.output_files return files + def make_binned_hist(workflow, trig_file, veto_file, veto_name, out_dir, bank_file, exclude=None, require=None, tags=None): @@ -478,6 +502,7 @@ def make_binned_hist(workflow, trig_file, veto_file, veto_name, files += node.output_files return files + def make_singles_plot(workflow, trig_files, bank_file, veto_file, veto_name, out_dir, exclude=None, require=None, tags=None): tags = [] if tags is None else tags @@ -505,71 +530,18 @@ def make_singles_plot(workflow, trig_files, bank_file, veto_file, veto_name, files += node.output_files return files -def make_dq_timeseries_trigger_rate_plot(workflow, dq_files, out_dir, tags=None): - tags = [] if tags is None else tags - makedir(out_dir) - files = FileList([]) - for dq_file in dq_files: - if workflow.cp.has_option_tags('bin_trigger_rates_dq', - 'background-bins', tags=tags): - background_bins = \ - workflow.cp.get_opt_tags('bin_trigger_rates_dq', - 'background-bins', tags=tags) - bin_names = [tuple(bbin.split(':'))[0] for bbin - in background_bins.split(' ')] - else: bin_names = ['all_bin'] - for bbin in bin_names: - plot_tags = [bbin] + tags - node = PlotExecutable(workflow.cp, 'plot_dq_likelihood_vs_time', - ifos=dq_file.ifo, - out_dir=out_dir, - tags=plot_tags).create_node() - node.add_opt('--ifo', dq_file.ifo) - node.add_opt('--background-bin', bbin) - node.add_input_opt('--dq-file', dq_file) - node.new_output_file_opt(dq_file.segment, '.png', '--output-file') - workflow += node - files += node.output_files - return files -def make_dq_percentile_plot(workflow, dq_files, out_dir, tags=None): +def make_dq_flag_trigger_rate_plot(workflow, dq_file, dq_label, out_dir, tags=None): tags = [] if tags is None else tags makedir(out_dir) files = FileList([]) - for dq_file in dq_files: - if workflow.cp.has_option_tags('bin_trigger_rates_dq', - 'background-bins', tags=tags): - background_bins = \ - workflow.cp.get_opt_tags('bin_trigger_rates_dq', - 'background-bins', tags=tags) - bin_names = [tuple(bbin.split(':'))[0] for bbin - in background_bins.split(' ')] - else: bin_names = ['all_bin'] - for bbin in bin_names: - plot_tags = [bbin] + tags - node = PlotExecutable(workflow.cp, 'plot_dq_percentiles', - ifos=dq_file.ifo, - out_dir=out_dir, - tags=plot_tags).create_node() - node.add_opt('--ifo', dq_file.ifo) - node.add_opt('--background-bin', bbin) - node.add_input_opt('--dq-file', dq_file) - node.new_output_file_opt(dq_file.segment, '.png', '--output-file') - workflow += node - files += node.output_files - return files - -def make_dq_flag_trigger_rate_plot(workflow, dq_files, out_dir, tags=None): - tags = [] if tags is None else tags - makedir(out_dir) - files = FileList([]) - for dq_file in dq_files: - node = PlotExecutable(workflow.cp, 'plot_dq_flag_likelihood', - ifos=dq_file.ifo, out_dir=out_dir, - tags=tags).create_node() - node.add_input_opt('--dq-file', dq_file) - node.add_opt('--ifo', dq_file.ifo) - node.new_output_file_opt(dq_file.segment, '.png', '--output-file') - workflow += node - files += node.output_files + node = PlotExecutable(workflow.cp, 'plot_dq_flag_likelihood', + ifos=dq_file.ifo, out_dir=out_dir, + tags=tags).create_node() + node.add_input_opt('--dq-file', dq_file) + node.add_opt('--dq-label', dq_label) + node.add_opt('--ifo', dq_file.ifo) + node.new_output_file_opt(dq_file.segment, '.png', '--output-file') + workflow += node + files += node.output_files return files diff --git a/setup.py b/setup.py index 35be1165e91..ee5ccc22515 100755 --- a/setup.py +++ b/setup.py @@ -119,7 +119,7 @@ def __getattr__(self, attr): vinfo = _version_helper.generate_git_version_info() except: vinfo = vdummy() - vinfo.version = '2.3.3' + vinfo.version = '2.3.4' vinfo.release = 'True' version_script = f"""# coding: utf-8