{"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"}],"dockerImageVersionId":30588,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import gc\nimport numpy as np\nimport pandas as pd\nimport polars as pl\nfrom tqdm.auto import tqdm\nfrom matplotlib import pyplot as plt\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.metrics import mean_squared_error\nfrom scipy.signal import find_peaks\n\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader\nfrom torch.optim.lr_scheduler import CosineAnnealingLR\nfrom torchinfo import summary","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-12-08T06:50:15.43219Z","iopub.execute_input":"2023-12-08T06:50:15.432645Z","iopub.status.idle":"2023-12-08T06:50:19.705103Z","shell.execute_reply.started":"2023-12-08T06:50:15.432622Z","shell.execute_reply":"2023-12-08T06:50:19.704155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"INPUT_DIR = \"/kaggle/input/child-mind-institute-detect-sleep-states/\"","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:50:19.706623Z","iopub.execute_input":"2023-12-08T06:50:19.707079Z","iopub.status.idle":"2023-12-08T06:50:19.710687Z","shell.execute_reply.started":"2023-12-08T06:50:19.707055Z","shell.execute_reply":"2023-12-08T06:50:19.709916Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1. Data Preparation","metadata":{}},{"cell_type":"code","source":"df_series = pl.read_parquet(INPUT_DIR + \"train_series.parquet\")\ndf_events = pl.read_csv(INPUT_DIR + \"train_events.csv\")\ndf_events = df_events.with_columns(\n    pl.col(\"event\").map_dict({\"wakeup\": 1.0, \"onset\": -1.0}, return_dtype=pl.Float32)\n)","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:50:19.711705Z","iopub.execute_input":"2023-12-08T06:50:19.712127Z","iopub.status.idle":"2023-12-08T06:50:32.344772Z","shell.execute_reply.started":"2023-12-08T06:50:19.712105Z","shell.execute_reply":"2023-12-08T06:50:32.343821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_groupby = df_series.group_by(\"series_id\", maintain_order=True)\n\ndict_valid_ratio = dict()\nlist_feature_array = []\nlist_df_1min = []\nfor series_id, df in tqdm(df_groupby, total=df_series.get_column(\"series_id\").n_unique()):\n    df = (\n        df.join(\n            df_events.filter(pl.col(\"series_id\") == series_id).select(\"timestamp\", \"event\"),\n            on=\"timestamp\", how=\"left\"\n        ).with_columns(\n            pl.col(\"timestamp\").str.to_datetime(),\n            pl.col(\"event\").fill_null(0.0),\n        ).with_columns(\n            pl.col(\"timestamp\").dt.date().cast(str).alias(\"date\"),\n            pl.col(\"timestamp\").dt.time().cast(str).alias(\"time\"),\n        )\n    ).to_pandas()\n    \n    df[\"timestamp\"] = df[\"timestamp\"].dt.tz_localize(None)\n    dup_count = df.groupby([\"anglez\", \"enmo\", \"time\"])[\"step\"].transform(\"count\")\n    df[\"valid_flag\"] = (dup_count == 1).astype(\"float32\")\n    dict_valid_ratio[series_id] = df[\"valid_flag\"].mean()\n    \n    list_feature_array_tmp = []\n    df[\"log_anglez_std\"] = np.log(df[\"anglez\"].rolling(25, min_periods=1, center=True).std() + 1).astype(\"float32\")\n    df[\"log_enmo\"] = np.log(df[\"enmo\"] + 0.01).astype(\"float32\")\n    for feature in [\"log_anglez_std\", \"log_enmo\", \"valid_flag\"]:\n        df_pivot = df.pivot(index=[\"series_id\", \"date\"], columns=\"time\", values=feature)\n        feature_array = df_pivot.fillna(0).values\n        feature_array_1day_bedore = df_pivot.shift(1).fillna(0).values\n        feature_array_1day_after = df_pivot.shift(-1).fillna(0).values\n        feature_array = np.concatenate([feature_array_1day_bedore[:, -180*12:], feature_array, feature_array_1day_after[:, :180*12]], axis=1)\n        list_feature_array_tmp.append(feature_array)\n    list_feature_array.append(np.stack(list_feature_array_tmp, axis=1))\n    \n    dict_agg = {\"series_id\": \"first\", \"date\": \"first\", \"time\": \"first\", \"step\": \"mean\", \"event\": \"sum\", \"valid_flag\": \"max\"}\n    df_1min = df.resample(\"1min\", on=\"timestamp\").agg(dict_agg).reset_index()\n    df_1min[\"step\"] = df_1min[\"step\"].astype(\"int32\")\n    values_event = df_1min[\"event\"].values\n    values_target = values_event.copy()\n    for j in range(30):\n        weight = np.exp(-j/2.8)\n        values_target[:-(j+1)] += weight*values_event[(j+1):]  # shift backward\n        if j > 0:\n            values_target[j:] += weight*values_event[:-j]  # shift forward\n    df_1min[\"target\"] = values_target\n    list_df_1min.append(df_1min)\n\ndel df_series, df_groupby","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:50:32.346837Z","iopub.execute_input":"2023-12-08T06:50:32.347136Z","iopub.status.idle":"2023-12-08T06:55:49.773288Z","shell.execute_reply.started":"2023-12-08T06:50:32.347112Z","shell.execute_reply":"2023-12-08T06:55:49.772263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12, 3))\nplt.plot(df_1min.iloc[2300:2350][\"step\"], df_1min.iloc[2300:2350][\"target\"], marker=\".\", label=\"target\")\nplt.plot(df_1min.iloc[2300:2350][\"step\"], df_1min.iloc[2300:2350][\"event\"], drawstyle=\"steps-pre\", label=\"event\")\nplt.title(\"visualization example of the target\")\nplt.xlabel(\"step\")\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:55:49.774468Z","iopub.execute_input":"2023-12-08T06:55:49.774798Z","iopub.status.idle":"2023-12-08T06:55:50.01123Z","shell.execute_reply.started":"2023-12-08T06:55:49.774777Z","shell.execute_reply":"2023-12-08T06:55:50.010344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = np.concatenate(list_feature_array)\ndel list_feature_array\nX = (X - X.min(axis=(0, 2), keepdims=True)) / (X.max(axis=(0, 2), keepdims=True) - X.min(axis=(0, 2), keepdims=True))\n\nX.shape","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:55:50.012591Z","iopub.execute_input":"2023-12-08T06:55:50.013167Z","iopub.status.idle":"2023-12-08T06:55:52.491544Z","shell.execute_reply.started":"2023-12-08T06:55:50.013137Z","shell.execute_reply":"2023-12-08T06:55:52.490638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_1min = pd.concat(list_df_1min)\ndel list_df_1min\n\ndf_1min","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:55:52.492531Z","iopub.execute_input":"2023-12-08T06:55:52.492807Z","iopub.status.idle":"2023-12-08T06:55:53.073258Z","shell.execute_reply.started":"2023-12-08T06:55:52.492785Z","shell.execute_reply":"2023-12-08T06:55:53.072341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:55:53.074167Z","iopub.execute_input":"2023-12-08T06:55:53.074441Z","iopub.status.idle":"2023-12-08T06:55:53.616381Z","shell.execute_reply.started":"2023-12-08T06:55:53.074416Z","shell.execute_reply":"2023-12-08T06:55:53.615481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_y = df_1min.pivot(index=[\"series_id\", \"date\"], columns=\"time\", values=\"target\").fillna(0)\ndf_mask = df_1min.pivot(index=[\"series_id\", \"date\"], columns=\"time\", values=\"valid_flag\").fillna(0)\ndf_mask.to_parquet(\"df_mask.parquet\")\n\ndf_y","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:55:53.617415Z","iopub.execute_input":"2023-12-08T06:55:53.617662Z","iopub.status.idle":"2023-12-08T06:56:05.09907Z","shell.execute_reply.started":"2023-12-08T06:55:53.617641Z","shell.execute_reply":"2023-12-08T06:56:05.097959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2. Regression Model (1D-CNN UNet)","metadata":{}},{"cell_type":"markdown","source":"The architecture of my NN model is based on the following notebook written by [@K_mat](https://www.kaggle.com/kmat2019) in the other competition.  \nhttps://www.kaggle.com/code/kmat2019/u-net-1d-cnn-with-keras\n\nThanks and congrats on 2nd Place!","metadata":{}},{"cell_type":"markdown","source":"### 2.1. definition of classes and functions","metadata":{}},{"cell_type":"code","source":"class MyDataset(Dataset):\n    def __init__(self, X, Y, flag):\n        self.X = torch.FloatTensor(X)\n        if Y is not None:\n            self.Y = torch.FloatTensor(Y)\n        self.flag = torch.FloatTensor(flag)\n    \n    def __len__(self):\n        return self.X.shape[0]\n    \n    def __getitem__(self, idx):\n        if \"Y\" in dir(self):\n            return (self.X[idx], self.Y[idx], self.flag[idx])\n        else:\n            return (self.X[idx], torch.Tensor(), self.flag[idx])","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:56:05.103228Z","iopub.execute_input":"2023-12-08T06:56:05.103484Z","iopub.status.idle":"2023-12-08T06:56:05.109998Z","shell.execute_reply.started":"2023-12-08T06:56:05.103461Z","shell.execute_reply":"2023-12-08T06:56:05.109042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ConvBNReLU(nn.Module):\n    def __init__(self, in_channels, out_channels, kernel_size, stride, groups=1):\n        super().__init__()\n        \n        if stride == 1:\n            padding = \"same\"\n        else:\n            padding = (kernel_size - stride) // 2\n        self.layers = nn.Sequential(\n            nn.Conv1d(in_channels, out_channels, kernel_size=kernel_size, stride=stride, padding=padding, groups=groups),\n            nn.BatchNorm1d(out_channels),\n            nn.ReLU()\n        )\n    \n    def forward(self, x):\n        x_out = self.layers(x)\n        return x_out\n\n\nclass SEBlock(nn.Module):\n    def __init__(self, n_channels, se_ratio):\n        super().__init__()\n        \n        self.layers = nn.Sequential(\n            nn.AdaptiveAvgPool1d(output_size=1),  #  Global Average Pooling\n            nn.Conv1d(n_channels, n_channels//se_ratio, kernel_size=1),\n            nn.ReLU(),\n            nn.Conv1d(n_channels//se_ratio, n_channels, kernel_size=1),\n            nn.Sigmoid()\n        )\n\n    def forward(self, x):\n        x_out = torch.mul(x, self.layers(x))\n        return x_out\n\n\nclass ResBlock(nn.Module):\n    def __init__(self, n_channels, kernel_size, se_ratio):\n        super().__init__()\n        \n        self.layers = nn.Sequential(\n            ConvBNReLU(n_channels, n_channels, kernel_size, stride=1),\n            ConvBNReLU(n_channels, n_channels, kernel_size, stride=1),\n            SEBlock(n_channels, se_ratio)\n        )\n    \n    def forward(self, x):\n        x_re = self.layers(x)\n        x_out = x + x_re\n        return x_out\n    \n\nclass UNet1d(nn.Module):\n    def __init__(self, input_channels, initial_channels, initial_kernel_size,\n                 down_channels, down_kernel_size, down_stride, res_depth, res_kernel_size, se_ratio, out_kernel_size):\n        super().__init__()\n        self.down_kernel_size = down_kernel_size\n        self.down_stride = down_stride\n        \n        self.initial_layers = ConvBNReLU(input_channels, initial_channels, initial_kernel_size, stride=1, groups=input_channels)\n        \n        self.down_layers = nn.ModuleList()\n        for i in range(len(down_channels)):\n            if i == 0:\n                in_channels = initial_channels\n            else:\n                in_channels = down_channels[i-1] + input_channels\n            out_channels = down_channels[i]\n            kernel_size = down_kernel_size[i]\n            stride = down_stride[i]\n            \n            block = []\n            block.append(ConvBNReLU(in_channels, out_channels, kernel_size, stride))\n            for j in range(res_depth):\n                block.append(ResBlock(out_channels, res_kernel_size, se_ratio))\n            self.down_layers.append(nn.Sequential(*block))\n        \n        self.up_layers = nn.ModuleList()\n        for i in range(len(down_channels)-1, 0, -1):\n            in_channels = out_channels + down_channels[i]\n            out_channels = down_channels[i]\n            kernel_size = down_kernel_size[i]\n            self.up_layers.append(ConvBNReLU(in_channels, out_channels, kernel_size, stride=1))\n        \n        self.out_layers = nn.Conv1d(down_channels[1], 1, out_kernel_size, padding=\"same\")\n    \n    def forward(self, x):\n        outs = []\n        x_avg = x\n        x = self.initial_layers(x)\n        \n        for i in range(len(self.down_layers)):\n            x_out = self.down_layers[i](x)\n            if i == len(self.down_layers) - 1:\n                x = x_out\n            else:\n                outs.append(x_out)\n                kernel_size = self.down_kernel_size[i]\n                stride = self.down_stride[i]\n                padding = (kernel_size - stride) // 2\n                x_avg = F.avg_pool1d(x_avg, kernel_size, stride, padding)\n                x = torch.cat([x_out, x_avg], dim=1)\n        \n        for i in range(len(self.up_layers)):\n            scale_factor = self.down_stride[-i-1]\n            x = F.interpolate(x, scale_factor=scale_factor, mode=\"linear\")\n            x = torch.cat([x, outs[-i-1]], dim=1)\n            x = self.up_layers[i](x)\n        \n        x_out = self.out_layers(x)\n        x_out = x_out[:, 0, 180:-180]\n        \n        return x_out","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:56:05.111291Z","iopub.execute_input":"2023-12-08T06:56:05.11212Z","iopub.status.idle":"2023-12-08T06:56:05.128503Z","shell.execute_reply.started":"2023-12-08T06:56:05.112094Z","shell.execute_reply":"2023-12-08T06:56:05.127573Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train(model, data_loader, optimizer, criterion, device):\n    model.train()\n    \n    for batch in data_loader:\n        X = batch[0].to(device)\n        Y = batch[1].to(device)\n        mask = batch[2].to(device)\n        \n        preds = model(X) * mask\n        loss = criterion(preds, Y)\n\n        optimizer.zero_grad()\n        loss.backward()\n        optimizer.step()\n\n\ndef evaluate(model, data_loader, criterion, device):\n    model.eval()\n    \n    n = 0\n    total_loss = 0.0\n    for batch in data_loader:\n        X = batch[0].to(device)\n        Y = batch[1].to(device)\n        mask = batch[2].to(device)\n        \n        with torch.no_grad():\n            preds = model(X) * mask\n        \n        loss = criterion(preds, Y)\n        total_loss += loss.item() * X.shape[0]\n        n += X.shape[0]\n    \n    avg_loss = total_loss / n\n    \n    return avg_loss\n\n\ndef predict(model, data_loader, device):\n    model.eval()\n    \n    preds_all = []\n    for batch in data_loader:\n        X = batch[0].to(device)\n        mask = batch[2].to(device)\n        \n        with torch.no_grad():\n            preds = model(X) * mask\n        preds = preds.cpu().numpy()\n        preds_all.append(preds)\n    \n    preds_all = np.concatenate(preds_all)\n        \n    return preds_all","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:56:05.129606Z","iopub.execute_input":"2023-12-08T06:56:05.129835Z","iopub.status.idle":"2023-12-08T06:56:05.143032Z","shell.execute_reply.started":"2023-12-08T06:56:05.129814Z","shell.execute_reply":"2023-12-08T06:56:05.142281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.2. train model","metadata":{}},{"cell_type":"code","source":"n_splits = 10\nn_epochs = 15\nbatch_size = 16\n    \nY = df_y.values\nmask = df_mask.values\ngroups = df_y.index.get_level_values(\"series_id\").tolist()\n\npreds_valid = np.zeros_like(Y)\ngkf = GroupKFold(n_splits=n_splits)\nfor fold, (idx_train, idx_valid) in enumerate(gkf.split(X, Y, groups=groups)):   \n    print(f\"fold: {fold}\")\n    X_train = X[idx_train]\n    Y_train = Y[idx_train]\n    mask_train = mask[idx_train]\n\n    X_valid = X[idx_valid]\n    Y_valid = Y[idx_valid]\n    mask_valid = mask[idx_valid]\n\n    # dataset\n    ds_train = MyDataset(X_train, Y_train, mask_train)\n    ds_valid = MyDataset(X_valid, Y_valid, mask_valid)\n\n    # dataloader\n    dl_train = DataLoader(ds_train, batch_size=batch_size, shuffle=True, \n                          num_workers=0, pin_memory=True, drop_last=True)\n    dl_valid = DataLoader(ds_valid, batch_size=batch_size, shuffle=False, \n                          num_workers=0, pin_memory=True, drop_last=False)\n\n    # build model\n    torch.manual_seed(0)\n    device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n    model = UNet1d(\n        input_channels=X.shape[1],\n        initial_channels=72,\n        initial_kernel_size=15,\n        down_channels=(72, 72, 72),\n        down_kernel_size=(12, 15, 15),\n        down_stride=(12, 9, 5),  # first element must be 12\n        res_depth=3,\n        res_kernel_size=15,\n        se_ratio=4,\n        out_kernel_size=21,\n    )\n    if fold == 0:\n        print(summary(\n            model=model,\n            input_size=(batch_size, X.shape[1], X.shape[2]),\n            col_names=[\"input_size\", \"output_size\", \"num_params\"],\n            col_width=20\n        ))\n    model.to(device)\n    optimizer = optim.Adam(model.parameters(), lr=1e-3)\n    scheduler = CosineAnnealingLR(optimizer, T_max=n_epochs)\n    criterion = nn.MSELoss()\n\n    best_loss = np.inf\n    for epoch in range(n_epochs):\n        train(model, dl_train, optimizer, criterion, device)\n        loss = evaluate(model, dl_valid, criterion, device)\n        scheduler.step()\n        if loss < best_loss:\n            best_loss = loss\n            torch.save(model.state_dict(), f\"model_{fold}.pth\")\n            print(f\"epoch: {epoch}\\tvalid-loss: {loss}\\tbest!\")\n        else:\n            print(f\"epoch: {epoch}\\tvalid-loss: {loss}\")\n\n    with torch.no_grad():\n        model.load_state_dict(torch.load(f\"model_{fold}.pth\"))\n        preds_valid[idx_valid] = predict(model, dl_valid, device)\n    print(f\"best_loss: {best_loss}\")\n    print()\n\nrmse = mean_squared_error(Y, preds_valid, squared=False)\nprint(f\"RMSE: {rmse:.5f}\")","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:56:05.143914Z","iopub.execute_input":"2023-12-08T06:56:05.144187Z","iopub.status.idle":"2023-12-08T06:59:47.384576Z","shell.execute_reply.started":"2023-12-08T06:56:05.144167Z","shell.execute_reply":"2023-12-08T06:59:47.383469Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12, 3))\nplt.plot(np.linspace(0, 24, 60*24), preds_valid.mean(axis=0), label=\"mean of prediction\")\nplt.plot(np.linspace(0, 24, 60*24), df_y.mean().values, label=\"mean of target\")\nplt.xlim(0, 24)\nplt.xlabel(\"hour (UTC)\")\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.385446Z","iopub.status.idle":"2023-12-08T06:59:47.385752Z","shell.execute_reply.started":"2023-12-08T06:59:47.385609Z","shell.execute_reply":"2023-12-08T06:59:47.385625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_pred = pd.DataFrame(preds_valid, index=df_y.index, columns=df_y.columns)\ndf_pred = df_pred.stack().reset_index(name=\"score\")\ndf_pred = pd.merge(\n    df_1min[[\"series_id\", \"date\", \"time\", \"step\", \"event\"]],\n    df_pred,\n    on=[\"series_id\", \"date\", \"time\"],\n    how=\"inner\"\n)\ndf_pred","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.386877Z","iopub.status.idle":"2023-12-08T06:59:47.387178Z","shell.execute_reply.started":"2023-12-08T06:59:47.38704Z","shell.execute_reply":"2023-12-08T06:59:47.387056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i, (series_id, df) in enumerate(df_pred.groupby(\"series_id\")):\n    if i == 10:\n        break\n    plt.figure(figsize=(12, 3))\n    plt.plot(df[\"step\"], df[\"score\"], linewidth=2, label=\"predicted score\")\n    plt.plot(df[\"step\"], df[\"event\"], linewidth=1, label=\"ground truth\")\n    plt.legend()\n    plt.title(series_id)\n    plt.show()\n    plt.close()","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.38792Z","iopub.status.idle":"2023-12-08T06:59:47.388236Z","shell.execute_reply.started":"2023-12-08T06:59:47.388099Z","shell.execute_reply":"2023-12-08T06:59:47.388116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Make Submission","metadata":{}},{"cell_type":"markdown","source":"### 3.1. without post-prosessing (only peak detection)","metadata":{}},{"cell_type":"code","source":"list_df = []\nfor series_id, df in tqdm(df_pred.groupby(\"series_id\")):\n    for event in [\"onset\", \"wakeup\"]:\n        values_step = df[\"step\"].values\n        if event == \"onset\":\n            values_score = -df[\"score\"].values\n        else:\n            values_score = df[\"score\"].values\n\n        peak_idx = find_peaks(values_score, height=0.0, distance=8)[0]\n        df_peak = pd.DataFrame(values_step[peak_idx], columns=[\"step\"])\n        df_peak[\"series_id\"] = series_id\n        df_peak[\"event\"] = event\n        df_peak[\"score\"] = values_score[peak_idx]\n        list_df.append(df_peak)\n\ndf_sub = pd.concat(list_df)\ndf_sub = df_sub.sort_values(\"score\", ascending=False).groupby(\"event\").head(100000)  # avoid Submission Scoring Error\ndf_sub = df_sub.sort_values([\"series_id\", \"step\"]).reset_index(drop=True)\ndf_sub = df_sub[[\"series_id\", \"step\", \"event\", \"score\"]].reset_index(names=\"row_id\")\ndf_sub.to_csv(\"df_sub_without_pp.csv\", index=False)\ndf_sub","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.389935Z","iopub.status.idle":"2023-12-08T06:59:47.390758Z","shell.execute_reply.started":"2023-12-08T06:59:47.390596Z","shell.execute_reply":"2023-12-08T06:59:47.390614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.2 with post-processing","metadata":{}},{"cell_type":"code","source":"df_events = pd.read_csv(INPUT_DIR + \"train_events.csv\").dropna()\ndf_events[\"timestamp\"] = pd.to_datetime(df_events[\"timestamp\"], utc=True).dt.tz_localize(None)\ndf_events[\"time\"] = df_events[\"timestamp\"].dt.time.astype(str)\ndf_events[\"minute_mod15\"] = df_events[\"timestamp\"].dt.minute % 15\n\ndf_agg = df_events.groupby([\"time\", \"event\"], as_index=False).size()\ndf_agg[\"rate\"] = df_agg[\"size\"] / df_agg.groupby(\"event\")[\"size\"].transform(\"sum\") * (60*24)\ndf_time = df_agg.pivot(index=\"time\", columns=\"event\", values=\"rate\").fillna(0).reset_index()\ndf_time = df_time.merge(df_pred[[\"time\"]].drop_duplicates(), how=\"right\").fillna(0)\ndf_time = pd.concat([df_time]*3, ignore_index=True)\ndf_time[\"onset\"] = df_time[\"onset\"].rolling(60, center=True).mean()\ndf_time[\"wakeup\"] = df_time[\"wakeup\"].rolling(60, center=True).mean()\ndf_time = df_time.iloc[60*24:-60*24].reset_index(drop=True)\n\ndf_agg = df_events.groupby([\"minute_mod15\", \"event\"], as_index=False).size()\ndf_agg[\"rate\"] = df_agg[\"size\"] / df_agg.groupby(\"event\")[\"size\"].transform(\"sum\") * 15\ndf_minute = df_agg.pivot(index=\"minute_mod15\", columns=\"event\", values=\"rate\").reset_index()\n\nplt.figure(figsize=(12, 3))\nplt.plot(np.linspace(0, 24, 60*24), df_time[\"onset\"], label=\"onset\")\nplt.plot(np.linspace(0, 24, 60*24), df_time[\"wakeup\"], label=\"wakeup\")\nplt.xlabel(\"hour (UTC)\")\nplt.legend()\nplt.show()\n\nplt.figure(figsize=(12, 3))\nplt.bar(np.linspace(0, 14, 15)-0.1, df_minute[\"onset\"], width=0.2, label=\"onset\")\nplt.bar(np.linspace(0, 14, 15)+0.1, df_minute[\"wakeup\"], width=0.2, label=\"wakeup\")\nplt.xlabel(\"minute mod 15\")\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.391838Z","iopub.status.idle":"2023-12-08T06:59:47.392148Z","shell.execute_reply.started":"2023-12-08T06:59:47.391999Z","shell.execute_reply":"2023-12-08T06:59:47.392014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_time[[\"onset\", \"wakeup\"]] = df_time[[\"onset\", \"wakeup\"]].clip(0.1, 1.1) ** 0.13\ndf_minute[[\"onset\", \"wakeup\"]] = df_minute[[\"onset\", \"wakeup\"]].clip(0.5, 1.3) ** 0.06\n\nplt.figure(figsize=(12, 3))\nplt.plot(np.linspace(0, 24, 60*24), df_time[\"onset\"], label=\"onset\")\nplt.plot(np.linspace(0, 24, 60*24), df_time[\"wakeup\"], label=\"wakeup\")\nplt.ylim(0, 1.1)\nplt.xlabel(\"hour (UTC)\")\nplt.legend()\nplt.show()\n\nplt.figure(figsize=(12, 3))\nplt.bar(np.linspace(0, 14, 15)-0.1, df_minute[\"onset\"], width=0.2, label=\"onset\")\nplt.bar(np.linspace(0, 14, 15)+0.1, df_minute[\"wakeup\"], width=0.2, label=\"wakeup\")\nplt.ylim(0, 1.1)\nplt.xlabel(\"minute mod 15\")\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.393469Z","iopub.status.idle":"2023-12-08T06:59:47.394334Z","shell.execute_reply.started":"2023-12-08T06:59:47.394129Z","shell.execute_reply":"2023-12-08T06:59:47.394149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_pred[\"minute_mod15\"] = df_pred[\"time\"].str[3:5].astype(int) % 15","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.395466Z","iopub.status.idle":"2023-12-08T06:59:47.396103Z","shell.execute_reply.started":"2023-12-08T06:59:47.39587Z","shell.execute_reply":"2023-12-08T06:59:47.39589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_df = []\nfor series_id, df in tqdm(df_pred.groupby(\"series_id\")):\n    df = df.merge(df_time, how=\"left\", on=\"time\")\n    df = df.merge(df_minute, how=\"left\", on=\"minute_mod15\")\n    \n    df_tmp = df.copy()\n    df_tmp[\"score\"] = df_tmp[\"score\"].replace(0.0, np.nan)\n    df_tmp = df_tmp.groupby(\"time\")[\"score\"].mean()\n    df_tmp = pd.concat([df_tmp]*3).rolling(90, center=True, min_periods=1).mean()\n    df_tmp = df_tmp.iloc[60*24:-60*24].reset_index().rename({\"score\": \"score_mean\"}, axis=1)\n    df = df.merge(df_tmp, on=\"time\", how=\"left\")\n    \n    df[\"score\"] = 0.9*df[\"score\"] + 0.1*df[\"score_mean\"]\n    df[\"score\"] *= np.where(df[\"score\"]>0, df[\"wakeup_x\"], df[\"onset_x\"])\n    df[\"score\"] *= np.where(df[\"score\"]>0, df[\"wakeup_y\"], df[\"onset_y\"])\n    valid_ratio = dict_valid_ratio[series_id]\n    \n    for event in [\"onset\", \"wakeup\"]:\n        values_step = df[\"step\"].values\n        if event == \"onset\":\n            values_score = -df[\"score\"].values\n        else:\n            values_score = df[\"score\"].values\n        \n        # measure peaks\n        peak_idx = find_peaks(values_score, height=0.04, distance=60*16)[0]  # at least 16 hours interval\n        df_measure_peak = pd.DataFrame(values_step[peak_idx], columns=[\"step\"])\n        df_measure_peak[\"series_id\"] = series_id\n        df_measure_peak[\"event\"] = event\n        df_measure_peak[\"score\"] = values_score[peak_idx] * 4 * valid_ratio**0.15\n        \n        # minor peaks\n        peak_idx = find_peaks(values_score, height=0.0, distance=6)[0]\n        df_minor_peak = pd.DataFrame(values_step[peak_idx], columns=[\"step\"])\n        df_minor_peak[\"series_id\"] = series_id\n        df_minor_peak[\"event\"] = event\n        df_minor_peak[\"score\"] = values_score[peak_idx]\n        \n        df_peak = pd.concat([df_measure_peak, df_minor_peak]).drop_duplicates(subset=[\"step\"])\n        list_df.append(df_peak)\n        \ndf_sub = pd.concat(list_df)\ndf_sub = df_sub.sort_values(\"score\", ascending=False).groupby(\"event\").head(100000)  # avoid Submission Scoring Error\ndf_sub = df_sub.sort_values([\"series_id\", \"step\"]).reset_index(drop=True)\ndf_sub = df_sub[[\"series_id\", \"step\", \"event\", \"score\"]].reset_index(names=\"row_id\")\ndf_sub.to_csv(\"df_sub_with_pp.csv\", index=False)\ndf_sub","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.397422Z","iopub.status.idle":"2023-12-08T06:59:47.397815Z","shell.execute_reply.started":"2023-12-08T06:59:47.397623Z","shell.execute_reply":"2023-12-08T06:59:47.397641Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Scoring","metadata":{}},{"cell_type":"markdown","source":"The fast evaluation functions implemented by [@213tubo](https://www.kaggle.com/tubotubo) were very beneficial in my experiments and saved me a lot of time. Thanks\\!  \nhttps://github.com/tubo213/kaggle-child-mind-institute-detect-sleep-states/blob/main/src/utils/metrics.py","metadata":{}},{"cell_type":"code","source":"from bisect import bisect_left\nfrom typing import Dict, List, Tuple\n\n\nclass ParticipantVisibleError(Exception):\n    pass\n\n\n# Set some placeholders for global parameters\nseries_id_column_name = \"series_id\"\ntime_column_name = \"step\"\nevent_column_name = \"event\"\nscore_column_name = \"score\"\nuse_scoring_intervals = False\ntolerances = {\n    \"onset\": [12, 36, 60, 90, 120, 150, 180, 240, 300, 360],\n    \"wakeup\": [12, 36, 60, 90, 120, 150, 180, 240, 300, 360],\n}\n\n\ndef score(\n    solution: pd.DataFrame,\n    submission: pd.DataFrame,\n    tolerances: Dict[str, List[float]],\n    series_id_column_name: str,\n    time_column_name: str,\n    event_column_name: str,\n    score_column_name: str,\n    use_scoring_intervals: bool = False,\n    verbose: bool = True,\n) -> Tuple[float, pd.DataFrame, pd.DataFrame]:\n    # Validate metric parameters\n    assert len(tolerances) > 0, \"Events must have defined tolerances.\"\n    assert set(tolerances.keys()) == set(solution[event_column_name]).difference(\n        {\"start\", \"end\"}\n    ), (\n        f\"Solution column {event_column_name} must contain the same events \"\n        \"as defined in tolerances.\"\n    )\n    assert pd.api.types.is_numeric_dtype(\n        solution[time_column_name]\n    ), f\"Solution column {time_column_name} must be of numeric type.\"\n\n    # Validate submission format\n    for column_name in [\n        series_id_column_name,\n        time_column_name,\n        event_column_name,\n        score_column_name,\n    ]:\n        if column_name not in submission.columns:\n            raise ParticipantVisibleError(\n                f\"Submission must have column '{column_name}'.\"\n            )\n\n    if not pd.api.types.is_numeric_dtype(submission[time_column_name]):\n        raise ParticipantVisibleError(\n            f\"Submission column '{time_column_name}' must be of numeric type.\"\n        )\n    if not pd.api.types.is_numeric_dtype(submission[score_column_name]):\n        raise ParticipantVisibleError(\n            f\"Submission column '{score_column_name}' must be of numeric type.\"\n        )\n\n    # Set these globally to avoid passing around a bunch of arguments\n    globals()[\"series_id_column_name\"] = series_id_column_name\n    globals()[\"time_column_name\"] = time_column_name\n    globals()[\"event_column_name\"] = event_column_name\n    globals()[\"score_column_name\"] = score_column_name\n    globals()[\"use_scoring_intervals\"] = use_scoring_intervals\n\n    return event_detection_ap(solution, submission, tolerances, verbose=verbose)\n\n\ndef find_nearest(xs: np.ndarray, value):\n    \"\"\"\n    Find the index of the closest value to x in the array xs.\n    \"\"\"\n    idx = np.searchsorted(xs, value, side=\"left\")\n    best_idx = None\n    best_error = float(\"inf\")\n    best_diff = float(\"inf\")\n\n    range_min = max(0, idx - 1)\n    range_max = min(len(xs), idx + 2)\n    for check_idx in range(\n        range_min, range_max\n    ):  # Check the exact, one before, and one after\n        error = abs(xs[check_idx] - value)\n        if error < best_error:\n            best_error = error\n            best_idx = check_idx\n            best_diff = xs[check_idx] - value\n\n    return best_idx, best_error, best_diff\n\n\ndef find_nearest_time_idx(sorted_gt_times, det_time, excluded_indices: set):\n    \"\"\"\n    search index of gt_times closest to det_time.\n\n    assumes gt_times is sorted in ascending order.\n    \"\"\"\n    # e.g. if gt_times = [0, 1, 2, 3, 4, 5] and det_time = 2.5, then idx = 3\n    sorted_gt_times = np.asarray(sorted_gt_times)\n    available_indices = np.asarray(\n        sorted(set(range(len(sorted_gt_times))) - excluded_indices), dtype=int\n    )\n    sorted_gt_times = sorted_gt_times[available_indices]\n    idx, error, diff = find_nearest(sorted_gt_times, det_time)\n    best_idx = available_indices[idx] if idx is not None else None\n\n    return best_idx, error, diff\n\n\ndef match_detections(\n    tolerance: float, ground_truths: pd.DataFrame, detections: pd.DataFrame\n) -> pd.DataFrame:\n    detections_sorted = detections.sort_values(\n        score_column_name, ascending=False\n    ).dropna()\n    is_matched = np.full_like(detections_sorted[event_column_name], False, dtype=bool)\n    diffs = np.full_like(\n        detections_sorted[event_column_name], float(\"inf\"), dtype=float\n    )\n    ground_truths_times = ground_truths.sort_values(time_column_name)[\n        time_column_name\n    ].to_list()\n    matched_gt_indices: set[int] = set()\n\n    for i, det in enumerate(detections_sorted.itertuples(index=False)):\n        det_time = getattr(det, time_column_name)\n\n        best_idx, best_error, best_diff = find_nearest_time_idx(\n            ground_truths_times, det_time, matched_gt_indices\n        )\n\n        if (best_idx is not None) and (best_error < tolerance):\n            is_matched[i] = True\n            diffs[i] = best_diff\n            matched_gt_indices.add(best_idx)\n\n    detections_sorted[\"matched\"] = is_matched\n    detections_sorted[\"diff\"] = diffs\n    return detections_sorted\n\n\ndef precision_recall_curve(\n    matches: np.ndarray, scores: np.ndarray, p: int\n) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:\n    if len(matches) == 0:\n        return [1], [0], []\n\n    # Sort matches by decreasing confidence\n    idxs = np.argsort(scores, kind=\"stable\")[::-1]\n    scores = scores[idxs]\n    matches = matches[idxs]\n\n    distinct_value_indices = np.where(np.diff(scores))[0]\n    threshold_idxs = np.r_[distinct_value_indices, matches.size - 1]\n    thresholds = scores[threshold_idxs]\n\n    # Matches become TPs and non-matches FPs as confidence threshold decreases\n    tps = np.cumsum(matches)[threshold_idxs]\n    fps = np.cumsum(~matches)[threshold_idxs]\n\n    precision = tps / (tps + fps)\n    precision[np.isnan(precision)] = 0\n    recall = (\n        tps / p\n    )  # total number of ground truths might be different than total number of matches\n\n    # Stop when full recall attained and reverse the outputs so recall is non-increasing.\n    last_ind = tps.searchsorted(tps[-1])\n    sl = slice(last_ind, None, -1)\n\n    # Final precision is 1 and final recall is 0\n    return np.r_[precision[sl], 1], np.r_[recall[sl], 0], thresholds[sl]\n\n\ndef average_precision_score(matches: np.ndarray, scores: np.ndarray, p: int) -> float:\n    precision, recall, _ = precision_recall_curve(matches, scores, p)\n    # Compute step integral\n    return -np.sum(np.diff(recall) * np.array(precision)[:-1])\n\n\ndef event_detection_ap(\n    solution: pd.DataFrame,\n    submission: pd.DataFrame,\n    tolerances: Dict[str, List[float]] = tolerances,\n    progress_bar: bool = True,\n    verbose: bool = True,\n) -> Tuple[float, pd.DataFrame, pd.DataFrame]:\n    # Ensure solution and submission are sorted properly\n    solution = solution.sort_values([series_id_column_name, time_column_name])\n    submission = submission.sort_values([series_id_column_name, time_column_name])\n\n    # Extract scoring intervals.\n    if use_scoring_intervals:\n        raise NotImplementedError(\"Scoring intervals not implemented.\")\n\n    # Extract ground-truth events.\n    ground_truths = solution.query(\"event not in ['start', 'end']\").reset_index(\n        drop=True\n    )\n\n    # Map each event class to its prevalence (needed for recall calculation)\n    class_counts = ground_truths.value_counts(event_column_name).to_dict()\n\n    # Create table for detections with a column indicating a match to a ground-truth event\n    detections = submission.assign(matched=False)\n\n    # Remove detections outside of scoring intervals\n    if use_scoring_intervals:\n        raise NotImplementedError(\"Scoring intervals not implemented.\")\n    else:\n        detections_filtered = detections\n\n    # Create table of event-class x tolerance x series_id values\n    aggregation_keys = pd.DataFrame(\n        [\n            (ev, tol, vid)\n            for ev in tolerances.keys()\n            for tol in tolerances[ev]\n            for vid in ground_truths[series_id_column_name].unique()\n        ],\n        columns=[event_column_name, \"tolerance\", series_id_column_name],\n    )\n\n    # Create match evaluation groups: event-class x tolerance x series_id\n    detections_grouped = aggregation_keys.merge(\n        detections_filtered, on=[event_column_name, series_id_column_name], how=\"left\"\n    ).groupby([event_column_name, \"tolerance\", series_id_column_name])\n    ground_truths_grouped = aggregation_keys.merge(\n        ground_truths, on=[event_column_name, series_id_column_name], how=\"left\"\n    ).groupby([event_column_name, \"tolerance\", series_id_column_name])\n\n    # Match detections to ground truth events by evaluation group\n    pbars = aggregation_keys.itertuples(index=False)\n    if progress_bar:\n        pbars = tqdm(pbars, total=len(aggregation_keys), desc=\"Matching detections\")\n    detections_matched = []\n    for key in pbars:\n        dets = detections_grouped.get_group(key)\n        gts = ground_truths_grouped.get_group(key)\n        detections_matched.append(\n            match_detections(dets[\"tolerance\"].iloc[0], gts, dets)\n        )\n    detections_matched = pd.concat(detections_matched)\n\n    # Compute AP per event x tolerance group\n    event_classes = ground_truths[event_column_name].unique()\n    ap_table = (\n        detections_matched.query(\"event in @event_classes\")\n        .groupby([event_column_name, \"tolerance\"])\n        .apply(\n            lambda group: average_precision_score(\n                group[\"matched\"].to_numpy(),\n                group[score_column_name].to_numpy(),\n                class_counts[group[event_column_name].iat[0]],\n            )\n        )\n        .reset_index()\n        .pivot(index=\"tolerance\", columns=\"event\", values=0)\n    )\n    if verbose:\n        display(ap_table)\n    # Average over tolerances, then over event classes\n    mean_ap = ap_table.mean().mean()\n\n    return mean_ap, ap_table, detections_matched","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.399242Z","iopub.status.idle":"2023-12-08T06:59:47.399736Z","shell.execute_reply.started":"2023-12-08T06:59:47.399582Z","shell.execute_reply":"2023-12-08T06:59:47.399599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sub = pd.read_csv(\"df_sub_without_pp.csv\")\nscore_all, df_score, df_result = score(\n    solution=df_events,\n    submission=df_sub,\n    tolerances=tolerances,\n    series_id_column_name=series_id_column_name,\n    time_column_name=time_column_name,\n    event_column_name=event_column_name,\n    score_column_name=score_column_name\n)\nprint(score_all)","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.400643Z","iopub.status.idle":"2023-12-08T06:59:47.401486Z","shell.execute_reply.started":"2023-12-08T06:59:47.401332Z","shell.execute_reply":"2023-12-08T06:59:47.401348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sub = pd.read_csv(\"df_sub_with_pp.csv\")\nscore_all, df_score, df_result = score(\n    solution=df_events,\n    submission=df_sub,\n    tolerances=tolerances,\n    series_id_column_name=series_id_column_name,\n    time_column_name=time_column_name,\n    event_column_name=event_column_name,\n    score_column_name=score_column_name\n)\nprint(score_all)","metadata":{"execution":{"iopub.status.busy":"2023-12-08T06:59:47.402375Z","iopub.status.idle":"2023-12-08T06:59:47.40268Z","shell.execute_reply.started":"2023-12-08T06:59:47.402513Z","shell.execute_reply":"2023-12-08T06:59:47.402533Z"},"trusted":true},"execution_count":null,"outputs":[]}]}