Generating waveforms

Training data for Dingo consist of pairs of parameters \(\theta\) and corresponding simulated strain data sets \(d_I\), where \(I\) runs over the GW interferometers (L1, H1, V1, etc.). Additionally, when conditioning on detector noise properties, data also include noise context (the PSD \(S_{\text{n},I}\)). Strain data sets are of the form

\[ d_I = h_I(\theta) + n_I, \]

where \(h_I(\theta)\) is a signal waveform (provided by a waveform model) and \(n_I\) is a noise realization (stationary and Gaussian, consistent with \(S_{\text{n}, I}\)).

Data domain

At present, Dingo works entirely with frequency domain data. Although NPE is very flexible and could in principle learn to interpret data in any representation, FD data are especially convenient because (1) stationary Gaussian noise is independent in each frequency bin, so noise generation is straightforward, (2) time shifts take a simple form, enabling improved data augmentation, and (3) the noise context is already in FD. Other domains could be useful in the future, however, so the code is written in a way that the domain could be adapted.

The domain is specified by instantiating a UniformFrequencyDomain,

import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import lal

from dingo.gw.domains import UniformFrequencyDomain
domain = UniformFrequencyDomain(f_min=20.0, f_max=1024.0, delta_f=0.125)

Derived class properties include, e.g., the frequency grid. Frequency arrays run from 0 to f_max, as is standard for GW data analysis software.

domain.sample_frequencies

array([0.000000e+00, 1.250000e-01, 2.500000e-01, ..., 1.023750e+03,
       1.023875e+03, 1.024000e+03], shape=(8193,), dtype=float32)

Note

The standard deviation of white noise in each frequency bin is stored in domain.noise_std. In frequency domain, this is given by \(\sqrt{1/4\delta f}\). See data conditioning for details.

Various class methods also act on data, to perform operations such as zeroing below f_min, truncating above f_max, or applying a time shift:

class dingo.gw.domains.UniformFrequencyDomain(f_min: float, f_max: float, delta_f: float)

Defines the physical domain on which the data of interest live.

The frequency bins are assumed to be uniform between [0, f_max] with spacing delta_f.

Given a finite length of time domain data, the Fourier domain data starts at a frequency f_min and is zero below this frequency.

window_kwargs specify windowing used for FFT to obtain FD data from TD data in practice.

Parameters:

• f_min (float)

• f_max (float)

• delta_f (float)

property delta_f: float

The frequency spacing of the uniform grid [Hz].

property domain_dict

Enables to rebuild the domain via calling build_domain(domain_dict).

property duration: float

Waveform duration in seconds.

property f_max: float

The maximum frequency [Hz] is typically set to half the sampling rate.

property f_min: float

The minimum frequency [Hz].

property frequency_mask: ndarray

Mask which selects frequency bins greater than or equal to the starting frequency

property frequency_mask_length: int

Number of samples in the subdomain domain[frequency_mask].

get_sample_frequencies_astype(data: ndarray | Tensor) → ndarray | Tensor

Returns a 1D frequency array compatible with the last index of data array.

Decides whether array is numpy or torch tensor (and cuda vs cpu), and whether it contains the leading zeros below f_min.

Parameters:

data (Union[np.array, torch.Tensor]) – Sample data

Return type:

frequency array compatible with last index

property sampling_rate: float

The sampling rate of the data [Hz].

update(new_settings: dict)

Update the domain with new settings. This is only allowed if the new settings are “compatible” with the old ones. E.g., f_min should be larger than the existing f_min.

Parameters:

new_settings (dict) – Settings dictionary. Must contain a subset of the keys contained in domain_dict.

update_data(data: ndarray, axis: int = -1, low_value: float = 0.0)

Adjusts data to be compatible with the domain:

• Below f_min, it sets the data to low_value (typically 0.0 for a waveform, but for a PSD this might be a large value).

• Above f_max, it truncates the data array.

Parameters:

• data (np.ndarray) – Data array

• axis (int) – Which data axis to apply the adjustment along.

• low_value (float) – Below f_min, set the data to this value.

Returns:

The new data array.

Return type:

np.ndarray

Waveform generator

