{"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":"**What are you trying to do in this notebook?**\n\nMy goal for this competition is to find continuous gravitational-wave signals. I will develop a model sensitive enough to detect weak yet long-lasting signals emitted by rapidly-spinning neutron stars within noisy data.\n\nMy work will help scientists detect something new: a second class of gravitational waves! The first gravitational wave discoveries earned a Nobel Prize. Further study of these waves may enable scientists to learn about the structure of the most extreme stars in our universe.\n\nWhen scientists detected the first class of gravitational waves in 2015, they expected the discoveries to continue. There are four classes, yet at present only signals from merging black holes and neutron stars have been detected. Among those remaining are continuous gravitational-wave signals. These are weak yet long-lasting signals emitted by rapidly-spinning neutron stars. Imagine the mass of our Sun but condensed into a ball the size of a city and spinning over 1,000 times a second. The extreme compactness of these stars, composed of the densest material in the universe, could allow continuous waves to be emitted and then detected on Earth. There are potentially many continuous signals from neutron stars in our own galaxy and the current challenge for scientists is to make the first detection, and hopefully data science can help with this mission.\n\n**Why are you trying it?**\n\nIn this challenge I'll enable scientists to improve their sensitivity, leading to new discoveries in the field. As a result, scientists could learn more about the structure of the most extreme stars in our universe with the help of G2Net.\n\nG2Net aims to create a broad network of scientists.\n\nFrom four different areas of expertise, namely GW physics, Geophysics, Computing Science and Robotics, these scientists have agreed on a common goal of tackling challenges in data analysis and noise characterization for GW detectors.\n\nBy helping G2Net in this challenge I'll enable scientists to improve their sensitivity, leading to new discoveries in the field. As a result, scientists could learn more about the structure of the most extreme stars in our universe.\n\nIn this competition, I'm provided with a training set containing time-frequency data from two gravitational-wave interferometers (LIGO Hanford & LIGO Livingston). Each data sample contains either real or simulated noise and possibly a simulated continuous gravitational-wave signal (CW). The task is to identify when a signal is present in the data (target=1).\n\n","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:24:53.907468Z","iopub.execute_input":"2022-11-21T12:24:53.90837Z","iopub.status.idle":"2022-11-21T12:24:56.069496Z","shell.execute_reply.started":"2022-11-21T12:24:53.908273Z","shell.execute_reply":"2022-11-21T12:24:56.065894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" COLAB = False\n\nif COLAB == True:\n    from google.colab import drive\n    drive.mount('/content/drive')\n    %cd '/content/drive/MyDrive/Colab Notebooks/kaggle/G2Net2022/code'","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:24:56.07115Z","iopub.execute_input":"2022-11-21T12:24:56.071484Z","iopub.status.idle":"2022-11-21T12:24:56.084757Z","shell.execute_reply.started":"2022-11-21T12:24:56.071447Z","shell.execute_reply":"2022-11-21T12:24:56.083864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"! pip3 install timm -q","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:24:56.08642Z","iopub.execute_input":"2022-11-21T12:24:56.087182Z","iopub.status.idle":"2022-11-21T12:25:08.643103Z","shell.execute_reply.started":"2022-11-21T12:24:56.087145Z","shell.execute_reply":"2022-11-21T12:25:08.64198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport time\nimport h5py\nimport timm\nimport torch\nimport torch.nn as nn\nimport torchaudio\nimport torchvision.transforms as TF\n\n\nfrom tqdm.auto import tqdm\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import roc_auc_score\nfrom timm.scheduler import CosineLRScheduler\n\ndevice = torch.device('cuda')\ncriterion = nn.BCEWithLogitsLoss()\n\ndi = '../input/g2net-detecting-continuous-gravitational-waves'\ndf = pd.read_csv(di + '/train_labels.csv')\ndf = df[df.target >= 0]  ","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:08.647599Z","iopub.execute_input":"2022-11-21T12:25:08.647936Z","iopub.status.idle":"2022-11-21T12:25:11.826837Z","shell.execute_reply.started":"2022-11-21T12:25:08.647905Z","shell.execute_reply":"2022-11-21T12:25:11.825717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"transforms_time_mask = nn.Sequential(\n                torchaudio.transforms.TimeMasking(time_mask_param=10),\n            )\n\ntransforms_freq_mask = nn.Sequential(\n                torchaudio.transforms.FrequencyMasking(freq_mask_param=10),\n            )\n\nflip_rate = 0.0 \nfre_shift_rate = 0.0 \n\ntime_mask_num = 0 \nfreq_mask_num = 0 ","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:11.829117Z","iopub.execute_input":"2022-11-21T12:25:11.829802Z","iopub.status.idle":"2022-11-21T12:25:11.836859Z","shell.execute_reply.started":"2022-11-21T12:25:11.82976Z","shell.execute_reply":"2022-11-21T12:25:11.835584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Dataset(torch.utils.data.Dataset):\n    \"\"\"\n    dataset = Dataset(data_type, df)\n\n    img, y = dataset[i]\n      img (np.float32): 2 x 360 x 128\n      y (np.float32): label 0 or 1\n    \"\"\"\n    def __init__(self, data_type, df, tfms=False):\n        self.data_type = data_type\n        self.df = df\n        self.tfms = tfms\n\n    def __len__(self):\n        return len(self.df)\n\n    def __getitem__(self, i):\n        \"\"\"\n        i (int): get ith data\n        \"\"\"\n        r = self.df.iloc[i]\n        y = np.float32(r.target)\n        file_id = r.id\n\n        img = np.empty((2, 360, 128), dtype=np.float32)\n\n        filename = '%s/%s/%s.hdf5' % (di, self.data_type, file_id)\n        with h5py.File(filename, 'r') as f:\n            g = f[file_id]\n\n            for ch, s in enumerate(['H1', 'L1']):\n                a = g[s]['SFTs'][:, :4096] * 1e22  \n\n                p = a.real**2 + a.imag**2  \n                p /= np.mean(p)  \n                p = np.mean(p.reshape(360, 128, 32), axis=2)  \n                img[ch] = p\n\n        if self.tfms:\n            if np.random.rand() <= flip_rate: \n                img = np.flip(img, axis=1).copy()\n            if np.random.rand() <= flip_rate: \n                img = np.flip(img, axis=2).copy()\n            if np.random.rand() <= fre_shift_rate: \n                img = np.roll(img, np.random.randint(low=0, high=img.shape[1]), axis=1)\n            \n            img = torch.from_numpy(img)\n\n            for _ in range(time_mask_num): \n                img = transforms_time_mask(img)\n            for _ in range(freq_mask_num): \n                img = transforms_freq_mask(img)\n        \n        else:\n            img = torch.from_numpy(img)\n                \n        return img, y","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:11.838332Z","iopub.execute_input":"2022-11-21T12:25:11.838681Z","iopub.status.idle":"2022-11-21T12:25:11.852867Z","shell.execute_reply.started":"2022-11-21T12:25:11.838646Z","shell.execute_reply":"2022-11-21T12:25:11.851931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dataset = Dataset('train', df, tfms=False)\nimg, y = dataset[10]\n\n\nplt.figure(figsize=(8, 3))\nplt.title('Spectrogram')\nplt.xlabel('time')\nplt.ylabel('frequency')\nplt.imshow(img[0, 0:360])\nplt.colorbar()\nplt.show()\n\n\nflip_rate = 1.0 \ndataset = Dataset('train', df, tfms=True)\nimg, y = dataset[10]\n\nplt.figure(figsize=(8, 3))\nplt.title('Spectrogram')\nplt.xlabel('time')\nplt.ylabel('frequency')\nplt.imshow(img[0, 0:360])\nplt.colorbar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:11.854473Z","iopub.execute_input":"2022-11-21T12:25:11.854828Z","iopub.status.idle":"2022-11-21T12:25:13.202798Z","shell.execute_reply.started":"2022-11-21T12:25:11.854782Z","shell.execute_reply":"2022-11-21T12:25:13.201772Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dataset = Dataset('train', df, tfms=False)\nimg, y = dataset[10]\n\n\nplt.figure(figsize=(8, 3))\nplt.title('Spectrogram')\nplt.xlabel('time')\nplt.ylabel('frequency')\nplt.imshow(img[0, 0:360])\nplt.colorbar()\nplt.show()\n\n\nflip_rate = 0.0  \nfre_shift_rate = 1.0 \n\ndataset = Dataset('train', df, tfms=True)\nimg, y = dataset[10]\n\nplt.figure(figsize=(8, 3))\nplt.title('Spectrogram')\nplt.xlabel('time')\nplt.ylabel('frequency')\nplt.imshow(img[0, 0:360])\nplt.colorbar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:13.204385Z","iopub.execute_input":"2022-11-21T12:25:13.204721Z","iopub.status.idle":"2022-11-21T12:25:14.189721Z","shell.execute_reply.started":"2022-11-21T12:25:13.204687Z","shell.execute_reply":"2022-11-21T12:25:14.188811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dataset = Dataset('train', df, tfms=False)\nimg, y = dataset[10]\n\n\nplt.figure(figsize=(8, 3))\nplt.title('Spectrogram')\nplt.xlabel('time')\nplt.ylabel('frequency')\nplt.imshow(img[0, 0:360])\nplt.colorbar()\nplt.show()\n\n\nflip_rate = 0.0  \nfre_shift_rate = 0.0 \ntime_mask_num = 3 \n\ndataset = Dataset('train', df, tfms=True)\nimg, y = dataset[10]\n\nplt.figure(figsize=(8, 3))\nplt.title('Spectrogram')\nplt.xlabel('time')\nplt.ylabel('frequency')\nplt.imshow(img[0, 0:360])\nplt.colorbar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:14.191106Z","iopub.execute_input":"2022-11-21T12:25:14.19174Z","iopub.status.idle":"2022-11-21T12:25:15.202451Z","shell.execute_reply.started":"2022-11-21T12:25:14.191704Z","shell.execute_reply":"2022-11-21T12:25:15.201342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dataset = Dataset('train', df, tfms=False)\nimg, y = dataset[10]\n\n\nplt.figure(figsize=(8, 3))\nplt.title('Spectrogram')\nplt.xlabel('time')\nplt.ylabel('frequency')\nplt.imshow(img[0, 0:360])\nplt.colorbar()\nplt.show()\n\n\nflip_rate = 0.0 \nfre_shift_rate = 0.0 \ntime_mask_num = 0 \nfreq_mask_num = 3 \n\ndataset = Dataset('train', df, tfms=True)\nimg, y = dataset[10]\n\nplt.figure(figsize=(8, 3))\nplt.title('Spectrogram')\nplt.xlabel('time')\nplt.ylabel('frequency')\nplt.imshow(img[0, 0:360])\nplt.colorbar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:15.205366Z","iopub.execute_input":"2022-11-21T12:25:15.205647Z","iopub.status.idle":"2022-11-21T12:25:16.343742Z","shell.execute_reply.started":"2022-11-21T12:25:15.205621Z","shell.execute_reply":"2022-11-21T12:25:16.342746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Model(nn.Module):\n    def __init__(self, name, *, pretrained=False):\n        \"\"\"\n        name (str): timm model name, e.g. tf_efficientnet_b2_ns\n        \"\"\"\n        super().__init__()\n\n        # Use timm\n        model = timm.create_model(name, pretrained=pretrained, in_chans=2)\n\n        clsf = model.default_cfg['classifier']\n        n_features = model._modules[clsf].in_features\n        model._modules[clsf] = nn.Identity()\n\n        self.fc = nn.Linear(n_features, 1)\n        self.model = model\n\n    def forward(self, x):\n        x = self.model(x)\n        x = self.fc(x)\n        return x","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:16.345227Z","iopub.execute_input":"2022-11-21T12:25:16.345565Z","iopub.status.idle":"2022-11-21T12:25:16.352684Z","shell.execute_reply.started":"2022-11-21T12:25:16.34553Z","shell.execute_reply":"2022-11-21T12:25:16.35138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def evaluate(model, loader_val, *, compute_score=True, pbar=None):\n    \"\"\"\n    Predict and compute loss and score\n    \"\"\"\n    tb = time.time()\n    was_training = model.training\n    model.eval()\n\n    loss_sum = 0.0\n    n_sum = 0\n    y_all = []\n    y_pred_all = []\n\n    if pbar is not None:\n        pbar = tqdm(desc='Predict', nrows=78, total=pbar)\n\n    for img, y in loader_val:\n        n = y.size(0)\n        img = img.to(device)\n        y = y.to(device)\n\n        with torch.no_grad():\n                y_pred = model(img.to(device))\n\n        loss = criterion(y_pred.view(-1), y)\n\n        n_sum += n\n        loss_sum += n * loss.item()\n\n        y_all.append(y.cpu().detach().numpy())\n        y_pred_all.append(y_pred.sigmoid().squeeze().cpu().detach().numpy())\n\n        if pbar is not None:\n            pbar.update(len(img))\n        \n        del loss, y_pred, img, y\n\n    loss_val = loss_sum / n_sum\n\n    y = np.concatenate(y_all)\n    y_pred = np.concatenate(y_pred_all)\n\n    score = roc_auc_score(y, y_pred) if compute_score else None\n\n    ret = {'loss': loss_val,\n           'score': score,\n           'y': y,\n           'y_pred': y_pred,\n           'time': time.time() - tb}\n    \n    model.train(was_training)  \n\n    return ret","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:16.354475Z","iopub.execute_input":"2022-11-21T12:25:16.354918Z","iopub.status.idle":"2022-11-21T12:25:16.367494Z","shell.execute_reply.started":"2022-11-21T12:25:16.354883Z","shell.execute_reply":"2022-11-21T12:25:16.366569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_name = 'tf_efficientnet_b6_ns'\nnfold = 5\nkfold = KFold(n_splits=nfold, random_state=42, shuffle=True)\n\nepochs = 25\nbatch_size = 32\nnum_workers = 2\nweight_decay = 1e-6\nmax_grad_norm = 1000\n\nlr_max = 4e-4\nepochs_warmup = 1.0\n \nflip_rate = 0.5 \nfre_shift_rate = 1.0 \ntime_mask_num = 1 \nfreq_mask_num = 2 \n\nfor ifold, (idx_train, idx_test) in enumerate(kfold.split(df)):\n    print('Fold %d/%d' % (ifold, nfold))\n    torch.manual_seed(42 + ifold + 1)\n\n    dataset_train = Dataset('train', df.iloc[idx_train], tfms=True)\n    dataset_val = Dataset('train', df.iloc[idx_test])\n\n    loader_train = torch.utils.data.DataLoader(dataset_train, batch_size=batch_size,\n                     num_workers=num_workers, pin_memory=True, shuffle=True, drop_last=True)\n    loader_val = torch.utils.data.DataLoader(dataset_val, batch_size=batch_size,\n                     num_workers=num_workers, pin_memory=True)\n\n    \n    model = Model(model_name, pretrained=True)\n    model.to(device)\n    model.train()\n\n    optimizer = torch.optim.AdamW(model.parameters(), lr=lr_max, weight_decay=weight_decay)\n\n    \n    nbatch = len(loader_train)\n    warmup = epochs_warmup * nbatch  \n    nsteps = epochs * nbatch        \n\n    scheduler = CosineLRScheduler(optimizer,\n                  warmup_t=warmup, warmup_lr_init=0.0, warmup_prefix=True, \n                  t_initial=(nsteps - warmup), lr_min=1e-6)                \n    \n    time_val = 0.0\n    lrs = []\n\n    tb = time.time()\n    print('Epoch   loss          score   lr')\n    for iepoch in range(epochs):\n        loss_sum = 0.0\n        n_sum = 0\n\n        # Train\n        for ibatch, (img, y) in enumerate(loader_train):\n            n = y.size(0)\n            img = img.to(device)\n            y = y.to(device)\n\n            optimizer.zero_grad()\n\n            y_pred = model(img)\n            loss = criterion(y_pred.view(-1), y)\n\n            loss_train = loss.item()\n            loss_sum += n * loss_train\n            n_sum += n\n\n            loss.backward()\n\n            grad_norm = torch.nn.utils.clip_grad_norm_(model.parameters(),\n                                                       max_grad_norm)\n            optimizer.step()\n            \n            scheduler.step(iepoch * nbatch + ibatch + 1)\n            lrs.append(optimizer.param_groups[0]['lr'])            \n\n        # Evaluate\n        val = evaluate(model, loader_val)\n        time_val += val['time']\n        loss_train = loss_sum / n_sum\n        lr_now = optimizer.param_groups[0]['lr']\n        dt = (time.time() - tb) / 60\n        print('Epoch %d %.4f %.4f %.4f  %.2e  %.2f min' %\n              (iepoch + 1, loss_train, val['loss'], val['score'], lr_now, dt))\n\n    dt = time.time() - tb\n    print('Training done %.2f min total, %.2f min val' % (dt / 60, time_val / 60))\n\n   \n    ofilename = 'model%d.pytorch' % ifold\n    torch.save(model.state_dict(), ofilename)\n    print(ofilename, 'written')","metadata":{"execution":{"iopub.status.busy":"2022-11-21T12:25:16.369021Z","iopub.execute_input":"2022-11-21T12:25:16.369366Z","iopub.status.idle":"2022-11-21T16:28:54.521447Z","shell.execute_reply.started":"2022-11-21T12:25:16.369332Z","shell.execute_reply":"2022-11-21T16:28:54.520152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.title('LR Schedule: Cosine with linear warmup')\nplt.xlabel('steps')\nplt.ylabel('learning rate')\nplt.plot(lrs)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-21T16:29:34.355041Z","iopub.execute_input":"2022-11-21T16:29:34.355454Z","iopub.status.idle":"2022-11-21T16:29:34.591891Z","shell.execute_reply.started":"2022-11-21T16:29:34.355419Z","shell.execute_reply":"2022-11-21T16:29:34.590969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submit = pd.read_csv(di + '/sample_submission.csv')\nif COLAB == False:\n    \n    submit['target'] = 0\n    for i in range(5):\n        model = Model(model_name, pretrained=False)\n        filename = f'model{i}.pytorch'\n        model.to(device)\n        model.load_state_dict(torch.load(filename, map_location=device))\n        model.eval()\n\n        dataset_test = Dataset('test', submit)\n        loader_test = torch.utils.data.DataLoader(dataset_test, batch_size=64,\n                                                num_workers=num_workers, pin_memory=True)\n\n        test = evaluate(model, loader_test, compute_score=False, pbar=len(submit))\n\n        submit['target'] += test['y_pred']/5\nsubmit.to_csv('submission-5folds.csv', index=False)\nprint('target range [%.2f, %.2f]' % (submit['target'].min(), submit['target'].max()))","metadata":{"execution":{"iopub.status.busy":"2022-11-21T16:30:03.603941Z","iopub.execute_input":"2022-11-21T16:30:03.604421Z","iopub.status.idle":"2022-11-21T18:27:10.685839Z","shell.execute_reply.started":"2022-11-21T16:30:03.604381Z","shell.execute_reply":"2022-11-21T18:27:10.683556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Did it work?**\n\nEach sample is comprised of a set of Short-time Fourier Transforms (SFTs) and corresponding GPS time stamps for each interferometer. The SFTs are not always contiguous in time, since the interferometers are not continuously online.\n\nThe simulated signals are present throughout the entire duration of the set of SFTs in both detectors. The signals are characterised by the location and orientation of the hypothetical astrophysical source as well as two intrinsic parameters: frequency and spin-down. In total there are eight parameters which have all been randomised. These are not provided as part of the data. The typical amplitudes of the resulting signals are one or two orders of magnitude lower than the amplitude of the detector noise.\n\n**What did you not understand about it?**\n\nWell, everything provides in the competition data page. I've no problem while working on it. My work will help scientists detect something new: a second class of gravitational waves! The first gravitational wave discoveries earned a Nobel Prize. Further study of these waves may enable scientists to learn about the structure of the most extreme stars in our universe.\n\n**Share your feedback to accomplish the goals!!!**","metadata":{}}]}