diff --git a/bilby/gw/__init__.py b/bilby/gw/__init__.py index cd09bc6f6..a996f853c 100644 --- a/bilby/gw/__init__.py +++ b/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 diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py index b55d49379..c2d3950d0 100644 --- a/bilby/gw/conversion.py +++ b/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( diff --git a/bilby/gw/detector/__init__.py b/bilby/gw/detector/__init__.py index 51ff98279..96423347c 100644 --- a/bilby/gw/detector/__init__.py +++ b/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 * diff --git a/bilby/gw/detector/acf.py b/bilby/gw/detector/acf.py new file mode 100644 index 000000000..fad7b5c4d --- /dev/null +++ b/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"] diff --git a/bilby/gw/detector/condition_data.py b/bilby/gw/detector/condition_data.py new file mode 100644 index 000000000..c7718bc5b --- /dev/null +++ b/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"] diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index f8eda202a..b41e6aabe 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -14,6 +14,8 @@ from .calibration import Recalibrate from .geometry import InterferometerGeometry from .strain_data import InterferometerStrainData +from ..time_domain_utils import ( + gohberg_semencul_product, gohberg_semencul_vectors, place_time_domain_signal, times_from_merger) from ..conversion import generate_all_bbh_parameters @@ -93,6 +95,11 @@ def __init__(self, name, power_spectral_density, minimum_frequency, maximum_freq maximum_frequency=maximum_frequency) self.meta_data = dict(name=name) self.reference_time = None + self.analysis_start_time = None + self.analysis_duration = None + self.analysis_cut_time = None + self._autocovariance_function = None + self.gohberg_semencul_vectors = None def __eq__(self, other): if self.name == other.name and \ @@ -263,6 +270,145 @@ def set_strain_data_from_zero_noise( sampling_frequency=sampling_frequency, duration=duration, start_time=start_time) + def set_analysis_window(self, start_time, duration, cut_time=None): + """ Set the analysis window used by time-domain likelihoods + + Setting it removes the autocovariance function, which must be set after it. + + Parameters + ========== + start_time: float + The GPS start time of the window + duration: float + The duration of the window (in s) + cut_time: float, optional + The GPS time of the inspiral/post-inspiral cut, if any + """ + self.analysis_start_time = float(start_time) + self.analysis_duration = float(duration) + self.analysis_cut_time = cut_time + self.meta_data["analysis_window"] = dict( + start_time=self.analysis_start_time, duration=self.analysis_duration, cut_time=cut_time) + self.autocovariance_function = None + + def condition_strain_data(self, sampling_frequency, trim=0.25, taper_and_trim=True, + remove_mean=True): + """ Downsample the time-domain strain data around the analysis window + + See :code:`bilby.gw.detector.condition_data.condition_strain_data`; the sample + nearest the cut time (or the window start) is kept. + + Parameters + ========== + sampling_frequency: float + The new sampling frequency (in Hz); must divide the current one + 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 .condition_data import condition_strain_data + if self.analysis_duration is None: + raise ValueError(f"{self.name}: set the analysis window first (set_analysis_window)") + condition_strain_data(self.strain_data, sampling_frequency, self.analysis_start_time, + self.analysis_duration, preserve_time=self.analysis_cut_time, trim=trim, + taper_and_trim=taper_and_trim, remove_mean=remove_mean) + + @property + def analysis_slice(self): + """ The slice of the time-domain strain data in the analysis window """ + if self.analysis_duration is None: + raise ValueError(f"{self.name}: set the analysis window first (set_analysis_window)") + first = int(np.argmin(abs(self.time_array - self.analysis_start_time))) + last = first + int(round(self.analysis_duration * self.sampling_frequency)) + if (last > len(self.time_array) + or abs(self.time_array[first] - self.analysis_start_time) > 1 / self.sampling_frequency): + raise ValueError( + f"{self.name}: the analysis window ({self.analysis_start_time:.6f} + {self.analysis_duration:g} s) " + f"is not inside the data ({self.start_time:.6f} + {self.duration:g} s)") + return slice(first, last) + + @property + def analysis_data(self): + """ The time-domain strain data in the analysis window """ + return self.time_domain_strain[self.analysis_slice] + + @property + def autocovariance_function(self): + """ The noise autocovariance function (:code:`bilby.gw.detector.AutoCovarianceFunction`) + + Setting it computes :code:`gohberg_semencul_vectors`, which give products with + the inverse noise covariance matrix of the analysis window. + """ + return self._autocovariance_function + + @autocovariance_function.setter + def autocovariance_function(self, autocovariance_function): + self._autocovariance_function = autocovariance_function + self.gohberg_semencul_vectors = None + self.meta_data.pop("autocovariance_function", None) + if autocovariance_function is None: + return + if self.analysis_duration is None: + raise ValueError(f"{self.name}: set the analysis window first (set_analysis_window)") + if self.strain_data.sampling_frequency is not None and not np.isclose( + autocovariance_function.sampling_frequency, self.strain_data.sampling_frequency, + rtol=1e-9, atol=0): + raise ValueError( + f"{self.name}: the sampling frequency of the autocovariance function " + f"({autocovariance_function.sampling_frequency:g} Hz) differs from that of the strain data " + f"({self.strain_data.sampling_frequency:g} Hz)") + self.meta_data["autocovariance_function"] = autocovariance_function.meta_data + self.gohberg_semencul_vectors = gohberg_semencul_vectors(autocovariance_function.get_acf_array( + int(round(self.analysis_duration * autocovariance_function.sampling_frequency)))) + + def set_autocovariance_function_from_power_spectral_density( + self, minimum_frequency=None, maximum_frequency=None, sampling_frequency=None, + duration_factor=16, fill_value=None, fill_multiplier=1e4): + """ Set the autocovariance function from the power spectral density + + See :code:`bilby.gw.detector.AutoCovarianceFunction.from_power_spectral_density_array`. + The power spectral density is evaluated on the frequency grid of the next + power of two at or above `duration_factor` x the number of samples in the + analysis window, so that the ACF has the sampling frequency of the data (as + in ringdown). + + Parameters + ========== + minimum_frequency: float, optional + The lower edge of the frequency band (in Hz), default: :code:`self.minimum_frequency` + maximum_frequency: float, optional + The upper edge of the frequency band (in Hz), default: :code:`self.maximum_frequency` + sampling_frequency: float, optional + The sampling frequency (in Hz), default: that of the strain data + duration_factor: float + The minimum duration of the autocovariance function divided by that of the analysis window + 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 + """ + from .acf import AutoCovarianceFunction + if self.analysis_duration is None: + raise ValueError(f"{self.name}: set the analysis window first (set_analysis_window)") + if sampling_frequency is None: + sampling_frequency = self.strain_data.sampling_frequency or 2048 + if minimum_frequency is None: + minimum_frequency = self.minimum_frequency + if maximum_frequency is None: + maximum_frequency = self.maximum_frequency + number_of_samples = int(round(self.analysis_duration * sampling_frequency)) + fft_length = 2 ** int(np.ceil(np.log2(duration_factor * number_of_samples))) + frequency_array = np.fft.rfftfreq(fft_length, 1 / sampling_frequency) + with np.errstate(invalid="ignore", divide="ignore"): + psd_array = self.power_spectral_density.get_power_spectral_density_array(frequency_array) + self.autocovariance_function = AutoCovarianceFunction.from_power_spectral_density_array( + frequency_array, psd_array, minimum_frequency=minimum_frequency, + maximum_frequency=maximum_frequency, sampling_frequency=sampling_frequency, + fill_value=fill_value, fill_multiplier=fill_multiplier) + def antenna_response(self, ra, dec, time, psi, mode): """ Calculate the antenna response function for a given sky location @@ -366,6 +512,199 @@ def get_detector_response(self, waveform_polarizations, parameters, frequencies= return signal_ifo + def get_time_domain_detector_response(self, waveform_polarizations, parameters, + placement="nearest", times=None): + """ Get the time-domain detector response for a particular waveform + + Parameters + ========== + waveform_polarizations: dict + Output of a time-domain source model, e.g. + :code:`bilby.gw.source.lal_binary_black_hole_time_domain` + parameters: dict + Parameters describing position and time of arrival of the signal + placement: str + How the waveform is placed on the data samples: "nearest" or "fd_shift" + (see :code:`bilby.gw.time_domain_utils.place_time_domain_signal`), or + "model_time" for a model-time source model, e.g. + :code:`bilby.gw.source.lal_binary_black_hole_model_time` + times: array_like, optional + The data sample times, default: those of the analysis window + + Returns + ======= + array_like: The detector response at `times` (None if a model-time waveform failed) + + Notes + ===== + If the :code:`reference_time` attribute is not :code:`None`, this is + used to set the time at which the antenna response is evaluated, + otherwise the provided :code:`parameters["geocent_time"]` is used. + """ + if times is None: + times = self.time_array[self.analysis_slice] + model_time = "polarizations_at" in waveform_polarizations + if model_time != (placement == "model_time"): + raise ValueError( + 'placement="model_time" needs a model-time source model (e.g. ' + 'bilby.gw.source.lal_binary_black_hole_model_time), and a model-time source ' + f'model needs placement="model_time"; got placement={placement!r}') + + if self.reference_time is None: + antenna_time = parameters["geocent_time"] + else: + antenna_time = self.reference_time + + time_shift = self.time_delay_from_geocenter( + parameters['ra'], parameters['dec'], parameters['geocent_time']) + + if model_time: + polarizations = waveform_polarizations["polarizations_at"](times_from_merger( + times[0], len(times), self.sampling_frequency, parameters['geocent_time'], time_shift)) + if polarizations is None: + return None + else: + polarizations = {mode: value for mode, value in waveform_polarizations.items() + if mode != "epoch"} + + signal_ifo = 0 + for mode, polarization in polarizations.items(): + det_response = self.antenna_response( + parameters['ra'], parameters['dec'], antenna_time, parameters['psi'], mode) + signal_ifo = signal_ifo + polarization * det_response + + if model_time: + return signal_ifo + + # Be careful to first subtract the two GPS times which are ~1e9 sec. + # And then add the time_shift which varies at ~1e-5 sec + dt_geocent = parameters['geocent_time'] - times[0] + return place_time_domain_signal( + signal_ifo, merger_index=-waveform_polarizations["epoch"] * self.sampling_frequency, + arrival_index=(dt_geocent + time_shift) * self.sampling_frequency, + number_of_samples=len(times), placement=placement) + + def inject_signal_time_domain(self, parameters, injection_polarizations=None, + waveform_generator=None, placement="nearest"): + """ General time-domain signal injection method. + Provide the injection parameters and either the injection polarizations + or the waveform generator to inject a signal into the time-domain strain + data. Defaults to the injection polarizations if both are given. + + Parameters + ========== + parameters: dict + Parameters of the injection. + injection_polarizations: dict, optional + Polarizations of waveform to inject, output of + `waveform_generator.time_domain_strain()`. If + `waveform_generator` is also given, the injection_polarizations will + be calculated directly and this argument can be ignored. + waveform_generator: bilby.gw.waveform_generator.WaveformGenerator, optional + A WaveformGenerator instance using the time-domain source model to + inject. If `injection_polarizations` is given, this will be ignored. + placement: str + "nearest", "fd_shift" or "model_time"; see + :code:`get_time_domain_detector_response` + + Returns + ======= + injection_polarizations: dict + The injected polarizations. This is the same as the injection_polarizations parameters + if it was passed in. Otherwise it is the return value of waveform_generator.time_domain_strain(). + """ + if injection_polarizations is None and waveform_generator is None: + raise ValueError( + "inject_signal_time_domain needs one of waveform_generator or " + "injection_polarizations.") + elif injection_polarizations is not None: + self.inject_signal_time_domain_from_waveform_polarizations( + parameters=parameters, injection_polarizations=injection_polarizations, placement=placement) + elif waveform_generator is not None: + injection_polarizations = self.inject_signal_time_domain_from_waveform_generator( + parameters=parameters, waveform_generator=waveform_generator, placement=placement) + return injection_polarizations + + def inject_signal_time_domain_from_waveform_generator(self, parameters, waveform_generator, + placement="nearest"): + """ Inject a time-domain signal using a waveform generator and a set of parameters. + Alternative to `inject_signal_time_domain` and + `inject_signal_time_domain_from_waveform_polarizations` + + Parameters + ========== + parameters: dict + Parameters of the injection. + waveform_generator: bilby.gw.waveform_generator.WaveformGenerator + A WaveformGenerator instance using the time-domain source model to inject. + placement: str + "nearest", "fd_shift" or "model_time" + + Returns + ======= + injection_polarizations: dict + The internally generated injection parameters + """ + if waveform_generator.sampling_frequency is not None and not np.isclose( + waveform_generator.sampling_frequency, self.sampling_frequency, rtol=1e-9, atol=0): + raise ValueError( + f"{self.name}: the sampling frequency of the waveform generator " + f"({waveform_generator.sampling_frequency:g} Hz) differs from that of the strain data " + f"({self.sampling_frequency:g} Hz)") + injection_polarizations = \ + waveform_generator.time_domain_strain(parameters) + self.inject_signal_time_domain_from_waveform_polarizations( + parameters=parameters, injection_polarizations=injection_polarizations, placement=placement) + return injection_polarizations + + def inject_signal_time_domain_from_waveform_polarizations(self, parameters, injection_polarizations, + placement="nearest"): + """ Inject a time-domain signal into the detector from a dict of waveform polarizations. + Alternative to `inject_signal_time_domain` and + `inject_signal_time_domain_from_waveform_generator`. + + Parameters + ========== + parameters: dict + Parameters of the injection. + injection_polarizations: dict + Polarizations of waveform to inject, output of + `waveform_generator.time_domain_strain()`. + placement: str + "nearest", "fd_shift" or "model_time" + """ + if "epoch" not in injection_polarizations and "polarizations_at" not in injection_polarizations: + raise ValueError("injection_polarizations must come from a time-domain source model " + "(with an 'epoch' or 'polarizations_at' entry)") + if self.gohberg_semencul_vectors is None: + raise ValueError(f"{self.name}: set the autocovariance function before injecting") + if not self.strain_data.time_within_data(parameters['geocent_time']): + logger.warning( + 'Injecting signal outside segment, start_time={}, merger time={}.' + .format(self.strain_data.start_time, parameters['geocent_time'])) + + signal_ifo = self.get_time_domain_detector_response( + injection_polarizations, parameters, placement=placement, times=self.time_array) + if signal_ifo is None: + raise ValueError(f"{self.name}: the injection waveform failed") + self.strain_data._time_domain_strain = self.time_domain_strain + signal_ifo + self.strain_data._frequency_domain_strain = None + + signal = signal_ifo[self.analysis_slice] + inverse_covariance_signal = gohberg_semencul_product(*self.gohberg_semencul_vectors, signal) + self.meta_data['optimal_SNR'] = float(np.dot(signal, inverse_covariance_signal)) ** 0.5 + self.meta_data['matched_filter_SNR'] = float( + np.dot(self.analysis_data, inverse_covariance_signal)) / self.meta_data['optimal_SNR'] + self.meta_data['parameters'] = parameters + self.meta_data['injection_domain'] = "time" + self.meta_data['injection_placement'] = placement + + logger.info("Injected time-domain signal in {}:".format(self.name)) + logger.info(" optimal SNR = {:.2f}".format(self.meta_data['optimal_SNR'])) + logger.info(" matched filter SNR = {:.2f}".format(self.meta_data['matched_filter_SNR'])) + for key in parameters: + logger.info(' {} = {}'.format(key, parameters[key])) + def check_signal_duration(self, parameters, raise_error=True): """ Check that the signal with the given parameters fits in the data diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py index dddc2f364..0f4034362 100644 --- a/bilby/gw/detector/networks.py +++ b/bilby/gw/detector/networks.py @@ -52,12 +52,18 @@ def reference_time(self, time): def _check_interferometers(self): """Verify IFOs 'duration', 'start_time', 'sampling_frequency' are the same. + If every interferometer has a time-domain analysis window, 'start_time' + is not checked. + If the above attributes are not the same, then the attributes are checked to see if they are the same up to 5 decimal places. If both checks fail, then a ValueError is raised. """ consistent_attributes = ["duration", "start_time", "sampling_frequency"] + if all(interferometer.analysis_duration is not None for interferometer in self): + # time-domain analysis windows can start at different times in each detector + consistent_attributes = ["duration", "sampling_frequency"] for attribute in consistent_attributes: x = [ getattr(interferometer.strain_data, attribute) @@ -195,6 +201,113 @@ def inject_signal( return all_injection_polarizations + def set_analysis_windows(self, duration, start_time=None, analysis_segment="imr", + segment_cut_time=None, cut_reference_parameters=None): + """ Set the analysis window of each interferometer + + For "imr", every interferometer gets the window [`start_time`, + `start_time` + `duration`). For "inspiral" ("post_inspiral"), the window of + each interferometer ends (starts) at its cut time: `segment_cut_time` after + the arrival of the merger, at the `geocent_time`, `ra` and `dec` in + `cut_reference_parameters`. + + Parameters + ========== + duration: float + The duration of the windows (in s) + start_time: float, optional + The GPS start time of the windows; required for "imr" + analysis_segment: str + "imr", "inspiral" or "post_inspiral" + segment_cut_time: float, optional + The cut time relative to the merger (in s); required for "inspiral" and "post_inspiral" + cut_reference_parameters: dict, optional + Dictionary with `geocent_time`, `ra` and `dec`; required for "inspiral" and "post_inspiral" + """ + if analysis_segment not in ("imr", "inspiral", "post_inspiral"): + raise ValueError( + "analysis_segment must be 'imr', 'inspiral' or 'post_inspiral', " + f"not {analysis_segment!r}") + for interferometer in self: + window_start = start_time + cut_time = None + if analysis_segment != "imr": + cut_time = cut_reference_parameters["geocent_time"] + segment_cut_time + \ + interferometer.time_delay_from_geocenter( + cut_reference_parameters["ra"], cut_reference_parameters["dec"], + cut_reference_parameters["geocent_time"]) + window_start = cut_time - duration if analysis_segment == "inspiral" else cut_time + interferometer.set_analysis_window(window_start, duration, cut_time) + interferometer.meta_data["analysis_window"].update( + analysis_segment=analysis_segment, segment_cut_time=segment_cut_time) + + def set_autocovariance_functions_from_power_spectral_densities(self, **kwargs): + """ Set the autocovariance function of each interferometer from its power spectral density + + See :code:`bilby.gw.detector.Interferometer.set_autocovariance_function_from_power_spectral_density` + for the keyword arguments. + """ + for interferometer in self: + interferometer.set_autocovariance_function_from_power_spectral_density(**kwargs) + + def inject_signal_time_domain( + self, + parameters=None, + injection_polarizations=None, + waveform_generator=None, + placement="nearest", + ): + """ Inject a time-domain signal into the time-domain strain data of each interferometer. + + Parameters + ========== + parameters: dict + Parameters of the injection. + injection_polarizations: dict + Polarizations of waveform to inject, output of + `waveform_generator.time_domain_strain()`. If + `waveform_generator` is also given, the injection_polarizations will + be calculated directly and this argument can be ignored. + waveform_generator: bilby.gw.waveform_generator.WaveformGenerator + A WaveformGenerator instance using the time-domain source model to + inject. If `injection_polarizations` is given, this will be ignored. + placement: str + "nearest", "fd_shift" or "model_time"; see + :code:`bilby.gw.detector.Interferometer.get_time_domain_detector_response` + + Returns + ======= + injection_polarizations: dict + """ + if injection_polarizations is None: + if waveform_generator is not None: + if waveform_generator.sampling_frequency is not None and not np.isclose( + waveform_generator.sampling_frequency, self.sampling_frequency, rtol=1e-9, atol=0): + raise ValueError( + f"The sampling frequency of the waveform generator " + f"({waveform_generator.sampling_frequency:g} Hz) differs from that of the strain " + f"data ({self.sampling_frequency:g} Hz)") + injection_polarizations = waveform_generator.time_domain_strain( + parameters + ) + else: + raise ValueError( + "inject_signal_time_domain needs one of waveform_generator or " + "injection_polarizations." + ) + + all_injection_polarizations = list() + for interferometer in self: + all_injection_polarizations.append( + interferometer.inject_signal_time_domain( + parameters=parameters, + injection_polarizations=injection_polarizations, + placement=placement, + ) + ) + + return all_injection_polarizations + def save_data(self, outdir, label=None): """Creates a save file for the data in plain text format diff --git a/bilby/gw/likelihood/__init__.py b/bilby/gw/likelihood/__init__.py index d752f77a9..360dff513 100644 --- a/bilby/gw/likelihood/__init__.py +++ b/bilby/gw/likelihood/__init__.py @@ -3,6 +3,8 @@ from .roq import BilbyROQParamsRangeError, ROQGravitationalWaveTransient from .multiband import MBGravitationalWaveTransient from .relative import RelativeBinningGravitationalWaveTransient +from .time_domain import TimeDomainGravitationalWaveTransient +from .time_domain_heterodyning import TimeDomainHeterodyningGravitationalWaveTransient from ..source import lal_binary_black_hole from ..waveform_generator import WaveformGenerator diff --git a/bilby/gw/likelihood/time_domain.py b/bilby/gw/likelihood/time_domain.py new file mode 100644 index 000000000..436b7a018 --- /dev/null +++ b/bilby/gw/likelihood/time_domain.py @@ -0,0 +1,307 @@ +import copy + +import numpy as np + +from ...core.utils import logger +from ..time_domain_utils import gohberg_semencul_product +from .base import GravitationalWaveTransient + + +class TimeDomainGravitationalWaveTransient(GravitationalWaveTransient): + """ A gravitational-wave transient likelihood object evaluated in the time domain + + It computes the log-likelihood of the data in the analysis window of each + interferometer, assuming stationary colored Gaussian noise described by an + autocovariance function, as in time-domain ringdown analyses (Carullo et al. + arXiv:1902.07527, Isi & Farr arXiv:2107.05609). Products with the inverse + noise covariance matrix use the Gohberg-Semencul formula. See Sharma et al. + arXiv:2601.11239. + + Parameters + ========== + interferometers: list, bilby.gw.detector.InterferometerList + A list of `bilby.detector.Interferometer` instances - contains the + detector data, analysis windows and autocovariance functions + waveform_generator: `bilby.waveform_generator.WaveformGenerator` + An object which computes the time-domain strain of the signal, + given some set of parameters + distance_marginalization: bool, optional + If true, marginalize over distance in the likelihood. + This uses a look up table calculated at run time. + The distance prior is set to be a delta function at the minimum + distance allowed in the prior being marginalised over. + priors: dict, optional + If given, used in the distance marginalization. + distance_marginalization_lookup_table: (dict, str), optional + If a dict, dictionary containing the lookup_table, distance_array, + (distance) prior_array, and reference_distance used to construct + the table. + If a string the name of a file containing these quantities. + The lookup table is stored after construction in either the + provided string or a default location: + '.distance_marginalization_lookup_dmin{}_dmax{}_n{}.npz' + reference_frame: (str, bilby.gw.detector.InterferometerList, list), optional + Definition of the reference frame for the sky location. + + - :code:`sky`: sample in RA/dec, this is the default + - e.g., :code:`"H1L1", ["H1", "L1"], InterferometerList(["H1", "L1"])`: + sample in azimuth and zenith, `azimuth` and `zenith` defined in the + frame where the z-axis is aligned the the vector connecting H1 + and L1. + + time_reference: str, optional + Name of the reference for the sampled time parameter. + + - :code:`geocent`/:code:`geocenter`: sample in the time at the + Earth's center, this is the default + - e.g., :code:`H1`: sample in the time of arrival at H1 + + placement: str, optional + How the waveform is placed on the data samples: "nearest" (default), + "fd_shift" or "model_time"; see + :code:`bilby.gw.detector.Interferometer.get_time_domain_detector_response` + time_marginalization, phase_marginalization, calibration_marginalization: bool, optional + Not implemented for this likelihood; must be False + + Returns + ======= + Likelihood: `bilby.core.likelihood.Likelihood` + A likelihood object, able to compute the likelihood of the data given + some model parameters + """ + + likelihood_domain = "time" + + def __init__( + self, interferometers, waveform_generator, distance_marginalization=False, + priors=None, distance_marginalization_lookup_table=None, reference_frame="sky", + time_reference="geocenter", placement="nearest", time_marginalization=False, + phase_marginalization=False, calibration_marginalization=False, + ): + for name, value in dict( + time_marginalization=time_marginalization, + phase_marginalization=phase_marginalization, + calibration_marginalization=calibration_marginalization).items(): + if value: + raise NotImplementedError( + f"{name} is not implemented for {self.__class__.__name__}") + if getattr(waveform_generator, "time_domain_source_model", None) is None: + raise ValueError( + f"{self.__class__.__name__} needs a waveform generator with a " + "time_domain_source_model, e.g. bilby.gw.source.lal_binary_black_hole_time_domain") + if placement not in ("nearest", "fd_shift", "model_time"): + raise ValueError( + f"placement must be 'nearest', 'fd_shift' or 'model_time', not {placement!r}") + + super(TimeDomainGravitationalWaveTransient, self).__init__( + interferometers=interferometers, + waveform_generator=waveform_generator, + time_marginalization=False, + distance_marginalization=distance_marginalization, + phase_marginalization=False, + calibration_marginalization=False, + priors=priors, + distance_marginalization_lookup_table=distance_marginalization_lookup_table, + jitter_time=False, + reference_frame=reference_frame, + time_reference=time_reference, + ) + self.placement = placement + self._setup_noise_model() + + def __repr__(self): + return self.__class__.__name__ + ( + '(interferometers={},\n\twaveform_generator={},\n\tdistance_marginalization={}, ' + 'priors={}, placement={})'.format( + self.interferometers, self.waveform_generator, self.distance_marginalization, + self.priors, self.placement)) + + def _setup_noise_model(self): + """ Compute and store C^-1 d and d^T C^-1 d for each interferometer """ + self._data_times_inverse_covariance = dict() + self._data_inverse_covariance_data = dict() + for interferometer in self.interferometers: + acf = interferometer.autocovariance_function + if acf is None or acf.acf_array is None: + raise ValueError(f"{interferometer.name}: no autocovariance_function set") + if not np.isclose(acf.sampling_frequency, interferometer.sampling_frequency, + rtol=1e-9, atol=0): + raise ValueError( + f"{interferometer.name}: autocovariance function sampling frequency " + f"{acf.sampling_frequency:g} Hz differs from the data's " + f"{interferometer.sampling_frequency:g} Hz") + nyquist_frequency = interferometer.sampling_frequency / 2 + if interferometer.maximum_frequency < nyquist_frequency and not np.isclose( + interferometer.maximum_frequency, nyquist_frequency, rtol=1e-9, atol=0): + logger.warning( + f"{interferometer.name}: {interferometer.maximum_frequency:g} - {nyquist_frequency:g} Hz " + "is outside the band and filled in the autocovariance function (consider downsampling " + f"the data to {2 * interferometer.maximum_frequency:g} Hz with condition_strain_data)") + data = interferometer.analysis_data + data_times_inverse_covariance = gohberg_semencul_product( + *interferometer.gohberg_semencul_vectors, data) + self._data_times_inverse_covariance[interferometer.name] = data_times_inverse_covariance + self._data_inverse_covariance_data[interferometer.name] = float( + np.dot(data_times_inverse_covariance, data)) + + def _compute_full_waveform(self, signal_polarizations, interferometer, parameters): + """ Project the waveform polarizations onto the interferometer and place them + on the samples of its analysis window. + + Parameters + ========== + signal_polarizations: dict + Output of the time-domain source model + interferometer: bilby.gw.detector.Interferometer + Interferometer to compute the response with respect to. + parameters: dict + The parameters of the signal + """ + return interferometer.get_time_domain_detector_response( + signal_polarizations, parameters, placement=self.placement) + + def calculate_snrs(self, waveform_polarizations, interferometer, *, return_array=True, parameters): + """ Compute the snrs + + Parameters + ========== + waveform_polarizations: dict + A dictionary of waveform polarizations and the corresponding array + interferometer: bilby.gw.detector.Interferometer + The bilby interferometer object + return_array: bool + Not used; kept for the :code:`GravitationalWaveTransient` interface + parameters: dict + The parameters of the signal + + Returns + ======= + calculated_snrs: _CalculatedSNRs + An object containing the SNR quantities (None if a model-time waveform failed) + """ + signal = self._compute_full_waveform( + signal_polarizations=waveform_polarizations, + interferometer=interferometer, + parameters=parameters, + ) + if signal is None: + return None + d_inner_h = float(np.dot(self._data_times_inverse_covariance[interferometer.name], signal)) + optimal_snr_squared = float(np.dot( + signal, gohberg_semencul_product(*interferometer.gohberg_semencul_vectors, signal))) + if optimal_snr_squared > 0: + matched_filter_snr = d_inner_h / optimal_snr_squared ** 0.5 + else: + matched_filter_snr = 0.0 + return self._CalculatedSNRs( + d_inner_h=d_inner_h, + optimal_snr_squared=optimal_snr_squared, + complex_matched_filter_snr=matched_filter_snr, + d_inner_h_array=None, + optimal_snr_squared_array=None, + ) + + def _calculate_noise_log_likelihood(self): + return -0.5 * sum(self._data_inverse_covariance_data.values()) + + def log_likelihood_ratio(self, parameters): + parameters = copy.deepcopy(parameters) + parameters.update(self.get_sky_frame_parameters(parameters)) + waveform_polarizations = self.waveform_generator.time_domain_strain(parameters) + if waveform_polarizations is None: + return np.nan_to_num(-np.inf) + + total_snrs = self._CalculatedSNRs() + for interferometer in self.interferometers: + snrs = self.calculate_snrs( + waveform_polarizations=waveform_polarizations, + interferometer=interferometer, + parameters=parameters, + ) + if snrs is None: + return np.nan_to_num(-np.inf) + total_snrs += snrs + log_l = self.compute_log_likelihood_from_snrs(total_snrs, parameters=parameters) + return float(np.real(log_l)) + + def compute_per_detector_log_likelihood(self, parameters): + parameters.update(self.get_sky_frame_parameters(parameters)) + waveform_polarizations = self.waveform_generator.time_domain_strain(parameters) + for interferometer in self.interferometers: + per_detector_snr = self.calculate_snrs( + waveform_polarizations=waveform_polarizations, + interferometer=interferometer, + parameters=parameters, + ) + parameters['{}_log_likelihood'.format(interferometer.name)] = ( + np.nan_to_num(-np.inf) if per_detector_snr is None else + self.compute_log_likelihood_from_snrs(per_detector_snr, parameters=parameters)) + return parameters.copy() + + def generate_posterior_sample_from_marginalized_likelihood(self, parameters): + """ Reconstruct the distance posterior from a run which used a likelihood + which explicitly marginalised over distance. + + See Eq. (C29-C32) of https://arxiv.org/abs/1809.02293 + + Returns + ======= + sample: dict + Returns the parameters with new samples. + """ + if not self.distance_marginalization: + return parameters + parameters['luminosity_distance'] = self.generate_distance_sample_from_marginalized_likelihood( + parameters=parameters) + return parameters.copy() + + def generate_distance_sample_from_marginalized_likelihood( + self, signal_polarizations=None, *, parameters): + """ Generate a single sample from the posterior distribution for luminosity + distance when using a likelihood which explicitly marginalises over + distance. + + See Eq. (C29-C32) of https://arxiv.org/abs/1809.02293 + + Parameters + ========== + signal_polarizations: dict, optional + Polarizations modes of the template. + Note: These are rescaled in place after the distance sample is + generated to allow further parameter reconstruction to occur. + + Returns + ======= + new_distance: float + Sample from the distance posterior. + """ + if signal_polarizations is None: + parameters.update(self.get_sky_frame_parameters(parameters)) + signal_polarizations = copy.deepcopy( + self.waveform_generator.time_domain_strain(parameters)) + return super(TimeDomainGravitationalWaveTransient, self).\ + generate_distance_sample_from_marginalized_likelihood( + signal_polarizations=signal_polarizations, parameters=parameters) + + def _rescale_signal(self, signal, new_distance): + factor = self._ref_dist / new_distance + if "polarizations_at" in signal: + polarizations_at = signal["polarizations_at"] + + def rescaled(times): + polarizations = polarizations_at(times) + if polarizations is None: + return None + return {mode: value * factor for mode, value in polarizations.items()} + + signal["polarizations_at"] = rescaled + return + for mode in signal: + if mode != "epoch": + signal[mode] *= factor + + @property + def meta_data(self): + meta_data = super(TimeDomainGravitationalWaveTransient, self).meta_data + meta_data.update(likelihood_domain=self.likelihood_domain, placement=self.placement) + return meta_data diff --git a/bilby/gw/likelihood/time_domain_heterodyning.py b/bilby/gw/likelihood/time_domain_heterodyning.py new file mode 100644 index 000000000..d4a49c0ff --- /dev/null +++ b/bilby/gw/likelihood/time_domain_heterodyning.py @@ -0,0 +1,658 @@ +import copy + +import numpy as np +from scipy.optimize import differential_evolution + +from ...core.utils import gravitational_constant, logger, solar_mass, speed_of_light +from ..conversion import component_masses_to_chirp_mass +from ..time_domain_utils import gohberg_semencul_product, times_from_merger +from .relative import RelativeBinningGravitationalWaveTransient +from .time_domain import TimeDomainGravitationalWaveTransient + + +class TimeDomainHeterodyningGravitationalWaveTransient(TimeDomainGravitationalWaveTransient): + """ A time-domain gravitational-wave transient likelihood object which uses + heterodyning (relative binning) to calculate a fast likelihood, for the + IMRPhenomT, IMRPhenomTHM, IMRPhenomTP and IMRPhenomTPHM waveform models. + + The waveform is evaluated only at the edges of time bins, and the + likelihood is computed from summary data built from a fiducial waveform. + See Sharma et al. arXiv:2601.11239, and Cornish arXiv:1007.4820 and + Zackay et al. arXiv:1806.08792 for the frequency-domain method. + + Parameters + ========== + interferometers: list, bilby.gw.detector.InterferometerList + A list of `bilby.detector.Interferometer` instances - contains the + detector data, analysis windows ("imr", "inspiral" or "post_inspiral") + and autocovariance functions + waveform_generator: `bilby.waveform_generator.WaveformGenerator` + A WaveformGenerator with + :code:`bilby.gw.source.lal_binary_black_hole_time_domain_heterodyning` + as its time-domain source model + fiducial_parameters: dict, optional + A starting guess for initial parameters of the event for finding the + maximum likelihood (fiducial) waveform. These should be specified in + the same parameter basis as the one that sampling is carried out in. + If no fiducial parameters are passed, a random point will be drawn from + the prior. + maximization_kwargs: dict, optional + Keyword arguments for :code:`scipy.optimize.differential_evolution` + update_fiducial_parameters: bool + Whether to try to update the fiducial parameters by finding the maximum + likelihood point using scipy differential evolution. Default is + :code:`False`. + parameters_to_be_updated: list, optional + The parameters over which the likelihood is maximized, within their + prior bounds. Default is :code:`["chirp_mass", "mass_ratio"]`. + distance_marginalization: bool, optional + If true, marginalize over distance in the likelihood. + This uses a look up table calculated at run time. + The distance prior is set to be a delta function at the minimum + distance allowed in the prior being marginalised over. + priors: dict, optional + If given, used in the distance marginalization and to find the + fiducial parameters. + distance_marginalization_lookup_table: (dict, str), optional + If a dict, dictionary containing the lookup_table, distance_array, + (distance) prior_array, and reference_distance used to construct + the table. + If a string the name of a file containing these quantities. + reference_frame: (str, bilby.gw.detector.InterferometerList, list), optional + Definition of the reference frame for the sky location, as in + :code:`GravitationalWaveTransient`. + time_reference: str, optional + Name of the reference for the sampled time parameter, as in + :code:`GravitationalWaveTransient`. + epsilon: tuple, optional + The bin width in the phase proxy before and after `time_split`; + smaller values give more bins (see :code:`time_domain_bin_edges`) + time_split: float, optional + The time relative to the merger (in s) at which `epsilon` changes + spacing: tuple, optional + The number of samples per bin after the merger: before + `spacing_times[0]`, before `spacing_times[1]` and after + spacing_times: tuple, optional + The times relative to the merger (in s) at which `spacing` changes + chi: float, optional + Tunable parameter of the phase proxy. See https://arxiv.org/abs/1806.08792. + fiducial_minimum_frequency: float, optional + The minimum frequency (in Hz) of the fiducial waveform, default: chosen + so that it starts before the analysis window + time_marginalization, phase_marginalization, calibration_marginalization: bool, optional + Not implemented for this likelihood; must be False + use_float32: bool, optional + If true, the summary data are stored and multiplied in single precision. + Default is :code:`False`. + + Returns + ======= + Likelihood: `bilby.core.likelihood.Likelihood` + A likelihood object, able to compute the likelihood of the data given + some model parameters. + + Notes + ===== + The time-domain heterodyning likelihood does not currently support time, + phase or calibration marginalization. + """ + + def __init__(self, interferometers, waveform_generator, fiducial_parameters=None, + maximization_kwargs=None, update_fiducial_parameters=False, + parameters_to_be_updated=("chirp_mass", "mass_ratio"), + distance_marginalization=False, priors=None, + distance_marginalization_lookup_table=None, reference_frame="sky", + time_reference="geocenter", epsilon=(0.5, 0.1), time_split=-0.02, + spacing=(2, 10, 100), spacing_times=(0.02, 0.1), chi=1, + fiducial_minimum_frequency=None, + time_marginalization=False, phase_marginalization=False, + calibration_marginalization=False, use_float32=False): + super(TimeDomainHeterodyningGravitationalWaveTransient, self).__init__( + interferometers=interferometers, + waveform_generator=waveform_generator, + distance_marginalization=distance_marginalization, + priors=priors, + distance_marginalization_lookup_table=distance_marginalization_lookup_table, + reference_frame=reference_frame, + time_reference=time_reference, + placement="model_time", + time_marginalization=time_marginalization, + phase_marginalization=phase_marginalization, + calibration_marginalization=calibration_marginalization, + ) + self.epsilon = epsilon + self.time_split = time_split + self.spacing = spacing + self.spacing_times = spacing_times + self.chi = chi + self.use_float32 = use_float32 + # this likelihood sets its own entries in waveform_arguments; copy the dictionary so that + # they do not reach other waveform generators that share it + self.waveform_generator.waveform_arguments = dict(self.waveform_generator.waveform_arguments) + self.requested_fiducial_minimum_frequency = fiducial_minimum_frequency + + if fiducial_parameters is None: + logger.info("Drawing fiducial parameters from prior.") + fiducial_parameters = priors.sample() + self.fiducial_parameters = dict(fiducial_parameters) + + if update_fiducial_parameters: + logger.info("Using scipy optimization to find maximum likelihood parameters.") + self.parameters_to_be_updated = list(parameters_to_be_updated) + logger.info(f"Parameters over which likelihood is maximized: {self.parameters_to_be_updated}") + self.parameter_bounds = self.get_bounds_from_priors(self.priors) + self.fiducial_parameters = self.find_maximum_likelihood_parameters( + self.parameter_bounds, maximization_kwargs=maximization_kwargs) + + if self.distance_marginalization: + self.fiducial_parameters["luminosity_distance"] = self._ref_dist + self.set_fiducial_waveforms(self.fiducial_parameters) + logger.info("Fiducial waveforms set up") + self.setup_bins() + self.compute_summary_data() + logger.info("Summary Data Obtained") + logger.info(f"Fiducial likelihood: {self.log_likelihood_ratio(self.fiducial_parameters):.2f}") + + def _model_times(self, interferometer, indices, parameters): + """ Times (in s) of the analysis-window samples `indices` relative to the arrival of the merger """ + # the time of the first sample, plus the samples `indices` (not all N samples) + first = times_from_merger( + interferometer.time_array[interferometer.analysis_slice][0], 1, + interferometer.sampling_frequency, parameters["geocent_time"], + interferometer.time_delay_from_geocenter(parameters["ra"], parameters["dec"], parameters["geocent_time"])) + return first[0] + np.asarray(indices) / interferometer.sampling_frequency + + def _modes_and_harmonics_at_combined_bin_edges(self, signal, model_times): + """ Modes and harmonic factors of `signal` at the model times of each + interferometer, from one waveform call (None if the waveform failed) + """ + all_modes_and_harmonics = signal["heterodyning_terms_at"](np.concatenate(model_times)) + if all_modes_and_harmonics is None: + return None + all_modes, all_harmonics = all_modes_and_harmonics + per_interferometer = list() + first = 0 + for interferometer_model_times in model_times: + interferometer_part = slice(first, first + len(interferometer_model_times)) + first += len(interferometer_model_times) + modes = {key: mode[interferometer_part] for key, mode in all_modes.items()} + # (Z, Zbar): arrays for the precessing models, constants otherwise + harmonics = { + key: tuple(np.asarray(factor)[interferometer_part] if np.ndim(factor) else factor + for factor in factors) + for key, factors in all_harmonics.items()} + per_interferometer.append((modes, harmonics)) + return per_interferometer + + def _r_tilde_coefficients(self, waveform_polarizations, interferometers, parameters): + """ The heterodyning coefficients of each interferometer in each bin (None if the waveform failed) """ + modes_and_harmonics = self._modes_and_harmonics_at_combined_bin_edges( + waveform_polarizations, + [self._model_times(interferometer, self.bin_edges[interferometer.name], parameters) + for interferometer in interferometers]) + if modes_and_harmonics is None: + return None + r_tilde = list() + for interferometer, (modes, harmonics) in zip(interferometers, modes_and_harmonics): + # as in Interferometer.get_time_domain_detector_response + antenna_time = parameters["geocent_time"] if interferometer.reference_time is None \ + else interferometer.reference_time + antenna_factor = ( # F = F_+ + i F_x + interferometer.antenna_response( + parameters["ra"], parameters["dec"], antenna_time, parameters["psi"], "plus") + + 1j * interferometer.antenna_response( + parameters["ra"], parameters["dec"], antenna_time, parameters["psi"], "cross")) + # r_tilde_lm = h_lm (F Z_lm + (F Zbar_lm)^*) / h^0_lm at the bin edges + r_tilde_at_edges = np.array([ + modes[mode] * (antenna_factor * harmonics[mode][0] + np.conj(antenna_factor * harmonics[mode][1])) + if mode in modes else np.zeros(len(self.bin_edges[interferometer.name])) + for mode in self.mode_keys]) * self._inverse_h0_at_edges[interferometer.name] + r_tilde_0 = (r_tilde_at_edges[:, 1:] + r_tilde_at_edges[:, :-1]) / 2 + r_tilde_1 = (r_tilde_at_edges[:, 1:] - r_tilde_at_edges[:, :-1]) / self.bin_widths[interferometer.name] + r_tilde.append(np.stack([r_tilde_0, r_tilde_1], axis=1)) + return r_tilde + + def set_fiducial_waveforms(self, parameters): + """ Compute the modes of the fiducial waveform at the samples of the analysis + window of each interferometer + + Parameters + ========== + parameters: dict + The fiducial parameters + """ + parameters = parameters.copy() + parameters.update(self.get_sky_frame_parameters(parameters=parameters)) + self._fiducial_times = { + interferometer.name: self._model_times( + interferometer, np.arange(len(self._data_times_inverse_covariance[interferometer.name])), parameters) + for interferometer in self.interferometers} + window_start = -min(times[0] for times in self._fiducial_times.values()) + # leading-order (2, 2) frequency 1.5 window_start + 0.1 s before the merger, minus 1 Hz + converted = self.waveform_generator.parameter_conversion(parameters.copy())[0] + chirp_mass_time = component_masses_to_chirp_mass(converted["mass_1"], converted["mass_2"]) \ + * solar_mass * gravitational_constant / speed_of_light ** 3 + frequency = (5 / 256 * chirp_mass_time ** (-5 / 3) / (1.5 * max(window_start, 0) + 0.1)) ** (3 / 8) / np.pi + derived = min(max(frequency - 1, 1.0), + self.waveform_generator.waveform_arguments.get("minimum_frequency", 20.0)) + requested = self.requested_fiducial_minimum_frequency + if requested is not None and requested > derived: + raise ValueError( + f"fiducial_minimum_frequency ({requested:g} Hz) must not be above {derived:.2f} Hz " + "(the lower of the waveforms' minimum_frequency and the frequency derived from the fiducial chirp mass " + "and the start of the analysis window)") + self.fiducial_minimum_frequency = derived if requested is None else requested + while True: + self._set_fiducial() + modes_and_harmonics = self._modes_and_harmonics_at_combined_bin_edges( + self.waveform_generator.time_domain_strain(parameters), + [self._fiducial_times[interferometer.name] for interferometer in self.interferometers]) + self._unset_fiducial() + if modes_and_harmonics is None: + raise ValueError(f"Cannot compute fiducial waveforms for {parameters}") + if all(np.any(modes[(2, 2)][:1]) for modes, _ in modes_and_harmonics): + break + if requested is not None or self.fiducial_minimum_frequency <= 1: + raise ValueError( + f"The fiducial waveform from {self.fiducial_minimum_frequency:.2f} Hz does not start " + f"before the analysis window for {parameters}") + self.fiducial_minimum_frequency = max(self.fiducial_minimum_frequency - 1, 1) + logger.info(f"Fiducial waveform from {self.fiducial_minimum_frequency:.2f} Hz") + self.fiducial_modes = dict() + self._fiducial_end = dict() + for interferometer, (modes, _) in zip(self.interferometers, modes_and_harmonics): + self.fiducial_modes[interferometer.name] = dict() + ends = list() + for key, mode in modes.items(): + nonzero = np.flatnonzero(mode) + if len(nonzero) == 0: # mode not in the fiducial waveform + continue + mode[nonzero[-1] + 1:] = mode[nonzero[-1]] # hold the last non-zero value after the ringdown + self.fiducial_modes[interferometer.name][key] = mode + ends.append(nonzero[-1]) + self._fiducial_end[interferometer.name] = max(ends) + self.mode_keys = sorted({key for modes in self.fiducial_modes.values() for key in modes}) + logger.debug(f"Modes of the heterodyning basis: {self.mode_keys}") + + @staticmethod + def time_domain_bin_edges(times, epsilon=(0.5, 0.1), time_split=-0.02, spacing=(2, 10, 100), + spacing_times=(0.02, 0.1), chi=1): + """ Bin edges for the time-domain heterodyning likelihood + + Before the merger, the bins are equally spaced in a phase proxy as in + Zackay et al. arXiv:1806.08792, with powers of time in place of powers of + frequency; after the merger, they have a fixed number of samples. See + Sharma et al. arXiv:2601.11239. + + Parameters + ========== + times: array_like + Increasing times relative to the merger (in s) + epsilon, time_split, spacing, spacing_times, chi: + See :code:`TimeDomainHeterodyningGravitationalWaveTransient` + + Returns + ======= + array_like: The indices of the bin edges in `times`, first and last included + """ + gamma = np.array([5 / 8, 3 / 8, 1 / 4, -3 / 8]) # powers of |t| in the phase proxy + edges = [0, len(times) - 1] + inspiral = np.flatnonzero(times < 0) + if len(inspiral) > 1: + t = np.abs(times[inspiral]) + d_alpha = 2 * np.pi * chi / np.where(gamma < 0, t.min(), t.max()) ** gamma + proxy = np.sum(np.sign(gamma) * d_alpha * t[:, None] ** gamma, axis=1) + for region, region_epsilon in zip((times[inspiral] <= time_split, times[inspiral] > time_split), epsilon): + if np.sum(region) < 2: + continue + phase = np.abs(proxy[region] - proxy[region][0]) + number_of_bins = max(int(phase[-1] / region_epsilon), 1) + targets = np.linspace(0, phase[-1], number_of_bins + 1) + edges.extend(inspiral[region][np.searchsorted(phase, targets, side="right") - 1]) + merger = inspiral[-1] if len(inspiral) else 0 + boundaries = [merger] + [max(merger, np.searchsorted(times, x)) for x in spacing_times] + [len(times)] + for start, stop, step in zip(boundaries[:-1], boundaries[1:], spacing): + edges.extend(range(start, stop, step)) + return np.unique(edges) + + def setup_bins(self): + """ Set up the bin edges of each interferometer from its fiducial waveform """ + self.bin_edges = dict() + self.bin_widths = dict() + self.number_of_bins = dict() + for interferometer in self.interferometers: + name = interferometer.name + times = self._fiducial_times[name] + end = self._fiducial_end[name] + edges = self.time_domain_bin_edges( + times[:end + 1], + epsilon=self.epsilon, + time_split=self.time_split, + spacing=self.spacing, + spacing_times=self.spacing_times, + chi=self.chi) + if end < len(times) - 1: + edges = np.append(edges, len(times) - 1) + self.bin_edges[name] = edges + self.bin_widths[name] = np.diff(times[edges]) + self.number_of_bins[name] = len(edges) - 1 + logger.info( + f"{name}: {self.number_of_bins[name]} bins between {times[edges[0]]:.3f} s " + f"and {times[edges[-1]]:.3f} s from the merger") + + def compute_summary_data(self): + """ Compute the summary data of each interferometer from its fiducial waveform """ + self.summary_data = dict() + self._inverse_h0_at_edges = dict() + for interferometer in self.interferometers: + name = interferometer.name + edges = self.bin_edges[name] # sample indices of the bin edges: first and last sample included + t = self._fiducial_times[name] # t_j + number_of_samples, number_of_bins = len(t), len(edges) - 1 + samples_per_bin = np.diff(edges) + samples_per_bin[-1] += 1 + bin_of_sample = np.repeat(np.arange(number_of_bins), samples_per_bin) + bin_centres = (t[edges[:-1]] + t[edges[1:]]) / 2 # t_b + time_from_bin_centre = t - bin_centres[bin_of_sample] # t_j - t_b + + h0 = np.array([self.fiducial_modes[name][mode] for mode in self.mode_keys]) # h^0_lm(t_j) + self._inverse_h0_at_edges[name] = 1 / h0[:, edges] + + # basis functions h^0_lm (t_j - t_b)^p as real functions, in the order (Re/-Im, lm, p); + # basis[i] restricted to bin b is the column (i, b) of V + h0_times_power_of_time = np.stack([h0, h0 * time_from_bin_centre], axis=1) + basis = np.concatenate([h0_times_power_of_time.real, -h0_times_power_of_time.imag]) + basis = basis.reshape(-1, number_of_samples) + number_of_basis = len(basis) + + # a = V^T C^-1 d: basis times C^-1 d, summed over the samples of each bin + a = np.add.reduceat(basis * self._data_times_inverse_covariance[name], edges[:-1], axis=1) + + # S = V^T C^-1 V without the terms quadratic in the slopes r_tilde_1 (p = p' = 1, odd i): + # C^-1 times each column (i, b) of V, then times basis[j], summed over the samples of + # each bin b'; S is symmetric, so only j >= i is computed + x, y = interferometer.gohberg_semencul_vectors + S = np.zeros((number_of_basis, number_of_bins, number_of_basis, number_of_bins)) + for i in range(number_of_basis): + columns = np.zeros((number_of_bins, number_of_samples)) + columns[bin_of_sample, np.arange(number_of_samples)] = basis[i] + C_inv_columns = np.array([gohberg_semencul_product(x, y, column) for column in columns]).T + for j in range(i, number_of_basis): + block = np.add.reduceat(basis[j][:, None] * C_inv_columns, edges[:-1], axis=0) + S[j, :, i, :] = block + S[i, :, j, :] = block.T + + # S = sum_n lambda_n e_n e_n^T, keeping the eigenvalues above 1e-12 of the largest + size = number_of_basis * number_of_bins + eigenvalues, eigenvectors = np.linalg.eigh(S.reshape(size, size)) + keep = np.abs(eigenvalues) > 1e-12 * np.max(np.abs(eigenvalues)) + + self.summary_data[name] = dict( + a=a.ravel(), S_eigenvalues=eigenvalues[keep], + S_eigenvectors=np.ascontiguousarray( + eigenvectors[:, keep].T, dtype=np.float32 if self.use_float32 else np.float64), + bin_of_sample=bin_of_sample, time_from_bin_centre=time_from_bin_centre) + logger.info(f"{name}: {len(self.mode_keys)} modes, {number_of_bins} bins, " + f"S of size {size} and rank {np.sum(keep)}") + + def calculate_snrs(self, waveform_polarizations, interferometer, *, return_array=True, parameters, + r_tilde=None): + """ Compute the snrs + + Parameters + ========== + waveform_polarizations: dict + The output of the heterodyning source model + interferometer: bilby.gw.detector.Interferometer + The bilby interferometer object + return_array: bool + Not used; kept for the :code:`GravitationalWaveTransient` interface + parameters: dict + The parameters of the signal + r_tilde: array_like, optional + The heterodyning coefficients of this interferometer, if already computed + + Returns + ======= + calculated_snrs: _CalculatedSNRs + An object containing the SNR quantities (None if the waveform failed) + """ + if r_tilde is None: + r_tilde = self._r_tilde_coefficients(waveform_polarizations, [interferometer], parameters) + if r_tilde is None: + return None + r_tilde = r_tilde[0] + summary = self.summary_data[interferometer.name] + v = np.concatenate([r_tilde.real.ravel(), r_tilde.imag.ravel()]) + d_inner_s = float(summary["a"] @ v) + # e_n . v in the precision of the stored eigenvectors (see use_float32) + projections = summary["S_eigenvectors"] @ v.astype(summary["S_eigenvectors"].dtype, copy=False) + s_inner_s = float(summary["S_eigenvalues"] @ projections.astype(np.float64) ** 2) + return self._CalculatedSNRs( + d_inner_h=d_inner_s, + optimal_snr_squared=s_inner_s, + complex_matched_filter_snr=d_inner_s / s_inner_s ** 0.5 if s_inner_s > 0 else 0.0, + d_inner_h_array=None, + optimal_snr_squared_array=None, + ) + + def log_likelihood_ratio(self, parameters): + parameters = copy.deepcopy(parameters) + parameters.update(self.get_sky_frame_parameters(parameters)) + signal = self.waveform_generator.time_domain_strain(parameters) + r_tilde = None if signal is None else self._r_tilde_coefficients(signal, self.interferometers, parameters) + if r_tilde is None: + return np.nan_to_num(-np.inf) + total_snrs = self._CalculatedSNRs() + for interferometer, interferometer_r_tilde in zip(self.interferometers, r_tilde): + total_snrs += self.calculate_snrs( + signal, interferometer, parameters=parameters, r_tilde=interferometer_r_tilde) + log_l = self.compute_log_likelihood_from_snrs(total_snrs, parameters=parameters) + return float(np.real(log_l)) + + def _compute_full_waveform(self, signal_polarizations, interferometer, parameters): + """ The heterodyned template at the samples of the analysis window """ + r_tilde = self._r_tilde_coefficients(signal_polarizations, [interferometer], parameters) + if r_tilde is None: + return None + summary = self.summary_data[interferometer.name] + r_tilde_at_samples = r_tilde[0][:, :, summary["bin_of_sample"]] + h0 = np.array([self.fiducial_modes[interferometer.name][mode] for mode in self.mode_keys]) + return np.sum(((r_tilde_at_samples[:, 0] + r_tilde_at_samples[:, 1] * summary["time_from_bin_centre"]) + * h0).real, axis=0) + + def _rescale_signal(self, signal, new_distance): + factor = self._ref_dist / new_distance + terms_at = signal["heterodyning_terms_at"] + + def rescaled(times): + terms = terms_at(times) + if terms is None: + return None + modes, harmonics = terms + return {key: mode * factor for key, mode in modes.items()}, harmonics + + signal["heterodyning_terms_at"] = rescaled + + def find_maximum_likelihood_parameters(self, parameter_bounds, maximization_kwargs=None): + """ Find the maximum of the time-domain likelihood over + :code:`parameters_to_be_updated`, using scipy differential evolution + + Parameters + ========== + parameter_bounds: list + The bounds of the parameters + maximization_kwargs: dict, optional + Keyword arguments for :code:`scipy.optimize.differential_evolution` + + Returns + ======= + dict: The updated fiducial parameters + """ + if maximization_kwargs is None: + maximization_kwargs = dict() + exact_likelihood = self.exact_likelihood() + + def negative_log_likelihood_ratio(parameter_list): + return -exact_likelihood.log_likelihood_ratio(self.get_parameter_dictionary_from_list(parameter_list)) + + start = self.get_parameter_list_from_dictionary(self.fiducial_parameters) + logger.info(f"Initial ln likelihood ratio: {-negative_log_likelihood_ratio(start):.2f}") + output = differential_evolution( + negative_log_likelihood_ratio, bounds=parameter_bounds, x0=start, **maximization_kwargs) + updated_parameters = self.get_parameter_dictionary_from_list(output["x"]) + logger.info(f"Maximum ln likelihood ratio: {-output['fun']:.2f} at " + f"{dict(zip(self.parameters_to_be_updated, output['x']))}") + return updated_parameters + + def exact_likelihood(self): + """ The exact time-domain likelihood (model time) with the same data, waveform + model, priors and settings as this likelihood + + Returns + ======= + bilby.gw.likelihood.TimeDomainGravitationalWaveTransient + """ + from ..source import lal_binary_black_hole_model_time + + waveform_arguments = { + key: value for key, value in self.waveform_generator.waveform_arguments.items() + if key not in ("fiducial", "fiducial_minimum_frequency")} + return TimeDomainGravitationalWaveTransient( + interferometers=self.interferometers, + waveform_generator=self.waveform_generator.__class__( + duration=self.waveform_generator.duration, + sampling_frequency=self.waveform_generator.sampling_frequency, + start_time=self.waveform_generator.start_time, + time_domain_source_model=lal_binary_black_hole_model_time, + parameter_conversion=self.waveform_generator.parameter_conversion, + waveform_arguments=waveform_arguments), + distance_marginalization=self.distance_marginalization, + priors=self.priors, + distance_marginalization_lookup_table=( + self.cached_lookup_table_filename if self.distance_marginalization else None), + reference_frame=self.reference_frame, + time_reference=self.time_reference, + placement="model_time") + + def compare_with_exact_likelihood(self, samples, outdir=".", label="heterodyning", n_samples=None): + """ Compare the heterodyned likelihood with the exact time-domain likelihood + (model time) at a set of samples, e.g. the posterior of a run + + Saves two figures in `outdir`: `{label}_heterodyning_vs_exact.png` (both log + likelihood ratios against the chirp mass, heterodyned against exact, and their + relative difference against the chirp mass) and + `{label}_heterodyning_differences.png` (heterodyned minus exact against the + other parameters). + + Parameters + ========== + samples: pandas.DataFrame or dict + The samples, e.g. :code:`result.posterior` + outdir: str + The directory for the figures + label: str + The label used in the file names + n_samples: int, optional + The number of samples drawn at random from `samples` (default: all) + + Returns + ======= + dict: The samples used and the exact and heterodyned log likelihood ratios + """ + import matplotlib.pyplot as plt + import pandas as pd + + from ...core.utils import safe_save_figure + + samples = pd.DataFrame(samples) + if n_samples is not None and n_samples < len(samples): + samples = samples.sample(n_samples) + exact_likelihood = self.exact_likelihood() + exact, heterodyned = np.array([ + (exact_likelihood.log_likelihood_ratio(dict(sample)), self.log_likelihood_ratio(dict(sample))) + for _, sample in samples.iterrows()]).T + difference = heterodyned - exact + logger.info(f"Heterodyned - exact ln likelihood ratio at {len(samples)} samples: median " + f"{np.median(difference):.3g}, maximum absolute value {np.max(np.abs(difference)):.3g}") + if "chirp_mass" in samples: + chirp_mass = samples["chirp_mass"].to_numpy() + else: + chirp_mass = component_masses_to_chirp_mass(samples["mass_1"], samples["mass_2"]).to_numpy() + + fig, axes = plt.subplots(1, 3, figsize=(16, 4.5)) + axes[0].plot(chirp_mass, exact, ".", label="exact (model time)") + axes[0].plot(chirp_mass, heterodyned, ".", alpha=0.5, label="heterodyned") + axes[0].set_xlabel("chirp_mass") + axes[0].set_ylabel("log likelihood ratio") + axes[0].legend() + limits = [min(exact.min(), heterodyned.min()), max(exact.max(), heterodyned.max())] + axes[1].plot(exact, heterodyned, ".") + axes[1].plot(limits, limits, "k--", label="y = x") + axes[1].set_xlabel("exact log likelihood ratio") + axes[1].set_ylabel("heterodyned log likelihood ratio") + axes[1].legend() + axes[2].plot(chirp_mass, difference / np.abs(exact), ".") + axes[2].axhline(0, color="k", ls="--") + axes[2].set_xlabel("chirp_mass") + axes[2].set_ylabel("(heterodyned - exact) / |exact|") + fig.tight_layout() + safe_save_figure(fig, f"{outdir}/{label}_heterodyning_vs_exact.png") + plt.close(fig) + + keys = ["mass_ratio", "chi_1", "chi_2", "a_1", "a_2", "tilt_1", "tilt_2", "theta_jn", "geocent_time"] + if not self.distance_marginalization: + keys.append("luminosity_distance") + keys = [key for key in keys if key in samples and samples[key].nunique() > 1] + if len(keys) > 0: + columns = min(len(keys), 4) + rows = int(np.ceil(len(keys) / columns)) + fig, axes = plt.subplots(rows, columns, figsize=(4 * columns, 3.5 * rows), squeeze=False) + for ax, key in zip(axes.flat, keys): + if key == "geocent_time": + ax.plot(samples[key] - self.fiducial_parameters["geocent_time"], difference, ".") + ax.set_xlabel("geocent_time - fiducial geocent_time (s)") + else: + ax.plot(samples[key], difference, ".") + ax.set_xlabel(key) + ax.axhline(0, color="k", ls="--") + ax.set_ylabel("heterodyned - exact") + for ax in list(axes.flat)[len(keys):]: + ax.set_visible(False) + fig.tight_layout() + safe_save_figure(fig, f"{outdir}/{label}_heterodyning_differences.png") + plt.close(fig) + return dict(samples=samples, exact=exact, heterodyned=heterodyned) + + def _set_fiducial(self): + """ Make the waveform generator generate the fiducial waveform """ + self.waveform_generator.waveform_arguments["fiducial"] = 1 + self.waveform_generator.waveform_arguments["fiducial_minimum_frequency"] = self.fiducial_minimum_frequency + self.waveform_generator._cache["parameters"] = None + + def _unset_fiducial(self): + """ Make the waveform generator generate waveforms at sampled parameters """ + self.waveform_generator.waveform_arguments["fiducial"] = 0 + self.waveform_generator._cache["parameters"] = None + + get_parameter_dictionary_from_list = RelativeBinningGravitationalWaveTransient.get_parameter_dictionary_from_list + get_parameter_list_from_dictionary = RelativeBinningGravitationalWaveTransient.get_parameter_list_from_dictionary + get_bounds_from_priors = RelativeBinningGravitationalWaveTransient.get_bounds_from_priors + + @property + def meta_data(self): + meta_data = super(TimeDomainHeterodyningGravitationalWaveTransient, self).meta_data + meta_data.update( + fiducial_parameters=self.fiducial_parameters, + parameters_to_be_updated=getattr(self, "parameters_to_be_updated", None), + fiducial_minimum_frequency=self.fiducial_minimum_frequency, + mode_keys=self.mode_keys, + epsilon=self.epsilon, + time_split=self.time_split, + spacing=self.spacing, + spacing_times=self.spacing_times, + chi=self.chi, + number_of_bins=self.number_of_bins, + bin_edges=self.bin_edges, + bin_edge_times={name: self._fiducial_times[name][edges] for name, edges in self.bin_edges.items()}) + return meta_data diff --git a/bilby/gw/model_time_waveforms.py b/bilby/gw/model_time_waveforms.py new file mode 100644 index 000000000..8f6fdfa3b --- /dev/null +++ b/bilby/gw/model_time_waveforms.py @@ -0,0 +1,182 @@ +""" +IMRPhenomT-family waveforms evaluated at arbitrary times. + +Wrappers of the lalsimulation functions that evaluate IMRPhenomT +(Estellés et al. arXiv:2004.08302), IMRPhenomTHM (arXiv:2012.11923), +IMRPhenomTP and IMRPhenomTPHM (arXiv:2105.05872) at given model times +t / M, where M is the total mass in seconds and t = 0 at the peak of the (2, 2) +amplitude. These functions require a lalsuite build that provides them. +""" +import numpy as np + +APPROXIMANTS = ("IMRPhenomT", "IMRPhenomTHM", "IMRPhenomTP", "IMRPhenomTPHM") +PRECESSING_APPROXIMANTS = ("IMRPhenomTP", "IMRPhenomTPHM") + +# lalsimulation function names, e.g. SimIMRPhenomTPHMModelTimes(Truncated) +_POLARIZATION_FUNCTION = "Sim{approximant}ModelTimes" +_MODE_FUNCTION = "Sim{approximant}ModesModelTimes" +_COPRECESSING_HARMONICS_FUNCTION = "Sim{approximant}CoprecModesHarmonicsModelTimes" + + +def _get_approximant(template, approximant, truncate): + import lalsimulation as lalsim + name = template.format(approximant=approximant) + ("Truncated" if truncate else "") + try: + return getattr(lalsim, name) + except AttributeError: + raise ImportError( + f"lalsimulation has no {name}: model-time waveforms need a lalsuite build " + "with the time-array IMRPhenomT-family functions") from None + + +def _time_sequence(model_times): + import lal + model_times = np.asarray(model_times, dtype=float) + sequence = lal.CreateREAL8Sequence(len(model_times)) + sequence.data = model_times + return sequence + + +def spin_arguments(approximant, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z): + """ The spin arguments of the lalsimulation model-time functions: all six + components for IMRPhenomTP and IMRPhenomTPHM, the z components for + IMRPhenomT and IMRPhenomTHM + + Raises a ValueError for other approximants, or for in-plane spins with a + non-precessing approximant. + """ + if approximant not in APPROXIMANTS: + raise ValueError( + f"model-time waveforms are available for {', '.join(APPROXIMANTS)}, not {approximant!r}") + if approximant in PRECESSING_APPROXIMANTS: + return (spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z) + if spin_1x != 0 or spin_1y != 0 or spin_2x != 0 or spin_2y != 0: + raise ValueError( + f"Non-zero transverse spins were given, but {approximant} is a non-precessing approximant.") + return (spin_1z, spin_2z) + + +def polarizations_at_model_times( + model_times, approximant, mass_1, mass_2, spin_1x, spin_1y, spin_1z, + spin_2x, spin_2y, spin_2z, luminosity_distance, inclination, phase, + delta_time, minimum_frequency, reference_frequency, waveform_dictionary, + truncate=True): + """ The plus and cross polarizations at given model times + + Parameters + ========== + model_times: array_like + The times t / M, with t = 0 at the peak of the (2, 2) amplitude + approximant: str + IMRPhenomT, IMRPhenomTHM, IMRPhenomTP or IMRPhenomTPHM + mass_1, mass_2: float + The component masses (in kg) + spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z: float + The dimensionless spin components in the lalsimulation frame + luminosity_distance: float + The luminosity distance (in m) + inclination: float + The inclination (in rad) + phase: float + The reference phase (in rad) + delta_time: float + The sampling interval (in s) of the waveform on which the precession angles are computed + minimum_frequency: float + The minimum frequency (in Hz) + reference_frequency: float + The reference frequency (in Hz) + waveform_dictionary: lal.Dict + The lal waveform dictionary + truncate: bool + If True, the waveform is zero at times outside its span, from the + minimum frequency to the end of the ringdown + + Returns + ======= + dict: The plus and cross polarizations + """ + spins = spin_arguments(approximant, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z) + function = _get_approximant(_POLARIZATION_FUNCTION, approximant, truncate) + hplus, hcross = function( + mass_1, mass_2, *spins, luminosity_distance, inclination, delta_time, + minimum_frequency, reference_frequency, phase, _time_sequence(model_times), + waveform_dictionary) + return dict(plus=hplus.data.data, cross=hcross.data.data) + + +def modes_at_model_times( + model_times, approximant, mass_1, mass_2, spin_1x, spin_1y, spin_1z, + spin_2x, spin_2y, spin_2z, luminosity_distance, phase, delta_time, + minimum_frequency, reference_frequency, waveform_dictionary, truncate=True): + """ The spherical-harmonic modes at given model times; for IMRPhenomTP and + IMRPhenomTPHM, in the frame of the orbital angular momentum at the reference + frequency + + Parameters are as in :code:`polarizations_at_model_times`, without the inclination. + + Returns + ======= + dict: The modes, with keys (l, m) + """ + spins = spin_arguments(approximant, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z) + function = _get_approximant(_MODE_FUNCTION, approximant, truncate) + hlms = function( + mass_1, mass_2, *spins, luminosity_distance, delta_time, minimum_frequency, + reference_frequency, phase, _time_sequence(model_times), waveform_dictionary) + return _modes_dict(hlms) + + +def _modes_dict(hlms, nonnegative=False): + modes = dict() + while hlms is not None: + if not nonnegative or hlms.m >= 0: + modes[(hlms.l, hlms.m)] = hlms.mode.data.data.copy() + hlms = hlms.next + return modes + + +def heterodyning_terms_at_model_times( + model_times, approximant, mass_1, mass_2, spin_1x, spin_1y, spin_1z, + spin_2x, spin_2y, spin_2z, luminosity_distance, inclination, phase, + delta_time, minimum_frequency, reference_frequency, waveform_dictionary, + truncate=True): + """ The modes and harmonic factors used by + :code:`bilby.gw.likelihood.TimeDomainHeterodyningGravitationalWaveTransient` + at given model times: the modes with m > 0 for IMRPhenomT and IMRPhenomTHM, + the co-precessing modes with m >= 0 for IMRPhenomTP and IMRPhenomTPHM. See + Sharma et al. arXiv:2601.11239. + + Parameters are as in :code:`polarizations_at_model_times`. + + Returns + ======= + modes: dict + The modes, with keys (l, m) + harmonics: dict + The harmonic factors of each mode and of its complex conjugate, with keys (l, m) + """ + spins = spin_arguments(approximant, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z) + if approximant not in PRECESSING_APPROXIMANTS: + import lal + modes = modes_at_model_times( + model_times, approximant, mass_1, mass_2, spin_1x, spin_1y, spin_1z, + spin_2x, spin_2y, spin_2z, luminosity_distance, phase, delta_time, + minimum_frequency, reference_frequency, waveform_dictionary, truncate=truncate) + modes = {(ell, m): mode for (ell, m), mode in modes.items() if m > 0} + # Y_lm = -2Y_lm(inclination, pi / 2 - phase) + harmonics = { + (ell, m): (lal.SpinWeightedSphericalHarmonic(inclination, np.pi / 2 - phase, -2, ell, m), + (-1) ** ell * lal.SpinWeightedSphericalHarmonic(inclination, np.pi / 2 - phase, -2, ell, -m)) + for ell, m in modes} + return modes, harmonics + + function = _get_approximant(_COPRECESSING_HARMONICS_FUNCTION, approximant, truncate) + hlms, factors = function( + mass_1, mass_2, *spins, luminosity_distance, inclination, delta_time, minimum_frequency, + reference_frequency, phase, _time_sequence(model_times), waveform_dictionary) + modes = {key: mode for key, mode in _modes_dict(hlms, nonnegative=True).items() if np.any(mode)} + factors = _modes_dict(factors) + harmonics = { + (ell, m): (factors[(ell, m)], (-1) ** ell * factors[(ell, -m)] if m > 0 else 0) + for ell, m in modes} + return modes, harmonics diff --git a/bilby/gw/source.py b/bilby/gw/source.py index 97ad51b51..a937cb584 100644 --- a/bilby/gw/source.py +++ b/bilby/gw/source.py @@ -8,6 +8,7 @@ lalsim_SimInspiralFD, lalsim_SimInspiralChooseFDWaveform, lalsim_SimInspiralChooseFDWaveformSequence, + lalsim_SimInspiralChooseTDWaveform, safe_cast_mode_to_int) UNUSED_KWARGS_MESSAGE = """There are unused waveform kwargs. @@ -690,6 +691,367 @@ def _base_lal_cbc_fd_waveform( return dict(plus=h_plus, cross=h_cross) +def lal_binary_black_hole_time_domain( + time_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, + phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs): + """ A time-domain Binary Black Hole waveform model using lalsimulation, for + :code:`bilby.gw.likelihood.TimeDomainGravitationalWaveTransient` + + Parameters + ========== + time_array: array_like + The times of the data (only the sample spacing is used) + mass_1: float + The mass of the heavier object in solar masses + mass_2: float + The mass of the lighter object in solar masses + luminosity_distance: float + The luminosity distance in megaparsec + a_1: float + Dimensionless primary spin magnitude + tilt_1: float + Primary tilt angle + phi_12: float + Azimuthal angle between the two component spins + a_2: float + Dimensionless secondary spin magnitude + tilt_2: float + Secondary tilt angle + phi_jl: float + Azimuthal angle between the total binary angular momentum and the + orbital angular momentum + theta_jn: float + Angle between the total binary angular momentum and the line of sight + phase: float + The phase at reference frequency or peak amplitude (depends on waveform) + kwargs: dict + Optional keyword arguments + Supported arguments: + + - waveform_approximant + - reference_frequency + - minimum_frequency + - catch_waveform_errors + - pn_spin_order + - pn_tidal_order + - pn_phase_order + - pn_amplitude_order + - mode_array: see :code:`lal_binary_black_hole` + - lal_waveform_dictionary + - turn_on_window: the alpha of the one-sided Tukey window applied to the + start of the waveform, for placement "nearest" or "fd_shift" (default 0: no window) + + Returns + ======= + dict: A dictionary with the plus and cross polarisation strain modes on the + waveform's own time grid, and the epoch, the time of the first sample relative + to the merger (in s) + """ + waveform_kwargs = dict( + waveform_approximant='IMRPhenomTP', reference_frequency=50.0, + minimum_frequency=20.0, catch_waveform_errors=False, + pn_spin_order=-1, pn_tidal_order=-1, pn_phase_order=-1, + pn_amplitude_order=0, turn_on_window=0) + waveform_kwargs.update(kwargs) + return _base_lal_cbc_td_waveform( + time_array=time_array, mass_1=mass_1, mass_2=mass_2, + luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase, + a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_12=phi_12, + phi_jl=phi_jl, **waveform_kwargs) + + +def lal_binary_neutron_star_time_domain( + time_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, + phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, lambda_1, lambda_2, + **kwargs): + """ A time-domain Binary Neutron Star waveform model using lalsimulation, for + :code:`bilby.gw.likelihood.TimeDomainGravitationalWaveTransient` + + Parameters + ========== + time_array: array_like + The times of the data (only the sample spacing is used) + mass_1: float + The mass of the heavier object in solar masses + mass_2: float + The mass of the lighter object in solar masses + luminosity_distance: float + The luminosity distance in megaparsec + a_1: float + Dimensionless primary spin magnitude + tilt_1: float + Primary tilt angle + phi_12: float + Azimuthal angle between the two component spins + a_2: float + Dimensionless secondary spin magnitude + tilt_2: float + Secondary tilt angle + phi_jl: float + Azimuthal angle between the total binary angular momentum and the + orbital angular momentum + theta_jn: float + Angle between the total binary angular momentum and the line of sight + phase: float + The phase at reference frequency or peak amplitude (depends on waveform) + lambda_1: float + Dimensionless tidal deformability of mass_1 + lambda_2: float + Dimensionless tidal deformability of mass_2 + kwargs: dict + Optional keyword arguments, as for :code:`lal_binary_black_hole_time_domain` + + Returns + ======= + dict: A dictionary with the plus and cross polarisation strain modes on the + waveform's own time grid, and the epoch + """ + waveform_kwargs = dict( + waveform_approximant='IMRPhenomPv2_NRTidal', reference_frequency=50.0, + minimum_frequency=20.0, catch_waveform_errors=False, + pn_spin_order=-1, pn_tidal_order=-1, pn_phase_order=-1, + pn_amplitude_order=0, turn_on_window=0) + waveform_kwargs.update(kwargs) + return _base_lal_cbc_td_waveform( + time_array=time_array, mass_1=mass_1, mass_2=mass_2, + luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase, + a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_12=phi_12, + phi_jl=phi_jl, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs) + + +def _base_lal_cbc_td_waveform( + time_array, mass_1, mass_2, luminosity_distance, theta_jn, phase, + a_1=0.0, a_2=0.0, tilt_1=0.0, tilt_2=0.0, phi_12=0.0, phi_jl=0.0, + lambda_1=0.0, lambda_2=0.0, eccentricity=0.0, **waveform_kwargs): + """ Generate a time-domain cbc waveform model using lalsimulation + + Parameters are as in :code:`_base_lal_cbc_fd_waveform`, with `time_array` + in place of `frequency_array`. + + Returns + ======= + dict: A dictionary with the plus and cross polarisation strain modes on the + waveform's own time grid, and the epoch, the time of the first sample relative + to the merger (in s) + """ + import lalsimulation as lalsim + from scipy.signal.windows import tukey + + waveform_approximant = waveform_kwargs.pop('waveform_approximant') + reference_frequency = waveform_kwargs.pop('reference_frequency') + minimum_frequency = waveform_kwargs.pop('minimum_frequency') + catch_waveform_errors = waveform_kwargs.pop('catch_waveform_errors') + turn_on_window = waveform_kwargs.pop('turn_on_window') + pn_amplitude_order = waveform_kwargs['pn_amplitude_order'] + + waveform_dictionary = set_waveform_dictionary(waveform_kwargs, lambda_1, lambda_2) + approximant = lalsim_GetApproximantFromString(waveform_approximant) + + if pn_amplitude_order != 0: + start_frequency = lalsim.SimInspiralfLow2fStart( + float(minimum_frequency), int(pn_amplitude_order), approximant + ) + else: + start_frequency = minimum_frequency + + delta_time = float(time_array[1] - time_array[0]) + + luminosity_distance = luminosity_distance * 1e6 * utils.parsec + mass_1 = mass_1 * utils.solar_mass + mass_2 = mass_2 * utils.solar_mass + + iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = bilby_to_lalsimulation_spins( + theta_jn=theta_jn, phi_jl=phi_jl, tilt_1=tilt_1, tilt_2=tilt_2, + phi_12=phi_12, a_1=a_1, a_2=a_2, mass_1=mass_1, mass_2=mass_2, + reference_frequency=reference_frequency, phase=phase) + + longitude_ascending_nodes = 0.0 + mean_per_ano = 0.0 + + try: + hplus, hcross = lalsim_SimInspiralChooseTDWaveform( + mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, + spin_2z, luminosity_distance, iota, phase, + longitude_ascending_nodes, eccentricity, mean_per_ano, delta_time, + start_frequency, reference_frequency, + waveform_dictionary, approximant) + except Exception as e: + if not catch_waveform_errors: + raise + else: + EDOM = (e.args[0] == 'Internal function call failed: Input domain error') + if EDOM: + failed_parameters = dict(mass_1=mass_1, mass_2=mass_2, + spin_1=(spin_1x, spin_1y, spin_1z), + spin_2=(spin_2x, spin_2y, spin_2z), + luminosity_distance=luminosity_distance, + iota=iota, phase=phase, + eccentricity=eccentricity, + start_frequency=start_frequency) + logger.warning("Evaluating the waveform failed with error: {}\n".format(e) + + "The parameters were {}\n".format(failed_parameters) + + "Likelihood will be set to -inf.") + return None + else: + raise + + if len(waveform_kwargs) > 0: + raise ValueError(UNUSED_KWARGS_MESSAGE.format(waveform_kwargs=waveform_kwargs)) + + window = tukey(hplus.data.length, turn_on_window) + window[hplus.data.length // 2:] = 1 + epoch = hplus.epoch.gpsSeconds + hplus.epoch.gpsNanoSeconds * 1e-9 + return dict(plus=window * hplus.data.data, cross=window * hcross.data.data, epoch=epoch) + + +def lal_binary_black_hole_model_time( + time_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, + phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs): + """ A time-domain Binary Black Hole waveform model for the IMRPhenomT family, + evaluated at the data samples of each interferometer, for + :code:`bilby.gw.likelihood.TimeDomainGravitationalWaveTransient` with + :code:`placement="model_time"`. See :code:`bilby.gw.model_time_waveforms`. + + Parameters + ========== + time_array: array_like + The times of the data (only the sample spacing is used) + mass_1: float + The mass of the heavier object in solar masses + mass_2: float + The mass of the lighter object in solar masses + luminosity_distance: float + The luminosity distance in megaparsec + a_1: float + Dimensionless primary spin magnitude + tilt_1: float + Primary tilt angle + phi_12: float + Azimuthal angle between the two component spins + a_2: float + Dimensionless secondary spin magnitude + tilt_2: float + Secondary tilt angle + phi_jl: float + Azimuthal angle between the total binary angular momentum and the + orbital angular momentum + theta_jn: float + Angle between the total binary angular momentum and the line of sight + phase: float + The phase at reference frequency or peak amplitude (depends on waveform) + kwargs: dict + Optional keyword arguments + Supported arguments: + + - waveform_approximant: IMRPhenomT, IMRPhenomTHM, IMRPhenomTP or IMRPhenomTPHM + - reference_frequency + - minimum_frequency + - catch_waveform_errors + - mode_array: see :code:`lal_binary_black_hole` + - lal_waveform_dictionary + + Returns + ======= + dict: polarizations_at, a function of the times relative to the merger (in s) + which returns the plus and cross polarizations (None if the waveform failed) + """ + from .model_time_waveforms import polarizations_at_model_times + + return dict(polarizations_at=_model_time_function( + polarizations_at_model_times, time_array, mass_1, mass_2, luminosity_distance, + a_1, tilt_1, phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs)) + + +def lal_binary_black_hole_time_domain_heterodyning( + time_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, + phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs): + """ A time-domain Binary Black Hole waveform model for the IMRPhenomT family, for + :code:`bilby.gw.likelihood.TimeDomainHeterodyningGravitationalWaveTransient`. See + :code:`bilby.gw.model_time_waveforms.heterodyning_terms_at_model_times`. + + Parameters and keyword arguments are as in + :code:`lal_binary_black_hole_model_time`, and: + + - fiducial: if 1, the waveform starts at fiducial_minimum_frequency and + polarizations_at is also returned (as in + :code:`lal_binary_black_hole_model_time`), so that the same waveform + generator can be used for an injection; otherwise only heterodyning_terms_at, + starting at minimum_frequency. Set by the likelihood. + - fiducial_minimum_frequency (default: minimum_frequency), set by the likelihood + + Returns + ======= + dict: heterodyning_terms_at, a function of the times relative to the merger + (in s) which returns the modes and harmonic factors (None if the waveform failed), + and, if fiducial is 1, polarizations_at + """ + from .model_time_waveforms import heterodyning_terms_at_model_times, polarizations_at_model_times + + fiducial = kwargs.pop("fiducial", 0) + fiducial_minimum_frequency = kwargs.pop("fiducial_minimum_frequency", None) + if fiducial and fiducial_minimum_frequency is not None: + kwargs["minimum_frequency"] = fiducial_minimum_frequency + arguments = (time_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, phi_12, a_2, + tilt_2, phi_jl, theta_jn, phase) + waveform = dict(heterodyning_terms_at=_model_time_function( + heterodyning_terms_at_model_times, *arguments, **kwargs)) + if fiducial: + waveform["polarizations_at"] = _model_time_function( + polarizations_at_model_times, *arguments, **kwargs) + return waveform + + +def _model_time_function( + function, time_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, + phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs): + """ The function of the times relative to the merger (in s) which evaluates + `function` of :code:`bilby.gw.model_time_waveforms` with these parameters + """ + import lal + from .model_time_waveforms import spin_arguments + + waveform_kwargs = dict( + waveform_approximant='IMRPhenomTPHM', reference_frequency=50.0, + minimum_frequency=20.0, catch_waveform_errors=False, + pn_spin_order=-1, pn_tidal_order=-1, pn_phase_order=-1, + pn_amplitude_order=0) + waveform_kwargs.update(kwargs) + waveform_approximant = waveform_kwargs.pop('waveform_approximant') + reference_frequency = waveform_kwargs.pop('reference_frequency') + minimum_frequency = waveform_kwargs.pop('minimum_frequency') + catch_waveform_errors = waveform_kwargs.pop('catch_waveform_errors') + waveform_dictionary = set_waveform_dictionary(waveform_kwargs) + if len(waveform_kwargs) > 0: + raise ValueError(UNUSED_KWARGS_MESSAGE.format(waveform_kwargs=waveform_kwargs)) + + delta_time = float(time_array[1] - time_array[0]) + luminosity_distance = luminosity_distance * 1e6 * utils.parsec + mass_1 = mass_1 * utils.solar_mass + mass_2 = mass_2 * utils.solar_mass + iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = bilby_to_lalsimulation_spins( + theta_jn=theta_jn, phi_jl=phi_jl, tilt_1=tilt_1, tilt_2=tilt_2, + phi_12=phi_12, a_1=a_1, a_2=a_2, mass_1=mass_1, mass_2=mass_2, + reference_frequency=reference_frequency, phase=phase) + spin_arguments(waveform_approximant, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z) + total_mass = (mass_1 / lal.MSUN_SI + mass_2 / lal.MSUN_SI) * lal.MTSUN_SI # M in s, as lalsimulation computes it + + def evaluate(times): + try: + return function( + np.asarray(times) / total_mass, waveform_approximant, mass_1, mass_2, + spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, + luminosity_distance, iota, phase, delta_time, minimum_frequency, + reference_frequency, waveform_dictionary, truncate=True) + except Exception as e: + if not catch_waveform_errors or e.args[0] != 'Internal function call failed: Input domain error': + raise + logger.warning("Evaluating the waveform failed with error: {}\n".format(e) + + "Likelihood will be set to -inf.") + return None + + return evaluate + + def binary_black_hole_roq( frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **waveform_arguments): diff --git a/bilby/gw/time_domain_utils.py b/bilby/gw/time_domain_utils.py new file mode 100644 index 000000000..4087f4ff3 --- /dev/null +++ b/bilby/gw/time_domain_utils.py @@ -0,0 +1,163 @@ +""" +Utilities for time-domain gravitational-wave inference: products with the +inverse of a Toeplitz noise covariance matrix, and placement of a time-domain +waveform on the data samples. +""" +import numpy as np +from scipy.fft import next_fast_len +from scipy.linalg import solve_toeplitz + + +def gohberg_semencul_vectors(acf): + """ The Gohberg-Semencul vectors of the inverse of the symmetric Toeplitz matrix + with first row `acf` (Gohberg & Semencul 1972), computed with the Levinson + recursion + + Parameters + ========== + acf: array_like + The first row of the covariance matrix + + Returns + ======= + x, y: array_like + The first column of the inverse matrix, and the same reversed + """ + acf = np.asarray(acf, dtype=float) + unit_vector = np.zeros(len(acf)) + unit_vector[0] = 1.0 + x = solve_toeplitz((acf, acf), unit_vector) + return x, x[::-1] + + +def gohberg_semencul_product(x, y, vector): + """ The product of the inverse of a symmetric Toeplitz matrix with a vector, + computed from its Gohberg-Semencul vectors with FFTs, without forming the + inverse matrix (Gohberg & Semencul 1972) + + Parameters + ========== + x, y: array_like + The Gohberg-Semencul vectors, from :code:`gohberg_semencul_vectors` + vector: array_like + The vector, of the same length as `x` + + Returns + ======= + array_like: The product of the inverse matrix with `vector` + """ + x = np.asarray(x, dtype=float) + number_of_samples = len(x) + fft_length = next_fast_len(2 * number_of_samples - 1, real=True) + x_fft = np.fft.rfft(x, n=fft_length) + zy_fft = np.fft.rfft(np.concatenate(([0.0], np.asarray(y, dtype=float)[:-1])), n=fft_length) + + def lower(w_fft, vector_fft): + """ The product of the lower-triangular Toeplitz matrix with first column w and v, from their FFTs """ + return np.fft.irfft(w_fft * vector_fft, n=fft_length)[:number_of_samples] + + # L(w)^T v = reverse(L(w) reverse(v)) + reversed_fft = np.fft.rfft(np.asarray(vector, dtype=float)[::-1], n=fft_length) + first = lower(x_fft, np.fft.rfft(lower(x_fft, reversed_fft)[::-1], n=fft_length)) + second = lower(zy_fft, np.fft.rfft(lower(zy_fft, reversed_fft)[::-1], n=fft_length)) + return (first - second) / x[0] + + +def fractional_time_shift(signal, fraction, pad=32): + """ Delay a time series by a fraction of a sample with a Fourier phase shift + + Parameters + ========== + signal: array_like + The time series + fraction: float + The delay in samples (positive values delay the signal) + pad: int + The number of zeros added before and after the signal, so that nothing wraps around + + Returns + ======= + array_like: The delayed time series, of length len(signal) + 2 * pad + """ + signal = np.asarray(signal, dtype=float) + length = next_fast_len(len(signal) + 2 * pad) + buffer = np.zeros(length) + buffer[pad:pad + len(signal)] = signal + frequencies = np.fft.rfftfreq(length) + shifted = np.fft.irfft( + np.fft.rfft(buffer) * np.exp(-2j * np.pi * frequencies * fraction), n=length) + return shifted[:len(signal) + 2 * pad] + + +def _copy_into_segment(signal, offset, number_of_samples): + """ Put `signal` at data indices `offset`, `offset` + 1, ... of a segment of `number_of_samples` zeros """ + segment = np.zeros(number_of_samples) + first = max(0, -offset) + last = min(len(signal), number_of_samples - offset) + if last > first: + segment[offset + first:offset + last] = signal[first:last] + return segment + + +def times_from_merger(start_time, number_of_samples, sampling_frequency, geocent_time, time_delay): + """ The times of data samples relative to the arrival of the merger at a detector + + Parameters + ========== + start_time: float + The GPS time of the first data sample + number_of_samples: int + The number of data samples + sampling_frequency: float + The sampling frequency (in Hz) + geocent_time: float + The GPS time of the merger at the geocenter + time_delay: float + The time delay of the detector from the geocenter (in s) + + Returns + ======= + array_like: The times (in s) + """ + # subtract the two GPS times (~1e9 s) first, then add the small terms + return (start_time - geocent_time) - time_delay + np.arange(number_of_samples) / sampling_frequency + + +def place_time_domain_signal(signal, merger_index, arrival_index, number_of_samples, + placement="nearest"): + """ Place a time-domain signal, sampled at the sampling frequency of the data, + on the data samples + + Signal samples outside the data are dropped and data samples not covered by + the signal are zero. + + Parameters + ========== + signal: array_like + The signal on its own time grid + merger_index: float + The index of the merger in `signal` + arrival_index: float + The index in the data at which the merger arrives + number_of_samples: int + The number of data samples + placement: str + "nearest": shift by the whole number of samples nearest to the delay; + "fd_shift": shift by whole samples and the remaining fraction of a + sample with :code:`fractional_time_shift` + + Returns + ======= + array_like: The signal on the data samples + """ + position = arrival_index - merger_index + offset = int(np.floor(position + 0.5)) + if placement == "nearest": + return _copy_into_segment(signal, offset, number_of_samples) + elif placement == "fd_shift": + fraction = position - offset + pad = 32 + shifted = fractional_time_shift(signal, fraction, pad=pad) + return _copy_into_segment(shifted, offset - pad, number_of_samples) + else: + raise ValueError(f"placement must be 'nearest' or 'fd_shift', not {placement!r}") diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py index 1a55de975..e8ef14abe 100644 --- a/bilby/gw/utils.py +++ b/bilby/gw/utils.py @@ -684,6 +684,50 @@ def lalsim_SimInspiralChooseFDWaveform( return SimInspiralChooseFDWaveform(*args, waveform_dictionary, approximant) +def lalsim_SimInspiralChooseTDWaveform( + mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, + spin_2z, luminosity_distance, iota, phase, + longitude_ascending_nodes, eccentricity, mean_per_ano, delta_time, + minimum_frequency, reference_frequency, waveform_dictionary, + approximant): + """ + Safely call lalsimulation.SimInspiralChooseTDWaveform + + Parameters + ========== + mass_1: float, int + mass_2: float, int + spin_1x: float, int + spin_1y: float, int + spin_1z: float, int + spin_2x: float, int + spin_2y: float, int + spin_2z: float, int + luminosity_distance: float, int + iota: float, int + phase: float, int + longitude_ascending_nodes: float, int + eccentricity: float, int + mean_per_ano: float, int + delta_time: float, int + minimum_frequency: float, int + reference_frequency: float, int + waveform_dictionary: None, lal.Dict + approximant: int, str + """ + from lalsimulation import SimInspiralChooseTDWaveform + + args = convert_args_list_to_float( + mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, + luminosity_distance, iota, phase, longitude_ascending_nodes, + eccentricity, mean_per_ano, delta_time, minimum_frequency, + reference_frequency) + + approximant = _get_lalsim_approximant(approximant) + + return SimInspiralChooseTDWaveform(*args, waveform_dictionary, approximant) + + @lru_cache(maxsize=10) def _get_lalsim_approximant(approximant): if isinstance(approximant, int): diff --git a/examples/gw_examples/data_examples/GW150914_time_domain.py b/examples/gw_examples/data_examples/GW150914_time_domain.py new file mode 100644 index 000000000..8f918cc1f --- /dev/null +++ b/examples/gw_examples/data_examples/GW150914_time_domain.py @@ -0,0 +1,124 @@ +#!/usr/bin/env python +""" +Tutorial to demonstrate running time-domain parameter estimation on GW150914 + +This is the GW150914.py analysis done in the time domain. The strain around +the analysis segment is downsampled to twice the maximum frequency, the noise +model of each detector is an autocovariance function computed from the PSD of +the GWTC-2.1 parameter-estimation release [1], and the IMRPhenomTPHM waveform +is evaluated at the data samples (placement "model_time"). Time and phase +marginalization are not available for the time-domain likelihood. + +[1] https://gwosc.org/eventapi/html/GWTC-2.1-confident/GW150914/v4/ +""" +import bilby +import h5py +from gwpy.timeseries import TimeSeries + +logger = bilby.core.utils.logger +outdir = "outdir" +label = "GW150914_time_domain" + +# Note you can get trigger times using the gwosc package, e.g.: +# > from gwosc import datasets +# > datasets.event_gps("GW150914") +trigger_time = 1126259462.4 +detectors = ["H1", "L1"] +maximum_frequency = 512 +minimum_frequency = 20 +sampling_frequency = 2 * maximum_frequency # The data are downsampled to this rate +duration = 4 # Analysis segment duration +post_trigger_duration = 2 # Time between trigger time and end of segment +end_time = trigger_time + post_trigger_duration +start_time = end_time - duration + +# Data downloaded on each side of the analysis segment, removed by the +# downsampling (condition_strain_data, which needs duration / 2 on each side) +padding = duration / 2 + +# The PSDs of the GWTC-2.1 parameter-estimation release, from GWOSC [1]: +# https://zenodo.org/records/6513631/files/IGWN-GWTC2p1-v2-GW150914_095045_PEDataRelease_mixed_cosmo.h5 +psd_file = "IGWN-GWTC2p1-v2-GW150914_095045_PEDataRelease_mixed_cosmo.h5" +psd_label = "C01:IMRPhenomXPHM" +with h5py.File(psd_file, "r") as ff: + psds = {det: ff[psd_label]["psds"][det][:] for det in detectors} + +ifo_list = bilby.gw.detector.InterferometerList([]) +for det in detectors: + logger.info("Downloading analysis data for ifo {}".format(det)) + ifo = bilby.gw.detector.get_empty_interferometer(det) + ifo.set_analysis_window(start_time, duration) + data = TimeSeries.fetch_open_data(det, start_time - padding, end_time + padding) + ifo.strain_data.set_from_gwpy_timeseries(data) + ifo.condition_strain_data(sampling_frequency) + + ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity( + frequency_array=psds[det][:, 0], psd_array=psds[det][:, 1] + ) + ifo.maximum_frequency = maximum_frequency + ifo.minimum_frequency = minimum_frequency + ifo.set_autocovariance_function_from_power_spectral_density() + ifo_list.append(ifo) + +logger.info("Saving data plots to {}".format(outdir)) +bilby.core.utils.check_directory_exists_and_if_not_mkdir(outdir) +ifo_list.plot_data(outdir=outdir, label=label) + +priors = bilby.gw.prior.BBHPriorDict(filename="GW150914.prior") + +priors["geocent_time"] = bilby.core.prior.Uniform( + trigger_time - 0.1, trigger_time + 0.1, name="geocent_time" +) + +""" +Fixed arguments passed into the source model. + +There are two different minimum frequencies: +- ifo.minimum_frequency (20 Hz): the lowest frequency used in the analysis. +- minimum_frequency below (10 Hz): the frequency at which the waveform starts. + +The waveform starts at 10 Hz, below 20 Hz, for two reasons: +1. Higher modes. Higher modes oscillate faster than the main (2, 2) mode + (the (4, 4) mode at twice its frequency). If the (2, 2) mode starts at + 10 Hz, all modes are already present at 20 Hz. +2. The start of the waveform must not be in the data. The waveform switches + on suddenly, and this sudden start disturbs the analysis even above 20 Hz. + Starting at 10 Hz, every waveform in the prior begins before the analysis + window, so its start is never in the data. With a wider prior or a longer + window, use a lower value. +""" + +waveform_generator = bilby.gw.WaveformGenerator( + time_domain_source_model=bilby.gw.source.lal_binary_black_hole_model_time, + parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters, + waveform_arguments={ + "waveform_approximant": "IMRPhenomTPHM", + "reference_frequency": 20, + "minimum_frequency": 10, + }, +) + +likelihood = bilby.gw.likelihood.TimeDomainGravitationalWaveTransient( + ifo_list, + waveform_generator, + priors=priors, + distance_marginalization=False, + placement="model_time", +) + +result = bilby.run_sampler( + likelihood, + priors, + sampler="dynesty", + sample="acceptance-walk", + nlive=500, + naccept=60, + outdir=outdir, + label=label, + check_point_delta_t=600, + check_point_plot=True, + npool=1, + conversion_function=bilby.gw.conversion.generate_all_bbh_parameters, + result_class=bilby.gw.result.CBCResult, +) +result.plot_corner() diff --git a/examples/gw_examples/injection_examples/fast_tutorial_time_domain.py b/examples/gw_examples/injection_examples/fast_tutorial_time_domain.py new file mode 100644 index 000000000..db86ebb12 --- /dev/null +++ b/examples/gw_examples/injection_examples/fast_tutorial_time_domain.py @@ -0,0 +1,137 @@ +#!/usr/bin/env python +""" +Tutorial to demonstrate running time-domain parameter estimation on a reduced +parameter space for an injected precessing signal. +""" + +import bilby + +# Set the duration and sampling frequency of the data segment that we're +# going to inject the signal into. The analysis window is the whole segment. +duration = 4.0 +sampling_frequency = 2048.0 +minimum_frequency = 20 + +# Specify the output directory and the name of the simulation. +outdir = "outdir" +label = "fast_tutorial_time_domain" +bilby.core.utils.setup_logger(outdir=outdir, label=label, log_level="DEBUG") + +# Set up a random seed for result reproducibility. This is optional! +bilby.core.utils.random.seed(88170235) + +injection_parameters = dict( + mass_1=36.0, + mass_2=29.0, + a_1=0.4, + a_2=0.3, + tilt_1=0.5, + tilt_2=1.0, + phi_12=1.7, + phi_jl=0.3, + luminosity_distance=2000.0, + theta_jn=0.4, + psi=2.659, + phase=1.3, + geocent_time=1126259642.413, + ra=1.375, + dec=-1.2108, +) +start_time = injection_parameters["geocent_time"] - 2 + +""" +Fixed arguments passed into the source model. + +There are two different minimum frequencies: +- ifo.minimum_frequency (20 Hz): the lowest frequency used in the analysis. +- minimum_frequency below (10 Hz): the frequency at which the waveform starts. + +The waveform starts at 10 Hz, below 20 Hz, for two reasons: +1. Higher modes. Higher modes oscillate faster than the main (2, 2) mode + (the (4, 4) mode at twice its frequency). If the (2, 2) mode starts at + 10 Hz, all modes are already present at 20 Hz. +2. The start of the waveform must not be in the data. The waveform switches + on suddenly, and this sudden start disturbs the analysis even above 20 Hz. + Starting at 10 Hz, every waveform in the prior begins before the analysis + window, so its start is never in the data. With a wider prior or a longer + window, use a lower value. +""" + +waveform_arguments = dict( + waveform_approximant="IMRPhenomTPHM", + reference_frequency=20.0, + minimum_frequency=10.0, +) + +# Create the waveform_generator using the model-time LAL BinaryBlackHole source +# function, at the sampling frequency of the data +waveform_generator = bilby.gw.WaveformGenerator( + duration=duration, + sampling_frequency=sampling_frequency, + start_time=start_time, + time_domain_source_model=bilby.gw.source.lal_binary_black_hole_model_time, + parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters, + waveform_arguments=waveform_arguments, +) + +# Set up interferometers. In this case we'll use two interferometers +# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design +# sensitivity. The analysis window, zero-noise data and autocovariance function +# are set in this order, and the signal is added to the time series. +ifos = bilby.gw.detector.InterferometerList(["H1", "L1"]) +for ifo in ifos: + ifo.minimum_frequency = minimum_frequency +ifos.set_analysis_windows(duration=duration, start_time=start_time) +ifos.set_strain_data_from_zero_noise( + sampling_frequency=sampling_frequency, + duration=duration, + start_time=start_time, +) +ifos.set_autocovariance_functions_from_power_spectral_densities() +ifos.inject_signal_time_domain( + waveform_generator=waveform_generator, + parameters=injection_parameters, + placement="model_time", +) + +priors = bilby.gw.prior.BBHPriorDict() +for key in [ + "a_1", + "a_2", + "tilt_1", + "tilt_2", + "phi_12", + "phi_jl", + "psi", + "ra", + "dec", + "geocent_time", + "phase", +]: + priors[key] = injection_parameters[key] + +priors["chirp_mass"] = bilby.gw.prior.UniformInComponentsChirpMass( + name="chirp_mass", minimum=25, maximum=30, unit="$M_{\\odot}$" +) + +# Perform a check that the prior does not extend to a parameter space longer than the data +priors.validate_prior(duration, minimum_frequency) + +likelihood = bilby.gw.likelihood.TimeDomainGravitationalWaveTransient( + interferometers=ifos, waveform_generator=waveform_generator, placement="model_time" +) + +result = bilby.run_sampler( + likelihood=likelihood, + priors=priors, + sampler="dynesty", + sample="acceptance-walk", + nlive=500, + naccept=60, + injection_parameters=injection_parameters, + outdir=outdir, + label=label, + result_class=bilby.gw.result.CBCResult, +) + +result.plot_corner() diff --git a/examples/gw_examples/injection_examples/time_domain_heterodyning.py b/examples/gw_examples/injection_examples/time_domain_heterodyning.py new file mode 100644 index 000000000..c5f25d4d7 --- /dev/null +++ b/examples/gw_examples/injection_examples/time_domain_heterodyning.py @@ -0,0 +1,159 @@ +#!/usr/bin/env python +""" +Tutorial to demonstrate running time-domain parameter estimation on a reduced +parameter space for an injected precessing signal, using the time-domain +heterodyned likelihood. + +This is the heterodyning analysis done in the time domain, with zero +noise and the IMRPhenomTPHM waveform. The heterodyned likelihood evaluates the +waveform only at the edges of time bins; the posterior samples are then +compared with the exact time-domain likelihood (placement "model_time"). +""" + +import bilby +from bilby.core.utils import random + +random.seed(123) + +# Set the duration and sampling frequency of the data segment that we're +# going to inject the signal into. The analysis window is the whole segment. +duration = 4.0 +sampling_frequency = 2048.0 +minimum_frequency = 20 + +# Specify the output directory and the name of the simulation. +outdir = "outdir" +label = "time_domain_heterodyning" +bilby.core.utils.setup_logger(outdir=outdir, label=label) + +injection_parameters = dict( + mass_1=36.0, + mass_2=29.0, + a_1=0.4, + a_2=0.3, + tilt_1=0.5, + tilt_2=1.0, + phi_12=1.7, + phi_jl=0.3, + luminosity_distance=2000.0, + theta_jn=0.4, + psi=2.659, + phase=1.3, + geocent_time=1126259642.413, + ra=1.375, + dec=-1.2108, +) +start_time = injection_parameters["geocent_time"] - 2 + +""" +Fixed arguments passed into the source model. + +There are two different minimum frequencies: +- ifo.minimum_frequency (20 Hz): the lowest frequency used in the analysis. +- minimum_frequency below (10 Hz): the frequency at which the waveform starts. + +The waveform starts at 10 Hz, below 20 Hz, for two reasons: +1. Higher modes. Higher modes oscillate faster than the main (2, 2) mode + (the (4, 4) mode at twice its frequency). If the (2, 2) mode starts at + 10 Hz, all modes are already present at 20 Hz. +2. The start of the waveform must not be in the data. The waveform switches + on suddenly, and this sudden start disturbs the analysis even above 20 Hz. + Starting at 10 Hz, every waveform in the prior begins before the analysis + window, so its start is never in the data. With a wider prior or a longer + window, use a lower value. +""" + +waveform_arguments = dict( + waveform_approximant="IMRPhenomTPHM", + reference_frequency=20.0, + minimum_frequency=10.0, + fiducial=1, +) + +# Create the waveform_generator using the heterodyning LAL BinaryBlackHole +# source function, at the sampling frequency of the data +waveform_generator = bilby.gw.WaveformGenerator( + duration=duration, + sampling_frequency=sampling_frequency, + start_time=start_time, + time_domain_source_model=bilby.gw.source.lal_binary_black_hole_time_domain_heterodyning, + parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters, + waveform_arguments=waveform_arguments, +) + +# Set up interferometers. In this case we'll use two interferometers +# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design +# sensitivity. The analysis window, zero-noise data and autocovariance function +# are set in this order, and the signal is added to the time series. +ifos = bilby.gw.detector.InterferometerList(["H1", "L1"]) +for ifo in ifos: + ifo.minimum_frequency = minimum_frequency +ifos.set_analysis_windows(duration=duration, start_time=start_time) +ifos.set_strain_data_from_zero_noise( + sampling_frequency=sampling_frequency, + duration=duration, + start_time=start_time, +) +ifos.set_autocovariance_functions_from_power_spectral_densities() +ifos.inject_signal_time_domain( + waveform_generator=waveform_generator, + parameters=injection_parameters, + placement="model_time", +) + +priors = bilby.gw.prior.BBHPriorDict() +for key in [ + "a_1", + "a_2", + "tilt_1", + "tilt_2", + "phi_12", + "phi_jl", + "psi", + "ra", + "dec", + "geocent_time", + "phase", +]: + priors[key] = injection_parameters[key] + +priors["chirp_mass"] = bilby.gw.prior.UniformInComponentsChirpMass( + name="chirp_mass", minimum=25, maximum=30, unit="$M_{\\odot}$" +) + +# Perform a check that the prior does not extend to a parameter space longer than the data +priors.validate_prior(duration, minimum_frequency) + +fiducial_parameters = injection_parameters.copy() +m1 = fiducial_parameters.pop("mass_1") +m2 = fiducial_parameters.pop("mass_2") +fiducial_parameters["chirp_mass"] = bilby.gw.conversion.component_masses_to_chirp_mass( + m1, m2 +) +fiducial_parameters["mass_ratio"] = m2 / m1 + +likelihood = bilby.gw.likelihood.TimeDomainHeterodyningGravitationalWaveTransient( + interferometers=ifos, + waveform_generator=waveform_generator, + priors=priors, + distance_marginalization=False, + fiducial_parameters=fiducial_parameters, +) + +result = bilby.run_sampler( + likelihood=likelihood, + priors=priors, + sampler="dynesty", + sample="acceptance-walk", + nlive=500, + naccept=60, + injection_parameters=injection_parameters, + outdir=outdir, + label=label, + conversion_function=bilby.gw.conversion.generate_all_bbh_parameters, + result_class=bilby.gw.result.CBCResult, +) + +result.plot_corner() + +likelihood.compare_with_exact_likelihood(result.posterior, outdir=outdir, label=label) diff --git a/examples/gw_examples/injection_examples/time_domain_inspiral.py b/examples/gw_examples/injection_examples/time_domain_inspiral.py new file mode 100644 index 000000000..7cf3e8290 --- /dev/null +++ b/examples/gw_examples/injection_examples/time_domain_inspiral.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python +""" +Tutorial to demonstrate running time-domain parameter estimation on the +inspiral of an injected precessing signal. + +The analysis window ends at segment_cut_time relative to the merger, in each +detector (the arrival of the cut at that detector). The data are zero noise and +the IMRPhenomTPHM waveform is evaluated at the data samples (placement +"model_time"). See time_domain_post_inspiral.py for the post-inspiral analysis. +""" + +import bilby +import lal + +# The data segment covers the analysis window (the 2 s before the cut, in each +# detector), from 3 s before the merger to 1 s after it. +duration = 4.0 +sampling_frequency = 2048.0 +minimum_frequency = 20 + +# The analysis window is the 2 s before the cut, in each detector. +inspiral_analysis_duration = 2.0 + +# Specify the output directory and the name of the simulation. +outdir = "outdir" +label = "time_domain_inspiral" +bilby.core.utils.setup_logger(outdir=outdir, label=label) + +bilby.core.utils.random.seed(88170235) + +injection_parameters = dict( + mass_1=36.0, + mass_2=29.0, + a_1=0.4, + a_2=0.3, + tilt_1=0.5, + tilt_2=1.0, + phi_12=1.7, + phi_jl=0.3, + luminosity_distance=2000.0, + theta_jn=0.4, + psi=2.659, + phase=1.3, + geocent_time=1126259642.413, + ra=1.375, + dec=-1.2108, +) +start_time = injection_parameters["geocent_time"] - 3 + +# The cut is 20 M before the merger (geocent_time, where the (2, 2) mode of +# the waveform peaks), with M the total mass of the injection in seconds. +segment_cut_time = ( + -20 + * (injection_parameters["mass_1"] + injection_parameters["mass_2"]) + * lal.MTSUN_SI +) + +""" +Fixed arguments passed into the source model. + +There are two different minimum frequencies: +- ifo.minimum_frequency (20 Hz): the lowest frequency used in the analysis. +- minimum_frequency below (10 Hz): the frequency at which the waveform starts. + +The waveform starts at 10 Hz, below 20 Hz, for two reasons: +1. Higher modes. Higher modes oscillate faster than the main (2, 2) mode + (the (4, 4) mode at twice its frequency). If the (2, 2) mode starts at + 10 Hz, all modes are already present at 20 Hz. +2. The start of the waveform must not be in the data. The waveform switches + on suddenly, and this sudden start disturbs the analysis even above 20 Hz. + Starting at 10 Hz, every waveform in the prior begins before the analysis + window, so its start is never in the data. With a wider prior or a longer + window, use a lower value. +""" + +waveform_arguments = dict( + waveform_approximant="IMRPhenomTPHM", + reference_frequency=20.0, + minimum_frequency=10.0, +) + +# Create the waveform_generator using the model-time LAL BinaryBlackHole source +# function, at the sampling frequency of the data +waveform_generator = bilby.gw.WaveformGenerator( + duration=duration, + sampling_frequency=sampling_frequency, + start_time=start_time, + time_domain_source_model=bilby.gw.source.lal_binary_black_hole_model_time, + parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters, + waveform_arguments=waveform_arguments, +) + +priors = bilby.gw.prior.BBHPriorDict() +for key in [ + "a_1", + "a_2", + "tilt_1", + "tilt_2", + "phi_12", + "phi_jl", + "psi", + "ra", + "dec", + "geocent_time", + "phase", +]: + priors[key] = injection_parameters[key] + +priors["chirp_mass"] = bilby.gw.prior.UniformInComponentsChirpMass( + name="chirp_mass", minimum=25, maximum=30, unit="$M_{\\odot}$" +) + +# Perform a check that the prior does not extend to a parameter space longer than the data +priors.validate_prior(duration, minimum_frequency) + +# Set up interferometers. In this case we'll use two interferometers +# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design +# sensitivity. The analysis window (cut at the arrival of segment_cut_time in +# each detector), zero-noise data and autocovariance function are set in this +# order, and the signal is added to the time series. +ifos = bilby.gw.detector.InterferometerList(["H1", "L1"]) +for ifo in ifos: + ifo.minimum_frequency = minimum_frequency +ifos.set_analysis_windows( + duration=inspiral_analysis_duration, + analysis_segment="inspiral", + segment_cut_time=segment_cut_time, + cut_reference_parameters=injection_parameters, +) +ifos.set_strain_data_from_zero_noise( + sampling_frequency=sampling_frequency, + duration=duration, + start_time=start_time, +) +ifos.set_autocovariance_functions_from_power_spectral_densities() +ifos.inject_signal_time_domain( + waveform_generator=waveform_generator, + parameters=injection_parameters, + placement="model_time", +) + +likelihood = bilby.gw.likelihood.TimeDomainGravitationalWaveTransient( + interferometers=ifos, waveform_generator=waveform_generator, placement="model_time" +) +result = bilby.run_sampler( + likelihood=likelihood, + priors=priors, + sampler="dynesty", + sample="acceptance-walk", + nlive=500, + naccept=60, + injection_parameters=injection_parameters, + outdir=outdir, + label=label, + result_class=bilby.gw.result.CBCResult, +) + +result.plot_corner() diff --git a/examples/gw_examples/injection_examples/time_domain_post_inspiral.py b/examples/gw_examples/injection_examples/time_domain_post_inspiral.py new file mode 100644 index 000000000..9faa86e74 --- /dev/null +++ b/examples/gw_examples/injection_examples/time_domain_post_inspiral.py @@ -0,0 +1,160 @@ +#!/usr/bin/env python +""" +Tutorial to demonstrate running time-domain parameter estimation on the +post-inspiral of an injected precessing signal. + +The analysis window starts at segment_cut_time relative to the merger, in each +detector (the arrival of the cut at that detector). The data are zero noise and +the IMRPhenomTPHM waveform is evaluated at the data samples (placement +"model_time"). See time_domain_inspiral.py for the inspiral analysis. + +""" + +import bilby +import lal + +# The data segment covers the analysis window (the 0.5 s after the cut, in each +# detector), from 2 s before the merger to 2 s after it. +duration = 4.0 +sampling_frequency = 2048.0 +minimum_frequency = 20 + +# The analysis window is the 0.5 s after the cut, in each detector. +postinspiral_analysis_duration = 0.5 + +# Specify the output directory and the name of the simulation. +outdir = "outdir" +label = "time_domain_post_inspiral" +bilby.core.utils.setup_logger(outdir=outdir, label=label) + +bilby.core.utils.random.seed(88170235) + +injection_parameters = dict( + mass_1=36.0, + mass_2=29.0, + a_1=0.4, + a_2=0.3, + tilt_1=0.5, + tilt_2=1.0, + phi_12=1.7, + phi_jl=0.3, + luminosity_distance=2000.0, + theta_jn=0.4, + psi=2.659, + phase=1.3, + geocent_time=1126259642.413, + ra=1.375, + dec=-1.2108, +) +start_time = injection_parameters["geocent_time"] - 2 + +# The cut is 20 M before the merger (geocent_time, where the (2, 2) mode of +# the waveform peaks), with M the total mass of the injection in seconds. +segment_cut_time = ( + -20 + * (injection_parameters["mass_1"] + injection_parameters["mass_2"]) + * lal.MTSUN_SI +) + +""" +Fixed arguments passed into the source model. + +There are two different minimum frequencies: +- ifo.minimum_frequency (20 Hz): the lowest frequency used in the analysis. +- minimum_frequency below (10 Hz): the frequency at which the waveform starts. + +The waveform starts at 10 Hz, below 20 Hz, for two reasons: +1. Higher modes. Higher modes oscillate faster than the main (2, 2) mode + (the (4, 4) mode at twice its frequency). If the (2, 2) mode starts at + 10 Hz, all modes are already present at 20 Hz. +2. The start of the waveform must not be in the data. The waveform switches + on suddenly, and this sudden start disturbs the analysis even above 20 Hz. + Starting at 10 Hz, every waveform in the prior begins before the analysis + window, so its start is never in the data. With a wider prior or a longer + window, use a lower value. +""" + +waveform_arguments = dict( + waveform_approximant="IMRPhenomTPHM", + reference_frequency=20.0, + minimum_frequency=10.0, +) + +# Create the waveform_generator using the model-time LAL BinaryBlackHole source +# function, at the sampling frequency of the data +waveform_generator = bilby.gw.WaveformGenerator( + duration=duration, + sampling_frequency=sampling_frequency, + start_time=start_time, + time_domain_source_model=bilby.gw.source.lal_binary_black_hole_model_time, + parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters, + waveform_arguments=waveform_arguments, +) + +priors = bilby.gw.prior.BBHPriorDict() +for key in [ + "a_1", + "a_2", + "tilt_1", + "tilt_2", + "phi_12", + "phi_jl", + "psi", + "ra", + "dec", + "geocent_time", + "phase", +]: + priors[key] = injection_parameters[key] + +priors["chirp_mass"] = bilby.gw.prior.UniformInComponentsChirpMass( + name="chirp_mass", minimum=25, maximum=30, unit="$M_{\\odot}$" +) + +# Perform a check that the prior does not extend to a parameter space longer than the data +priors.validate_prior(duration, minimum_frequency) + +# Set up interferometers. In this case we'll use two interferometers +# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design +# sensitivity. The analysis window (cut at the arrival of segment_cut_time in +# each detector), zero-noise data and autocovariance function are set in this +# order, and the signal is added to the time series. +ifos = bilby.gw.detector.InterferometerList(["H1", "L1"]) +for ifo in ifos: + ifo.minimum_frequency = minimum_frequency +ifos.set_analysis_windows( + duration=postinspiral_analysis_duration, + analysis_segment="post_inspiral", + segment_cut_time=segment_cut_time, + cut_reference_parameters=injection_parameters, +) +ifos.set_strain_data_from_zero_noise( + sampling_frequency=sampling_frequency, + duration=duration, + start_time=start_time, +) +ifos.set_autocovariance_functions_from_power_spectral_densities() +ifos.inject_signal_time_domain( + waveform_generator=waveform_generator, + parameters=injection_parameters, + placement="model_time", +) + +likelihood = bilby.gw.likelihood.TimeDomainGravitationalWaveTransient( + interferometers=ifos, waveform_generator=waveform_generator, placement="model_time" +) + +result = bilby.run_sampler( + likelihood=likelihood, + priors=priors, + sampler="dynesty", + sample="acceptance-walk", + nlive=500, + naccept=60, + injection_parameters=injection_parameters, + outdir=outdir, + label=label, + result_class=bilby.gw.result.CBCResult, +) + +result.plot_corner()