Waveforms are generated using the WaveformGenerator class (or its subclass NewInterfaceWaveformGenerator, for employing the new LIGO waveform interface, needed for some approximants). This depends on a Domain as well as a waveform approximant and a reference frequency f_ref. In the backend, the WaveformGenerator class calls LALSimulation functions (typically SimInspiralFD) via the SWIG-Python interface. For time domain waveforms, SimInspiralFD takes care of FFTing to frequency domain. The NewInterfaceWaveformGenerator class calls the gwsignal module, a Python interface recently implemented in LALSimulation, which is needed for employing some of the latest waveform approximants, as the SEOBNRv5HM and SEOBNRv5PHM.

from dingo.gw.waveform_generator import WaveformGenerator #, NewInterfaceWaveformGenerator

wfg = WaveformGenerator(approximant='IMRPhenomXPHM', domain=domain, f_ref=20.0)
# wfg = NewInterfaceWaveformGenerator(approximant='SEOBNRv5PHM', domain=domain, f_ref=20.0)

Setting spin_conversion_phase = None. Using phase parameter for conversion to cartesian spins.

To generate a waveform we first need to choose parameters. Here we sample parameters from a bilby.core.prior.PriorDict. We use the default Dingo intrinsic prior.

from bilby.core.prior import PriorDict
from dingo.gw.prior import default_intrinsic_dict

prior = PriorDict(default_intrinsic_dict)
prior

{'mass_1': Constraint(minimum=10.0, maximum=80.0, name='mass_1', latex_label='$m_1$', unit=None),
 'mass_2': Constraint(minimum=10.0, maximum=80.0, name='mass_2', latex_label='$m_2$', unit=None),
 'mass_ratio': bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0, name='mass_ratio', latex_label='$q$', unit=None, boundary=None, equal_mass=False),
 'chirp_mass': bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0, name='chirp_mass', latex_label='$\\mathcal{M}$', unit=None, boundary=None),
 'luminosity_distance': DeltaFunction(peak=1000.0, name=None, latex_label=None, unit=None),
 'theta_jn': Sine(minimum=0.0, maximum=3.141592653589793, name='theta_jn', latex_label='$\\theta_{JN}$', unit=None, boundary=None),
 'phase': Uniform(minimum=0.0, maximum=6.283185307179586, name='phase', latex_label='$\\phi$', unit=None, boundary='periodic'),
 'a_1': Uniform(minimum=0.0, maximum=0.99, name='a_1', latex_label='$a_1$', unit=None, boundary=None),
 'a_2': Uniform(minimum=0.0, maximum=0.99, name='a_2', latex_label='$a_2$', unit=None, boundary=None),
 'tilt_1': Sine(minimum=0.0, maximum=3.141592653589793, name='tilt_1', latex_label='$\\theta_1$', unit=None, boundary=None),
 'tilt_2': Sine(minimum=0.0, maximum=3.141592653589793, name='tilt_2', latex_label='$\\theta_2$', unit=None, boundary=None),
 'phi_12': Uniform(minimum=0.0, maximum=6.283185307179586, name='phi_12', latex_label='$\\Delta\\phi$', unit=None, boundary='periodic'),
 'phi_jl': Uniform(minimum=0.0, maximum=6.283185307179586, name='phi_jl', latex_label='$\\phi_{JL}$', unit=None, boundary='periodic'),
 'geocent_time': DeltaFunction(peak=0.0, name=None, latex_label=None, unit=None)}

p = prior.sample()
p

{'mass_ratio': 0.7494256666503227,
 'chirp_mass': 98.08572768991085,
 'luminosity_distance': 1000.0,
 'theta_jn': np.float64(1.3287182465211944),
 'phase': 0.3888401594037996,
 'a_1': 0.14340135010437627,
 'a_2': 0.6549630187410872,
 'tilt_1': np.float64(2.3718796322236075),
 'tilt_2': np.float64(1.6161703419232794),
 'phi_12': 0.9218708216349888,
 'phi_jl': 2.147703087420897,
 'geocent_time': 0.0}

Finally, we generate the waveform. This is returned as a dictionary, with entries for each polarization. This way of representing a sample is used throughout Dingo, and will be very convenient when applying transforms (to apply extrinsic parameters, add noise, etc.).

