add ability to cache detector response to marginalize time model and …

…us… (#4806)

* add ability to cache detect response to marginalize time model and use alternate sampler rate than data

* cc

* update

* update example

* update example

examples/inference/margtime/margtime.ini

@@ -3,8 +3,11 @@
 name = marginalized_time
 low-frequency-cutoff = 30.0
 
+# This is the sample rate used for the model and marginalization
+sample_rate = 4096
+
 marginalize_vector_params = tc, ra, dec, polarization
-marginalize_vector_samples = 500
+marginalize_vector_samples = 2000
 
 ; You shouldn't use phase marginalization if the approximant has
 ; higher-order modes

examples/inference/margtime/run.sh

@@ -1,6 +1,6 @@
 OMP_NUM_THREADS=1 pycbc_inference \
 --config-file `dirname "$0"`/margtime.ini \
---nprocesses 2 \
+--nprocesses 1 \
 --processing-scheme mkl \
 --output-file marg_150914.hdf \
 --seed 0 \
@@ -23,4 +23,5 @@ pycbc_inference_plot_posterior \
  "primary_mass(mass1, mass2) / (1 + redshift(distance)):srcmass1" \
  "secondary_mass(mass1, mass2) / (1 + redshift(distance)):srcmass2" \
  ra dec tc inclination coa_phase polarization distance \
+--vmin 23.2 \
 --z-arg snr

pycbc/inference/models/marginalized_gaussian_noise.py

@@ -19,6 +19,7 @@
 """
 
 import itertools
+import logging
 import numpy
 from scipy import special
 
@@ -207,8 +208,10 @@ class MarginalizedTime(DistMarg, BaseGaussianNoise):
     def __init__(self, variable_params,
                  data, low_frequency_cutoff, psds=None,
                  high_frequency_cutoff=None, normalize=False,
+                 sample_rate=None,
                  **kwargs):
 
+        self.sample_rate = float(sample_rate)
         self.kwargs = kwargs
         variable_params, kwargs = self.setup_marginalization(
                                variable_params,
@@ -241,6 +244,14 @@ def __init__(self, variable_params,
 
         self.dets = {}
 
+        if sample_rate is not None:
+            for ifo in self.data:
+                if self.sample_rate < self.data[ifo].sample_rate:
+                    raise ValueError("Model sample rate was set less than the"
+                                     " data. ")
+            logging.info("Using %s sample rate for marginalization",
+                         sample_rate)
+
     def _nowaveform_loglr(self):
         """Convenience function to set loglr values if no waveform generated.
         """
@@ -296,8 +307,15 @@ def _loglr(self):
             hp[self._kmin[det]:kmax] *= self._weight[det][slc]
             hc[self._kmin[det]:kmax] *= self._weight[det][slc]
 
-            hp.resize(len(self._whitened_data[det]))
-            hc.resize(len(self._whitened_data[det]))
+            # Use a higher sample rate if requested
+            if self.sample_rate is not None:
+                tlen = int(round(self.sample_rate *
+                           self.whitened_data[det].duration))
+                flen = tlen // 2 + 1
+                hp.resize(flen)
+                hc.resize(flen)
+                self._whitened_data[det].resize(flen)
+
             cplx_hpd[det], _, _ = matched_filter_core(
                                  hp,
                                  self._whitened_data[det],
@@ -325,15 +343,20 @@ def _loglr(self):
         for det in wfs:
             if det not in self.dets:
                 self.dets[det] = Detector(det)
-            fp, fc = self.dets[det].antenna_pattern(
-                                    params['ra'],
-                                    params['dec'],
-                                    params['polarization'],
-                                    params['tc'])
-            dt = self.dets[det].time_delay_from_earth_center(params['ra'],
-                                                             params['dec'],
-                                                             params['tc'])
+
+            if self.precalc_antenna_factors:
+                fp, fc, dt = self.get_precalc_antenna_factors(det)
+            else:
+                fp, fc = self.dets[det].antenna_pattern(
+                                        params['ra'],
+                                        params['dec'],
+                                        params['polarization'],
+                                        params['tc'])
+                dt = self.dets[det].time_delay_from_earth_center(params['ra'],
+                                                                 params['dec'],
+                                                                 params['tc'])
             dtc = params['tc'] + dt
+
             cplx_hd = fp * cplx_hpd[det].at_time(dtc,
                                                  interpolate='quadratic')
             cplx_hd += fc * cplx_hcd[det].at_time(dtc,