{'fraction_non_alphanumeric': 0.02708350575188281, 'fraction_numerical': 0.02033849209633369, 'mean_word_length': 4.751189909566873, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, '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}}
Using Deep Learning to Localize Gravitational Wave Sources

In this paper, we report on the construction of a deep Artificial Neural Network (ANN) to localize simulated gravitational wave signals in the sky with high accuracy. We have modelled the sky as a sphere and have considered cases where the sphere is divided into 18, 50, 128, 1024, 2048 and 4096 sectors. The sky direction of the gravitational wave source is estimated by classifying the signal into one of these sectors based on it's right ascension and declination values for each of these cases. In order to do this, we have injected simulated binary black hole gravitational wave signals of component masses sampled uniformly between 30-80 solar mass into Gaussian noise and used the whitened strain values to obtain the input features for training our ANN. We input features such as the delays in arrival times, phase differences and amplitude ratios at each of the three detectors Hanford, Livingston and Virgo, from the raw time-domain strain values as well as from analytical versions of these signals, obtained through Hilbert transformation. We show that our model is able to classify gravitational wave samples, not used in the training process, into their correct sectors with very high accuracy (>90%) for coarse angular resolution using 18, 50 and 128 sectors. We also test our localization on test samples with injection parameters of the published LIGO binary black hole merger events GW150914, GW170818 and GW170823 for 1024, 2048 and 4096 sectors and compare the result with that from BAYESTAR and Parameter Estimation (PE). In addition, we report that the time taken by our model to localize one GW signal is around 0.018 secs on 14 Intel Xeon CPU cores.


I. INTRODUCTION
The Laser Interferometer Gravitational Wave Observatory (LIGO) [1,2] and its European counterpart, Virgo [3], have made, at the time of writing, more than twenty reported detections of gravitational waves (GW) since their inception [4][5][6][7][8][9][10]. These GW signals, produced by the coalescence of binary black hole [4][5][6][7] and binary neutron star systems (BNS) [8][9][10], have made Gravitational Wave Astronomy a reality and it is now a major component of multi-messenger astronomy, alongside neutrino and cosmic ray [11][12][13] observations. With more detectors and better sensitivities, new GW sources will be detected on a regular basis and faster follow-up of electromagnetic (EM) counterparts of these mergers will provide detailed information about the origin and Machine Learning algorithms, like deep learning [35][36][37][38][39], can be used to address these issues since the only computationally intensive step is the one-time training phase that happens prior to the analysis of actual data. This enables low latency detection, as well as fast parameter estimation, making this approach possibly orders of magnitude faster than other techniques.
The important parameters used for localizing the GW sources are its Right Ascension (RA) and Declination (Dec) values which are its coordinates in the sky. It is necessary to obtain RA and Decs with minimum uncertainty so that EM telescopes can be pointed in the precise direction, in real time so that extremely directional and transient counterparts like kilonova and short gamma ray bursts following binary neutron star mergers can be observed. In O3 for example, for BNS events, the median sky localisation accuracy in terms of the 90% credible area is 120-180 deg 2 . 12-21% of BNS mergers are expected to be localised to less than 20 deg 2 [40].
In this paper we present a new method of localizing GW sources using Artificial Neural Network (ANN) [41], a very popular Deep Learning algorithm. We have modelled the sky as a sphere divided into many sectors and used our network to classify simulated GW signals injected in Gaussian noise into one of the sectors of the sphere. As the number of sectors is increased, the area of the individual sectors decreases, thereby potentially helping improve the angular resolution. This, coupled with the low computational costs and high speed of machine learning algorithms makes our approach a viable option for future localization endeavours. Section II gives a brief overview of the applications and successes of machine learning and in particular, deep learning in GW research. Section III describes the method of triangulation of three detectors which helps constrain the sky direction of the source of GW signals. Section IV discusses briefly the method of ANNs and their applications. Section V describes the sample generation code used to create the GW signals and noise samples used for training and testing. Section VI describes in detail the architecture of our ANN model, the rationale behind our chosen input features and the signal processing techniques used to extract those features. In Section VII we describe the details of the experiments carried out for GW signals with and without noise and discuss the results and accuracy obtained for different number of sectors and SNR ranges. In Section VIII we evaluate the performance of our classification algorithm and describe the metrics used for the evaluation. We also highlight ways to improve the accuracy in our future work. In Section IX we use our model to test samples simulated with published BH parameters for events GW150914, GW170818 and GW170823 injected into Gaussian noise and with advanced LIGO PSD. Lastly, we conclude in Section X with discussions and plans for future work.

II. MACHINE LEARNING IN GRAVITATIONAL WAVE RESEARCH
Machine Learning algorithms, in particular deep learning, have already found many applications in GW research and is quickly gaining popularity as a viable alternative to the more computationally intensive, template based matched filtering algorithm [42]. The reliability of deep learning as a powerful classification and regression tool is firmly established by its varied applications in myriad sectors. Apart from GW, deep learning has been applied in image processing and classification [43,44], self-driving cars, medical diagnostics etc.
In GW research, deep learning algorithms, in particular convolutional neural networks (CNN) [38] have been used to match the accuracy of matched filtering based search pipelines in signal detection and parameter estimation [45][46][47][48][49][50][51]. It has also been extensively applied in glitch classification [52][53][54][55][56][57], signal denoising [58][59][60][61] and has been used with numerical relativity simulations to detect higher order multipole waveforms for eccentric binary mergers [62]. The preferred machine learning algorithms of choice are mostly different types of neural networks, with the most popular being CNNs, used by several groups for detection, parameter estimation and glitch classification. Autoencoders, another popular deep learning algorithm has been used in denoising GW signals from both parametric and real LIGO noise, with very promising results.

III. LOCALIZATION OF GRAVITATIONAL WAVES
As mentioned in section I, it is of pivotal importance to have accurate estimates of the sky direction of GW signals in the sky the EM telescopes have small field of view and need an accurate GW localization. It is also possible to have EM counterparts from BBH mergers [63][64][65][66][67]. In this work however, we use BBHs as a test case as the signals are short and are fast to compute, but will apply our method in the future work for BNS or NSBH events. It is well known that each interferometer is more or less an all-sky monitor with very little directional information. Therefore localization isn't possible with a single detector. Having two detectors however does give us some directional information. Using wave arrival information and simple trigonometric calculations, we can get a ring of possible source directions in the sky for a detected GW signal. It is further possible to narrow down this range of possibilities by using relative amplitude and phase information from the detected Figure 1. Accuracy of localization by our ANN model, described in Section VIII. The dots here are the exact sky directions of a few of our injected GW test samples. Orange dots (three of them in this case) are injections that belong to the colored sectors but were classified into immediate neighbouring sectors. Sky blue dots represent injections classified into the two sectors. Sky directions of injections lying outside the two sectors can be further correctly classified by using a multi-labelling scheme.

signals.
The detectors are most sensitive to waves perpendicular to the plane formed by the two arms of the interferometers. As the incidence angle becomes less than the perpendicular, the sensitivity drops. The curvature of the Earth causes an angle difference of around 27 degrees between the zenith direction of LIGO Hanford and Livingston that creates the aforementioned amplitude and phase inconsistencies. Therefore, it is possible to constrain the sky direction of the source to sufficiently small portions of the sky using only the time delay, phase difference and relative amplitude information.
The localization can be further improved by using more detectors with large separation baselines between each individual interferometers. For three detectors, LIGO Hanford, LIGO Livingston and Virgo, the signals are localized via a triangulation of the detectors based on observed time delays of the signals. This gives two sky locations that are mirror images with respect the plane of the detectors. The degeneracy between the two sky directions can be broken by using the amplitude and phase information mentioned previously. There are several available algorithms for localising GW signals [68][69][70][71][72][73][74][75][76]. In this work, we use different signal processing techniques to obtain the time delays, amplitude ratios, and phase lags between the three detector signals and use them as input features of our ANN model described later in detail.

IV. ARTIFICIAL NEURAL NETWORKS
As discussed in section I, ANNs [41] are machine learning algorithms that are used to perform various tasks, including classification and regression. For this algorithm to work, the network must be 'trained' with sufficient labelled examples so that it learns the complex relationships between the input data and the output. The performance of the model is then tested on a set of data which the model has never seen before. A well-trained model will then be able to generalize well on unseen data and give high test accuracy.
The basic processing unit of an ANN is a node or a 'neuron', modelled after the neurons in our brain. The ANN is composed of stacked arrays of these neurons that take a vector of inputs ( x) and performs a linear operation on it with some weights ( w) and bias parameters (b), which are optimised during the course of training. The output of the transformation operation f ( x) = w · x + b from each layer is then passed through a non-linear activation function that restricts the neuron output to a certain range. The commonly used activation functions include rectified linear unit (ReLU), the hyperbolic tangent (tanh), softmax and the sigmoid functions. In our work, we have used the ReLU activation function, mathematically expressed as max(0, x).
The ANN consists of an input layer of neurons, followed by one or more hidden layers and an output layer as shown in Figure 2. Each neuron in these layers is connected to all the neurons in the next layer. In a classification problem, the output at each neuron in the output layer is the probability of an input data to belong to one of the several classes involved in the classification.
The flow of information from the input layer to the output layer is unidirectional and happens through a process called 'forward propagation'. The values obtained at the output nodes are then compared with their true outputs and the losses or mean squared deviations from the true outputs are calculated. To minimise these losses, the errors are back-propagated through the network and the weights and biases are adjusted at each iteration of the training in order to reduce the overall error. This is achieved through algorithms called 'gradient descent'. The weights and biases are initialized to small values and are then adjusted through backpropagation over several steps called 'epochs', until the errors reach a minimum.
There are several advantages of using ANNs or neural networks in general, over other machine learning algorithms. ANNs have the ability to learn complex nonlinear relationships between input data and output and has been found to be able to infer unseen relationships on unseen data. Also, ANNs do not impose any restrictions on the distribution of input variables which is a great advantage over other prediction techniques.

V. SAMPLE GENERATION
In this section, we give a short description of the sample generation code we used to create the training and testing datasets. We have used a modified version of the code developed by Gebhard et.al. [77] for this purpose. While Gebhard et.al. generated GW strains for only LIGO Hanford and Livingston interferometers, we have edited the code to include Virgo as well for better localization of the generated samples. Figure 3 shows an example of a gravitational wave strain sample generated by the code for the three detectors, Hanford, Livingston and Virgo. Here, the pure GW signals, shown in red,

Input Layer
Hidden Layers Output Layer were injected into random Gaussian noise and the strain signal, shown in blue was obtained.
To simulate the GW waveforms, we use a joint distribution of a set of relevant parameters as enlisted in Table I. The component black hole masses, spins, polarization angle, cosine of inclination, COA phase angle and injection SNRs are all uniformly sampled from the specified ranges. A brief explanation for each of these parameters and their significance are given in the Appendix in [77].
We concentrate on binary black hole mergers in our classification problem whose waveforms are modelled using the effective-one-body (EOB) approximant SEOBNRv4 in the time domain [78]. The waveforms are produced with a sampling rate of 2048 Hz. The samples are all uniformly distributed throughout the sky with the RA and Decs assigned as follows: RA = 2πu and Dec = sin −1 (1 − 2v), where u and v are uniformly sampled between [0,1]. The cosine of the inclination angle is uniformly sampled within the range [-1,1].
These parameters are then fed to routines in LALSUITE that performs the actual simulation and returns a tuple with the two polarization modes of the simulated waves (h + and h × ), resized to a chosen waveform length. These These 'raw' waveforms are then projected onto the antenna patterns of the detectors [79][80][81] based on their individual RA, Decs and polarization angles using PyCBC [82] functions and are time-corrected to ensure that the signals at the three detectors have the correct time offset on the basis of their relative positions.
This generates the pure signals observed by the detectors without any noise. In this work, we do not take into consideration the effect of the rotation of the Earth on it's axis. We assume that the Earth is stationary and the GW signals are uniformly distributed over the entire sky and approach the detectors from all directions. Since the sensitivity of the network detectors has a functional dependence on the angle of incidence of the GW, this effect is not entirely insignificant. We shall take this into account in our future work.

A. Background noise
After the pure signals are generated, they need to be injected into suitable background noise. In this work, we concentrate on simulated Gaussian noise whose frequency distribution is shown in Figure 4. We have used the advanced LIGO PSD, created using the Py-CBC module:aLIGOZeroDetHighPower. The simulated waveform is then added to the noise to get the dummy strain. The PSD estimate is obtained from the dummy strain and is then used to calculate the optimal matched filtering signal-to-noise ratio (SNR) for each of the detectors for a particular injection. From this, the net-  work optimal matched filtering SNR (NOMF-SNR) is calculated which is given by : The NOMF-SNR is used to calculate a rescaling factor that is multiplied to the waveforms to ensure that the resulting waveforms have the chosen injection SNR for the particular sample.
The resulting strain is then whitened and bandpassed, with the lower limit set to 18 Hz and upper limit to 500 Hz. We reject higher frequency components, since for this work, we concentrate on high and intermediate mass black holes. The bandpassed strain is then cut to a length of 0.25 secs, the choice of which is motivated by the estimated inspiral-merger-ringdown length for black hole mergers of our chosen masses and from system memory considerations as well. The distance of all the sources from the Earth is kept fixed at 1000 Mpc. Although the distance scales with the amplitudes and therefore the SNRs of the samples, it is an irrelevant parameter in this context because the waveforms are rescaled to match our chosen injection SNRs, as explained above. More details of the sample generation process can be found in [77].

B. Labelling of generated samples
Since we are treating the localization of GW signals as a classification problem, we need to first create the classes to which the signals need to be classified into. In order to do this, we model the sky as a sphere and divide it into several sectors based on ranges of RA and sine of Dec. We consider cases where the sphere is divided into 18, 50, 128, 1024, 2048 and 4096 sectors. Figure 5 shows an example where the sphere has been divided into 18 sectors and labelled according to our defined convention. We have followed a similar labelling convention for all the cases we have considered. For each of these cases, we assign labels to the GW samples according to the sectors they belong to. After the model is trained, we test its classification accuracy on a set of GW samples that the model has not been trained with. Each GW sample is assigned only one label. We shall consider the case for multilabelling in a later section. We have used the Pandas library [83] for the purpose of handling data frames and assigning the labels as explained above.
Here we report the training and test accuracy of classification for the coarse sky resolution (18, 50 and 128 sectors) and use finer resolution (1024, 2048 and 4096 sectors) for testing our model's localization with the parameters of the BBH merger events, GW150914, GW170818 and GW170823. We thereby report and comment on how well our model can localize these events with Gaussian noise and advanced LIGO PSD for three detectors.

VI. OUR METHOD
In this section, we describe the architecture of our ANN and the input features we have used to train the network in the context of localization of GW signals. We explain the method to extract these features from the raw time series data generated by the code described earlier. We consider two different scenarios, pure signals without noise; and signals injected into Gaussian noise. The architectures of our ANN and the input features are tailored for these two cases, which we describe in detail in the following sections. We have written our code in Python 3.6 and have used Keras [84] with Tensorflow backend [85] to implement our ANN model.

A. Model Architecture
In this section, we describe the architecture of our ANN models. Tables II and III list the components of the two networks for pure signals and signals plus Gaussian noise respectively. The tables consist of three columns. The first column lists the type of layer. 'Dense' refers to a fully connected layer, in which each node is connected to all other nodes in the next layer. 'Batch Norm' refers to Batch Normalization layer which normalizes the outputs of the previous activation layers by subtracting the batch mean and dividing by the batch standard deviation. It is applied to prevent the network from becoming too sensitive to the initial weights. By including this layer, we reduce the chances of getting very different losses on the test set at each run. The 'Dropout' layers ensure that the model does not 'overfit' on the training data. Overfitting is a common problem in machine learning, which happens when the model is trained too heavily on the training set, giving a very high training accuracy, but a low test set accuracy. This happens when the model picks up the particular set of details and random fluctuations in the training data too well, but is unable to generalize when presented with new data where these details and fluctuations do not apply. Dropout helps to prevent this over-training by randomly dropping a fraction of the total number of weights in that particular layer. We have used a dropout fraction of 0.2 after each hidden layer in our ANN model for signal plus Gaussian noise.
Columns 2 in Table I and II lists the output shape of the particular layer. All the layers in our network consist of vectors of dimension equal to the number of nodes in that layer. Column 3 shows the total number of parameters (weights and biases) corresponding to the particular layer. The activation functions we have used in our models are ReLU in all the hidden layers and softmax layer in the output layer.

B. Training Process
The model is trained after extracting the input vector of features from the time series strain data using several signal processing techniques as described in the next section. Before they are fed into the network, we normalize the features using (x − µ)/σ, where µ and σ are the mean and standard deviations for all the samples used in the training.
In order to ensure that the model does not overfit, we split our training set into two parts: the training and validation set. The model is evaluated on the validation set at every epoch and the loss on this set is  computed to keep track of the generalization efficiency of our model. We also use Keras' EarlyStopping module for this purpose, which monitors the loss on both the training and the validation sets and forcefully ends training when the loss on the validation fails to improve after a certain number of epochs. The weights were initialized randomly and during training, were optimized using the stochastic gradient algorithm, Adam [86]. The loss function used was categorical cross entropy. The calculated loss is backpropagated through the network to updates the weights and biases in order to minimize the loss. We had generated 12000 samples for pure GW signals without noise. We split it into a training and test set consisting of 10800 samples and 1200 samples respectively. For signals plus Gaussian noise, we have trained the network with 160000 samples and used a validation set of 40000 samples and a test set of 4000 samples. The choice of the number of samples were motivated by both memory considerations and the complexity of the problem. Since it is difficult to make models learn from noisy signals, we used a significantly larger dataset for signal plus Gaussian noise. During the training, we had used a mini-batch size of 250 for pure signals and 2000 for signals plus Gaussian noise. We had experimented with several batch sizes and found these yielded the best results.
We considered three different SNR ranges in our experiments with signal plus Gaussian noise : [50][51][52][53][54][55], [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35] and . The third case was studied by applying curriculum learning [87] in which the training begins with very high SNRs and then as the training progresses, the SNRs are gradually decreased. This is a form of transfer learning, in which a network uses the tuned weights and biases from earlier training, and applies it to the current training session. This approach helped improve the classification performance at low SNRs while retaining its efficiency at high SNRs. Curriculum learning has already been applied in classification, regression and denoising of GW signals previously [46,48,50,58,61,[88][89][90]. We describe the results of the classification for all the three cases in later sections.

C. Input Features
The most crucial step involved in this work is the choice of robust and effective input features that our ANN has to learn from in order to classify GW signals to their correct sectors. Our choice is directly motivated by existing methods of triangulation of three detectors, as explained earlier. The three most important requirements for accurate localization of GW signals are the arrival time delays of the signals at the three detectors, the phase lags and the amplitudes observed at the three detectors. While it is easy to extract this information from signals without noise, the problem becomes much more difficult for noisy signals because information from the actual signal could be lost due to the random fluctuations from terrestrial and other sources. Therefore, on top of the raw time-series strain data, we applied several robust signal processing techniques to extract the features of the signals buried in noise. This method can be generalized to a variety of sources with different black hole masses, spins etc. The techniques are described as follows:

Cross-correlations of time domain signals
A common technique in signal processing used to measure the similarity between two signals u(t) and v(t) as a function of the time lag between the them is called cross-correlations, which effectively is a sliding inner product of the two vectors. The cross-correlation between the signals u(t) and v(t) is given by: where u * (τ ) denotes the complex conjugate of u(τ ) and t is the time lag between the two signals. Calculating the cross-correlation gives a peak at the time lag where u(τ ) best matches v(τ + t).
Since we are dealing with three signals at the three detectors: Hanford, Livingston and Virgo, with different arrival times relative to each other and different amplitudes, we calculate three cross-correlations: between signals at Hanford (H1) and Livingston (L1), Livingston and Virgo (V1) and Hanford and Virgo. The time corresponding to the maximum cross-correlation values give the arrival time delay of one signal relative to the other.
We use the maximum cross-correlation values themselves as the amplitude information necessary to break the degeneracy between two possible sky directions of the source, as explained earlier. Therefore, we use direct cross-correlations of the time series signals to obtain a total of six feature vectors : • Arrival time delays : H1-L1, L1-V1, H1-V1.

Hilbert transforms
Since the arrival time delays and maximum crosscorrelation values of the time-domain signals are sufficient to estimate the sources of GW signals, we get a very high classification accuracy by running our ANN on pure signals using these features. But for noisy signals, the accuracy drops significantly because of the distortion of the signal due to the randomness of the noise. Therefore, besides the cross-correlations, we make use of the analytic representations obtained through a Hilbert transformation [91,92]. The Hilbert transformation converts a time domain signal into a vector of complex numbers and allows us to easily extract quantities like the instantaneous amplitudes, instantaneous phases and phase lags between two signals, from noisy signals in the time-domain.
The analytic signal of s(t) with its Fourier transform S(f ) is given by : where sgn(f ) is the signum function and F −1 is the inverse Fourier transform. This can be expressed as: where * denotes convolution operation. From here we get, where H[s(t)] is the Hilbert transform of s(t).
The analytic forms of the detector signals (H1 a , L1 a , V1 a ) are therefore arrays of 512 complex numbers (all the GW samples generated are 0.25 s long with a sampling frequency of 2048 Hz) whose real part gives the instantaneous amplitude and complex part, the instantaneous phase.
We used cross-correlations of the analytic signals to obtain the arrival time lags and amplitude information from the three pairs of detectors, the same way as we had done for the time domain signals.
We calculate the average instantaneous phase around the merger time of the three signals from the complex part of the Hilbert transform and find the difference of these average phases for each pair of detectors. We first find the instantaneous amplitudes around the merger region. By searching for ten highest peaks from each GW sample for the three detectors and record their instantaneous phases. Then we calculate the average of the ten instantaneous phases. This gives us the average instantaneous phase around the merger of each signal. We then calculate the differences of these instantaneous phases for each pair of detectors to get the instantaneous phase lags.
Having already found the instantaneous amplitudes around merger for each signal, we calculate their average and take the ratios for each pair of detectors. This serves as an additional feature for the ANN. The importance of using features around merger time is that since the signal is loudest in this region, the effect of noise is minimal, especially for high SNRs. Hence the ANN is able to learn the characteristics of the pure signal much better than while using the entire signal.
We therefore obtain 12 more input features from this method of using analytic signals via Hilbert transforms: • Arrival time delays: H1 a -L1 a , L1 a -V1 a , H1 a -V1 a • Maximum cross-correlations: H1 a ⊗L1 a , L1 a ⊗V1 a , H1 a ⊗V1 a .
Note, to compute the analytic signals using the Hilbert transforms, we used Scipy's [83] signal.hilbert method.

Complex correlation coefficients
We discussed earlier that for the purpose of localising the GW signals, the three interferometers are triangulated, which means while operating together, they form a two dimensional plane in space. GW signals may arrive at the interferometers at any angle with respect to that plane. However, since the interferometers themselves are located some distance from each other, the curvature of the Earth will introduce slight differences in the angles between the signals received at the three detectors. This quantity can therefore also be thought of as a measure of the phase difference of the waves. We calculate the angles between the signal vectors and use them as the final input feature for our ANN. To calculate the angles, we compute the complex correlation coefficient of each pair of signals in the following way: For complex-valued sinusoids, for example, of the form x(t) = Ae j(ωt+ψ) and y(t) = Be j(ωt+φ) , the angle θ between them is given by: where x and y have N samples each. We perform this operation on H1, L1 and V1 signals using Scipy's [83] libraries and get the final set of input features for signals + Gaussian noise: • Complex correlation coefficients = H1 * L1, L1 * V1, H1 * V1.

A. Signal without noise
Our first experiment was with pure signals. We used this as a sanity check for the validity of our approach, since the ANN is certainly not expected to classify noisy signals correctly if it fails for signals without noise. Since the GW signals are not contaminated by noise, we only use input features from the actual signals rather than their analytic representations. These features are the following: • Arrival time difference at each pair of detectors.
• Relative phase differences between each pair of signals.
• Ratios of amplitude at merger time. (taken as the time corresponding to highest amplitude of the signal).
We generated a sample of 12000 GW signals, distributed uniformly over the sky and split it into the training and test sets having 10800 and 1200 samples respectively. We used fixed masses, spins, COA phase angles, inclinations and polarization angles for these experiments. The architecture of the ANN we used is shown in Table II. We used a mini-batch size of 250 and trained the ANN for 300 epochs. We repeated this experiment for three different sector numbers: 18, 50 and 128. The training and test accuracy for the three cases are shown in Table IV.
We see that both the training and test accuracy are very high for all the three cases, implying our model is able to learn the relationships of the input vector of features with the individual sectors. It is also able to generalize well on unseen data from the test set, resulting in the high test accuracy. We observe that the test accuracy drops as the number of sectors is increased. This is because as the sectors become finer, it gets more and more difficult for the model to pick the correct sector from the ones in its immediate vicinity since the arrival times lags, phase lags and relative amplitudes of the GW signals from the correct sector and those around it become more and more similar in value for closely separated regions of the sky. Besides, the angular resolution is also limited by the SNRs and timing resolution. It must also be noted that this test accuracy corresponds to the number of GW samples that are classified 'exactly' correctly. In the next section, for signals plus Gaussian noise, we consider a second measure of test accuracy where samples that are classified to the immediate neighbouring sectors are counted as correctly classified. We find that the test accuracy we obtain by implementing this scheme is above 90% for all our experiments.  Table III shows the architecture of our ANN for signal plus Gaussian noise. The network here is deeper with greater number of nodes in each hidden layer to allow the ANN to learn characteristics of the actual signal from the noisy background. To make the training better, we generated 200000 GW samples of which we used 160000 samples for training and 20000 samples for validation. The ANN was tested on a final, unseen, 4000 samples. We generated our samples with variable component masses, spins etc, randomly sampled within the range as shown in Table I. We considered two different SNR ranges : [50][51][52][53][54][55] and [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35]. For each of these SNR ranges, we have evaluated our model on 18, 50 and 128 sectors each.
In addition, in order to improve the performance of the algorithm at low SNRs, we used the curriculum learning method, explained earlier, where we begin training with high SNR GW samples and then gradually decrease the SNR to 10. We evaluated the performance of our ANN, trained with the curriculum learning approach, using a test set of 4000 samples of SNR in the range of [20][21][22][23][24][25][26][27][28][29][30] and found a 3% increase in accuracy for 128 sectors.
We have used the full set of 21 input features for our experiments with noisy signals, which are the following: • Arrival time delays of signals.
• Maximum cross-correlation values of signals.
• Arrival time delays of analytic (signal plus Hilbert transform of the signal) signal.
• Maximum cross-correlation values of analytic signals.
• Ratios of average instantaneous amplitudes around merger.
• Average phase lags around merger.
• Complex correlation coefficients between the three signals.
As explained earlier, the test accuracy in all the cases improves once we consider the GWs that are classified into the immediate neighbouring sectors are counted as correctly classified. This proves that our model is indeed reliable. Figures 6 and 7 show the accuracy and loss curves in the training of the network for 18, 50 and 128 sectors for SNR ranges of 50-55 and 20-35 respectively. From the training and validation loss curves, we can conclude that our model is a good fit for the data for all the cases that we have considered. Both the training and validation loss curves in each of the Figures 6 and 7 decreases to a point of stability after the first few epochs after which they remain almost constant. The fact that our model does not overfit on the training data is confirmed by the validation loss curve remaining constant until the last epoch and it does not show an inflection point. If the model had started to overfit, the losses in the validation set would have climbed after reaching a minimum. Note the usual or preferable scenario is to always have a training loss curve lower than the validation loss curve. The fact that our validation loss is lower than the training loss curve in our plots suggests that we have an unrepresentative validation set, meaning that the validation dataset does not provide sufficient information to evaluate the ability of the model to generalize. This is expected because of the large disparity in the volume of the training and validation sets. While we have 160000 training examples, there are only 40000 validation samples. However since we have established that the model does not overfit, we can safely ignore this behaviour.

VIII. EVALUATION OF PERFORMANCE
We observe that the training and the test set accuracy of our ANN is high for both pure and noisy signals, indicating that our model is robust and is able to generalize well on examples not used in training. However, as the number of sectors are increased, the accuracy for both the training and test sets decrease. This is because as the area of the sectors gets smaller, it becomes more and more difficult for the model to differentiate between GW signals coming from adjacent sectors. Signals from adjacent sectors, would have nearly the same arrival time delays, phase lags and relative amplitudes when the individual sectors are small. Therefore the model is confused between the 'correct' sector and the ones lying adjacent to it, thereby decreasing the accuracy. In order to confirm this, we picked those incorrectly classified GW signals from our test set which were located at the sectors adjacent to the correct one and included them with our correctly classified samples. The revised test accuracy we obtained in doing so was greater than 90% for all our test cases, thereby proving our hypothesis and establishing the efficiency of our ANN algorithm. The last column of Table V shows this revised test accuracy. Fig 1 shows the performance of the ANN network on a test set comprising 4000 simulated gravitational wave samples generated with SNRs uniformly sampled between 20 and 35, injected into Gaussian noise. These samples are distributed uniformly over the entire sky and has been classified into 50 sectors based on their RA and Dec values, according the convention described earlier. The colored patches in the sphere represent two of those 50 sectors. The dots are the actual sky directions of the GWs plotted at their exact RA and Decs. The blue dots represent GW samples that the ANN classifies into the coloured sectors. Therefore, the blue dots lying inside the coloured patches represent samples correctly classified by the network and those outside it are incorrectly classified. We observe that the actual sky directions of almost all of these incorrectly classified signals are very close to the sectors the ANN classifies them into, thereby indicating that our model is indeed very accurate and is therefore a confirmation of our high 'Revised Test Accuracy Scores' in Table V.
The orange dots represent samples which belong to the coloured sectors, but were classified to the immediate neighbouring sectors by our ANN. The incorrect sectors assigned to these GW samples by our ANN lie just adjacent to the sectors they actually belong to, reaffirming our ANN's high classification accuracy.
Instead of counting the number of misclassified GW signals lying in the sectors adjacent to the correct one, we can directly implement Keras's 'multilabelling' scheme in our code for the same purpose. The idea behind multilabelling is to assign more than one labels or in this case, sectors, to each of our GW samples. Therefore each GW sample can be labelled with the 'exactly correct' sector ID as well as with the IDs of all the sectors adjacent to it. This will not only result in high test accuracy, as we have demonstrated, but high training accuracy as well.

IX. TEST RESULT FOR GW150914, GW170818 AND GW170823-LIKE SOURCES
In Section VIII we have established that our model is able to correctly classify GW samples not used in training with more than 90% accuracy for SNRs ≥ 20 and number of sectors ≤ 128. In this section, we investigate how well our model is able to localize the O1 event GW150914 and the O2 events GW170818 and GW170823. Note our training samples are generated with advanced LIGO PSD for three detectors since the purpose of this investigation is a sanity check to test the feasibility of our approach and to confirm whether in an ideal situation, our method gives good results.  [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35] In order to test this, we first simulate the GW signals with the published BBH parameters as listed in Table VI for GW150914 and Tables VII and VIII for GW170818 and GW170823 events. Next, we consider a large number of sectors: 1024, 2048, 4096 and train our model with samples of SNR in the range of [20][21][22][23][24][25] for the GW150914 test sample and with SNR in the range of [10][11][12][13][14][15] for GW170818 and GW170823 test samples. We then test our model's classification on them and use the probability assigned to each sector by our ANN for these three test samples as ranking statistics from which we obtain the cumulative frequency distribution of the probabilities and get the 90% and 50% contours in the same way as in [93]. show probability heatmaps and confidence intervals respectively, for our test samples with 2048 sectors. They have been plotted using the Python package 'healpy', that is based on the Hierarchical Equal Area isoLatitude Pixelization (HEALPix) library. For plots generated using healpy, the sky is divided into 'N' pixels. The size and number of pixels depend on our preferred resolution. Higher resolution skymaps have finer and greater number of pixels. Each of these pixels have a specific RA and Dec and a probability associated with it for the event to be located in that particular pixel or coordinate in the sky. Since the sectors in our definition have a finite area and is not simply a point or pixel in the skymap, we assign all the pixels within our sectors the same value of probability equal to the probability that was assigned to the sector by our model. Heatmaps are then obtained from these probabilities, with darker colours denoting higher probability. The 50% and 90% contour intervals are then plotted from the cumulative distribution of these probability values. The plots were made using ligo.skymaps software developed by Leo P. Singer [93].
We can see from Figures 8 and 9 that our model has been able to localize the signal quite well, implying . The reason for this large difference, we believe, is because of insufficient training data in the sectors around the source location. It must also be kept in mind that our localization is limited by our choice of sampling rate of 2048 Hz for our generated signals. With higher sampling rate, we can obtain the time delay with greater precision and therefore localize our source better. We shall take into consideration the effect of higher sampling rate in our future work. Another probable reason may be the unreliability of some of our extracted features at low SNRs. For example, we assume that information like the amplitude ratios around merger regions that we use as one of our feature vectors becomes ambiguous at SNRs < 12. However, we also find that at the same SNR, the 90% area of our GW170823 test sample is 429 deg 2 in comparison with Bayestar's area of 2145 deg 2 and PE's area of 1651 deg 2 [95]. One of the major advantages of our model, and also of all ML algorithms in general is that once the model is trained, it can perform classifications or any other operation on the test set at very high speeds, which is of particular importance for EM follow-up, as discussed earlier. We find that our trained model is able to perform the classification on a single test sample in 0.0003 secs while running on 14 Intel Xeon CPU cores.
The total time of execution, taking into consideration the time to extract the input features from the test sample's raw strain data and classifying it into one of 8192 sectors is 0.0189 secs on 14 CPU cores. Therefore our model is possibly orders of magnitude faster than any other existing localization techniques, since Bayestar takes 10 secs -1 min and other offline PE algorithms take hours to perform the localization.
We have also performed a convergence test of our localization method on the GW150914, GW170818 and GW170823 test samples for SNR in the range of 15 to 50. Figure 10 shows our convergence plot of 90% contour area vs. SNR for the GW150914 test sample. This area is expected to decrease as the SNR is increased by following an inverse square relation, as expected from diffraction limit [94]. The curve has been fitted with the function A/x 2 . We find that the 90% areas of our test sample does follow this curve as the SNR is increased. We have performed this convergence test for SNR in the range of 10 to 50 for our GW170818 and GW170823 test samples as well and we find that the 90% areas follow a similar trend.
For all three of our test samples, the 90% areas decreased as we increased the number of sectors from 1024 to 2048 but slightly increased when we increased the sector number to 4096. We suspect that this is due to the fact that with finer sectors, our training samples are distributed over a larger number of sectors and hence there aren't enough training examples in the sectors around the event to enable the model to accurately learn the features from that part of the sky. As discussed previously, the angular resolution of the source is also limited by our choice of timing resolution of 1/2048 s. Besides, in a classification problem, as the number of classes are increased, the number of error values that are backpropagated through the network during the training also increases and the weights and biases are updated every time. This might lead to instability if the number of classes (sectors in this case) are too large. Since the network may not be tuned to handle so many optimization operations in each epoch of the training, it may affect its performance.

X. DISCUSSION AND CONCLUSION
We have attempted to localize GW signals injected into Gaussian noise using deep learning. We have built an ANN and have approached it like a classification problem, where the sky, modelled as a sphere has been divided into several sectors and we train our ANN model so that it classifies simulated GW signals into their correct sector based on it's RA and Dec's. We report high (> 90%) classification accuracy for all the cases with coarse angular resolution down to a few hundredths of the sky. To achieve finer angular resolution, we calculate confidence levels using ML probability as the detection statistics. We tested our model's localization accuracy by using test samples with injection parameters of GW150914, GW170818 and GW170823 and conclude that our model gives feasible results for Gaussian noise and advanced LIGO PSD. We demonstrate that the 90% contour area for our test samples decrease in an inverse proportion to the square of SNRs, as expected from the diffraction limit. The major advantage of our approach as we have established in this work is that it is orders of magnitude faster than any current localization technique. We would like to emphasize that this work is meant to be a feasibility test of our method. We plan to try with real LIGO noise and realistic PSD and compare our results with Bayestar and PE in our future work. We also plan to work with higher numbers of sectors and ensure that the model has enough training examples for higher sector numbers by having a larger and more uniform training dataset. Another potential improvement of the localization can be done by increasing the data sampling rate. We will explore the performance of our network for binary neutron star mergers and lower mass black holes in the future.