h = wfg.generate_hplus_hcross(p)
h

{'h_plus': array([0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j], shape=(8193,)),
 'h_cross': array([0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j], shape=(8193,))}

import matplotlib.pyplot as plt
plt.plot(domain.sample_frequencies, h['h_plus'].real, label='real')
plt.plot(domain.sample_frequencies, h['h_plus'].imag, label='imag')
plt.xlim((10,1024))
plt.xscale('log')
plt.legend()
plt.xlabel('f')
plt.ylabel(r'$h_+$')
plt.show()

Note that the waveform is nonzero slightly below f_min. This simply arises from the model implementation in LALSimulation. When training networks, input data will be truncated below f_min.

The complete specification of the WaveformGenerator class is given as

class dingo.gw.waveform_generator.WaveformGenerator(approximant: str, domain: Domain, f_ref: float, f_start: float | None = None, mode_list: List[Tuple] | None = None, transform=None, spin_conversion_phase=None, **kwargs)

Generate polarizations using LALSimulation routines in the specified domain for a single GW coalescence given a set of waveform parameters.

Parameters:

• approximant (str) – Waveform “approximant” string understood by lalsimulation This is defines which waveform model is used.

• domain (Domain) – Domain object that specifies on which physical domain the waveform polarizations will be generated, e.g. Fourier domain, time domain.

• f_ref (float) – Reference frequency for the waveforms

• f_start (float) – Starting frequency for waveform generation. This is optional, and if not included, the starting frequency will be set to f_min. This exists so that EOB waveforms can be generated starting from a lower frequency than f_min.

• mode_list (List[Tuple]) – A list of waveform (ell, m) modes to include when generating the polarizations.

• spin_conversion_phase (float = None) – Value for phiRef when computing cartesian spins from bilby spins via bilby_to_lalsimulation_spins. The common convention is to use the value of the phase parameter here, which is also used in the spherical harmonics when combining the different modes. If spin_conversion_phase = None, this default behavior is adapted. For dingo, this convention for the phase parameter makes it impossible to treat the phase as an extrinsic parameter, since we can only account for the change of phase in the spherical harmonics when changing the phase (in order to also change the cartesian spins – specifically, to rotate the spins by phase in the sx-sy plane – one would need to recompute the modes, which is expensive). By setting spin_conversion_phase != None, we impose the convention to always use phase = spin_conversion_phase when computing the cartesian spins.

generate_FD_modes_LO(parameters)

Generate FD modes in the L0 frame.

Parameters:

parameters (dict) – Dictionary of parameters for the waveform. For details see see self.generate_hplus_hcross.

Returns:

• hlm_fd (dict) – Dictionary with (l,m) as keys and the corresponding FD modes in lal format as values.

• iota (float)

generate_FD_waveform(parameters_lal: Tuple, target_function: Callable) → Dict[str, ndarray]

Generate Fourier domain GW polarizations (h_plus, h_cross).

Parameters:

• parameters_lal – A tuple of parameters for the lalsimulation waveform generator

• target_function – Lalsimulation function for waveform generation.

Returns:

A dictionary of generated waveform polarizations

Return type:

pol_dict

generate_TD_modes_L0(parameters)

Generate TD modes in the L0 frame.

Parameters:

parameters (dict) – Dictionary of parameters for the waveform. For details see see self.generate_hplus_hcross.

Returns:

• hlm_td (dict) – Dictionary with (l,m) as keys and the corresponding TD modes in lal format as values.

• iota (float)

generate_TD_waveform(parameters_lal: Tuple) → Dict[str, ndarray]

Generate time domain GW polarizations (h_plus, h_cross)

Parameters:

parameters_lal – A tuple of parameters for the lalsimulation waveform generator

Returns:

A dictionary of generated waveform polarizations

Return type:

pol_dict

generate_hplus_hcross(parameters: Dict[str, float], catch_waveform_errors=True) → Dict[str, ndarray]

Generate GW polarizations (h_plus, h_cross).

If the generation of the lalsimulation waveform fails with an “Input domain error”, we return NaN polarizations.

Use the domain, approximant, and mode_list specified in the constructor along with the waveform parameters to generate the waveform polarizations.

Parameters:

