diff --git a/.github/PULL_REQUEST_TEMPLATE/standard_template.md b/.github/PULL_REQUEST_TEMPLATE/standard_template.md new file mode 100644 index 00000000000..79a06d213f5 --- /dev/null +++ b/.github/PULL_REQUEST_TEMPLATE/standard_template.md @@ -0,0 +1,52 @@ + + + + +- [ ] The author of this pull request confirms they will adhere to the [code of conduct](https://github.com/gwastro/pycbc/blob/master/CODE_OF_CONDUCT.md) + +## Standard information about the request + + +This is a: bug fix, new feature, efficiency update, other (please describe) + + +This change affects: the offline search, the live search, inference, PyGRB + + +This change changes: documentation, result presentation / plotting, scientific output + + +This change: has appropriate unit tests, follows style guidelines (See e.g. [PEP8](https://peps.python.org/pep-0008/)), has been proposed using the [contribution guidelines](https://github.com/gwastro/pycbc/blob/master/CONTRIBUTING.md) + + +This change will: break current functionality, require additional dependencies, require a new release, other (please describe) + +## Motivation + + +## Contents + + +## Links to any issues or associated PRs + + +## Testing performed + + +## Additional notes + diff --git a/.github/workflows/inference-workflow.yml b/.github/workflows/inference-workflow.yml index 16e5b9cc65d..d95e1c5b714 100644 --- a/.github/workflows/inference-workflow.yml +++ b/.github/workflows/inference-workflow.yml @@ -25,7 +25,7 @@ jobs: wget -qO - https://download.pegasus.isi.edu/pegasus/gpg.txt | sudo apt-key add - echo "deb https://download.pegasus.isi.edu/pegasus/ubuntu bionic main" | sudo tee -a /etc/apt/sources.list sudo apt-get -o Acquire::Retries=3 update - sudo apt-get -o Acquire::Retries=3 install pegasus + sudo apt-get -o Acquire::Retries=3 install pegasus=5.0.6-1+ubuntu18 - run: sudo apt-get -o Acquire::Retries=3 install *fftw3* intel-mkl* - name: Install pycbc run: | diff --git a/.github/workflows/search-workflow.yml b/.github/workflows/search-workflow.yml index 18219411078..9383c5066cb 100644 --- a/.github/workflows/search-workflow.yml +++ b/.github/workflows/search-workflow.yml @@ -30,7 +30,7 @@ jobs: wget -qO - https://download.pegasus.isi.edu/pegasus/gpg.txt | sudo apt-key add - echo "deb https://download.pegasus.isi.edu/pegasus/ubuntu bionic main" | sudo tee -a /etc/apt/sources.list sudo apt-get -o Acquire::Retries=3 update - sudo apt-get -o Acquire::Retries=3 install pegasus + sudo apt-get -o Acquire::Retries=3 install pegasus=5.0.6-1+ubuntu18 - run: sudo apt-get -o Acquire::Retries=3 install *fftw3* intel-mkl* - name: Install pycbc run: | diff --git a/.github/workflows/tmpltbank-workflow.yml b/.github/workflows/tmpltbank-workflow.yml index bf29c9051d0..66a4f1ba3b2 100644 --- a/.github/workflows/tmpltbank-workflow.yml +++ b/.github/workflows/tmpltbank-workflow.yml @@ -29,7 +29,7 @@ jobs: wget -qO - https://download.pegasus.isi.edu/pegasus/gpg.txt | sudo apt-key add - echo "deb https://download.pegasus.isi.edu/pegasus/ubuntu bionic main" | sudo tee -a /etc/apt/sources.list sudo apt-get -o Acquire::Retries=3 update - sudo apt-get -o Acquire::Retries=3 install pegasus + sudo apt-get -o Acquire::Retries=3 install pegasus=5.0.6-1+ubuntu18 - run: sudo apt-get -o Acquire::Retries=3 install *fftw3* intel-mkl* - name: Install pycbc run: | diff --git a/.github/workflows/workflow-tests.yml b/.github/workflows/workflow-tests.yml index 754e25f8f11..86410066aff 100644 --- a/.github/workflows/workflow-tests.yml +++ b/.github/workflows/workflow-tests.yml @@ -34,7 +34,7 @@ jobs: wget -qO - https://download.pegasus.isi.edu/pegasus/gpg.txt | sudo apt-key add - echo "deb https://download.pegasus.isi.edu/pegasus/ubuntu bionic main" | sudo tee -a /etc/apt/sources.list sudo apt-get -o Acquire::Retries=3 update - sudo apt-get -o Acquire::Retries=3 install pegasus + sudo apt-get -o Acquire::Retries=3 install pegasus=5.0.6-1+ubuntu18 - run: sudo apt-get -o Acquire::Retries=3 install *fftw3* intel-mkl* - name: Install pycbc run: | diff --git a/bin/all_sky_search/pycbc_bin_trigger_rates_dq b/bin/all_sky_search/pycbc_bin_trigger_rates_dq index 25234a687f0..ef069d9a6b1 100644 --- a/bin/all_sky_search/pycbc_bin_trigger_rates_dq +++ b/bin/all_sky_search/pycbc_bin_trigger_rates_dq @@ -46,6 +46,18 @@ logging.info('Start') ifo, flag_name = args.flag_name.split(':') +if args.gating_windows: + gate_times = [] + with h5.File(args.trig_file, 'r') as trig_file: + 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') + trigs = SingleDetTriggers( args.trig_file, ifo, @@ -134,6 +146,7 @@ with h5.File(args.output_file, 'w') as f: frac_dt = abs(segs) / livetime dq_rates[state] = frac_eff / frac_dt bin_grp['dq_rates'] = dq_rates + bin_grp['num_triggers'] = len(trig_times_bin) # save dq state segments for dq_state, segs in dq_state_segs_dict.items(): @@ -142,7 +155,10 @@ with h5.File(args.output_file, 'w') as f: starts, ends = segments_to_start_end(segs) dq_grp['segment_starts'] = starts dq_grp['segment_ends'] = ends + dq_grp['livetime'] = abs(segs) f.attrs['stat'] = f'{ifo}-dq_stat_info' + f.attrs['sngl_ranking'] = args.sngl_ranking + f.attrs['sngl_ranking_threshold'] = args.stat_threshold logging.info('Done!') diff --git a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam index d4b2dac0f76..33e063203b3 100755 --- a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam +++ b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam @@ -174,18 +174,14 @@ parser = argparse.ArgumentParser(usage="", pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action=pycbc.version.Version) -parser.add_argument("--template-fit-file", +parser.add_argument("--template-fit-file", required=True, help="hdf5 file containing fit coefficients for each" " individual template. Required") -parser.add_argument("--bank-file", default=None, - help="hdf file containing template parameters. Required " - "unless reading param from template fit file") +parser.add_argument("--bank-file", required=True, + help="hdf file containing template parameters. Required") parser.add_argument("--output", required=True, help="Location for output file containing smoothed fit " - "coefficients. Required") -parser.add_argument("--use-template-fit-param", action="store_true", - help="Use parameter values stored in the template fit " - "file as template_param for smoothing.") + "coefficients. Required") parser.add_argument("--fit-param", nargs='+', help="Parameter(s) over which to regress the background " "fit coefficients. Required. Either read from " @@ -196,20 +192,19 @@ parser.add_argument("--fit-param", nargs='+', "multiple parameters, provide them as a list.") parser.add_argument("--approximant", default="SEOBNRv4", help="Approximant for template duration. Default SEOBNRv4") -parser.add_argument("--f-lower", type=float, default=0., - help="Starting frequency for calculating template " - "duration, if not reading from the template fit file") +parser.add_argument("--f-lower", type=float, + help="Start frequency for calculating template duration.") parser.add_argument("--min-duration", type=float, default=0., help="Fudge factor for templates with tiny or negative " "values of template_duration: add to duration values" " before fitting. Units seconds.") parser.add_argument("--log-param", nargs='+', - help="Take the log of the fit param before smoothing.") + help="Take the log of the fit param before smoothing. " + "Must be a list corresponding to fit params.") parser.add_argument("--smoothing-width", type=float, nargs='+', required=True, - help="Distance in the space of fit param values (or the " - "logs of them) to smooth over. Required. " - "This must be a list corresponding to the smoothing " - "parameters.") + help="Distance in the space of fit param values (or their" + " logs) to smooth over. Required. Must be a list " + "corresponding to fit params.") parser.add_argument("--smoothing-method", default="smooth_tophat", choices = _smooth_dist_func.keys(), help="Method used to smooth the fit parameters; " @@ -220,15 +215,14 @@ parser.add_argument("--smoothing-method", default="smooth_tophat", "the smoothing until 500 triggers are reached. " "'distance_weighted' weights the closest templates " "with a normal distribution of width smoothing-width " - "trucated at three smoothing-widths.") + "truncated at three smoothing-widths.") parser.add_argument("--smoothing-keywords", nargs='*', help="Keywords for the smoothing function, supplied " "as key:value pairs, e.g. total_trigs:500 to define " - "the number of templates in the n_closest smoothing " - "method") + "the number of templates for n_closest smoothing.") parser.add_argument("--output-fits-by-template", action='store_true', help="If given, will output the input file fits to " - "fit_by_template group") + "fit_by_template group.") args = parser.parse_args() if args.smoothing_keywords: diff --git a/bin/all_sky_search/pycbc_prepare_xml_for_gracedb b/bin/all_sky_search/pycbc_prepare_xml_for_gracedb index 0baf6cb19dc..f915954884a 100755 --- a/bin/all_sky_search/pycbc_prepare_xml_for_gracedb +++ b/bin/all_sky_search/pycbc_prepare_xml_for_gracedb @@ -96,7 +96,7 @@ def read_psds(psd_files): 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 + return psds psds = read_psds(args.psd_files) diff --git a/bin/bank/pycbc_aligned_bank_cat b/bin/bank/pycbc_aligned_bank_cat index e09b6f274d9..f564c3b9bad 100644 --- a/bin/bank/pycbc_aligned_bank_cat +++ b/bin/bank/pycbc_aligned_bank_cat @@ -48,9 +48,8 @@ __program__ = "pycbc_aligned_bank_cat" parser = argparse.ArgumentParser(description=__doc__, formatter_class=tmpltbank.IndentedHelpFormatterWithNL) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=__version__) -parser.add_argument("--verbose", action="store_true", default=False, - help="verbose output") parser.add_argument("-i", "--input-glob", help="file glob the list of paramters") parser.add_argument("-I", "--input-files", nargs='+', @@ -76,11 +75,7 @@ tmpltbank.insert_base_bank_options(parser, match_req=False) options = parser.parse_args() -if options.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -logging.basicConfig(format='%(asctime)s %(message)s', level=log_level) +pycbc.init_logging(options.verbose) # Sanity check options if not options.output_file: diff --git a/bin/bank/pycbc_aligned_stoch_bank b/bin/bank/pycbc_aligned_stoch_bank index 9b9a473d3fa..e615fd75415 100644 --- a/bin/bank/pycbc_aligned_stoch_bank +++ b/bin/bank/pycbc_aligned_stoch_bank @@ -23,6 +23,7 @@ Stochastic aligned spin bank generator. import argparse import numpy import logging + import pycbc import pycbc.version from pycbc import tmpltbank @@ -45,8 +46,7 @@ parser = argparse.ArgumentParser(description=_desc, # Begin with code specific options parser.add_argument("--version", action="version", version=__version__) -parser.add_argument("--verbose", action="store_true", default=False, - help="verbose output") +pycbc.add_common_pycbc_options(parser) parser.add_argument("-V", "--vary-fupper", action="store_true", default=False, help="Use a variable upper frequency cutoff in laying " "out the bank. OPTIONAL.") @@ -96,12 +96,7 @@ tmpltbank.insert_ethinca_metric_options(parser) opts = parser.parse_args() -if opts.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -log_format='%(asctime)s %(message)s' -logging.basicConfig(format=log_format, level=log_level) +pycbc.init_logging(opts.verbose) # delete defaults for redundant options if not varying fupper if not opts.vary_fupper: diff --git a/bin/bank/pycbc_bank_verification b/bin/bank/pycbc_bank_verification index 9c76e28aaf5..6758daf9556 100644 --- a/bin/bank/pycbc_bank_verification +++ b/bin/bank/pycbc_bank_verification @@ -26,14 +26,16 @@ import argparse import logging import numpy import h5py +import matplotlib +matplotlib.use('Agg') +import pylab + from ligo.lw import lsctables, utils as ligolw_utils + import pycbc import pycbc.version from pycbc import tmpltbank, psd, strain from pycbc.io.ligolw import LIGOLWContentHandler -import matplotlib -matplotlib.use('Agg') -import pylab __author__ = "Ian Harry " @@ -48,8 +50,7 @@ parser = argparse.ArgumentParser(description=__doc__, # Begin with code specific options parser.add_argument("--version", action="version", version=__version__) -parser.add_argument("--verbose", action="store_true", default=False, - help="verbose output") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--histogram-output-file", action="store", default=None, help="Output a histogram of fitting factors to the " "supplied file. If not given no histogram is produced.") @@ -111,11 +112,7 @@ pycbc.strain.insert_strain_option_group(parser) opts = parser.parse_args() -if opts.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -logging.basicConfig(format='%(asctime)s %(message)s', level=log_level) +pycbc.init_logging(opts.verbose) # Sanity check options tmpltbank.verify_metric_calculation_options(opts, parser) diff --git a/bin/bank/pycbc_brute_bank b/bin/bank/pycbc_brute_bank index e887be710a1..583793b43c9 100644 --- a/bin/bank/pycbc_brute_bank +++ b/bin/bank/pycbc_brute_bank @@ -19,16 +19,21 @@ """Generate a bank of templates using a brute force stochastic method. """ -import numpy, h5py, logging, argparse, numpy.random +import numpy +import h5py +import logging +import argparse +import numpy.random +from scipy.stats import gaussian_kde + import pycbc.waveform, pycbc.filter, pycbc.types, pycbc.psd, pycbc.fft, pycbc.conversions from pycbc import transforms from pycbc.waveform.spa_tmplt import spa_length_in_time from pycbc.distributions import read_params_from_config from pycbc.distributions.utils import draw_samples_from_config, prior_from_config -from scipy.stats import gaussian_kde parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument('--verbose', action='store_true') +pycbc.add_common_pycbc_options(parser) parser.add_argument('--output-file', required=True, help='Output file name for template bank.') parser.add_argument('--input-file', @@ -71,10 +76,10 @@ parser.add_argument('--tau0-end', type=float) parser.add_argument('--tau0-cutoff-frequency', type=float, default=15.0) pycbc.psd.insert_psd_option_group(parser) args = parser.parse_args() + pycbc.init_logging(args.verbose) numpy.random.seed(args.seed) - # Read the .ini file if it's in the input. if args.input_config is not None: config_parser = pycbc.types.config.InterpolatingConfigParser() diff --git a/bin/bank/pycbc_coinc_bank2hdf b/bin/bank/pycbc_coinc_bank2hdf index fb03caea604..db35f48d4be 100644 --- a/bin/bank/pycbc_coinc_bank2hdf +++ b/bin/bank/pycbc_coinc_bank2hdf @@ -23,6 +23,7 @@ with their template. import argparse import logging import numpy + import pycbc from pycbc.waveform import bank @@ -56,6 +57,7 @@ def parse_parameters(parameters): parser = argparse.ArgumentParser() parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--bank-file', required=True, help="The bank file to load. Must end in '.xml[.gz]' " "and must contain a SnglInspiral table or must end " @@ -82,8 +84,6 @@ bank.add_approximant_arg(parser, parser.add_argument("--force", action="store_true", default=False, help="Overwrite the given hdf file if it exists. " "Otherwise, an error is raised.") -parser.add_argument("--verbose", action="store_true", default=False, - help="Be verbose.") args = parser.parse_args() pycbc.init_logging(args.verbose) diff --git a/bin/bank/pycbc_geom_aligned_2dstack b/bin/bank/pycbc_geom_aligned_2dstack index 179c74f8216..215440e76c1 100644 --- a/bin/bank/pycbc_geom_aligned_2dstack +++ b/bin/bank/pycbc_geom_aligned_2dstack @@ -31,6 +31,7 @@ import copy import numpy import logging import h5py + import pycbc.tmpltbank import pycbc.version from pycbc import pnutils @@ -47,10 +48,9 @@ usage = """usage: %prog [options]""" _desc = __doc__[1:] parser = argparse.ArgumentParser(usage, description=_desc, formatter_class=pycbc.tmpltbank.IndentedHelpFormatterWithNL) +pycbc.add_common_pycbc_options(parser) # Code specific options parser.add_argument('--version', action='version', version=__version__) -parser.add_argument("--verbose", action="store_true", default=False,\ - help="verbose output") parser.add_argument("--pn-order", action="store", type=str,\ default=None,\ help="Determines the PN order to use. Note that if you "+\ @@ -106,11 +106,7 @@ opts = parser.parse_args() opts.eval_vec4_depth = not opts.skip_vec4_depth opts.eval_vec5_depth = not opts.skip_vec5_depth -if opts.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -logging.basicConfig(format='%(asctime)s %(message)s', level=log_level) +pycbc.init_logging(opts.verbose) # Sanity check options if not opts.pn_order: diff --git a/bin/bank/pycbc_geom_aligned_bank b/bin/bank/pycbc_geom_aligned_bank index 97da6a5674c..d18c3490f91 100644 --- a/bin/bank/pycbc_geom_aligned_bank +++ b/bin/bank/pycbc_geom_aligned_bank @@ -31,8 +31,10 @@ import logging import h5py import configparser from scipy import spatial + from ligo.lw import ligolw from ligo.lw import utils as ligolw_utils + import pycbc import pycbc.psd import pycbc.strain @@ -175,9 +177,8 @@ parser = argparse.ArgumentParser(description=_desc, formatter_class=pycbc.tmpltbank.IndentedHelpFormatterWithNL) # Begin with code specific options +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=__version__) -parser.add_argument("-V", "--verbose", action="store_true", default=False, - help="verbose output") parser.add_argument("-s", "--stack-distance", action="store", type=float,\ default=0.2, help="Minimum metric spacing before we "+\ "stack.") @@ -249,12 +250,7 @@ outdoc.appendChild(ligolw.LIGO_LW()) create_process_table(outdoc, __program__, options=vars(opts)) ligolw_utils.write_filename(outdoc, opts.metadata_file) -if opts.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -log_format='%(asctime)s %(message)s' -logging.basicConfig(format=log_format, level=log_level) +pycbc.init_logging(opts.verbose) opts.max_mismatch = 1 - opts.min_match diff --git a/bin/bank/pycbc_geom_nonspinbank b/bin/bank/pycbc_geom_nonspinbank index 79298eff547..dcf5aa71890 100644 --- a/bin/bank/pycbc_geom_nonspinbank +++ b/bin/bank/pycbc_geom_nonspinbank @@ -25,6 +25,7 @@ import math import copy import numpy import logging + import pycbc import pycbc.version import pycbc.psd @@ -48,8 +49,7 @@ parser = argparse.ArgumentParser( # Begin with code specific options parser.add_argument('--version', action='version', version=__version__) -parser.add_argument("-V", "--verbose", action="store_true", default=False, - help="verbose output") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--random-seed", action="store", type=int, default=None, help="""Random seed to use when calling numpy.random @@ -78,11 +78,7 @@ tmpltbank.insert_ethinca_metric_options(parser) opts = parser.parse_args() -if opts.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -logging.basicConfig(format='%(asctime)s %(message)s', level=log_level) +pycbc.init_logging(opts.verbose) opts.max_mismatch = 1 - opts.min_match tmpltbank.verify_metric_calculation_options(opts, parser) diff --git a/bin/bank/pycbc_tmpltbank_to_chi_params b/bin/bank/pycbc_tmpltbank_to_chi_params index df0201ccc0f..ad7429b5f86 100644 --- a/bin/bank/pycbc_tmpltbank_to_chi_params +++ b/bin/bank/pycbc_tmpltbank_to_chi_params @@ -33,7 +33,9 @@ import argparse import logging import numpy import h5py + from ligo.lw import lsctables, utils as ligolw_utils + from pycbc import tmpltbank, psd, strain from pycbc.io.ligolw import LIGOLWContentHandler @@ -44,8 +46,7 @@ parser = argparse.ArgumentParser(description=__doc__, # Begin with code specific options parser.add_argument("--version", action="version", version=__version__) -parser.add_argument("--verbose", action="store_true", default=False, - help="verbose output") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--input-bank", action="store", required=True, help="The template bank to use an input.") parser.add_argument("--output-file", action="store", required=True, @@ -88,11 +89,7 @@ pycbc.strain.insert_strain_option_group(parser) opts = parser.parse_args() -if opts.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -logging.basicConfig(format='%(asctime)s %(message)s', level=log_level) +pycbc.init_logging(opts.verbose) # Sanity check options tmpltbank.verify_metric_calculation_options(opts, parser) diff --git a/bin/inference/pycbc_inference b/bin/inference/pycbc_inference index a4c8dd50649..d0cfd7b953a 100644 --- a/bin/inference/pycbc_inference +++ b/bin/inference/pycbc_inference @@ -40,10 +40,9 @@ from pycbc.workflow import configuration # command line usage parser = argparse.ArgumentParser(usage=__file__ + " [--options]", description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=__version__, help="Prints version information.") -parser.add_argument("--verbose", action="store_true", default=False, - help="Print logging messages.") # output options parser.add_argument("--output-file", type=str, required=True, help="Output file path.") diff --git a/bin/inference/pycbc_inference_extract_samples b/bin/inference/pycbc_inference_extract_samples index a404c6b5846..ffbf89ee996 100644 --- a/bin/inference/pycbc_inference_extract_samples +++ b/bin/inference/pycbc_inference_extract_samples @@ -71,7 +71,7 @@ def isthesame(current_val, val): parser = ResultsArgumentParser(defaultparams='all', autoparamlabels=False, description=__doc__) - +pycbc.add_common_pycbc_options(parser) parser.add_argument("--output-file", type=str, required=True, help="Output file to create.") parser.add_argument("--force", action="store_true", default=False, @@ -85,8 +85,6 @@ parser.add_argument("--skip-groups", default=None, nargs="+", "to write all groups if only one file is provided, " "and all groups from the first file except " "sampler_info if multiple files are provided.") -parser.add_argument("--verbose", action="store_true", default=False, - help="Be verbose") opts = parser.parse_args() diff --git a/bin/inference/pycbc_inference_model_stats b/bin/inference/pycbc_inference_model_stats index 8a8512cb4b8..9abf8c6dba8 100644 --- a/bin/inference/pycbc_inference_model_stats +++ b/bin/inference/pycbc_inference_model_stats @@ -27,6 +27,7 @@ import shutil import logging import time import numpy +import tqdm import pycbc from pycbc.pool import (use_mpi, choose_pool) @@ -37,6 +38,7 @@ from pycbc.workflow import WorkflowConfigParser parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--input-file", type=str, required=True, help="Input HDF file to read. Can be either an inference " "file or a posterior file.") @@ -51,8 +53,6 @@ parser.add_argument("--nprocesses", type=int, default=1, parser.add_argument("--config-file", nargs="+", type=str, default=None, help="Override the config file stored in the input file " "with the given file(s).") -parser.add_argument("--verbose", action="count", default=0, - help="Print logging messages.") parser.add_argument('--reconstruct-parameters', action="store_true", help="Reconstruct marginalized parameters") @@ -117,7 +117,8 @@ logging.info('Getting samples') with loadfile(opts.input_file, 'r') as fp: # we'll need the shape; all the arrays in the samples group should have the # same shape - shape = fp[fp.samples_group]['loglikelihood'].shape + pick = list(fp[fp.samples_group].keys())[0] + shape = fp[fp.samples_group][pick].shape samples = {} for p in model.variable_params: samples[p] = fp[fp.samples_group][p][()].flatten() @@ -128,7 +129,8 @@ logging.info("Loaded %i samples", samples.size) # get the stats logging.info("Calculating stats") -data = list(pool.map(model_call, enumerate(samples))) +data = list(tqdm.tqdm(pool.imap(model_call, enumerate(samples)), + total=len(samples))) stats = [x[0] for x in data] rec = [x[1] for x in data] diff --git a/bin/inference/pycbc_inference_monitor b/bin/inference/pycbc_inference_monitor index b4311762ca2..090bdf6a9ff 100644 --- a/bin/inference/pycbc_inference_monitor +++ b/bin/inference/pycbc_inference_monitor @@ -6,13 +6,13 @@ import os, sys, time, glob, shutil, logging, argparse, pycbc import pycbc.workflow as wf parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument('--inference-file', help="The name of the inference file") parser.add_argument('--check-interval', default=10, type=int, help="Polling interval to check for file changes (s)") parser.add_argument('--output-dir', help="Output plots / results directory") parser.add_argument('--allow-failure', action='store_true', help="Allow for a failure in plot generation") -parser.add_argument('--verbose', action='store_true') wf.add_workflow_command_line_group(parser) args = parser.parse_args() diff --git a/bin/inference/pycbc_inference_plot_acceptance_rate b/bin/inference/pycbc_inference_plot_acceptance_rate index b3f1e9aafac..fe0b81555b4 100644 --- a/bin/inference/pycbc_inference_plot_acceptance_rate +++ b/bin/inference/pycbc_inference_plot_acceptance_rate @@ -34,6 +34,7 @@ parser = argparse.ArgumentParser( usage="pycbc_inference_plot_acceptance_rate [--options]", description="Plots histogram of the fractions of steps " "accepted by walkers.") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--input-file", type=str, required=True, help="Path to input HDF file.") parser.add_argument('--version', action='version', version=__version__, @@ -56,9 +57,6 @@ parser.add_argument("--temps", type=int, nargs="+", default=None, # add number of bins for histogram parser.add_argument("--bins", type=int, default=10, help="Specify number of bins for the histogram plot.") -# verbose option -parser.add_argument("--verbose", action="store_true", default=False, - help="") # parse the command line opts = parser.parse_args() diff --git a/bin/inference/pycbc_inference_plot_acf b/bin/inference/pycbc_inference_plot_acf index d38645be49c..41ab706395a 100644 --- a/bin/inference/pycbc_inference_plot_acf +++ b/bin/inference/pycbc_inference_plot_acf @@ -38,9 +38,7 @@ from pycbc.inference.sampler import samplers parser = io.ResultsArgumentParser(skip_args='thin-interval', description="Plots autocorrelation function " "from inference samples.") -# verbose option -parser.add_argument("--verbose", action="store_true", default=False, - help="Print logging info.") +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=__version__, help='show version number and exit') # output plot options diff --git a/bin/inference/pycbc_inference_plot_acl b/bin/inference/pycbc_inference_plot_acl index 348dbb10763..860a9810713 100644 --- a/bin/inference/pycbc_inference_plot_acl +++ b/bin/inference/pycbc_inference_plot_acl @@ -39,8 +39,7 @@ parser = io.ResultsArgumentParser(skip_args=['thin-interval', 'temps'], description="Histograms autocorrelation " "length per walker from an MCMC " "sampler.") -parser.add_argument("--verbose", action="store_true", default=False, - help="Print logging info.") +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=__version__, help='show version number and exit') # output plot options diff --git a/bin/inference/pycbc_inference_plot_dynesty_run b/bin/inference/pycbc_inference_plot_dynesty_run index 04a1aaf0371..6e2c8d03afb 100644 --- a/bin/inference/pycbc_inference_plot_dynesty_run +++ b/bin/inference/pycbc_inference_plot_dynesty_run @@ -34,6 +34,7 @@ parser = argparse.ArgumentParser( usage="pycbc_inference_plot_dynesty_run [--options]", description="Plots various figures showing evolution of " "dynesty run.") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--input-file", type=str, required=True, help="Path to input HDF file.") parser.add_argument('--version', action='version', version=__version__, @@ -42,9 +43,6 @@ parser.add_argument('--version', action='version', version=__version__, parser.add_argument("--output-file", type=str, required=True, help="Path to output plot.") -parser.add_argument("--verbose", action="store_true", default=False, - help="") - # parse the command line opts = parser.parse_args() diff --git a/bin/inference/pycbc_inference_plot_dynesty_traceplot b/bin/inference/pycbc_inference_plot_dynesty_traceplot index c328313ec94..c7e352f91b7 100644 --- a/bin/inference/pycbc_inference_plot_dynesty_traceplot +++ b/bin/inference/pycbc_inference_plot_dynesty_traceplot @@ -36,6 +36,7 @@ parser = argparse.ArgumentParser( description="Plots various figures showing evolution of " "dynesty run.") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--input-file", type=str, required=True, help="Path to input HDF file.") parser.add_argument('--version', action='version', version=__version__, @@ -44,8 +45,6 @@ parser.add_argument('--version', action='version', version=__version__, parser.add_argument("--output-file", type=str, required=True, help="Path to output plot.") -parser.add_argument("--verbose", action="store_true", default=False, - help="") # parse the command line opts = parser.parse_args() diff --git a/bin/inference/pycbc_inference_plot_geweke b/bin/inference/pycbc_inference_plot_geweke index 16c04463932..755b2f3268a 100644 --- a/bin/inference/pycbc_inference_plot_geweke +++ b/bin/inference/pycbc_inference_plot_geweke @@ -32,11 +32,9 @@ from pycbc.inference import (io, geweke, option_utils) # add options to command line parser = io.ResultsArgumentParser(skip_args=['walkers']) -# program-specific +pycbc.add_common_pycbc_options(parser) -# verbose option -parser.add_argument("--verbose", action="store_true", default=False, - help="Print logging info.") +# program-specific parser.add_argument('--version', action='version', version=__version__, help='show version number and exit') diff --git a/bin/inference/pycbc_inference_plot_inj_recovery b/bin/inference/pycbc_inference_plot_inj_recovery index b65dd77c4ac..fc51655bdd8 100644 --- a/bin/inference/pycbc_inference_plot_inj_recovery +++ b/bin/inference/pycbc_inference_plot_inj_recovery @@ -21,12 +21,11 @@ from pycbc.results import save_fig_with_metadata # parse command line parser = io.ResultsArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=__version__, help="Prints version information.") parser.add_argument("--output-file", required=True, type=str, help="Path to save output plot.") -parser.add_argument("--verbose", action="store_true", - help="Allows print statements.") parser.add_argument("--percentiles", nargs=2, type=float, default=[5, 95], help="Percentiles to use as limits.") option_utils.add_scatter_option_group(parser) diff --git a/bin/inference/pycbc_inference_plot_movie b/bin/inference/pycbc_inference_plot_movie index edc04be8c00..1f3a4b83c88 100644 --- a/bin/inference/pycbc_inference_plot_movie +++ b/bin/inference/pycbc_inference_plot_movie @@ -110,6 +110,7 @@ def integer_logspace(start, end, num): # the frame number/step options skip_args = ['thin-start', 'thin-interval', 'thin-end', 'iteration'] parser = io.ResultsArgumentParser(description=__doc__, skip_args=skip_args) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=__version__, help="show version number and exit") # make frame number and frame step mutually exclusive @@ -135,7 +136,6 @@ parser.add_argument("--log-steps", action="store_true", default=False, parser.add_argument("--output-prefix", type=str, required=True, help="Output path and prefix for the frame files " "(without extension).") -parser.add_argument('--verbose', action='store_true') parser.add_argument('--dpi', type=int, default=200, help='Set the dpi for each frame; default is 200') parser.add_argument("--nprocesses", type=int, default=None, diff --git a/bin/inference/pycbc_inference_plot_posterior b/bin/inference/pycbc_inference_plot_posterior index e9df4e79048..7d507a7ab92 100644 --- a/bin/inference/pycbc_inference_plot_posterior +++ b/bin/inference/pycbc_inference_plot_posterior @@ -52,13 +52,12 @@ use('agg') # add options to command line parser = io.ResultsArgumentParser() +pycbc.add_common_pycbc_options(parser) # program-specific parser.add_argument("--version", action="version", version=__version__, help="Prints version information.") parser.add_argument("--output-file", type=str, required=True, help="Output plot path.") -parser.add_argument("--verbose", action="store_true", default=False, - help="Be verbose") parser.add_argument("--plot-prior", nargs="+", type=str, help="Plot the prior on the 1D marginal plots using the " "given config file(s).") diff --git a/bin/inference/pycbc_inference_plot_pp b/bin/inference/pycbc_inference_plot_pp index f976a1f5317..5e6490499c8 100644 --- a/bin/inference/pycbc_inference_plot_pp +++ b/bin/inference/pycbc_inference_plot_pp @@ -42,12 +42,11 @@ from pycbc import __version__ # parse command line parser = io.ResultsArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=__version__, help='show version number and exit') parser.add_argument("--output-file", required=True, type=str, help="Path to save output plot.") -parser.add_argument("--verbose", action="store_true", - help="Allows print statements.") parser.add_argument("--injection-hdf-group", default="injections", help="HDF group that contains injection values. " "Default is 'injections'.") diff --git a/bin/inference/pycbc_inference_plot_prior b/bin/inference/pycbc_inference_plot_prior index 73d39b2d481..436d0c6f0e6 100644 --- a/bin/inference/pycbc_inference_plot_prior +++ b/bin/inference/pycbc_inference_plot_prior @@ -46,6 +46,7 @@ def cartesian(arrays): parser = argparse.ArgumentParser( usage="pycbc_inference_plot_prior [--options]", description="Plots prior distributions.") +pycbc.add_common_pycbc_options(parser) # add input options parser.add_argument("--config-files", type=str, nargs="+", required=True, help="A file parsable by " @@ -75,9 +76,6 @@ parser.add_argument("--nsamples", type=int, default=10000, "plotting. Default is 10000.") parser.add_argument("--output-file", type=str, required=True, help="Path to output plot.") -# verbose option -parser.add_argument("--verbose", action="store_true", default=False, - help="") parser.add_argument("--version", action="version", version=__version__, help="show version number and exit") diff --git a/bin/inference/pycbc_inference_plot_samples b/bin/inference/pycbc_inference_plot_samples index 4ece1a64ac0..651c79e5178 100644 --- a/bin/inference/pycbc_inference_plot_samples +++ b/bin/inference/pycbc_inference_plot_samples @@ -34,8 +34,7 @@ import sys # command line usage parser = argparse.parser = io.ResultsArgumentParser( skip_args=['chains', 'iteration']) -parser.add_argument("--verbose", action="store_true", default=False, - help="Print logging info.") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=__version__, help="show version number and exit") parser.add_argument("--chains", nargs='+', default=None, diff --git a/bin/inference/pycbc_inference_pp_table_summary b/bin/inference/pycbc_inference_pp_table_summary index 25d41539d79..23195eb8495 100644 --- a/bin/inference/pycbc_inference_pp_table_summary +++ b/bin/inference/pycbc_inference_pp_table_summary @@ -31,8 +31,7 @@ from pycbc.inference.io import (ResultsArgumentParser, results_from_cli, injections_from_cli) parser = ResultsArgumentParser(description=__doc__) -parser.add_argument("--verbose", action="store_true", default=False, - help="Print logging info.") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--output-file", type=str, required=True, help="Path to output plot.") add_injsamples_map_opt(parser) diff --git a/bin/inference/pycbc_inference_table_summary b/bin/inference/pycbc_inference_table_summary index 261715e35da..6e29fb636ff 100644 --- a/bin/inference/pycbc_inference_table_summary +++ b/bin/inference/pycbc_inference_table_summary @@ -29,8 +29,7 @@ from pycbc.inference.io import ResultsArgumentParser, results_from_cli parser = ResultsArgumentParser( description="Makes a table of posterior results.") -parser.add_argument("--verbose", action="store_true", default=False, - help="Print logging info.") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--output-file", type=str, required=True, help="Path to output plot.") parser.add_argument("--print-metadata", nargs="+", metavar="PARAM[:LABEL]", diff --git a/bin/live/pycbc_live_combine_single_fits b/bin/live/pycbc_live_combine_single_fits index a0638ecb036..c2847a1e774 100644 --- a/bin/live/pycbc_live_combine_single_fits +++ b/bin/live/pycbc_live_combine_single_fits @@ -84,10 +84,11 @@ for f in args.trfits_files: same_conf = (fits_f.attrs['sngl_ranking'] == sngl_rank and fits_f.attrs['fit_threshold'] == fit_thresh and fits_f.attrs['fit_function'] == fit_func + and fits_f['bins_lower'].size == bl.size and all(fits_f['bins_lower'][:] == bl) and all(fits_f['bins_upper'][:] == bu)) if not same_conf: - logging.warn( + logging.warning( "Found a change in the fit configuration, skipping %s", f ) @@ -119,8 +120,8 @@ for f in args.trfits_files: counts_all[ifo].append(ffi['counts'][:]) alphas_all[ifo].append(ffi['fit_coeff'][:]) if any(np.isnan(ffi['fit_coeff'][:])): - logging.warn("nan in %s, %s", f, ifo) - logging.warn(ffi['fit_coeff'][:]) + logging.warning("nan in %s, %s", f, ifo) + logging.warning(ffi['fit_coeff'][:]) # Set up the date array, this is stored as an offset from the first trigger time of # the first file to the last trigger of the file diff --git a/bin/plotting/pycbc_page_dq_table b/bin/plotting/pycbc_page_dq_table new file mode 100644 index 00000000000..83b39b01569 --- /dev/null +++ b/bin/plotting/pycbc_page_dq_table @@ -0,0 +1,54 @@ +#!/usr/bin/env python +""" Make a table of dq state information +""" +import sys +import argparse +import h5py as h5 +import numpy as np + +import pycbc +import pycbc.results +from pycbc.version import git_verbose_msg as version + +parser = argparse.ArgumentParser() +parser.add_argument("--version", action="version", version=version) +parser.add_argument('--ifo', required=True) +parser.add_argument('--dq-file', required=True) +parser.add_argument('--verbose', action='count') +parser.add_argument('--output-file') +args = parser.parse_args() + +dq_states = { + 'dq_state_0': 'Clean', + 'dq_state_1': 'DQ Flag', + 'dq_state_2': 'Autogating', +} + +f = h5.File(args.dq_file, 'r') +grp = f[args.ifo]['dq_segments'] + +livetimes = [] +total_livetime = 0 +for dq_state in dq_states: + livetime = grp[dq_state]['livetime'][()] + livetimes.append(livetime) + total_livetime += livetime +livetimes.append(total_livetime) + +frac_livetimes = [lt / total_livetime for lt in livetimes] +state_names = list(dq_states.values()) + ['Total'] +columns = [state_names, livetimes, frac_livetimes] +columns = [np.array(c) for c in columns] +col_names = ['DQ State', 'Livetime', '% of Livetime'] + +format_strings = [None, '0.0', '0.00%'] + +html_table = pycbc.results.html_table(columns, col_names, + page_size=len(state_names), + format_strings=format_strings) +title = f'{args.ifo} DQ State Livetimes' +caption = 'Table of DQ state livetimes' + +pycbc.results.save_fig_with_metadata( + str(html_table), args.output_file, title=title, + caption=caption, cmd=' '.join(sys.argv)) diff --git a/bin/plotting/pycbc_page_template_bin_table b/bin/plotting/pycbc_page_template_bin_table new file mode 100644 index 00000000000..ef2f7ffb2aa --- /dev/null +++ b/bin/plotting/pycbc_page_template_bin_table @@ -0,0 +1,75 @@ +#!/usr/bin/env python +""" Make a table of template bin information +""" +import sys +import argparse +import h5py as h5 +import numpy as np + +import pycbc +import pycbc.results +from pycbc.version import git_verbose_msg as version + +parser = argparse.ArgumentParser() +parser.add_argument("--version", action="version", version=version) +parser.add_argument('--ifo', required=True) +parser.add_argument('--dq-file', required=True) +parser.add_argument('--verbose', action='count') +parser.add_argument('--output-file') +args = parser.parse_args() + +f = h5.File(args.dq_file, 'r') +grp = f[args.ifo]['bins'] +bin_names = list(grp.keys()) + +sngl_ranking = f.attrs['sngl_ranking'] +sngl_thresh = f.attrs['sngl_ranking_threshold'] + +livetime = 0 +seg_grp = f[args.ifo]['dq_segments'] +for k in seg_grp.keys(): + livetime += seg_grp[k]['livetime'][()] + +num_templates = [] +num_triggers = [] +total_templates = 0 +total_triggers = 0 +for bin_name in bin_names: + bin_grp = grp[bin_name] + + n_tmp = len(bin_grp['tids'][:]) + num_templates.append(n_tmp) + total_templates += n_tmp + + n_trig = bin_grp['num_triggers'][()] + num_triggers.append(n_trig) + total_triggers += n_trig + +bin_names.append('Total') +num_triggers.append(total_triggers) +num_templates.append(total_templates) +frac_triggers = [n / total_triggers for n in num_triggers] +frac_templates = [n / total_templates for n in num_templates] +trigger_rate = [n / livetime for n in num_triggers] + +col_names = ['Template Bin', 'Number of Templates', '% of Templates', + 'Number of Loud Triggers', '% of Loud Triggers', + 'Loud Trigger Rate (Hz)'] +columns = [bin_names, num_templates, frac_templates, + num_triggers, frac_triggers, trigger_rate] +columns = [np.array(c) for c in columns] + +format_strings = [None, '#', '0.000%', '#', '0.0%', '0.00E0'] + +html_table = pycbc.results.html_table(columns, col_names, + page_size=len(bin_names), + format_strings=format_strings) +title = f'{args.ifo} DQ Template Bin Information' +caption = 'Table of information about template bins ' \ + + 'used for DQ trigger rate calculations. ' \ + + 'Loud triggers are defined as those with ' \ + + f'{sngl_ranking} > {sngl_thresh}.' + +pycbc.results.save_fig_with_metadata( + str(html_table), args.output_file, title=title, + caption=caption, cmd=' '.join(sys.argv)) diff --git a/bin/plotting/pycbc_plot_dq_flag_likelihood b/bin/plotting/pycbc_plot_dq_flag_likelihood index a1d7beed5c1..ed4a6a1231f 100644 --- a/bin/plotting/pycbc_plot_dq_flag_likelihood +++ b/bin/plotting/pycbc_plot_dq_flag_likelihood @@ -77,7 +77,10 @@ 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 Trigger Rates' +plot_title = f'{ifo}:{args.dq_label} DQ Trigger Rates' + +ax.set_title(plot_title) + plot_caption = 'The log likelihood correction \ during during each dq state for each template bin.' pycbc.results.save_fig_with_metadata( diff --git a/bin/pycbc_banksim b/bin/pycbc_banksim index a7598f43226..66cf7e9f531 100644 --- a/bin/pycbc_banksim +++ b/bin/pycbc_banksim @@ -21,8 +21,11 @@ import logging from tqdm import tqdm from numpy import complex64, array from argparse import ArgumentParser +from math import ceil, log + from ligo.lw import utils as ligolw_utils from ligo.lw import lsctables + from pycbc.pnutils import mass1_mass2_to_mchirp_eta from pycbc.pnutils import mass1_mass2_to_tau0_tau3 from pycbc.pnutils import nearest_larger_binary_number @@ -32,7 +35,6 @@ from pycbc import DYN_RANGE_FAC from pycbc.types import FrequencySeries, TimeSeries, zeros, complex_same_precision_as from pycbc.filter import match, sigmasq from pycbc.io.ligolw import LIGOLWContentHandler -from math import ceil, log import pycbc.psd, pycbc.scheme, pycbc.fft, pycbc.strain, pycbc.version from pycbc.detector import overhead_antenna_pattern as generate_fplus_fcross from pycbc.waveform import TemplateBank @@ -144,8 +146,7 @@ aprs = sorted(list(set(td_approximants() + fd_approximants()))) parser = ArgumentParser(description=__doc__) parser.add_argument("--match-file", dest="out_file", metavar="FILE", required=True, help="File to output match results") -parser.add_argument("--verbose", action='store_true', default=False, - help="Print verbose statements") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) diff --git a/bin/pycbc_banksim_combine_banks b/bin/pycbc_banksim_combine_banks index 3c3570b27f4..56f382f19ee 100644 --- a/bin/pycbc_banksim_combine_banks +++ b/bin/pycbc_banksim_combine_banks @@ -26,6 +26,7 @@ from os.path import isfile import argparse import logging from numpy import * + import pycbc import pycbc.version @@ -37,8 +38,7 @@ _desc = __doc__[1:] parser = argparse.ArgumentParser(description=_desc) parser.add_argument('--version', action=pycbc.version.Version) -parser.add_argument("--verbose", action="store_true", default=False, - help="verbose output") +pycbc.add_common_pycbc_options(parser) parser.add_argument("-I", "--input-files", nargs='+', help="Explicit list of input files.") parser.add_argument("-o", "--output-file", required=True, @@ -46,6 +46,8 @@ parser.add_argument("-o", "--output-file", required=True, options = parser.parse_args() +pycbc.init_logging(options.verbose) + dtypef = dtype([('match', float64), ('bank', unicode_, 256), ('bank_i', int32), ('sim', unicode_, 256), ('sim_i', int32), ('sigmasq', float64)]) diff --git a/bin/pycbc_banksim_match_combine b/bin/pycbc_banksim_match_combine index f2008744887..1692504e4c1 100644 --- a/bin/pycbc_banksim_match_combine +++ b/bin/pycbc_banksim_match_combine @@ -26,7 +26,9 @@ import imp import argparse import numpy as np import h5py + from ligo.lw import utils, table + import pycbc from pycbc import pnutils from pycbc.waveform import TemplateBank @@ -43,8 +45,7 @@ __program__ = "pycbc_banksim_match_combine" parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--version", action="version", version=__version__) -parser.add_argument("--verbose", action="store_true", default=False, - help="verbose output") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--match-files", nargs='+', help="Explicit list of match files.") #parser.add_argument("--inj-files", nargs='+', @@ -62,6 +63,8 @@ parser.add_argument("--filter-func-file", default=None, "as numpy arrays.") options = parser.parse_args() +pycbc.init_logging(options.verbose) + dtypem = np.dtype([('match', np.float64), ('bank', np.unicode_, 256), ('bank_i', np.int32), ('sim', np.unicode_, 256), ('sim_i', np.int32), ('sigmasq', np.float64)]) diff --git a/bin/pycbc_banksim_skymax b/bin/pycbc_banksim_skymax index 40a208b9132..91e0973ca4f 100644 --- a/bin/pycbc_banksim_skymax +++ b/bin/pycbc_banksim_skymax @@ -21,8 +21,12 @@ import sys import logging from numpy import complex64, float32, sqrt, argmax, real, array from argparse import ArgumentParser +from math import ceil, log +from tqdm import tqdm + from ligo.lw import utils as ligolw_utils from ligo.lw import lsctables + from pycbc.pnutils import mass1_mass2_to_mchirp_eta, f_SchwarzISCO from pycbc.pnutils import nearest_larger_binary_number from pycbc.pnutils import mass1_mass2_to_tau0_tau3 @@ -35,11 +39,9 @@ from pycbc.filter import overlap_cplx, matched_filter from pycbc.filter import compute_max_snr_over_sky_loc_stat from pycbc.filter import compute_max_snr_over_sky_loc_stat_no_phase from pycbc.io.ligolw import LIGOLWContentHandler -from math import ceil, log import pycbc.psd, pycbc.scheme, pycbc.fft, pycbc.strain, pycbc.version from pycbc.detector import overhead_antenna_pattern as generate_fplus_fcross from pycbc.waveform import TemplateBank -from tqdm import tqdm ## Remove the need for these functions ######################################## @@ -150,8 +152,7 @@ aprs = sorted(list(set(td_approximants() + fd_approximants()))) parser = ArgumentParser(description=__doc__) parser.add_argument("--match-file", dest="out_file", metavar="FILE", required=True, help="File to output match results") -parser.add_argument("--verbose", action='store_true', default=False, - help="Print verbose statements") +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) diff --git a/bin/pycbc_coinc_time b/bin/pycbc_coinc_time index 7f0803fa5c4..34e4687ea63 100644 --- a/bin/pycbc_coinc_time +++ b/bin/pycbc_coinc_time @@ -1,7 +1,13 @@ #!/bin/env python -import argparse, logging, pycbc.version, ligo.segments, numpy +import argparse +import logging +import numpy from dqsegdb.apicalls import dqsegdbQueryTimes as query + +import ligo.segments + #from pycbc.workflow.segment import cat_to_veto_def_cat as convert_cat +import pycbc.version def sane(seg_list): """ Convert list of len two lists containing strs to segment list """ @@ -85,7 +91,7 @@ def get_vetoes(veto_def, ifo, server, veto_name, start_default, end_default): parser = argparse.ArgumentParser() parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') +pycbc.add_common_pycbc_options(parser) parser.add_argument('--gps-start-time', type=int, required=True, help="integer gps start time") diff --git a/bin/pycbc_compress_bank b/bin/pycbc_compress_bank index 851e6bd78b9..0d98bbffc75 100755 --- a/bin/pycbc_compress_bank +++ b/bin/pycbc_compress_bank @@ -30,6 +30,7 @@ import argparse import numpy import h5py import logging + import pycbc from pycbc import psd, DYN_RANGE_FAC from pycbc.waveform import compress @@ -40,6 +41,7 @@ from pycbc import filter parser = argparse.ArgumentParser(description=__description__) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--bank-file", type=str, required=True, help="Bank hdf file to load.") parser.add_argument("--output", type=str, required=True, @@ -96,7 +98,6 @@ parser.add_argument("--precision", type=str, choices=["double", "single"], parser.add_argument("--force", action="store_true", default=False, help="Overwrite the given hdf file if it exists. " "Otherwise, an error is raised.") -parser.add_argument("--verbose", action="store_true", default=False) # Insert the PSD options pycbc.psd.insert_psd_option_group(parser, include_data_options=False) diff --git a/bin/pycbc_convertinjfiletohdf b/bin/pycbc_convertinjfiletohdf index 9c703d51f98..0998024f30a 100755 --- a/bin/pycbc_convertinjfiletohdf +++ b/bin/pycbc_convertinjfiletohdf @@ -23,6 +23,7 @@ or the LIGO HDF format) into a PyCBC hdf injection format import argparse import numpy import h5py +import shutil import pycbc from pycbc.inject import InjectionSet, CBCHDFInjectionSet @@ -150,11 +151,16 @@ class LVKNewStyleInjectionSet(object): elif lvk_name == 't_co_gps_add': continue else: - data[pycbc_name] = self.inj_file[f'{self.subdir}/{lvk_name}'][:] + lvk_file_dset = self.inj_file[f'{self.subdir}/{lvk_name}'] + if lvk_file_dset.dtype.char in ['U', 'O']: + data[pycbc_name] = lvk_file_dset[:].astype('S') + else: + data[pycbc_name] = lvk_file_dset[:] return data parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg) parser.add_argument('--injection-file', required=True, @@ -162,8 +168,6 @@ parser.add_argument('--injection-file', required=True, "and must contain a SimInspiral table") parser.add_argument('--output-file', required=True, help="The ouput file name. Must end in '.hdf'.") -parser.add_argument("--verbose", action="store_true", default=False, - help="Be verbose.") args = parser.parse_args() pycbc.init_logging(args.verbose) @@ -179,14 +183,11 @@ else: injclass = LVKNewStyleInjectionSet(args.injection_file) inj_file.close() -if injclass is None: - # Assume this is a PyCBC format - injclass = InjectionSet(args.injection_file) - data = {} - for key in xinj.table[0].__slots__: - data[str(key)] = numpy.array([getattr(t, key) for t in xinj.table]) -else: +if injclass is not None: + # The injection file is of known type data = injclass.pack_data_into_pycbc_format_input() - -samples = FieldArray.from_kwargs(**data) -CBCHDFInjectionSet.write(args.output_file, samples) + samples = FieldArray.from_kwargs(**data) + CBCHDFInjectionSet.write(args.output_file, samples) +else: + # Assume the injection is a HDF file (PyCBC format), only make a copy + shutil.copy(args.injection_file, args.output_file) diff --git a/bin/pycbc_create_injections b/bin/pycbc_create_injections index a9616808a94..3fac44c3d66 100644 --- a/bin/pycbc_create_injections +++ b/bin/pycbc_create_injections @@ -108,10 +108,14 @@ along with any numpy ufunc. So, in the above example, mass ratio (q) is constrained to be <= 4 by using a function from the conversions module. """ -import os, sys +import os +import sys import argparse import logging import numpy +import h5py +from numpy.random import uniform + import pycbc import pycbc.version from pycbc.inject import InjectionSet @@ -122,12 +126,11 @@ from pycbc.io import record from pycbc.waveform import parameters from pycbc.workflow import configuration from pycbc.workflow import WorkflowConfigParser -from numpy.random import uniform -import h5py parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) configuration.add_workflow_command_line_group(parser) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg, help='Prints version information.') @@ -163,8 +166,6 @@ parser.add_argument('--static-params-section', default="static_params", parser.add_argument("--force", action="store_true", default=False, help="If the output-file already exists, overwrite it. " "Otherwise, an OSError is raised.") -parser.add_argument("--verbose", action="store_true", default=False, - help="Print logging messages.") opts = parser.parse_args() pycbc.init_logging(opts.verbose) diff --git a/bin/pycbc_data_store b/bin/pycbc_data_store index eb501e709a3..a6b17bac67f 100755 --- a/bin/pycbc_data_store +++ b/bin/pycbc_data_store @@ -1,16 +1,23 @@ #!/usr/bin/env python """ Create HDF strain cache file """ -import logging, argparse, numpy, h5py -import pycbc, pycbc.strain, pycbc.dq +import logging +import argparse +import numpy +import h5py + +import pycbc +import pycbc.strain +import pycbc.dq from pycbc.version import git_verbose_msg as version from pycbc.fft.fftw import set_measure_level from pycbc.events.veto import segments_to_start_end + 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("--science-name", help="Science flag definition") parser.add_argument("--segment-server") parser.add_argument("--veto-definer-file") diff --git a/bin/pycbc_faithsim b/bin/pycbc_faithsim index c1405316909..e1bf979a15a 100644 --- a/bin/pycbc_faithsim +++ b/bin/pycbc_faithsim @@ -27,6 +27,7 @@ import logging from numpy import complex64 import argparse import sys + from ligo.lw import utils as ligolw_utils from ligo.lw import lsctables @@ -84,8 +85,7 @@ parser = argparse.ArgumentParser(usage='', description="Calculate faithfulness for a set of waveforms.") 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) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--param-file", dest="bank_file", metavar="FILE", help="Sngl or Sim Inspiral Table containing waveform " "parameters.") diff --git a/bin/pycbc_faithsim_collect_results b/bin/pycbc_faithsim_collect_results index 7d98ac02b22..a574e5e0415 100755 --- a/bin/pycbc_faithsim_collect_results +++ b/bin/pycbc_faithsim_collect_results @@ -5,12 +5,16 @@ Program for collecting the results of pycbc_faithsim comparing two approximants computing the match between them and creating a .dat file with the results. """ +import argparse import numpy as np + from ligo.lw import utils, table, lsctables + +from pycbc import add_common_pycbc_options, init_logging from pycbc.io.ligolw import LIGOLWContentHandler -import argparse parser = argparse.ArgumentParser(description=__doc__) +add_common_pycbc_options(parser) parser.add_argument( "--match-inputs", action="append", @@ -26,6 +30,8 @@ parser.add_argument( parser.add_argument("--output", required=True, help="name of the output .dat file") args = parser.parse_args() +init_logging(options.verbose) + mfields = ("match", "overlap", "time_offset", "sigma1", "sigma2") bfields = ( "match", diff --git a/bin/pycbc_fit_sngl_trigs b/bin/pycbc_fit_sngl_trigs index cef76488606..9f0c347656d 100644 --- a/bin/pycbc_fit_sngl_trigs +++ b/bin/pycbc_fit_sngl_trigs @@ -55,10 +55,8 @@ def get_bins(opt, pmin, pmax): 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("--inputs", nargs="+", help="Input file or space-separated list of input files " "containing single triggers. Currently .xml(.gz) " @@ -134,11 +132,7 @@ binopt.add_argument("--irregular-bins", opt = parser.parse_args() -if opt.verbose: - log_level = logging.DEBUG -else: - log_level = logging.WARN -logging.basicConfig(format='%(asctime)s : %(message)s', level=log_level) +pycbc.init_logging(opt.verbose) statname = opt.sngl_stat.replace("_", " ") paramname = opt.fit_param.replace("_", " ") diff --git a/bin/pycbc_gwosc_segment_query b/bin/pycbc_gwosc_segment_query index 4403af29168..0ac34512a33 100644 --- a/bin/pycbc_gwosc_segment_query +++ b/bin/pycbc_gwosc_segment_query @@ -6,14 +6,11 @@ import json import argparse import shutil from urllib.request import urlopen -import ligo.segments -from pycbc.workflow import SegFile +import ligo.segments -# Logging formatting from pycbc optimal snr -log_fmt = '%(asctime)s %(message)s' -log_date_fmt = '%Y-%m-%d %H:%M:%S' -logging.basicConfig(level=logging.INFO, format=log_fmt, datefmt=log_date_fmt) +import pycbc +from pycbc.workflow import SegFile # Function to query json segment data from GWOSC @@ -64,6 +61,7 @@ def write_xml_file(ifo, summary_segment, segments, filename): parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument('--gps-start-time', type=int, required=True) parser.add_argument('--gps-end-time', type=int, required=True) parser.add_argument('--include-segments', type=str, required=True) @@ -71,6 +69,8 @@ parser.add_argument('--output-file', type=str, required=True) parser.add_argument('--protract-hw-inj', type=int, default=0) args = parser.parse_args() +pycbc.init_logging(args.verbose) + gps_start_time = args.gps_start_time gps_end_time = args.gps_end_time duration = gps_end_time - gps_start_time diff --git a/bin/pycbc_hdf5_splitbank b/bin/pycbc_hdf5_splitbank index e86e649c28b..d080d5774bb 100755 --- a/bin/pycbc_hdf5_splitbank +++ b/bin/pycbc_hdf5_splitbank @@ -25,13 +25,15 @@ import argparse import numpy import h5py import logging +from numpy import random + import pycbc, pycbc.version from pycbc.waveform import bank -from numpy import random __author__ = "Soumi De " parser = argparse.ArgumentParser(description=__doc__[1:]) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) parser.add_argument("--bank-file", type=str, @@ -61,11 +63,11 @@ parser.add_argument("--random-seed", type=int, parser.add_argument("--force", action="store_true", default=False, help="Overwrite the given hdf file if it exists. " "Otherwise, an error is raised.") -parser.add_argument("--verbose", action="store_true", default=False) - args = parser.parse_args() +pycbc.init_logging(args.verbose) + # input checks if args.mchirp_sort and (args.random_seed is not None): raise RuntimeError("Can't supply a random seed if not sorting randomly!") diff --git a/bin/pycbc_hdf_splitinj b/bin/pycbc_hdf_splitinj index 0d255e3850a..6a6050f5abd 100644 --- a/bin/pycbc_hdf_splitinj +++ b/bin/pycbc_hdf_splitinj @@ -6,14 +6,16 @@ Split sets are organized to maximize time between injections. """ import argparse -from pycbc.inject import InjectionSet import h5py import numpy as np + +from pycbc.inject import InjectionSet import pycbc.version # Parse command line 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("-f", "--output-files", nargs='*', required=True, @@ -23,6 +25,8 @@ parser.add_argument("-i", "--input-file", required=True, args = parser.parse_args() +pycbc.init_logging(args.verbose) + # Read in input file as both an hdf file and an InjectionSet object inj_file = h5py.File(args.input_file, 'r') inj_set = InjectionSet(args.input_file) diff --git a/bin/pycbc_inj_cut b/bin/pycbc_inj_cut index 0622d21f0d2..1ca70bcb5ca 100644 --- a/bin/pycbc_inj_cut +++ b/bin/pycbc_inj_cut @@ -29,14 +29,17 @@ __email__ = "stephen.fairhurst@ligo.org" import argparse import logging import numpy + +from ligo.lw import utils as ligolw_utils +from ligo.lw import lsctables + import pycbc import pycbc.inject -from ligo.lw import utils as ligolw_utils from pycbc.types import MultiDetOptionAction -from ligo.lw import lsctables import pycbc.version parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action=pycbc.version.Version) parser.add_argument('--input', dest='inj_xml', required=True, help='Input LIGOLW injections file.') parser.add_argument('--output-missed', dest='output_missed', required=False, @@ -53,20 +56,13 @@ parser.add_argument('--snr-columns', nargs='+', action=MultiDetOptionAction, ' no useful data in it, good candidates are usually' \ ' alpha1, alpha2 etc.') parser.add_argument("-z", "--write-compress", action="store_true", help="Write compressed xml.gz files.") -parser.add_argument("-v", "--verbose", action="store_true", default=False, - help="Extended standard output.") opts = parser.parse_args() -log_fmt = '%(asctime)s %(message)s' -log_date_fmt = '%Y-%m-%d %H:%M:%S' -logging.basicConfig(level=logging.INFO, format=log_fmt, datefmt=log_date_fmt) +pycbc.init_logging(opts.verbose) -if opts.verbose: - logging.info("Loading injections") injections = pycbc.inject.InjectionSet(opts.inj_xml) - # injections that we should do: table and name of the file that will store it output_inject = opts.output_inject if opts.write_compress: @@ -89,11 +85,20 @@ for i, inj in enumerate(injections.table): else: out_sim_inspiral_inject.append(inj) -logging.info('Found %d/%d (potentially) found injections. Storing them to %s.', len(out_sim_inspiral_inject), len(injections.table), output_inject) -logging.info('Found %d/%d missed injections.', len(out_sim_inspiral_missed), len(injections.table)) - -if opts.verbose: - logging.info('Writing output') +logging.info( + 'Found %d/%d (potentially) found injections. Storing them to %s.', + len(out_sim_inspiral_inject), + len(injections.table), + output_inject +) + +logging.info( + 'Found %d/%d missed injections.', + len(out_sim_inspiral_missed), + len(injections.table) +) + +logging.info('Writing output') llw_doc = injections.indoc llw_root = llw_doc.childNodes[0] llw_root.removeChild(injections.table) diff --git a/bin/pycbc_inspiral b/bin/pycbc_inspiral index 4d5e6e9a585..3d4924c3e0b 100644 --- a/bin/pycbc_inspiral +++ b/bin/pycbc_inspiral @@ -16,19 +16,24 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -import sys, os, copy -import logging, argparse, numpy, itertools +import sys +import os +import copy +import logging +import argparse +import numpy +import itertools +import time +from multiprocessing import Pool + import pycbc import pycbc.version from pycbc import vetoes, psd, waveform, strain, scheme, fft, DYN_RANGE_FAC, events from pycbc.vetoes.sgchisq import SingleDetSGChisq from pycbc.filter import MatchedFilterControl, make_frequency_series, qtransform from pycbc.types import TimeSeries, FrequencySeries, zeros, float32, complex64 -from multiprocessing import Pool -import pycbc.version import pycbc.opt import pycbc.inject -import time last_progress_update = -1.0 @@ -48,9 +53,8 @@ tstart = time.time() parser = argparse.ArgumentParser(usage='', description="Find single detector gravitational-wave triggers.") +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("--update-progress", help="updates a file 'progress.txt' with a value 0 .. 1.0 when this amount of (filtering) progress was made", type=float, default=0) diff --git a/bin/pycbc_inspiral_skymax b/bin/pycbc_inspiral_skymax index f082f7d3841..2f2c477f72d 100755 --- a/bin/pycbc_inspiral_skymax +++ b/bin/pycbc_inspiral_skymax @@ -16,7 +16,11 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -import logging, argparse, sys, numpy +import logging +import argparse +import sys +import numpy + from pycbc import vetoes, psd, waveform, events, strain, scheme, fft from pycbc.vetoes.sgchisq import SingleDetSGChisq from pycbc import DYN_RANGE_FAC @@ -28,8 +32,7 @@ import pycbc.opt parser = argparse.ArgumentParser(usage='', description="Find single detector precessing gravitational-wave triggers.") -parser.add_argument("-V", "--verbose", action="store_true", - help="print extra debugging information", default=False ) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--output", type=str, help="FIXME: ADD") parser.add_argument("--bank-file", type=str, help="FIXME: ADD") parser.add_argument("--snr-threshold", diff --git a/bin/pycbc_live b/bin/pycbc_live index 646d4bebf40..7cfb3af8290 100755 --- a/bin/pycbc_live +++ b/bin/pycbc_live @@ -82,7 +82,13 @@ class LiveEventManager(object): def __init__(self, args, bank): self.low_frequency_cutoff = args.low_frequency_cutoff self.bank = bank - self.skymap_only_ifos = [] if args.skymap_only_ifos is None else list(set(args.skymap_only_ifos)) + # all interferometers involved in the analysis, whether for generating + # candidates or doing localization only + self.ifos = set(args.channel_name.keys()) + # interferometers used for localization only + self.skymap_only_ifos = set(args.skymap_only_ifos or []) + # subset of the interferometers allowed to produce candidates + self.trigg_ifos = self.ifos - self.skymap_only_ifos # Figure out what we are supposed to process within the pool of MPI processes self.comm = mpi.COMM_WORLD @@ -106,6 +112,7 @@ class LiveEventManager(object): self.enable_gracedb_upload = args.enable_gracedb_upload self.run_snr_optimization = args.run_snr_optimization self.snr_opt_label = args.snr_opt_label + self.snr_opt_options = args.snr_opt_extra_opts self.gracedb = None # Keep track of which events have been uploaded @@ -116,12 +123,11 @@ class LiveEventManager(object): if self.run_snr_optimization: # preestimate the number of CPU cores that we can afford giving # to followup processes without slowing down the main search - bg_cores = len(tuple(itertools.combinations(ifos, 2))) + bg_cores = len(tuple(itertools.combinations(self.trigg_ifos, 2))) analysis_cores = 1 + bg_cores if platform.system() != 'Darwin': available_cores = len(os.sched_getaffinity(0)) self.fu_cores = available_cores - analysis_cores - self.optimizer = args.snr_opt_method if self.fu_cores <= 0: logging.warning( 'Insufficient number of CPU cores (%d) to ' @@ -133,10 +139,6 @@ class LiveEventManager(object): else: # To enable mac testing, this is just set to 1 self.fu_cores = 1 - # Convert SNR optimizer options into a string - self.snr_opt_options = snr_optimizer.args_to_string(args) - else: - self.snr_opt_options = None if args.enable_embright_has_massgap: if args.embright_massgap_max < self.mc_area_args['mass_bdary']['ns_max']: @@ -506,7 +508,7 @@ class LiveEventManager(object): logging.info('computing followup data for coinc') coinc_ifos = coinc_results['foreground/type'].split('-') followup_ifos = set(ifos) - set(coinc_ifos) - followup_ifos = list(followup_ifos | set(self.skymap_only_ifos)) + followup_ifos = list(followup_ifos | self.skymap_only_ifos) double_ifar = coinc_results['foreground/ifar'] if double_ifar < args.ifar_double_followup_threshold: @@ -540,14 +542,24 @@ class LiveEventManager(object): ) logging.info('Coincident candidate! Saving as %s', fname) + # Which IFOs were active? + live_ifos = [ifo for ifo in sld if 'snr_series' in sld[ifo]] + # Verbally explain some details not obvious from the other info - comment = ('Trigger produced as a {} coincidence. ' - 'FAR is based on all listed detectors.
' - 'Two-detector ranking statistic: {}
' - 'Followup detectors: {}') - comment = comment.format(ppdets(coinc_ifos), - args.ranking_statistic, - ppdets(followup_ifos)) + comment = ( + 'Trigger produced as a {} coincidence.
' + 'Two-detector ranking statistic: {}
' + 'Detectors used for FAR calculation: {}.
' + 'Detectors used for localization: {}.
' + 'Detectors used only for localization: {}.' + ) + comment = comment.format( + ppdets(coinc_ifos), + args.ranking_statistic, + ppdets(set(ifos) - self.skymap_only_ifos), + ppdets(live_ifos), + ppdets(set(live_ifos) & self.skymap_only_ifos) + ) ifar = coinc_results['foreground/ifar'] upload_checks = self.enable_gracedb_upload and self.ifar_upload_threshold < ifar @@ -568,9 +580,6 @@ class LiveEventManager(object): # even if not running it - do this before the thread so no # data buffers move on in a possible interim period - # Which IFOs were active? - live_ifos = [ifo for ifo in sld if 'snr_series' in sld[ifo]] - # Tell SNR optimized event about p_terr if hasattr(event, 'p_terr') and event.p_terr is not None: coinc_results['p_terr'] = event.p_terr @@ -608,7 +617,7 @@ class LiveEventManager(object): logging.info(f'Found {ifo} single with ifar {sifar}') followup_ifos = [i for i in active if i is not ifo] - followup_ifos = list(set(followup_ifos) | set(self.skymap_only_ifos)) + followup_ifos = list(set(followup_ifos) | self.skymap_only_ifos) # Don't recompute ifar considering other ifos sld = self.compute_followup_data( [ifo], @@ -639,11 +648,21 @@ class LiveEventManager(object): ) logging.info('Single-detector candidate! Saving as %s', fname) + # Which IFOs were active? + live_ifos = [ifo for ifo in sld if 'snr_series' in sld[ifo]] + # Verbally explain some details not obvious from the other info - comment = ('Trigger produced as a {0} single. ' - 'FAR is based on {0} only.
' - 'Followup detectors: {1}') - comment = comment.format(ifo, ppdets(followup_ifos)) + comment = ( + 'Trigger produced as a {0} single.
' + 'Detectors used for FAR calculation: {0}.
' + 'Detectors used for localization: {1}.
' + 'Detectors used only for localization: {2}.' + ) + comment = comment.format( + ifo, + ppdets(live_ifos), + ppdets(set(live_ifos) & self.skymap_only_ifos) + ) # Has a coinc event at this time been uploaded recently? # If so, skip upload - Note that this means that we _always_ @@ -680,9 +699,6 @@ class LiveEventManager(object): # where there is already a coinc continue - # Which IFOs were active? - live_ifos = [ifo for ifo in sld if 'snr_series' in sld[ifo]] - # Tell SNR optimized event about p_terr if hasattr(event, 'p_terr') and event.p_terr is not None: single['p_terr'] = event.p_terr @@ -987,6 +1003,11 @@ parser.add_argument('--snr-opt-timeout', type=int, default=400, metavar='SECONDS help='Maximum allowed duration of followup process to maximize SNR') parser.add_argument('--snr-opt-label', default='SNR_OPTIMIZED', help='Label to apply to snr-optimized GraceDB uploads') +parser.add_argument('--snr-opt-extra-opts', + help='Extra options to pass to the optimizer subprocess. Example: ' + '--snr-opt-extra-opts "--snr-opt-method differential_evolution ' + '--snr-opt-di-maxiter 50 --snr-opt-di-popsize 100 ' + '--snr-opt-seed 42 --snr-opt-include-candidate "') parser.add_argument('--enable-embright-has-massgap', action='store_true', default=False, help='Estimate HasMassGap probability for EMBright info. Lower limit ' @@ -1010,15 +1031,12 @@ Coincer.insert_args(parser) SingleDetSGChisq.insert_option_group(parser) mchirp_area.insert_args(parser) livepau.insert_live_pastro_option_group(parser) -snr_optimizer.insert_snr_optimizer_options(parser) args = parser.parse_args() scheme.verify_processing_options(args, parser) fft.verify_fft_options(args, parser) Coincer.verify_args(args, parser) -ifos = set(args.channel_name.keys()) -analyze_singles = LiveSingle.verify_args(args, parser, ifos) if args.output_background is not None and len(args.output_background) != 2: parser.error('--output-background takes two parameters: period and path') @@ -1045,14 +1063,16 @@ if bank.min_f_lower < args.low_frequency_cutoff: 'minimum f_lower across all templates ' '({} Hz)'.format(args.low_frequency_cutoff, bank.min_f_lower)) -logging.info('Analyzing data from detectors %s', ppdets(ifos)) - evnt = LiveEventManager(args, bank) -logging.info('Detectors that only aid in the sky localization %s', ppdets(evnt.skymap_only_ifos)) + +logging.info('Analyzing data from detectors %s', ppdets(evnt.ifos)) +logging.info('Using %s for localization only', ppdets(evnt.skymap_only_ifos)) + +analyze_singles = LiveSingle.verify_args(args, parser, evnt.trigg_ifos) # include MPI rank and functional description into proctitle task_name = 'root' if evnt.rank == 0 else 'filtering' -setproctitle('PyCBC Live rank {:d} [{}]'.format(evnt.rank, task_name)) +setproctitle(f'PyCBC Live rank {evnt.rank:d} [{task_name}]') sg_chisq = SingleDetSGChisq.from_cli(args, bank, args.chisq_bins) @@ -1111,23 +1131,27 @@ with ctx: args.round_start_time logging.info('Starting from: %s', args.start_time) - # initialize the data readers for all detectors + # Initialize the data readers for all detectors. For rank 0, we need data + # from all detectors, including the localization-only ones. For higher + # ranks, we only need the detectors that can generate candidates. if args.max_length is not None: maxlen = args.max_length maxlen = int(maxlen) - data_reader = {ifo: StrainBuffer.from_cli(ifo, args, maxlen) - for ifo in ifos} + data_reader = { + ifo: StrainBuffer.from_cli(ifo, args, maxlen) + for ifo in (evnt.ifos if evnt.rank == 0 else evnt.trigg_ifos) + } evnt.data_readers = data_reader # create single-detector background "estimators" if analyze_singles and evnt.rank == 0: sngl_estimator = {ifo: LiveSingle.from_cli(args, ifo) - for ifo in ifos} + for ifo in evnt.trigg_ifos} - # Create double coincident background estimator for every combo + # Create double coincident background estimator + # for every pair of triggering interferometers if args.enable_background_estimation and evnt.rank == 0: - trigg_ifos = [ifo for ifo in ifos if ifo not in evnt.skymap_only_ifos] - ifo_combos = itertools.combinations(trigg_ifos, 2) + ifo_combos = itertools.combinations(evnt.trigg_ifos, 2) estimators = [] for combo in ifo_combos: logging.info('Will calculate %s background', ppdets(combo, "-")) @@ -1167,12 +1191,14 @@ with ctx: # main analysis loop data_end = lambda: data_reader[tuple(data_reader.keys())[0]].end_time last_bg_dump_time = int(data_end()) - psd_count = {ifo:0 for ifo in ifos} + psd_count = {ifo:0 for ifo in evnt.ifos} # Create dicts to track whether the psd has been recalculated and to hold # psd variation filters - psd_recalculated = {ifo: True for ifo in ifos} - psd_var_filts = {ifo: None for ifo in ifos} + psd_recalculated = { + ifo: True for ifo in (evnt.ifos if evnt.rank == 0 else evnt.trigg_ifos) + } + psd_var_filts = {ifo: None for ifo in evnt.trigg_ifos} while data_end() < args.end_time: t1 = pycbc.gps_now() @@ -1181,7 +1207,7 @@ with ctx: results = {} evnt.live_detectors = set() - for ifo in ifos: + for ifo in (evnt.ifos if evnt.rank == 0 else evnt.trigg_ifos): results[ifo] = False status = data_reader[ifo].advance( valid_pad, @@ -1259,23 +1285,27 @@ with ctx: # Calculate and add the psd variation for the results if args.psd_variation: - for ifo in results: - logging.info(f"Calculating PSD Variation Statistic for {ifo}") + logging.info("Calculating PSD Variation Statistic for %s", ifo) # A new filter is needed if the PSD has been recalculated if psd_recalculated[ifo] is True: - psd_var_filts[ifo] = variation.live_create_filter(data_reader[ifo].psd, - args.psd_segment_length, - int(args.sample_rate)) + psd_var_filts[ifo] = variation.live_create_filter( + data_reader[ifo].psd, + args.psd_segment_length, + int(args.sample_rate) + ) psd_recalculated[ifo] = False - psd_var_ts = variation.live_calc_psd_variation(data_reader[ifo].strain, - psd_var_filts[ifo], - args.increment) + psd_var_ts = variation.live_calc_psd_variation( + data_reader[ifo].strain, + psd_var_filts[ifo], + args.increment + ) - psd_var_vals = variation.live_find_var_value(results[ifo], - psd_var_ts) + psd_var_vals = variation.live_find_var_value( + results[ifo], psd_var_ts + ) results[ifo]['psd_var_val'] = psd_var_vals @@ -1295,7 +1325,7 @@ with ctx: gates = {ifo: data_reader[ifo].gate_params for ifo in data_reader} # map the results file to an hdf file - prefix = '{}-{}-{}-{}'.format(''.join(sorted(ifos)), + prefix = '{}-{}-{}-{}'.format(''.join(sorted(evnt.ifos)), args.file_prefix, data_end() - args.analysis_chunk, valid_pad) @@ -1310,8 +1340,9 @@ with ctx: data_end() - last_bg_dump_time > float(args.output_background[0]): last_bg_dump_time = int(data_end()) bg_dists = coinc_pool.broadcast(output_background, None) - bg_fn = '{}-LIVE_BACKGROUND-{}.hdf'.format(''.join(sorted(ifos)), - last_bg_dump_time) + bg_fn = '{}-LIVE_BACKGROUND-{}.hdf'.format( + ''.join(sorted(evnt.trigg_ifos)), last_bg_dump_time + ) bg_fn = os.path.join(args.output_background[1], bg_fn) with h5py.File(bg_fn, 'w') as bgf: for bg_ifos, bg_data, bg_time in bg_dists: diff --git a/bin/pycbc_make_banksim b/bin/pycbc_make_banksim index cfbc8985bca..09b0efbad5c 100644 --- a/bin/pycbc_make_banksim +++ b/bin/pycbc_make_banksim @@ -9,6 +9,7 @@ import tempfile from optparse import OptionParser from glue.pipeline import CondorDAGJob, CondorDAGNode, CondorDAG, CondorJob +from pycbc import init_logging class BaseJob(CondorDAGJob, CondorJob): def __init__(self, log_dir, executable, cp, section, gpu=False, @@ -131,7 +132,8 @@ parser.add_option('--config', type=str) if options.config is None: raise ValueError("Config file is required") -logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO) +# logging INFO +init_logging(1) confs = ConfigParser.ConfigParser() confs.read(options.config) diff --git a/bin/pycbc_make_faithsim b/bin/pycbc_make_faithsim index 8c1f33bbc3e..c27934745c3 100644 --- a/bin/pycbc_make_faithsim +++ b/bin/pycbc_make_faithsim @@ -9,6 +9,7 @@ import tempfile from optparse import OptionParser from glue.pipeline import CondorDAGJob, CondorDAGNode, CondorDAG, CondorJob +from pycbc import init_logging class BaseJob(CondorDAGJob, CondorJob): def __init__(self, log_dir, executable, cp, section, accounting_group=None): @@ -65,7 +66,8 @@ parser.add_option('--config', type=str) if options.config is None: raise ValueError("Config file is required") -logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO) +# logging INFO +init_logging(1) confs = ConfigParser.ConfigParser() confs.read(options.config) diff --git a/bin/pycbc_make_html_page b/bin/pycbc_make_html_page index 1ed5612dd48..f4ea0f2e809 100644 --- a/bin/pycbc_make_html_page +++ b/bin/pycbc_make_html_page @@ -17,14 +17,17 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. import argparse -import os, stat -import pycbc.results +import os +import stat import random import shutil import zipfile import codecs -from ligo import segments from jinja2 import Environment, FileSystemLoader + +from ligo import segments + +import pycbc.results from pycbc.results.render import get_embedded_config, render_workflow_html_template, setup_template_render from pycbc.workflow import segment import pycbc.version @@ -166,6 +169,7 @@ parser = argparse.ArgumentParser(usage='pycbc_make_html_page \ description="Create static html pages of a filesystem's content.") parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) +pycbc.add_common_pycbc_options(parser) parser.add_argument('-f', '--template-file', type=str, help='Template file to use for skeleton html page.') parser.add_argument('-b', '--output-path', type=str, @@ -181,12 +185,12 @@ parser.add_argument('-s', '--analysis-subtitle', type=str, parser.add_argument('-l', '--logo', type=str, default=default_logo_location, help='File location of logo to include at top of page') -parser.add_argument('-v', '--verbose', action='store_true', - help='Print extra debugging information.', default=False) parser.add_argument('-P', '--protect-results', action='store_true', help='Make the output web pages read only.', default=False) opts = parser.parse_args() +pycbc.init_logging(opts.verbose) + # edit command line options analysis_title = opts.analysis_title.strip('"').rstrip('"') analysis_subtitle = opts.analysis_subtitle.strip('"').rstrip('"') diff --git a/bin/pycbc_make_sky_grid b/bin/pycbc_make_sky_grid index 8a1809441a4..e6e2bf85ddd 100644 --- a/bin/pycbc_make_sky_grid +++ b/bin/pycbc_make_sky_grid @@ -10,10 +10,13 @@ The grid is constructed following the method described in Section V of https://a import numpy as np import h5py import argparse -from pycbc.detector import Detector import itertools from scipy.spatial.transform import Rotation as R +import pycbc +from pycbc.detector import Detector + + def spher_to_cart(sky_points): """Convert spherical coordinates to cartesian coordinates. """ @@ -32,16 +35,29 @@ def cart_to_spher(sky_points): return spher parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument('--ra', type=float, help="Right ascension (in rad) of the center of the external trigger error box") -parser.add_argument('--dec', type=float, help="Declination (in rad) of the center of the external trigger error box") -parser.add_argument('--instruments', nargs="+", type=str, required=True, help="List of instruments to analyze.") -parser.add_argument('--sky-error', type=float, required=True, help="3-sigma confidence radius (in rad) of the external trigger error box") -parser.add_argument('--trigger-time', type=int, required=True, help="Time (in s) of the external trigger") -parser.add_argument('--timing-uncertainty', type=float, default=0.0001, help="Timing uncertainty (in s) we are willing to accept") -parser.add_argument('--output', type=str, required=True, help="Name of the sky grid") +pycbc.add_common_pycbc_options(parser) +parser.add_argument('--ra', type=float, + help="Right ascension (in rad) of the center of the external trigger " + "error box") +parser.add_argument('--dec', type=float, + help="Declination (in rad) of the center of the external trigger " + "error box") +parser.add_argument('--instruments', nargs="+", type=str, required=True, + help="List of instruments to analyze.") +parser.add_argument('--sky-error', type=float, required=True, + help="3-sigma confidence radius (in rad) of the external trigger error " + "box") +parser.add_argument('--trigger-time', type=int, required=True, + help="Time (in s) of the external trigger") +parser.add_argument('--timing-uncertainty', type=float, default=0.0001, + help="Timing uncertainty (in s) we are willing to accept") +parser.add_argument('--output', type=str, required=True, + help="Name of the sky grid") args = parser.parse_args() +pycbc.init_logging(args.verbose) + if len(args.instruments) == 1: parser.error('Can not make a sky grid for only one detector.') diff --git a/bin/pycbc_make_skymap b/bin/pycbc_make_skymap index 4b011831011..89e5933db88 100755 --- a/bin/pycbc_make_skymap +++ b/bin/pycbc_make_skymap @@ -19,6 +19,9 @@ mpl_use_backend('agg') import numpy as np import h5py + +from ligo.gracedb.rest import GraceDb + import pycbc from pycbc.filter import sigmasq from pycbc.io import live, WaveformArray @@ -29,7 +32,6 @@ from pycbc.pnutils import nearest_larger_binary_number from pycbc.waveform.spa_tmplt import spa_length_in_time from pycbc import frame, waveform, dq from pycbc.psd import interpolate -from ligo.gracedb.rest import GraceDb def default_frame_type(time, ifo): @@ -476,9 +478,9 @@ def main(trig_time, mass1, mass2, spin1z, spin2z, f_low, f_upper, sample_rate, shutil.rmtree(tmpdir) if __name__ == '__main__': - pycbc.init_logging(True) parser = argparse.ArgumentParser(description=__doc__) + pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action=pycbc.version.Version) # note that I am not using a MultiDetOptionAction for --trig-time as I # explicitly want to handle cases like `--trig-time 1234` and @@ -563,6 +565,8 @@ if __name__ == '__main__': opt = parser.parse_args() + pycbc.init_logging(opt.verbose) + frame_type_dict = {f.split(':')[0]: f.split(':')[1] for f in opt.frame_type} \ if opt.frame_type is not None else None chan_name_dict = {f.split(':')[0]: f for f in opt.channel_name} \ diff --git a/bin/pycbc_merge_inj_hdf b/bin/pycbc_merge_inj_hdf index fe0d698e935..9b208738d47 100755 --- a/bin/pycbc_merge_inj_hdf +++ b/bin/pycbc_merge_inj_hdf @@ -24,6 +24,7 @@ import logging import argparse import numpy as np import h5py + import pycbc import pycbc.inject import pycbc.version @@ -43,20 +44,16 @@ def get_gc_end_time(injection): if __name__ == '__main__': parser = argparse.ArgumentParser(description=__doc__) + pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action=pycbc.version.Version) parser.add_argument('--injection-files', '-i', dest='injection_file', required=True, nargs='+', help='Input HDF5 files defining injections') parser.add_argument('--output-file', '-o', dest='out_file', required=True, help='Output HDF5 file') - parser.add_argument("--verbose", action="store_true", default=False, - help="Print logging messages") opts = parser.parse_args() - log_level = logging.INFO if opts.verbose else logging.WARN - logging.basicConfig(level=log_level, - format='%(asctime)s %(message)s', - datefmt='%Y-%m-%d %H:%M:%S') + pycbc.init_logging(opts.verbose) logging.info("Loading injections") diff --git a/bin/pycbc_multi_inspiral b/bin/pycbc_multi_inspiral index c0e8189c57e..3d9967f42d0 100755 --- a/bin/pycbc_multi_inspiral +++ b/bin/pycbc_multi_inspiral @@ -31,9 +31,10 @@ from collections import defaultdict import argparse import numpy as np import h5py + from pycbc import ( detector, fft, init_logging, inject, opt, psd, scheme, strain, vetoes, - waveform, DYN_RANGE_FAC + waveform, DYN_RANGE_FAC, add_common_pycbc_options ) from pycbc.events import ranking, coherent as coh, EventManagerCoherent from pycbc.filter import MatchedFilterControl @@ -46,9 +47,7 @@ from pycbc.vetoes import sgchisq # pycbc_multi_inspiral executable. time_init = time.time() parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument("-V", "--verbose", action="store_true", - help="prints extra debugging information during runtime", - default=False) +add_common_pycbc_options(parser) parser.add_argument("--output", type=str) parser.add_argument("--instruments", nargs="+", type=str, required=True, help="List of instruments to analyze.") @@ -352,10 +351,16 @@ with ctx: } network_names = sorted(network_out_vals.keys()) event_mgr = EventManagerCoherent( - args, args.instruments, ifo_names, - [ifo_out_types[n] for n in ifo_names], network_names, + args, + args.instruments, + ifo_names, + [ifo_out_types[n] for n in ifo_names], + network_names, [network_out_types[n] for n in network_names], - segments=segments, time_slides=time_slides) + segments=segments, + time_slides=time_slides, + gating_info={det: strain_dict[det].gating_info for det in strain_dict} + ) # Template bank: filtering and thinning logging.info("Read in template bank") @@ -697,6 +702,9 @@ with ctx: event_mgr.cluster_template_network_events( 'time_index', 'reweighted_snr', cluster_window, slide=slide) event_mgr.finalize_template_events() + +logging.info('Writing output') event_mgr.write_events(args.output) + logging.info("Finished") logging.info("Time to complete analysis: %d", int(time.time() - time_init)) diff --git a/bin/pycbc_optimal_snr b/bin/pycbc_optimal_snr index 7e99bed39f4..5511c06880f 100644 --- a/bin/pycbc_optimal_snr +++ b/bin/pycbc_optimal_snr @@ -27,12 +27,14 @@ import argparse import multiprocessing import numpy as np import h5py + +from ligo.lw import utils as ligolw_utils +from ligo.lw import lsctables + import pycbc import pycbc.inject import pycbc.psd import pycbc.version -from ligo.lw import utils as ligolw_utils -from ligo.lw import lsctables from pycbc.filter import sigma, make_frequency_series from pycbc.types import TimeSeries, FrequencySeries, zeros, float32, \ MultiDetOptionAction, load_frequencyseries @@ -117,6 +119,7 @@ def get_gc_end_time(injection): if __name__ == '__main__': parser = argparse.ArgumentParser(description=__doc__) + pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action=pycbc.version.Version) parser.add_argument('--input-file', '-i', dest='injection_file', required=True, @@ -159,11 +162,6 @@ if __name__ == '__main__': parser.add_argument('--ignore-waveform-errors', action='store_true', help='Ignore errors in waveform generation and keep ' 'the corresponding column unchanged') - parser.add_argument('--verbose', action='store_true', - help='Print which injection is going to be attempted ' - 'next and when it completes. Use this to debug ' - 'misbehaving approximants which crash or quit ' - 'the Python interpreter') parser.add_argument('--progress', action='store_true', help='Show a progress bar (requires tqdm)') psd_group = pycbc.psd.insert_psd_option_group_multi_ifo(parser) @@ -187,10 +185,7 @@ if __name__ == '__main__': if not opts.time_varying_psds: pycbc.psd.verify_psd_options_multi_ifo(opts, parser, detectors) - log_level = logging.INFO if opts.verbose else logging.WARN - logging.basicConfig(level=log_level, - format='%(asctime)s %(message)s', - datefmt='%Y-%m-%d %H:%M:%S') + pycbc.init_logging(opts.verbose) seg_len = opts.seg_length sample_rate = opts.sample_rate @@ -240,20 +235,27 @@ if __name__ == '__main__': psd = psds[det](injection_time) if psd is None: continue - logging.info('Trying injection %s at %s', inj.simulation_id, det) + logging.debug('Trying injection %s at %s', inj.simulation_id, det) try: wave = get_injection(injections, det, injection_time, simulation_id=inj.simulation_id) except Exception as e: if opts.ignore_waveform_errors: - logging.warn('%s: waveform generation failed, skipping (%s)', - inj.simulation_id, e) + logging.debug( + '%s: waveform generation failed, skipping (%s)', + inj.simulation_id, + e + ) continue else: logging.error('%s: waveform generation failed with the ' 'following exception', inj.simulation_id) raise - logging.info('Injection %s at %s completed', inj.simulation_id, det) + logging.debug( + 'Injection %s at %s completed', + inj.simulation_id, + det + ) sval = sigma(wave, psd=psd, low_frequency_cutoff=f_low) if ligolw: setattr(inj, column, sval) diff --git a/bin/pycbc_optimize_snr b/bin/pycbc_optimize_snr index 5094a02dc62..fb36ba6c721 100755 --- a/bin/pycbc_optimize_snr +++ b/bin/pycbc_optimize_snr @@ -5,14 +5,14 @@ import os import argparse import logging +import numpy +import h5py # we will make plots on a likely headless machine, so make sure matplotlib's # backend is set appropriately from matplotlib import use as mpl_use_backend mpl_use_backend('agg') -import numpy -import h5py import pycbc from pycbc import ( fft, scheme, version @@ -27,9 +27,9 @@ from pycbc.live import snr_optimizer parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action='version', version=version.git_verbose_msg) -parser.add_argument('--verbose', action='store_true') parser.add_argument('--params-file', required=True, help='Location of the attributes file created by PyCBC ' 'Live') @@ -45,12 +45,12 @@ parser.add_argument('--psd-files', type=str, nargs='+', 'by PyCBC Live.') parser.add_argument('--approximant', required=True, help='Waveform approximant string.') -parser.add_argument('--snr-threshold', default=4.0, +parser.add_argument('--snr-threshold', type=float, default=4.0, help='If the SNR in ifo X is below this threshold do not ' 'consider it part of the coincidence. Not implemented') -parser.add_argument('--chirp-time-f-lower', default=20., +parser.add_argument('--chirp-time-f-lower', type=float, default=20., help='Starting frequency for chirp time window (Hz).') -parser.add_argument('--chirp-time-window', default=2., +parser.add_argument('--chirp-time-window', type=float, default=2., help='Chirp time window (s).') parser.add_argument('--gracedb-server', metavar='URL', help='URL of GraceDB server API for uploading events. ' diff --git a/bin/pycbc_process_sngls b/bin/pycbc_process_sngls index c740415ecdd..3e5909e449d 100644 --- a/bin/pycbc_process_sngls +++ b/bin/pycbc_process_sngls @@ -6,12 +6,14 @@ import argparse import logging import numpy import h5py + +import pycbc from pycbc.io import SingleDetTriggers from pycbc.events import stat, coinc, veto parser = argparse.ArgumentParser(description=__doc__) - +pycbc.add_common_pycbc_options(parser) parser.add_argument('--single-trig-file', required=True, help='Path to file containing single-detector triggers in ' 'HDF5 format. Required') @@ -41,7 +43,7 @@ stat.insert_statistic_option_group(parser) args = parser.parse_args() -logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO) +pycbc.init_logging(args.verbose) # Munge together SNR cut and any other filter specified snr_filter = 'self.snr>%f' % (args.min_snr) if args.min_snr > 0. else None diff --git a/bin/pycbc_randomize_inj_dist_by_optsnr b/bin/pycbc_randomize_inj_dist_by_optsnr index 1bc4f507c27..c3b3a334209 100644 --- a/bin/pycbc_randomize_inj_dist_by_optsnr +++ b/bin/pycbc_randomize_inj_dist_by_optsnr @@ -14,6 +14,8 @@ from argparse import ArgumentParser from ligo.lw import lsctables from ligo.lw import utils as ligolw_utils from ligo.lw import table + +import pycbc from pycbc.io.ligolw import LIGOLWContentHandler @@ -50,6 +52,7 @@ def get_weights(distribution, r, r1, r2): return min_vol, weight parser = ArgumentParser(description = __doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--xml-file', required = True, help = 'Input LIGOLW file defining injections.') parser.add_argument('--output-file', required=True, @@ -90,10 +93,11 @@ parser.add_argument('--scale-by-chirp-dist', action='store_true', parser.add_argument('--seed', type = int, default = int((time.time()*100)% 1e6), help = 'Set the seed to use for the random number generator. ' \ 'If none specified, will use the current time.') -parser.add_argument('--verbose', action = 'store_true', help = 'Be verbose.') opts = parser.parse_args() +pycbc.init_logging(opts.verbose) + # parse the distributions # beta is the only special one if opts.distribution in ['volume', 'distance', 'logdist']: diff --git a/bin/pycbc_single_template b/bin/pycbc_single_template index 062f2c11815..ee78facdcc1 100755 --- a/bin/pycbc_single_template +++ b/bin/pycbc_single_template @@ -17,7 +17,13 @@ """ Calculate the SNR and CHISQ timeseries for either a chosen template, or a specific Nth loudest coincident event. """ -import sys, logging, argparse, numpy, pycbc, h5py +import sys +import logging +import argparse +import numpy +import h5py + +import pycbc from pycbc import vetoes, psd, waveform, strain, scheme, fft, filter from pycbc.io import WaveformArray from pycbc import events @@ -86,11 +92,10 @@ def select_segments(fname, anal_name, data_name, ifo, time, pad_data): parser = argparse.ArgumentParser(usage='', description="Single template gravitational-wave followup") +pycbc.add_common_pycbc_options(parser) parser.add_argument('--version', action=pycbc.version.Version) parser.add_argument('--output-file', required=True) parser.add_argument('--subtract-template', action='store_true') -parser.add_argument("-V", "--verbose", action="store_true", - help="print extra debugging information", default=False ) parser.add_argument("--low-frequency-cutoff", type=float, help="The low frequency cutoff to use for filtering (Hz)") parser.add_argument("--high-frequency-cutoff", type=float, @@ -182,6 +187,8 @@ scheme.insert_processing_option_group(parser) fft.insert_fft_option_group(parser) opt = parser.parse_args() +pycbc.init_logging(opt.verbose) + # Exclude invalid options if opt.window and opt.trigger_time is None: raise RuntimeError("Can't use --window option without a valid trigger time!") @@ -210,7 +217,6 @@ strain.verify_strain_options(opt, parser) strain.StrainSegments.verify_segment_options(opt, parser) scheme.verify_processing_options(opt, parser) fft.verify_fft_options(opt,parser) -pycbc.init_logging(opt.verbose) ctx = scheme.from_cli(opt) gwstrain = strain.from_cli(opt, pycbc.DYN_RANGE_FAC) diff --git a/bin/pycbc_source_probability_offline b/bin/pycbc_source_probability_offline index 2fc8eaceb84..f12447603ef 100755 --- a/bin/pycbc_source_probability_offline +++ b/bin/pycbc_source_probability_offline @@ -9,12 +9,14 @@ import tqdm import argparse import logging import numpy as np + import pycbc from pycbc.io import hdf from pycbc.pnutils import mass1_mass2_to_mchirp_eta from pycbc import mchirp_area parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument('--trigger-file', required=True) parser.add_argument('--bank-file', required=True) parser.add_argument('--single-detector-triggers', nargs='+', required=True) @@ -27,7 +29,6 @@ parser.add_argument('--ifar-threshold', type=float, default=None, 'above threshold.') parser.add_argument('--include-mass-gap', action='store_true', help='Option to include the Mass Gap region.') -parser.add_argument('--verbose', action='count') parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) mchirp_area.insert_args(parser) diff --git a/bin/pycbc_split_inspinj b/bin/pycbc_split_inspinj index ce38299eda1..28e5ee9c1ad 100644 --- a/bin/pycbc_split_inspinj +++ b/bin/pycbc_split_inspinj @@ -5,12 +5,14 @@ import argparse from ligo.lw import utils as ligolw_utils from ligo.lw import lsctables from itertools import cycle + import pycbc.version from pycbc.io.ligolw import LIGOLWContentHandler, get_table_columns # Parse command line parser = argparse.ArgumentParser() +pycbc.add_common_pycbc_options(parser) parser.add_argument("--version", action="version", version=pycbc.version.git_verbose_msg) group = parser.add_mutually_exclusive_group(required=True) @@ -30,7 +32,10 @@ if args.output_files and args.output_dir: # Read in input file xmldoc = ligolw_utils.load_filename( - args.input_file, verbose=True, contenthandler=LIGOLWContentHandler) + args.input_file, + verbose=args.verbose, + contenthandler=LIGOLWContentHandler +) tabletype = lsctables.SimInspiralTable allinjs = tabletype.get_table(xmldoc) diff --git a/bin/pycbc_splitbank b/bin/pycbc_splitbank index 42d2e6ba437..8919c31ac2f 100644 --- a/bin/pycbc_splitbank +++ b/bin/pycbc_splitbank @@ -29,9 +29,12 @@ import argparse from numpy import random, ceil + from ligo.lw import ligolw from ligo.lw import lsctables from ligo.lw import utils as ligolw_utils + +import pycbc from pycbc import version from pycbc.io.ligolw import LIGOLWContentHandler, create_process_table from pycbc.conversions import mchirp_from_mass1_mass2 @@ -46,6 +49,7 @@ __program__ = "pycbc_splitbank" parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--version', action='version', version=version.git_verbose_msg) +pycbc.add_common_pycbc_options(parser) group = parser.add_mutually_exclusive_group(required=True) group.add_argument('--templates-per-bank', metavar='SAMPLES', help='number of templates in the output banks', type=int) @@ -61,8 +65,6 @@ group.add_argument("-O", "--output-filenames", nargs='*', default=None, parser.add_argument("-o", "--output-prefix", default=None, help="Prefix to add to the template bank name (name becomes output#.xml[.gz])" ) -parser.add_argument("-V", "--verbose", action="store_true", - help="Print extra debugging information", default=False ) parser.add_argument("-t", "--bank-file", metavar='INPUT_FILE', help='Template bank to split', required=True) parser.add_argument("--sort-frequency-cutoff", diff --git a/bin/pycbc_upload_xml_to_gracedb b/bin/pycbc_upload_xml_to_gracedb index a922d0923cb..df9fc57848a 100755 --- a/bin/pycbc_upload_xml_to_gracedb +++ b/bin/pycbc_upload_xml_to_gracedb @@ -23,24 +23,26 @@ Take a coinc xml file containing multiple events and upload to gracedb. import os import argparse import logging -from ligo.gracedb.rest import GraceDb import numpy as np import h5py import matplotlib matplotlib.use('agg') -import pycbc +from ligo.gracedb.rest import GraceDb import lal import lal.series from ligo.lw import lsctables from ligo.lw import utils as ligolw_utils from ligo.segments import segment, segmentlist + +import pycbc from pycbc.io.live import gracedb_tag_with_version from pycbc.io.ligolw import LIGOLWContentHandler from pycbc.psd import interpolate from pycbc.types import FrequencySeries from pycbc.results import generate_asd_plot + def check_gracedb_for_event(gdb_handle, query, far): """ Check if there is an event in gracedb with the queried string @@ -68,6 +70,7 @@ logging.basicConfig(format='%(asctime)s:%(levelname)s : %(message)s', level=logging.INFO) parser = argparse.ArgumentParser(description=__doc__) +pycbc.add_common_pycbc_options(parser) parser.add_argument("--psd-files", nargs='+', required=True, help='HDF file(s) containing the PSDs to upload') parser.add_argument('--input-file', required=True, type=str, @@ -109,6 +112,8 @@ parser.add_argument('--generate-plots', action='store_true', args = parser.parse_args() +pycbc.init_logging(args.verbose) + if args.production_server: gracedb = GraceDb() else: diff --git a/bin/pygrb/pycbc_grb_inj_finder b/bin/pygrb/pycbc_grb_inj_finder index 623e45c231e..230e09d3bc9 100644 --- a/bin/pygrb/pycbc_grb_inj_finder +++ b/bin/pygrb/pycbc_grb_inj_finder @@ -42,6 +42,7 @@ from ligo.segments.utils import fromsegwizard from pycbc import __version__ from pycbc.inject import InjectionSet from pycbc.io import FieldArray +from pycbc.results.pygrb_postprocessing_utils import template_hash_to_id __author__ = "Duncan Macleod " @@ -96,6 +97,8 @@ def read_hdf5_triggers(inputfiles, verbose=False): elif name.startswith("network") and \ NETWORK_IFO_EVENT_ID_REGEX.match(name): return + elif "template_id" in name: + return else: shape = obj.shape dtype = obj.dtype @@ -126,7 +129,7 @@ def read_hdf5_triggers(inputfiles, verbose=False): # copy dataset contents for filename in inputfiles: with h5py.File(filename, 'r') as h5in: - # Skipping emtpy files + # Skipping empty files if len(h5in['network']['end_time_gc'][:]): for dset in datasets: data = h5in[dset][:] @@ -174,6 +177,11 @@ parser.add_argument( nargs="+", help="path(s) to input injection file(s)", ) +parser.add_argument( + "--bank-file", + required=True, + help="full template bank file for mapping hashes to ids", +) parser.add_argument( "-o", "--output-dir", @@ -329,12 +337,21 @@ with h5py.File(outfilename, "w") as h5out: for x in out_trigs: xgroup = h5out.create_group(x) for col in out_trigs[x]: + if "template_id" == col: + continue xgroup.create_dataset( col, data=out_trigs[x][col], compression="gzip", compression_opts=9, ) - + # map and write template ids + ids = template_hash_to_id(h5out, args.bank_file) + h5out.create_dataset( + "network/template_id", + data=ids, + compression="gzip", + compression_opts=9, + ) if args.verbose: print("Found/missed written to {}".format(outfilename)) diff --git a/bin/pygrb/pycbc_grb_trig_combiner b/bin/pygrb/pycbc_grb_trig_combiner index 84cdf0b4ca5..e1687a3c237 100644 --- a/bin/pygrb/pycbc_grb_trig_combiner +++ b/bin/pygrb/pycbc_grb_trig_combiner @@ -36,6 +36,7 @@ from ligo import segments from ligo.segments.utils import fromsegwizard from pycbc import __version__ +from pycbc.results.pygrb_postprocessing_utils import template_hash_to_id __author__ = "Duncan Macleod " @@ -55,6 +56,8 @@ def _init_combined_file(h5file, attributes, datasets, **compression_kw): h5file.attrs.update(attributes) # create datasets for dset, (shape, dtype) in datasets.items(): + if "template_id" in dset: + continue h5file.create_dataset(dset, shape=shape, dtype=dtype, **compression_kw) @@ -207,6 +210,11 @@ def merge_hdf5_files(inputfiles, outputfile, verbose=False, **compression_kw): else: ifo_event_id_pos = None + # Template ID becomes unhelpful in the merged file. + # This will be handled separately. + if "template_id" in dset: + continue + # merge dataset position[dset] += _merge_dataset( h5in[dset], h5out[dset], position[dset], @@ -218,6 +226,12 @@ def merge_hdf5_files(inputfiles, outputfile, verbose=False, **compression_kw): while search_datasets: dataset_names.remove(search_datasets.pop()) + template_ids = template_hash_to_id(trigger_file=h5out, bank_path=args.bank_file) + h5out.create_dataset("network/template_id", + data=template_ids, + compression="gzip", + compression_opts=9) + if verbose: print("Merged triggers written to {}".format(outputfile)) return outputfile @@ -423,6 +437,11 @@ parser.add_argument( metavar="TRIGGER FILE", help="read in listed trigger files", ) +parser.add_argument( + "--bank-file", + required=True, + help="full template bank file for mapping hashes to ids", +) parser.add_argument( "-o", "--output-dir", diff --git a/bin/pygrb/pycbc_make_offline_grb_workflow b/bin/pygrb/pycbc_make_offline_grb_workflow index 3dd1ee9d831..59160f4f79c 100644 --- a/bin/pygrb/pycbc_make_offline_grb_workflow +++ b/bin/pygrb/pycbc_make_offline_grb_workflow @@ -273,9 +273,14 @@ all_files.extend(datafind_veto_files) # TEMPLATE BANK AND SPLIT BANK bank_files = _workflow.setup_tmpltbank_workflow(wflow, sciSegs, datafind_files, df_dir) +# Check there are not multiple banks +if len(bank_files) > 1: + raise NotImplementedError("Multiple banks not supported") +full_bank_file = bank_files[0] +# Note: setup_splittable_workflow requires a FileList as input splitbank_files = _workflow.setup_splittable_workflow(wflow, bank_files, df_dir, tags=["inspiral"]) -all_files.extend(bank_files) +all_files.append(full_bank_file) all_files.extend(splitbank_files) # INJECTIONS @@ -411,7 +416,7 @@ if post_proc_method == "PYGRB_OFFLINE": # The content and structure of pp_files are described in the definition # of _workflow.setup_pygrb_pp_workflow pp_files = _workflow.setup_pygrb_pp_workflow(wflow, pp_dir, seg_dir, - sciSegs[ifo][0], inspiral_files, + sciSegs[ifo][0], full_bank_file, inspiral_files, injs, inj_insp_files, inj_tags) sec_name = 'workflow-pygrb_pp_workflow' if not wflow.cp.has_section(sec_name): diff --git a/bin/pygrb/pycbc_pygrb_efficiency b/bin/pygrb/pycbc_pygrb_efficiency index 9c5b48e4660..b3755c228c0 100644 --- a/bin/pygrb/pycbc_pygrb_efficiency +++ b/bin/pygrb/pycbc_pygrb_efficiency @@ -24,6 +24,7 @@ import sys import os import logging +import h5py import matplotlib.pyplot as plt from matplotlib import rc import numpy as np @@ -34,10 +35,7 @@ from pycbc.detector import Detector from pycbc import init_logging from pycbc.results import save_fig_with_metadata from pycbc.results import pygrb_postprocessing_utils as ppu -try: - from glue.ligolw import lsctables -except ImportError: - pass +from pycbc.conversions import mchirp_from_mass1_mass2 plt.switch_backend('Agg') rc("image") @@ -98,18 +96,26 @@ parser.add_argument("-g", "--glitch-check-factor", action="store", "exclusion efficiencies this value is multiplied " + "to the offsource around the injection trigger to " + "determine if it is just a loud glitch.") +parser.add_argument("--found-missed-file", action="store", type=str, + default=None, help="Location of found/missed " + "injections and trigger file") parser.add_argument("--injection-set-name", action="store", type=str, default="", help="Name of the injection set to be " + "used in the plot title.") parser.add_argument("-C", "--cluster-window", action="store", type=float, default=0.1, help="The cluster window used " + "to cluster triggers in time.") +parser.add_argument("--bank-file", action="store", type=str, + help="Location of the full template bank used.") ppu.pygrb_add_injmc_opts(parser) ppu.pygrb_add_bestnr_opts(parser) opts = parser.parse_args() init_logging(opts.verbose, format="%(asctime)s: %(levelname)s: %(message)s") +# Load bank file +bank_file = h5py.File(opts.bank_file, 'r') + # Check options do_injections = opts.found_missed_file is not None @@ -165,8 +171,8 @@ ifos, vetoes = ppu.extract_ifos_and_vetoes(trig_file, opts.veto_files, # Load triggers, time-slides, and segment dictionary logging.info("Loading triggers.") -trigs = ppu.load_xml_table(trig_file, lsctables.MultiInspiralTable.tableName) -logging.info("%d triggers loaded.", len(trigs)) +trigs = ppu.load_triggers(trig_file, vetoes) +logging.info("%d triggers loaded.", trigs['network/event_id'].len()) logging.info("Loading timeslides.") slide_dict = ppu.load_time_slides(trig_file) logging.info("Loading segments.") @@ -181,7 +187,7 @@ logging.info("%d trials generated.", total_trials) # Extract basic trigger properties and store as dictionaries trig_time, trig_snr, trig_bestnr = \ - ppu.extract_basic_trig_properties(trial_dict, trigs, slide_dict, segment_dict, opts) + ppu.extract_basic_trig_properties(trial_dict, trigs, slide_dict, segment_dict) # Calculate BestNR values and maximum time_veto_max_bestnr = {} @@ -219,6 +225,17 @@ offsource_trigs.reverse() max_bestnr, _, full_time_veto_max_bestnr =\ ppu.max_median_stat(slide_dict, time_veto_max_bestnr, trig_bestnr, total_trials) +# ========================== +# Calculate template chirp masses from bank +# ========================== +logging.info('Reading template chirp masses') +with h5py.File(opts.bank_file, 'r') as bank_file: + template_mchirps = mchirp_from_mass1_mass2( + bank_file['mass1'][:], + bank_file['mass2'][:] + ) + + # ======================= # Load on source triggers @@ -226,47 +243,38 @@ max_bestnr, _, full_time_veto_max_bestnr =\ if onsource_file: # Get trigs - on_trigs = ppu.load_xml_table(onsource_file, "multi_inspiral") - logging.info("%d onsource triggers loaded.", len(on_trigs)) + on_trigs = ppu.load_triggers(onsource_file, vetoes) + logging.info("%d onsource triggers loaded.", on_trigs['network/event_id'].len()) - # Separate off chirp mass column - on_mchirp = on_trigs.get_column('mchirp') + # Calculate chirp mass values + on_mchirp = template_mchirps[on_trigs['network/template_id']] # Set loudest event arrays - #loud_on_bestnr_trigs = None loud_on_bestnr = 0 - #loud_on_fap = 1 - # Calculate BestNRs and record loudest trig by BestNR + # Retrieve BestNRs and record loudest trig by BestNR. + # Get indices of loudest triggers and pick the loudest. if on_trigs: - on_trigs_bestnrs = ppu.get_bestnrs(on_trigs, - q=chisq_index, n=chisq_nhigh, - null_thresh=null_thresh, - snr_threshold=snr_thresh, - sngl_snr_threshold=sngl_snr_thresh, - chisq_threshold=new_snr_thresh, - null_grad_thresh=null_grad_thresh, - null_grad_val=null_grad_val) - loud_on_bestnr_trigs, loud_on_bestnr = \ - np.asarray([[x, y] for y, x in sorted(zip(on_trigs_bestnrs, on_trigs), - key=lambda elem: elem[0], - reverse=True)])[0] + loud_on_bestnr_idx = np.argmax(on_trigs['network/reweighted_snr']) + loud_on_bestnr = np.max(on_trigs['network/reweighted_snr']) + + # Convert to float for output + loud_on_bestnr = float(loud_on_bestnr) # If the loudest event has bestnr = 0, there is no event at all! if loud_on_bestnr == 0: - loud_on_bestnr_trigs = None + loud_on_bestnr_idx = None loud_on_fap = 1 logging.info("Onsource analysed.") - if loud_on_bestnr_trigs: - trig = loud_on_bestnr_trigs + if loud_on_bestnr_idx is not None: num_trials_louder = 0 tot_off_snr = np.array([]) for slide_id in slide_dict: - num_trials_louder += sum(time_veto_max_bestnr[slide_id] > \ + num_trials_louder += sum(time_veto_max_bestnr[slide_id] > loud_on_bestnr) - tot_off_snr = np.concatenate([tot_off_snr,\ + tot_off_snr = np.concatenate([tot_off_snr, time_veto_max_bestnr[slide_id]]) #fap_test = sum(tot_off_snr > loud_on_bestnr)/total_trials loud_on_fap = num_trials_louder/total_trials @@ -274,7 +282,7 @@ if onsource_file: else: tot_off_snr = np.array([]) for slide_id in slide_dict: - tot_off_snr = np.concatenate([tot_off_snr,\ + tot_off_snr = np.concatenate([tot_off_snr, time_veto_max_bestnr[slide_id]]) med_snr = np.median(tot_off_snr) #loud_on_fap = sum(tot_off_snr > med_snr)/total_trials @@ -287,55 +295,38 @@ if do_injections: sites = [ifo[0] for ifo in ifos] # Triggers and injections recovered in some form: discard ones at vetoed times - found_trigs = ppu.load_injections(found_file, vetoes) - found_injs = ppu.load_injections(found_file, vetoes, sim_table=True) + # This injs object contains the information about found/missed injections AND triggers + injs = ppu.load_triggers(found_missed_file, vetoes) logging.info("Missed/found injections/triggers loaded.") - # Extract columns of found injections - found_inj = {} - found_inj['ra'] = np.asarray(found_injs.get_column('longitude')) - found_inj['dec'] = np.asarray(found_injs.get_column('latitude')) - found_inj['time'] = np.asarray(found_injs.get_column('geocent_end_time')) +\ - np.asarray(found_injs.get_column('geocent_end_time_ns')*\ - 10**-9) - found_inj['dist'] = np.asarray(found_injs.get_column('distance')) - - # Extract columns of triggers - found_trig = {} - found_trig['mchirp'] = np.asarray(found_trigs.get_column('mchirp')) - found_trig['bestnr'] = ppu.get_bestnrs(found_trigs, q=chisq_index, - n=chisq_nhigh, - null_thresh=null_thresh, - snr_threshold=snr_thresh, - sngl_snr_threshold=sngl_snr_thresh, - chisq_threshold=new_snr_thresh, - null_grad_thresh=null_grad_thresh, - null_grad_val=null_grad_val) + # Calculate quantities not included in trigger files, such as chirp mass + found_trig_mchirp = template_mchirps[injs['network/template_id']] + # Construct conditions for injection: # 1) found louder than background, - zero_fap = np.zeros(len(found_injs)).astype(bool) - zero_fap_cut = found_trig['bestnr'] > max_bestnr + zero_fap = np.zeros(len(injs['network/end_time_gc'])).astype(bool) + zero_fap_cut = injs['network/reweighted_snr'][:] > max_bestnr zero_fap = zero_fap | (zero_fap_cut) # 2) found (bestnr > 0) but not louder than background (non-zero FAP) - nonzero_fap = ~zero_fap & (found_trig['bestnr'] != 0) + nonzero_fap = ~zero_fap & (injs['network/reweighted_snr'] != 0) # 3) missed after being recovered (i.e., vetoed) are not used here # missed = (~zero_fap) & (~nonzero_fap) # Non-zero FAP triggers (g_ifar) g_ifar = {} - g_ifar['bestnr'] = found_trig['bestnr'][nonzero_fap] + g_ifar['bestnr'] = injs['network/reweighted_snr'][nonzero_fap] g_ifar['stat'] = np.zeros([len(g_ifar['bestnr'])]) for ix, (mc, bestnr) in \ - enumerate(zip(found_trig['mchirp'][nonzero_fap], g_ifar['bestnr'])): + enumerate(zip(found_trig_mchirp[nonzero_fap], g_ifar['bestnr'])): g_ifar['stat'][ix] = (full_time_veto_max_bestnr > bestnr).sum() g_ifar['stat'] = g_ifar['stat'] / total_trials # Set the sigma values - inj_sigma = found_trigs.get_sigmasqs() + inj_sigma = {ifo: injs[f'{ifo}/sigmasq'][:] for ifo in ifos} # If the sigmasqs are not populated, we can still do calibration errors, # but only in the 1-detector case for ifo in ifos: @@ -354,10 +345,10 @@ if do_injections: for ifo in ifos: antenna = Detector(ifo) f_resp[ifo] = ppu.get_antenna_responses(antenna, - found_inj['ra'], found_inj['dec'], - found_inj['time']) + injs['found/ra'][:], injs['found/dec'][:], + injs['found/tc'][:]) - inj_sigma_mult = (np.asarray(list(inj_sigma.values())) *\ + inj_sigma_mult = (np.asarray(list(inj_sigma.values())) * np.asarray(list(f_resp.values()))) inj_sigma_tot = inj_sigma_mult[0, :] @@ -368,19 +359,14 @@ if do_injections: for ifo in ifos: inj_sigma_mean[ifo] = ((inj_sigma[ifo]*f_resp[ifo])/inj_sigma_tot).mean() - logging.info("%d found injections analysed.", len(found_injs)) - - # Missed injections (ones not recovered at all) - missed_injs = ppu.load_injections(found_missed_file, vetoes, sim_table=True) + logging.info("%d found injections analysed.", len(injs['found/tc'])) - # Process missed injections 'missed_inj' - missed_inj = {} - missed_inj['dist'] = np.asarray(missed_injs.get_column('distance')) + # Process missed injections (injs['missed']) - logging.info("%d missed injections analysed.", len(missed_injs)) + logging.info("%d missed injections analysed.", len(injs['missed/tc'])) # Create new set of injections for efficiency calculations - total_injs = len(found_injs) + len(missed_injs) + total_injs = len(injs['found/distance']) + len(injs['missed/distance']) long_inj = {} long_inj['dist'] = stats.uniform.rvs(size=total_injs) * (upper_dist-lower_dist) +\ upper_dist @@ -388,10 +374,10 @@ if do_injections: logging.info("%d long distance injections created.", total_injs) # Set distance bins and data arrays - dist_bins = zip(np.arange(lower_dist, upper_dist + (upper_dist-lower_dist),\ - (upper_dist-lower_dist)/num_bins),\ - np.arange(lower_dist, upper_dist + (upper_dist-lower_dist),\ - (upper_dist-lower_dist)/num_bins) +\ + dist_bins = zip(np.arange(lower_dist, upper_dist + (upper_dist-lower_dist), + (upper_dist-lower_dist)/num_bins), + np.arange(lower_dist, upper_dist + (upper_dist-lower_dist), + (upper_dist-lower_dist)/num_bins) + (upper_dist-lower_dist)/num_bins) dist_bins = list(dist_bins) @@ -405,7 +391,7 @@ if do_injections: found_on_bestnr[key] = np.zeros(num_dist_bins_plus_one) # Construct FAP list for all found injections - inj_fap = np.zeros(len(found_injs)) + inj_fap = np.zeros(len(injs['found/distance'])) inj_fap[nonzero_fap] = g_ifar['stat'] # Calculate the amplitude error @@ -428,18 +414,17 @@ if do_injections: # NOTE: the loop on num_mc_injs would fill up the *_inj['dist_mc']'s at the # same time, so filling them up sequentially will vary the numbers a little # (this is an MC, order of operations matters!) - found_inj['dist_mc'] = ppu.mc_cal_wf_errs(num_mc_injs, found_inj['dist'], + found_inj_dist_mc = ppu.mc_cal_wf_errs(num_mc_injs, injs['found/distance'], cal_error, wav_err, max_dc_cal_error) - missed_inj['dist_mc'] = ppu.mc_cal_wf_errs(num_mc_injs, missed_inj['dist'], + missed_inj_dist_mc = ppu.mc_cal_wf_errs(num_mc_injs, injs['missed/distance'], cal_error, wav_err, max_dc_cal_error) long_inj['dist_mc'] = ppu.mc_cal_wf_errs(num_mc_injs, long_inj['dist'], cal_error, wav_err, max_dc_cal_error) - logging.info("MC injection set distributed with %d iterations.",\ + logging.info("MC injection set distributed with %d iterations.", num_mc_injs) # Check injections against on source - more_sig_than_onsource = np.ndarray(len(found_injs)) if onsource_file: more_sig_than_onsource = (inj_fap <= loud_on_fap) else: @@ -447,29 +432,29 @@ if do_injections: distance_count = np.zeros(len(dist_bins)) - found_trig['max_bestnr'] = np.empty(len(found_trig['mchirp'])) - found_trig['max_bestnr'].fill(max_bestnr) + found_trig_max_bestnr = np.empty(len(injs['network/event_id'])) + found_trig_max_bestnr.fill(max_bestnr) - max_bestnr_cut = (found_trig['bestnr'] > found_trig['max_bestnr']) + max_bestnr_cut = (injs['network/reweighted_snr'] > found_trig_max_bestnr) # Check louder than on source - found_trig['loud_on_bestnr'] = np.empty(len(found_trig['mchirp'])) + found_trig_loud_on_bestnr = np.empty(len(injs['network/event_id'])) if onsource_file: - found_trig['loud_on_bestnr'].fill(loud_on_bestnr) + found_trig_loud_on_bestnr.fill(loud_on_bestnr) else: - found_trig['loud_on_bestnr'].fill(med_snr) - on_bestnr_cut = found_trig['bestnr'] > found_trig['loud_on_bestnr'] + found_trig_loud_on_bestnr.fill(med_snr) + on_bestnr_cut = injs['network/reweighted_snr'] > found_trig_loud_on_bestnr # Check whether injection is found for the purposes of exclusion # distance calculation. # Found: if louder than all on source # Missed: if not louder than loudest on source found_excl = on_bestnr_cut & (more_sig_than_onsource) & \ - (found_trig['bestnr'] != 0) + (injs['network/reweighted_snr'] != 0) # If not missed, double check bestnr against nearby triggers near_test = np.zeros((found_excl).sum()).astype(bool) - for j, (t, bestnr) in enumerate(zip(found_inj['time'][found_excl],\ - found_trig['bestnr'][found_excl])): + for j, (t, bestnr) in enumerate(zip(injs['found/tc'][found_excl], + injs['network/reweighted_snr'][found_excl])): # 0 is the zero-lag timeslide near_bestnr = trig_bestnr[0]\ [np.abs(trig_time[0]-t) < cluster_window] @@ -486,10 +471,10 @@ if do_injections: # Loop over the distance bins for j, dist_bin in enumerate(dist_bins): # Construct distance cut - found_dist_cut = (dist_bin[0] <= found_inj['dist_mc'][k, :]) &\ - (found_inj['dist_mc'][k, :] < dist_bin[1]) - missed_dist_cut = (dist_bin[0] <= missed_inj['dist_mc'][k, :]) &\ - (missed_inj['dist_mc'][k, :] < dist_bin[1]) + found_dist_cut = (dist_bin[0] <= found_inj_dist_mc[k, :]) &\ + (found_inj_dist_mc[k, :] < dist_bin[1]) + missed_dist_cut = (dist_bin[0] <= missed_inj_dist_mc[k, :]) &\ + (missed_inj_dist_mc[k, :] < dist_bin[1]) long_dist_cut = (dist_bin[0] <= long_inj['dist_mc'][k, :]) &\ (long_inj['dist_mc'][k, :] < dist_bin[1]) @@ -530,7 +515,7 @@ if do_injections: # Calculate error bars for efficiency/distance plots and datafiles # using max BestNR of background - yerr_low_mc, yerr_high_mc, fraction_mc = efficiency_with_errs(\ + yerr_low_mc, yerr_high_mc, fraction_mc = efficiency_with_errs( found_max_bestnr['mc'], num_injections['mc'], num_mc_injs=num_mc_injs) yerr_low_no_mc, yerr_high_no_mc, fraction_no_mc = efficiency_with_errs(\ found_max_bestnr['no_mc'], num_injections['no_mc']) diff --git a/bin/pygrb/pycbc_pygrb_plot_injs_results b/bin/pygrb/pycbc_pygrb_plot_injs_results index 22915b945cd..e79494f23e3 100644 --- a/bin/pygrb/pycbc_pygrb_plot_injs_results +++ b/bin/pygrb/pycbc_pygrb_plot_injs_results @@ -100,7 +100,7 @@ def load_mass_data(input_file_handle, key, tag): return data -def load_sky_error_data(input_file_handle, key, tag, trig_data): +def load_sky_error_data(input_file_handle, key, tag): """Extract data related to sky_error from raw injection and trigger data""" if tag == 'missed': # Missed injections are assigned null values @@ -110,8 +110,8 @@ def load_sky_error_data(input_file_handle, key, tag, trig_data): inj['ra'] = input_file_handle[tag+'/ra'][:] inj['dec'] = input_file_handle[tag+'/dec'][:] trig = {} - trig['ra'] = np.rad2deg(trig_data['network/longitude'][:]) - trig['dec'] = np.rad2deg(trig_data['network/latitude'][:]) + trig['ra'] = input_file_handle['network/ra'][:] + trig['dec'] = input_file_handle['network/dec'][:] data = np.arccos(np.cos(inj['dec'] - trig['dec']) - np.cos(inj['dec']) * np.cos(trig['dec']) * (1 - np.cos(inj['ra'] - trig['ra']))) @@ -128,7 +128,7 @@ easy_keys = ['distance', 'mass1', 'mass2', 'polarization', 'spin2_azimuthal', 'spin2_polar', 'dec', 'ra', 'phi_ref'] # Function to extract desired data from an injection file -def load_data(input_file_handle, keys, tag, trig_data): +def load_data(input_file_handle, keys, tag): """Create a dictionary containing the data specified by the list of keys extracted from an injection file""" @@ -148,7 +148,7 @@ def load_data(input_file_handle, keys, tag, trig_data): data_dict[key] = load_incl_data(input_file_handle, key, tag) elif key == 'sky_error': data_dict[key] = load_sky_error_data(input_file_handle, key, - tag, trig_data) + tag) except KeyError: #raise NotImplemented(key+' not allowed: returning empty entry') logging.info(key+' not allowed yet') @@ -247,10 +247,10 @@ max_reweighted_snr = max(trig_data['network/reweighted_snr'][:]) # Post-process injections # ======================= # Extract the necessary data from the missed injections for the plot -missed_inj = load_data(f, [x_qty, y_qty], 'missed', trig_data) +missed_inj = load_data(f, [x_qty, y_qty], 'missed') # Extract the necessary data from the found injections for the plot -found_inj = load_data(f, [x_qty, y_qty], 'found', trig_data) +found_inj = load_data(f, [x_qty, y_qty], 'found') # Injections found surviving vetoes found_after_vetoes = found_inj if 'found_after_vetoes' not in f.keys() else f['found_after_vetoes'] diff --git a/bin/workflows/pycbc_make_offline_search_workflow b/bin/workflows/pycbc_make_offline_search_workflow index a7fc7daf5c8..25fc1e627ba 100755 --- a/bin/workflows/pycbc_make_offline_search_workflow +++ b/bin/workflows/pycbc_make_offline_search_workflow @@ -57,6 +57,92 @@ def ifo_combos(ifos): for ifocomb in combinations: yield ifocomb +def finalize(container, workflow, finalize_workflow): + # Create the final log file + log_file_html = wf.File(workflow.ifos, 'WORKFLOW-LOG', workflow.analysis_time, + extension='.html', directory=rdir['workflow']) + + gen_file_html = wf.File(workflow.ifos, 'WORKFLOW-GEN', workflow.analysis_time, + extension='.html', directory=rdir['workflow']) + + # Create a page to contain a dashboard link + dashboard_file = wf.File(workflow.ifos, 'DASHBOARD', workflow.analysis_time, + extension='.html', directory=rdir['workflow']) + dashboard_str = """

