{"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":"![](https://www.sciencenews.org/wp-content/uploads/2016/02/gwaves_feat-new.jpg)\n\n\nIn this competition you are provided with a training set of time series data containing simulated gravitational wave measurements from a network of 3 gravitational wave interferometers (LIGO Hanford, LIGO Livingston, and Virgo). Each time series contains either detector noise or detector noise plus a simulated gravitational wave signal. The task is to identify when a signal is present in the data (target=1).\n\nSo we need to use the training data along with the target value to build our model and make predictions on the test IDs in form of probability that the target exists for that ID.\n\nSo basically data-science helping here by building models to filter out this noises from data-streams (which includes both noise frequencies and Gravitational Waves frequencies) so we can single out frequencies for Gravitational-Waves. This is very well-explained by by Professor Rana Adhikari of Caltech and a member of the LIGO team, who were the first to measure gravitational waves. See his [interview here](https://youtu.be/1D2j8nTjOZ4?t=1946)\n\n\n# Basic Description of the Data Provided\n\nWe are provided with a train and test set of time series data containing simulated gravitational wave measurements from a network of 3 gravitational wave interferometers:\n\n- LIGO Hanford\n\n- LIGO Livingston\n\n- Virgo\n\nEach time series contains either detector noise or detector noise plus a simulated gravitational wave signal.\n\nThe task is to identify when a signal is present in the data (target=1).\n\nEach .npy data file contains 3 time series (1 coming for each detector) and each spans 2 sec and is sampled at 2,048 Hz.\n\nAnd we have a total of 5,60,000 files, each file of dimension of 3 * 4096, this turns out to be a huge time series\n\n\n# What are Gravitational Waves?\n\nGravitational waves are 'ripples' in space-time caused by some of the most violent and energetic processes in the Universe. Albert Einstein predicted the existence of gravitational waves in 1916 in his general theory of relativity. Einstein's mathematics showed that massive accelerating objects (such as neutron stars or black holes orbiting each other) would disrupt space-time in such a way that 'waves' of undulating space-time would propagate in all directions away from the source. These cosmic ripples would travel at the speed of light, carrying with them information about their origins, as well as clues to the nature of gravity itself.\n\nThe strongest gravitational waves are produced by cataclysmic events such as colliding black holes, supernovae (massive stars exploding at the end of their lifetimes), and colliding neutron stars. Other waves are predicted to be caused by the rotation of neutron stars that are not perfect spheres, and possibly even the remnants of gravitational radiation created by the Big Bang.\n\n[Source](https://www.ligo.caltech.edu/page/what-are-gw)\n\nOn September 14, 2015 for the very first time, LIGO (Laser Interferometer Gravitational-Wave Observatory) physically sensed the undulations in spacetime caused by gravitational waves generated by two colliding black holes 1.3 billion light-years away. LIGO's discovery will go down in history as one of humanity's greatest scientific achievements.\n\nWhile the processes that generate gravitational waves can be extremely violent and destructive, by the time the waves reach Earth they are thousands of billions of times smaller! In fact, by the time gravitational waves from LIGO's first detection reached us, the amount of space-time wobbling they generated was a 1000 times smaller than the nucleus of an atom! Such inconceivably small measurements are what LIGO was designed to make.\n\n# Time Domain vs Frequency Domain Analysis ?\n\nA Time domain analysis is an analysis of physical signals, mathematical functions, or time series of economic or environmental data, in reference to time. Also, in the time domain, the signal or function's value is understood for all real numbers at various separate instances in the case of discrete-time or the case of continuous-time. Furthermore, an oscilloscope is a tool commonly used to see real-world signals in the time domain.\n\nMoreover, a time-domain graph can show how a signal changes with time.\n\nIn Frequency domain your model/system is analyzed according to it's response for different frequencies. How much of the signal lie in different frequency range. Theoretically signals are composed of many sinusoidal signals with different frequencies (Fourier series), like triangle signal, its actually composed of infinite sinusoidal signal (fundamental and odd harmonics frequencies).\n\nWe can move from time domain to frequency domain with the help of Fourier transform.\n\n# Sources of Noise in Gravitational Waves Signal data and how Data Science can help\n\nNoise of non-astrophysical origin contaminates data taken by LIGO. The sensitivity of advanced gravitational-wave detectors will be limited by multiple sources of noise from the hardware subsystems and the environment. The low frequency sensitivity of the detectors (10 Hz) will be limited by the effects of seismic noise. Thermal noise due to Brownian motion will be the most dominant noise source in the most sensitive frequency range of the instruments. At frequencies higher than ∼150 Hz, shot noise, due to quantum uncertainties in the laser light, is expected to be the dominant noise source. Instrumental and environmental disturbances can also produce non-astrophysical triggers in science data, so called ‘glitches,’ as well as increasing the false alarm rate of searches and producing a decrease in the detectors’ duty cycles. The non-Gaussian and nonstationary nature of advanced detector noise may produce glitches, which could affect the sensitivity of searches and be mistaken as gravitational-wave detections, in particular for unmodelled sources.","metadata":{}},{"cell_type":"code","source":"!python -m pip install gwpy\n# !pip install astropy==4.2.1\n!pip install astropy\n!pip install nnAudio\n\n# from google.colab import drive\n# drive.mount('/content/gdrive')\n# import os\n\n# One time ONLY - Installation of Kaggle package\n# !pip install --upgrade --force-reinstall --no-deps kaggle\n\n# TO IMPORT DATA - ONE TIME CODE\n# os.environ['KAGGLE_CONFIG_DIR']='/content/gdrive/MyDrive/Kaggle_Datasets'\n# os.chdir('/content/gdrive/MyDrive/G2Net_Gravitational_Waves')\n# !kaggle competitions download -c g2net-gravitational-wave-detection\n\n#Unzip Data\n# !unzip /content/gdrive/MyDrive/G2Net_Gravitational_Waves/g2net-gravitational-wave-detection.zip","metadata":{"execution":{"iopub.status.busy":"2021-08-23T15:57:35.140539Z","iopub.execute_input":"2021-08-23T15:57:35.140931Z","iopub.status.idle":"2021-08-23T15:58:03.840735Z","shell.execute_reply.started":"2021-08-23T15:57:35.140851Z","shell.execute_reply":"2021-08-23T15:58:03.839779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport seaborn as sns\nfrom scipy import signal\nfrom gwpy.timeseries import TimeSeries\nfrom gwpy.plot import Plot\nimport numpy as np\nfrom sklearn.preprocessing import MinMaxScaler\nfrom PIL import Image\nfrom glob import glob\nfrom matplotlib import pyplot as plt\nimport random\nfrom colorama import Fore, Back, Style\nplt.style.use('ggplot')\nimport warnings\nwarnings.filterwarnings('ignore')\nfrom sklearn.model_selection import train_test_split\n\nfrom tensorflow.keras.utils import Sequence\n\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import Dense, Dropout, Flatten, Conv1D, MaxPool1D, BatchNormalization\nfrom tensorflow.keras.optimizers import RMSprop, Adam\n\nimport torch\nfrom nnAudio.Spectrogram import CQT1992v2","metadata":{"execution":{"iopub.status.busy":"2021-08-23T15:58:03.844297Z","iopub.execute_input":"2021-08-23T15:58:03.844586Z","iopub.status.idle":"2021-08-23T15:58:12.213373Z","shell.execute_reply.started":"2021-08-23T15:58:03.844556Z","shell.execute_reply":"2021-08-23T15:58:12.21241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some basic definition\n\n### Wha are .npy files?\n\nIt is a standard binary file format for persisting a single arbitrary NumPy array on disk. The format stores all of the shape and data type information necessary to reconstruct the array correctly even on another machine with a different architecture. The format is designed to be as simple as possible while achieving its limited goals. The implementation is intended to be pure Python and distributed as part of the main numpy package.\n\n**This file format makes incredibly fast reading speed enhancement over reading from plain text or CSV files.**\n\n### Sampling frequency in a TimeSeries\n\nThe sampling frequency, or sample rate, is the number of equal-spaced samples per unit of time. For instance, if you have 96 equally spaced observation per day, then you sampling rate is 96/day, or 96/24/3600=0.0011 Hz. Hz, which means per second, is widely used for sample rate.","metadata":{}},{"cell_type":"code","source":"# Checking the contents of one file\nroot_dir = \"../input/g2net-gravitational-wave-detection/\"\nfile = root_dir + 'train/0/0/0/000a5b6e5c.npy'\ndata = np.load(file)\nprint(data.shape)\nprint(data)\n# print(data[0, :].shape)\n# print(data[1, :].shape)\n# print(data[2, :].shape)\nprint(\"data[0, :] is \", data[0, :])\n# data_1","metadata":{"execution":{"iopub.status.busy":"2021-08-23T15:58:12.215143Z","iopub.execute_input":"2021-08-23T15:58:12.215488Z","iopub.status.idle":"2021-08-23T15:58:12.239819Z","shell.execute_reply.started":"2021-08-23T15:58:12.215452Z","shell.execute_reply":"2021-08-23T15:58:12.238899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The gwpy's TimeSeries function expects array-like input data array as its first argument\n# and sample_rate : float, Quantity, optional the rate of samples per second (Hertz)\ndef get_tseries_from_file(file_name):\n  t_data = np.load(file_name)\n  tseries1 = TimeSeries(t_data[0,:], sample_rate=2048)\n  tseries2 = TimeSeries(t_data[1,:], sample_rate=2048)\n  tseries3 = TimeSeries(t_data[2,:], sample_rate=2048)\n  return tseries1, tseries2, tseries3\n\n''' Multi-data plots with gwpy.plot\n\nhttps://gwpy.github.io/docs/latest/plot/index.html#multi-data-plots - \n\nGWpy enables trivial generation of plots with multiple datasets. The Plot constructor will accept an arbitrary collection of data objects and will build a figure with the required geometry automatically. By default, a flat set of objects are shown on the same axes: \n\nseparate=True can be given to show each dataset on a separate Axes\n\nThe returned object is a Plot, a sub-class of matplotlib.figure. Figure adapted for GPS time-stamped data. Customisations of the figure or the underlying Axes can be done using standard matplotlib methods. Hence I can use .gca()\n\n.gac() - Get the current Axes instance on the current figure matching the given keyword args, or create one. The plt.gca() function gets the current axes so that you can draw on it directly. \n\nhttps://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.gca.html\n\n'''\n\ndef plot_tseries(t1, t2, t3):\n  plot = Plot(t1, t2, t3, separate=True, sharex=True, figsize=[20, 12])\n  ax = plot.gca()\n  ax.set_xlim(0, 2)\n  ax.set_xlabel('Time [s]')\n  plt.show()\n  \nfile_1 = root_dir + 'train/0/0/0/000a5b6e5c.npy'\n  \ntseries1, tseries2, tseries3 = get_tseries_from_file(file_1)\n\n# Plotting the 3 TimeSeries\nplot_tseries(tseries1, tseries2, tseries3)\n  ","metadata":{"execution":{"iopub.status.busy":"2021-08-23T15:58:12.241518Z","iopub.execute_input":"2021-08-23T15:58:12.241884Z","iopub.status.idle":"2021-08-23T15:58:16.087043Z","shell.execute_reply.started":"2021-08-23T15:58:12.241834Z","shell.execute_reply":"2021-08-23T15:58:16.08611Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load the .npy files from all the nested folder-structure and get the ids from file names","metadata":{}},{"cell_type":"code","source":"train_labels = pd.read_csv(root_dir + \"/training_labels.csv\")\nprint(Fore.YELLOW + 'Dataset has ', Style.RESET_ALL + \"{} Observations\".format(train_labels.shape[0]) )\n\nprint(Fore.MAGENTA + \"Printing first 5 Labels: \", Style.RESET_ALL )\ndisplay(train_labels.head())","metadata":{"execution":{"iopub.status.busy":"2021-08-23T15:58:16.088277Z","iopub.execute_input":"2021-08-23T15:58:16.088619Z","iopub.status.idle":"2021-08-23T15:58:16.488766Z","shell.execute_reply.started":"2021-08-23T15:58:16.088585Z","shell.execute_reply":"2021-08-23T15:58:16.487967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Construct a Training dataframe for all the available .npy files \n\n# Get all the file file path from all 4-labels of nested folder structure\nfiles_paths = glob(root_dir + '/train/*/*/*/*')\n''' The glob module finds all the pathnames matching a specified pattern according to the rules \nused by the Unix shell, although results are returned in arbitrary order. \nNo tilde expansion is done, but *, ?, and character ranges expressed with [] will be correctly matched. \nWe can use glob to search for a specific file pattern, or perhaps more usefully, search for files where the \nfilename matches a certain pattern by using wildcard characters.\n\n'''\n\n# get the list of ids from the .npy files\nids_from_npy_files = [path.split(\"/\")[-1].split(\".\")[0] for path in files_paths]\n# [-1] means the last element in a sequence,\n# print(ids_from_npy_files)\n\n# get a dataframe with paths and ids of those .npy files\ndf_path_id = pd.DataFrame({'path': files_paths, 'id':ids_from_npy_files})\ndf_path_id.head()\n\n# merging that above df with the target\ndf_train = pd.merge(left=train_labels, right=df_path_id, on='id')\ndisplay(df_train.head())\n\n# verifying the shape of the merged df has 5,60,000 rows and 3 columns\ndf_train.shape\n","metadata":{"execution":{"iopub.status.busy":"2021-08-23T15:58:16.490098Z","iopub.execute_input":"2021-08-23T15:58:16.490467Z","iopub.status.idle":"2021-08-23T16:00:24.011294Z","shell.execute_reply.started":"2021-08-23T15:58:16.490413Z","shell.execute_reply":"2021-08-23T16:00:24.01053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Classify the the 2 classes of targets of 1 and 0\ntarget_1_df_train = df_train[df_train.target == 1]\ntarget_0_df_train = df_train[df_train.target == 0]\nprint(\"Class distribution of Target: \\n \", train_labels.target.value_counts())\ndisplay(target_1_df_train.head())","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:24.01263Z","iopub.execute_input":"2021-08-23T16:00:24.012985Z","iopub.status.idle":"2021-08-23T16:00:24.142317Z","shell.execute_reply.started":"2021-08-23T16:00:24.012948Z","shell.execute_reply":"2021-08-23T16:00:24.141495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.countplot(x = 'target' , data=train_labels)\nplt.title('Target Class Distribution')","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:24.145478Z","iopub.execute_input":"2021-08-23T16:00:24.145789Z","iopub.status.idle":"2021-08-23T16:00:24.396682Z","shell.execute_reply.started":"2021-08-23T16:00:24.145759Z","shell.execute_reply":"2021-08-23T16:00:24.395611Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Basic CQT EDA and A quick Note on nnAudio and CQT\n\n### Q-Transform\n\nThe constant quality factor transform (CQT), introduced by J.C. Brown in 1988, is an interesting alternative to the windowed Fourier transform (STFT / Short Time Fourier Transform) or wavelets, for time-frequency analysis.\n\nThe constant-Q transform transforms a data series to the frequency domain. It is related to the Fourier transform.\nIn general, the transform is well suited to musical data and proves useful where frequencies span several octaves.It is more useful in the identification of instruments.\n\nUnlike the Fourier transform, but similar to the mel scale, the constant-Q transform (Wikipedia) uses a logarithmically spaced frequency axis. The original paper reference below\n\n[Judith C. Brown, \"Calculation of a constant Q spectral transform,\" J. Acoust. Soc. Am., 89(1):425–434, 1991.](http://academics.wellesley.edu/Physics/brown/pubs/cq1stPaper.pdf)\n\n[From Wikipedia](https://en.wikipedia.org/wiki/Constant-Q_transform) - In mathematics and signal processing, the constant-Q transform, simply known as CQT transforms a data series to the frequency domain. It is related to the Fourier transform[1] and very closely related to the complex Morlet wavelet transform. In general, the transform is well suited to musical data, and this can be seen in some of its advantages compared to the fast Fourier transform. As the output of the transform is effectively amplitude/phase against log frequency, fewer frequency bins are required to cover a given range effectively, and this proves useful where frequencies span several octaves. As the range of human hearing covers approximately ten octaves from 20 Hz to around 20 kHz, this reduction in output data is significant.\n\nA Constant Q transform is a variation on the Discrete Fourier Transform (DFT). In other words, it is a type of wavelet transform.\n\nI only have a casual understanding of both types of transforms myself, so take what I'm saying with a grain of salt.\n\nA standard DFT uses a constant window size throughout all frequencies. This typically leads to a pretty consistent, fully continuous transform. However, the constant bin size for all frequencies leads to some problems when you map frequency on a logarithmic scale. Specifically, peaks on the lower end are incredibly wide (sometimes up to half an octave), lacking any sort of detail.\n\nThis is an issue for emulating human perception because humans perceive frequency on a logarithmic scale.\n\nA Constant Q transform seeks to solve this problem by increasing the window size for lower frequencies, and alleviate some of the computational strain caused by this by reducing the window size used for high frequencies. It's pretty effective at this, but has a few drawbacks.\n\nThe computational complexity of a Constant Q transform is only slightly larger than that of a standard DFT, but because the window size changes per frequency, it is impossible to apply the typical optimizations of the FFT to a Constant Q transform.\n\nIn other words, a Constant Q transform will yield better results where low frequencies and logarithmic frequency mapping are concerned.\n\nThe transform exhibits a reduction in frequency resolution with higher frequency bins, which is desirable for auditory applications. The transform mirrors the human auditory system, whereby at lower-frequencies spectral resolution is better, whereas temporal resolution improves at higher frequencies.\n\nCQT refers to a time-frequency representation where the frequency bins are geometrically spaced\nand the Q-factors (ratios of the center frequencies to bandwidths) of all bins are equal.\n\n\n### nnAudio\n\n**nnAudio** is an audio processing toolbox using PyTorch convolutional neural network as its backend. By doing so, spectrograms can be generated from audio on-the-fly during neural network training and the Fourier kernels (e.g. or CQT kernels) can be trained. **Kapre** and **torch-stft** have a similar concept in which they also use 1D convolution from Keras and PyTorch to do the waveforms to spectrogram conversions. Other GPU audio processing tools are torchaudio and tf.signal. But they are not using the neural network approach, and hence the Fourier basis can not be trained.\n\n\nFrom [this](https://www.readcube.com/articles/10.1109%2Faccess.2020.3019084) Paper\n\n\"nnAudio, a new neural network-based audio processing framework with graphics processing unit (GPU) support that leverages 1D convolutional neural networks to perform time domain to frequency domain conversion. It allows on-the-fly spectrogram extraction due to its fast speed, without the need to store any spectrograms on the disk. Moreover, this approach also allows back-propagation on the waveforms-to-spectrograms transformation layer, and hence, the transformation process can be made trainable, further optimizing the waveform-to-spectrogram transformation for the specific task that the neural network is trained on. All spectrogram implementations scale as Big-O of linear time with respect to the input length. nnAudio, however, leverages the compute unified device architecture (CUDA) of 1D convolutional neural network from PyTorch, its short-time Fourier transform (STFT), Mel spectrogram, and constant-Q transform (CQT) implementations are an order of magnitude faster than other implementations using only the central processing unit (CPU).\"\n\n","metadata":{}},{"cell_type":"code","source":"# First quickly check the shape of my df_train that I earlier defined\ndf_train.head()","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:24.398712Z","iopub.execute_input":"2021-08-23T16:00:24.399067Z","iopub.status.idle":"2021-08-23T16:00:24.40925Z","shell.execute_reply.started":"2021-08-23T16:00:24.39903Z","shell.execute_reply":"2021-08-23T16:00:24.408184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"waves_from_each_file = np.load(df_train.loc[2, 'path'])\nprint(\"waves_from_each_file shape: \", waves_from_each_file.shape)\nwaves_stacked = np.hstack(waves_from_each_file)\nprint(\"HStacked waves_from_each_file shape: \", waves_stacked.shape)\n","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:24.410589Z","iopub.execute_input":"2021-08-23T16:00:24.411025Z","iopub.status.idle":"2021-08-23T16:00:24.449378Z","shell.execute_reply.started":"2021-08-23T16:00:24.410989Z","shell.execute_reply":"2021-08-23T16:00:24.448546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to show the CQT Spectogram\ndef get_cqt_spectrogram(\n    waves_from_each_file,\n    transform=CQT1992v2(sr=2048, fmin=20, fmax=1024, hop_length=64),\n):\n    stacked_waves_from_each_file = np.hstack(waves_from_each_file)\n    stacked_waves_from_each_file = stacked_waves_from_each_file / np.max(\n        stacked_waves_from_each_file\n    )\n    stacked_waves_from_each_file = torch.from_numpy(stacked_waves_from_each_file).float()\n    cqt_image = transform(stacked_waves_from_each_file)\n    return cqt_image\n\n\nfor i in range(5):\n  waves = np.load(df_train.loc[i, 'path'])\n  cqt_image = get_cqt_spectrogram(waves)\n  # print(\"cqt_image \", i, cqt_image)\n  target = df_train.loc[i, 'target']\n  plt.imshow(cqt_image[0])\n  plt.title(f\"target: {target}\")\n  plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:24.450696Z","iopub.execute_input":"2021-08-23T16:00:24.451063Z","iopub.status.idle":"2021-08-23T16:00:26.039475Z","shell.execute_reply.started":"2021-08-23T16:00:24.451026Z","shell.execute_reply":"2021-08-23T16:00:26.038491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In above **sr** represents Sampling rate.\nSampling rate or Sampling frequency is the rate at which we are capturing amplitudes. In other words, it is just the number of data points recorded per second. These amplitudes can be electric or magnetic fileds in the case of light waves, current or voltage in the case of digital signals and pressure or displacement values for the case of sound waves.","metadata":{}},{"cell_type":"markdown","source":"# Plotting (KDE, Boxplot, and Timeseries)\n\nAs we are provided with a train and test set of time series data from a network of 3 gravitational wave interferometers:\n\n- LIGO Hanford\n\n- LIGO Livingston\n\n- Virgo\n","metadata":{}},{"cell_type":"code","source":"# Defining a multi-plot function\n# And I am going to call the 3 series as Series-1, Series-2 and Series-3\ndef multi_plot(series, plot_type, target):\n    if plot_type == 'box' or plot_type == 'kde':\n      plt.figure(figsize=(20, 2))\n    else:\n      plt.figure(figsize=(15,12))\n\n    for idx in range(3):\n      if plot_type == 'box':\n          plt.subplot(1, 3, idx+1)\n          sns.boxplot(series[idx:idx+1], color = 'b')\n\n      elif plot_type == 'kde':\n            plt.subplot(1, 3, idx+1)\n            sns.kdeplot(series[idx], color = 'g', shade=True, lw=2, alpha=0.5)\n      else:\n        plt.subplot(3, 1, idx+1) # As in the case of 'time' based plot, I wan 3 rows\n        plt.plot(series[idx: idx+1].T, color= 'r')\n        plt.title(\"\\nSeries-\" + str(idx+1))\n\n    if plot_type == 'box':\n      plt.suptitle(\"Box Plots(target = \" + target + \")\")\n    elif plot_type == 'kde':\n      plt.suptitle(\"Probability Distribution Plots(target = \" + target + \")\")\n    else:\n      plt.suptitle(\"Time Distribution of Signals - Spans 2 sec, Sampled at 2,048 Hz(target = \" + target + \")\")\n        \n\n  ","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:26.040927Z","iopub.execute_input":"2021-08-23T16:00:26.041263Z","iopub.status.idle":"2021-08-23T16:00:26.051028Z","shell.execute_reply.started":"2021-08-23T16:00:26.041227Z","shell.execute_reply":"2021-08-23T16:00:26.049812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extract 1 random element from the target_1_df_train\n# using random_state to ensure the reproducibility of the examples.\ntarget_1_rand_sample_file = target_1_df_train.sample(1).path.values[0]\nprint(target_1_rand_sample_file)\n\nselected_rand_series_1 = np.load(target_1_rand_sample_file)\n\nmulti_plot(selected_rand_series_1, 'time', \"1\")\n","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:26.052551Z","iopub.execute_input":"2021-08-23T16:00:26.052931Z","iopub.status.idle":"2021-08-23T16:00:26.69194Z","shell.execute_reply.started":"2021-08-23T16:00:26.052896Z","shell.execute_reply":"2021-08-23T16:00:26.691097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extract 1 random element from the target_0_df_train\n# using random_state to ensure the reproducibility of the examples.\ntarget_0_rand_sample_file = target_0_df_train.sample(1).path.values[0]\nprint(target_0_rand_sample_file)\n\nselected_rand_series_0 = np.load(target_0_rand_sample_file)\n\nmulti_plot(selected_rand_series_0, 'time', \"0\")\n","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:26.693207Z","iopub.execute_input":"2021-08-23T16:00:26.69369Z","iopub.status.idle":"2021-08-23T16:00:27.359332Z","shell.execute_reply.started":"2021-08-23T16:00:26.693649Z","shell.execute_reply":"2021-08-23T16:00:27.358384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"multi_plot(selected_rand_series_1, \"box\", \"1\")","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:27.360651Z","iopub.execute_input":"2021-08-23T16:00:27.361113Z","iopub.status.idle":"2021-08-23T16:00:28.073817Z","shell.execute_reply.started":"2021-08-23T16:00:27.361069Z","shell.execute_reply":"2021-08-23T16:00:28.072808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"multi_plot(selected_rand_series_0, 'box', \"0\")","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:28.075292Z","iopub.execute_input":"2021-08-23T16:00:28.075706Z","iopub.status.idle":"2021-08-23T16:00:28.548184Z","shell.execute_reply.started":"2021-08-23T16:00:28.075663Z","shell.execute_reply":"2021-08-23T16:00:28.547229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Basic observation from Box Plots:**\n1. The 3 series from the 3 Interferometer have fairly similar distribution for both the class types.","metadata":{}},{"cell_type":"code","source":"multi_plot(selected_rand_series_1, 'kde', \"1\")","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:28.549708Z","iopub.execute_input":"2021-08-23T16:00:28.550058Z","iopub.status.idle":"2021-08-23T16:00:29.224894Z","shell.execute_reply.started":"2021-08-23T16:00:28.550021Z","shell.execute_reply":"2021-08-23T16:00:29.224042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"multi_plot(selected_rand_series_0, 'kde', \"0\")","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:29.226398Z","iopub.execute_input":"2021-08-23T16:00:29.226824Z","iopub.status.idle":"2021-08-23T16:00:29.859552Z","shell.execute_reply.started":"2021-08-23T16:00:29.226785Z","shell.execute_reply":"2021-08-23T16:00:29.858559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# Note on Keras Sequential model\n\nThere are two ways to build Keras models: sequential and functional.\n\n**The sequential API** allows you to create models layer-by-layer for most problems. It is limited in that it does not allow you to create models that share layers or have multiple inputs or outputs.\nIn short, you create a sequential model where you can easily add layers, and each layer can have convolution, max pooling, activation, drop-­out, and batch normalization.\n\nAlternatively, **the functional API** allows you to create models that have a lot more flexibility as you can easily define models where layers connect to more than just the previous and next layers. In fact, you can connect layers to (literally) any other layer. As a result, creating complex networks such as siamese networks and residual networks become possible.\n\nFrom the definition of Keras documentation the Sequential model is a linear stack of layers. You can create a Sequential model by passing a list of layer instances to the constructor. The common architecture of ConvNets is a sequential architecture. However, some architectures are not linear stacks. For example, siamese networks are two parallel neural networks with some shared layers.\n\nWith the Sequential Models, you need to ensure the\ninput layer has the right number of inputs. Assume that you have 3,072\ninput variables; then you need to create the first hidden layer with 512\nnodes/neurons. In the second hidden layer, you have 120 nodes/neurons.\nFinally, you have ten nodes in the output layer. For example, an image\nmaps onto ten nodes that shows the probability of being label1 (airplane),\nlabel2 (automobile), label3 (cat), ..., label10 (truck). The node of highest\nprobability is the predicted class/label.\n\n\n### Some common layer types in Keras are as follows:\n\n**Dense**: This is a fully connected layer in which all the nodes of the layer are\ndirectly connected to all the inputs and all the outputs. ANNs for classification or\nregression tasks on tabular data usually have a large percentage of their layers with\nthis type in the architecture.\n\n**Convolutional**: This layer type creates a convolutional kernel that is convolved with\nthe input layer to produce a tensor of outputs. This convolution can occur in one\nor multiple dimensions. ANNs for the classification of images usually feature one or\nmore convolutional layers in their architecture.\n\n**Pooling**: This type of layer is used to reduce the dimensionality of an input layer.\nCommon types of pooling include max pooling, in which the maximum value of\na given window is passed through to the output, or average pooling, in which\nthe average value of a window is passed through. These layers are often used\nin conjunction with a convolutional layer, and their purpose is to reduce the\ndimensions of the subsequent layers, allowing for fewer training parameters to be\nlearned with little information loss.\n\n**Recurrent**: Recurrent layers learn patterns from sequences, so each output is\ndependent on the results from the previous step. ANNs that model sequential data\nsuch as natural language or time-series data often feature one or more recurrent\nlayer types.\n\n\nThere are other layer types in Keras; however, these are the most common types when\nit comes to building models using Keras.\n\n\n## Why I would need a Custom Data Generator Function for Keras Sequential Model building\n\nThe key reason is to be able to handle large data with batching, so the RAM/CPU/GPU does not need to handle the full data at once, which will anyway not be possible for this 72GB dataset.\n\nSo basically, since our code will in most cases be multicore-friendly, so we focus on doing more complex operations (e.g. computations from source files) without worrying about data generation becoming a bottleneck in the training process.\n\n\n**DataGenerator(Sequence)** => Now, let's go through the details of how to set the Python class DataGenerator, which will be used for real-time data feeding to your Keras model. We make `DataGenerator` inherit the properties of `keras.utils.Sequence` so that we can leverage nice functionalities such as multiprocessing.\n\nWhile we have built-in Data Generator like [ImageDataGenerator](https://www.tensorflow.org/api_docs/python/tf/keras/preprocessing/image/ImageDataGenerator), we still need a plethora of custom Generator function. Because, Model training is not limited to a single type of input and target. There are times when a model is fed with multiple types of inputs at once. For example, say in a multi-modal classification problem which needs to process text and image data simultaneously. Here, obviously we cannot use ImageDataGenerator. Hence, we need a custom data generator.\n\nAccording to [Keras Documentation](https://www.tensorflow.org/api_docs/python/tf/keras/utils/Sequence) - Every Sequence must implement the `__getitem__` and the `__len__` methods. If you want to modify your dataset between epochs you may implement on_epoch_end. The method __getitem__ should return a complete batch.\n\n ### A note on `yield` function with respect to the custom data generator here for Sequence API of Keras\n\nHere I am using the Sequence API, which works a bit different than plain generators. In a generator function, you would use the `yield` keyword to perform iteration inside a while True: loop, so each time Keras calls the generator, it gets a batch of data and it automatically wraps around the end of the data.\n\nBut in a Sequence-API, there is an index parameter to the `__getitem__` function, so no iteration or `yield` is required, this is performed by Keras for you. This is made so the sequence can run in parallel using multiprocessing, which is not possible with old generator functions.\n\n---\n\n### Difference between fit() and fit_generator() in Keras\n\nIn Keras, using fit() and predict() is fine for smaller datasets which can be loaded into memory. But in practice, for most practical-use cases, almost all datasets are large and cannot be loaded into memory at once. The solution is to use fit_generator() and predict_generator() with custom data generator functions which can load data to memory during training or predicting.\n\nIn keras, fit() is much similar to sklearn's fit method, where you pass array of features as x values and target as y values. You pass your whole dataset at once in fit method. Also, use it if you can load whole data into your memory (small dataset).\n\nIn fit_generator(), you don't pass the x and y directly, instead they come from a generator. As it is written in keras documentation, generator is used when you want to avoid duplicate data when using multiprocessing. This is for practical purpose, when you have large dataset.\n\n# Final Model Building","metadata":{}},{"cell_type":"code","source":"\"\"\" First, we define the constructor to initialize the configuration of the generator.\nNote that here, we assume the path to the data is in a dataframe column.\n\n\"\"\"\n\nclass DataGenerator(Sequence):\n\n    # For this dataset the list_IDs are the value of the ids\n    # for each of the time-series file\n    # i.e. for Train data => values of column 'id' from training_labels.csv\n\n    # Also Note we have earlier defined our labels to be the below\n    # labels = pd.read_csv(root_dir + \"training_labels.csv\")\n    # and the argument \"data\" is that label here.\n    def __init__(self, path, list_IDs, data, batch_size):\n        self.path = path\n        self.list_IDs = list_IDs\n        self.data = data\n        self.batch_size = batch_size\n        self.indexes = np.arange(len(self.list_IDs))\n\n    \"\"\" __len__ essentially returns the number of steps in an epoch, using the samples and the batch size.\n        Each call requests a batch index between 0 and the total number of batches, where the latter is specified in the __len__ method.\n        A common practice is to set this value to (samples / batch size)\n        so that the model sees the training samples at most once per epoch.\n        Now, when the batch corresponding to a given index is called, the generator executes the __getitem__ method to generate it.\n    \"\"\"\n\n    def __len__(self):\n        len_ = int(len(self.list_IDs)/self.batch_size)\n        if len_ * self.batch_size < len(self.list_IDs):\n            len_ += 1\n        return len_\n\n    \"\"\"  __getitem__ method is called with the batch number as an argument to obtain a given batch of data.\n\n    \"\"\"\n    def __getitem__(self, index):\n        # get the range to to feed to keras for each epoch\n        # incrementing by +1 the bath_size\n        indexes = self.indexes[index * self.batch_size : (index + 1) * self.batch_size]\n        list_IDs_temp = [self.list_IDs[k] for k in indexes]\n        X, y = self.__data_generation(list_IDs_temp)\n        return X, y\n\n    \"\"\" And finally the core method which will actually produce batches of data. This private method __data_generation \"\"\"\n\n    def __data_generation(self, list_IDs_temp):\n        # We have 5,60,000 files, each with dimension of 3 * 4096\n        X = np.zeros((self.batch_size, 3, 4096))\n        y = np.zeros((self.batch_size, 1))\n        for i, ID in enumerate(list_IDs_temp):\n            id_ = self.data.loc[ID, \"id\"]\n            file = id_ + \".npy\"  # build the file name\n            path_in = \"/\".join([self.path, id_[0], id_[1], id_[2]]) + \"/\"\n            # there are three nesting labels inside train/ or test/\n            data_array = np.load(path_in + file)            \n            data_array = (data_array - data_array.mean())/data_array.std()\n            X[i, ] = data_array\n            y[i, ] = self.data.loc[ID, 'target']\n        # print(X)\n        return X, y","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:29.861068Z","iopub.execute_input":"2021-08-23T16:00:29.861419Z","iopub.status.idle":"2021-08-23T16:00:29.873888Z","shell.execute_reply.started":"2021-08-23T16:00:29.861381Z","shell.execute_reply":"2021-08-23T16:00:29.872796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission = pd.read_csv(root_dir +  'sample_submission.csv')\n# print(len(train_labels)) # 5,60,000\n# print(len(sample_submission)) # 2,26,000\ntrain_ids = train_labels['id'].values\n# train_ids # ['00000e74ad', '00001f4945', '0000661522' ... ]\ny = train_labels['target'].values\ntest_ids = sample_submission['id'].values","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:29.875455Z","iopub.execute_input":"2021-08-23T16:00:29.876217Z","iopub.status.idle":"2021-08-23T16:00:30.049385Z","shell.execute_reply.started":"2021-08-23T16:00:29.876172Z","shell.execute_reply":"2021-08-23T16:00:30.04852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_labels = pd.read_csv(root_dir + \"training_labels.csv\", nrows=1000)\n# Now I shall genereate train indices, validation indices and test indices\n# Which are just the values from the 0-based indices\ntrain_indices, validation_indices = train_test_split(list(train_labels.index), test_size=0.33, random_state=2021)\n# print(len(train_indices))\nprint(len(validation_indices))\ntest_indices = list(sample_submission.index)\n# test_indices","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:30.051256Z","iopub.execute_input":"2021-08-23T16:00:30.051882Z","iopub.status.idle":"2021-08-23T16:00:30.31196Z","shell.execute_reply.started":"2021-08-23T16:00:30.051843Z","shell.execute_reply":"2021-08-23T16:00:30.310896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ntrain_generator_for_seq_model = DataGenerator( root_dir +  'train/', train_indices, train_labels, 64)\n# print(train_generator_for_seq_model)\nvalidation_generator_for_seq_model = DataGenerator( root_dir + 'train/', validation_indices, train_labels, 64)\ntest_generator_for_seq_model = DataGenerator( root_dir + 'test/', test_indices, sample_submission, 64)","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:30.313497Z","iopub.execute_input":"2021-08-23T16:00:30.313872Z","iopub.status.idle":"2021-08-23T16:00:30.319903Z","shell.execute_reply.started":"2021-08-23T16:00:30.313833Z","shell.execute_reply":"2021-08-23T16:00:30.318881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_keras_seq = Sequential()\nmodel_keras_seq.add(Conv1D(64, input_shape=(3, 4096), kernel_size=3, activation='relu'))\nmodel_keras_seq.add(BatchNormalization())\nmodel_keras_seq.add(Flatten())\nmodel_keras_seq.add(Dense(64, activation='relu'))\nmodel_keras_seq.add(Dense(1, activation='sigmoid'))\n\nmodel_keras_seq.compile(optimizer= Adam(lr=2e-4), loss='binary_crossentropy', metrics=['acc'])\nmodel_keras_seq.summary()","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:30.323816Z","iopub.execute_input":"2021-08-23T16:00:30.324205Z","iopub.status.idle":"2021-08-23T16:00:36.211881Z","shell.execute_reply.started":"2021-08-23T16:00:30.324174Z","shell.execute_reply":"2021-08-23T16:00:36.210332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history = model_keras_seq.fit_generator(generator=train_generator_for_seq_model, validation_data=validation_generator_for_seq_model, epochs = 1, workers=-1)\n# 1 Epoch took around 2 and half hours\n\npredicted_test_seq_keras = model_keras_seq.predict_generator(test_generator_for_seq_model, verbose=1)\n\nsample_submission['target'] = predicted_test_seq_keras[:len(sample_submission)]\n\nsample_submission.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-08-23T16:00:36.213327Z","iopub.execute_input":"2021-08-23T16:00:36.213708Z"},"trusted":true},"execution_count":null,"outputs":[]}]}