• parameters (Dict[str, float]) –

  A dictionary of parameter names and scalar values. The parameter dictionary must include the following keys. For masses, spins, and distance there are multiple options.

  Mass: (mass_1, mass_2) or a pair of quantities from

  ((chirp_mass, total_mass), (mass_ratio, symmetric_mass_ratio))

  Spin:

  (a_1, a_2, tilt_1, tilt_2, phi_12, phi_jl) if precessing binary or (chi_1, chi_2) if the binary has aligned spins

  Reference frequency: f_ref at which spin vectors are defined Extrinsic:

  Distance: one of (luminosity_distance, redshift, comoving_distance) Inclination: theta_jn Reference phase: phase Geocentric time: geocent_time (GPS time)

  The following parameters are not required:

  Sky location: ra, dec, Polarization angle: psi

  Units:

  Masses should be given in units of solar masses. Distance should be given in megaparsecs (Mpc). Frequencies should be given in Hz and time in seconds. Spins should be dimensionless. Angles should be in radians.

• catch_waveform_errors (bool) – Whether to catch lalsimulation errors

Returns:

A dictionary of generated waveform polarizations

Return type:

wf_dict

generate_hplus_hcross_m(parameters: Dict[str, float]) → Dict[tuple, Dict[str, ndarray]]

Generate GW polarizations (h_plus, h_cross), separated into contributions from the different modes. This method is identical to self.generate_hplus_hcross, except that it generates the individual contributions of the modes to the polarizations and sorts these according to their transformation behavior (see below), instead of returning the overall sum.

This is useful in order to treat the phase as an extrinsic parameter. Instead of {“h_plus”: hp, “h_cross”: hc}, this method returns a dict in the form of {m: {“h_plus”: hp_m, “h_cross”: hc_m} for m in [-l_max,…,0,…,l_max]}. Each key m contains the contribution to the polarization that transforms according to exp(-1j * m * phase) under phase transformations (due to the spherical harmonics).

Note:

• pol_m[m] contains contributions of the m modes and and the -m modes. This is because the frequency domain (FD) modes have a positive frequency part which transforms as exp(-1j * m * phase), while the negative frequency part transforms as exp(+1j * m * phase). Typically, one of these dominates [e.g., the (2,2) mode is dominated by the negative frequency part and the (-2,2) mode is dominated by the positive frequency part] such that the sum of (l,|m|) and (l,-|m|) modes transforms approximately as exp(1j * |m| * phase), which is e.g. used for phase marginalization in bilby/lalinference. However, this is not exact. In this method we account for this effect, such that each contribution pol_m[m] transforms exactly as exp(-1j * m * phase).

• Phase shifts contribute in two ways: Firstly via the spherical harmonics, which we account for with the exp(-1j * m * phase) transformation. Secondly, the phase determines how the PE spins transform to cartesian spins, by rotating (sx,sy) by phase. This is not accounted for in this function. Instead, the phase for computing the cartesian spins is fixed to self.spin_conversion_phase (if not None). This effectively changes the PE parameters {phi_jl, phi_12} to parameters {phi_jl_prime, phi_12_prime}. For parameter estimation, a postprocessing operation can be applied to account for this, {phi_jl_prime, phi_12_prime} -> {phi_jl, phi_12}. See also documentation of __init__ method for more information on self.spin_conversion_phase.

Differences to self.generate_hplus_hcross: - We don’t catch errors yet TODO - We don’t apply transforms yet TODO

Parameters:

parameters (dict) – Dictionary of parameters for the waveform. For details see see self.generate_hplus_hcross.

Returns:

pol_m – Dictionary with contributions to h_plus and h_cross, sorted by their transformation behaviour under phase shifts: {m: {“h_plus”: hp_m, “h_cross”: hc_m} for m in [-l_max,…,0,…,l_max]} Each contribution h_m transforms as exp(-1j * m * phase) under phase shifts (for fixed self.spin_conversion_phase, see above).

Return type:

dict

setup_mode_array(mode_list: List[Tuple]) → Dict

Define a mode array to select waveform modes to include in the polarizations from a list of modes.

Parameters:

mode_list (a list of (ell, m) modes)

Returns:

A lal parameter dictionary

Return type:

lal_params

Waveform modes

Add later.