Pegasus Dashboard Page

""" + kwds = {'title': "Pegasus Dashboard", + 'caption': "Link to Pegasus Dashboard", + 'cmd': "PYCBC_SUBMIT_DAX_ARGV", } + save_fig_with_metadata(dashboard_str, dashboard_file.storage_path, **kwds) + + # Create pages for the submission script to write data + wf.makedir(rdir['workflow/dax']) + wf.makedir(rdir['workflow/input_map']) + wf.makedir(rdir['workflow/output_map']) + wf.makedir(rdir['workflow/planning']) + + wf.make_results_web_page(finalize_workflow, os.path.join(os.getcwd(), + rdir.base)) + + container += workflow + container += finalize_workflow + + container.add_subworkflow_dependancy(workflow, finalize_workflow) + + container.save() + + open_box_result_path = rdir['open_box_result'] + if os.path.exists(open_box_result_path): + os.chmod(open_box_result_path, 0o0700) + else: + pass + + logging.info("Written dax.") + + # Close the log and flush to the html file + logging.shutdown() + with open(wf_log_file.storage_path, "r") as logfile: + logdata = logfile.read() + log_str = """ +

Workflow generation script created workflow in output directory: %s

+

Workflow name is: %s

+

Workflow generation script run on host: %s

+ 
%s
+ """ % (os.getcwd(), args.workflow_name, socket.gethostname(), logdata) + kwds = {'title': 'Workflow Generation Log', + 'caption': "Log of the workflow script %s" % sys.argv[0], + 'cmd': ' '.join(sys.argv), } + save_fig_with_metadata(log_str, log_file_html.storage_path, **kwds) + + # Add the command line used to a specific file + args_to_output = [sys.argv[0]] + for arg in sys.argv[1:]: + if arg.startswith('--'): + # This is an option, add tab + args_to_output.append(' ' + arg) + else: + # This is a parameter, add two tabs + args_to_output.append(' ' + arg) + + gen_str = '
' + ' \\\n'.join(args_to_output) + '
' + kwds = {'title': 'Workflow Generation Command', + 'caption': "Command used to generate the workflow.", + 'cmd': ' '.join(sys.argv), } + save_fig_with_metadata(gen_str, gen_file_html.storage_path, **kwds) + layout.single_layout(rdir['workflow'], ([dashboard_file, gen_file_html, log_file_html])) + sys.exit(0) + +def check_stop(job_name, container, workflow, finalize_workflow): + #This function will finalize the workflow and stop it at the job specified. + #Under [workflow] supply the option stop-after = and the job name you want the workflow to stop at. Current options are [inspiral, hdf_trigger_merge, statmap] + if workflow.cp.has_option('workflow', 'stop-after'): + stop_after = workflow.cp.get('workflow', 'stop-after') + if stop_after == job_name: + logging.info("Search worflow will be stopped after " + str(job_name)) + finalize(container, workflow, finalize_workflow) + else: + pass + + parser = argparse.ArgumentParser(description=__doc__[1:]) parser.add_argument('--version', action='version', version=__version__) parser.add_argument('--verbose', action='count', @@ -183,10 +269,15 @@ ind_insps = insps = wf.setup_matchedfltr_workflow(workflow, analyzable_segs, datafind_files, splitbank_files_fd, output_dir, tags=['full_data']) +#check to see if workflow should stop at inspiral jobs. +check_stop('inspiral', container, workflow, finalize_workflow) + insps = wf.merge_single_detector_hdf_files(workflow, hdfbank, insps, output_dir, tags=['full_data']) +check_stop('hdf_trigger_merge', container, workflow, finalize_workflow) + # setup sngl trigger distribution fitting jobs # 'statfiles' is list of files used in calculating statistic # 'dqfiles' is the subset of files containing data quality information @@ -264,6 +355,9 @@ for ifo in ifo_ids.keys(): workflow, hdfbank, inspcomb, statfiles, final_veto_file, final_veto_name, output_dir, tags=ctagsngl) +#check to see if workflow should stop at statmap jobs. +check_stop('statmap', container, workflow, finalize_workflow) + ifo_sets = list(ifo_combos(ifo_ids.keys())) if analyze_singles: ifo_sets += [(ifo,) for ifo in ifo_ids.keys()] @@ -561,19 +655,25 @@ for key in final_bg_files: # DQ log likelihood plots for dqf, dql in zip(dqfiles, dqfile_labels): + ifo = dqf.ifo + outdir = rdir[f'data_quality/{dqf.ifo}_DQ_results'] + # 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 -ifo_bins = workflow.cp.get_subsections('bin_templates') -for ifo in ifo_bins: + # make table of dq segment info + wf.make_dq_segment_table(workflow, dqf, outdir, tags=[dql]) + + # plot background bins background_bins = workflow.cp.get_opt_tags( 'bin_templates', 'background-bins', tags=[ifo]) - wf.make_template_plot(workflow, hdfbank, rdir['data_quality'], + wf.make_template_plot(workflow, hdfbank, outdir, bins=background_bins, tags=[ifo, 'dq_bins'] + bank_tags) + # make table of template bin info + wf.make_template_bin_table(workflow, dqf, outdir, tags=[ifo]) + ############################## Setup the injection runs ####################### splitbank_files_inj = wf.setup_splittable_workflow(workflow, [hdfbank], @@ -830,72 +930,4 @@ wf.make_versioning_page( ) ############################ Finalization #################################### - -# Create the final log file -log_file_html = wf.File(workflow.ifos, 'WORKFLOW-LOG', workflow.analysis_time, - extension='.html', directory=rdir['workflow']) - -gen_file_html = wf.File(workflow.ifos, 'WORKFLOW-GEN', workflow.analysis_time, - extension='.html', directory=rdir['workflow']) - -# Create a page to contain a dashboard link -dashboard_file = wf.File(workflow.ifos, 'DASHBOARD', workflow.analysis_time, - extension='.html', directory=rdir['workflow']) -dashboard_str = """

