Time-domain likelihood (with heterodyning)

bilby/gw/__init__.py

@@ -1,7 +1,7 @@
 from . import (conversion, cosmology, detector, eos, likelihood, prior,
                result, source, utils, waveform_generator)
 from .waveform_generator import WaveformGenerator, LALCBCWaveformGenerator
-from .likelihood import GravitationalWaveTransient
+from .likelihood import GravitationalWaveTransient, TimeDomainGravitationalWaveTransient
 from .detector import calibration
 from . import compat
 

bilby/gw/conversion.py

@@ -2227,6 +2227,14 @@ def generate_source_frame_parameters(sample):
     return output_sample
 
 
+def _signal_polarizations_for_likelihood(likelihood, parameters):
+    """Return the time-domain waveform for the time-domain likelihood, and
+    the frequency-domain waveform for all other likelihoods."""
+    if getattr(likelihood, "likelihood_domain", "frequency") == "time":
+        return likelihood.waveform_generator.time_domain_strain(parameters)
+    return likelihood.waveform_generator.frequency_domain_strain(parameters)
+
+
 def compute_snrs(sample, likelihood, npool=1, pool=None):
     """
     Compute the optimal and matched filter snrs of all posterior samples
@@ -2242,7 +2250,7 @@ def compute_snrs(sample, likelihood, npool=1, pool=None):
     """
     if likelihood is not None:
         if isinstance(sample, dict):
-            signal_polarizations = likelihood.waveform_generator.frequency_domain_strain(sample.copy())
+            signal_polarizations = _signal_polarizations_for_likelihood(likelihood, sample.copy())
             for ifo in likelihood.interferometers:
                 per_detector_snr = likelihood.calculate_snrs(signal_polarizations, ifo, parameters=sample)
                 sample['{}_matched_filter_snr'.format(ifo.name)] =\
@@ -2280,9 +2288,7 @@ def _compute_snrs(args):
     likelihood = sampling_convenience_dump.likelihood
     ii, sample = args
     sample = dict(sample).copy()
