Hello community, here is the log from the commit of package python-PyCBC for openSUSE:Factory checked in at 2020-06-08 23:58:54 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Comparing /work/SRC/openSUSE:Factory/python-PyCBC (Old) and /work/SRC/openSUSE:Factory/.python-PyCBC.new.3606 (New) ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Package is "python-PyCBC" Mon Jun 8 23:58:54 2020 rev:2 rq:810052 version:1.16.1 Changes: -------- --- /work/SRC/openSUSE:Factory/python-PyCBC/python-PyCBC.changes 2020-05-26 17:17:09.727695685 +0200 +++ /work/SRC/openSUSE:Factory/.python-PyCBC.new.3606/python-PyCBC.changes 2020-06-09 00:01:52.468696842 +0200 @@ -1,0 +2,7 @@ +Wed May 27 20:51:44 UTC 2020 - Atri Bhattacharya <[email protected]> + +- Update to version 0.16.1: + * Backward incompatible Changes to the offline search analysis, + notably in the naming of ranking statistics + +------------------------------------------------------------------- Old: ---- PyCBC-1.15.6.tar.gz New: ---- PyCBC-1.16.1.tar.gz ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Other differences: ------------------ ++++++ python-PyCBC.spec ++++++ --- /var/tmp/diff_new_pack.2mkdSA/_old 2020-06-09 00:01:54.024702363 +0200 +++ /var/tmp/diff_new_pack.2mkdSA/_new 2020-06-09 00:01:54.028702376 +0200 @@ -1,7 +1,7 @@ # # spec file for package python-PyCBC # -# Copyright (c) 2019 SUSE LINUX GmbH, Nuernberg, Germany. +# Copyright (c) 2020 SUSE LLC # # All modifications and additions to the file contributed by third parties # remain the property of their copyright owners, unless otherwise agreed @@ -12,24 +12,26 @@ # license that conforms to the Open Source Definition (Version 1.9) # published by the Open Source Initiative. -# Please submit bugfixes or comments via http://bugs.opensuse.org/ +# Please submit bugfixes or comments via https://bugs.opensuse.org/ +# + Name: python-PyCBC -Version: 1.15.6 +Version: 1.16.1 Release: 0 -License: GPL-3.0-or-later Summary: Core library to analyze gravitational-wave data -Url: http://www.pycbc.org/ +License: GPL-3.0-or-later Group: Development/Languages/Python +URL: http://www.pycbc.org/ Source: https://files.pythonhosted.org/packages/source/p/pycbc/PyCBC-%{version}.tar.gz -BuildRequires: gcc-c++ -BuildRequires: python-rpm-macros -BuildRequires: %{python_module devel} BuildRequires: %{python_module Cython} -BuildRequires: %{python_module setuptools} -BuildRequires: fdupes +BuildRequires: %{python_module devel} BuildRequires: %{python_module numpy >= 1.16.0} BuildRequires: %{python_module numpy-devel >= 1.16.0} +BuildRequires: %{python_module setuptools} +BuildRequires: fdupes +BuildRequires: gcc-c++ +BuildRequires: python-rpm-macros Requires: python-numpy >= 1.16.0 ExclusiveArch: %ix86 x86_64 ++++++ PyCBC-1.15.6.tar.gz -> PyCBC-1.16.1.tar.gz ++++++ diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/PKG-INFO new/PyCBC-1.16.1/PKG-INFO --- old/PyCBC-1.15.6/PKG-INFO 2020-04-17 11:07:50.000000000 +0200 +++ new/PyCBC-1.16.1/PKG-INFO 2020-04-28 20:25:45.000000000 +0200 @@ -1,12 +1,12 @@ Metadata-Version: 2.1 Name: PyCBC -Version: 1.15.6 +Version: 1.16.1 Summary: Core library to analyze gravitational-wave data, find signals, and study their parameters. Home-page: http://www.pycbc.org/ Author: The PyCBC team Author-email: [email protected] License: UNKNOWN -Download-URL: https://github.com/gwastro/pycbc/tarball/v1.15.6 +Download-URL: https://github.com/gwastro/pycbc/tarball/v1.16.1 Description: `PyCBC <http://pycbc.org>`_ is a software package used to explore astrophysical sources of gravitational waves. It contains algorithms to analyze gravitational-wave data from the LIGO and Virgo detectors, detect coalescing compact binaries, and measure the astrophysical parameters of detected sources. PyCBC was used in the `first direct detection of gravitational waves <https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.116.061102>`_ and is used in the flagship analysis of LIGO and Virgo data. PyCBC is developed collaboratively and lead by a team of LIGO scientists with the aim to build accessible tools for gravitational-wave data analysis. One of the easiest ways to get a full software environment is by `downloading one of our docker images. <http://pycbc.org/pycbc/latest/html/docker.html>`_ diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/PyCBC.egg-info/PKG-INFO new/PyCBC-1.16.1/PyCBC.egg-info/PKG-INFO --- old/PyCBC-1.15.6/PyCBC.egg-info/PKG-INFO 2020-04-17 11:07:50.000000000 +0200 +++ new/PyCBC-1.16.1/PyCBC.egg-info/PKG-INFO 2020-04-28 20:25:44.000000000 +0200 @@ -1,12 +1,12 @@ Metadata-Version: 2.1 Name: PyCBC -Version: 1.15.6 +Version: 1.16.1 Summary: Core library to analyze gravitational-wave data, find signals, and study their parameters. Home-page: http://www.pycbc.org/ Author: The PyCBC team Author-email: [email protected] License: UNKNOWN -Download-URL: https://github.com/gwastro/pycbc/tarball/v1.15.6 +Download-URL: https://github.com/gwastro/pycbc/tarball/v1.16.1 Description: `PyCBC <http://pycbc.org>`_ is a software package used to explore astrophysical sources of gravitational waves. It contains algorithms to analyze gravitational-wave data from the LIGO and Virgo detectors, detect coalescing compact binaries, and measure the astrophysical parameters of detected sources. PyCBC was used in the `first direct detection of gravitational waves <https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.116.061102>`_ and is used in the flagship analysis of LIGO and Virgo data. PyCBC is developed collaboratively and lead by a team of LIGO scientists with the aim to build accessible tools for gravitational-wave data analysis. One of the easiest ways to get a full software environment is by `downloading one of our docker images. <http://pycbc.org/pycbc/latest/html/docker.html>`_ diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/PyCBC.egg-info/SOURCES.txt new/PyCBC-1.16.1/PyCBC.egg-info/SOURCES.txt --- old/PyCBC-1.15.6/PyCBC.egg-info/SOURCES.txt 2020-04-17 11:07:50.000000000 +0200 +++ new/PyCBC-1.16.1/PyCBC.egg-info/SOURCES.txt 2020-04-28 20:25:44.000000000 +0200 @@ -220,6 +220,9 @@ examples/inference/gw150914/plot.sh examples/inference/gw150914/run.sh examples/inference/gw150914/run_test.sh +examples/inference/relative/get.sh +examples/inference/relative/plot.sh +examples/inference/relative/run.sh examples/inference/samplers/run.sh examples/inference/single/get.sh examples/inference/single/plot.sh diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/bin/hdfcoinc/pycbc_page_coinc_snrchi new/PyCBC-1.16.1/bin/hdfcoinc/pycbc_page_coinc_snrchi --- old/PyCBC-1.15.6/bin/hdfcoinc/pycbc_page_coinc_snrchi 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/bin/hdfcoinc/pycbc_page_coinc_snrchi 2020-04-28 19:55:41.000000000 +0200 @@ -50,9 +50,15 @@ # Add the background triggers f = h5py.File(args.coinc_statistic_file, 'r') b_tids = {} -# FIXME: rewrite for multi-ifo case -b_tids[f.attrs['detector_1']] = f['background_exc/trigger_id1'] -b_tids[f.attrs['detector_2']] = f['background_exc/trigger_id2'] +old_style_format_flag = 'detector_1' in f.attrs +if old_style_format_flag: + ifos = [f.attrs['detector_1'], f.attrs['detector_2']] + b_tids[f.attrs['detector_1']] = f['background_exc/trigger_id1'] + b_tids[f.attrs['detector_2']] = f['background_exc/trigger_id2'] +else: + ifos = f.attrs['ifos'].split(' ') + for tmpifo in ifos: + b_tids[tmpifo] = f['background_exc/{}/trigger_id'.format(tmpifo)] f = h5py.File(args.single_trigger_file, 'r') ifo = tuple(f.keys())[0] @@ -74,9 +80,13 @@ # Add the found injection points f = h5py.File(args.found_injection_file, 'r') inj_tids = {} -# FIXME : make compatible with multi-ifo case -inj_tids[f.attrs['detector_1']] = f['found_after_vetoes/trigger_id1'] -inj_tids[f.attrs['detector_2']] = f['found_after_vetoes/trigger_id2'] +if old_style_format_flag: + inj_tids[f.attrs['detector_1']] = f['found_after_vetoes/trigger_id1'] + inj_tids[f.attrs['detector_2']] = f['found_after_vetoes/trigger_id2'] +else: + for tmpifo in ifos: + inj_tids[tmpifo] = f['found_after_vetoes/{}/trigger_id'.format(tmpifo)] + inj_idx = f['found_after_vetoes/injection_index'][:] eff_dist = f['injections'][eff_map[ifo]][:][inj_idx] @@ -165,13 +175,14 @@ pylab.grid(which='major', ls='solid', alpha=0.7, linewidth=.5) pylab.grid(which='minor', ls='solid', alpha=0.7, linewidth=.1) -title = '%s %s chisq vs SNR. Background with injections %s' \ - % (ifo.upper(), args.chisq_choice, +title = '%s %s chisq vs SNR. %s background with injections %s' \ + % (ifo.upper(), args.chisq_choice, ''.join(ifos).upper(), 'behind' if args.background_front else 'ontop') caption = """Distribution of SNR and %s chi-squared veto for single detector -triggers. Black points are background triggers. Triangles are injection +triggers. Black points are %s background triggers. Triangles are injection triggers colored by %s of the injection. Dashed lines show contours of -constant NewSNR.""" % (args.chisq_choice, coloring[args.colorbar_choice][1]) +constant NewSNR.""" % (args.chisq_choice, ''.join(ifos).upper(), + coloring[args.colorbar_choice][1]) pycbc.results.save_fig_with_metadata(fig, args.output_file, title=title, diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/bin/hdfcoinc/pycbc_page_ifar new/PyCBC-1.16.1/bin/hdfcoinc/pycbc_page_ifar --- old/PyCBC-1.15.6/bin/hdfcoinc/pycbc_page_ifar 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/bin/hdfcoinc/pycbc_page_ifar 2020-04-28 19:55:41.000000000 +0200 @@ -24,18 +24,20 @@ import pylab import pycbc.results import pycbc.version +import copy +from pycbc.events import veto from ligo import segments -def calculate_time_slide_duration(ifo1_segments, ifo2_segments, offset=0): +def calculate_time_slide_duration(pifo_segments, fifo_segments, offset=0): ''' Returns the amount of coincident time between two segmentlists. ''' # dertermine how much coincident time - overlap = ifo1_segments.coalesce() & ifo2_segments.coalesce().shift(offset) + overlap = pifo_segments.coalesce() & fifo_segments.coalesce().shift(offset) duration = abs(overlap.coalesce()) # unshift the shifted segmentlist - ifo2_segments.shift(-offset) + fifo_segments.shift(-offset) return duration @@ -89,6 +91,11 @@ # read file fp = h5py.File(opts.trigger_file, 'r') +if 'ifos' in fp.attrs: + legacy_smap_file = False +else: + legacy_smap_file = True + # Parse which inclusive background to use for the plotting h_inc_back_num = opts.use_hierarchical_level @@ -110,17 +117,17 @@ ax.set_ylim(0, 1) output_message = "No more foreground events louder than all background\n" \ - "at this removal level.\n" \ - "Attempted to show " + str(h_inc_back_num) + " removal(s),\n" \ - "but only " + str(h_iterations) + " removal(s) done." + "at this removal level.\nAttempted to show %d " \ + "removal(s),\nbut only %d removal(s) done." % \ + (h_inc_back_num, h_iterations) ax.text(0.5, 0.5, output_message, horizontalalignment='center', verticalalignment='center') pycbc.results.save_fig_with_metadata(fig, opts.output_file, - title='Cumulative Number vs. IFAR', - caption=output_message, - cmd=' '.join(sys.argv)) + title='Cumulative Number vs. IFAR', + caption=output_message, + cmd=' '.join(sys.argv)) # Exit the code successfully and bypass the rest of the plotting code. sys.exit(0) @@ -184,23 +191,36 @@ back_tsid = fp['background/timeslide_id'][:] back_ifar = fp['background/ifar'][:] -# get IFO segments -h1_segments = segments.segmentlist([segments.segment(s,e) for s,e in zip(fp['segments']['H1']['start'][:],fp['segments']['H1']['end'][:])]) -l1_segments = segments.segmentlist([segments.segment(s,e) for s,e in zip(fp['segments']['L1']['start'][:],fp['segments']['L1']['end'][:])]) - # make figure fig = pylab.figure(1) # get a unique list of timeslide_ids and loop over them interval = fp.attrs['timeslide_interval'] -ifo1, ifo2 = fp.attrs['detector_1'], fp.attrs['detector_2'] -min_tsid = (fp['segments'][ifo1]['start'][:].min() - \ - fp['segments'][ifo2]['end'][:].max()) / interval +if legacy_smap_file: + pifo, fifo = fp.attrs['detector_1'], fp.attrs['detector_2'] + p_starts = fp['segments'][pifo]['start'][:] + p_ends = fp['segments'][pifo]['end'][:] + pifo_segments = veto.start_end_to_segments(p_starts, p_ends) + f_starts = fp['segments'][pifo]['start'][:] + f_ends = fp['segments'][pifo]['end'][:] + fifo_segments = veto.start_end_to_segments(f_starts, f_ends) + min_tsid = (p_starts.min() - f_ends.max()) / interval + max_tsid = (p_ends.max() - f_starts.min()) / interval +else: + pifo, fifo = fp.attrs['pivot'], fp.attrs['fixed'] + ifo_joined = fp.attrs['ifos'].replace(' ','') + p_starts = fp['segments'][ifo_joined]['start'][:] + p_ends = fp['segments'][ifo_joined]['end'][:] + pifo_segments = veto.start_end_to_segments(p_starts, p_ends) + fifo_segments = copy.deepcopy(pifo_segments) + min_tsid = (p_starts.min() - p_ends.max()) / interval + max_tsid = -min_tsid + min_tsid = numpy.rint(min_tsid) // opts.decimation_factor -max_tsid = (fp['segments'][ifo1]['end'][:].max() - fp['segments'][ifo2]['start'][:].min()) / interval max_tsid = numpy.rint(max_tsid) // opts.decimation_factor -tsids = numpy.arange(min_tsid, max_tsid, 1).astype(numpy.int64) * opts.decimation_factor +tsids = numpy.arange(min_tsid, max_tsid, 1).astype(numpy.int64) * \ + opts.decimation_factor allbkg_dur = 0 allbkg_ifars = [] @@ -216,9 +236,11 @@ # calculate the amount of coincident time in this time slide offset = tsid*fp.attrs['timeslide_interval'] - back_dur = calculate_time_slide_duration(h1_segments, l1_segments, offset=offset) + back_dur = calculate_time_slide_duration(pifo_segments, fifo_segments, + offset=offset) if back_dur == 0: continue + plotted_slide_count += 1 allbkg_dur = allbkg_dur + back_dur # Limit how far back the decimated background gets plotted @@ -229,7 +251,7 @@ # apply the correction factor for this time slide to its IFAR # you need a correction factor because the analyzed time of the time slide # is not the same as the analyzed time of the foreground - ts_ifar = back_ifar[ts_indx] * ( fp.attrs['foreground_time'] / back_dur ) + ts_ifar = back_ifar[ts_indx] * (fp.attrs['foreground_time'] / back_dur) ts_ifar.sort() # get the cumulative number of triggers in this time slide @@ -251,30 +273,36 @@ label="All decimated background") # plot the expected background -pylab.loglog(expected_ifar, expected_cumnum, linestyle='--', linewidth=2, color='black', label='Expected Background') +pylab.loglog(expected_ifar, expected_cumnum, linestyle='--', linewidth=2, + color='black', label='Expected Background') # plot the counting error error_plus = expected_cumnum + numpy.sqrt(expected_cumnum) error_minus = expected_cumnum - numpy.sqrt(expected_cumnum) error_minus = numpy.where(error_minus<=0, 1e-5, error_minus) -pylab.fill_between(expected_ifar, error_minus, error_plus, facecolor='y', alpha=0.4, label='$N^{1/2}$ Errors') +pylab.fill_between(expected_ifar, error_minus, error_plus, facecolor='y', + alpha=0.4, label='$N^{1/2}$ Errors') # plot the counting error error_plus = expected_cumnum + 2 * numpy.sqrt(expected_cumnum) error_minus = expected_cumnum - 2*numpy.sqrt(expected_cumnum) error_minus = numpy.where(error_minus<=0, 1e-5, error_minus) -pylab.fill_between(expected_ifar, error_minus, error_plus, facecolor='y', alpha=0.2, label='$2N^{1/2}$ Errors') +pylab.fill_between(expected_ifar, error_minus, error_plus, facecolor='y', + alpha=0.2, label='$2N^{1/2}$ Errors') # plot the foreground triggers if opts.open_box: - pylab.loglog(fore_ifar, fore_cumnum, linestyle='None', color='blue', marker='^', label='Foreground') + pylab.loglog(fore_ifar, fore_cumnum, linestyle='None', color='blue', + marker='^', label='Foreground') if h_inc_back_num > 0: - pylab.loglog(h_rm_ifar, h_rm_cumnum, linestyle='None', color='#b66dff', marker='v', label='Hierarchically Removed Foreground') + pylab.loglog(h_rm_ifar, h_rm_cumnum, linestyle='None', color='#b66dff', + marker='v', label='Hierarchically Removed Foreground') # format plot if opts.open_box and len(fore_cumnum) > 100: - # If we have > 100 foreground triggers, scale the plot around the foreground + # If we have > 100 foreground triggers, scale the plot around the + # foreground pylab.ylim(0.8, 1.1 * len(fore_cumnum)) pylab.xlim(0.9 * min(fore_ifar)) elif len(allbkg_cumnum) > 0: @@ -288,13 +316,16 @@ pylab.xlabel('Inverse False Alarm Rate (yr)') # save -caption = 'This is a cumulative histogram of triggers. The blue triangles represent ' \ - + 'coincident foreground triggers. The dashed line represents the expected background ' \ - + 'given the analysis time. The shaded regions represent ' \ - + 'counting errors. The gray lines are time slides treated as zero lag, here there are %d time ' %(plotted_slide_count) \ - + 'slides plotted. Gray dots are time slides with only one event. ' \ - + '%d of the plotted slides have zero events. ' %(empty_slide_count) \ - + 'The green line represents all decimated time slides rescaled to the analysis time.' +caption = 'This is a cumulative histogram of triggers. The blue triangles ' \ + 'represent coincident foreground triggers. The dashed line ' \ + 'represents the expected background given the analysis time. The ' \ + 'shaded regions represent counting errors. The gray lines are ' \ + 'time slides treated as zero lag, here there are %d time slides ' \ + 'plotted. Gray dots are time slides with only one event. %d of ' \ + 'the plotted slides have zero events. The green line represents ' \ + 'all decimated time slides rescaled to the analysis time.' \ + % (plotted_slide_count, empty_slide_count) + pycbc.results.save_fig_with_metadata(fig, opts.output_file, title='Cumulative Number vs. IFAR', caption=caption, diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/bin/plotting/pycbc_ifar_catalog new/PyCBC-1.16.1/bin/plotting/pycbc_ifar_catalog --- old/PyCBC-1.15.6/bin/plotting/pycbc_ifar_catalog 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/bin/plotting/pycbc_ifar_catalog 2020-04-28 19:55:41.000000000 +0200 @@ -145,7 +145,7 @@ # get expected (noise-only) foreground IFAR values and cumulative number for # each IFAR value -expected_ifar = numpy.logspace(-4., numpy.log10(opts.truncate_threshold) or 5, +expected_ifar = numpy.logspace(-4., numpy.log10(opts.truncate_threshold or 10000), num=1000, base=10.0) fg_time = 0 for f in trigf: diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/bin/plotting/pycbc_plot_qscan new/PyCBC-1.16.1/bin/plotting/pycbc_plot_qscan --- old/PyCBC-1.15.6/bin/plotting/pycbc_plot_qscan 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/bin/plotting/pycbc_plot_qscan 2020-04-28 19:55:41.000000000 +0200 @@ -98,13 +98,24 @@ help="If given, use this as the plot caption") parser.add_argument('--colormap', choices=plt.colormaps(), default='viridis', help='Colormap to use (default %(default)s)') +parser.add_argument('--mass1', type=float, + help='Provide masses to plot the waveform track on top of ' + 'the qscan.') +parser.add_argument('--mass2', type=float, + help='Provide masses to plot the waveform track on top of ' + 'the qscan.') +parser.add_argument('--spin1z', type=float, default=0, + help='If masses provided, provide spins to plot the ' + 'waveform track on top of the qscan (default=0).') +parser.add_argument('--spin2z', type=float, default=0, + help='If masses provided, provide spins to plot the ' + 'waveform track on top of the qscan (default=0).') pycbc.strain.insert_strain_option_group(parser) opts = parser.parse_args() logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO) - if opts.center_time is None: center_time = (opts.gps_start_time + opts.gps_end_time) / 2. else: @@ -200,6 +211,17 @@ cb = fig.colorbar(im, ax=(axes.ravel().tolist() if len(wins) > 1 else axes)) cb.set_label('Normalized power') +if opts.mass1: + if not opts.mass2: + raise ValueError('mass1 provided but mass2 missing') + from pycbc.pnutils import get_inspiral_tf + f_low = 20. + approximant = 'SPAtmplt' if opts.mass1+opts.mass2<4 else 'SEOBNRv4_ROM' + track_t, track_f = get_inspiral_tf(opts.center_time, opts.mass1, + opts.mass2, opts.spin1z, opts.spin2z, + f_low, approximant=approximant) + ax.plot(track_t - opts.center_time, track_f, 'r-', lw=1.5) + if opts.plot_title is None: opts.plot_title = 'Q-transform plot around {:.3f}'.format(opts.center_time) if opts.plot_caption is None: diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/bin/workflows/pycbc_create_offline_search_workflow new/PyCBC-1.16.1/bin/workflows/pycbc_create_offline_search_workflow --- old/PyCBC-1.15.6/bin/workflows/pycbc_create_offline_search_workflow 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/bin/workflows/pycbc_create_offline_search_workflow 2020-04-28 19:55:41.000000000 +0200 @@ -424,7 +424,7 @@ table = wf.make_foreground_table(workflow, combined_bg_file, hdfbank, rdir['open_box_result'], singles=insps, extension='.html', - tags=["html"]) + tags=combined_bg_file.tags) #symlink_result(snrifar, 'open_box_result/significance') #symlink_result(ratehist, 'open_box_result/significance') @@ -473,21 +473,18 @@ snrifar_ifar = wf.make_snrifar_plot(workflow, bg_file, open_dir, cumulative=False, tags=bg_file.tags + ['ifar']) - #ifar_ob = wf.make_ifar_plot(workflow, bg_file, open_dir, - # tags=bg_file.tags + ['open_box']) - #ifar_cb = wf.make_ifar_plot(workflow, bg_file, closed_dir, - # tags=bg_file.tags + ['closed_box']) - #table = wf.make_foreground_table(workflow, bg_file, hdfbank, open_dir, - # singles=insps, extension='.html', - # tags=["html"]) - - #detailed_page = [(snrifar, ratehist), (snrifar_ifar, ifar_ob), (table,)] - # FIXME: NEED TO ADD NOT WORKING PLOTS, AS ABOVE - detailed_page = [(snrifar, ratehist), (snrifar_ifar,)] + ifar_ob = wf.make_ifar_plot(workflow, bg_file, open_dir, + tags=bg_file.tags + ['open_box']) + ifar_cb = wf.make_ifar_plot(workflow, bg_file, closed_dir, + tags=bg_file.tags + ['closed_box']) + table = wf.make_foreground_table(workflow, bg_file, hdfbank, open_dir, + singles=insps, extension='.html', + tags=bg_file.tags) + + detailed_page = [(snrifar, ratehist), (snrifar_ifar, ifar_ob), (table,)] layout.two_column_layout(open_dir, detailed_page) - #closed_page = [(snrifar_cb, ifar_cb)] - closed_page = [(snrifar_cb,)] + closed_page = [(snrifar_cb, ifar_cb)] snrifar_summ += closed_page layout.two_column_layout(closed_dir, closed_page) @@ -651,9 +648,9 @@ for inj_insp, trig_insp in zip(inspcomb, fdinspcmb): final_bg_file = final_bg_files[ordered_ifo_list] curr_tags = [tag, ordered_ifo_list] - #f += wf.make_coinc_snrchi_plot(workflow, found_inj, inj_insp, - # final_bg_file, trig_insp, - # injdir, tags=curr_tags) + f += wf.make_coinc_snrchi_plot(workflow, found_inj, inj_insp, + final_bg_file, trig_insp, + injdir, tags=curr_tags) # Make pages from plots inj_layout = list(layout.grouper(f, 2)) + [(found_table,), (missed_table,)] diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/examples/inference/relative/get.sh new/PyCBC-1.16.1/examples/inference/relative/get.sh --- old/PyCBC-1.15.6/examples/inference/relative/get.sh 1970-01-01 01:00:00.000000000 +0100 +++ new/PyCBC-1.16.1/examples/inference/relative/get.sh 2020-04-28 19:55:41.000000000 +0200 @@ -0,0 +1,3 @@ +wget -nc https://dcc.ligo.org/public/0146/P1700349/001/H-H1_LOSC_CLN_4_V1-1187007040-2048.gwf +wget -nc https://dcc.ligo.org/public/0146/P1700349/001/L-L1_LOSC_CLN_4_V1-1187007040-2048.gwf +wget -nc https://dcc.ligo.org/public/0146/P1700349/001/V-V1_LOSC_CLN_4_V1-1187007040-2048.gwf diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/examples/inference/relative/plot.sh new/PyCBC-1.16.1/examples/inference/relative/plot.sh --- old/PyCBC-1.15.6/examples/inference/relative/plot.sh 1970-01-01 01:00:00.000000000 +0100 +++ new/PyCBC-1.16.1/examples/inference/relative/plot.sh 2020-04-28 19:55:41.000000000 +0200 @@ -0,0 +1,4 @@ +pycbc_inference_plot_posterior \ +--input-file relative.hdf \ +--output-file relative.png \ +--z-arg snr diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/examples/inference/relative/run.sh new/PyCBC-1.16.1/examples/inference/relative/run.sh --- old/PyCBC-1.15.6/examples/inference/relative/run.sh 1970-01-01 01:00:00.000000000 +0100 +++ new/PyCBC-1.16.1/examples/inference/relative/run.sh 2020-04-28 19:55:41.000000000 +0200 @@ -0,0 +1,7 @@ +pycbc_inference \ +--config-file `dirname "$0"`/relative.ini \ +--nprocesses=4 \ +--output-file relative.hdf \ +--seed 0 \ +--force \ +--verbose diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/pycbc/events/stat.py new/PyCBC-1.16.1/pycbc/events/stat.py --- old/PyCBC-1.15.6/pycbc/events/stat.py 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/pycbc/events/stat.py 2020-04-28 19:55:41.000000000 +0200 @@ -274,14 +274,19 @@ # Assign attribute so that it can be replaced with other functions self.get_newsnr = ranking.get_newsnr - self.hist = self.hist_ifos = None + self.has_hist = False + self.hist_ifos = None self.ref_snr = 5.0 self.relsense = {} self.swidth = self.pwidth = self.twidth = None + self.srbmin = self.srbmax = None self.max_penalty = None self.pdtype = [] self.weights = {} self.param_bin = {} + self.two_det_flag = (len(ifos) == 2) + self.two_det_weights = {} + self.pb_int_size = None def get_hist(self, ifos=None): """Read in a signal density file for the ifo combination""" @@ -310,16 +315,30 @@ logging.info("Using signal histogram %s for ifos %s", selected, ifos) histfile = self.files[selected] - - # This order matters, we need to retrieve the order used to - # generate the histogram as the first ifos if the reference self.hist_ifos = histfile.attrs['ifos'] + n_ifos = len(self.hist_ifos) + + # Histogram bin attributes + self.twidth = histfile.attrs['twidth'] + self.pwidth = histfile.attrs['pwidth'] + self.swidth = histfile.attrs['swidth'] + self.srbmin = histfile.attrs['srbmin'] + self.srbmax = histfile.attrs['srbmax'] + bin_volume = (self.twidth * self.pwidth * self.swidth) ** (n_ifos - 1) + + # Read histogram for each ifo, to use if that ifo has smallest SNR in + # the coinc for ifo in self.hist_ifos: - self.weights[ifo] = histfile[ifo]['weights'][:] + + weights = histfile[ifo]['weights'][:] + # renormalise to PDF + self.weights[ifo] = weights / (weights.sum() * bin_volume) + param = histfile[ifo]['param_bin'][:] + if param.dtype == numpy.int8: - # Old style, incorrectly sorted histogram file + # Older style, incorrectly sorted histogram file ncol = param.shape[1] self.pdtype = [('c%s' % i, param.dtype) for i in range(ncol)] self.param_bin[ifo] = numpy.zeros(len(self.weights[ifo]), @@ -335,19 +354,47 @@ # param bin and weights have already been sorted self.param_bin[ifo] = param self.pdtype = self.param_bin[ifo].dtype - self.max_penalty = self.weights[ifo].min() - self.hist = {} + # Max_penalty is a small number to assigned to any bins without + # histogram entries. All histograms in a given file have the same + # min entry by design, so use the min of the last one read in. + self.max_penalty = self.weights[ifo].min() - # Bin boundaries are stored in the hdf file - self.twidth = histfile.attrs['twidth'] - self.pwidth = histfile.attrs['pwidth'] - self.swidth = histfile.attrs['swidth'] + if self.two_det_flag: + # The density of signals is computed as a function of 3 binned + # parameters: time difference (t), phase difference (p) and + # SNR ratio (s). These are computed for each combination of + # detectors, so for detectors 6 differences are needed. However + # many combinations of these parameters are highly unlikely and + # no instances of these combinations occurred when generating + # the statistic files. Rather than storing a bunch of 0s, these + # values are just not stored at all. This reduces the size of + # the statistic file, but means we have to identify the correct + # value to read for every trigger. For 2 detectors we can + # expand the weights lookup table here, basically adding in all + # the "0" values. This makes looking up a value in the + # "weights" table a O(N) rather than O(NlogN) operation. It + # sacrifices RAM to do this, so is a good tradeoff for 2 + # detectors, but not for 3! + pb_iinfo = numpy.iinfo(self.param_bin[ifo]['c0'].dtype) + self.pb_int_size = pb_iinfo.max - pb_iinfo.min + 1 + + array_size = [self.pb_int_size, self.pb_int_size, + self.pb_int_size] + dtypec = self.weights[ifo].dtype + self.two_det_weights[ifo] = \ + numpy.zeros(array_size, dtype=dtypec) + self.max_penalty + id0 = self.param_bin[ifo]['c0'] + self.pb_int_size // 2 + id1 = self.param_bin[ifo]['c1'] + self.pb_int_size // 2 + id2 = self.param_bin[ifo]['c2'] + self.pb_int_size // 2 + self.two_det_weights[ifo][id0, id1, id2] = self.weights[ifo] relfac = histfile.attrs['sensitivity_ratios'] for ifo, sense in zip(self.hist_ifos, relfac): self.relsense[ifo] = sense + self.has_hist = True + def single(self, trigs): """Calculate the single detector statistic & assemble other parameters @@ -377,21 +424,35 @@ return self.logsignalrate_multiifo(stats, shift, to_shift) def logsignalrate_multiifo(self, stats, shift, to_shift): - """Calculate the normalized log rate density of signals via lookup""" - # Convert to dict as hist ifos and self.ifos may not be in same - # order + """Calculate the normalized log rate density of signals via lookup + + Parameters + ---------- + stats: list of dicts giving single-ifo quantities, ordered as + self.ifos + shift: numpy array of float, size of the time shift vector for each + coinc to be ranked + to_shift: list of int, multiple of the time shift to apply ordered + as self.ifos + + Returns + ------- + value: log of coinc signal rate density for the given single-ifo + triggers and time shifts + """ + # Convert time shift vector to dict, as hist ifos and self.ifos may + # not be in same order to_shift = {ifo: s for ifo, s in zip(self.ifos, to_shift)} - # does not require ifos to be specified, only 1 p/t/a file - if self.hist is None: + if not self.has_hist: self.get_hist() - # figure out which ifo has the smallest SNR of the contributing ifos - # Store a list 'rtypes' by ifo of which triggers that reference ifo - # should handle. This corresponds to which histogram will be used + # Figure out which ifo of the contributing ifos has the smallest SNR, + # to use as reference for choosing the signal histogram. snrs = numpy.array([numpy.array(stats[ifo]['snr'], ndmin=1) for ifo in self.ifos]) smin = numpy.argmin(snrs, axis=0) + # Store a list of the triggers using each ifo as reference rtypes = {ifo: numpy.where(smin == j)[0] for j, ifo in enumerate(self.ifos)} @@ -431,20 +492,30 @@ sbin = (sdif / self.swidth).astype(numpy.int) binned += [tbin, pbin, sbin] - # convert binned to same dtype as stored in hist + # Convert binned to same dtype as stored in hist nbinned = numpy.zeros(len(pbin), dtype=self.pdtype) for i, b in enumerate(binned): nbinned['c%s' % i] = b # Read signal weight from precalculated histogram - loc = numpy.searchsorted(self.param_bin[ref_ifo], nbinned) - loc[loc == len(self.weights[ref_ifo])] = 0 - rate[rtype] = self.weights[ref_ifo][loc] - - # These weren't in our histogram so give them max penalty instead - # of random value - missed = numpy.where(self.param_bin[ref_ifo][loc] != nbinned)[0] - rate[rtype][missed] = self.max_penalty + if self.two_det_flag: + # High-RAM, low-CPU option for two-det + id0 = nbinned['c0'] + self.pb_int_size // 2 + id1 = nbinned['c1'] + self.pb_int_size // 2 + id2 = nbinned['c2'] + self.pb_int_size // 2 + rate[rtype] = self.two_det_weights[ref_ifo][id0, id1, id2] + else: + # Low[er]-RAM, high[er]-CPU option for >two det + loc = numpy.searchsorted(self.param_bin[ref_ifo], nbinned) + loc[loc == len(self.weights[ref_ifo])] = 0 + rate[rtype] = self.weights[ref_ifo][loc] + + # These weren't in our histogram so give them max penalty + # instead of random value + missed = numpy.where( + self.param_bin[ref_ifo][loc] != nbinned + )[0] + rate[rtype][missed] = self.max_penalty # Scale by signal population SNR rate[rtype] *= (sref / self.ref_snr) ** -4.0 @@ -1091,7 +1162,7 @@ return loglr -class ExpFitSGFgBgRateNewStatistic(PhaseTDNewStatistic, +class ExpFitSGFgBgNormNewStatistic(PhaseTDNewStatistic, ExpFitSGBgRateStatistic): def __init__(self, files=None, ifos=None, **kwargs): @@ -1138,7 +1209,7 @@ tnum = trigs['template_id'] # exists for SingleDetTriggers # Should only be one ifo fit file provided assert len(self.ifos) == 1 - # store benchmark log volume as single-ifo information since the coinc + # Store benchmark log volume as single-ifo information since the coinc # method does not have access to template id singles['benchmark_logvol'] = self.benchmark_logvol[tnum] return numpy.array(singles, ndmin=1) @@ -1168,32 +1239,56 @@ axis=0) # Volume \propto sigma^3 or sigmasq^1.5 network_logvol = 1.5 * numpy.log(network_sigmasq) - # Get benchmark log volume as single-ifo information - # NB benchmark logvol for a given template is not ifo-dependent - # - choose the first ifo for convenience + # Get benchmark log volume as single-ifo information : + # benchmark_logvol for a given template is not ifo-dependent, so + # choose the first ifo for convenience benchmark_logvol = s[0][1]['benchmark_logvol'] network_logvol -= benchmark_logvol + # Use prior histogram to get Bayes factor for signal vs noise + # given the time, phase and SNR differences between IFOs + + # First get signal PDF logr_s stat = {ifo: st for ifo, st in s} logr_s = self.logsignalrate_multiifo(stat, slide * step, to_shift) - loglr = logr_s + network_logvol - ln_noise_rate + # Find total volume of phase-time-amplitude space occupied by noise + # coincs + # Extent of time-difference space occupied + noise_twindow = coinc_rate.multiifo_noise_coincident_area( + self.hist_ifos, kwargs['time_addition']) + # Volume is the allowed time difference window, multiplied by 2pi for + # each phase difference dimension and by allowed range of SNR ratio + # for each SNR ratio dimension : there are (n_ifos - 1) dimensions + # for both phase and SNR + n_ifos = len(self.hist_ifos) + hist_vol = noise_twindow * \ + (2 * numpy.pi * (self.srbmax - self.srbmin) * self.swidth) ** \ + (n_ifos - 1) + # Noise PDF is 1/volume, assuming a uniform distribution of noise + # coincs + logr_n = - numpy.log(hist_vol) + + # Combine to get final statistic: log of + # ((rate of signals / rate of noise) * PTA Bayes factor) + loglr = network_logvol - ln_noise_rate + logr_s - logr_n + # cut off underflowing and very small values loglr[loglr < -30.] = -30. return loglr -class TwoOGCStatistic(ExpFitSGFgBgRateNewStatistic): +class ExpFitSGPSDFgBgNormStatistic(ExpFitSGFgBgNormNewStatistic): def __init__(self, files=None, ifos=None, **kwargs): - ExpFitSGFgBgRateNewStatistic.__init__(self, files=files, ifos=ifos, + ExpFitSGFgBgNormNewStatistic.__init__(self, files=files, ifos=ifos, **kwargs) self.get_newsnr = ranking.get_newsnr_sgveto_psdvar_scaled -class TwoOGCBBHStatistic(ExpFitSGFgBgRateNewStatistic): +class ExpFitSGPSDFgBgNormBBHStatistic(ExpFitSGFgBgNormNewStatistic): def __init__(self, files=None, ifos=None, max_chirp_mass=None, **kwargs): - ExpFitSGFgBgRateNewStatistic.__init__(self, files=files, ifos=ifos, + ExpFitSGFgBgNormNewStatistic.__init__(self, files=files, ifos=ifos, **kwargs) self.get_newsnr = ranking.get_newsnr_sgveto_psdvar_scaled_threshold self.mcm = max_chirp_mass @@ -1206,12 +1301,12 @@ if self.mcm is not None: # Careful - input might be a str, so cast to float self.curr_mchirp = min(self.curr_mchirp, float(self.mcm)) - return ExpFitSGFgBgRateNewStatistic.single(self, trigs) + return ExpFitSGFgBgNormNewStatistic.single(self, trigs) def logsignalrate_multiifo(self, stats, shift, to_shift): # model signal rate as uniform over chirp mass, background rate is # proportional to mchirp^(-11/3) due to density of templates - logr_s = ExpFitSGFgBgRateNewStatistic.logsignalrate_multiifo( + logr_s = ExpFitSGFgBgNormNewStatistic.logsignalrate_multiifo( self, stats, shift, to_shift) logr_s += numpy.log((self.curr_mchirp / 20.0) ** (11./3.0)) return logr_s @@ -1238,9 +1333,11 @@ PhaseTDExpFitSGPSDScaledStatistic, 'exp_fit_sg_bg_rate': ExpFitSGBgRateStatistic, 'exp_fit_sg_fgbg_rate': ExpFitSGFgBgRateStatistic, - 'exp_fit_sg_fgbg_rate_new': ExpFitSGFgBgRateNewStatistic, - '2ogc': TwoOGCStatistic, - '2ogcbbh': TwoOGCBBHStatistic, + 'exp_fit_sg_fgbg_norm_new': ExpFitSGFgBgNormNewStatistic, + '2ogc': ExpFitSGPSDFgBgNormStatistic, # backwards compatible + '2ogcbbh': ExpFitSGPSDFgBgNormBBHStatistic, # backwards compatible + 'exp_fit_sg_fgbg_norm_psdvar': ExpFitSGPSDFgBgNormStatistic, + 'exp_fit_sg_fgbg_norm_psdvar_bbh': ExpFitSGPSDFgBgNormBBHStatistic } sngl_statistic_dict = { diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/pycbc/filter/matchedfilter.py new/PyCBC-1.16.1/pycbc/filter/matchedfilter.py --- old/PyCBC-1.15.6/pycbc/filter/matchedfilter.py 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/pycbc/filter/matchedfilter.py 2020-04-28 19:55:41.000000000 +0200 @@ -1088,7 +1088,7 @@ if psd: try: numpy.testing.assert_almost_equal(ht.delta_f, psd.delta_f) - except: + except AssertionError: raise ValueError('Waveform does not have same delta_f as psd') if psd is None: @@ -1234,10 +1234,11 @@ if psd is not None: if isinstance(psd, FrequencySeries): - if psd.delta_f == stilde.delta_f : - qtilde[kmin:kmax] /= psd[kmin:kmax] - else: - raise TypeError("PSD delta_f does not match data") + try: + numpy.testing.assert_almost_equal(stilde.delta_f, psd.delta_f) + except AssertionError: + raise ValueError("PSD delta_f does not match data") + qtilde[kmin:kmax] /= psd[kmin:kmax] else: raise TypeError("PSD must be a FrequencySeries") @@ -1902,4 +1903,3 @@ 'compute_u_val_for_sky_loc_stat_no_phase', 'compute_u_val_for_sky_loc_stat', 'followup_event_significance'] - diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/pycbc/inference/models/__init__.py new/PyCBC-1.16.1/pycbc/inference/models/__init__.py --- old/PyCBC-1.15.6/pycbc/inference/models/__init__.py 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/pycbc/inference/models/__init__.py 2020-04-28 19:55:41.000000000 +0200 @@ -25,7 +25,7 @@ from .gaussian_noise import GaussianNoise from .marginalized_gaussian_noise import MarginalizedPhaseGaussianNoise from .single_template import SingleTemplate -from .relbin import RelativeSPA +from .relbin import Relative # Used to manage a model instance across multiple cores or MPI @@ -183,5 +183,5 @@ GaussianNoise, MarginalizedPhaseGaussianNoise, SingleTemplate, - RelativeSPA + Relative )} diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/pycbc/inference/models/relbin.py new/PyCBC-1.16.1/pycbc/inference/models/relbin.py --- old/PyCBC-1.15.6/pycbc/inference/models/relbin.py 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/pycbc/inference/models/relbin.py 2020-04-28 19:55:41.000000000 +0200 @@ -31,8 +31,9 @@ from scipy.interpolate import interp1d from scipy import special -from pycbc.waveform.spa_tmplt import spa_tmplt +from pycbc.waveform import get_fd_waveform_sequence from pycbc.detector import Detector +from pycbc.types import Array from .gaussian_noise import BaseGaussianNoise @@ -81,14 +82,15 @@ # frequency grid points fbin = dphi2f(dphi_grid) # indices of frequency grid points in the FFT array - fbin_ind = numpy.array([numpy.argmin(numpy.absolute(f_full - ff)) for - ff in fbin]) + fbin_ind = numpy.unique([numpy.argmin(numpy.absolute(f_full - ff)) for + ff in fbin]) # make sure grid points are precise fbin = numpy.array([f_full[i] for i in fbin_ind]) + nbin = len(fbin) return nbin, fbin, fbin_ind -class RelativeSPA(BaseGaussianNoise): +class Relative(BaseGaussianNoise): r"""Model that assumes the likelihood in a region around the peak is slowly varying such that a linear approximation can be made, and likelihoods can be calculated at a coarser frequency resolution. For @@ -140,14 +142,11 @@ All other keyword arguments are passed to :py:class:`BaseGaussianNoise`. """ - name = "relative_spa" + name = "relative" def __init__(self, variable_params, data, low_frequency_cutoff, - mass1_ref, mass2_ref, spin1z_ref, spin2z_ref, - ra_ref, dec_ref, tc_ref, - epsilon=0.5, - **kwargs): - super(RelativeSPA, self).__init__( + fiducial_params=None, epsilon=0.5, **kwargs): + super(Relative, self).__init__( variable_params, data, low_frequency_cutoff, **kwargs) # check that all of the frequency cutoffs are the same # FIXME: this can probably be loosened at some point @@ -168,32 +167,39 @@ self.comp_data = {ifo: d.numpy() for ifo, d in self.data.items()} self.comp_psds = {ifo: p.numpy() for ifo, p in self.psds.items()} # store fiducial waveform params - self.mass1_ref = float(mass1_ref) - self.mass2_ref = float(mass2_ref) - self.spin1z_ref = float(spin1z_ref) - self.spin2z_ref = float(spin2z_ref) - self.ra_ref = float(ra_ref) - self.dec_ref = float(dec_ref) - self.tc_ref = float(tc_ref) + self.fid_params = fiducial_params # get detector-specific arrival times relative to end of data dt = {ifo: - self.det[ifo].time_delay_from_earth_center(self.ra_ref, - self.dec_ref, - self.tc_ref) + self.det[ifo].time_delay_from_earth_center( + self.fid_params['ra'], self.fid_params['dec'], + self.fid_params['tc']) for ifo in self.data} - self.ta = {ifo: self.tc_ref + dt[ifo] - self.end_time + self.ta = {ifo: self.fid_params['tc'] + dt[ifo] - self.end_time for ifo in self.data} # generate fiducial waveform - logging.info("Generating fiducial waveform") - hp = spa_tmplt(f_lower=kmins[0]*self.df, f_upper=(kmaxs[0]+1)*self.df, - delta_f=self.df, mass1=self.mass1_ref, - mass2=self.mass2_ref, spin1z=self.spin1z_ref, - spin2z=self.spin2z_ref, distance=1., - spin_order=-1, phase_order=-1) - hp.resize(len(self.f)) - self.h00 = numpy.array(hp) + f_lo = kmins[0] * self.df + f_hi = kmaxs[0] * self.df + logging.info("Generating fiducial waveform from %s to %s Hz", + f_lo, f_hi) + # prune low frequency samples to avoid waveform errors + nbelow = sum(self.f < 10) + fpoints = Array(self.f.astype(numpy.float64))[nbelow:] + approx = self.static_params['approximant'] + fid_hp, fid_hc = get_fd_waveform_sequence(approximant=approx, + sample_points=fpoints, + **self.fid_params) + self.h00 = {} + for ifo in self.data: + # make copy of fiducial wfs, adding back in low frequencies + hp0 = numpy.concatenate([[0j] * nbelow, fid_hp.copy()]) + hc0 = numpy.concatenate([[0j] * nbelow, fid_hc.copy()]) + fp, fc = self.det[ifo].antenna_pattern( + self.fid_params['ra'], self.fid_params['dec'], + self.fid_params['polarization'], self.fid_params['tc']) + tshift = numpy.exp(-2.0j * numpy.pi * self.f * self.ta[ifo]) + self.h00[ifo] = numpy.array(hp0 * fp + hc0 * fc) * tshift # compute frequency bins logging.info("Computing frequency bins") @@ -203,13 +209,14 @@ logging.info("Using %s bins for this model", nbin) # store bins and edges in sample and frequency space self.edges = fbin_ind - self.fedges = numpy.array(fbin).astype(numpy.float32) + self.fedges = numpy.array(fbin).astype(numpy.float64) self.bins = numpy.array([(self.edges[i], self.edges[i+1]) for i in range(len(self.edges) - 1)]) self.fbins = numpy.array([(fbin[i], fbin[i+1]) for i in range(len(fbin) - 1)]) # store low res copy of fiducial waveform - self.h00_sparse = self.h00.copy().take(self.edges) + self.h00_sparse = {ifo: self.h00[ifo].copy().take(self.edges) for ifo + in self.h00} # compute summary data logging.info("Calculating summary data at frequency resolution %s Hz", @@ -227,53 +234,30 @@ containing bin coefficients a0, b0, a1, b1, for each frequency bin. """ - # timeshift the fiducial waveform for each detector - shift = {ifo: numpy.exp(-2.0j * numpy.pi * self.f * self.ta[ifo]) for - ifo in self.data} - h0 = {ifo: self.h00.copy() * shift[ifo] for ifo in self.data} # calculate coefficients sdat = {} for ifo in self.data: - hd = numpy.conjugate(self.comp_data[ifo]) * h0[ifo] + hd = numpy.conjugate(self.comp_data[ifo]) * self.h00[ifo] hd /= self.comp_psds[ifo] - hh = (numpy.absolute(h0[ifo]) ** 2.0) / self.comp_psds[ifo] + hh = (numpy.absolute(self.h00[ifo]) ** 2.0) / self.comp_psds[ifo] # constant terms a0 = numpy.array([4. * self.df * numpy.sum(hd[l:h]) for l, h in self.bins]) b0 = numpy.array([4. * self.df * numpy.sum(hh[l:h]) for l, h in self.bins]) # linear terms - bin_centers = [0.5 * (fl + fh) for fl, fh in self.fbins] + bin_lefts = [fl for fl, fh in self.fbins] a1 = numpy.array([4. * self.df - * numpy.sum(hd[l:h] * (self.f[l:h] - bc)) for - (l, h), bc in zip(self.bins, bin_centers)]) + * numpy.sum(hd[l:h] * (self.f[l:h] - bl)) for + (l, h), bl in zip(self.bins, bin_lefts)]) b1 = numpy.array([4. * self.df - * numpy.sum(hh[l:h] * (self.f[l:h] - bc)) for - (l, h), bc in zip(self.bins, bin_centers)]) + * numpy.sum(hh[l:h] * (self.f[l:h] - bl)) for + (l, h), bl in zip(self.bins, bin_lefts)]) sdat[ifo] = {'a0': a0, 'a1': a1, 'b0': b0, 'b1': b1} return sdat - def waveform_ratio(self, p, htf, dtc=0.0): - """Calculate waveform ratio between template and fiducial - waveforms. - """ - # generate template - hp = spa_tmplt(sample_points=self.fedges, - mass1=p['mass1'], mass2=p['mass2'], - spin1z=p['spin1z'], spin2z=p['spin2z'], - distance=1., spin_order=-1, phase_order=-1) - htarget = numpy.array(hp) - # apply antenna pattern, inclination, and distance factor - htarget *= htf - # compute waveform ratio and timeshift - shift = numpy.exp(-2.0j * numpy.pi * self.fedges * dtc) - r = htarget / self.h00_sparse * shift - r0 = 0.5 * (r[:-1] + r[1:]) - r1 = (r[1:] - r[:-1]) / (self.fedges[1:] - self.fedges[:-1]) - return numpy.array([r0, r1], dtype=numpy.complex128) - def _loglr(self): r"""Computes the log likelihood ratio, @@ -294,32 +278,35 @@ p = self.current_params.copy() p.update(self.static_params) - llr = 0. + hh = 0. + hd = 0j for ifo in self.data: # get detector antenna pattern fp, fc = self.det[ifo].antenna_pattern(p['ra'], p['dec'], p['polarization'], p['tc']) - ip = numpy.cos(p['inclination']) - ic = 0.5 * (1.0 + ip * ip) - htf = (fp * ip + 1.0j * fc * ic) / p['distance'] - # get timeshift relative to fiducial waveform + # get timeshift relative to end of data dt = self.det[ifo].time_delay_from_earth_center(p['ra'], p['dec'], p['tc']) - dtc = p['tc'] + dt - self.end_time - self.ta[ifo] + dtc = p['tc'] + dt - self.end_time + tshift = numpy.exp(-2.0j * numpy.pi * self.fedges * dtc) # generate template and calculate waveform ratio - r0, r1 = self.waveform_ratio(p, htf, dtc=dtc) - # <h, d> is real part of sum over bins of A0r0 + A1r1 - hd = numpy.sum(self.sdat[ifo]['a0'] * r0 - + self.sdat[ifo]['a1'] * r1).real - # marginalize over phase - hd = numpy.log(special.i0e(hd)) + abs(hd) - # <h, h> is real part of sum over bins of B0|r0|^2 + 2B1Re(r1r0*) - hh = numpy.sum(self.sdat[ifo]['b0'] * numpy.absolute(r0) ** 2. - + 2. * self.sdat[ifo]['b1'] - * (r1 * numpy.conjugate(r0)).real).real - # increment loglr - llr += (hd - 0.5 * hh) + hp, hc = get_fd_waveform_sequence(sample_points=Array(self.fedges), + **p) + htilde = numpy.array(fp * hp + fc * hc) * tshift + r = (htilde / self.h00_sparse[ifo]).astype(numpy.complex128) + r0 = r[:-1] + r1 = (r[1:] - r[:-1]) / (self.fedges[1:] - self.fedges[:-1]) + + # <h, d> is sum over bins of A0r0 + A1r1 + hd += numpy.sum(self.sdat[ifo]['a0'] * r0 + + self.sdat[ifo]['a1'] * r1) + # <h, h> is sum over bins of B0|r0|^2 + 2B1Re(r1r0*) + hh += numpy.sum(self.sdat[ifo]['b0'] * numpy.absolute(r0) ** 2. + + 2. * self.sdat[ifo]['b1'] + * (r1 * numpy.conjugate(r0)).real) + hd = abs(hd) + llr = numpy.log(special.i0e(hd)) + hd - 0.5 * hh return float(llr) def write_metadata(self, fp): @@ -330,12 +317,27 @@ fp : pycbc.inference.io.BaseInferenceFile instance The inference file to write to. """ - super(RelativeSPA, self).write_metadata(fp) - fp.attrs['mass1_ref'] = self.mass1_ref - fp.attrs['mass2_ref'] = self.mass2_ref - fp.attrs['spin1z_ref'] = self.spin1z_ref - fp.attrs['spin2z_ref'] = self.spin2z_ref - fp.attrs['ra_ref'] = self.ra_ref - fp.attrs['dec_ref'] = self.dec_ref - fp.attrs['tc_ref'] = self.tc_ref + super(Relative, self).write_metadata(fp) fp.attrs['epsilon'] = self.epsilon + for p, v in self.fid_params.items(): + fp.attrs['{}_ref'.format(p)] = v + + @staticmethod + def extra_args_from_config(cp, section, skip_args=None, dtypes=None): + """Adds reading fiducial waveform parameters from config file. + """ + # add fiducial params to skip list + skip_args += [option for option in cp.options(section) if + option.endswith('_ref')] + args = super(Relative, Relative).extra_args_from_config( + cp, section, skip_args=skip_args, dtypes=dtypes) + # get fiducial params from config + fid_params = {p.replace('_ref', ''): float(cp.get('model', p)) for p + in cp.options('model') if p.endswith('_ref')} + # add optional params with default values if not specified + opt_params = {'ra': numpy.pi, 'dec': 0.0, 'inclination': 0.0, + 'polarization': numpy.pi} + fid_params.update({p: opt_params[p] for p in opt_params if p + not in fid_params}) + args.update({'fiducial_params': fid_params}) + return args diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/pycbc/version.py new/PyCBC-1.16.1/pycbc/version.py --- old/PyCBC-1.15.6/pycbc/version.py 2020-04-17 11:07:50.000000000 +0200 +++ new/PyCBC-1.16.1/pycbc/version.py 2020-04-28 20:25:44.000000000 +0200 @@ -1,24 +1,24 @@ # coding: utf-8 -# Generated by setup.py for PyCBC on 2020-04-17 09:07:50 +0000. +# Generated by setup.py for PyCBC on 2020-04-28 18:25:44 +0000. -version = '1.15.6' -date = '2020-04-17 08:16:41 +0000' +version = '1.16.1' +date = '2020-04-28 17:24:30 +0000' release = True -last_release = '1.15.6' +last_release = '1.16.1' -git_hash = '693210e77b2d7d24332f345dc4d8f1aa9bc54b99' +git_hash = 'af7a10948ed6d0d82a231bb6a0fe576cb880128a' git_branch = 'None' -git_tag = 'v1.15.6' -git_author = 'Ian Harry <[email protected]>' +git_tag = 'v1.16.1' +git_author = 'Duncan Brown <[email protected]>' git_committer = 'GitHub <[email protected]>' git_status = 'CLEAN: All modifications committed' git_builder = 'Travis CI User <[email protected]>' -git_build_date = '2020-04-17 09:07:50 +0000' -git_verbose_msg = """Version: 1.15.6 +git_build_date = '2020-04-28 18:25:44 +0000' +git_verbose_msg = """Version: 1.16.1 Branch: None -Tag: v1.15.6 -Id: 693210e77b2d7d24332f345dc4d8f1aa9bc54b99 +Tag: v1.16.1 +Id: af7a10948ed6d0d82a231bb6a0fe576cb880128a Builder: Travis CI User <[email protected]> -Build date: 2020-04-17 09:07:50 +0000 +Build date: 2020-04-28 18:25:44 +0000 Repository status is CLEAN: All modifications committed""" from pycbc._version import * diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/setup.py new/PyCBC-1.16.1/setup.py --- old/PyCBC-1.15.6/setup.py 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/setup.py 2020-04-28 19:55:41.000000000 +0200 @@ -124,8 +124,8 @@ vinfo = _version_helper.generate_git_version_info() except: vinfo = vdummy() - vinfo.version = '1.15.6' - vinfo.release = 'False' + vinfo.version = '1.16.1' + vinfo.release = 'True' with open('pycbc/version.py', 'w') as f: f.write("# coding: utf-8\n") diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/PyCBC-1.15.6/test/test_matchedfilter.py new/PyCBC-1.16.1/test/test_matchedfilter.py --- old/PyCBC-1.15.6/test/test_matchedfilter.py 2020-04-17 10:45:44.000000000 +0200 +++ new/PyCBC-1.16.1/test/test_matchedfilter.py 2020-04-28 19:55:41.000000000 +0200 @@ -134,7 +134,7 @@ #Check that an incompatible psd produces an error self.assertRaises(TypeError,match,self.filt,self.filt,psd=self.filt) psd = FrequencySeries(zeros(len(self.filt)/2+1),delta_f=100000) - self.assertRaises(TypeError,match,self.filt,self.filt,psd=psd) + self.assertRaises(ValueError,match,self.filt,self.filt,psd=psd) #Check that only TimeSeries or FrequencySeries are accepted self.assertRaises(TypeError,match,zeros(10),zeros(10))