Pegasus Dashboard Page

""" -kwds = { 'title' : "Pegasus Dashboard", - 'caption' : "Link to Pegasus Dashboard", - 'cmd' : "PYCBC_SUBMIT_DAX_ARGV", } -save_fig_with_metadata(dashboard_str, dashboard_file.storage_path, **kwds) - -# Create pages for the submission script to write data -wf.makedir(rdir['workflow/dax']) -wf.makedir(rdir['workflow/input_map']) -wf.makedir(rdir['workflow/output_map']) -wf.makedir(rdir['workflow/planning']) - -wf.make_results_web_page(finalize_workflow, os.path.join(os.getcwd(), - rdir.base)) - -container += workflow -container += finalize_workflow - -container.add_subworkflow_dependancy(workflow, finalize_workflow) - -container.save() - -# Protect the open box results folder -os.chmod(rdir['open_box_result'], 0o0700) - -logging.info("Written dax.") - -# Close the log and flush to the html file -logging.shutdown() -with open (wf_log_file.storage_path, "r") as logfile: - logdata=logfile.read() -log_str = """ -

Workflow generation script created workflow in output directory: %s

-

Workflow name is: %s

-

Workflow generation script run on host: %s

-
%s
-""" % (os.getcwd(), args.workflow_name, socket.gethostname(), logdata) -kwds = { 'title' : 'Workflow Generation Log', - 'caption' : "Log of the workflow script %s" % sys.argv[0], - 'cmd' :' '.join(sys.argv), } -save_fig_with_metadata(log_str, log_file_html.storage_path, **kwds) - -# Add the command line used to a specific file -args_to_output = [sys.argv[0]] -for arg in sys.argv[1:]: - if arg.startswith('--'): - # This is an option, add tab - args_to_output.append(' ' + arg) - else: - # This is a parameter, add two tabs - args_to_output.append(' ' + arg) - -gen_str = '
' + ' \\\n'.join(args_to_output) + '
' -kwds = { 'title' : 'Workflow Generation Command', - 'caption' : "Command used to generate the workflow.", - 'cmd' :' '.join(sys.argv), } -save_fig_with_metadata(gen_str, gen_file_html.storage_path, **kwds) -layout.single_layout(rdir['workflow'], ([dashboard_file, gen_file_html, log_file_html])) +finalize(container, workflow, finalize_workflow) diff --git a/companion.txt b/companion.txt index 19ea6b0bb73..d49f4121e9e 100644 --- a/companion.txt +++ b/companion.txt @@ -18,6 +18,7 @@ https://github.com/willvousden/ptemcee/archive/master.tar.gz --extra-index-url https://download.pytorch.org/whl/cpu torch nessai>=0.11.0 +snowline # useful to look at PyCBC Live with htop setproctitle diff --git a/docs/frame.rst b/docs/frame.rst index 974f3c81261..0966c1abafc 100644 --- a/docs/frame.rst +++ b/docs/frame.rst @@ -33,7 +33,7 @@ Reading a frame file The ``pycbc.frame`` module provides methods for reading these files into ``TimeSeries`` objects as follows:: >>> from pycbc import frame - >>> data = frame.read_frame('G-G1_RDS_C01_L3-1049587200-60.gwf', 'G1:DER_DATA_H') + >>> data = frame.read_frame('G-G1_RDS_C01_L3-1049587200-60.gwf', 'G1:DER_DATA_H', 1049587200, 1049587200 + 60) Here the first argument is the path to the frame file of interest, while the second lists the `data channel` of interest whose data exist within the file. diff --git a/docs/inference/examples/sampler_platter.rst b/docs/inference/examples/sampler_platter.rst index 963232188f1..35f4158da99 100644 --- a/docs/inference/examples/sampler_platter.rst +++ b/docs/inference/examples/sampler_platter.rst @@ -79,6 +79,21 @@ or an external package (multinest). .. literalinclude:: ../../../examples/inference/samplers/multinest_stub.ini :language: ini +============================================================ +`Snowline `_ +============================================================ + +.. literalinclude:: ../../../examples/inference/samplers/snowline_stub.ini + :language: ini + + +============================================================ +`nessai `_ +============================================================ + +.. literalinclude:: ../../../examples/inference/samplers/nessai_stub.ini + :language: ini + If we run these samplers, we create the following plot: .. image:: ../../_include/sample.png diff --git a/examples/inference/samplers/run.sh b/examples/inference/samplers/run.sh index 91e41e7fbb7..bdeaefc8d6a 100755 --- a/examples/inference/samplers/run.sh +++ b/examples/inference/samplers/run.sh @@ -1,5 +1,5 @@ #!/bin/sh -for f in cpnest_stub.ini emcee_stub.ini emcee_pt_stub.ini dynesty_stub.ini ultranest_stub.ini epsie_stub.ini nessai_stub.ini; do +for f in cpnest_stub.ini emcee_stub.ini emcee_pt_stub.ini dynesty_stub.ini ultranest_stub.ini epsie_stub.ini nessai_stub.ini snowline_stub.ini; do echo $f pycbc_inference \ --config-files `dirname $0`/simp.ini `dirname $0`/$f \ @@ -17,6 +17,7 @@ ultranest_stub.ini.hdf:ultranest \ epsie_stub.ini.hdf:espie \ cpnest_stub.ini.hdf:cpnest \ nessai_stub.ini.hdf:nessai \ +snowline_stub.ini.hdf:snowline \ --output-file sample.png \ --plot-contours \ --plot-marginal \ diff --git a/examples/inference/samplers/snowline_stub.ini b/examples/inference/samplers/snowline_stub.ini new file mode 100644 index 00000000000..ac17879fe74 --- /dev/null +++ b/examples/inference/samplers/snowline_stub.ini @@ -0,0 +1,6 @@ +[sampler] +name = snowline + +num_gauss_samples = 500 +min_ess = 500 +max_improvement_loops = 3 diff --git a/examples/live/run.sh b/examples/live/run.sh index 82b79db5a24..8857c5d996e 100755 --- a/examples/live/run.sh +++ b/examples/live/run.sh @@ -195,10 +195,11 @@ python -m mpi4py `which pycbc_live` \ --src-class-eff-to-lum-distance 0.74899 \ --src-class-lum-distance-to-delta -0.51557 -0.32195 \ --run-snr-optimization \ ---snr-opt-di-maxiter 50 \ ---snr-opt-di-popsize 100 \ ---snr-opt-include-candidate \ ---snr-opt-seed 42 \ +--snr-opt-extra-opts \ + "--snr-opt-method differential_evolution \ + --snr-opt-di-maxiter 50 \ + --snr-opt-di-popsize 100 \ + --snr-opt-include-candidate " \ --sngl-ifar-est-dist conservative \ --single-newsnr-threshold 9 \ --single-duration-threshold 7 \ @@ -210,11 +211,14 @@ python -m mpi4py `which pycbc_live` \ # If you would like to use the pso optimizer, change --optimizer to pso # and include these arguments while removing other optimizer args. # You will need to install the pyswarms package into your environment. -# --snr-opt-pso-iters 5 \ -# --snr-opt-pso-particles 250 \ -# --snr-opt-pso-c1 0.5 \ -# --snr-opt-pso-c2 2.0 \ -# --snr-opt-pso-w 0.01 \ +# --snr-opt-extra-opts \ +# "--snr-opt-method pso \ +# --snr-opt-pso-iters 5 \ +# --snr-opt-pso-particles 250 \ +# --snr-opt-pso-c1 0.5 \ +# --snr-opt-pso-c2 2.0 \ +# --snr-opt-pso-w 0.01 \ +# --snr-opt-include-candidate " \ # note that, at this point, some SNR optimization processes may still be # running, so the checks below may ignore their results diff --git a/pycbc/__init__.py b/pycbc/__init__.py index 97acf32e7a9..e716caf4a82 100644 --- a/pycbc/__init__.py +++ b/pycbc/__init__.py @@ -79,11 +79,11 @@ def add_common_pycbc_options(parser): title="PyCBC common options", description="Common options for PyCBC executables.", ) - group.add_argument('-V', '--verbose', action='count', default=0, + 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 ' + '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, diff --git a/pycbc/_version.py b/pycbc/_version.py index 31ebf7f4fea..617cf4a523e 100644 --- a/pycbc/_version.py +++ b/pycbc/_version.py @@ -13,27 +13,40 @@ # 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. + """ This modules contains a function to provide an argparse action that reports extremely verbose version information for PyCBC, lal, and lalsimulation. """ -import os, sys +import os +import sys +import glob import argparse import inspect import subprocess + def print_link(library): - err_msg = "Could not execute runtime linker to determine\n" + \ - "shared library paths for library:\n " + library + "\n" - FNULL = open(os.devnull, 'w') + err_msg = ( + "Could not execute runtime linker to determine\n" + f"shared library paths for library:\n {library}\n" + ) try: - link = subprocess.check_output(['ldd', library], - stderr=FNULL) + # Linux + link = subprocess.check_output( + ['ldd', library], + stderr=subprocess.DEVNULL, + text=True + ) except OSError: try: - link = subprocess.check_output(['otool', '-L', library], - stderr=FNULL) + # macOS + link = subprocess.check_output( + ['otool', '-L', library], + stderr=subprocess.DEVNULL, + text=True + ) except: link = err_msg except: @@ -41,42 +54,58 @@ def print_link(library): return link +def get_lal_info(module, lib_glob): + """Return a string reporting the version and runtime library information + for a LAL Python import. + """ + module_path = inspect.getfile(module) + version_str = ( + module.git_version.verbose_msg + + "\n\nImported from: " + module_path + + "\n\nRuntime libraries:\n" + ) + possible_lib_paths = glob.glob( + os.path.join(os.path.dirname(module_path), lib_glob) + ) + for lib_path in possible_lib_paths: + version_str += print_link(lib_path) + return version_str + + class Version(argparse.Action): - """ print the pycbc, lal and lalsimulation versions """ + """Subclass of argparse.Action that prints version information for PyCBC, + LAL and LALSimulation. + """ def __init__(self, nargs=0, **kw): super(Version, self).__init__(nargs=nargs, **kw) - def __call__(self, parser, namespace, values, option_string=None): - import pycbc - version_str="--- PyCBC Version --------------------------\n" + \ - pycbc.version.git_verbose_msg + \ + + version_str = ( + "--- PyCBC Version --------------------------\n" + + pycbc.version.git_verbose_msg + "\n\nImported from: " + inspect.getfile(pycbc) + ) version_str += "\n\n--- LAL Version ----------------------------\n" try: import lal.git_version - lal_module = inspect.getfile(lal) - lal_library = os.path.join( os.path.dirname(lal_module), - '_lal.so') - version_str += lal.git_version.verbose_msg + \ - "\n\nImported from: " + lal_module + \ - "\n\nRuntime libraries:\n" + print_link(lal_library) except ImportError: version_str += "\nLAL not installed in environment\n" + else: + version_str += get_lal_info(lal, '_lal*.so') version_str += "\n\n--- LALSimulation Version-------------------\n" try: import lalsimulation.git_version - lalsimulation_module = inspect.getfile(lalsimulation) - lalsimulation_library = os.path.join( os.path.dirname(lalsimulation_module), - '_lalsimulation.so') - version_str += lalsimulation.git_version.verbose_msg + \ - "\n\nImported from: " + lalsimulation_module + \ - "\n\nRuntime libraries:\n" + print_link(lalsimulation_library) except ImportError: version_str += "\nLALSimulation not installed in environment\n" + else: + version_str += get_lal_info(lalsimulation, '_lalsimulation*.so') print(version_str) sys.exit(0) + + +__all__ = ['Version'] diff --git a/pycbc/_version_helper.py b/pycbc/_version_helper.py index 694895a4b93..4dde5936e5c 100644 --- a/pycbc/_version_helper.py +++ b/pycbc/_version_helper.py @@ -64,30 +64,29 @@ def call(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, if p.returncode != 0 and on_error == 'raise': raise GitInvocationError('Failed to run "%s"' % " ".join(command)) - out = out.decode('utf-8') + out = out.decode('utf-8').strip() if returncode: - return out.strip(), p.returncode - else: - return out.strip() + return out, p.returncode + return out def get_build_name(git_path='git'): """Find the username of the current builder """ - name,retcode = call(('git', 'config', 'user.name'), returncode=True) + name, retcode = call(('git', 'config', 'user.name'), returncode=True) if retcode: name = "Unknown User" - email,retcode = call(('git', 'config', 'user.email'), returncode=True) + email, retcode = call(('git', 'config', 'user.email'), returncode=True) if retcode: email = "" - return "%s <%s>" % (name, email) + return f"{name} <{email}>" def get_build_date(): """Returns the current datetime as the git build date """ - return time.strftime('%Y-%m-%d %H:%M:%S +0000', time.gmtime()) + return time.strftime(r'%Y-%m-%d %H:%M:%S +0000', time.gmtime()) def get_last_commit(git_path='git'): @@ -96,12 +95,15 @@ def get_last_commit(git_path='git'): Returns a tuple (hash, date, author name, author e-mail, committer name, committer e-mail). """ - hash_, udate, aname, amail, cname, cmail = ( - call((git_path, 'log', '-1', - '--pretty=format:%H,%ct,%an,%ae,%cn,%ce')).split(",")) - date = time.strftime('%Y-%m-%d %H:%M:%S +0000', time.gmtime(float(udate))) - author = '%s <%s>' % (aname, amail) - committer = '%s <%s>' % (cname, cmail) + hash_, udate, aname, amail, cname, cmail = call(( + git_path, + 'log', + '-1', + r'--pretty=format:%H,%ct,%an,%ae,%cn,%ce' + )).split(",") + date = time.strftime(r'%Y-%m-%d %H:%M:%S +0000', time.gmtime(float(udate))) + author = f'{aname} <{amail}>' + committer = f'{cname} <{cmail}>' return hash_, date, author, committer @@ -111,8 +113,7 @@ def get_git_branch(git_path='git'): branch_match = call((git_path, 'rev-parse', '--symbolic-full-name', 'HEAD')) if branch_match == "HEAD": return None - else: - return os.path.basename(branch_match) + return os.path.basename(branch_match) def get_git_tag(hash_, git_path='git'): @@ -122,26 +123,26 @@ def get_git_tag(hash_, git_path='git'): '--tags', hash_), returncode=True) if status == 0: return tag - else: - return None + return None + def get_num_commits(): return call(('git', 'rev-list', '--count', 'HEAD')) + def get_git_status(git_path='git'): """Returns the state of the git working copy """ status_output = subprocess.call((git_path, 'diff-files', '--quiet')) if status_output != 0: return 'UNCLEAN: Modified working tree' - else: - # check index for changes - status_output = subprocess.call((git_path, 'diff-index', '--cached', - '--quiet', 'HEAD')) - if status_output != 0: - return 'UNCLEAN: Modified index' - else: - return 'CLEAN: All modifications committed' + # check index for changes + status_output = subprocess.call((git_path, 'diff-index', '--cached', + '--quiet', 'HEAD')) + if status_output != 0: + return 'UNCLEAN: Modified index' + return 'CLEAN: All modifications committed' + def determine_latest_release_version(): """Query the git repository for the last released version of the code. @@ -152,22 +153,21 @@ def determine_latest_release_version(): tag_list = call((git_path, 'tag')).split('\n') # Reduce to only versions - tag_list = [t[1:] for t in tag_list if t.startswith('v')] + re_magic = re.compile(r"v\d+\.\d+\.\d+$") + tag_list = [t[1:] for t in tag_list if re_magic.match(t)] - # Determine if indeed a tag and store largest + # find latest version latest_version = None latest_version_string = None - re_magic = re.compile("\d+\.\d+\.\d+$") for tag in tag_list: - # Is this a version string - if re_magic.match(tag): - curr_version = distutils.version.StrictVersion(tag) - if latest_version is None or curr_version > latest_version: - latest_version = curr_version - latest_version_string = tag + curr_version = distutils.version.StrictVersion(tag) + if latest_version is None or curr_version > latest_version: + latest_version = curr_version + latest_version_string = tag return latest_version_string + def generate_git_version_info(): """Query the git repository information to generate a version module. """ diff --git a/pycbc/events/eventmgr.py b/pycbc/events/eventmgr.py index ec065900414..80e4d2c8a64 100644 --- a/pycbc/events/eventmgr.py +++ b/pycbc/events/eventmgr.py @@ -24,10 +24,13 @@ """This modules defines functions for clustering and thresholding timeseries to produces event triggers """ -import numpy, copy, os.path +import os.path +import copy +import itertools import logging -import h5py import pickle +import numpy +import h5py from pycbc.types import Array from pycbc.scheme import schemed @@ -177,6 +180,23 @@ def cluster_reduce(idx, snr, window_size): return idx.take(ind), snr.take(ind) +class H5FileSyntSugar(object): + """Convenience class that adds some syntactic sugar to h5py.File. + """ + def __init__(self, name, prefix=''): + self.f = h5py.File(name, 'w') + self.prefix = prefix + + def __setitem__(self, name, data): + self.f.create_dataset( + self.prefix + '/' + name, + data=data, + compression='gzip', + compression_opts=9, + shuffle=True + ) + + class EventManager(object): def __init__(self, opt, column, column_types, **kwds): self.opt = opt @@ -407,33 +427,21 @@ def save_performance(self, ncores, nfilters, ntemplates, run_time, self.write_performance = True def write_events(self, outname): - """ Write the found events to a sngl inspiral table + """Write the found events to a file. The only currently supported + format is HDF5, indicated by an .hdf or .h5 extension. """ self.make_output_dir(outname) - - if '.hdf' in outname: + if outname.endswith(('.hdf', '.h5')): self.write_to_hdf(outname) - else: - raise ValueError('Cannot write to this format') + return + raise ValueError('Unsupported event output file format') def write_to_hdf(self, outname): - class fw(object): - def __init__(self, name, prefix): - self.f = h5py.File(name, 'w') - self.prefix = prefix - - def __setitem__(self, name, data): - col = self.prefix + '/' + name - self.f.create_dataset(col, data=data, - compression='gzip', - compression_opts=9, - shuffle=True) - self.events.sort(order='template_id') th = numpy.array([p['tmplt'].template_hash for p in self.template_params]) tid = self.events['template_id'] - f = fw(outname, self.opt.channel_name[0:2]) + f = H5FileSyntSugar(outname, self.opt.channel_name[0:2]) if len(self.events): f['snr'] = abs(self.events['snr']) @@ -471,6 +479,10 @@ def __setitem__(self, name, data): # Not precessing f['sigmasq'] = self.events['sigmasq'] + # Template durations should ideally be stored in the bank file. + # At present, however, a few plotting/visualization codes + # downstream in the offline search workflow rely on durations being + # stored in the trigger files instead. template_durations = [p['tmplt'].template_duration for p in self.template_params] f['template_duration'] = numpy.array(template_durations, @@ -595,6 +607,31 @@ def add_template_events_to_ifo(self, ifo, columns, vectors): self.template_event_dict[ifo] = self.template_events self.template_events = None + def write_gating_info_to_hdf(self, hf): + """Write per-detector gating information to an h5py file object. + The information is laid out according to the following groups and + datasets: `//gating/{file, auto}/{time, width, pad}` where + "file" and "auto" indicate respectively externally-provided gates and + internally-generated gates (autogating), and "time", "width" and "pad" + indicate the gate center times, total durations and padding durations + in seconds respectively. + """ + if 'gating_info' not in self.global_params: + return + gates = self.global_params['gating_info'] + for ifo, gate_type in itertools.product(self.ifos, ['file', 'auto']): + if gate_type not in gates[ifo]: + continue + hf[f'{ifo}/gating/{gate_type}/time'] = numpy.array( + [float(g[0]) for g in gates[ifo][gate_type]] + ) + hf[f'{ifo}/gating/{gate_type}/width'] = numpy.array( + [g[1] for g in gates[ifo][gate_type]] + ) + hf[f'{ifo}/gating/{gate_type}/pad'] = numpy.array( + [g[2] for g in gates[ifo][gate_type]] + ) + class EventManagerCoherent(EventManagerMultiDetBase): def __init__(self, opt, ifos, column, column_types, network_column, @@ -656,7 +693,6 @@ def add_template_network_events(self, columns, vectors): """ Add a vector indexed """ # initialize with zeros - since vectors can be None, look for the # longest one that isn't - new_events = None new_events = numpy.zeros( max([len(v) for v in vectors if v is not None]), dtype=self.network_event_dtype @@ -681,19 +717,11 @@ def add_template_events_to_network(self, columns, vectors): self.template_events = None def write_to_hdf(self, outname): - class fw(object): - def __init__(self, name): - self.f = h5py.File(name, 'w') - - def __setitem__(self, name, data): - col = self.prefix + '/' + name - self.f.create_dataset( - col, data=data, compression='gzip', compression_opts=9, - shuffle=True) self.events.sort(order='template_id') th = numpy.array( [p['tmplt'].template_hash for p in self.template_params]) - f = fw(outname) + f = H5FileSyntSugar(outname) + self.write_gating_info_to_hdf(f) # Output network stuff f.prefix = 'network' network_events = numpy.array( @@ -721,8 +749,7 @@ def __setitem__(self, name, data): ifo_events = numpy.array([e for e in self.events if e['ifo'] == self.ifo_dict[ifo]], dtype=self.event_dtype) if len(ifo_events): - ifo_str = ifo.lower()[0] if ifo != 'H1' else ifo.lower() - f['snr_%s' % ifo_str] = abs(ifo_events['snr']) + f['snr'] = abs(ifo_events['snr']) f['event_id'] = ifo_events['event_id'] try: # Precessing @@ -736,8 +763,8 @@ def __setitem__(self, name, data): f['bank_chisq_dof'] = ifo_events['bank_chisq_dof'] f['cont_chisq'] = ifo_events['cont_chisq'] f['end_time'] = ifo_events['time_index'] / \ - float(self.opt.sample_rate[ifo_str]) + \ - self.opt.gps_start_time[ifo_str] + float(self.opt.sample_rate[ifo]) + \ + self.opt.gps_start_time[ifo] f['time_index'] = ifo_events['time_index'] f['slide_id'] = ifo_events['slide_id'] try: @@ -764,11 +791,6 @@ def __setitem__(self, name, data): dtype=numpy.float32) f['sigmasq'] = template_sigmasq[tid] - template_durations = [p['tmplt'].template_duration for p in - self.template_params] - f['template_duration'] = numpy.array(template_durations, - dtype=numpy.float32)[tid] - # FIXME: Can we get this value from the autochisq instance? # cont_dof = self.opt.autochi_number_points # if self.opt.autochi_onesided is None: @@ -820,17 +842,6 @@ def __setitem__(self, name, data): f['search/setup_time_fraction'] = \ numpy.array([float(self.setup_time) / float(self.run_time)]) - if 'gating_info' in self.global_params: - gating_info = self.global_params['gating_info'] - for gate_type in ['file', 'auto']: - if gate_type in gating_info: - f['gating/' + gate_type + '/time'] = numpy.array( - [float(g[0]) for g in gating_info[gate_type]]) - f['gating/' + gate_type + '/width'] = numpy.array( - [g[1] for g in gating_info[gate_type]]) - f['gating/' + gate_type + '/pad'] = numpy.array( - [g[2] for g in gating_info[gate_type]]) - def finalize_template_events(self): # Check that none of the template events have the same time index as an # existing event in events. I.e. don't list the same ifo event multiple @@ -965,41 +976,20 @@ def finalize_template_events(self, perform_coincidence=True, self.template_event_dict[ifo] = numpy.array([], dtype=self.event_dtype) - def write_events(self, outname): - """ Write the found events to a sngl inspiral table - """ - self.make_output_dir(outname) - - if '.hdf' in outname: - self.write_to_hdf(outname) - else: - raise ValueError('Cannot write to this format') - def write_to_hdf(self, outname): - class fw(object): - def __init__(self, name): - self.f = h5py.File(name, 'w') - - def __setitem__(self, name, data): - col = self.prefix + '/' + name - self.f.create_dataset(col, data=data, - compression='gzip', - compression_opts=9, - shuffle=True) - self.events.sort(order='template_id') th = numpy.array([p['tmplt'].template_hash for p in self.template_params]) tid = self.events['template_id'] - f = fw(outname) + f = H5FileSyntSugar(outname) + self.write_gating_info_to_hdf(f) for ifo in self.ifos: f.prefix = ifo ifo_events = numpy.array([e for e in self.events if e['ifo'] == self.ifo_dict[ifo]], dtype=self.event_dtype) if len(ifo_events): - ifo_str = ifo.lower()[0] if ifo != 'H1' else ifo.lower() - f['snr_%s' % ifo_str] = abs(ifo_events['snr']) + f['snr'] = abs(ifo_events['snr']) try: # Precessing f['u_vals'] = ifo_events['u_vals'] @@ -1012,8 +1002,8 @@ def __setitem__(self, name, data): f['bank_chisq_dof'] = ifo_events['bank_chisq_dof'] f['cont_chisq'] = ifo_events['cont_chisq'] f['end_time'] = ifo_events['time_index'] / \ - float(self.opt.sample_rate[ifo_str]) + \ - self.opt.gps_start_time[ifo_str] + float(self.opt.sample_rate[ifo]) + \ + self.opt.gps_start_time[ifo] try: # Precessing template_sigmasq_plus = numpy.array([t['sigmasq_plus'] for @@ -1034,11 +1024,6 @@ def __setitem__(self, name, data): dtype=numpy.float32) f['sigmasq'] = template_sigmasq[tid] - template_durations = [p['tmplt'].template_duration for p in - self.template_params] - f['template_duration'] = \ - numpy.array(template_durations, dtype=numpy.float32)[tid] - # FIXME: Can we get this value from the autochisq instance? cont_dof = self.opt.autochi_number_points if self.opt.autochi_onesided is None: @@ -1093,17 +1078,6 @@ def __setitem__(self, name, data): f['search/setup_time_fraction'] = \ numpy.array([float(self.setup_time) / float(self.run_time)]) - if 'gating_info' in self.global_params: - gating_info = self.global_params['gating_info'] - for gate_type in ['file', 'auto']: - if gate_type in gating_info: - f['gating/' + gate_type + '/time'] = numpy.array( - [float(g[0]) for g in gating_info[gate_type]]) - f['gating/' + gate_type + '/width'] = numpy.array( - [g[1] for g in gating_info[gate_type]]) - f['gating/' + gate_type + '/pad'] = numpy.array( - [g[2] for g in gating_info[gate_type]]) - __all__ = ['threshold_and_cluster', 'findchirp_cluster_over_window', 'threshold', 'cluster_reduce', 'ThresholdCluster', diff --git a/pycbc/inference/io/__init__.py b/pycbc/inference/io/__init__.py index 4b3fd0ce909..3157f4408f6 100644 --- a/pycbc/inference/io/__init__.py +++ b/pycbc/inference/io/__init__.py @@ -37,6 +37,7 @@ from .multinest import MultinestFile from .dynesty import DynestyFile from .ultranest import UltranestFile +from .snowline import SnowlineFile from .nessai import NessaiFile from .posterior import PosteriorFile from .txt import InferenceTXTFile @@ -51,6 +52,7 @@ PosteriorFile.name: PosteriorFile, UltranestFile.name: UltranestFile, NessaiFile.name: NessaiFile, + SnowlineFile.name: SnowlineFile, } try: diff --git a/pycbc/inference/io/snowline.py b/pycbc/inference/io/snowline.py new file mode 100644 index 00000000000..4aa9edcbf13 --- /dev/null +++ b/pycbc/inference/io/snowline.py @@ -0,0 +1,32 @@ +# Copyright (C) 2024 Alex Nitz +# 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 +# self.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 +# +# ============================================================================= +# +"""Provides IO for the snowline sampler. +""" +from .posterior import PosteriorFile + + +class SnowlineFile(PosteriorFile): + """Class to handle file IO for the ``snowline`` sampler.""" + + name = 'snowline_file' diff --git a/pycbc/inference/sampler/__init__.py b/pycbc/inference/sampler/__init__.py index 41da16f39c6..a25573b7ee2 100644 --- a/pycbc/inference/sampler/__init__.py +++ b/pycbc/inference/sampler/__init__.py @@ -25,6 +25,7 @@ from .ultranest import UltranestSampler from .dummy import DummySampler from .refine import RefineSampler +from .snowline import SnowlineSampler # list of available samplers samplers = {cls.name: cls for cls in ( @@ -32,6 +33,7 @@ UltranestSampler, DummySampler, RefineSampler, + SnowlineSampler, )} try: diff --git a/pycbc/inference/sampler/snowline.py b/pycbc/inference/sampler/snowline.py new file mode 100644 index 00000000000..6a54825879a --- /dev/null +++ b/pycbc/inference/sampler/snowline.py @@ -0,0 +1,168 @@ +# Copyright (C) 2023 Alex Nitz +# 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 +# +# ============================================================================= +# +""" +This modules provides classes and functions for using the snowline sampler +packages for parameter estimation. +""" + +import sys +import logging + +from pycbc.inference.io.snowline import SnowlineFile +from pycbc.io.hdf import dump_state +from pycbc.pool import use_mpi +from .base import (BaseSampler, setup_output) +from .base_cube import setup_calls + + +# +# ============================================================================= +# +# Samplers +# +# ============================================================================= +# + +class SnowlineSampler(BaseSampler): + """This class is used to construct an Snowline sampler from the snowline + package. + + Parameters + ---------- + model : model + A model from ``pycbc.inference.models`` + """ + name = "snowline" + _io = SnowlineFile + + def __init__(self, model, **kwargs): + super().__init__(model) + + import snowline + log_likelihood_call, prior_call = setup_calls(model, copy_prior=True) + + self._sampler = snowline.ReactiveImportanceSampler( + list(self.model.variable_params), + log_likelihood_call, + transform=prior_call) + + do_mpi, _, rank = use_mpi() + self.main = (not do_mpi) or (rank == 0) + self.nlive = 0 + self.ndim = len(self.model.variable_params) + self.kwargs = kwargs + + def run(self): + self.result = self._sampler.run(**self.kwargs) + if not self.main: + sys.exit(0) + self._sampler.print_results() + + @property + def io(self): + return self._io + + @property + def niterations(self): + return self.result['niter'] + + @classmethod + def from_config(cls, cp, model, output_file=None, **kwds): + """ + Loads the sampler from the given config file. + """ + skeys = {} + opts = {'num_global_samples': int, + 'num_gauss_samples': int, + 'max_ncalls': int, + 'min_ess': int, + 'max_improvement_loops': int + } + for opt_name in opts: + if cp.has_option('sampler', opt_name): + value = cp.get('sampler', opt_name) + skeys[opt_name] = opts[opt_name](value) + inst = cls(model, **skeys) + + do_mpi, _, rank = use_mpi() + if not do_mpi or (rank == 0): + setup_output(inst, output_file) + return inst + + def checkpoint(self): + """ There is currently no checkpointing implemented""" + pass + + def resume_from_checkpoint(self): + """ There is currently no checkpointing implemented""" + pass + + def finalize(self): + logging.info("Writing samples to files") + for fn in [self.checkpoint_file, self.backup_file]: + self.write_results(fn) + + @property + def model_stats(self): + return {} + + @property + def samples(self): + samples = self.result['samples'] + params = list(self.model.variable_params) + samples_dict = {p: samples[:, i] for i, p in enumerate(params)} + return samples_dict + + def write_results(self, filename): + """Writes samples, model stats, acceptance fraction, and random state + to the given file. + + Parameters + ---------- + filename : str + The file to write to. The file is opened using the ``io`` class + in an an append state. + """ + with self.io(filename, 'a') as fp: + # write samples + fp.write_samples(self.samples, self.samples.keys()) + # write log evidence + fp.write_logevidence(self.logz, self.logz_err) + + # write full results + dump_state(self.result, fp, + path='sampler_info', + dsetname='presult') + + @property + def logz(self): + """Return bayesian evidence estimated by snowline sampler. + """ + return self.result['logz'] + + @property + def logz_err(self): + """Return error in bayesian evidence estimated by snowline sampler. + """ + return self.result['logzerr'] diff --git a/pycbc/live/snr_optimizer.py b/pycbc/live/snr_optimizer.py index 9e6b80481a3..34ab4cd2de2 100644 --- a/pycbc/live/snr_optimizer.py +++ b/pycbc/live/snr_optimizer.py @@ -328,13 +328,13 @@ def insert_snr_optimizer_options(parser): opt_opt_group.add_argument('--snr-opt-include-candidate', action='store_true', help='Include parameters of the candidate event in the initialized ' - 'array for the optimizer. Only relevant for --optimizer pso or ' - 'differential_evolution') + 'array for the optimizer. Only relevant for --snr-opt-method pso ' + 'or differential_evolution') opt_opt_group.add_argument('--snr-opt-seed', default='42', help='Seed to supply to the random generation of initial array to ' - 'pass to the optimizer. Only relevant for --optimizer pso or ' - 'differential_evolution. Set to ''random'' for a random seed') + 'pass to the optimizer. Only relevant for --snr-opt-method pso ' + 'or differential_evolution. Set to ''random'' for a random seed') # For each optimizer, add the possible options for optimizer, option_subdict in option_dict.items(): @@ -345,7 +345,7 @@ def insert_snr_optimizer_options(parser): option_name = f"--snr-opt-{optimizer_name}-{opt_name}" opt_opt_group.add_argument(option_name, type=float, - help=f'Only relevant for --optimizer {optimizer}: ' + + help=f'Only relevant for --snr-opt-method {optimizer}: ' + opt_help_default[0] + f' Default = {opt_help_default[1]}') @@ -381,24 +381,3 @@ def check_snr_optimizer_options(args, parser): key_name = f'snr_opt_{optimizer_name}_{key}' if not getattr(args, key_name): setattr(args, key_name, value[1]) - - -def args_to_string(args): - """ - Convert the supplied arguments for SNR optimization config into - a string - this is to be used when running subprocesses - """ - argstr = f'--snr-opt-method {args.snr_opt_method} ' - optimizer_name = args.snr_opt_method.replace('_', '-') - if optimizer_name == 'differential-evolution': - optimizer_name = 'di' - for opt in option_dict[args.snr_opt_method]: - key_name = f'snr_opt_{optimizer_name}_{opt}' - option_value = getattr(args, key_name) - # If the option is not given, don't pass it and use default - if option_value is None: - continue - option_fullname = f'--snr-opt-{optimizer_name}-{opt}' - argstr += f'{option_fullname} {option_value} ' - - return argstr diff --git a/pycbc/pool.py b/pycbc/pool.py index 9750b4f4c91..8f2095181ee 100644 --- a/pycbc/pool.py +++ b/pycbc/pool.py @@ -121,6 +121,12 @@ def broadcast(self, fcn, args): def map(self, f, items): return [f(a) for a in items] + # This is single core, so imap and map + # would not behave differently. This is defined + # so that the general pool interfaces can use + # imap irrespective of the pool type. + imap = map + def use_mpi(require_mpi=False, log=True): """ Get whether MPI is enabled and if so the current size and rank """ diff --git a/pycbc/psd/analytical_space.py b/pycbc/psd/analytical_space.py index 1f388ae59bd..4237a724124 100644 --- a/pycbc/psd/analytical_space.py +++ b/pycbc/psd/analytical_space.py @@ -29,7 +29,7 @@ """ import numpy as np -from astropy import constants +from astropy.constants import c from pycbc.psd.read import from_numpy_arrays @@ -49,11 +49,11 @@ def psd_lisa_acc_noise(f, acc_noise_level=3e-15): The PSD value or array for acceleration noise. Notes ----- - Pease see Eq.(11-13) in for more details. + Please see Eq.(11-13) in for more details. """ s_acc = acc_noise_level**2 * (1+(4e-4/f)**2)*(1+(f/8e-3)**4) s_acc_d = s_acc * (2*np.pi*f)**(-4) - s_acc_nu = (2*np.pi*f/constants.c.value)**2 * s_acc_d + s_acc_nu = (2*np.pi*f/c.value)**2 * s_acc_d return s_acc_nu @@ -74,10 +74,10 @@ def psd_lisa_oms_noise(f, oms_noise_level=15e-12): The PSD value or array for OMS noise. Notes ----- - Pease see Eq.(9-10) in for more details. + Please see Eq.(9-10) in for more details. """ s_oms_d = oms_noise_level**2 * (1+(2e-3/f)**4) - s_oms_nu = s_oms_d * (2*np.pi*f/constants.c.value)**2 + s_oms_nu = s_oms_d * (2*np.pi*f/c.value)**2 return s_oms_nu @@ -96,13 +96,13 @@ def lisa_psd_components(f, acc_noise_level=3e-15, oms_noise_level=15e-12): Returns ------- - [low_freq_component, high_freq_component] : list + low_freq_component, high_freq_component : The PSD value or array for acceleration and OMS noise. """ low_freq_component = psd_lisa_acc_noise(f, acc_noise_level) high_freq_component = psd_lisa_oms_noise(f, oms_noise_level) - return [low_freq_component, high_freq_component] + return low_freq_component, high_freq_component def omega_length(f, len_arm=2.5e9): @@ -120,7 +120,7 @@ def omega_length(f, len_arm=2.5e9): omega_len : float or numpy.array The value of 2*pi*f*LISA_arm_length. """ - omega_len = 2*np.pi*f * len_arm/constants.c.value + omega_len = 2*np.pi*f * len_arm/c.value return omega_len @@ -151,23 +151,18 @@ def analytical_psd_lisa_tdi_1p5_XYZ(length, delta_f, low_freq_cutoff, The TDI-1.5 PSD (X,Y,Z channel) for LISA. Notes ----- - Pease see Eq.(19) in for more details. + Please see Eq.(19) in for more details. """ len_arm = np.float64(len_arm) acc_noise_level = np.float64(acc_noise_level) oms_noise_level = np.float64(oms_noise_level) - psd = [] fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length) - for f in fr: - [s_acc_nu, s_oms_nu] = lisa_psd_components( - f, acc_noise_level, oms_noise_level) - omega_len = omega_length(f, len_arm) - psd.append(16*(np.sin(omega_len))**2 * - (s_oms_nu+s_acc_nu*(3+np.cos(omega_len)))) - fseries = from_numpy_arrays(fr, np.array(psd), - length, delta_f, low_freq_cutoff) + s_acc_nu, s_oms_nu = lisa_psd_components( + fr, acc_noise_level, oms_noise_level) + omega_len = omega_length(fr, len_arm) + psd = 16*(np.sin(omega_len))**2 * (s_oms_nu+s_acc_nu*(3+np.cos(omega_len))) - return fseries + return from_numpy_arrays(fr, psd, length, delta_f, low_freq_cutoff) def analytical_psd_lisa_tdi_2p0_XYZ(length, delta_f, low_freq_cutoff, @@ -196,23 +191,19 @@ def analytical_psd_lisa_tdi_2p0_XYZ(length, delta_f, low_freq_cutoff, The TDI-2.0 PSD (X,Y,Z channel) for LISA. Notes ----- - Pease see Eq.(20) in for more details. + Please see Eq.(20) in for more details. """ len_arm = np.float64(len_arm) acc_noise_level = np.float64(acc_noise_level) oms_noise_level = np.float64(oms_noise_level) - psd = [] fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length) - for f in fr: - [s_acc_nu, s_oms_nu] = lisa_psd_components( - f, acc_noise_level, oms_noise_level) - omega_len = omega_length(f, len_arm) - psd.append(64*(np.sin(omega_len))**2 * (np.sin(2*omega_len))**2 * - (s_oms_nu+s_acc_nu*(3+np.cos(2*omega_len)))) - fseries = from_numpy_arrays(fr, np.array(psd), - length, delta_f, low_freq_cutoff) + s_acc_nu, s_oms_nu = lisa_psd_components( + fr, acc_noise_level, oms_noise_level) + omega_len = omega_length(fr, len_arm) + psd = (64*(np.sin(omega_len))**2 * (np.sin(2*omega_len))**2 * + (s_oms_nu+s_acc_nu*(3+np.cos(2*omega_len)))) - return fseries + return from_numpy_arrays(fr, psd, length, delta_f, low_freq_cutoff) def analytical_csd_lisa_tdi_1p5_XY(length, delta_f, low_freq_cutoff, @@ -241,22 +232,19 @@ def analytical_csd_lisa_tdi_1p5_XY(length, delta_f, low_freq_cutoff, The CSD between LISA's TDI-1.5 channel X and Y. Notes ----- - Pease see Eq.(56) in for more details. + Please see Eq.(56) in for more details. """ len_arm = np.float64(len_arm) acc_noise_level = np.float64(acc_noise_level) oms_noise_level = np.float64(oms_noise_level) - csd = [] fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length) - for f in fr: - omega_len = omega_length(f, len_arm) - [s_acc_nu, s_oms_nu] = lisa_psd_components( - f, acc_noise_level, oms_noise_level) - csd.append(-8*np.sin(omega_len)**2 * np.cos(omega_len) * - (s_oms_nu+4*s_acc_nu)) - fseries = from_numpy_arrays(fr, np.array(csd), - length, delta_f, low_freq_cutoff) - return fseries + s_acc_nu, s_oms_nu = lisa_psd_components( + fr, acc_noise_level, oms_noise_level) + omega_len = omega_length(fr, len_arm) + csd = (-8*np.sin(omega_len)**2 * np.cos(omega_len) * + (s_oms_nu+4*s_acc_nu)) + + return from_numpy_arrays(fr, csd, length, delta_f, low_freq_cutoff) def analytical_psd_lisa_tdi_1p5_AE(length, delta_f, low_freq_cutoff, @@ -285,24 +273,20 @@ def analytical_psd_lisa_tdi_1p5_AE(length, delta_f, low_freq_cutoff, The PSD of LISA's TDI-1.5 channel A and E. Notes ----- - Pease see Eq.(58) in for more details. + Please see Eq.(58) in for more details. """ len_arm = np.float64(len_arm) acc_noise_level = np.float64(acc_noise_level) oms_noise_level = np.float64(oms_noise_level) - psd = [] fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length) - for f in fr: - [s_acc_nu, s_oms_nu] = lisa_psd_components( - f, acc_noise_level, oms_noise_level) - omega_len = omega_length(f, len_arm) - psd.append(8*(np.sin(omega_len))**2 * - (4*(1+np.cos(omega_len)+np.cos(omega_len)**2)*s_acc_nu + - (2+np.cos(omega_len))*s_oms_nu)) - fseries = from_numpy_arrays(fr, np.array(psd), - length, delta_f, low_freq_cutoff) + s_acc_nu, s_oms_nu = lisa_psd_components( + fr, acc_noise_level, oms_noise_level) + omega_len = omega_length(fr, len_arm) + psd = (8*(np.sin(omega_len))**2 * + (4*(1+np.cos(omega_len)+np.cos(omega_len)**2)*s_acc_nu + + (2+np.cos(omega_len))*s_oms_nu)) - return fseries + return from_numpy_arrays(fr, psd, length, delta_f, low_freq_cutoff) def analytical_psd_lisa_tdi_1p5_T(length, delta_f, low_freq_cutoff, @@ -331,23 +315,19 @@ def analytical_psd_lisa_tdi_1p5_T(length, delta_f, low_freq_cutoff, The PSD of LISA's TDI-1.5 channel T. Notes ----- - Pease see Eq.(59) in for more details. + Please see Eq.(59) in for more details. """ len_arm = np.float64(len_arm) acc_noise_level = np.float64(acc_noise_level) oms_noise_level = np.float64(oms_noise_level) - psd = [] fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length) - for f in fr: - [s_acc_nu, s_oms_nu] = lisa_psd_components( - f, acc_noise_level, oms_noise_level) - omega_len = omega_length(f, len_arm) - psd.append(32*np.sin(omega_len)**2 * np.sin(omega_len/2)**2 * - (4*s_acc_nu*np.sin(omega_len/2)**2 + s_oms_nu)) - fseries = from_numpy_arrays(fr, np.array(psd), - length, delta_f, low_freq_cutoff) + s_acc_nu, s_oms_nu = lisa_psd_components( + fr, acc_noise_level, oms_noise_level) + omega_len = omega_length(fr, len_arm) + psd = (32*np.sin(omega_len)**2 * np.sin(omega_len/2)**2 * + (4*s_acc_nu*np.sin(omega_len/2)**2 + s_oms_nu)) - return fseries + return from_numpy_arrays(fr, psd, length, delta_f, low_freq_cutoff) def averaged_lisa_fplus_sq_approx(f, len_arm=2.5e9): @@ -367,7 +347,7 @@ def averaged_lisa_fplus_sq_approx(f, len_arm=2.5e9): The sky and polarization angle averaged squared antenna response. Notes ----- - Pease see Eq.(36) in for more details. + Please see Eq.(36) in for more details. """ from scipy.interpolate import interp1d from astropy.utils.data import download_file @@ -405,7 +385,7 @@ def averaged_response_lisa_tdi_1p5(f, len_arm=2.5e9): The sky and polarization angle averaged TDI-1.5 response to GW. Notes ----- - Pease see Eq.(39) in for more details. + Please see Eq.(39) in for more details. """ omega_len = omega_length(f, len_arm) ave_fp2 = averaged_lisa_fplus_sq_approx(f, len_arm) @@ -431,7 +411,7 @@ def averaged_response_lisa_tdi_2p0(f, len_arm=2.5e9): The sky and polarization angle averaged TDI-2.0 response to GW. Notes ----- - Pease see Eq.(40) in for more details. + Please see Eq.(40) in for more details. """ omega_len = omega_length(f, len_arm) response_tdi_1p5 = averaged_response_lisa_tdi_1p5(f, len_arm) @@ -469,21 +449,19 @@ def sensitivity_curve_lisa_semi_analytical(length, delta_f, low_freq_cutoff, LISA's sensitivity curve (6-links). Notes ----- - Pease see Eq.(42-43) in for more details. + Please see Eq.(42-43) in for more details. """ - sense_curve = [] len_arm = np.float64(len_arm) acc_noise_level = np.float64(acc_noise_level) oms_noise_level = np.float64(oms_noise_level) fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length) fp_sq = averaged_lisa_fplus_sq_approx(fr, len_arm) - for i in range(len(fr)): - [s_acc_nu, s_oms_nu] = lisa_psd_components( - fr[i], acc_noise_level, oms_noise_level) - omega_len = 2*np.pi*fr[i] * len_arm/constants.c.value - sense_curve.append((s_oms_nu + s_acc_nu*(3+np.cos(2*omega_len))) / - (omega_len**2*fp_sq[i])) - fseries = from_numpy_arrays(fr, np.array(sense_curve)/2, + s_acc_nu, s_oms_nu = lisa_psd_components( + fr, acc_noise_level, oms_noise_level) + omega_len = 2*np.pi*fr * len_arm/c.value + sense_curve = ((s_oms_nu + s_acc_nu*(3+np.cos(2*omega_len))) / + (omega_len**2*fp_sq)) + fseries = from_numpy_arrays(fr, sense_curve/2, length, delta_f, low_freq_cutoff) return fseries @@ -509,19 +487,15 @@ def sensitivity_curve_lisa_SciRD(length, delta_f, low_freq_cutoff): LISA's sensitivity curve in SciRD. Notes ----- - Pease see Eq.(114) in for more details. + Please see Eq.(114) in for more details. """ - sense_curve = [] fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length) - for f in fr: - s_I = 5.76e-48 * (1+(4e-4/f)**2) - s_II = 3.6e-41 - R = 1 + (f/2.5e-2)**2 - sense_curve.append(10/3 * (s_I/(2*np.pi*f)**4+s_II) * R) - fseries = from_numpy_arrays(fr, sense_curve, - length, delta_f, low_freq_cutoff) + s_I = 5.76e-48 * (1+(4e-4/fr)**2) + s_II = 3.6e-41 + R = 1 + (fr/2.5e-2)**2 + sense_curve = 10/3 * (s_I/(2*np.pi*fr)**4+s_II) * R - return fseries + return from_numpy_arrays(fr, sense_curve, length, delta_f, low_freq_cutoff) def sensitivity_curve_lisa_confusion(length, delta_f, low_freq_cutoff, @@ -557,7 +531,7 @@ def sensitivity_curve_lisa_confusion(length, delta_f, low_freq_cutoff, LISA's sensitivity curve with Galactic confusion noise. Notes ----- - Pease see Eq.(85-86) in for more details. + Please see Eq.(85-86) in for more details. """ if base_model == "semi": base_curve = sensitivity_curve_lisa_semi_analytical( @@ -571,13 +545,11 @@ def sensitivity_curve_lisa_confusion(length, delta_f, low_freq_cutoff, if duration < 0 or duration > 10: raise Exception("Must between 0 and 10.") fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length) - sh_confusion = [] f1 = 10**(-0.25*np.log10(duration)-2.7) fk = 10**(-0.27*np.log10(duration)-2.47) - for f in fr: - sh_confusion.append(0.5*1.14e-44*f**(-7/3)*np.exp(-(f/f1)**1.8) * - (1.0+np.tanh((fk-f)/(0.31e-3)))) - fseries_confusion = from_numpy_arrays(fr, np.array(sh_confusion), + sh_confusion = (0.5*1.14e-44*fr**(-7/3)*np.exp(-(fr/f1)**1.8) * + (1.0+np.tanh((fk-fr)/0.31e-3))) + fseries_confusion = from_numpy_arrays(fr, sh_confusion, length, delta_f, low_freq_cutoff) fseries = from_numpy_arrays(base_curve.sample_frequencies, base_curve+fseries_confusion, @@ -622,7 +594,7 @@ def sh_transformed_psd_lisa_tdi_XYZ(length, delta_f, low_freq_cutoff, noise, transformed from LISA sensitivity curve. Notes ----- - Pease see Eq.(7,41-43) in for more details. + Please see Eq.(7,41-43) in for more details. """ fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length) if tdi == "1.5": @@ -680,13 +652,11 @@ def semi_analytical_psd_lisa_confusion_noise(length, delta_f, low_freq_cutoff, raise Exception("The version of TDI, currently only for 1.5 or 2.0.") fseries_response = from_numpy_arrays(fr, np.array(response), length, delta_f, low_freq_cutoff) - sh_confusion = [] f1 = 10**(-0.25*np.log10(duration)-2.7) fk = 10**(-0.27*np.log10(duration)-2.47) - for f in fr: - sh_confusion.append(0.5*1.14e-44*f**(-7/3)*np.exp(-(f/f1)**1.8) * - (1.0+np.tanh((fk-f)/(0.31e-3)))) - fseries_confusion = from_numpy_arrays(fr, np.array(sh_confusion), + sh_confusion = (0.5*1.14e-44*fr**(-7/3)*np.exp(-(fr/f1)**1.8) * + (1.0+np.tanh((fk-fr)/0.31e-3))) + fseries_confusion = from_numpy_arrays(fr, sh_confusion, length, delta_f, low_freq_cutoff) psd_confusion = 2*fseries_confusion.data * fseries_response.data fseries = from_numpy_arrays(fseries_confusion.sample_frequencies, diff --git a/pycbc/results/pygrb_postprocessing_utils.py b/pycbc/results/pygrb_postprocessing_utils.py index 8760a025224..531f7124a40 100644 --- a/pycbc/results/pygrb_postprocessing_utils.py +++ b/pycbc/results/pygrb_postprocessing_utils.py @@ -310,16 +310,17 @@ def _get_id_numbers(ligolw_table, column): # ============================================================================= # Function to build a dictionary (indexed by ifo) of time-slid vetoes # ============================================================================= -def _slide_vetoes(vetoes, slide_dict_or_list, slide_id): +def _slide_vetoes(vetoes, slide_dict_or_list, slide_id, ifos): """Build a dictionary (indexed by ifo) of time-slid vetoes""" # Copy vetoes - slid_vetoes = copy.deepcopy(vetoes) - - # Slide them - ifos = vetoes.keys() - for ifo in ifos: - slid_vetoes[ifo].shift(-slide_dict_or_list[slide_id][ifo]) + if vetoes is not None: + slid_vetoes = copy.deepcopy(vetoes) + # Slide them + for ifo in ifos: + slid_vetoes[ifo].shift(-slide_dict_or_list[slide_id][ifo]) + else: + slid_vetoes = {ifo: segments.segmentlist() for ifo in ifos} return slid_vetoes @@ -461,36 +462,39 @@ def sort_trigs(trial_dict, trigs, slide_dict, seg_dict): sorted_trigs = {} # Begin by sorting the triggers into each slide - # New seems pretty slow, so run it once and then use deepcopy - tmp_table = glsctables.New(glsctables.MultiInspiralTable) for slide_id in slide_dict: - sorted_trigs[slide_id] = copy.deepcopy(tmp_table) - for trig in trigs: - sorted_trigs[int(trig.time_slide_id)].append(trig) + sorted_trigs[slide_id] = [] + for slide_id, event_id in zip(trigs['network/slide_id'], + trigs['network/event_id']): + sorted_trigs[slide_id].append(event_id) for slide_id in slide_dict: # These can only *reduce* the analysis time curr_seg_list = seg_dict[slide_id] # Check the triggers are all in the analysed segment lists - for trig in sorted_trigs[slide_id]: - if trig.end_time not in curr_seg_list: + for event_id in sorted_trigs[slide_id]: + index = numpy.flatnonzero(trigs['network/event_id'] == event_id)[0] + end_time = trigs['network/end_time_gc'][index] + if end_time not in curr_seg_list: # This can be raised if the trigger is on the segment boundary, # so check if the trigger is within 1/100 of a second within # the list - if trig.get_end() + 0.01 in curr_seg_list: + if end_time + 0.01 in curr_seg_list: continue - if trig.get_end() - 0.01 in curr_seg_list: + if end_time - 0.01 in curr_seg_list: continue err_msg = "Triggers found in input files not in the list of " err_msg += "analysed segments. This should not happen." raise RuntimeError(err_msg) # END OF CHECK # - # The below line works like the inverse of .veto and only returns trigs - # that are within the segment specified by trial_dict[slide_id] - sorted_trigs[slide_id] = \ - sorted_trigs[slide_id].vetoed(trial_dict[slide_id]) + # Keep triggers that are in trial_dict + sorted_trigs[slide_id] = [event_id for event_id in + sorted_trigs[slide_id] + if trigs['network/end_time_gc'][ + trigs['network/event_id'] == event_id][0] + in trial_dict[slide_id]] return sorted_trigs @@ -498,8 +502,7 @@ def sort_trigs(trial_dict, trigs, slide_dict, seg_dict): # ============================================================================= # Extract basic trigger properties and store them as dictionaries # ============================================================================= -def extract_basic_trig_properties(trial_dict, trigs, slide_dict, seg_dict, - opts): +def extract_basic_trig_properties(trial_dict, trigs, slide_dict, seg_dict): """Extract and store as dictionaries time, SNR, and BestNR of time-slid triggers""" @@ -507,16 +510,6 @@ def extract_basic_trig_properties(trial_dict, trigs, slide_dict, seg_dict, sorted_trigs = sort_trigs(trial_dict, trigs, slide_dict, seg_dict) logging.info("Triggers sorted.") - # Local copies of variables entering the BestNR definition - chisq_index = opts.chisq_index - chisq_nhigh = opts.chisq_nhigh - null_thresh = list(map(float, opts.null_snr_threshold.split(','))) - snr_thresh = opts.snr_threshold - sngl_snr_thresh = opts.sngl_snr_threshold - new_snr_thresh = opts.newsnr_threshold - null_grad_thresh = opts.null_grad_thresh - null_grad_val = opts.null_grad_val - # Build the 3 dictionaries trig_time = {} trig_snr = {} @@ -524,21 +517,12 @@ def extract_basic_trig_properties(trial_dict, trigs, slide_dict, seg_dict, for slide_id in slide_dict: slide_trigs = sorted_trigs[slide_id] if slide_trigs: - trig_time[slide_id] = numpy.asarray(slide_trigs.get_end()).\ - astype(float) - trig_snr[slide_id] = numpy.asarray(slide_trigs.get_column('snr')) + trig_time[slide_id] = trigs['network/end_time_gc'][slide_trigs] + trig_snr[slide_id] = trigs['network/coherent_snr'][slide_trigs] else: trig_time[slide_id] = numpy.asarray([]) trig_snr[slide_id] = numpy.asarray([]) - trig_bestnr[slide_id] = get_bestnrs(slide_trigs, - q=chisq_index, - n=chisq_nhigh, - null_thresh=null_thresh, - snr_threshold=snr_thresh, - sngl_snr_threshold=sngl_snr_thresh, - chisq_threshold=new_snr_thresh, - null_grad_thresh=null_grad_thresh, - null_grad_val=null_grad_val) + trig_bestnr[slide_id] = trigs['network/reweighted_snr'][slide_trigs] logging.info("Time, SNR, and BestNR of triggers extracted.") return trig_time, trig_snr, trig_bestnr @@ -699,7 +683,7 @@ def construct_trials(seg_files, seg_dict, ifos, slide_dict, vetoes): seg_buffer.coalesce() # Construct the ifo-indexed dictionary of slid veteoes - slid_vetoes = _slide_vetoes(vetoes, slide_dict, slide_id) + slid_vetoes = _slide_vetoes(vetoes, slide_dict, slide_id, ifos) # Construct trial list and check against buffer trial_dict[slide_id] = segments.segmentlist() @@ -804,3 +788,24 @@ def get_coinc_snr(trigs_or_injs, ifos): coinc_snr = numpy.sqrt(snr_sum_square) return coinc_snr + + +def template_hash_to_id(trigger_file, bank_path): + """ + This function converts the template hashes from a trigger file + into 'template_id's that represent indices of the + templates within the bank. + Parameters + ---------- + trigger_file: h5py File object for trigger file + bank_file: filepath for template bank + """ + with h5py.File(bank_path, "r") as bank: + hashes = bank['template_hash'][:] + ifos = [k for k in trigger_file.keys() if k != 'network'] + trig_hashes = trigger_file[f'{ifos[0]}/template_hash'][:] + trig_ids = numpy.zeros(trig_hashes.shape[0], dtype=int) + for idx, t_hash in enumerate(hashes): + matches = numpy.where(trig_hashes == t_hash) + trig_ids[matches] = idx + return trig_ids diff --git a/pycbc/results/render.py b/pycbc/results/render.py index d10c4989e6d..3f58c4bf6fb 100644 --- a/pycbc/results/render.py +++ b/pycbc/results/render.py @@ -1,5 +1,3 @@ -#!/usr/bin/python - # Copyright (C) 2015 Christopher M. Biwer # # This program is free software; you can redistribute it and/or modify it diff --git a/pycbc/results/versioning.py b/pycbc/results/versioning.py index 64536377d22..529b39eeb91 100644 --- a/pycbc/results/versioning.py +++ b/pycbc/results/versioning.py @@ -1,5 +1,3 @@ -#!/usr/bin/python - # Copyright (C) 2015 Ian Harry # # This program is free software; you can redistribute it and/or modify it @@ -20,8 +18,9 @@ import subprocess import urllib.parse -import lal, lalframe -import pycbc.version, glue.git_version +import lal +import lalframe +import pycbc.version def get_library_version_info(): """This will return a list of dictionaries containing versioning @@ -88,19 +87,6 @@ def add_info_new_version(info_dct, curr_module, extra_str): pass library_list.append(lalsimulationinfo) - glueinfo = {} - glueinfo['Name'] = 'LSCSoft-Glue' - glueinfo['ID'] = glue.git_version.id - glueinfo['Status'] = glue.git_version.status - glueinfo['Version'] = glue.git_version.version - glueinfo['Tag'] = glue.git_version.tag - glueinfo['Author'] = glue.git_version.author - glueinfo['Builder'] = glue.git_version.builder - glueinfo['Branch'] = glue.git_version.branch - glueinfo['Committer'] = glue.git_version.committer - glueinfo['Date'] = glue.git_version.date - library_list.append(glueinfo) - pycbcinfo = {} pycbcinfo['Name'] = 'PyCBC' pycbcinfo['ID'] = pycbc.version.version diff --git a/pycbc/tmpltbank/bank_conversions.py b/pycbc/tmpltbank/bank_conversions.py index 9859d4e7a0a..f987d1df0eb 100644 --- a/pycbc/tmpltbank/bank_conversions.py +++ b/pycbc/tmpltbank/bank_conversions.py @@ -26,10 +26,14 @@ specific values from PyCBC template banks. """ +import logging import numpy as np + from pycbc import conversions as conv from pycbc import pnutils +logger = logging.getLogger('pycbc.tmpltbank.bank_conversions') + # Convert from parameter name to helper function # some multiple names are used for the same function conversion_options = ['mass1', 'mass2', 'spin1z', 'spin2z', 'duration', diff --git a/pycbc/tmpltbank/bank_output_utils.py b/pycbc/tmpltbank/bank_output_utils.py index 2fae646fb1e..4d5e281ad16 100644 --- a/pycbc/tmpltbank/bank_output_utils.py +++ b/pycbc/tmpltbank/bank_output_utils.py @@ -1,7 +1,10 @@ +import logging import numpy import h5py + from lal import PI, MTSUN_SI, TWOPI, GAMMA from ligo.lw import ligolw, lsctables, utils as ligolw_utils + from pycbc import pnutils from pycbc.tmpltbank.lambda_mapping import ethinca_order_from_string from pycbc.io.ligolw import ( @@ -10,6 +13,8 @@ from pycbc.waveform import get_waveform_filter_length_in_time as gwflit +logger = logging.getLogger('pycbc.tmpltbank.bank_output_utils') + def convert_to_sngl_inspiral_table(params, proc_id): ''' Convert a list of m1,m2,spin1z,spin2z values into a basic sngl_inspiral diff --git a/pycbc/tmpltbank/brute_force_methods.py b/pycbc/tmpltbank/brute_force_methods.py index c5ff02858c0..34ca473c5df 100644 --- a/pycbc/tmpltbank/brute_force_methods.py +++ b/pycbc/tmpltbank/brute_force_methods.py @@ -1,6 +1,10 @@ +import logging import numpy + from pycbc.tmpltbank.coord_utils import get_cov_params +logger = logging.getLogger('pycbc.tmpltbank.brute_force_methods') + def get_physical_covaried_masses(xis, bestMasses, bestXis, req_match, massRangeParams, metricParams, fUpper, diff --git a/pycbc/tmpltbank/calc_moments.py b/pycbc/tmpltbank/calc_moments.py index 305fb8834e6..c472dfacd91 100644 --- a/pycbc/tmpltbank/calc_moments.py +++ b/pycbc/tmpltbank/calc_moments.py @@ -14,9 +14,13 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +import logging import numpy + from pycbc.tmpltbank.lambda_mapping import generate_mapping +logger = logging.getLogger('pycbc.tmpltbank.calc_moments') + def determine_eigen_directions(metricParams, preserveMoments=False, vary_fmax=False, vary_density=None): diff --git a/pycbc/tmpltbank/coord_utils.py b/pycbc/tmpltbank/coord_utils.py index c1bb38a7f62..0d0e4db509d 100644 --- a/pycbc/tmpltbank/coord_utils.py +++ b/pycbc/tmpltbank/coord_utils.py @@ -16,11 +16,15 @@ import logging import numpy + from pycbc.tmpltbank.lambda_mapping import get_chirp_params from pycbc import conversions from pycbc import pnutils from pycbc.neutron_stars import load_ns_sequence +logger = logging.getLogger('pycbc.tmpltbank.coord_utils') + + def estimate_mass_range(numPoints, massRangeParams, metricParams, fUpper,\ covary=True): """ @@ -235,7 +239,7 @@ def get_random_mass(numPoints, massRangeParams, eos='2H'): warn_msg += "(%s). " %(max_ns_g_mass) warn_msg += "The code will proceed using the latter value " warn_msg += "as the boundary mass." - logging.warn(warn_msg) + logger.warn(warn_msg) boundary_mass = max_ns_g_mass # Empty arrays to store points that pass all cuts @@ -711,7 +715,7 @@ def find_max_and_min_frequencies(name, mass_range_params, freqs): warn_msg += "for the metric: %s Hz. " %(freqs.min()) warn_msg += "Distances for these waveforms will be calculated at " warn_msg += "the lowest available metric frequency." - logging.warn(warn_msg) + logger.warn(warn_msg) if upper_f_cutoff > freqs.max(): warn_msg = "WARNING: " warn_msg += "Highest frequency cutoff is %s Hz " %(upper_f_cutoff,) @@ -719,7 +723,7 @@ def find_max_and_min_frequencies(name, mass_range_params, freqs): warn_msg += "for the metric: %s Hz. " %(freqs.max()) warn_msg += "Distances for these waveforms will be calculated at " warn_msg += "the largest available metric frequency." - logging.warn(warn_msg) + logger.warn(warn_msg) return find_closest_calculated_frequencies(cutoffs, freqs) diff --git a/pycbc/tmpltbank/lambda_mapping.py b/pycbc/tmpltbank/lambda_mapping.py index ed0121ed912..3e3516dbc96 100644 --- a/pycbc/tmpltbank/lambda_mapping.py +++ b/pycbc/tmpltbank/lambda_mapping.py @@ -14,12 +14,17 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. import re +import logging import numpy -import pycbc.libutils + from lal import MTSUN_SI, PI, CreateREAL8Vector +import pycbc.libutils + lalsimulation = pycbc.libutils.import_optional('lalsimulation') +logger = logging.getLogger('pycbc.tmpltbank.lambda_mapping') + # PLEASE ENSURE THESE ARE KEPT UP TO DATE WITH THE REST OF THIS FILE pycbcValidTmpltbankOrders = ['zeroPN','onePN','onePointFivePN','twoPN',\ 'twoPointFivePN','threePN','threePointFivePN'] diff --git a/pycbc/tmpltbank/lattice_utils.py b/pycbc/tmpltbank/lattice_utils.py index c938a8b3612..d6db4f98c86 100644 --- a/pycbc/tmpltbank/lattice_utils.py +++ b/pycbc/tmpltbank/lattice_utils.py @@ -14,10 +14,14 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +import logging import copy import numpy + import lal +logger = logging.getLogger('pycbc.tmpltbank.lattice_utils') + def generate_hexagonal_lattice(maxv1, minv1, maxv2, minv2, mindist): """ This function generates a 2-dimensional lattice of points using a hexagonal diff --git a/pycbc/tmpltbank/option_utils.py b/pycbc/tmpltbank/option_utils.py index 722c250b205..e9451a7e716 100644 --- a/pycbc/tmpltbank/option_utils.py +++ b/pycbc/tmpltbank/option_utils.py @@ -19,11 +19,15 @@ import textwrap import numpy import os + from pycbc.tmpltbank.lambda_mapping import get_ethinca_orders, pycbcValidOrdersHelpDescriptions from pycbc import pnutils from pycbc.neutron_stars import load_ns_sequence from pycbc.types import positive_float, nonnegative_float +logger = logging.getLogger('pycbc.tmpltbank.option_utils') + + class IndentedHelpFormatterWithNL(argparse.ArgumentDefaultsHelpFormatter): """ This class taken from @@ -589,21 +593,21 @@ def verify_mass_range_options(opts, parser, nonSpin=False): # inform him/her about the caveats involved in this. if hasattr(opts, 'remnant_mass_threshold') \ and opts.remnant_mass_threshold is not None: - logging.info("""You have asked to exclude EM dim NS-BH systems from the - target parameter space. The script will assume that m1 is - the BH and m2 is the NS: make sure that your settings - respect this convention. The script will also treat the - NS as non-spinning: use NS spins in the template bank - at your own risk!""") + logger.info("""You have asked to exclude EM dim NS-BH systems from the + target parameter space. The script will assume that m1 + is the BH and m2 is the NS: make sure that your settings + respect this convention. The script will also treat the + NS as non-spinning: use NS spins in the template bank + at your own risk!""") if opts.use_eos_max_ns_mass: - logging.info("""You have asked to take into account the maximum NS + logger.info("""You have asked to take into account the maximum NS mass value for the EOS in use.""") # Find out if the EM constraint surface data already exists or not # and inform user whether this will be read from file or generated. # This is the minumum eta as a function of BH spin and NS mass # required to produce an EM counterpart if os.path.isfile('constraint_em_bright.npz'): - logging.info("""The constraint surface for EM bright binaries + logger.info("""The constraint surface for EM bright binaries will be read in from constraint_em_bright.npz.""") # Assign min/max total mass from mass1, mass2 if not specified diff --git a/pycbc/tmpltbank/partitioned_bank.py b/pycbc/tmpltbank/partitioned_bank.py index 29505b2a03e..9b22658aaa8 100644 --- a/pycbc/tmpltbank/partitioned_bank.py +++ b/pycbc/tmpltbank/partitioned_bank.py @@ -17,8 +17,11 @@ import copy import numpy import logging + from pycbc.tmpltbank import coord_utils +logger = logging.getLogger('pycbc.tmpltbank.partitioned_bank') + class PartitionedTmpltbank(object): """ This class is used to hold a template bank partitioned into numerous bins @@ -524,7 +527,7 @@ def add_point_by_masses(self, mass1, mass2, spin1z, spin2z, warn_msg += "convention is that mass1 > mass2. Swapping mass1 " warn_msg += "and mass2 and adding point to bank. This message " warn_msg += "will not be repeated." - logging.warn(warn_msg) + logger.warn(warn_msg) self.spin_warning_given = True # These that masses obey the restrictions of mass_range_params diff --git a/pycbc/transforms.py b/pycbc/transforms.py index c8fa44a9a9b..487de5ab103 100644 --- a/pycbc/transforms.py +++ b/pycbc/transforms.py @@ -318,6 +318,98 @@ def from_config(cls, cp, section, outputs): return cls(inputs, outputs, transform_functions, jacobian=jacobian) +class CustomTransformMultiOutputs(CustomTransform): + """Allows for any transform to be defined. Based on CustomTransform, + but also supports multi-returning value functions. + + Parameters + ---------- + input_args : (list of) str + The names of the input parameters. + output_args : (list of) str + The names of the output parameters. + transform_functions : dict + Dictionary mapping input args to a string giving a function call; + e.g., ``{'q': 'q_from_mass1_mass2(mass1, mass2)'}``. + jacobian : str, optional + String giving a jacobian function. The function must be in terms of + the input arguments. + """ + + name = "custom_multi" + + def __init__(self, input_args, output_args, transform_functions, + jacobian=None): + super(CustomTransformMultiOutputs, self).__init__( + input_args, output_args, transform_functions, jacobian) + + def transform(self, maps): + """Applies the transform functions to the given maps object. + Parameters + ---------- + maps : dict, or FieldArray + Returns + ------- + dict or FieldArray + A map object containing the transformed variables, along with the + original variables. The type of the output will be the same as the + input. + """ + if self.transform_functions is None: + raise NotImplementedError("no transform function(s) provided") + # copy values to scratch + self._copytoscratch(maps) + # ensure that we return the same data type in each dict + getslice = self._getslice(maps) + # evaluate the functions + # func[0] is the function itself, func[1] is the index, + # this supports multiple returning values function + out = { + p: self._scratch[func[0]][func[1]][getslice] if + len(self._scratch[func[0]]) > 1 else + self._scratch[func[0]][getslice] + for p, func in self.transform_functions.items() + } + return self.format_output(maps, out) + + @classmethod + def from_config(cls, cp, section, outputs): + """Loads a CustomTransformMultiOutputs from the given config file. + + Example section: + + .. code-block:: ini + + [{section}-outvar1+outvar2] + name = custom_multi + inputs = inputvar1, inputvar2 + outvar1, outvar2 = func1(inputs) + jacobian = func2(inputs) + """ + tag = outputs + outputs = list(outputs.split(VARARGS_DELIM)) + all_vars = ", ".join(outputs) + inputs = map(str.strip, + cp.get_opt_tag(section, "inputs", tag).split(",")) + # get the functions for each output + transform_functions = {} + output_index = slice(None, None, None) + for var in outputs: + # check if option can be cast as a float + try: + func = cp.get_opt_tag(section, var, tag) + except Exception: + func = cp.get_opt_tag(section, all_vars, tag) + output_index = slice(outputs.index(var), outputs.index(var)+1) + transform_functions[var] = [func, output_index] + s = "-".join([section, tag]) + if cp.has_option(s, "jacobian"): + jacobian = cp.get_opt_tag(section, "jacobian", tag) + else: + jacobian = None + return cls(inputs, outputs, transform_functions, jacobian=jacobian) + + # # ============================================================================= # @@ -2725,6 +2817,7 @@ def from_config(cls, cp, section, outputs, # dictionary of all transforms transforms = { CustomTransform.name: CustomTransform, + CustomTransformMultiOutputs.name: CustomTransformMultiOutputs, MchirpQToMass1Mass2.name: MchirpQToMass1Mass2, Mass1Mass2ToMchirpQ.name: Mass1Mass2ToMchirpQ, MchirpEtaToMass1Mass2.name: MchirpEtaToMass1Mass2, diff --git a/pycbc/types/timeseries.py b/pycbc/types/timeseries.py index 3bf94065473..70d8802cb23 100644 --- a/pycbc/types/timeseries.py +++ b/pycbc/types/timeseries.py @@ -188,10 +188,10 @@ def time_slice(self, start, end, mode='floor'): start_idx = float(start - self.start_time) * self.sample_rate end_idx = float(end - self.start_time) * self.sample_rate - if _numpy.isclose(start_idx, round(start_idx)): + if _numpy.isclose(start_idx, round(start_idx), rtol=0, atol=1E-3): start_idx = round(start_idx) - if _numpy.isclose(end_idx, round(end_idx)): + if _numpy.isclose(end_idx, round(end_idx), rtol=0, atol=1E-3): end_idx = round(end_idx) if mode == 'floor': diff --git a/pycbc/workflow/core.py b/pycbc/workflow/core.py index a455b37c81e..2853e50d528 100644 --- a/pycbc/workflow/core.py +++ b/pycbc/workflow/core.py @@ -36,7 +36,6 @@ import lal import lal.utils import Pegasus.api # Try and move this into pegasus_workflow -from glue import lal as gluelal from ligo import segments from ligo.lw import lsctables, ligolw from ligo.lw import utils as ligolw_utils @@ -1572,6 +1571,8 @@ def convert_to_lal_cache(self): """ Return all files in this object as a glue.lal.Cache object """ + from glue import lal as gluelal + lal_cache = gluelal.Cache([]) for entry in self: try: diff --git a/pycbc/workflow/datafind.py b/pycbc/workflow/datafind.py index 9dc898effcd..42055d25949 100644 --- a/pycbc/workflow/datafind.py +++ b/pycbc/workflow/datafind.py @@ -34,7 +34,6 @@ import urllib.parse from ligo import segments from ligo.lw import utils, table -from glue import lal from gwdatafind import find_urls as find_frame_urls from pycbc.workflow.core import SegFile, File, FileList, make_analysis_dir from pycbc.io.ligolw import LIGOLWContentHandler @@ -694,6 +693,8 @@ def setup_datafind_from_pregenerated_lcf_files(cp, ifos, outputDir, tags=None): datafindOuts : pycbc.workflow.core.FileList List of all the datafind output files for use later in the pipeline. """ + from glue import lal + if tags is None: tags = [] @@ -821,6 +822,8 @@ def get_missing_segs_from_frame_file_cache(datafindcaches): missingFrames: Dict. of ifo keyed lal.Cache instances The list of missing frames """ + from glue import lal + missingFrameSegs = {} missingFrames = {} for cache in datafindcaches: @@ -947,6 +950,8 @@ def run_datafind_instance(cp, outputDir, observatory, frameType, Cache file listing all of the datafind output files for use later in the pipeline. """ + from glue import lal + if tags is None: tags = [] diff --git a/pycbc/workflow/grb_utils.py b/pycbc/workflow/grb_utils.py index 8f96f793d39..2d2113516b3 100644 --- a/pycbc/workflow/grb_utils.py +++ b/pycbc/workflow/grb_utils.py @@ -233,8 +233,8 @@ def get_sky_grid_scale( return out -def setup_pygrb_pp_workflow(wf, pp_dir, seg_dir, segment, insp_files, - inj_files, inj_insp_files, inj_tags): +def setup_pygrb_pp_workflow(wf, pp_dir, seg_dir, segment, bank_file, + insp_files, inj_files, inj_insp_files, inj_tags): """ Generate post-processing section of PyGRB offline workflow """ @@ -258,7 +258,7 @@ def setup_pygrb_pp_workflow(wf, pp_dir, seg_dir, segment, insp_files, job_instance = exe_class(wf.cp, "trig_combiner") # Create node for coherent no injections jobs node, trig_files = job_instance.create_node(wf.ifos, seg_dir, segment, - insp_files, pp_dir) + insp_files, pp_dir, bank_file) wf.add_node(node) pp_outs.append(trig_files) @@ -285,7 +285,7 @@ def setup_pygrb_pp_workflow(wf, pp_dir, seg_dir, segment, insp_files, if inj_tag in f.tags[1]]) node, inj_find_file = job_instance.create_node( tag_inj_files, tag_insp_files, - pp_dir) + bank_file, pp_dir) wf.add_node(node) inj_find_files.append(inj_find_file) pp_outs.append(inj_find_files) @@ -321,7 +321,7 @@ def __init__(self, cp, name): self.num_trials = int(cp.get('trig_combiner', 'num-trials')) def create_node(self, ifo_tag, seg_dir, segment, insp_files, - out_dir, tags=None): + out_dir, bank_file, tags=None): node = Node(self) node.add_opt('--verbose') node.add_opt("--ifo-tag", ifo_tag) @@ -331,6 +331,7 @@ def create_node(self, ifo_tag, seg_dir, segment, insp_files, node.add_input_list_opt("--input-files", insp_files) node.add_opt("--user-tag", "PYGRB") node.add_opt("--num-trials", self.num_trials) + node.add_input_opt("--bank-file", bank_file) # Prepare output file tag user_tag = f"PYGRB_GRB{self.trigger_name}" if tags: @@ -385,13 +386,14 @@ class PycbcGrbInjFinderExecutable(Executable): def __init__(self, cp, exe_name): super().__init__(cp=cp, name=exe_name) - def create_node(self, inj_files, inj_insp_files, + def create_node(self, inj_files, inj_insp_files, bank_file, out_dir, tags=None): if tags is None: tags = [] node = Node(self) node.add_input_list_opt('--input-files', inj_insp_files) node.add_input_list_opt('--inj-files', inj_files) + node.add_input_opt('--bank-file', bank_file) ifo_tag, desc, segment = filename_metadata(inj_files[0].name) desc = '_'.join(desc.split('_')[:-1]) out_name = "{}-{}_FOUNDMISSED-{}-{}.h5".format( diff --git a/pycbc/workflow/pegasus_workflow.py b/pycbc/workflow/pegasus_workflow.py index e442072e28a..9a80ff09fcd 100644 --- a/pycbc/workflow/pegasus_workflow.py +++ b/pycbc/workflow/pegasus_workflow.py @@ -75,7 +75,7 @@ def set_num_retries(self, number): self.add_profile("dagman", "retry", number) def set_execution_site(self, site): - self.add_profile("selector", "execution_site", site) + self.add_profile("selector", "execution.site", site) class Executable(ProfileShortcuts): diff --git a/pycbc/workflow/plotting.py b/pycbc/workflow/plotting.py index 165786c5715..cdf71b4119a 100644 --- a/pycbc/workflow/plotting.py +++ b/pycbc/workflow/plotting.py @@ -534,7 +534,6 @@ def make_singles_plot(workflow, trig_files, bank_file, veto_file, veto_name, 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([]) node = PlotExecutable(workflow.cp, 'plot_dq_flag_likelihood', ifos=dq_file.ifo, out_dir=out_dir, tags=tags).create_node() @@ -543,5 +542,29 @@ def make_dq_flag_trigger_rate_plot(workflow, dq_file, dq_label, out_dir, tags=No 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 + return node.output_files[0] + + +def make_dq_segment_table(workflow, dq_file, out_dir, tags=None): + tags = [] if tags is None else tags + makedir(out_dir) + node = PlotExecutable(workflow.cp, 'page_dq_table', 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, '.html', '--output-file') + workflow += node + return node.output_files[0] + + +def make_template_bin_table(workflow, dq_file, out_dir, tags=None): + tags = [] if tags is None else tags + makedir(out_dir) + node = PlotExecutable(workflow.cp, 'page_template_bin_table', + 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, '.html', '--output-file') + workflow += node + return node.output_files[0] diff --git a/tools/pycbc_test_suite.sh b/tools/pycbc_test_suite.sh index 08ef955dc49..6e1f745194e 100755 --- a/tools/pycbc_test_suite.sh +++ b/tools/pycbc_test_suite.sh @@ -31,7 +31,7 @@ fi if [ "$PYCBC_TEST_TYPE" = "help" ] || [ -z ${PYCBC_TEST_TYPE+x} ]; then # check that all executables that do not require # special environments can return a help message - for prog in `find ${PATH//:/ } -maxdepth 1 -name 'pycbc*' -print 2>/dev/null | egrep -v '(pycbc_live_nagios_monitor|pycbc_make_offline_grb_workflow|pycbc_mvsc_get_features|pycbc_upload_xml_to_gracedb|pycbc_coinc_time)'` + for prog in `find ${PATH//:/ } -maxdepth 1 -name 'pycbc*' -print 2>/dev/null | egrep -v '(pycbc_live_nagios_monitor|pycbc_make_offline_grb_workflow|pycbc_mvsc_get_features|pycbc_upload_xml_to_gracedb|pycbc_coinc_time)' | sort | uniq` do echo -e ">> [`date`] running $prog --help" $prog --help &> $LOG_FILE @@ -45,6 +45,18 @@ if [ "$PYCBC_TEST_TYPE" = "help" ] || [ -z ${PYCBC_TEST_TYPE+x} ]; then echo -e " Pass." fi done + # also check that --version works for one executable + echo -e ">> [`date`] running pycbc_inspiral --version" + pycbc_inspiral --version &> $LOG_FILE + if test $? -ne 0 ; then + RESULT=1 + echo -e " FAILED!" + echo -e "---------------------------------------------------------" + cat $LOG_FILE + echo -e "---------------------------------------------------------" + else + echo -e " Pass." + fi fi if [ "$PYCBC_TEST_TYPE" = "search" ] || [ -z ${PYCBC_TEST_TYPE+x} ]; then