{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Exploratory Data Analysis G2Net Gravitational Wave Detection 2021\n\nThis notebook covers basic EDA for the G2Net Gravitational Wave Detection competition. Also included: spectrogram and power spectral density (PSD) transformations and their visualizations using [GWpy library](https://gwpy.github.io/docs/stable/). Brief definitions of important concepts are also provided. References such as notebooks or papers used are listed below.\n\n**Contents**\n\n* [1. Libraries](#section-1)\n* [2. Count labels in training set](#section-2)\n* [3. Import and visualize training data as time series](#section-3)\n    - [3.1 Time series example, target=1](#subsection-3_1)\n    - [3.2 Time series example, target=0](#subsection-3_2)\n* [4. Represent training data as spectrograms](#section-4)\n    - [4.1 Spectrogram example, target=1](#subsection-4_1)\n    - [4.2 Spectrogram example, target=0](#subsection-4_2)\n* [5. Represent data as power spectral density (PSD)](#section-5)\n    - [5.1 PSD example, target=1](#subsection-5_1)\n    - [5.2 PSD example, target=0](#subsection-5_2)\n\n\n\n\n**References**\n\nNotebooks: \n\n(1) [ [EDA] G2Net GW Detection by Mukharbek Organokov ](https://www.kaggle.com/muhakabartay/eda-g2net-gw-detection).\n\n(2) [G2Net Keras baseline [Training] by Ammar Alhaj Ali](https://www.kaggle.com/ammarnassanalhajali/g2net-keras-baseline-training).\n\nPapers:\n\n[1] [Bailes et al. *Gravitational-wave physics and astronomy in the 2020s and 2030s*, Nature (2021) ](https://www.nature.com/articles/s42254-021-00303-8).\n\nLibraries:\n\n{1} [*GWpy*: Duncan Macleod, Alex L. Urban, Scott Coughlin, Thomas Massinger, Matt Pitkin, rngeorge, … The Gitter Badger. (2021, April 27). gwpy/gwpy: 2.0.4 (Version v2.0.4). Zenodo. http://doi.org/10.5281/zenodo.472331](https://gwpy.github.io/docs/stable/).","metadata":{}},{"cell_type":"markdown","source":"<a id=\"section-1\"></a>\n## 1. Libraries \n\nFirst we install [GWpy](https://gwpy.github.io/docs/latest/index.html) which \"is a collaboration-driven Python package providing tools for studying data from ground-based gravitational-wave detectors\" and then load other libraries.","metadata":{}},{"cell_type":"code","source":"!python -m pip install gwpy ","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-07-14T21:35:36.519390Z","iopub.execute_input":"2021-07-14T21:35:36.519693Z","iopub.status.idle":"2021-07-14T21:35:51.575638Z","shell.execute_reply.started":"2021-07-14T21:35:36.519629Z","shell.execute_reply":"2021-07-14T21:35:51.574565Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport random\n\nimport numpy as np\nimport pandas as pd\n\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\nfrom gwpy.timeseries import TimeSeries, TimeSeriesDict\nfrom gwpy.plot import Plot","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:35:51.577077Z","iopub.execute_input":"2021-07-14T21:35:51.577354Z","iopub.status.idle":"2021-07-14T21:35:53.566413Z","shell.execute_reply.started":"2021-07-14T21:35:51.577322Z","shell.execute_reply":"2021-07-14T21:35:53.564664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section-2\"></a>\n## 2. Count labels in training set","metadata":{}},{"cell_type":"code","source":"PATH='../input/g2net-gravitational-wave-detection/'\n\ndef get_train_file_path(event_id):\n    \"\"\"From an event id returns the path to the corresponding file\"\"\"\n    file=\"train/{}/{}/{}/{}.npy\".format(event_id[0],event_id[1],event_id[2],event_id) \n    return file\n\ntrain_labels = pd.read_csv(os.path.join(PATH,'training_labels.csv'))\n\ntrain_labels['file_path'] = [get_train_file_path(event_id=event) for event in train_labels['id']] \n\ntrain_labels.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:35:53.567818Z","iopub.execute_input":"2021-07-14T21:35:53.568167Z","iopub.status.idle":"2021-07-14T21:35:54.261348Z","shell.execute_reply.started":"2021-07-14T21:35:53.568140Z","shell.execute_reply":"2021-07-14T21:35:54.260604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(7,6))\nax = sns.countplot( x = train_labels['target'] )\n","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:35:54.262590Z","iopub.execute_input":"2021-07-14T21:35:54.262817Z","iopub.status.idle":"2021-07-14T21:35:54.607706Z","shell.execute_reply.started":"2021-07-14T21:35:54.262792Z","shell.execute_reply":"2021-07-14T21:35:54.606640Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get random event paths for target 0 and 1\n\ndef get_random_sample(df,value):\n    \"Gets a random event id with target=value=(0,1)\"\n    ids=df.loc[df['target'] == value]['id'].tolist()\n    return random.sample(ids,1)[0]\n\n# get a random sample with target = 1\nevent_1 = get_random_sample(df=train_labels,value=1)\npath_1=train_labels.loc[train_labels['id'] == event_1]['file_path'].to_string(index=False,buf=None)\nprint('random_event(target=1): ', event_1)\nprint(path_1, '\\n')\n\n# get a random sample with target = 0\nevent_0 = get_random_sample(df=train_labels,value=0)\npath_0=train_labels.loc[train_labels['id'] == event_0]['file_path'].to_string(index=False)\nprint('random_event(target=0): ', event_0)\nprint(path_0, '\\n')\n\n","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:35:54.610088Z","iopub.execute_input":"2021-07-14T21:35:54.610290Z","iopub.status.idle":"2021-07-14T21:35:54.783920Z","shell.execute_reply.started":"2021-07-14T21:35:54.610269Z","shell.execute_reply":"2021-07-14T21:35:54.782928Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section-3\"></a>\n## 3. Import and visualize training data as time series","metadata":{}},{"cell_type":"markdown","source":"From the data description: \"Each data sample (npy file) contains 3 time series (1 for each detector) and each spans 2 sec and is sampled at 2,048 Hz\". This means that each time series is a ndarray with shape (3, 4096): one time series per GW interferometer containing 2 x 2048 = 4096 time points. We will check this with an example. The amplitude is measured by a dimensionless parameter called *strain*, a full description can be found in the paper [1].","metadata":{}},{"cell_type":"code","source":"# functions\n\ndef get_npy_data(file):\n    \"\"\"Get .npy file content.\"\"\"\n    data_npy = np.load(file)\n    Hanford    = TimeSeries(data_npy[0,:], sample_rate=2048)\n    Livingston = TimeSeries(data_npy[1,:], sample_rate=2048)\n    Virgo      = TimeSeries(data_npy[2,:], sample_rate=2048)\n    return Hanford, Livingston, Virgo\n\ndef plot_time_series(Hanford, Livingston, Virgo):\n    \"\"\"Plot time series. Separate subplots.\n    Detectors order: LIGO Hanford, LIGO Livingston, and Virgo\n    \"\"\"\n    plot = Plot(Hanford, Livingston, Virgo, \n                separate=True, \n                sharex=True,\n                ylabel='Strain',\n                figsize=[18, 12])\n    ax = plot.gca()\n    ax.set_xlim(0,2)\n    ax.set_xlabel('Time [s]')\n    plot.show()\n    \ndef plot_time_series_all(Hanford, Livingston, Virgo):\n    \"\"\"Plot time series. All detectors together.\"\"\"\n    plot = Plot(figsize=(18, 4))\n    ax = plot.add_subplot()\n    ax.plot(Hanford, color='gwpy:ligo-hanford', label='LIGO-Hanford')\n    ax.plot(Livingston, color='gwpy:ligo-livingston', label='LIGO-Livingston')\n    ax.plot(Virgo, color='gwpy:virgo', label='Virgo')\n    ax.set_ylabel('Strain')\n    ax.set_xlim(0, 2)\n    ax.set_ylim(-2e-20, 2e-20)\n    ax.set_xlabel('Time [s]')\n    ax.legend()\n    \ndef plot_time_series_LIGO(Hanford, Livingston):\n    \"\"\"Plot time series. LIGO detectors.\"\"\"\n    plot = Plot(figsize=(18, 4))\n    ax = plot.add_subplot()\n    ax.plot(Hanford, color='gwpy:ligo-hanford', label='LIGO-Hanford')\n    ax.plot(Livingston, color='gwpy:ligo-livingston', label='LIGO-Livingston')\n    ax.set_ylabel('Strain')\n    ax.set_xlim(0, 2)\n    ax.set_ylim(-2e-20, 2e-20)\n    ax.set_xlabel('Time [s]')\n    ax.legend()","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:35:54.784895Z","iopub.execute_input":"2021-07-14T21:35:54.785126Z","iopub.status.idle":"2021-07-14T21:35:54.793896Z","shell.execute_reply.started":"2021-07-14T21:35:54.785102Z","shell.execute_reply":"2021-07-14T21:35:54.793374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"subsection-3_1\"></a>\n### 3.1 Time series example, target=1 : ","metadata":{}},{"cell_type":"code","source":"# data sample target 1\nfile_npy_ev_1 = os.path.join(PATH,path_1)\nprint( file_npy_ev_1, '\\n' )\nprint( np.load(file_npy_ev_1).shape , '\\n')\n# plot data sample\nHanford, Livingston, Virgo = get_npy_data(file_npy_ev_1)\nplot_time_series(Hanford, Livingston, Virgo)\nplot_time_series_LIGO(Hanford, Livingston)\nplot_time_series_all(Hanford, Livingston, Virgo)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:35:54.794663Z","iopub.execute_input":"2021-07-14T21:35:54.794956Z","iopub.status.idle":"2021-07-14T21:35:58.999033Z","shell.execute_reply.started":"2021-07-14T21:35:54.794927Z","shell.execute_reply":"2021-07-14T21:35:58.998261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"subsection-3_2\"></a>\n### 3.2 Time series example, target=0 : ","metadata":{}},{"cell_type":"code","source":"# data sample target 1\nfile_npy_ev_0 = os.path.join(PATH,path_0)\nprint( file_npy_ev_0, '\\n' )\nprint( np.load(file_npy_ev_0).shape , '\\n')\n# plot data sample\nHanford, Livingston, Virgo = get_npy_data(file_npy_ev_0)\nplot_time_series(Hanford, Livingston, Virgo)\nplot_time_series_LIGO(Hanford, Livingston)\nplot_time_series_all(Hanford, Livingston, Virgo)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:35:59.001297Z","iopub.execute_input":"2021-07-14T21:35:59.001834Z","iopub.status.idle":"2021-07-14T21:36:03.354238Z","shell.execute_reply.started":"2021-07-14T21:35:59.001797Z","shell.execute_reply":"2021-07-14T21:36:03.353308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section-4\"></a>\n## 4. Represent training data as spectrograms\nA [spectrogram ](https://en.wikipedia.org/wiki/Spectrogram) \"is a visual representation of the spectrum of frequencies of a signal as it varies with time\". GWpy provides a way to generate a spectrogram from a time series using a [Q-transform or Constant-Q transform](https://en.wikipedia.org/wiki/Constant-Q_transform), which basically \"transforms a data series to the frequency domain\". This Q-transform is used to \"to produce high-resolution time-frequency maps of transient noise (glitches) and potential gravitational-wave signals\", [see GWpy documentation on this topic](https://gwpy.github.io/docs/latest/examples/timeseries/qscan.html).","metadata":{}},{"cell_type":"code","source":"def plot_spectrograms(detector_data, detector_name, event, target, grid):\n    \"\"\"Plot spectrograms for specific detector's event.\"\"\"\n    qspecgram = detector_data.q_transform(outseg=(0.0, 2.0))\n    plot = qspecgram.plot(figsize=[6, 4])\n    ax = plot.gca()\n    ax.set_yscale('log')\n    ax.set_title('{}. Event: {}. Target: {}'.format(detector_name, event, target))\n    ax.set_xlabel('Time [s]')\n    ax.set_ylabel('Frequency [Hz]')\n    ax.grid(grid)\n    ax.colorbar(cmap='viridis', label='Normalized energy')\n    plot.show()\n\ndef plot_all_spectrograms(detectors, names, event, target, grid):\n    \"\"\"From a list o data from detectors and a list of \n    their names plots the spectrograms of the three \n    time series, it also requires the event id, known \n    target and an option to add a grid (True/False)\"\"\"\n    for detector, name in zip(detectors,names):\n        plot_spectrograms(detector, name, event, target, grid)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:36:03.355310Z","iopub.execute_input":"2021-07-14T21:36:03.355513Z","iopub.status.idle":"2021-07-14T21:36:03.362070Z","shell.execute_reply.started":"2021-07-14T21:36:03.355492Z","shell.execute_reply":"2021-07-14T21:36:03.361015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"subsection-4_1\"></a>\n### 4.1 Spectrogram example, target=1 : ","metadata":{}},{"cell_type":"code","source":"# get the data (same data sample as in the previous section)\nHanford, Livingston, Virgo = get_npy_data(file_npy_ev_1)\ninterferometers = [Hanford, Livingston, Virgo] \ninternames = ['LIGO Hanford', 'LIGO Livingston', 'Virgo']\ntarget_sample = train_labels[train_labels['id'] == event_1]['target'].iloc[0]\n# plot the spectrograms\nplot_all_spectrograms(detectors=interferometers,names=internames,event=event_1,target=target_sample, grid=False)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:36:03.362966Z","iopub.execute_input":"2021-07-14T21:36:03.363156Z","iopub.status.idle":"2021-07-14T21:36:08.075879Z","shell.execute_reply.started":"2021-07-14T21:36:03.363135Z","shell.execute_reply":"2021-07-14T21:36:08.075040Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"subsection-4_2\"></a>\n### 4.2 Spectrogram example, target=0 : ","metadata":{}},{"cell_type":"code","source":"# get the data (same data sample as in the previous section)\nHanford, Livingston, Virgo = get_npy_data(file_npy_ev_0)\ninterferometers = [Hanford, Livingston, Virgo] \ninternames = ['LIGO Hanford', 'LIGO Livingston', 'Virgo']\ntarget_sample = train_labels[train_labels['id'] == event_0]['target'].iloc[0]\n# plot the spectrograms\nplot_all_spectrograms(detectors=interferometers,names=internames,event=event_0,target=target_sample, grid=False)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:36:08.076956Z","iopub.execute_input":"2021-07-14T21:36:08.077227Z","iopub.status.idle":"2021-07-14T21:36:13.220429Z","shell.execute_reply.started":"2021-07-14T21:36:08.077199Z","shell.execute_reply":"2021-07-14T21:36:13.219093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section-5\"></a>\n## 5. Represent data as power spectral density (PSD)\nPower spectral density (PSD) is used to represent [spectral density](https://en.wikipedia.org/wiki/Spectral_density), which \"describes the distribution of power into frequency components composing that signal\", PSD \"refers to the spectral energy distribution that would be found per unit time\". GWpy provides a method for [spectral density estimation from a time series called psd()](https://gwpy.github.io/docs/latest/signal/index.html). \n","metadata":{}},{"cell_type":"code","source":"def plot_PSD(detector_data, detector_name, event, target, grid, clr='black'):\n    \"\"\"Plot spectrograms for specific detector's event.\"\"\"\n    PSD = detector_data.psd()\n    plot = PSD.plot(figsize=[10, 8], color=clr)\n    ax = plot.gca()\n    ax.set_title('{}. Event: {}. Target: {}'.format(detector_name, event, target))\n    ax.set_xlabel('Frequency [Hz]')\n    ax.set_ylabel(r'PSD [strain$/\\sqrt{\\mathrm{Hz}}$]')\n    ax.grid(grid)\n    plot.show()\n\ndef plot_all_PSD(detectors, names, event, target, grid):\n    \"\"\"From a list o data from detectors and a list of \n    their names plots the spectrograms of the three \n    time series, it also requires the event id, known \n    target and an option to add a grid (True/False)\"\"\"\n    colors=['red','blue','purple']\n    for detector, name, colr in zip(detectors,names,colors):\n        plot_PSD(detector, name, event, target, grid, clr=colr)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:36:13.221448Z","iopub.execute_input":"2021-07-14T21:36:13.221695Z","iopub.status.idle":"2021-07-14T21:36:13.231677Z","shell.execute_reply.started":"2021-07-14T21:36:13.221673Z","shell.execute_reply":"2021-07-14T21:36:13.228367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"subsection-5_1\"></a>\n### 5.1 PSD example, target=1 : ","metadata":{}},{"cell_type":"code","source":"# plot all PSD for the three detectors\nHanford, Livingston, Virgo = get_npy_data(file_npy_ev_1)\ninterferometers = [Hanford, Livingston, Virgo] \ninternames = ['LIGO Hanford', 'LIGO Livingston', 'Virgo']\ntarget_sample = train_labels[train_labels['id'] == event_1]['target'].iloc[0]\nplot_all_PSD(detectors=interferometers,names=internames,event=event_1,target=target_sample, grid=True)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:36:13.233544Z","iopub.execute_input":"2021-07-14T21:36:13.233783Z","iopub.status.idle":"2021-07-14T21:36:15.008087Z","shell.execute_reply.started":"2021-07-14T21:36:13.233760Z","shell.execute_reply":"2021-07-14T21:36:15.007582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"subsection-5_2\"></a>\n### 5.2 PSD example, target=0 : ","metadata":{}},{"cell_type":"code","source":"# plot all PSD for the three detectors\nHanford, Livingston, Virgo = get_npy_data(file_npy_ev_0)\ninterferometers = [Hanford, Livingston, Virgo] \ninternames = ['LIGO Hanford', 'LIGO Livingston', 'Virgo']\ntarget_sample = train_labels[train_labels['id'] == event_0]['target'].iloc[0]\nplot_all_PSD(detectors=interferometers,names=internames,event=event_0,target=target_sample, grid=True)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T21:36:15.008926Z","iopub.execute_input":"2021-07-14T21:36:15.009229Z","iopub.status.idle":"2021-07-14T21:36:16.564421Z","shell.execute_reply.started":"2021-07-14T21:36:15.009208Z","shell.execute_reply":"2021-07-14T21:36:16.563491Z"},"trusted":true},"execution_count":null,"outputs":[]}]}