{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":53666,"databundleVersionId":6589269,"sourceType":"competition"},{"sourceId":6864546,"sourceType":"datasetVersion","datasetId":3945246},{"sourceId":8402711,"sourceType":"datasetVersion","datasetId":4999796},{"sourceId":149045505,"sourceType":"kernelVersion"}],"dockerImageVersionId":30559,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Child Mind Sleep States\n\nI'll introduce this notebook as a way to use the KLDivLoss, as well as an architecture that is a mix of GRU cells and the UNET architecture.  It can be improved, and done rather easily.  I will discuss the architecture below.\n\nI have used the Kullback-Leibler divergence loss from torch.  \n\n`L(y_pred, y_true) = y_true * (log y_true - log y_pred)`\n\n\nMy reasoning:\n\n1. It is easily interpretable as a probability distribution\n2. Ensembles are also easy to interpret, no matter how differently they are trained.\n\nWe are going to predict the the critical points, i.e. where the onset and wakeup events happen.  We will take as the target a gaussian of a tunable width centered around the actual point.\n\nThe training code:\n\nhttps://www.kaggle.com/code/danielphalen/cmss-grunet-train\n\n\nThis, of course, draws heavily from @werus23: \n\nhttps://www.kaggle.com/code/werus23/sleep-critical-point-train/notebook\n\nhttps://www.kaggle.com/code/werus23/sleep-critical-point-infer?scriptVersionId=147143158\n\nWhich is based on the wonderful discussion pose of @tolgadincer: \n\nhttps://www.kaggle.com/competitions/child-mind-institute-detect-sleep-states/discussion/441470\n\n\nThere could be bugs and other issues in here.  No promises!","metadata":{}},{"cell_type":"code","source":"from math import pi, sqrt, exp\nimport torch\nfrom torch import nn,Tensor\nimport torch.nn.functional as F\nfrom torch.utils.data import DataLoader, Dataset, SubsetRandomSampler\n\nimport pandas as pd\nimport numpy as np\nimport gc\nimport math\n\nimport datetime\nfrom tqdm.auto import tqdm\n\nimport ctypes\nimport torch\nfrom torch.utils.data import Dataset, DataLoader\n\nimport copy\n\nclass PATHS:\n    MAIN_DIR = \"/kaggle/input/child-mind-institute-detect-sleep-states/\"\n    # CSV FILES : \n    SUBMISSION = MAIN_DIR + \"sample_submission.csv\"\n    TRAIN_EVENTS = MAIN_DIR + \"train_events.csv\"\n    # PARQUET FILES:\n    TRAIN_SERIES = MAIN_DIR + \"train_series.parquet\"\n    TEST_SERIES = MAIN_DIR + \"test_series.parquet\"\n    # TEST_SERIES = MAIN_DIR + \"train_series.parquet\"\n    \n    WORKING_DIR = '/kaggle/working/'\n    \n    # Version 3\n    MODEL = ['/kaggle/input/models/grunet_lightweight_data0.pth',\n             ]\n    \n    @staticmethod\n    def get_series_filename(series_id):\n        f = f'{series_id}_test_series.parquet'\n        return PATHS.WORKING_DIR + f\n    \nclass CFG:\n    DEMO_MODE = False\n    VERBOSE = True\n    \n    DEVICE = 'cuda' if torch.cuda.is_available() else 'cpu'\n    \ndef clean_memory():\n    gc.collect()\n    ctypes.CDLL(\"libc.so.6\").malloc_trim(0)\n    torch.cuda.empty_cache()\n   \n    \narch = [(2, 8,  2, 17),\n        (8, 32,  2, 11),\n        (32, 64, 2, 7)]\n\nin_channels = 2\nhidden_size = arch[-1][1]\nkernel_size=25\nstride = arch[-1][0]\ndilation = 1\nn_layers = 5\ndconv_padding=5\nlen_mult = 2**len(arch)","metadata":{"execution":{"iopub.status.busy":"2024-05-13T16:06:01.212491Z","iopub.execute_input":"2024-05-13T16:06:01.212768Z","iopub.status.idle":"2024-05-13T16:06:05.319974Z","shell.execute_reply.started":"2024-05-13T16:06:01.212742Z","shell.execute_reply":"2024-05-13T16:06:05.318970Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load data and convert to individual frames for later processing","metadata":{}},{"cell_type":"code","source":"t1 = datetime.datetime.now()","metadata":{"execution":{"iopub.status.busy":"2024-05-13T16:06:05.321691Z","iopub.execute_input":"2024-05-13T16:06:05.322109Z","iopub.status.idle":"2024-05-13T16:06:05.326691Z","shell.execute_reply.started":"2024-05-13T16:06:05.322083Z","shell.execute_reply":"2024-05-13T16:06:05.325623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test = pd.read_parquet(PATHS.TEST_SERIES)\nseries_ids = list(test.series_id.unique())\nif len(series_ids) == 3:\n    print('Is test script')\n    test = pd.read_parquet(PATHS.TRAIN_SERIES)\n    series_ids = list(test.series_id.unique())[:10]\n    \nlen(series_ids)","metadata":{"execution":{"iopub.status.busy":"2024-05-13T16:06:05.328205Z","iopub.execute_input":"2024-05-13T16:06:05.328685Z","iopub.status.idle":"2024-05-13T16:06:59.629356Z","shell.execute_reply.started":"2024-05-13T16:06:05.328653Z","shell.execute_reply":"2024-05-13T16:06:59.628426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"series_ids","metadata":{"execution":{"iopub.status.busy":"2024-05-13T16:06:59.631923Z","iopub.execute_input":"2024-05-13T16:06:59.632650Z","iopub.status.idle":"2024-05-13T16:06:59.638355Z","shell.execute_reply.started":"2024-05-13T16:06:59.632612Z","shell.execute_reply":"2024-05-13T16:06:59.637468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Takes about 30 minutes for 277, and 28.2 GB of memory\nfor series_id in tqdm(series_ids):\n    s = test[test.series_id == series_id]\n    filename = PATHS.get_series_filename(series_id)\n    s.to_parquet(filename)","metadata":{"execution":{"iopub.status.busy":"2024-05-13T16:06:59.639696Z","iopub.execute_input":"2024-05-13T16:06:59.640045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del test\nclean_memory()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Functions to produce predictions and datasets","metadata":{}},{"cell_type":"code","source":"def get_predictions(res_df, target, SIGMA):\n    \"\"\"Function will take a predicted dataframe, and find local maxima to get event location.  The score is determined by the area under the curve\"\"\"\n    q = res_df[target].max() * 0.1\n    tmp = res_df.loc[res_df[target] > q].copy()\n    # print(f'Target max = {q}, len = {len(tmp)}')\n    tmp['gap'] = tmp['step'].diff()\n    tmp = tmp[tmp['gap'] > 5*5]\n    # print(f'Target max = {q}, len = {len(tmp)}')\n    res = []\n    for i in range(len(tmp) + 1):\n        start_i = 0 if i == 0 else tmp['step'].iloc[i-1]\n        end_i = tmp['step'].iloc[i] if i < len(tmp) else res_df['step'].max()\n        v = res_df.loc[(res_df['step'] > start_i) & (res_df['step'] < end_i)]\n        if v[target].max() > q:\n            # print('Locate in ', start_i, end_i)\n            idx = v.idxmax()[target]\n            step = v.loc[idx, 'step']\n            span = 3*SIGMA\n            score = res_df.loc[(res_df['step'] > step - span) & (res_df['step'] < step + span), target].sum()\n            res.append([step, target, score])\n            \n    return res","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_predictions(test_ds, i, net):\n    net.eval()\n    with torch.no_grad():\n        series_id = test_ds.series_ids[i]\n        data = test_ds.load_data(series_id)\n        X = test_ds[i]\n        pred = torch.zeros(X.shape).to(CFG.DEVICE, non_blocking=True)\n\n        h = None\n\n        seq_len = X.shape[0]\n        for j in range(0, seq_len, max_chunk_size):\n            X_chunk = X[j: j + max_chunk_size].float().to(CFG.DEVICE, non_blocking=True)\n            y_pred, h = net(X_chunk, h)\n            h = [hi.detach() for hi in h]\n            \n            # Trim edges where convolution not done well\n            v = torch.min(y_pred)\n            y_pred[:10, :] = v\n            y_pred[-10:, :] = v\n            \n            pred[j: j+max_chunk_size, :] = y_pred\n\n            del X_chunk, y_pred\n        clean_memory()\n    res_df = pd.DataFrame(torch.softmax(pred.cpu(), axis=0).numpy(), columns=['wakeup', 'onset'])   \n    res_df['step'] = data['step'].values\n    return res_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class SleepDatasetTrain(Dataset):\n    \"\"\"\n    We load the dataset, and select only the longest continuous chunk with events each night\n    \n    \"\"\"\n    def __init__(\n        self,\n        series_ids,\n        events,\n        len_mult,\n        continuous = None,\n        sigma = None\n    ):\n        self.series_ids = series_ids\n        self.continuous = continuous\n        self.len_mult = len_mult\n        if events is not None:\n            self.events = events\n            self.sigma = sigma\n        else:\n            self.events = None\n            self.sigma = None\n    \n    def load_data(self, series_id):\n        filename = PATHS.get_series_filename(series_id)\n        data = pd.read_parquet(filename)\n        if self.events is not None:\n            if self.continuous is not None:\n                start, end = self.continuous[series_id]\n            else:\n                start, end = 0, 1000000\n            gap = 6*60*12\n            tmp = self.events[(self.events.series_id == series_id) & (self.events.night >= start) & (self.events.night <= end)]\n            data = data[(data.step > (tmp.step.min() - gap)) & (data.step < (tmp.step.max() + gap))]\n            \n            data = data.set_index(['series_id', 'step']).join(tmp.set_index(['series_id', 'step'])[['event', 'night']]).reset_index()\n            norm = 1/ np.sqrt(pi / self.sigma)\n            for evt in ['wakeup', 'onset']:\n                steps = data[data.event == evt]['step'].values\n                col = f'{evt}_val'\n                data[col] = 0.0\n                for i in steps:\n                    x = 0.5*((data.step.astype(np.int64) - i)/self.sigma)**2\n                    data[col] += np.exp(-x)*norm\n                data[col] /= data[col].sum()\n                \n        n = int((len(data) // len_mult) * len_mult)\n        \n        return data.iloc[:n]\n        \n    def __len__(self):\n        return len(self.series_ids)\n\n    def __getitem__(self, index):\n        series_id = self.series_ids[index]\n        data = self.load_data(series_id)\n        X = data[['anglez','enmo']].values.astype(np.float32)\n        X = torch.from_numpy(X)\n        if self.sigma is not None:\n            Y = data[['wakeup_val', 'onset_val']].values.astype(np.float32)\n            Y = torch.from_numpy(Y)\n            return X, Y\n        else:\n            return X","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Network\n\nI decided that we are detecting edge points, similar to segmentation of images, the a UNET like architecture would be helpful.  Instead of downsampling to get indicators, we let the network learn the convolution.  Since we have wakeup events always after onset events, there is some time component, which is why the GRU cells are in the bottleneck of the UNET.\n\n\nPossible improvements:\n\n1. The encoder layers and decoder layers could better include more convolutions\n2. We can add skip connections.\n3. We could just put more indicators in the input.  \n4. We can add a time component to the input.\n5. We can increase the size of the hidden layers.\n\n\nThe Residual GRU is from \n\nhttps://www.kaggle.com/competitions/tlvmc-parkinsons-freezing-gait-prediction/discussion/416410","metadata":{}},{"cell_type":"code","source":"# Made using Draw.io\nfrom IPython.display import Image\nImage('/kaggle/input/grunet-diagram/grunet.jpg')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ResidualBiGRU(nn.Module):\n    def __init__(self, hidden_size, n_layers=1, bidir=True):\n        super(ResidualBiGRU, self).__init__()\n\n        self.hidden_size = hidden_size\n        self.n_layers = n_layers\n\n        self.gru = nn.GRU(\n            hidden_size,\n            hidden_size,\n            n_layers,\n            batch_first=True,\n            bidirectional=bidir,\n        )\n        dir_factor = 2 if bidir else 1\n        self.fc1 = nn.Linear(\n            hidden_size * dir_factor, hidden_size * dir_factor * 2\n        )\n        self.ln1 = nn.LayerNorm(hidden_size * dir_factor * 2)\n        self.fc2 = nn.Linear(hidden_size * dir_factor * 2, hidden_size)\n        self.ln2 = nn.LayerNorm(hidden_size)\n\n    def forward(self, x, h=None):\n        res, new_h = self.gru(x, h)\n        # res.shape = (batch_size, sequence_size, 2*hidden_size)\n\n        res = self.fc1(res)\n        res = self.ln1(res)\n        res = nn.functional.relu(res)\n\n        res = self.fc2(res)\n        res = self.ln2(res)\n        res = nn.functional.relu(res)\n\n        # skip connection\n        res = res + x\n\n        return res, new_h\n\nclass MultiResidualBiGRU(nn.Module):\n    def __init__(self, input_size, hidden_size, out_size, n_layers, bidir=True):\n        super(MultiResidualBiGRU, self).__init__()\n\n        self.input_size = input_size\n        self.hidden_size = hidden_size\n        self.out_size = out_size\n        self.n_layers = n_layers\n\n        self.fc_in = nn.Linear(input_size, hidden_size)\n        self.ln = nn.LayerNorm(hidden_size)\n        self.res_bigrus = nn.ModuleList(\n            [\n                ResidualBiGRU(hidden_size, n_layers=1, bidir=bidir)\n                for _ in range(n_layers)\n            ]\n        )\n        self.fc_out = nn.Linear(hidden_size, out_size)\n\n    def forward(self, x, h=None):\n        # if we are at the beginning of a sequence (no hidden state)\n        if h is None:\n            # (re)initialize the hidden state\n            h = [None for _ in range(self.n_layers)]\n\n        x = self.fc_in(x)\n        x = self.ln(x)\n        x = nn.functional.relu(x)\n\n        new_h = []\n        for i, res_bigru in enumerate(self.res_bigrus):\n            x, new_hi = res_bigru(x, h[i])\n            new_h.append(new_hi)\n\n        x = self.fc_out(x)\n#         x = F.normalize(x,dim=0)\n        return x, new_h  # log probabilities + hidden states\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class EncoderLayer(nn.Module):\n    def __init__(self, in_channels, hidden_size, kernel_size, stride, padding, dilation, use_layernorm, print_shape):\n        super(EncoderLayer, self).__init__()\n        self.conv = nn.Conv1d(in_channels, hidden_size, kernel_size, stride, padding=padding, dilation=1)\n        self.ln = nn.LayerNorm(hidden_size) if use_layernorm else None\n        self.print_shape = print_shape\n        \n    def forward(self, x):\n        x = self.conv(x.transpose(-1,-2))\n        if self.print_shape:\n            print('After Conv', x.shape)\n        if self.ln is not None:\n            x = self.ln(x.transpose(-1, -2))\n        else:\n            x = x.transpose(-1,-2)\n        if self.print_shape:\n            print('After Layernorm', x.shape)\n        x = nn.functional.relu(x)\n        return x","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class GRUNET(nn.Module):\n    def __init__(self, arch, out_channels, kernel_size, stride, dconv_padding, hidden_size, n_layers, bidir=True, print_shape=False):\n        super(GRUNET, self).__init__()\n\n        self.input_size = in_channels\n        self.kernel_size = kernel_size\n        self.stride = stride\n        self.hidden_size = hidden_size\n        # self.out_size = out_size\n        self.n_layers = n_layers\n        self.padding = kernel_size//2\n        self.print_shape = print_shape\n        self.arch = arch\n        self.dilation = 1\n        assert arch[-1][1] == hidden_size\n\n        self.conv = nn.Sequential(*[EncoderLayer(in_chan, out_chan, ksize, stride=stride, padding=ksize//2, dilation=self.dilation, use_layernorm=True, print_shape=print_shape) for in_chan, out_chan, stride, ksize in self.arch])\n        self.res_bigrus = nn.ModuleList(\n            [\n                ResidualBiGRU(hidden_size, n_layers=1, bidir=bidir)\n                for _ in range(n_layers)\n            ]\n        )\n        self.dconv = nn.Sequential(*sum([[nn.ConvTranspose1d(out_chan, in_chan, ksize, stride=stride, padding=ksize//2, dilation=self.dilation, output_padding=1), \n                                  nn.Conv1d(in_chan, in_chan, ksize, stride=1, padding=ksize//2, dilation=self.dilation), nn.ReLU(), \n                                  nn.Conv1d(in_chan, in_chan, ksize, stride=1, padding=ksize//2, dilation=self.dilation), nn.ReLU()] for in_chan, out_chan, stride, ksize in reversed(arch)], []))\n        self.output_layer = nn.Conv1d(2, 2, kernel_size=1, stride=1)\n        \n    def forward(self, x, h=None):\n        # if we are at the beginning of a sequence (no hidden state)\n        init_shape = x.shape\n        if h is None:\n            # (re)initialize the hidden state\n            h = [None for _ in range(self.n_layers)]\n\n        if self.print_shape:\n            print('In', x.shape)\n        x = self.conv(x)\n        if self.print_shape:\n            print('After EncoderLayer', x.shape)\n        new_h = []\n        for i, res_bigru in enumerate(self.res_bigrus):\n            x, new_hi = res_bigru(x, h[i])\n            new_h.append(new_hi)\n        if self.print_shape:\n            print('After GRU', x.shape)\n        x = self.dconv(x.transpose(-1, -2))\n        if self.print_shape:\n            print('After DConv', x.shape)\n            \n        x = self.output_layer(x)\n        x = x.transpose(-1,-2)\n            \n        if self.print_shape:\n            print('After SmoothConv', x.shape)\n        \n        return x, new_h  # probabilities + hidden states\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"net = GRUNET(arch=arch,out_channels=2, hidden_size=hidden_size, kernel_size=kernel_size, stride=stride, dconv_padding=dconv_padding, n_layers=n_layers, bidir=True, print_shape=True)\nX = torch.randn(1040, 2).float()\nZ, h = net(X)\nassert X.shape == Z.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Inference","metadata":{}},{"cell_type":"code","source":"max_chunk_size = 24*60*12  ## Process each day at a time with the network\nmin_interval = 30  ## minimum sleep interval\nSIGMA = 36","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_ds = SleepDatasetTrain(series_ids, events=None, len_mult=len_mult, sigma=None)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nets = []\nfor model in PATHS.MODEL:\n    net = GRUNET(arch=arch,out_channels=2, hidden_size=hidden_size, kernel_size=kernel_size, stride=stride, \n                 dconv_padding=dconv_padding, n_layers=n_layers, bidir=True, print_shape=False).to(CFG.DEVICE)\n    net.load_state_dict(torch.load(model, map_location=CFG.DEVICE))\n    net.eval()\n    nets.append(net)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_df = []\nwith torch.no_grad():\n    for i in tqdm(range(len(test_ds))):\n        series_id = test_ds.series_ids[i]\n        # print(series_id)\n        data = test_ds.load_data(series_id)\n        res_df = None\n        for net in nets:\n            tmp = create_predictions(test_ds, i, net)\n            if res_df is None:\n                res_df = tmp\n            else:\n                res_df['onset'] += tmp['onset']\n                res_df['wakeup'] += tmp['wakeup']\n                \n        res_df['onset'] /= len(nets)\n        res_df['wakeup'] /= len(nets)\n        onset_pred = get_predictions(res_df, target='onset', SIGMA=SIGMA)\n        wakeup_pred = get_predictions(res_df, target='wakeup', SIGMA=SIGMA)\n        pred_df = pd.DataFrame(wakeup_pred + onset_pred, columns=['step', 'event', 'score'])\n        pred_df['series_id'] = series_id\n        pred_df['row_id'] = pred_df.index\n        pred_df = pred_df.sort_values(by='step').drop_duplicates(subset='step').reset_index(drop=True)\n        min_step, max_step = res_df.step.min(), res_df.step.max()\n        pred_df = pred_df[(pred_df.step > min_step + 12*60) & (pred_df.step < max_step - 12*60)]\n\n        all_df.append(pred_df)\n        clean_memory()\n\n    pred_df = pd.concat(all_df).reset_index(drop=True)\n    pred_df['row_id'] = pred_df.index\n    pred_df = pred_df[['row_id', 'series_id', 'step', 'event', 'score']]\n    pred_df = pred_df.sort_values(by=['series_id', 'step'])\n    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred_df[['row_id', 'series_id', 'step', 'event', 'score']].to_csv('submission.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!head submission.csv","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t2 = datetime.datetime.now()\nprint(f'Ran in {t2-t1}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from plotly.offline import init_notebook_mode\ninit_notebook_mode()\nfrom IPython.display import Markdown\n\nimport dateutil.relativedelta as rd\nfrom plotly import __version__\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, iplot\nfrom plotly.graph_objs import *\nfrom plotly.tools import FigureFactory as FF","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"iplot({'data': [Scatter(x=res_df['step'], y=res_df['wakeup'], name='wakeup'), \n                Scatter(x=res_df['step'], y=res_df['onset'], name='onset')], \n       'layout': Layout(title='Event Probabilities')})","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}