{"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":"\n","metadata":{}},{"cell_type":"markdown","source":"# Welcome to (one more) Exploratory Data Analysis of the G2Net Dataset. \n<img src=\"https://i.ytimg.com/vi/TWqhUANNFXw/maxresdefault.jpg\" alt=\"title\">\nThe image shows a chirp pattern of gravitational waves detected by LIGO on September 14, 2015.\nCredit: LIGO (http://www.ligo.org)","metadata":{}},{"cell_type":"markdown","source":"\n\nUndoubtedly, one of the most bright breakthroughs in science in the recent decade was the detection of gravitational waves back in 2015. These waves are tiny ripples of the space-time fabric coming from collisions of some super-heavy objects, like black holes or neutron stars, predicted by Einstein 100 years before. The waves can travel billions of light years before hitting the ultra-sensitive instruments, called interferometers. Constructing such a device is the result of immense work by many scientists, engineers, and information science experts from whole over the world. \n\nThe first discovery of gravitational waves was widely outreached in the scientific community, as well as in general audience media. After the first event, more than 50 similar mergers or candidates have been reported. Now, when more and more data is collected, it is needed to more accurately detect GW signals, so that it can help build a more complete picture of our universe.\n\nIn the competition, you have a training set containing simulated time series data coming from three different locations (LIGO Hanford, LIGO Livingston, both in the US, and Virgo, in Italy). Each data file (in .npy format) represents either instrument noise or noise with a simulated gravitational wave signal, and is labeled 0 or 1, respectively. These labels are stored in a different .csv file. The task is to build a system capable of identification time series instances with GW signal present. So, this is a pure **binary classification problem** well known in the ML community. ","metadata":{}},{"cell_type":"markdown","source":"## Importing libraries\n\nFirst we need to import some libraries required to load and process the data.","metadata":{}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\npd.set_option(\"display.max_colwidth\", None) # setting the maximum width in characters when displaying pandas column. \"None\" value means unlimited.\n\nimport matplotlib.pyplot as plt  # plotting\nfrom glob import glob     # pathname management\n\nimport random    # generating (pseudo)-random numbers\n\nimport matplotlib.mlab as mlab  # some MATLAB commands\nfrom scipy.interpolate import interp1d  # interpolating a 1-D function","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:38:07.864118Z","iopub.execute_input":"2021-08-17T08:38:07.864572Z","iopub.status.idle":"2021-08-17T08:38:08.243546Z","shell.execute_reply.started":"2021-08-17T08:38:07.864477Z","shell.execute_reply":"2021-08-17T08:38:08.242581Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Importing data. General analysis","metadata":{}},{"cell_type":"code","source":"training_labels_path = '../input/g2net-gravitational-wave-detection/training_labels.csv'\ntraining_labels = pd.read_csv(training_labels_path)","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:38:17.583334Z","iopub.execute_input":"2021-08-17T08:38:17.583711Z","iopub.status.idle":"2021-08-17T08:38:18.020817Z","shell.execute_reply.started":"2021-08-17T08:38:17.583679Z","shell.execute_reply":"2021-08-17T08:38:18.019680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"training_labels.head(3)","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:38:18.979436Z","iopub.execute_input":"2021-08-17T08:38:18.979803Z","iopub.status.idle":"2021-08-17T08:38:19.004620Z","shell.execute_reply.started":"2021-08-17T08:38:18.979771Z","shell.execute_reply":"2021-08-17T08:38:19.003719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"training_labels['target'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:38:20.819279Z","iopub.execute_input":"2021-08-17T08:38:20.819669Z","iopub.status.idle":"2021-08-17T08:38:20.836262Z","shell.execute_reply.started":"2021-08-17T08:38:20.819629Z","shell.execute_reply":"2021-08-17T08:38:20.835208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can conclude that the training dataset is pretty much balanced.","metadata":{}},{"cell_type":"code","source":"training_paths = glob(\"../input/g2net-gravitational-wave-detection/train/*/*/*/*\")\nprint(\"The total number of files in the training set:\", len(training_paths))","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:38:22.571197Z","iopub.execute_input":"2021-08-17T08:38:22.571546Z","iopub.status.idle":"2021-08-17T08:39:50.909977Z","shell.execute_reply.started":"2021-08-17T08:38:22.571516Z","shell.execute_reply":"2021-08-17T08:39:50.908982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It turned out to be useful to merge labels and file paths based on their ids.","metadata":{}},{"cell_type":"code","source":"ids = [path.split(\"/\")[-1].split(\".\")[0] for path in training_paths]\npaths_df = pd.DataFrame({\"path\":training_paths, \"id\": ids})\ntrain_data = pd.merge(left=training_labels, right=paths_df, on=\"id\")","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:39:50.911674Z","iopub.execute_input":"2021-08-17T08:39:50.912070Z","iopub.status.idle":"2021-08-17T08:39:52.147444Z","shell.execute_reply.started":"2021-08-17T08:39:50.912027Z","shell.execute_reply":"2021-08-17T08:39:52.146448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data.head(3)","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:41:33.114495Z","iopub.execute_input":"2021-08-17T08:41:33.114846Z","iopub.status.idle":"2021-08-17T08:41:33.125001Z","shell.execute_reply.started":"2021-08-17T08:41:33.114814Z","shell.execute_reply":"2021-08-17T08:41:33.124186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To load a random data sample, we can make a helper function.","metadata":{}},{"cell_type":"code","source":"def load_random_file(signal = None):\n    \"\"\"Selecting a random file from the training dataset. \n    \n    Args:\n        signal: bool\n            optional flag defining whether to select pure detector \n            noise (False) or detector noise plus simulated signal (True).\n            If skipped, the flag is chosen randomly.\n    Returns:\n        file_id: str\n            unique id of the selected file\n        target: int\n            0 or 1, target value\n        data: numpy.ndarray\n            numpy array in the shape (3, 4096), where 3 is the number\n            of detectors, 4096 is number of data points (each time series\n            instance spans over 2 seconds and is sampled at 2048 Hz)\n        \n    \"\"\"    \n    if signal is None:\n        signal = random.choice([True, False])\n        \n    filtered = train_data[\"target\"]==signal   # filtering dataframe based on the target value\n    \n    index = random.choice(train_data[filtered].index)   # random index \n    \n    file_id = train_data['id'].at[index]\n    target = train_data['target'].at[index]\n    path = train_data['path'].at[index]\n    \n    data = np.load(path)\n    \n    return file_id, target, data","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:41:35.417658Z","iopub.execute_input":"2021-08-17T08:41:35.418251Z","iopub.status.idle":"2021-08-17T08:41:35.426946Z","shell.execute_reply.started":"2021-08-17T08:41:35.418194Z","shell.execute_reply":"2021-08-17T08:41:35.425819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plotting the raw data in time domain","metadata":{}},{"cell_type":"code","source":"file_id, target, data = load_random_file()\nylim = 1.1*np.max(data)\n\nplt.style.use('ggplot')\n\nfig, axs = plt.subplots(ncols=1, nrows=3, figsize=(10, 5))\n\nfor i in range(3):\n    ax = axs.ravel()[i]\n    ax.plot(data[i])\n    ax.margins(0)\n    axs[i].set_title(f\"Detector {i+1}\", loc='center')\n    ax.set_ylabel(f\"Amplitude\")\n    ax.set_ylim([-ylim, ylim])\n    \naxs[0].xaxis.set_visible(False)\naxs[1].xaxis.set_visible(False)\n\naxs[2].set_xlabel(\"Time stamp\")\nfig.suptitle(f\"Raw data visualization. ID: {file_id}. Target: {target}.\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-17T09:26:14.315494Z","iopub.execute_input":"2021-08-17T09:26:14.315845Z","iopub.status.idle":"2021-08-17T09:26:14.664789Z","shell.execute_reply.started":"2021-08-17T09:26:14.315813Z","shell.execute_reply":"2021-08-17T09:26:14.663650Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plotting the data in frequency domain\n\nOne of the ways to explore the frequency components of the data, is to plot the amplitude spectral density. To read more on this topic, please refer to the following link: https://www.gw-openscience.org/GW150914data/LOSC_Event_tutorial_GW150914.html#Whitening","metadata":{}},{"cell_type":"code","source":"fs = 2048      # sampling rate\nNFFT = 4*fs    # the Nyquist frequency \nf_min = 20.\nf_max = fs/2","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:41:45.641800Z","iopub.execute_input":"2021-08-17T08:41:45.642180Z","iopub.status.idle":"2021-08-17T08:41:45.646769Z","shell.execute_reply.started":"2021-08-17T08:41:45.642146Z","shell.execute_reply":"2021-08-17T08:41:45.645569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_, target, data = load_random_file(True)\n\nstrain1, strain2, strain3 = data[0], data[1], data[2]\n\nPxx_1, freqs = mlab.psd(strain1, Fs = fs, NFFT = NFFT)\nPxx_2, freqs = mlab.psd(strain2, Fs = fs, NFFT = NFFT)\nPxx_3, freqs = mlab.psd(strain3, Fs = fs, NFFT = NFFT)\n\npsd_1 = interp1d(freqs, Pxx_1)\npsd_2 = interp1d(freqs, Pxx_2)\npsd_3 = interp1d(freqs, Pxx_3)\n\nfig, ax = plt.subplots(ncols=1, nrows=1, figsize=(10, 5))\nax.loglog(freqs, np.sqrt(Pxx_1),\"g\",label=\"Detector 1\")\nax.loglog(freqs, np.sqrt(Pxx_2),\"r\",label=\"Detector 2\")\nax.loglog(freqs, np.sqrt(Pxx_3),\"b\",label=\"Detector 3\")\n\nax.set_xlim([f_min, f_max])\nax.set_ylabel(\"ASD (strain/$\\sqrt{Hz}$)\")\nax.set_xlabel(\"Frequency (Hz)\")\nax.legend()\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-17T09:26:05.869749Z","iopub.execute_input":"2021-08-17T09:26:05.870144Z","iopub.status.idle":"2021-08-17T09:26:07.002570Z","shell.execute_reply.started":"2021-08-17T09:26:05.870106Z","shell.execute_reply":"2021-08-17T09:26:07.001388Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Constant Q-Transform\n\nAnother very common way to visualize a GW signal is to perform a constant Q-transform (or CQT). This is a time-frequency representation widely used in processing musical data. To quickly perform Q-transform, we are going to use PyCBC library (the docs are available [here](http://pycbc.org/pycbc/latest/html/)).","metadata":{}},{"cell_type":"code","source":"!pip -q install pycbc\nimport pycbc","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:51:05.700696Z","iopub.execute_input":"2021-08-17T08:51:05.701136Z","iopub.status.idle":"2021-08-17T08:51:38.161544Z","shell.execute_reply.started":"2021-08-17T08:51:05.701098Z","shell.execute_reply":"2021-08-17T08:51:38.160436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can prepare some helper functions to generate and visualize Q-transforms. Some useful demos with PyCBC methods can be found [here](https://github.com/gwastro/PyCBC-Tutorials).","metadata":{}},{"cell_type":"code","source":"def generate_qtransform(data, fs):\n    \"\"\"Function for generating constant Q-transform. \n    \n    Args:\n        data: numpy.ndarray\n            numpy array in the shape (3, 4096), where 3 is the number\n            of detectors, 4096 is number of data points (each time series\n            instance spans over 2 seconds and is sampled at 2048 Hz)\n        fs: int\n            sampling frequency\n    Returns:\n        times: numpy.ndarray\n            array of time bins\n        freqs: numpy.ndarray\n            array of frequency bins\n        qplanes: list\n            list with 3 elements corresponding to each detector in the raw\n            data file. Each element is a 2-d vector of the power in each \n            time-frequency bin\n    \"\"\"    \n    \n    qplanes = []\n    for i in range(len(data)):\n        \n        # converting data into PyCBC Time Series format\n        ts = pycbc.types.TimeSeries(data[i, :], epoch=0, delta_t=1.0/fs)   \n        \n        # whitening the data within some frequency range\n        ts = ts.whiten(0.125, 0.125) \n        \n        # calculating CQT values\n        times, freqs, qplane = ts.qtransform(.002, logfsteps=100, qrange=(10, 10), frange=(20, 512))\n\n        qplanes.append(qplane)\n        \n    return times, freqs, qplanes ","metadata":{"execution":{"iopub.status.busy":"2021-08-17T08:55:39.537614Z","iopub.execute_input":"2021-08-17T08:55:39.538078Z","iopub.status.idle":"2021-08-17T08:55:39.546245Z","shell.execute_reply.started":"2021-08-17T08:55:39.538035Z","shell.execute_reply":"2021-08-17T08:55:39.545147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_qtransform(file_id, target, data):\n    \"\"\"Plotting constant Q-transform data.\n    \n    Args:\n        file_id: str\n            unique id of the selected file\n        target: int\n            0 or 1, target value\n        data: numpy.ndarray\n            numpy array in the shape (3, 4096), where 3 is the number\n            of detectors, 4096 is number of data points (each time series\n            instance spans over 2 seconds and is sampled at 2048 Hz)\n    \"\"\"\n    \n    times, freqs, qplanes = generate_qtransform(data, fs=fs)\n    \n    fig, axs = plt.subplots(ncols=1, nrows=3, figsize=(12, 8))\n\n    for i in range(3):\n\n        axs[i].pcolormesh(times, freqs, qplanes[i], shading = 'auto')\n        axs[i].set_yscale('log')\n        axs[i].set_ylabel('Frequency (Hz)')\n        axs[i].set_xlabel('Time (s)')\n        axs[i].set_title(f\"Detector {i+1}\", loc='left')\n        axs[i].grid(False)\n\n    axs[0].xaxis.set_visible(False)\n    axs[1].xaxis.set_visible(False)\n\n    fig.suptitle(f\"Q transform visualization. ID: {file_id}. Target: {target}.\", fontsize=16)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-17T09:20:18.128668Z","iopub.execute_input":"2021-08-17T09:20:18.129168Z","iopub.status.idle":"2021-08-17T09:20:18.137109Z","shell.execute_reply.started":"2021-08-17T09:20:18.129135Z","shell.execute_reply":"2021-08-17T09:20:18.136173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can select a random data file, perform CQT and plot the results.","metadata":{}},{"cell_type":"code","source":"file_id, target, data = load_random_file()\nplot_qtransform(file_id, target, data)","metadata":{"execution":{"iopub.status.busy":"2021-08-17T09:27:29.875385Z","iopub.execute_input":"2021-08-17T09:27:29.875776Z","iopub.status.idle":"2021-08-17T09:27:31.309546Z","shell.execute_reply.started":"2021-08-17T09:27:29.875743Z","shell.execute_reply":"2021-08-17T09:27:31.308466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we have a sample with a strong GW signal, characterized by a frequency chirp on the CQT spectrogram:","metadata":{}},{"cell_type":"code","source":"file_id = '7945e449f3'\ntarget = 1\ndata  = np.load(train_data[train_data['id']==file_id]['path'].values[0])\n\nplot_qtransform(file_id, target, data)","metadata":{"execution":{"iopub.status.busy":"2021-08-17T09:20:24.169901Z","iopub.execute_input":"2021-08-17T09:20:24.170253Z","iopub.status.idle":"2021-08-17T09:20:25.792189Z","shell.execute_reply.started":"2021-08-17T09:20:24.170223Z","shell.execute_reply":"2021-08-17T09:20:25.791168Z"},"trusted":true},"execution_count":null,"outputs":[]}]}