{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"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"},"papermill":{"default_parameters":{},"duration":202.068623,"end_time":"2023-11-03T00:47:08.842663","environment_variables":{},"exception":true,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2023-11-03T00:43:46.774040","version":"2.4.0"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":53666,"databundleVersionId":6589269,"sourceType":"competition"},{"sourceId":5611341,"sourceType":"datasetVersion","datasetId":3227383},{"sourceId":6809268,"sourceType":"datasetVersion","datasetId":3227387},{"sourceId":7128252,"sourceType":"datasetVersion","datasetId":4112330},{"sourceId":7129379,"sourceType":"datasetVersion","datasetId":4112439},{"sourceId":153570242,"sourceType":"kernelVersion"}],"dockerImageVersionId":30587,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install --no-deps \\\n    /kaggle/input/ruptures117/ruptures-1.1.8-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl \\\n    /kaggle/input/fcwt118/fCWT-0.1.18-cp310-cp310-linux_x86_64.whl","metadata":{"papermill":{"duration":36.0795,"end_time":"2023-11-03T00:44:26.841039","exception":false,"start_time":"2023-11-03T00:43:50.761539","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-12-05T13:30:34.875519Z","iopub.execute_input":"2023-12-05T13:30:34.875826Z","iopub.status.idle":"2023-12-05T13:30:57.113512Z","shell.execute_reply.started":"2023-12-05T13:30:34.875801Z","shell.execute_reply":"2023-12-05T13:30:57.111884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from pathlib import Path\n\nimport fcwt\nimport numpy as np\nimport pandas as pd\nimport polars as pl\nimport ruptures as rpt\nimport stumpy\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom lightgbm import Booster, LGBMClassifier\nfrom sklearn.metrics import confusion_matrix\nfrom sklearn.model_selection import GroupKFold, train_test_split\nfrom torch.optim.lr_scheduler import ReduceLROnPlateau\nfrom torch.utils.data import DataLoader, Dataset\nfrom tsfresh import extract_features","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":5.672184,"end_time":"2023-11-03T00:44:32.519635","exception":false,"start_time":"2023-11-03T00:44:26.847451","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-12-05T13:30:57.115126Z","iopub.execute_input":"2023-12-05T13:30:57.115428Z","iopub.status.idle":"2023-12-05T13:31:04.412755Z","shell.execute_reply.started":"2023-12-05T13:30:57.115402Z","shell.execute_reply":"2023-12-05T13:31:04.411540Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"_________\nLightGBM classification of breaking points","metadata":{}},{"cell_type":"code","source":"fc_params = {\n    \"absolute_sum_of_changes\": None,\n    \"approximate_entropy\": [{\"m\": 2, \"r\": 0.5}, {\"m\": 5, \"r\": 0.9}],\n    \"c3\": [{\"lag\": 1}, {\"lag\": 2}, {\"lag\": 3}],\n    \"cid_ce\": [{\"normalize\": False}],\n    \"kurtosis\": None,\n    \"lempel_ziv_complexity\": [{\"bins\": 3}, {\"bins\": 7}],\n    \"longest_strike_above_mean\": None,\n    \"longest_strike_below_mean\": None,\n    \"mean_abs_change\": None,\n    \"mean_change\": None,\n    \"mean_n_absolute_max\": [{\"number_of_maxima\": 7}],\n    \"mean_second_derivative_central\": None,\n    \"number_crossing_m\": [{\"m\": 0}, {\"m\": -0.5}, {\"m\": 0.5}],\n    \"number_cwt_peaks\": [{\"n\": 3}, {\"n\": 10}],\n    \"number_peaks\": [{\"n\": 3}, {\"n\": 10}, {\"n\": 50}],\n    \"percentage_of_reoccurring_datapoints_to_all_datapoints\": None,\n    \"percentage_of_reoccurring_values_to_all_values\": None,\n    \"permutation_entropy\": [{\"tau\": 1, \"dimension\": 3}],\n    \"ratio_value_number_to_time_series_length\": None,\n    \"spkt_welch_density\": [{\"coeff\": 2}, {\"coeff\": 5}, {\"coeff\": 8}],\n}","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:04.413913Z","iopub.execute_input":"2023-12-05T13:31:04.414446Z","iopub.status.idle":"2023-12-05T13:31:04.421643Z","shell.execute_reply.started":"2023-12-05T13:31:04.414422Z","shell.execute_reply":"2023-12-05T13:31:04.420351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def _calc_frac_day(df: pl.DataFrame):\n    # get local time and calculate fraction of the day\n    timelocal = df[\"timestamp\"].str.slice(0, 19).str.to_datetime()\n    frac_day = timelocal.dt.hour() * 3600\n    frac_day += timelocal.dt.minute() * 60\n    frac_day += timelocal.dt.second()\n    frac_day /= 24 * 6 * 6\n    return frac_day.cast(pl.Int8)\n\n\ndef _fast_wavelet_transform(df, wavelet, n_band=8, n_threads=1):\n    # wavelet transformation and filtering, keep the high-band signal combining enmo and anglez\n    scales = fcwt.Scales(wavelet, fcwt.FCWT_LINFREQS, 1, 0.001, 0.5, n_band)\n    fcwt_obj = fcwt.FCWT(wavelet, n_threads, False, False)\n\n    spectrogram = np.zeros((n_band, len(df)), dtype=np.complex64)\n    fcwt_obj.ccwt(\n        df[\"enmo\"].to_numpy()  # real\n        + np.exp(1j * df[\"anglez\"].to_numpy() * np.pi, dtype=np.complex64),  # imag\n        scales,\n        spectrogram,\n    )\n\n    hi_band = np.abs(spectrogram[: n_band - 1]).mean(axis=0)\n    hi_band /= hi_band.max()\n    return hi_band\n\n\ndef _calc_matrix_profile(df: pl.DataFrame, col=\"enmo\", m=50):\n    # calculate matrix profile for motif uniqueness and non-wearing detection\n    # sleeping periods usually have lower uniqueness than active periods.\n    # set default sequence length m = 50\n\n    vals = df[col].to_numpy().astype(np.float64)\n    approx = stumpy.scrump(\n        vals,\n        m=m,\n        percentage=0.01,\n        pre_scrump=True,\n    )\n    approx.update()\n\n    matrix_profile = np.zeros_like(vals)\n    matrix_profile[: len(approx.P_)] = approx.P_\n\n    return matrix_profile\n\n\ndef _calc_lids(df: pl.DataFrame, col=\"hi_band\"):\n    # Locomotor Inactivity During Sleep (LIDS)\n    # https://www.kaggle.com/competitions/child-mind-institute-detect-sleep-states/discussion/444007\n\n    lids = (\n        df.with_columns(loco_inactivity=(df[col] - 0.02).clip(0, 1))\n        .set_sorted([\"timestamp\"])\n        # .group_by_rolling(\n        .rolling(\n            index_column=\"timestamp\",\n            period=\"300s\",\n        )\n        .agg([pl.sum(\"loco_inactivity\").add(1).pow(-1)])\n        # .group_by_rolling(\n        .rolling(\n            index_column=\"timestamp\",\n            period=\"900s\",\n        )\n        .agg(pl.mean(\"loco_inactivity\"))\n    )\n\n    return lids[\"loco_inactivity\"].alias(\"lids\")\n\n\ndef _add_events(df, events, subject):\n    df = df.join(\n        events.groupby(\"series_id\")\n        .get_group(subject)\n        .dropna()\n        .set_index(\"timestamp\")[[\"night\", \"event\"]],\n        how=\"left\",\n    )\n    df.event = df.event.ffill().replace({\"onset\": 1, \"wakeup\": 0}).fillna(0)\n    df.night = df.night.fillna(0)\n    return df\n\n\ndef _calc_break_points(df: pl.DataFrame, cols, period=\"24h\", every=\"12h\", bp_pen=50):\n    # detect break points\n    break_points = set()\n    for i, d in df.set_sorted([\"timestamp\"]).group_by_dynamic(\n        \"timestamp\", period=period, every=every\n    ):\n        bp = rpt.KernelCPD(kernel=\"rbf\", min_size=15, jump=1).fit_predict(\n            d[cols].fill_null(0).to_numpy(), pen=bp_pen\n        )\n\n        bp = d[\"timestamp\"].to_numpy()[bp[:-1]]\n        break_points.update(bp[:-1])\n\n    break_points = pd.Series(sorted(break_points), name=\"timestamp\")\n    if len(break_points) > 0:\n        break_points = (\n            break_points.drop_duplicates()\n            .dt.tz_localize(tz=\"utc\")\n            .to_frame()\n            .assign(bp=1.0)\n            .set_index(\"timestamp\")\n        )\n\n    return break_points\n\n\ndef get_one_subject(\n    ts: pl.DataFrame,\n    subject: str,\n    events=None,\n    rate=\"30s\",\n    bp_pen=30,\n    pickle=False,\n) -> pd.DataFrame:\n    # load one subject\n    for s, d in ts.group_by(\"series_id\"):\n        if s == subject:\n            df = d\n            break\n\n    # convert anglez to -1, 1\n    df = df.with_columns(\n        anglez=df[\"anglez\"].map_elements(np.deg2rad) / np.pi * 2,\n        timestamp=df[\"timestamp\"].str.to_datetime(),\n        frac_day=_calc_frac_day(df),\n    )\n\n    # add anglez standard deviation (for break point detection)\n\n    wavelet = fcwt.Morlet(2.0)\n    hi_band = pl.Series(values=_fast_wavelet_transform(df, wavelet))\n\n    anglez_std = (\n        df.set_sorted([\"timestamp\"])\n        # .group_by_rolling(\n        .rolling(\n            index_column=\"timestamp\",\n            period=\"60s\",\n        )\n        .agg([pl.std(\"anglez\")])[\"anglez\"]\n    )\n\n    df = df.with_columns(\n        anglez_std=anglez_std.shift(6),\n        hi_band=hi_band,\n    ).fill_null(0)\n\n    df = (\n        df.set_sorted([\"timestamp\"])\n        .group_by_dynamic(\"timestamp\", every=rate)\n        .agg(\n            [\n                pl.col(\"enmo\").mean(),\n                pl.col(\"anglez\").median(),\n                pl.col(\"frac_day\").first(),\n                pl.col(\"step\").first(),\n                pl.col(\"hi_band\").mean(),\n                pl.col(\"anglez_std\").mean(),\n            ]\n        )\n    )\n\n    df = df.with_columns(\n        mp=pl.lit(_calc_matrix_profile(df, col=\"hi_band\")),\n        lids=_calc_lids(df),\n    )\n\n    break_points = _calc_break_points(df, cols=[\"anglez_std\", \"hi_band\"], bp_pen=30)\n\n    df = df.to_pandas().set_index(\"timestamp\")\n\n    if len(break_points) > 0:\n        df = df.join(break_points)\n        df.bp = df.bp.fillna(0).rolling(window=30, center=True).max()\n\n    # add events if given\n    if isinstance(events, pd.DataFrame):\n        df = _add_events(df, events, subject)\n\n    # save to pickle file\n    if pickle:\n        print(f\"saving: {subject}.pkl\")\n        df.to_pickle(f\"{subject}.pkl\")\n\n    return df\n\n\ndef calc_features(\n    df,\n    locs,\n    fc_params,\n    targets=None,\n    delta=\"1800s\",\n    cols=[\"hi_band\", \"mp\"],\n):\n    training = pd.DataFrame([])\n\n    for i in locs:\n        window = df.loc[\n            i - pd.Timedelta(delta) : i + pd.Timedelta(delta),\n            [\"lids\", \"hi_band\", \"anglez_std\", \"mp\", \"frac_day\", \"bp\", \"subj\"],\n        ].copy()\n        w_l = window.iloc[: len(window.index) // 2]\n        w_r = window.iloc[len(window.index) // 2 :]\n        feats = []\n        for c in cols:\n            for tag, w in zip(\n                (\"pre_half\", \"post_half\"),\n                (w_l, w_r),\n            ):\n                f = extract_features(\n                    w,\n                    column_value=c,\n                    column_id=\"subj\",\n                    default_fc_parameters=fc_params,\n                    disable_progressbar=True,\n                    n_jobs=1,\n                ).add_suffix(f\"_{delta}_{tag}\")\n                feats.append(f)\n\n        if targets is not None:\n            event = targets.loc[i]\n        else:\n            event = -1\n\n        feats = pd.concat(feats, axis=1).assign(\n            frac_day=df.loc[i].frac_day, event=event, step=df.loc[i].step\n        )\n\n        training = pd.concat([training, feats], axis=0)\n\n    training = training.reset_index().rename(columns={\"index\": \"subj\"})\n\n    return training","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:04.423923Z","iopub.execute_input":"2023-12-05T13:31:04.424261Z","iopub.status.idle":"2023-12-05T13:31:04.452743Z","shell.execute_reply.started":"2023-12-05T13:31:04.424235Z","shell.execute_reply":"2023-12-05T13:31:04.450953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = [\n    pd.read_pickle(d).assign(subj=d.stem)\n    for d in Path(\"/kaggle/input/data-prep-polars\").glob(\"*.pkl\")\n]\n\n# training = []\n\nfor df in data:\n    df.frac_day = df.frac_day.div(100).mul(2 * np.pi).apply(np.cos).add(1).div(2)\n    df.mp /= 5\n\n#     matched_night = df.night.fillna(0).rolling(window=2 * 60, center=True).max()\n#     locs = df[df.night > 0].index.tolist()\n#     bps = [i for i in df[df.bp.shift(15).diff() > 0].index if matched_night.loc[i] < 1]\n\n#     target = df.event.copy()\n#     target.loc[bps] = -1\n\n#     feat = calc_features(df, locs + bps, fc_params, targets=target)\n#     training.append(feat)\n\n# training = pd.concat(training)","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:04.453722Z","iopub.execute_input":"2023-12-05T13:31:04.454062Z","iopub.status.idle":"2023-12-05T13:31:20.902867Z","shell.execute_reply.started":"2023-12-05T13:31:04.454031Z","shell.execute_reply":"2023-12-05T13:31:20.901586Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# cv = GroupKFold(n_splits=2)\n# model = LGBMClassifier(is_unbalance=True)\n# drop_cols = [\"event\", \"subj\", \"step\"]\n\n# for t, v in cv.split(training, groups=training.subj):\n#     X = training.iloc[t]\n#     model.fit(\n#         X=X.fillna(0).drop(drop_cols, axis=1),\n#         y=X.event.astype(str),\n#     )\n\n#     Xv = training.iloc[v]\n#     y_prob = model.predict_proba(Xv.fillna(0).drop(drop_cols, axis=1))\n#     y_pred = y_prob.argmax(axis=1)\n#     y_true = Xv.event + 1\n#     print(confusion_matrix(y_pred=y_pred, y_true=y_true))\n\n# model.fit(X=training.fillna(0).drop(drop_cols, axis=1), y=training.event.astype(str))\n# model.booster_.save_model(\"lgb.model.txt\")\n\n# pd.Series(\n#     model.feature_importances_,\n#     index=training.drop(drop_cols, axis=1).columns[:],\n# ).sort_values().tail(100).plot(kind=\"barh\", figsize=(6, 16))","metadata":{"papermill":{"duration":144.007671,"end_time":"2023-11-03T00:46:57.087749","exception":false,"start_time":"2023-11-03T00:44:33.080078","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-12-05T13:31:20.904182Z","iopub.execute_input":"2023-12-05T13:31:20.904567Z","iopub.status.idle":"2023-12-05T13:31:20.910294Z","shell.execute_reply.started":"2023-12-05T13:31:20.904534Z","shell.execute_reply":"2023-12-05T13:31:20.909367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# model = Booster(model_file=\"/kaggle/input/cmi-breakpoint-lgbc/lgb.model.txt\")\n# drop_cols = [\"event\", \"subj\", \"step\"]\n# model.predict(calc_features(df, df[df.night > 0].index.tolist(), fc_params).drop(drop_cols, axis=1))","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:20.911180Z","iopub.execute_input":"2023-12-05T13:31:20.912216Z","iopub.status.idle":"2023-12-05T13:31:20.928098Z","shell.execute_reply.started":"2023-12-05T13:31:20.912168Z","shell.execute_reply":"2023-12-05T13:31:20.926762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"_____________________\n\n1D UNet segmentation (start/end points)","metadata":{}},{"cell_type":"code","source":"def dice_loss(inputv, target):\n    # AI generated code\n    # Flatten the inputv and target tensors\n    inputv = inputv.flatten()\n    target = target.flatten()\n\n    # Calculate the intersection and union of the input and target sets\n    intersection = (inputv * target).sum()\n    union = inputv.sum() + target.sum()\n\n    # Compute Dice loss\n    dice = (2.0 * intersection + 1e-12) / (union + 1e-12)\n\n    return 1 - dice\n\n\ndef jaccard_loss(inputv, target):\n    # AI generated code\n    # Flatten the inputv and target tensors\n    inputv = inputv.flatten()\n    target = target.flatten()\n\n    # Calculate the intersection and union of the input and target sets\n    intersection = (inputv * target).sum()\n    union = inputv.sum() + target.sum()\n\n    # Compute Jaccard loss\n    jaccard = (intersection + 1e-12) / (union - intersection + 1e-12)\n\n    return 1 - jaccard\n\n\ndef replace_segments(binary_list, value_to_replace, replacement_value, threshold):\n    # AI generated code\n    # Initialize variables\n    current_segment = []\n\n    # Iterate over binary list\n    for i, num in enumerate(binary_list):\n        # Check if number is the same as the value to replace\n        if num == value_to_replace:\n            current_segment.append(i)\n        else:\n            # If segment size is smaller than threshold, replace with replacement value\n            if len(current_segment) < threshold:\n                for index in current_segment:\n                    binary_list[index] = replacement_value\n            current_segment = []\n\n    # Handle last segment\n    if len(current_segment) < threshold:\n        for index in current_segment:\n            binary_list[index] = replacement_value\n\n    return binary_list\n\n\nclass sensor_dataset(Dataset):\n    def __init__(self, dfs, valid=False, window=2 * 60 * 24):\n        self.dfs = dfs\n        self.valid = valid\n        self.window = window\n\n    def __len__(self):\n        return len(self.dfs)\n\n    @staticmethod\n    def input_signal(df):\n        select_cols = [\n            \"lids\",\n            \"hi_band\",\n            \"anglez_std\",\n            \"mp\",\n        ]\n        dx = df[select_cols]\n        dx = dx.div(dx.sum(axis=1), axis=0).assign(\n            frac_day=df[\"frac_day\"], bp=(df[\"bp\"].fillna(0).shift(15) > 0).astype(\"int\")\n        )\n\n        signal = dx[\n            select_cols\n            + [\n                \"frac_day\",\n                \"bp\",\n            ]\n        ].values\n        return signal\n\n    def __getitem__(self, item):\n        df = self.dfs[item]\n        signal = self.input_signal(df)\n        mask = df[\"event\"].values\n        window = self.window\n\n        if self.valid and len(mask) <= self.window:\n            offset = 0\n        else:\n            offset = np.random.randint(len(mask) - window)\n\n        X = torch.from_numpy(signal[offset : offset + window, :]).permute(1, 0)\n        y = torch.from_numpy(mask[offset : offset + window]).unsqueeze(0)\n\n        return X, y\n\n\nclass UNet(nn.Module):\n    def __init__(self, p=0.5, in_channels=4):\n        super(UNet, self).__init__()\n\n        self.encoder1 = nn.Sequential(\n            nn.Conv1d(in_channels, 64, kernel_size=5, padding=2, dilation=1),\n            nn.BatchNorm1d(64),\n            nn.ReLU(inplace=True),\n            nn.Dropout(p),\n            nn.MaxPool1d(kernel_size=2, stride=2),\n        )\n\n        self.encoder2 = nn.Sequential(\n            nn.Conv1d(64, 128, kernel_size=5, padding=2, dilation=1),\n            nn.BatchNorm1d(128),\n            nn.ReLU(inplace=True),\n            nn.Dropout(p),\n            nn.MaxPool1d(kernel_size=2, stride=2),\n        )\n\n        self.encoder3 = nn.Sequential(\n            nn.Conv1d(128, 256, kernel_size=5, padding=2, dilation=1),\n            nn.BatchNorm1d(256),\n            nn.ReLU(inplace=True),\n            nn.Dropout(p),\n            nn.MaxPool1d(kernel_size=2, stride=2),\n        )\n\n        self.middle = nn.GRU(256, 256, batch_first=True)\n\n        self.decoder1 = nn.Sequential(\n            nn.Conv1d(512, 256, kernel_size=3, padding=1),\n            nn.ReLU(inplace=True),\n            nn.Dropout(p),\n            nn.ConvTranspose1d(256, 128, kernel_size=2, stride=2),\n        )\n\n        self.decoder2 = nn.Sequential(\n            nn.Conv1d(256, 128, kernel_size=3, padding=1),\n            nn.ReLU(inplace=True),\n            nn.Dropout(p),\n            nn.ConvTranspose1d(128, 64, kernel_size=2, stride=2),\n        )\n\n        self.decoder3 = nn.Sequential(\n            nn.Conv1d(128, 64, kernel_size=3, padding=1),\n            nn.ReLU(inplace=True),\n            nn.Dropout(p),\n            nn.ConvTranspose1d(64, 1, kernel_size=2, stride=2),\n            nn.Sigmoid(),\n        )\n\n    def forward(self, x):\n        enc1 = self.encoder1(x)\n        enc2 = self.encoder2(enc1)\n        enc3 = self.encoder3(enc2)\n\n        # GRU expects input of shape (batch_size, seq_len, num_features)\n        mid, _ = self.middle(enc3.permute(0, 2, 1))\n        # Permute back to (batch_size, num_features, seq_len)\n        mid = mid.permute(0, 2, 1)\n\n        dec1 = self.decoder1(torch.cat((mid, enc3), dim=1))\n        dec2 = self.decoder2(torch.cat((dec1, enc2), dim=1))\n        dec3 = self.decoder3(torch.cat((dec2, enc1), dim=1))\n        return dec3","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:20.929504Z","iopub.execute_input":"2023-12-05T13:31:20.929875Z","iopub.status.idle":"2023-12-05T13:31:20.954109Z","shell.execute_reply.started":"2023-12-05T13:31:20.929805Z","shell.execute_reply":"2023-12-05T13:31:20.952507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# device = \"cpu\"\n# total_epoch = 50\n# lr = 5e-4\n# batch_size = 4\n# window = 2 * 60 * 24\n\n# t, v = train_test_split(data, test_size=20)\n# v_all = [\n#     i for u in v for i in u.rolling(window=window, step=window) if len(i) == window\n# ]\n# tr_dataset = sensor_dataset(t, window=window)\n# vd_dataset = sensor_dataset(v_all, window=window, valid=True)\n\n# tr_loader = DataLoader(tr_dataset, batch_size=batch_size, shuffle=True, drop_last=True)\n# vd_loader = DataLoader(vd_dataset, batch_size=1)","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:20.956555Z","iopub.execute_input":"2023-12-05T13:31:20.957178Z","iopub.status.idle":"2023-12-05T13:31:20.971152Z","shell.execute_reply.started":"2023-12-05T13:31:20.957145Z","shell.execute_reply":"2023-12-05T13:31:20.969144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# n_channels = tr_dataset[0][0].shape[0]\n# model = UNet(p=0.5, in_channels=n_channels).double().to(device)\n\n# # optim = torch.optim.SGD(model.parameters(), lr=lr, weight_decay=1e-5)\n# optim = torch.optim.AdamW(model.parameters(), lr=lr, weight_decay=1e-4)\n# scheduler = ReduceLROnPlateau(optim, \"min\", patience=5, factor=0.8, min_lr=1e-8)\n\n# # loss_f = torch.nn.BCELoss()\n# # loss_f = jaccard_loss\n# loss_f = dice_loss","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:20.974503Z","iopub.execute_input":"2023-12-05T13:31:20.974953Z","iopub.status.idle":"2023-12-05T13:31:20.981010Z","shell.execute_reply.started":"2023-12-05T13:31:20.974923Z","shell.execute_reply":"2023-12-05T13:31:20.979585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# best_vloss = 1e10\n# for epoch in range(total_epoch):\n#     model.train()\n#     train_loss = 0\n#     tr_output_list = list()\n#     for idx, (X, y) in enumerate(tr_loader):\n#         optim.zero_grad()\n\n#         y_pred = model(X.to(device))\n#         loss = loss_f(y_pred, y.to(device))\n\n#         loss.backward()\n#         optim.step()\n\n#         train_loss += loss.cpu().item() / len(tr_loader)\n#         out_pred = F.softmax(y_pred, 1).detach().cpu().numpy().argmax(axis=1)\n#         tr_output_list.extend(list(out_pred))\n\n#     vd_loss = 0\n#     vd_jloss = 0\n#     for idx, (X, y) in enumerate(vd_loader):\n#         y_pred = model(X.to(device))\n#         out_pred = F.softmax(y_pred, 1).detach().cpu().numpy().argmax(axis=1)\n\n#         loss = loss_f(y_pred, y.to(device))\n#         vd_loss += loss.cpu().item() / len(vd_loader)\n\n#         jloss = jaccard_loss(y_pred, y.to(device))\n#         vd_jloss += jloss.cpu().item() / len(vd_loader)\n\n#     print(\n#         f\"Epoch {epoch} ----- Train loss: {train_loss:.3f}   Valid loss: {vd_loss:.3f} {vd_jloss:.3f}\"\n#     )\n\n#     if vd_loss < best_vloss:\n#         torch.jit.save(torch.jit.script(model), \"unet.model.pth\")\n#         best_vloss = vd_loss","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:20.982363Z","iopub.execute_input":"2023-12-05T13:31:20.982657Z","iopub.status.idle":"2023-12-05T13:31:20.992761Z","shell.execute_reply.started":"2023-12-05T13:31:20.982635Z","shell.execute_reply":"2023-12-05T13:31:20.991021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"-----------------","metadata":{}},{"cell_type":"code","source":"# feat\n# calc_features(vdf, locs=locs, fc_params=fc_params)\n# submission","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:25.486780Z","iopub.execute_input":"2023-12-05T13:31:25.487205Z","iopub.status.idle":"2023-12-05T13:31:25.492327Z","shell.execute_reply.started":"2023-12-05T13:31:25.487177Z","shell.execute_reply":"2023-12-05T13:31:25.490962Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# vdf = df\n# with torch.no_grad():\n#     preds = []\n#     for x in vdf.rolling(window=window, step=window // 2):\n#         if not len(x) == window:\n#             continue\n\n#         X = torch.from_numpy(sensor_dataset.input_signal(x)).permute(1, 0)[None, :]\n#         y_hat = unet_model(X).squeeze().numpy() > 0.5\n#         x = x.assign(y=y_hat)\n#         x.loc[\n#             x[\"mp\"].rolling(window=60, center=True).max().fillna(0) < 0.01, \"y\"\n#         ] = 0\n\n#         preds.append(x)\n\n#     preds = pd.concat(preds).reset_index().groupby(\"timestamp\").max().sort_index()\n#     preds[\"y\"] = replace_segments(\n#         replace_segments(\n#             preds[\"y\"].astype(int),\n#             0,\n#             1,\n#             threshold=2 * 15,  # remove awake periods smaller than 15 mins\n#         ),\n#         1,\n#         0,\n#         threshold=2 * 180,  # remove sleep periods smaller than 3h\n#     )\n\n#     drop_cols = [\"event\", \"subj\", \"step\"]\n\n\n#     locs = preds.index[preds.y.diff().abs() > 0].tolist()\n#     feat = calc_features(vdf, locs=locs, fc_params=fc_params)\n#     pred_prob = lgbc_model.predict(feat.drop(drop_cols, axis=1))[:, 1:]\n\n\n#     pred_event = (preds.y.diff().loc[locs] > 0).astype(\"int\")\n#     feat = feat.assign(event=pred_event, score=pred_prob.max(axis=1))\n\n#     submission = (\n#         feat.assign(event=pred_prob.argmax(axis=1), score=pred_prob.max(axis=1))\n#         .loc[:, [\"subj\", \"step\", \"event\", \"score\"]]\n#         .replace({\"event\": {1: \"onset\", 0: \"wakeup\"}})\n#         .rename(columns={\"subj\": \"series_id\"})\n#         .sort_values(by=['subj', 'step'])\n#     )","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:33.012537Z","iopub.execute_input":"2023-12-05T13:31:33.012947Z","iopub.status.idle":"2023-12-05T13:31:33.018810Z","shell.execute_reply.started":"2023-12-05T13:31:33.012916Z","shell.execute_reply":"2023-12-05T13:31:33.018016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def df_to_submit(vdf, unet_model, lgbc_model, fc_params, window=2 * 60 * 24, bp_only=False):\n    if bp_only:\n        locs = vdf.index[vdf.bp > 0].tolist()\n    else:\n        with torch.no_grad():\n            preds = []\n            for x in vdf.rolling(window=window, step=window // 2):\n                if not len(x) == window:\n                    continue\n\n                X = torch.from_numpy(sensor_dataset.input_signal(x)).permute(1, 0)[None, :]\n                y_hat = unet_model(X).squeeze().numpy() > 0.5\n                x = x.assign(y=y_hat)\n                x.loc[\n                    x[\"mp\"].rolling(window=60, center=True).max().fillna(0) < 0.01, \"y\"\n                ] = 0\n\n                preds.append(x)\n\n            preds = pd.concat(preds).reset_index().groupby(\"timestamp\").max().sort_index()\n            preds[\"y\"] = replace_segments(\n                replace_segments(\n                    preds[\"y\"].astype(int),\n                    0,\n                    1,\n                    threshold=2 * 15,  # remove awake periods smaller than 15 mins\n                ),\n                1,\n                0,\n                threshold=2 * 180,  # remove sleep periods smaller than 3h\n            )\n            locs = preds.index[preds.y.diff().abs() > 0].tolist()\n\n    feat = calc_features(vdf, locs=locs, fc_params=fc_params)\n    pred_prob = lgbc_model.predict(feat.drop([\"event\", \"subj\", \"step\"], axis=1))[:, 1:]\n    feat = feat.assign(event=pred_prob.argmax(axis=1), score=pred_prob.max(axis=1))\n\n    submission = (\n        feat.assign(event=pred_prob.argmax(axis=1), score=pred_prob.max(axis=1))\n        .loc[:, [\"subj\", \"step\", \"event\", \"score\"]]\n        .replace({\"event\": {1: \"onset\", 0: \"wakeup\"}})\n        .rename(columns={\"subj\": \"series_id\"})\n    )\n\n    return submission","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:34.779409Z","iopub.execute_input":"2023-12-05T13:31:34.779728Z","iopub.status.idle":"2023-12-05T13:31:34.790215Z","shell.execute_reply.started":"2023-12-05T13:31:34.779705Z","shell.execute_reply":"2023-12-05T13:31:34.789265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_ts = pl.read_parquet(\n    \"/kaggle/input/child-mind-institute-detect-sleep-states/train_series.parquet\"\n#     \"/kaggle/input/child-mind-institute-detect-sleep-states/test_series.parquet\"\n)\nlgbc_model = Booster(model_file=\"/kaggle/input/cmi-breakpoint-lgbc/lgb.model.txt\")\nunet_model = torch.jit.load(\"/kaggle/input/cmi-1d-unet/unet.model.pth\", map_location='cpu')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"window=2 * 60 * 24\nbreakpoint_only = True\n\nsubmissions = []\nfor i, subj in enumerate(test_ts[\"series_id\"].unique()):\n    df = get_one_subject(test_ts, subject=subj, rate=\"30s\").assign(subj=subj)\n\n    df.frac_day = df.frac_day.div(100).mul(2 * np.pi).apply(np.cos).add(1).div(2)\n    df.mp /= 5\n    df.bp = (df.bp.shift(15).diff() > 0).astype(\"int\")\n    df.loc[df.mp.rolling(window=30, center=True).max() < 0.002, 'bp'] = 0\n    \n    \n    try:\n        submissions.append(df_to_submit(df, unet_model, lgbc_model, fc_params, window, bp_only=breakpoint_only))\n    except:\n        print(f\"error processing {subj}\")\n        pass\n    \n","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:35:37.701007Z","iopub.execute_input":"2023-12-05T13:35:37.701389Z","iopub.status.idle":"2023-12-05T13:37:03.131738Z","shell.execute_reply.started":"2023-12-05T13:35:37.701358Z","shell.execute_reply":"2023-12-05T13:37:03.130589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submissions = (\n    pd.concat(submissions, ignore_index=True)\n    .reset_index()\n    .rename(columns={\"index\": \"row_id\"})\n).query('score > 0.01')\nsubmissions.to_csv(\"submission.csv\", index=False)\nsubmissions","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:37:03.133635Z","iopub.execute_input":"2023-12-05T13:37:03.133998Z","iopub.status.idle":"2023-12-05T13:37:03.155340Z","shell.execute_reply.started":"2023-12-05T13:37:03.133969Z","shell.execute_reply":"2023-12-05T13:37:03.154457Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# score(\n#     solution=solutions,\n#     submission=submissions,\n#     tolerances={\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#     series_id_column_name=\"series_id\",\n#     time_column_name=\"step\",\n#     event_column_name=\"event\",\n#     score_column_name=\"score\",\n# )","metadata":{"execution":{"iopub.status.busy":"2023-12-05T13:31:21.440479Z","iopub.status.idle":"2023-12-05T13:31:21.441130Z","shell.execute_reply.started":"2023-12-05T13:31:21.440971Z","shell.execute_reply":"2023-12-05T13:31:21.440986Z"},"trusted":true},"execution_count":null,"outputs":[]}]}