{'fraction_non_alphanumeric': 0.0318307157415178, 'fraction_numerical': 0.015016185109699078, 'mean_word_length': 4.704393913489485, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Rapid localization of gravitational wave sources from compact binary coalescences using deep learning

The mergers of neutron star-neutron star and neutron star-black hole binaries are the most promising gravitational wave events with electromagnetic counterparts. The rapid detection, localization and simultaneous multi-messenger follow-up of these sources is of primary importance in the upcoming science runs of the LIGO-Virgo-KAGRA Collaboration. While prompt electromagnetic counterparts during binary mergers can last less than two seconds, the time scales of existing localization methods that use Bayesian techniques, varies from seconds to days. In this paper, we propose the first deep learning-based approach for rapid and accurate sky localization of all types of binary coalescences, including neutron star-neutron star and neutron star-black hole binaries for the first time. Specifically, we train and test a normalizing flow model on matched-filtering output from gravitational wave searches. Our model produces sky direction posteriors in milliseconds using a single P100 GPU, which is three to six orders of magnitude faster than Bayesian techniques.

In 2015, the first detection of a gravitational wave (GW) signal from a merging binary black hole (BBH) [1] system was made by the LIGO Scientific Collaboration [2] and the Virgo Collaboration [3]. Since then, the KA-GRA [4] detector has joined the GW search effort and the total number of GW detections has grown, at the time of writing, to 90 [5]. This includes two confirmed detections from merging binary neutron stars (BNS) and at least four from mergers of neutron star-black hole (NSBH) binaries [5]. The first BNS detection in 2017 also involved a coincident observation of a short-gamma ray burst (GRB) originating from the same source [6][7][8], making it the first multi-messenger observation involving both GW and electromagnetic (EM) signals.
These prompt EM emissions, like the short GRBs, could last for a few seconds only, and therefore rapid localization of the GW sources will be essential to enable fast EM follow-up observations by ground and spacebased telescopes [9][10][11].
In GW data analysis, a full Bayesian framework, involving Markov Chain Monte Carlo (MCMC) and nested sampling, contained within the code LALINFERENCE [12] has been traditionally used to estimate optimal source parameters of detected GW events [5]. These tools are, however, very computationally intensive and time-consuming, since they probe the full 15-dimensional parameter space that typically characterises GW signals. The rapid, Bayesian, non-MCMC localization-only code called BAYESTAR [13] can generate posterior distributions of right ascension (α) and declination (δ) angles in a few seconds. Other fast parameter estimation methods involving the use of heterodyned likelihood [14] and excess power methods [15] also exist. Most of these techniques involve solving for the posterior distribution p(θ|d) of the source parameters θ, given detector strain data d, using Bayes' theorem: where p(d|θ) is the likelihood distribution of the observed strain data, assumed to consist of a GW signal h(θ) plus stationary Gaussian noise n. p(θ) represents our prior belief about the distribution of the source parameters and p(d) is a normalizing factor. On the other hand, deep learning techniques, in particular, neural posterior estimation (NPE) [16][17][18] models have proven to be efficient at estimating accurate posterior distribution of BBH source parameters [19][20][21], enabling sky direction inference at much faster speeds than conventional sampling algorithms like MCMC. However, there have been no published results, at the time of writing, on deep learning-based posterior estimation of BNS or NSBH sources, which are the primary compact binary coalescence (CBC) candidates with potential EM coun-terparts.
The main challenge involving the use of deep learning models for BNS and NSBH parameter estimation is that the GW waveforms of these sources are much longer than BBHs and can last up to a thousand seconds at detector design sensitivity. This means that for the same overall signal-to-noise ratio (SNR), BNS and NSBH GW signals will have lower signal strength in each individual time segment, compared to BBH. This makes feature extraction using deep learning models for BNS and NSBH strains extremely challenging and is one of the possible reasons why no parameter estimation results on BNS and NSBH using deep learning has been demonstrated yet.
While all previous deep learning-based methods for parameter estimation used the GW strain data directly for training and testing, we propose an alternative approach, in which we train our neural network on the signal-to-noise ratio (SNR) time series data, obtained by the template-based matched filtering technique used for detecting binary coalescence [22] in existing online and offline GW search pipelines [23][24][25][26][27].
The SNR time series of a GW signal can be thought of as a 'condensed' version of the GW strain data, in which all information about the GW source sky direction is preserved in its amplitude and a constant phase [28], and is typically less than 0.1 sec around the merger time. The top panel in Fig. 1 shows an example BNS waveform and the corresponding strain obtained after injecting the signal in stationary Gaussian noise with advanced LIGO Power Spectral Density (PSD). The bottom panel shows the optimal SNR time series obtained by matched filtering the strain data with the optimal signal template. The time difference of the SNR peaks between pairs of detectors equals the arrival time delays of the signal, which is essential information for source localisation based on the triangulation method [29]. The sensitivity of each detector to incoming GW signals is a function of the sky direction of the source. This information about the detector response further helps constrain the possible sky directions of the source and is also encoded in the amplitude and the constant phase of the SNR time series [28].
Since the SNR time series data is readily available from matched filtering-based GW search pipelines, we can directly apply our model in conjuction with a GW detection pipeline to perform rapid sky localization immediately after the detection of a GW signal.
Matched filtering involves cross-correlation of the GW strain data, d, with waveform templates h(θ) that depend on the intrinsic source parameters θ which include masses and spins. The output of matched filtering is the SNR time series ρ(t), defined as [28], where z(t) is the complex matched-filter defined as [28], Here d(f ) is the Fourier transform of the GW strain, h c (f ) and h s (f ) are the Fourier transforms of the cosine and sine components of the template used for matched filtering, and S n (f ) is the PSD of the noise. σ is the standard deviation of the real part of the matched filter, which is used to normalize ρ(t). In this work, we train and test our model using the 'optimal' SNR time series data, (ρ opt (t)), which is obtained by performing matched filtering of the GW strain with its optimal signal template, or the template with parameters exactly equal to that of the source components. Matched filtering of GW strains with the optimal signal template produces the highest peak in the SNR time series output, compared to when the filtering is done with non-optimal templates.
The dependence of ρ opt on important source parameters like α, δ, distance etc., can be found in eq. 7.181 in [22]. Since our model is trained and tested on very short (0.2 sec long) ρ opt (t) data, we are able to obtain highly accurate localization results for all CBC sources with a relatively simple neural network architecture, compared to deeper and more sophisticated neural networks used in similar analyses [19,20]. The use of short SNR time series data instead of strains for training and inference was the crucial reason for the success of our model on BNS and NSBH events, since machine learning models usually struggle to extract features from hundreds of seconds long input data, even after robust feature engineering, data compression or the use of extremely deep neural networks.
Our localization model is called CBC-SkyNet or Compact Binary Coalescence -Sky Localization Neural Network. CBC-SkyNet consists of three main components: two identical ResNet-34 networks [30] built using one-dimensional (1D) convolutional and max-pooling layers, to extract essential features from real and imaginary parts of ρ opt (t), and a Normalizing Flow [31] model, in particular, a Masked Autoregressive Flow (MAF) network [32,33], to approximate the true posterior distribution of α and δ (p true (α, δ|ρ opt )). The two 1D ResNet-34 models are used to extract essential features from the real and imaginary parts of ρ opt (t) which are then passed as an input feature vector to our normalizing flow model.
A normalizing flow is a chain of invertible transformations with simple Jacobian determinant that maps a simple base distribution to a more complex target distribution. The motivation behind using a Normalizing Flow is to produce flexible posterior distributions, since the shape of p true (α, δ|ρ opt ) cannot be accurately approximated using a single or a even a mixture of Gaussian distributions. Recent work using normalizing flows for full 15 dimensional parameter estimation of BBH sources has generated near-optimal results at only a few seconds [19].
If z is a random sample drawn from the base distribu-FIG. 1: A BNS waveform, the strain signal obtained after injecting the waveform in noise, and the corresponding SNR time series obtained by matched filtering the strain with the optimal waveform template is shown. Top panel: The whitened GW strain signal corresponding to a BNS merger is plotted in blue. This data was obtained after injecting the GW waveform (rescaled amplitude shown in orange) into stationary Gaussian noise coloured by advanced LIGO PSD. Bottom panel: The optimal SNR time series obtained by matched filtering the GW strain with the signal template shown in the above panel.
tion π(z) and f is the transformation induced by the Normalizing Flow, then the new random variable is x = f (z).
Since f is invertible, we can write z = f −1 (x). The probability distribution of the new random variable p(x) is then given by the change of variables formula for probability distributions: We use a multivariate Gaussian as the base distribution, and the MAF as the transformation f , to obtain our target posterior p φ (α, δ|ρ opt ). Since our target distribution p φ (α, δ|ρ opt ) is a conditional probability distribution, the transformation f is conditioned on features extracted from ρ opt (t) using the 1D ResNet-34 networks. The architectures of these networks are discussed in the Methods section.
The basic unit of the MAF is called a Masked Autoencoder for Distribution Estimation (MADE) network [34], which is a fully connected neural network with several connections that are masked or set to zero. This is done to ensure the network obeys the 'autoregressive' property which has shown substantial success in density estimation and sample generation. Fig. 2 shows the workflow of the CBC-SkyNet. The MAF takes as input random samples (z) drawn from our base distribution, which is a multivariate Gaussian with three dimensions, equal to the number of parameters we intend to predict. The α parameter is cyclic, i.e. for a fixed δ, α = 0 and α = 2π points to the same location in the sky. In order to make our model learn this cyclical nature of α, and prevent it from mistaking α = 0 with α = 2π and vice-versa, we decompose α into α x and α y , equal to the cosine and sine components of α respectively. Therefore, instead of predicting α and δ directly, we predict α x , α y and δ. We found that doing this trick improved the performance of the model on injection samples located at α = 0 and 2π. The final α posterior is obtained by taking the arc tangent of cos α and sin α. The Normalizing Flow model then scales and shifts the z samples to obtain p φ (α, δ|ρ opt ), the model's approximation of p true (α, δ|ρ opt ).
The objective of inference using NPE algorithms like Normalizing Flows is to approximate the true posterior distribution p true (θ|ρ opt ) of source parameters θ with a neural network without having to evaluate a likelihood. In order to do that, we optimize a loss function which in our case, is defined as the expectation over ρ opt of the cross entropy between p true (θ|ρ opt ) and the neural network predicted posterior p φ (θ|ρ opt ), L = dρ opt p true (ρ opt ) dθ p true (θ|ρ opt )log p φ (θ|ρ opt ), (5) where φ are the network parameters. Using Bayes' theorem and Monte Carlo approximation, we obtain the following approximate expression for the loss function: where N is the size of a mini-batch of training examples. The θ (i) and ρ (i) opt samples are drawn from the prior θ ∼ p true (θ) and likelihood ρ opt ∼ p true (ρ opt |θ). Equation (6) shows that the optimization step is free from costly likelihood and posterior evaluations. Instead, we only need to draw samples from the likelihood, which is readily available as matched filtering outputs from online GW search pipelines [23][24][25][26][27]. The optimization of L is done using stochastic gradient descent [35] to calculate the gradients with respect to φ.
Under the transformation given in eq. 4, the loss function, eq. 6 becomes: where n refers to the number of parameters/labels being predicted.
We train CBC-SkyNet with 10 5 (2 × 10 5 ) simulated SNR time series samples each, for three (two) detector BBH, BNS and NSBH events and perform injection run tests on 2000 (4000) samples each, drawn from the same prior distribution of source parameters as training (waveform approximants and prior distributions are defined in Extended Data Table I). We train and test on samples drawn uniformly from detector network SNR range 10 to 20 (10 to 30 for two detectors) for each CBC source, following realistic detection scenarios. We first generated the strain samples by injecting simulated waveforms in stationary Gaussian noise, sampled at 2048 Hz, coloured by advanced LIGO ZeroDetHighPower PSD and then performed matched filtering using the optimal template to obtain ρ opt (t).
We train two different versions of CBC-SkyNet for three detectors (Hanford (H1), Livingston (L1), Virgo (V1)), and two detector (H1 and L1) cases. The architectures of the three detector and two detector networks differ in only their choice of hyper-parameters and network depth. (see Methods section). For inference, we draw 5000 α and δ posterior samples for each event in our test set and construct P-P plots and cumulative histograms of the areas of 90% and 50% credible intervals, and of the searched areas, defined as the area of the smallest contour that contains the true sky direction of the events. The details of the sample generation process, training and testing are described in the Methods section.

A. Skymaps of test samples
We discuss our model's performance on simulated GW events and obtain skymaps for three and two detector cases. Fig. 3 shows CBC-SkyNet's localization performance on simulated BBH, BNS and NSBH events for both two (H1 and L1) and three detector (H1, L1 and V1). The true sky co-ordinates of all these events (indicated by the black cross) lie within CBC-SkyNet's predicted 50% credible interval (inner contour in Fig. 3). The shape of the sky map contours is determined by the signal arrival time delays, modulated by the phase lag and amplitude difference at different detectors. For a three detector network, the individual equal time delay contours intersect at two patches in the sky which are mirror images with respect to the plane defined by the detector network. This degeneracy can, in principle, be broken by considering the amplitude and phase consistency of the signals between detectors, allowing us to assign different significance to each path and potentially neglect one of them. For a two detector network, using only time delays, we obtain elongated ellipses and rings in the sky, which are only partially broken by including the amplitude and phase information. This explains our three detector event skymaps ( Fig. 3(a), (b) and (c)), with the two detector event skymaps ( Fig. 3 (d), (e) and (f)).

B. Localization of GW170817-like event
We describe results obtained by testing our model on a GW170817-like injection sample. Fig. 4, shows CBC-SkyNet's performance on an injection sample with parameters of the first BNS event detected, GW170817. The strain sample for the event was generated using source parameters from [6], and injected into Gaussian noise with advanced LIGO PSD. Matched filtering was then performed to obtain the corresponding ρ opt (t). Our model's posterior distribution, denoted by the heatmap in Fig. 4 (see inset), contains our injected true sky direction (denoted by blue marker) within the 50% credible interval. The area of the 90% and 50% intervals obtained are 48 deg 2 and 14 deg 2 respectively. The shape of our predicted 90% contour is broader than BAYESTAR [36] or LALINFERENCE [6] contours. This is possibly because our simulated Virgo sample has an SNR of 7.5, while the former analyses used real LIGO data, in which Virgo had SNR less than 3. For this reason, our predicted contour resembles a three detector localization contour, giving it FIG. 2: The workflow of our localization model, CBC-SkyNet, is shown here. The input data consists of 0.2 sec long SNR time series data obtained from matched filtering. The real and imaginary parts of the SNR time series is passed through identical ResNet-34 networks built using 1D convolutional layers. The features extracted from the two ResNet-34 model are then combined into a single 1D feature vector. This vector of features is then passed as a conditional input to the MAF that draws samples from the base distribution and transforms them to the desired posterior distribution. a broader shape, rather than a narrower and elongated contour obtained using real detector data. In the future, we plan to train and test our model on SNR time series obtained by injecting waveforms into real noise from O2 or O3. This would be a more robust test of our model's performance and make direct comparison with other analyses possible.

C. Injection run results
We have performed a full injection run for all CBC sources, for two and three detector cases each, and have plotted the cumulative distributions of the areas of different credible intervals for each run. Figs 5 (a), (b) and (c) show cumulative distributions of areas of the 90% credible interval, 50% credible intervals and searched area contours for 2000 simulated three detector BBH, BNS and NSBH events respectively, with network SNR be-tween 10 and 20. The median areas of the 90% contour we obtained for BBH is consistent with anticipated sky localization areas for O3 and O4, as shown in Table 5 of [37]. For BNS and NSBH, our median 90% contour areas are 9% and 4.3% smaller than the lower limits reported in [37] for these sources. One possible reason for this is the analysis in [37] considered O3 detector sensitivity for sample generation, which is much noisier than the advanced LIGO PSD we used in our analysis. There are also differences in the distribution of masses considered between our analysis and that of [37]. While we consider uniform mass distributions within certain ranges (Table I in Methods), [37] uses a power law distribution for BH masses and a Gaussian distribution for the NS masses with mean 1.33 and standard deviation 0.09. Although the masses of the binary components only affects the width of the SNR time series peak, and is not expected to have a very significant impact on localization accuracy, having similar source parameter distributions in future analyses will lead to better comparison of our results with other methods. In terms of inference time, we are again orders of magnitude faster, since it takes us a few milli-seconds to sample from our predicted posterior distribution, as opposed to hours for LALINFERENCE, as used in [37]. Figs 5 (d) -(f) show CDFs for 2000 two detector BBH, BNS and NSBH samples, involving H1 and L1 only. We compare our two detector BNS results with [38] and find that our median 90% credible interval area is larger than LALINFERENCE and BAYESTAR by 67% and 28% respectively (see Table 2 in [38]). This may be due to the fact that for the MAF, two detector localization posteriors are hard to approximate, since they have more complex shapes (elongated ellipses or broken rings) than three detector contours, which are usually single patches in the sky. However, it must also be noted that we use a more conservative lower SNR threshold of 10 in our injections, compared to 12 in [38]. Our upper SNR limit was chosen to be 30, whereas around 5% of the injections in [38] have SNR > 30.
Our two detector model also shows higher prediction uncertainty for BBH events, as evidenced by the broader 1-sigma uncertainty region in the 90% credible interval area CDF in Fig. 5 (d), compared to BNS and NSBH (Fig. 5 (e) and (f)). This possibly indicates that the two detector model needs further hyper-parameter tuning to make it robust against different CBC source data. Figs 6 (a), (b) and (c) show Probability-Probability (P-P) plots for the three detector BBH, BNS and NSBH injection runs respectively. The P-P plots show the cumulative histograms of the true sky directions of events lying in different credible intervals. It serves as a test of consistency of any Bayesian parameter estimation method with the frequentist interpretation (i.e., on average the true sky directions should lie within the 90% interval for 90% of the test cases). The distribution should therefore follow the dotted diagonal line (with small deviations due to finite sample size, indicated by the 90% error margin in grey). As the plots indicate, our pre-dicted P-P curves are consistent within the 90% error bounds of the expected diagonal line for all the sources. Figs 6 (d) -(f) show P-P plots for two detector BBH, BNS and NSBH respectively. For the two detector case, our predicted distributions show similar behaviour to the three detector scenario, although slight deviations from the diagonal line especially from the 60% credible interval are observed for all three sources. This further indicates that better optimization of our two detector network is needed.

II. DISCUSSION
In summary, we have reported the first deep learning based approach for sky localization of all CBC sources that enables orders of magnitude faster inference than Bayesian methods using limited computational resources. Using SNR time series data instead of the actual GW strain data for training and testing, we were able to achieve localization of BNS and NSBH sources for the first time using deep learning. The use of SNR time series data also significantly reduces computational overheads and memory, needed in most deep learning-based approaches. We plan to adapt this method for training and testing on real noise injections. This would require injecting simulated GW waveforms in real detector noise prior to performing matched filtering. Besides sky localization, this method can be adopted for inference of other important parameters essential for multi-messenger follow-up, like chirp mass and luminosity distance.

A. Sample generation
To train and test our model, we simulate nonprecessing BBH, BNS and NSBH waveforms and inject them into stationary Gaussian noise coloured by advanced LIGO PSD to get the detector strains. We also sample the desired network SNR of the events from the prior and rescale the amplitudes of our simulated strains to match the network SNR that we sample from the prior. We then perform matched filtering on the generated samples to obtain ρ opt (t) for each detector. We simulate h(θ) using source parameters θ drawn from the prior distribution of source parameters listed in Extended Data Table I. For BBH samples, we use the waveform approximant SEOBNRv4 [39] and SpinTaylorT4 [40] for BNS and NSBH. All the waveforms were simulated with a fixed coalescence GPS time of 1187008882.4. The generated waveforms are then injected into Gaussian noise and the resultant strains, sampled at 2048 Hz, are whitened with the advanced LIGO Zero Detuned High Power PSD using the Welch median average method [41]. After whitening, we rescale the amplitudes of the strains to match a desired network SNR. In order to do that, we first calcu-late the network optimal matched filtering SNR (NOMF-SNR) of the whitened waveforms using the formula: NOMF-SNR = SNR(H1) 2 + SNR(L1) 2 + SNR(V1) 2 (8) We then use the NOMF-SNR to calculate a scaling factor that we use to re-scale the strain amplitudes to match the injection SNR drawn from the prior network SNR distribution. After rescaling the strains, matched filtering is performed using the optimal templates to obtain ρ opt (t). We chop ρ opt (t) for all CBC sources to 0.2 secs around the merger peak. We train our model for 55 epochs each for BBH, BNS and NSBH. Because of the low dimensionality of our input data, training is very quick, requiring less than an hour on a NVIDIA Tesla P100 GPU. During testing, we sample from the predicted distribution to obtain posterior samples of α and δ and then evaluate a 2-D Kernel Density Estimator to generate the joint probability density of α and δ. Overall this sampling and evaluation step takes less than a second. Sample generation and matched filtering was implemented with a modified version of the code developed by [42] that uses PyCBC software [43]. CBC-SkyNet was written in Tensor-Flow 2.4 [44] using the Python language. The code is available in: GitHub.

B. Masked Autoregressive Flow
In this section, we review Masked Autoregressive Flows and the autoregressive density estimation technique. The MAF is a type of normalizing flow that maps a simple base distribution to a more complex distribution. In autoregressive density estimation techniques, a model learns a complex joint density by decomposing it into a product of one-dimensional conditional densities as follows [32,33]: In autoregressive models, these conditional densities are usually parameterized as a univariate Gaussian. The means and standard deviations of the Gaussians are predicted by neural networks that depend on the previous x 1:i−1 : where, and for i = 1,...,n.
f µi and f αi are the transformations induced by the MAF. If we sample u ∼ N (0, 1) n , then, defines a sample from p(x). This defines a transformation of the base distribution, which is a univariate normal, to the target distribution p(x). We can stack multiple autoregressive transformations to create a sufficiently expressive flow. It is easy to compute the inverse transformation u = f −1 (x) since it only requires us to reverse the scale and shift as follows: Therefore the forward and inverse pass of the flow both require only the forward evaluation of f µ and f α . The Jacobian in eq. 7 can be calculated as follows: det For parameterizing a normalizing flow, two conditions need to be necessarily fulfilled: the inverse of the transformation must exist and the Jacobian must be easy to compute. The above arguments suggest that the MAF satisfied both the requirements necessary for a normalizing flow. The autoregressive property is implemented in neural networks by 'masking' the weights of a fully connected neural network, or by leaving out only those connections which are allowed, with the rest being set to zero. This architecture is called the MADE or Masked Autoencoder for Density Estimation. Usually to obtain a sufficiently flexible model, a stack of several MADE networks are used. Between each pair of MADE networks, the order of the components are permuted so that if one network is unable to model a distribution well due to poor choice of variable ordering, a subsequent network with a different combination of weights and masking is able to do it. In our MAF implementation, we also use a batch normalization layer between each MADE unit to stabilize training and prevent the performance of the model from suffering due to spurious random weight initializations.

C. Network Architectures
In this section, we describe the architectures of our 1D ResNet-34 and the MAF models. We use two identical 1D ResNet-34 networks for the real and imaginary parts of the SNR time series input. Our ResNet-34 architecture is similar to [30] but with 1D convolutional and max-pooling layers, instead of the traditional 2D layers. The ResNet-34 is a deep network with 34 layers which has a special feature called the skip connections. The skip connection connects the output of one layer to another by skipping some layers in between. These layers, together with the skip connection forms a residual block. ResNets are made by stacking these residual blocks together. The advantage of using several skip connections is that the network starts to learn the features of the input data right from the start of training, even if some layers have not started learning, i.e., their weights are close to the values they were initialized with (usually close to zero). The skip connections make it possible for the signal to flow through the whole network and prevent training degradation and the vanishing gradients problem usually observed in very deep neural networks.

Three detectors
The first two layers of 1D ResNet-34 network are a 1D convolutional and with ReLU activation function and a batch normalization layer. We use 16 filters for the first convolutional layer with a kernel size of 7 and stride of 2. The stride determines the length by which the receptive field of a convolutional kernel shifts over the activations of the previous layer. We then use a max-pooling layer with a pool size of 3 to reduce the input dimensions by a factor of 3. This is followed by a stack of residual units. After every few residual units, the number of filters is doubled, starting from 16. Each residual unit is comprised of two convolutional layers with batch normalization, ReLU activation and kernel size of 3. Overall, our ResNet-34 implementation consists of 34 layers (only counting the convolutional and fully connected layers) containing 3 residual units with 16 filters, followed by 4 residual units with 32 filters, followed by 6 residual units with 164 filters, and then 3 residual units with 128 filters. Similar architecture is used for the imaginary input branch. The architecture of the ResNet networks is tabulated in Extended Data Table 2. The features extracted from the real and imaginary branches are passed through two fully connected layers each with 256 and 128 neurons respectively, following which the features from the two 1D ResNet branches are concatenated into a single vector of features. This is then passed as conditional input to the MAF.
The workflow of the MAF has already been described in the previous section. As described, we use the MADE network to implement the MAF. We use a stack of six MADE units, each consisting of three hidden layers with 256 neurons. In between each pair of MADE networks, we use batch normalization to stabilize training and produce a more flexible model. The same 1D ResNet and MAF architectures are used for BBH, BNS and NSBH sources.

Two detectors
For the two detector case, we used a 30 layer 1D ResNet model instead of the 34 layer 1D ResNet that we used in the three detector case. We also used 36 filters in the initial convolutional layers for the two detector case, compared to 16 for three detectors. At the end of the real and imaginary 1D ResNet branches, we added a 1D global average pooling layer which takes the two dimensional features generated from the previous layers and performs a downsampling operation by returning the average of the time or spatial dimension of the input. We then concatenate the final vector of features from the real and imaginary 1D ResNet branches and then perform another downsampling operation by passing the concatenated features through a fully connected layer of 256 neurons, which is then passed as a conditional input to the MAF. The 1D ResNets therefore consist of 3 residual units with 36 filters, followed by 4 residual units with 72 filters, followed by 4 residual units with 144 filters, and finally 3 residual units with 288 filters. Each residual unit is again comprised of two convolutional layers with batch normalization, ReLU activation and kernel size of 3. The architecture of the ResNet networks is tabulated in Extended Data Table 3.
We also tested with different MAF network architectures and found the best performance with 5 MAF hidden layers with 1024 neurons each for BBH and NSBH and 5 MAF hidden layers 2048 neurons each for BNS. Both the three and two detector networks were trained for 55 epochs with a batch size of 2048 using the Adam optimizer [35]. The learning rate was set to 10 −4 throughout the training. We arrived at this value after testing the model's performance on a range of initial learning rates and investigating the under-and over-fitting effects. The total training time was less than an hour with a P100 GPU.