{"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":"code","source":"fast_sub = True # set this to False to generate the whole dataset\ntrain = True # set this to True to generate the train set\ntest = True # set this to True to generate the test set","metadata":{"execution":{"iopub.status.busy":"2021-08-18T22:21:41.087611Z","iopub.execute_input":"2021-08-18T22:21:41.088205Z","iopub.status.idle":"2021-08-18T22:21:41.102764Z","shell.execute_reply.started":"2021-08-18T22:21:41.08816Z","shell.execute_reply":"2021-08-18T22:21:41.10151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Constant-Q Transform Preprocessing for G2Net Gravitational Wave Detection\n\nIf you find this notebook and dataset helpful, please, consider upvoting it. I took inspiration from the notebook of Geir Drange: https://www.kaggle.com/mistag/data-preprocessing-with-gwpy, so if give him an upvote too.\n\nIf you would like to use the pre-baked dataset I made, here are the links, but if you are here to learn how I did it, you are in the right place =)\n\n## Link to the dataset\n* **TRAIN DATASET PT1** https://www.kaggle.com/coldfir3/g2net-cqt-dataset-pt1-jpgrgb\n* **TRAIN DATASET PT2** https://www.kaggle.com/coldfir3/g2net-cqt-dataset-pt2-jpgrgb\n* **TRAIN DATASET PT3**: https://www.kaggle.com/coldfir3/g2net-cqt-dataset-pt3-jpgrgb\n* **TEST DATASET**: https://www.kaggle.com/coldfir3/g2net-cqt-dataset-test-jpgrgb\n\n## Explanation\nThis notebook inputs raw waveform data and outputs RGB images. For this, I employed the following pipeline:\n\n1. Read the .npy file and store it as a NumPy array.\n1. Convert each channel to a gwpy TimeSeries\n1. Whiten the signal and apply a Tukey window.\n1. (optional) Apply a bandpass filter. It is optional as `q_transform` will already do this. I also noticed that by not using this step the borders of the image are sharper.\n1. Apply the q-transform\n1. Convert to RGB:\n    1. Stack each channel on `dim = -1`\n    1. Clip the values to q_max and scale it from `[0, q_max]` to `[0, 255]`\n    1. Convert the array to unsigned 8 bits\n    1. Create the image and rotate it 90º, with frequency on the vertical axis and time on the horizontal axis.\n1. Save it as a .jpeg file\n\nThe main difference from Geir's notebook is the way I built the RGB images. Geir uses `MinMaxScaler` to normalize each channel to [0 - 1] before assembling the images in the original code. I think there are a couple of problems with this approach. 1) when using min-max coupled with the subsequent 8-bits discretization, outliers could wash out all the information contained in the image (i.e. the image would be mostly black). Secondly, if the signal is more potent in detector A when compared to detector B, you will lose this information by normalizing each channel independently.\n\nIt is worth noticing, though, that the default value I chose for the max_q was utterly arbitrary, and different values could lead to different results. The other defaults were ported from Geir's notebook, and I encourage you to experiment with different values and be kind to share your findings in the comment section below.","metadata":{}},{"cell_type":"markdown","source":"## Installing and loading the dependencies","metadata":{}},{"cell_type":"code","source":"%%capture\n!python -m pip install gwpy\n!pip install astropy==4.2.1","metadata":{"execution":{"iopub.status.busy":"2021-08-18T22:21:41.104436Z","iopub.execute_input":"2021-08-18T22:21:41.104728Z","iopub.status.idle":"2021-08-18T22:22:10.277623Z","shell.execute_reply.started":"2021-08-18T22:21:41.104701Z","shell.execute_reply":"2021-08-18T22:22:10.276237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from gwpy.timeseries import TimeSeries\nfrom gwpy.plot import Plot\nimport numpy as np\nfrom scipy import signal\nfrom PIL import Image\n\nimport os\nimport pandas as pd\nimport numpy as np\nfrom tqdm.auto import tqdm\nfrom joblib import Parallel, delayed\nfrom pathlib import Path","metadata":{"execution":{"iopub.status.busy":"2021-08-18T23:38:49.602745Z","iopub.execute_input":"2021-08-18T23:38:49.603182Z","iopub.status.idle":"2021-08-18T23:38:49.608451Z","shell.execute_reply.started":"2021-08-18T23:38:49.603148Z","shell.execute_reply":"2021-08-18T23:38:49.607557Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Main processing function","metadata":{}},{"cell_type":"code","source":"def sig2rgb(fname, whiten = True, window=0.2, bandpass=False, f_range = (30,400), q_range = (16,32), q_max = 10):\n    \n    # Load the file \n    data = np.load(fname)\n    # Split each chanel and convert to TimeSeries\n    data = map(lambda x: TimeSeries(x, sample_rate=2048), data)\n    # Whiten the signal and apply a tukey window\n    data = map(lambda x: x.whiten(window=(\"tukey\", window)), data)\n    # (optional) bandpass filter\n    if bandpass:\n        data = map(lambda x: x.bandpass(*f_range), data)\n    # Q-transform\n    data = map(lambda x: x.q_transform(qrange=q_range, frange=f_range, logf=True, whiten=False), data)\n    # Convert to RGB image\n    img = np.stack(list(data), axis = -1)\n    img = np.clip(img, 0, q_max)/q_max * 255\n    img = img.astype(np.uint8)\n    img = Image.fromarray(img).rotate(90, expand=1)\n    return img","metadata":{"execution":{"iopub.status.busy":"2021-08-18T23:01:40.401731Z","iopub.execute_input":"2021-08-18T23:01:40.402124Z","iopub.status.idle":"2021-08-18T23:01:40.411374Z","shell.execute_reply.started":"2021-08-18T23:01:40.402091Z","shell.execute_reply":"2021-08-18T23:01:40.410251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sig2rgb('../input/g2net-gravitational-wave-detection/train/0/0/0/000a5b6e5c.npy')","metadata":{"execution":{"iopub.status.busy":"2021-08-18T23:01:42.346736Z","iopub.execute_input":"2021-08-18T23:01:42.347602Z","iopub.status.idle":"2021-08-18T23:01:42.875964Z","shell.execute_reply.started":"2021-08-18T23:01:42.347554Z","shell.execute_reply":"2021-08-18T23:01:42.874892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Generating the dataset and ziping the files","metadata":{}},{"cell_type":"code","source":"def save_img(x, folder_out, **kwargs):\n    fname = Path('../input/g2net-gravitational-wave-detection/' + folder_out + '/' + '/'.join([x[0], x[1], x[2], x]) + '.npy')\n    file_out = folder_out + '/' + fname.with_suffix('.jpg').name\n    x = sig2rgb(fname, **kwargs)\n    x.save(file_out)","metadata":{"execution":{"iopub.status.busy":"2021-08-18T23:19:27.113751Z","iopub.execute_input":"2021-08-18T23:19:27.114377Z","iopub.status.idle":"2021-08-18T23:19:27.121255Z","shell.execute_reply.started":"2021-08-18T23:19:27.114326Z","shell.execute_reply":"2021-08-18T23:19:27.120393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from zipfile import ZipFile\nimport shutil\nimport os\ndef zip_folder(folder, rm_original = True):\n    # iterate over all the files in directory\n    for folderName, subfolders, filenames in os.walk(folder):\n        # create a ZipFile object\n        with ZipFile(folderName.split('/')[-1] + '.zip', 'w') as zipObj:\n            for filename in filenames:\n                # create complete filepath of file in directory\n                filePath = os.path.join(folderName, filename)\n                # add file to zip\n                zipObj.write(filePath, os.path.basename(filePath))\n                # delete the file to open space\n                if rm_original:\n                    os.remove(filePath)\n    if rm_original:\n        shutil.rmtree(folder)","metadata":{"execution":{"iopub.status.busy":"2021-08-18T23:40:20.026331Z","iopub.execute_input":"2021-08-18T23:40:20.027009Z","iopub.status.idle":"2021-08-18T23:40:20.034287Z","shell.execute_reply.started":"2021-08-18T23:40:20.026972Z","shell.execute_reply":"2021-08-18T23:40:20.033311Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test data","metadata":{}},{"cell_type":"code","source":"test_df = pd.read_csv('../input/g2net-gravitational-wave-detection/sample_submission.csv')\nif fast_sub: test_ids = test_df['id'][:100]\nelse: test_ids = test_df['id']\nif test:\n    os.makedirs('test', exist_ok = True)\n    o = Parallel(n_jobs=4)(delayed(save_img)(x, 'test') for x in tqdm(test_ids))\n    zip_folder('test')","metadata":{"execution":{"iopub.status.busy":"2021-08-18T23:39:16.573711Z","iopub.execute_input":"2021-08-18T23:39:16.574232Z","iopub.status.idle":"2021-08-18T23:39:32.493948Z","shell.execute_reply.started":"2021-08-18T23:39:16.574199Z","shell.execute_reply":"2021-08-18T23:39:32.492931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Train data","metadata":{}},{"cell_type":"code","source":"test_df = pd.read_csv('../input/g2net-gravitational-wave-detection/training_labels.csv')\nif fast_sub: test_ids = test_df['id'][:100]\nelse: test_ids = test_df['id']\nif train:\n    os.makedirs('train', exist_ok = True)\n    o = Parallel(n_jobs=4)(delayed(save_img)(x, 'train') for x in tqdm(test_ids))\n    zip_folder('train')","metadata":{"execution":{"iopub.status.busy":"2021-08-18T23:40:59.612815Z","iopub.execute_input":"2021-08-18T23:40:59.613197Z","iopub.status.idle":"2021-08-18T23:41:12.480633Z","shell.execute_reply.started":"2021-08-18T23:40:59.613166Z","shell.execute_reply":"2021-08-18T23:41:12.47942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}