Save overlapping PSD and read it to use while filtering

bin/pycbc_inspiral

@@ -34,7 +34,6 @@ from pycbc.filter import MatchedFilterControl, make_frequency_series, qtransform
 from pycbc.types import TimeSeries, FrequencySeries, zeros, float32, complex64
 import pycbc.opt
 import pycbc.inject
-from pycbc.io import HFile
 
 last_progress_update = -1.0
 
@@ -205,10 +204,6 @@ parser.add_argument("--multiprocessing-nprocesses", type=int,
                          "Used in conjunction with the option"
                          "--finalize-events-template-rate which should be set"
                          "to a multiple of the number of processes.")
-parser.add_argument("--associated-psd-output",
-                    help="(Optional) Write PSD to specified file in the format "
-                    "of calculate_psd output. --output-psd only saves psd of "
-                    "only one segment")
 
 # Add options groups
 psd.insert_psd_option_group(parser)
@@ -441,37 +436,9 @@ with ctx:
                                  maximal_value_dof=opt.autochi_max_valued_dof)
 
     logging.info("Overwhitening frequency-domain data segments")
-    if hasattr(opt, 'associated_psd_output') and opt.associated_psd_output:
-        logging.info("Saving PSDs for filter segments")
-        f = HFile(opt.associated_psd_output, 'w')
-        ifo = opt.channel_name[0:2]
-        psd_group = f.create_group(ifo + '/psds')
-        start, end = [], []
-
     for inc, seg in enumerate(segments):
         seg /= seg.psd
 
-
-        if hasattr(opt, 'associated_psd_output') and opt.associated_psd_output:
-            key = str(inc)
-            # start_time = gwstrain.start_time + seg.analyze.start/gwstrain.sample_rate
-            # end_time = gwstrain.start_time + seg.analyze.stop/gwstrain.sample_rate
-            # ic(start_time, end_time)
-            start.append(int(seg.start_time))
-            end.append(int(seg.end_time))
-            psd_group.create_dataset(key, data=seg.psd, compression='gzip',
-                                 compression_opts=9, shuffle=True)
-            psd_group[key].attrs['epoch'] = int(seg.start_time)
-            psd_group[key].attrs['delta_f'] = seg.psd.delta_f
-
-    if hasattr(opt, 'associated_psd_output') and opt.associated_psd_output:
-        f[ifo + '/start_time'] = numpy.array(start, dtype=numpy.uint32)
-        f[ifo + '/end_time'] = numpy.array(end, dtype=numpy.uint32)
-        f.attrs['low_frequency_cutoff'] = opt.low_frequency_cutoff
-        f.attrs['dynamic_range_factor'] = pycbc.DYN_RANGE_FAC
-        f.close()
-
-
     logging.info("Read in template bank")
     bank = waveform.FilterBank(opt.bank_file, flen, delta_f,
         low_frequency_cutoff=None if opt.enable_bank_start_frequency else flow,

pycbc/events/ranking.py

@@ -26,7 +26,7 @@ def newsnr(snr, reduced_x2, q=6., n=2.):
     reduced chi-squared values. See http://arxiv.org/abs/1208.3491 for
     definition. Previous implementation in glue/ligolw/lsctables.py
     """
-    nsnr = numpy.array(snr, ndmin=1, dtype=numpy.float64)
+    nsnr = numpy.array(numpy.abs(snr), ndmin=1, dtype=numpy.float64)
     reduced_x2 = numpy.array(reduced_x2, ndmin=1, dtype=numpy.float64)
 
     # newsnr is only different from snr if reduced chisq > 1

pycbc/psd/__init__.py

@@ -15,6 +15,9 @@
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 import copy
+import numpy
+import logging
+from pycbc.io import HFile
 from ligo import segments
 from pycbc.psd.read import *
 from pycbc.psd.analytical import *
@@ -27,6 +30,9 @@
 from pycbc.types import copy_opts_for_single_ifo
 from pycbc.types import required_opts, required_opts_multi_ifo
 from pycbc.types import ensure_one_opt, ensure_one_opt_multi_ifo
+from pycbc import DYN_RANGE_FAC
+
+logger = logging.getLogger('pycbc.psd')
 
 def from_cli(opt, length, delta_f, low_frequency_cutoff,
              strain=None, dyn_range_factor=1, precision=None):
@@ -272,7 +278,10 @@ def insert_psd_option_group(parser, output=True, include_data_options=True):
     if output:
         psd_options.add_argument("--psd-output",
                           help="(Optional) Write PSD to specified file")
-
+        psd_options.add_argument("--overlapping-psd-output",
+                    help="(Optional) Write PSD to specified file in the format "
+                    "of calculate_psd output. --output-psd only saves psd of "
+                    "only one segment")
     return psd_options
 
 def insert_psd_option_group_multi_ifo(parser):
@@ -505,6 +514,32 @@ def generate_overlapping_psds(opt, gwstrain, flen, delta_f, flow,
         psd = from_cli(opt, flen, delta_f, flow, strain=strain_part,
                        dyn_range_factor=dyn_range_factor, precision=precision)
         psds_and_times.append( (start_idx, end_idx, psd) )
+
+    if hasattr(opt, 'overlapping_psd_output') and opt.overlapping_psd_output:
+        logging.info(f"Saving overlapping PSDs for segments to "
+                f"{opt.overlapping_psd_output}")
+        f = HFile(opt.overlapping_psd_output, 'w')
+        ifo = opt.channel_name[0:2]
+        psd_group = f.create_group(ifo + '/psds')
+        start, end = [], []
+        for inc, pnt in enumerate(psds_and_times):
+            start_idx, end_idx, psd = pnt
+            key = str(inc)
+            start_time = gwstrain.start_time + start_idx/gwstrain.sample_rate
+            end_time = gwstrain.start_time + end_idx/gwstrain.sample_rate
+            start.append(int(start_time))
+            end.append(int(end_time))
+            psd_group.create_dataset(key, data=psd.numpy(), compression='gzip',
+                    compression_opts=9, shuffle=True)
+            psd_group[key].attrs['epoch'] = int(start_time)
+            psd_group[key].attrs['delta_f'] = psd.delta_f
+
+        f[ifo + '/start_time'] = numpy.array(start, dtype=numpy.uint32)
+        f[ifo + '/end_time'] = numpy.array(end, dtype=numpy.uint32)
+        f.attrs['low_frequency_cutoff'] = opt.low_frequency_cutoff
+        f.attrs['dynamic_range_factor'] = DYN_RANGE_FAC
+        f.close()
+
     return psds_and_times
 
 def associate_psds_to_segments(opt, fd_segments, gwstrain, flen, delta_f, flow,
@@ -541,6 +576,7 @@ def associate_psds_to_segments(opt, fd_segments, gwstrain, flen, delta_f, flow,
         not already in that precision.
     """
     if opt.precomputed_psd_file:
+        logging.info(f"Reading Precomuted PSDs from {opt.precomputed_psd_file}")
         tpsd = PrecomputedTimeVaryingPSD(opt,
                                          length=len(fd_segments[0].data),
                                          delta_f=fd_segments[0].delta_f,

pycbc/psd/read.py

@@ -230,7 +230,7 @@ def assosiate_psd_to_inspiral_segment(self, inp_seg, delta_f=None):
             sidx = numpy.argpartition(
                 numpy.abs(self.start_times - inp_seg[0]), 2)[:2]
         else:
-            sidx = np.array([0, 1])
+            sidx = numpy.argsort(numpy.abs(self.start_times - inp_seg[0]))
 
         nearest = segments.segment(
             self.start_times[sidx[0]], self.end_times[sidx[0]])