{"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":"# This Notebook\n\nThis notebook uses [Bayley et al.'s (2020)](https://arxiv.org/pdf/2007.08207.pdf) SOAP algorithm to extract features from spectrograms for the G2Net continuous gravitational waves task. The extracted features are intended to be used with [our other notebook](https://www.kaggle.com/code/hnsyprst/publish-denoised-data-soap-cnn), **which includes a visualisation of the features extracted by this notebook**. For speed, the code for generating these visualisations has note been replicated here; to see this code, please view our other notebook.\n\nThe input data to this notebook is a set of spectrograms gathered from Chirita's multiple synthetic datasets and the original G2Net data. These spectrograms have been passed through [laeyoung's stem](https://www.kaggle.com/code/laeyoung/g2net-large-kernel-inference) to improve their signal-to-noise ratio. The notebook used to gather and process the input data can be found [here](https://www.kaggle.com/code/hnsyprst/publish-consolidatedg2netdataset?scriptVersionId=115401063).\n\nSOAP is based on the Viterbi algorithm. The Viterbi algorithm finds the most likely sequence of states through a series of states with given transition probabilities. In SOAP's context, this means finding the most likely path through a given spectrogram (the path that will give the highest sum of FFT power). If a signal is present in a given spectrogram, this path should correspond with that signal. These paths, which Bayley et al. call 'Viterbi tracks', have been extracted for each of the noise reduced spectrograms described above. In addition to the Viterbi tracks, this notebook extracts 'Viterbi maps', which encode the probability that the signal (the Viterbi track) is in each frequency bin at each time across a given spectrogram.\n\nBayley, J., Messenger, C. and Woan, G. (2020) ‘A robust machine learning algorithm to search for continuous gravitational waves’, *Physical Review D*, 102(8), p. 083024. Available at: https://doi.org/10.1103/PhysRevD.102.083024.","metadata":{}},{"cell_type":"markdown","source":"# Example Output\n(with Original Data for Comparison)","metadata":{}},{"cell_type":"code","source":"!pip install timm\n\n!wget https://github.com/lscsoft/lalsuite/raw/master/lalpulsar/lib/earth00-19-DE405.dat.gz -P /kaggle/working/\n!wget https://github.com/lscsoft/lalsuite/raw/master/lalpulsar/lib/sun00-19-DE405.dat.gz -P /kaggle/working/\n\n!env LAL_DATA_PATH=/kaggle/working/\n\n!pip install soapcw","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-04T11:31:14.23027Z","iopub.execute_input":"2023-01-04T11:31:14.230797Z","iopub.status.idle":"2023-01-04T11:33:54.571804Z","shell.execute_reply.started":"2023-01-04T11:31:14.230755Z","shell.execute_reply":"2023-01-04T11:33:54.569883Z"},"jupyter":{"source_hidden":true,"outputs_hidden":true},"collapsed":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import shutil\nimport os\nimport gc, glob, os\nfrom concurrent.futures import ProcessPoolExecutor\n\nimport h5py\nimport numpy as np\nimport pandas as pd\nimport torch\nimport torch.nn as nn\nfrom scipy.stats import norm\nfrom timm import create_model\nfrom tqdm.notebook import tqdm\n\nimport matplotlib.pyplot as plt\nimport soapcw as soap\nfrom soapcw import cw\n\ntracks_save_dir = '/kaggle/working/vit_tracks/train'\nmaps_save_dir = '/kaggle/working/vit_maps/train'\n\nos.makedirs(tracks_save_dir + '/signals', exist_ok=True)\nos.makedirs(tracks_save_dir + '/noises', exist_ok=True)\nos.makedirs(maps_save_dir + '/signals', exist_ok=True)\nos.makedirs(maps_save_dir + '/noises', exist_ok=True)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-04T11:33:54.574602Z","iopub.execute_input":"2023-01-04T11:33:54.575065Z","iopub.status.idle":"2023-01-04T11:33:57.441302Z","shell.execute_reply.started":"2023-01-04T11:33:54.575023Z","shell.execute_reply":"2023-01-04T11:33:57.439587Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_track(file):\n    file = np.squeeze(file)\n    return file\n\ndef display_data(H1, L1, vitmap, track):\n    fig = plt.figure(figsize=(16,16), constrained_layout=True)\n    \n    subfigs = fig.subfigures(nrows=2, ncols=1)\n    subfigs[0].suptitle('Spectrograms before Feature Extraction', fontsize=16)\n    axs_0 = subfigs[0].subplots(nrows=2, ncols=1)\n    axs_0[0].imshow(H1.T,aspect=\"auto\",origin=\"lower\",cmap=\"YlGnBu\")\n    axs_0[0].set_title('H1 Spectrogram')\n    \n    axs_0[1].imshow(L1.T,aspect=\"auto\",origin=\"lower\",cmap=\"YlGnBu\")\n    axs_0[1].set_title('L1 Spectrogram')\n    \n    subfigs[1].suptitle('Features Extracted by this Notebook', fontsize=16)\n    axs_1 = subfigs[1].subplots(nrows=2, ncols=1)\n    axs_1[0].imshow(vitmap.T,aspect=\"auto\",origin=\"lower\",cmap=\"YlGnBu\")\n    axs_1[0].set_title('Viterbi Map')\n    \n    axs_1[1].scatter(np.arange(H1.shape[0]), np.argmax(track, axis=1), color=\"red\", s=3)\n    axs_1[1].set_ylim([0, H1.shape[1]])\n    axs_1[1].set_xlim([0, H1.shape[0]])\n    axs_1[1].set_title('Viterbi Track')\n    \nexample_H1 = np.load('/kaggle/input/g2net-examples/example_original_H1.npy')\nexample_L1 = np.load('/kaggle/input/g2net-examples/example_original_L1.npy')\nexample_vitmap = np.load('/kaggle/input/g2net-examples/example_vitmap.npy')\nexample_track = load_track(np.load('/kaggle/input/g2net-examples/example_track.npy'))\n\ndisplay_data(example_H1, example_L1, example_vitmap, example_track)","metadata":{"execution":{"iopub.status.busy":"2023-01-04T11:36:42.443976Z","iopub.execute_input":"2023-01-04T11:36:42.444532Z","iopub.status.idle":"2023-01-04T11:36:44.325308Z","shell.execute_reply.started":"2023-01-04T11:36:42.444492Z","shell.execute_reply":"2023-01-04T11:36:44.323524Z"},"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature Extraction","metadata":{}},{"cell_type":"code","source":"# change the number of elements to reduce generation time (500 -> 100)\npowers = np.linspace(1,200,500)\nv1d = soap.line_aware_stat.gen_lookup_python.LineAwareStatistic(powers,ndet=1,signal_prior_width=4.0,line_prior_width=5.0,noise_line_model_ratio=0.038)\nv1d.save_lookup(\"/kaggle/working/\")","metadata":{"execution":{"iopub.status.busy":"2023-01-02T21:44:16.778934Z","iopub.execute_input":"2023-01-02T21:44:16.779826Z","iopub.status.idle":"2023-01-02T21:44:18.230592Z","shell.execute_reply.started":"2023-01-02T21:44:16.779789Z","shell.execute_reply":"2023-01-02T21:44:18.229742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"di = '/kaggle/input/g2net-detecting-continuous-gravitational-waves'\nTRAIN_DF = pd.read_csv(di + '/train_labels.csv', index_col=False)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T21:45:42.962676Z","iopub.execute_input":"2023-01-02T21:45:42.963142Z","iopub.status.idle":"2023-01-02T21:45:42.974589Z","shell.execute_reply.started":"2023-01-02T21:45:42.963098Z","shell.execute_reply":"2023-01-02T21:45:42.973535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_id_and_spect(filepath):\n    # Load and preprocess denoised STFT\n    denoised = np.load(filepath)\n    \n    a = denoised[:, :720]\n\n    p = a.real**2 + a.imag**2  # power\n    p /= np.mean(p)  # normalize\n    p = np.sum(p.reshape(360, 90, 8), axis=2)\n\n    spect = p\n    spect = np.rot90(spect, 3)\n    \n    file_id = os.path.basename(filepath).replace('.npy', '')\n    \n    gc.collect()\n    return file_id, spect \n\ndef load_label(file_id):\n    return TRAIN_DF.loc[TRAIN_DF['id'] == file_id].target.item()\n\ndef train_dataload(filepath):\n    file_id, spect = load_id_and_spect(filepath)\n    #y = load_label(file_id)\n    \n    return file_id, spect\n\ndef test_dataload(filepath):\n    file_id, spect = load_id_and_spect(filepath)\n    \n    return file_id, spect","metadata":{"execution":{"iopub.status.busy":"2023-01-02T21:45:45.857435Z","iopub.execute_input":"2023-01-02T21:45:45.85786Z","iopub.status.idle":"2023-01-02T21:45:45.868088Z","shell.execute_reply.started":"2023-01-02T21:45:45.857828Z","shell.execute_reply":"2023-01-02T21:45:45.866718Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def vector_to_onehot(vector, vocab_length):\n    tensor = np.zeros((len(vector), 1, vocab_length))\n    for i, value in enumerate(vector):\n        tensor[i][0][value] = 1\n    return tensor\n\ndef make_vit(spect, vocab_length):\n    tr_1 = soap.tools.transition_matrix(1.0) \n    tracks_ng_ls = soap.single_detector(tr_1, spect, lookup_file='/kaggle/working/log_signoiseline_1det_96degfree_4.0_5.0_0.038.pkl')\n    \n    vit_track = vector_to_onehot(tracks_ng_ls.vit_track, vocab_length)\n    return vit_track, tracks_ng_ls.vitmap","metadata":{"execution":{"iopub.status.busy":"2023-01-02T21:48:42.298727Z","iopub.execute_input":"2023-01-02T21:48:42.299129Z","iopub.status.idle":"2023-01-02T21:48:42.306904Z","shell.execute_reply.started":"2023-01-02T21:48:42.299096Z","shell.execute_reply":"2023-01-02T21:48:42.305531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file_path = glob.glob(os.path.join(\"/kaggle/input/g2net-consolidated-reduced-noise/signals\", \"*.npy\"))\nvocab_length = 360\n\nwith ProcessPoolExecutor(4) as pool:\n    data = []\n    for file_id, spect in pool.map(train_dataload, sorted(file_path)):\n        vit_track, vit_map = make_vit(spect, vocab_length)\n        np.save('%s/signals/%s.npy' % (tracks_save_dir, file_id), vit_track)\n        np.save('%s/signals/%s.npy' % (maps_save_dir, file_id), vit_map)        \n        del(vit_track)\n        del(vit_map)\n        gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T21:48:43.480097Z","iopub.execute_input":"2023-01-02T21:48:43.480531Z","iopub.status.idle":"2023-01-02T21:54:56.101356Z","shell.execute_reply.started":"2023-01-02T21:48:43.480493Z","shell.execute_reply":"2023-01-02T21:54:56.099974Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file_path = glob.glob(os.path.join(\"/kaggle/input/g2net-consolidated-reduced-noise/noises\", \"*.npy\"))\nvocab_length = 360\n\nwith ProcessPoolExecutor(4) as pool:\n    data = []\n    for file_id, spect in pool.map(train_dataload, sorted(file_path)):\n        vit_track, vit_map = make_vit(spect, vocab_length)\n        np.save('%s/noises/%s.npy' % (tracks_save_dir, file_id), vit_track)\n        np.save('%s/noises/%s.npy' % (maps_save_dir, file_id), vit_map)        \n        del(vit_track)\n        del(vit_map)\n        gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T21:54:56.103487Z","iopub.execute_input":"2023-01-02T21:54:56.103908Z","iopub.status.idle":"2023-01-02T22:03:26.109577Z","shell.execute_reply.started":"2023-01-02T21:54:56.10387Z","shell.execute_reply":"2023-01-02T22:03:26.108173Z"},"trusted":true},"execution_count":null,"outputs":[]}]}