-    signal_polarizations = likelihood.waveform_generator.frequency_domain_strain(
-        sample.copy()
-    )
+    signal_polarizations = _signal_polarizations_for_likelihood(likelihood, sample.copy())
     snrs = list()
     for ifo in likelihood.interferometers:
         snrs.append(likelihood.calculate_snrs(

bilby/gw/detector/__init__.py

@@ -1,4 +1,6 @@
 from ..conversion import convert_to_lal_binary_black_hole_parameters
+from .acf import *
+from .condition_data import *
 from .calibration import *
 from .interferometer import *
 from .networks import *

bilby/gw/detector/acf.py

@@ -0,0 +1,153 @@
+import numpy as np
+
+from ...core.utils import logger
+
+
+class AutoCovarianceFunction:
+
+    def __init__(self, acf_array=None, sampling_frequency=None, minimum_frequency=None,
+                 maximum_frequency=None, fill_value=None, fill_multiplier=None):
+        """
+        Instantiate a new AutoCovarianceFunction object.
+
+        The autocovariance function (ACF) of stationary Gaussian noise gives the
+        noise covariance matrix of N consecutive samples, which is a Toeplitz
+        matrix with first row equal to the first N lags of the ACF.
+
+        Examples
+        ========
+        From a one-sided power spectral density:
+
+        .. code-block:: python
+
+            >>> from bilby.gw.detector import PowerSpectralDensity
+            >>> psd = PowerSpectralDensity.from_aligo()
+            >>> acf = AutoCovarianceFunction.from_power_spectral_density_array(
+            ...     psd.frequency_array, psd.psd_array, minimum_frequency=20,
+            ...     maximum_frequency=896, sampling_frequency=2048)
+
+        From a precomputed ACF:
+
+        .. code-block:: python
+
+            >>> acf_array = acf.get_acf_array(4096)
+            >>> acf = AutoCovarianceFunction(acf_array=acf_array, sampling_frequency=2048)
+
+        Parameters
+        ==========
+        acf_array: array_like
+            The ACF at lags 0, 1 / sampling_frequency, 2 / sampling_frequency, ...
+        sampling_frequency: float
+            The sampling frequency (in Hz) of the lags; must equal that of the data
+        minimum_frequency: float, optional
+            The lower edge (in Hz) of the frequency band used to compute the ACF
+        maximum_frequency: float, optional
+            The upper edge (in Hz) of the frequency band used to compute the ACF
+        fill_value: float, optional
+            The PSD value used outside the band
+        fill_multiplier: float, optional
+            The multiple of the maximum in-band PSD used outside the band
+
+        """
+        self.acf_array = None if acf_array is None else np.asarray(acf_array, dtype=float)
+        self.sampling_frequency = None if sampling_frequency is None else float(sampling_frequency)
+        self.minimum_frequency = minimum_frequency
+        self.maximum_frequency = maximum_frequency
+        self.fill_value = fill_value
+        self.fill_multiplier = fill_multiplier
+
+    @property
+    def meta_data(self):
+        return dict(
+            sampling_frequency=self.sampling_frequency,
+            minimum_frequency=self.minimum_frequency,
+            maximum_frequency=self.maximum_frequency,
+            fill_value=self.fill_value,
+            fill_multiplier=self.fill_multiplier,
+        )
+
+    @classmethod
+    def from_power_spectral_density_array(
+            cls, frequency_array, psd_array, minimum_frequency, maximum_frequency,
+            sampling_frequency=None, fill_value=None, fill_multiplier=1e4):
+        """ Compute the ACF from a one-sided power spectral density (PSD).
+
+        The PSD outside [`minimum_frequency`, `maximum_frequency`] is set to
+        `fill_value`, or to `fill_multiplier` times its maximum in the band, and
+        the ACF is its inverse Fourier transform (Wiener-Khinchin theorem).
+
+        Parameters
+        ==========
+        frequency_array: array_like
+            Uniformly spaced frequencies (in Hz) up to the Nyquist frequency
+        psd_array: array_like
+            The one-sided PSD at `frequency_array`
+        minimum_frequency: float
+            The lower edge of the frequency band (in Hz)
+        maximum_frequency: float
+            The upper edge of the frequency band (in Hz)
+        sampling_frequency: float, optional
+            The sampling frequency of the data (in Hz); a warning is logged if it
+            is not twice the highest frequency in `frequency_array`
+        fill_value: float, optional
+            The PSD value used outside the band; overrides `fill_multiplier`
+        fill_multiplier: float
+            The multiple of the maximum in-band PSD used outside the band
+
+        Returns
+        =======
+        AutoCovarianceFunction
+
+        """
+        frequency_array = np.asarray(frequency_array, dtype=float)
+        psd_array = np.asarray(psd_array, dtype=float)
+        delta_f = frequency_array[1] - frequency_array[0]
+        frequency_values_missing = int(round(frequency_array[0] / delta_f))
+        if frequency_values_missing > 0:
+            frequency_array = np.concatenate(
+                [np.arange(frequency_values_missing) * delta_f, frequency_array])
+            psd_array = np.concatenate([np.zeros(frequency_values_missing), psd_array])
+
+        psd_sampling_frequency = 2 * frequency_array[-1]
+        if sampling_frequency is not None and not np.isclose(
+                sampling_frequency, psd_sampling_frequency, rtol=1e-9, atol=0):
+            logger.warning(
+                f"The sampling frequency given ({sampling_frequency:g} Hz) is not equal to the one "
+                f"computed from the PSD ({psd_sampling_frequency:g} Hz); using "
+                f"{psd_sampling_frequency:g} Hz. An ACF computed with {sampling_frequency:g} Hz "
+                "would be wrongly scaled and its lags wrongly spaced, giving wrong likelihood "
+                "and SNR values.")
+
+        band = (frequency_array >= minimum_frequency) & (frequency_array <= maximum_frequency)
+        if fill_value is None:
+            fill_value = fill_multiplier * np.max(psd_array[band])
+        else:
+            fill_multiplier = None
+        psd_array = psd_array.copy()
+        psd_array[~band] = fill_value
+        acf_array = 0.5 * np.fft.irfft(psd_array) * psd_sampling_frequency
+        return cls(acf_array=acf_array, sampling_frequency=psd_sampling_frequency,
+                   minimum_frequency=minimum_frequency, maximum_frequency=maximum_frequency,
+                   fill_value=float(fill_value), fill_multiplier=fill_multiplier)
+
+    def get_acf_array(self, number_of_samples):
+        """ The first `number_of_samples` lags of the ACF
+
+        Parameters
+        ==========
+        number_of_samples: int
+            The number of samples in the analysis window
+
+        Returns
+        =======
+        array_like: The first row of the noise covariance matrix
+
+        """
+        if number_of_samples > len(self.acf_array):
+            raise ValueError(
+                f"The analysis segment ({number_of_samples} samples) is longer than the ACF "
+                f"({len(self.acf_array)} samples)")
+        return self.acf_array[:number_of_samples]
+
+
+__all__ = ["AutoCovarianceFunction"]

bilby/gw/detector/condition_data.py

@@ -0,0 +1,84 @@
+import numpy as np
+
+
+def condition_strain_data(strain_data, sampling_frequency, start_time, duration,
+                          preserve_time=None, trim=0.25, taper_and_trim=True,
+                          remove_mean=True):
+    """ Downsample time-domain strain data around an analysis window
+
+    By default, the data in the window plus a buffer of
+    `trim` / (1 - 2 `trim`) x `duration` on each side are tapered with a Tukey
+    window (alpha = `trim`), low-passed below the new Nyquist frequency,
+    downsampled, and the buffer is removed. The samples are aligned so that the
+    sample nearest `preserve_time` is kept. With `taper_and_trim` = False,
+    :code:`scipy.signal.decimate` is used instead. The result replaces the
+    data in `strain_data`.
+
+    Parameters
+    ==========
+    strain_data: bilby.gw.detector.InterferometerStrainData
+        The strain data, modified in place
+    sampling_frequency: float
+        The new sampling frequency (in Hz); must divide the current one
+    start_time: float
+        The GPS start time of the analysis window
+    duration: float
+        The duration of the analysis window (in s)
+    preserve_time: float, optional
+        GPS time whose nearest sample is kept (default: `start_time`)
+    trim: float
+        The fraction of the data removed at each end
+    taper_and_trim: bool
+        If True, taper, low-pass and trim; if False, use :code:`scipy.signal.decimate`
+    remove_mean: bool
+        If True, subtract the mean of the result
+
+    """
+    from scipy.signal import decimate
+    from scipy.signal.windows import tukey
+
+    factor = int(round(strain_data.sampling_frequency / sampling_frequency))
+    strain = np.asarray(strain_data.time_domain_strain, dtype=float)
+    times = np.asarray(strain_data.time_array, dtype=float)
+    if taper_and_trim:
+        # step 1: the window and a buffer on each side
+        buffer = duration * trim / (1 - 2 * trim)
+        first, last = start_time - buffer, start_time + duration + buffer
+        data_end = strain_data.start_time + strain_data.duration
+        if (strain_data.start_time > first + 1 / strain_data.sampling_frequency
+                or data_end < last - 1 / strain_data.sampling_frequency):
+            raise ValueError(
+                f"The data ({strain_data.start_time:.6f} - {data_end:.6f}) must cover {first:.6f} - {last:.6f}: "
+                f"the analysis window plus {buffer:g} s on each side")
+        keep = slice(max(int(round((first - strain_data.start_time) * strain_data.sampling_frequency)), 0),
+                     int(round((last - strain_data.start_time) * strain_data.sampling_frequency)))
+        strain, times = strain[keep], times[keep]
+
+    # step 2: keep the sample nearest preserve_time
+    shift = int(np.argmin(abs(times - (start_time if preserve_time is None else preserve_time)))) % factor
+    if taper_and_trim and trim > 0:  # the wrapped samples are at the end and are removed in step 4
+        strain, times = np.roll(strain, -shift), np.roll(times, -shift)
+    else:
+        strain, times = strain[shift:], times[shift:]
+
+    if taper_and_trim:
+        # step 3: taper, low-pass and downsample
+        strain_fd = np.fft.rfft(strain * tukey(len(strain), trim))
+        strain_fd[np.fft.rfftfreq(len(strain), 1 / strain_data.sampling_frequency) > sampling_frequency / 2] = 0
+        strain, times = np.fft.irfft(strain_fd, n=len(strain))[::factor], times[::factor]
+        # step 4: remove the fraction trim at each end, but always keep the analysis samples
+        window = int(np.argmin(abs(times - start_time)))  # first sample of the analysis window
+        keep = slice(min(int(round(trim * len(strain))), window),
+                     max(int(round((1 - trim) * len(strain))), window + int(round(duration * sampling_frequency))))
+        strain, times = strain[keep], times[keep]
+    else:
+        strain, times = decimate(strain, factor, zero_phase=True), times[::factor]
+
+    if remove_mean:
+        strain = strain - np.mean(strain)
+    strain_data.set_from_time_domain_strain(
+        strain, sampling_frequency=sampling_frequency, duration=len(strain) / sampling_frequency,
+        start_time=times[0])
+
+
+__all__ = ["condition_strain_data"]