{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.12"},"papermill":{"default_parameters":{},"duration":333.993192,"end_time":"2023-09-10T17:40:59.325820","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2023-09-10T17:35:25.332628","version":"2.4.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 💤 LightGBM + Feature Implementation from Academic Paper","metadata":{}},{"cell_type":"markdown","source":"This notebook presents modeling with LightGBM classifier, with features adapted from [1]. It scored 19%++ score in public LB.\n\nSome of the features implemented:\n- Median absolute deviation\n- Locomotor Inactivity During Sleep (LIDS)\n\nI have also added plots for feature importances and prediction signals.\n\nHeavily modified from impressive work of [Zzzs: Random Forest model starter](https://www.kaggle.com/code/carlmcbrideellis/zzzs-random-forest-model-starter).\n\n[1] [*Sleep classification from wrist-worn accelerometer data using random forests*](https://www.nature.com/articles/s41598-020-79217-x.pdf)","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport lightgbm as lgbm\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport gc\n\nfrom time import time\nfrom itertools import groupby, product\nfrom scipy.signal import savgol_filter","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":4.36326,"end_time":"2023-09-10T17:35:34.531335","exception":false,"start_time":"2023-09-10T17:35:30.168075","status":"completed"},"tags":[],"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-09-14T15:59:26.578861Z","iopub.execute_input":"2023-09-14T15:59:26.579170Z","iopub.status.idle":"2023-09-14T15:59:27.986078Z","shell.execute_reply.started":"2023-09-14T15:59:26.579145Z","shell.execute_reply":"2023-09-14T15:59:27.985211Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loading Data","metadata":{}},{"cell_type":"code","source":"FULL, MINI = 'full', 'mini'\n\n# adjust to FULL or MINI\n# FULL = whole data, MINI = only 3 series with least steps\nMODE = FULL\n\nseries_with_least_steps = ['349c5562ee2c', 'c535634d7dcd', 'a9e5f5314bcb']","metadata":{"execution":{"iopub.status.busy":"2023-09-14T15:59:27.987737Z","iopub.execute_input":"2023-09-14T15:59:27.987969Z","iopub.status.idle":"2023-09-14T15:59:27.993336Z","shell.execute_reply.started":"2023-09-14T15:59:27.987948Z","shell.execute_reply":"2023-09-14T15:59:27.991785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"st = time()\n\ntrain = pd.read_parquet(\"/kaggle/input/zzzs-lightweight-training-dataset-target/Zzzs_train.parquet\")\nevents  = pd.read_csv(\"/kaggle/input/child-mind-institute-detect-sleep-states/train_events.csv\")\n\nif MODE == MINI:\n    train = train[train.series_id.isin(series_with_least_steps)]\n    events = events[events.series_id.isin(series_with_least_steps)]\n\n# the bottleneck for the running time is in converting to UTC timestamp\n# FULL: ~100s, MINI: ~10s\ntrain[\"timestamp\"] = pd.to_datetime(train[\"timestamp\"]).apply(lambda t: t.tz_localize(None))\nevents[\"timestamp\"] = pd.to_datetime(events[\"timestamp\"]).apply(lambda t: t.tz_localize(None))\n\nen = time()\nprint(f'Done loading in {int(en - st)}s')\n\nprint(train.shape)\ntrain.head()","metadata":{"_kg_hide-input":true,"papermill":{"duration":115.808222,"end_time":"2023-09-10T17:37:30.346411","exception":false,"start_time":"2023-09-10T17:35:34.538189","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-09-14T15:59:27.994360Z","iopub.execute_input":"2023-09-14T15:59:27.994646Z","iopub.status.idle":"2023-09-14T16:01:46.404662Z","shell.execute_reply.started":"2023-09-14T15:59:27.994621Z","shell.execute_reply":"2023-09-14T16:01:46.403625Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Helper Functions","metadata":{}},{"cell_type":"code","source":"# Credits: https://www.kaggle.com/code/patrick0302/viz-of-sleeping-time-series\ndef plot_series_with_events(series_list, timestamp, events, events_pred=None, min_vals=None, ylim=None):\n    series_id = events.series_id.iloc[0]\n    if min_vals is None:\n        min_vals = [0] * len(series_list)\n    \n    # Filter event data based on the series_id\n    sample_onset = events.loc[events['event'] == 'onset', 'timestamp'].dropna()\n    sample_wakeup = events.loc[events['event'] == 'wakeup', 'timestamp'].dropna()\n    \n    if events_pred is not None:\n        onset_pred = events_pred.loc[events_pred['event'] == 'onset', 'timestamp'].dropna()\n        wakeup_pred = events_pred.loc[events_pred['event'] == 'wakeup', 'timestamp'].dropna()\n    \n    # Plot enmo and anglez\n    plt.figure(figsize=(15, 3))\n\n    for series, min_val in zip(series_list, min_vals):\n        label = series.name\n        \n        mi, ma =  series.min(), series.max()\n        if min_val == 0:\n            plt.plot(timestamp, (series - mi) / (ma - mi), label=label,\n                     linewidth=1)\n        else:\n            plt.plot(timestamp, 2 * (series - mi) / (ma - mi) - 1, label=label,\n                     linewidth=1)\n\n        for onset, wakeup in zip(sample_onset, sample_wakeup):\n            plt.fill_between(timestamp, min_val, 1, zorder=2,\n                             where=(timestamp >= onset) & (timestamp <= wakeup),\n                             color='red', alpha=0.2, label='Event Area')\n\n    if events_pred is not None:\n        for onset, wakeup in zip(onset_pred, wakeup_pred):\n            plt.fill_between(timestamp, min_val, 1, zorder=2,\n                             where=(timestamp >= onset) & (timestamp <= wakeup),\n                             color='yellow', alpha=0.2, label='Pred Event Area')\n\n    handles, labels = plt.gca().get_legend_handles_labels()\n    new_labels, new_handles = [], []\n    for handle, label_ in zip(handles, labels):\n        if label_ not in new_labels:\n            new_handles.append(handle)\n            new_labels.append(label_)\n            \n    if ylim != None:\n        plt.ylim(ylim)\n\n    plt.legend(new_handles, new_labels)\n    plt.xlabel('Timestamp')\n    plt.ylabel(f'Normalized value')\n    plt.title(f'Normalized data over time with event flags ({series_id})')\n    plt.show()","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2023-09-14T16:01:46.405873Z","iopub.execute_input":"2023-09-14T16:01:46.406140Z","iopub.status.idle":"2023-09-14T16:01:46.417788Z","shell.execute_reply.started":"2023-09-14T16:01:46.406118Z","shell.execute_reply":"2023-09-14T16:01:46.416775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature Engineering","metadata":{}},{"cell_type":"code","source":"signal_awake = dict(zip(range(1440), np.sin(np.linspace(0, np.pi, 1440) + 0.208 * np.pi) ** 24))\nsignal_onset = dict(zip(range(1440), np.sin(np.linspace(0, np.pi, 1440) + 0.555 * np.pi) ** 24))\n\ndef diff_feature(df: pd.Series, col: str, window_size: int = 60):\n    ret = df.groupby('series_id')[col].diff(periods=window_size)\n    ret = ret.fillna(ret.loc[ret.first_valid_index()])\n    gc.collect()\n    return ret\n\n\ndef rolling_feature(df: pd.Series, col: str, agg: str, window_size: int = 60):\n    ret = df.groupby('series_id')[col].rolling(window_size, center=True).agg(agg)\n    ret = ret.fillna(method='bfill').fillna(method='ffill')\n    gc.collect()\n    return ret\n\n\ndef create_features(df_, window_size=60, verbose=False):\n    start_time = time()\n    df = df_.copy()\n    \n    if 'awake' in df.columns:\n        df.awake = df.awake.astype(np.uint8)\n    \n    df['weekday'] = df['timestamp'].dt.weekday.astype(np.int8)\n    df['signal_onset'] = (df.timestamp.dt.hour * 60 + df.timestamp.dt.minute).map(signal_onset).astype(np.float16)\n    df['signal_awake'] = (df.timestamp.dt.hour * 60 + df.timestamp.dt.minute).map(signal_awake).astype(np.float16)\n    \n    df['anglez'] = abs(df['anglez'])\n    \n    df['enmo_rolling_mean']     = rolling_feature(df, col='enmo', agg='mean', window_size=window_size).values\n#     df['enmo_rolling_min']      = rolling_feature(df, col='enmo', agg='min', window_size=window_size).astype('float16').values\n    df['enmo_rolling_max']      = rolling_feature(df, col='enmo', agg='max', window_size=window_size).values\n    df['enmo_rolling_std']      = rolling_feature(df, col='enmo', agg='std', window_size=window_size).astype('float16').values\n    df['enmo_rolling_median']   = rolling_feature(df, col='enmo', agg='median', window_size=window_size).values\n    \n    df['anglez_rolling_mean']     = rolling_feature(df, col='anglez', agg='mean', window_size=window_size).values\n#     df['anglez_rolling_min']      = rolling_feature(df, col='anglez', agg='min', window_size=window_size).astype('float16').values\n    df['anglez_rolling_max']      = rolling_feature(df, col='anglez', agg='max', window_size=window_size).values\n    df['anglez_rolling_std']      = rolling_feature(df, col='anglez', agg='std', window_size=window_size).astype('float16').values\n    df['anglez_rolling_median']   = rolling_feature(df, col='anglez', agg='median', window_size=window_size).values\n\n    df['enmo_prev_diff_rolling_mean']  = diff_feature(df, col='enmo_rolling_mean', window_size=window_size).astype('float16').values\n    df['enmo_prev_diff_rolling_max']   = diff_feature(df, col='enmo_rolling_max', window_size=window_size).astype('float16').values\n    \n    df['anglez_prev_diff_rolling_mean']  = diff_feature(df, col='anglez_rolling_mean', window_size=window_size).astype('float16').values\n    df['anglez_prev_diff_rolling_max']   = diff_feature(df, col='anglez_rolling_max', window_size=window_size).astype('float16').values\n\n    gc.collect()\n    \n    x = (df.anglez - df.anglez_rolling_median).abs().to_frame(name='ad').assign(series_id=df.series_id)\n    x = x.groupby('series_id')['ad'].rolling(window_size, center=True).agg('median')\n    df['anglez_mad'] = x.fillna(method='bfill').fillna(method='ffill').astype('float16').values\n    \n    x = (df.enmo - df.enmo_rolling_median).abs().to_frame(name='ad').assign(series_id=df.series_id)\n    x = x.groupby('series_id')['ad'].rolling(window_size, center=True).agg('median')\n    df['enmo_mad'] = x.fillna(method='bfill').fillna(method='ffill').astype('float16').values\n\n    z = np.maximum(0., df.enmo - 0.02)\n    z = z.to_frame(name='lids').assign(series_id=df.series_id)\n    z = rolling_feature(z, col='lids', agg='sum', window_size=120)\n    z = 100 / (z + 1)\n    z = z.to_frame().reset_index()\n    z = rolling_feature(z, col='lids', agg='mean', window_size=360)\n    df['lids'] = z.astype('float16').values\n    \n    del (df['enmo_rolling_median'], df['anglez_rolling_median'])\n    \n    cols_astype_float16 = df.columns[df.columns.str.contains('enmo|anglez|signal|lids')].tolist()\n    for col in cols_astype_float16:\n        df[col] = df[col].astype('float16')\n    \n    gc.collect()    \n    \n    end_time = time()\n    eta = int(end_time - start_time)\n    if verbose:\n        print(f'Done in: {eta}s')\n    \n    return df","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.024916,"end_time":"2023-09-10T17:37:30.377942","exception":false,"start_time":"2023-09-10T17:37:30.353026","status":"completed"},"tags":[],"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2023-09-14T16:01:46.421450Z","iopub.execute_input":"2023-09-14T16:01:46.421820Z","iopub.status.idle":"2023-09-14T16:01:46.444082Z","shell.execute_reply.started":"2023-09-14T16:01:46.421787Z","shell.execute_reply":"2023-09-14T16:01:46.442614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"window_size = 36\n# 1 step = 5s, 10 steps = 50s, 20 steps = 100s\n# 1 min = 12 steps, 3 mins = 36 steps, 5 mins = 60 steps\n\ntrain = create_features(train, window_size=window_size, verbose=True)\ntrain.info()","metadata":{"papermill":{"duration":21.620544,"end_time":"2023-09-10T17:37:52.004913","exception":false,"start_time":"2023-09-10T17:37:30.384369","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-09-14T16:01:46.445406Z","iopub.execute_input":"2023-09-14T16:01:46.446458Z","iopub.status.idle":"2023-09-14T16:03:08.795562Z","shell.execute_reply.started":"2023-09-14T16:01:46.446396Z","shell.execute_reply":"2023-09-14T16:03:08.794684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features = train.columns[train.columns.str.contains('enmo|anglez|signal|lids|weekday')].tolist()\nprint(len(features))","metadata":{"execution":{"iopub.status.busy":"2023-09-14T16:03:08.796650Z","iopub.execute_input":"2023-09-14T16:03:08.797782Z","iopub.status.idle":"2023-09-14T16:03:08.804013Z","shell.execute_reply.started":"2023-09-14T16:03:08.797751Z","shell.execute_reply":"2023-09-14T16:03:08.802620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling","metadata":{}},{"cell_type":"code","source":"classifier = lgbm.LGBMClassifier(n_estimators=100,\n                                 min_samples_leaf=300,\n                                 random_state=42)\n\nst = time()\nclassifier.fit(train[features], train[\"awake\"])\nen = time()\n\neta = int(en - st)\nprint(f'Done training in {eta}s')","metadata":{"papermill":{"duration":63.139012,"end_time":"2023-09-10T17:38:55.150425","exception":false,"start_time":"2023-09-10T17:37:52.011413","status":"completed"},"tags":[],"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2023-09-14T16:03:08.805740Z","iopub.execute_input":"2023-09-14T16:03:08.806053Z","iopub.status.idle":"2023-09-14T16:04:58.175478Z","shell.execute_reply.started":"2023-09-14T16:03:08.806029Z","shell.execute_reply":"2023-09-14T16:04:58.174582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feat_imp = pd.Series(classifier.feature_importances_, index=features).sort_values()\nplt.figure(figsize=(6, 6))\nplt.barh(feat_imp.index, feat_imp.values)\nplt.grid(axis='x')\nplt.title('Feature Importance')\nplt.gca().set_axisbelow(True)\nplt.show()","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2023-09-14T16:04:58.176494Z","iopub.execute_input":"2023-09-14T16:04:58.176898Z","iopub.status.idle":"2023-09-14T16:04:58.441728Z","shell.execute_reply.started":"2023-09-14T16:04:58.176875Z","shell.execute_reply":"2023-09-14T16:04:58.440733Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot Predictions","metadata":{}},{"cell_type":"code","source":"# https://stackoverflow.com/questions/73777727/how-to-mark-start-end-of-a-series-of-non-null-and-non-0-values-in-a-column-of-a\ndef get_event(df):\n    lstCV = zip(df.series_id, df.label)\n    lstPOI = []\n    for (c, v), g in groupby(lstCV, lambda cv: \n                            (cv[0], cv[1]!=0 and not pd.isnull(cv[1]))):\n        llg = sum(1 for item in g)\n        if v is False: \n            lstPOI.extend([0]*llg)\n        else: \n            lstPOI.extend(['onset']+(llg-2)*[0]+['wakeup'] if llg > 1 else [0])\n    return lstPOI\n\n\ndef predict_events(classifier, test, threshold=0.5, smoothing_length=None):\n    test[\"score\"] = classifier.predict_proba(test.drop('series_id', axis=1))[:, 0]\n\n    # smoothing of the predictions\n    if smoothing_length != None:\n        test[\"smoothed_score\"] = test.groupby('series_id')[\"score\"] \\\n            .rolling(smoothing_length, center=True) \\\n            .mean().fillna(method=\"bfill\") \\\n            .fillna(method=\"ffill\") \\\n            .values\n\n    test[\"label\"] = 1 * (test.smoothed_score >= threshold)\n    test[\"event\"] = get_event(test)\n    return test","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.033982,"end_time":"2023-09-10T17:38:55.191572","exception":false,"start_time":"2023-09-10T17:38:55.157590","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-09-14T16:04:58.442592Z","iopub.execute_input":"2023-09-14T16:04:58.442911Z","iopub.status.idle":"2023-09-14T16:04:58.453120Z","shell.execute_reply.started":"2023-09-14T16:04:58.442890Z","shell.execute_reply":"2023-09-14T16:04:58.452285Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"smoothing_length = 720\n\nt1 = time()\n\ntrain_event_preds = predict_events(\n    classifier,\n    train[['series_id'] + features],\n    smoothing_length=smoothing_length\n)\ntrain_event_preds['timestamp'] = train['timestamp']\n\nt2 = time()\neta = int(t2 - t1)\nprint(f'Done predicting in {eta}s')\n\nfor chosen_series_id in series_with_least_steps:\n    sample = train[train.series_id == chosen_series_id]\n    sample_events = events[events.series_id == chosen_series_id]\n    sample_event_preds = train_event_preds.loc[\n        train_event_preds.series_id == chosen_series_id,\n        ['timestamp', 'series_id', 'event', 'score', 'smoothed_score']\n    ]\n\n    score  = sample_event_preds.score\n    sscore = sample_event_preds.smoothed_score\n    sample_event_preds.drop(['score', 'smoothed_score'], axis=1, inplace=True)\n\n    t1 = time()\n    plot_series_with_events(\n        [sample.lids, score, sscore],\n        sample.timestamp,\n        events=sample_events,\n        events_pred=sample_event_preds,\n        min_vals=[0., 0., 0.]\n    )\n    t2 = time()\n    eta = int(t2 - t1)\n    print(f'Done plot in {eta}s')\n\n# del(train, events) # turn this on to reduce memory usage during submission\ngc.collect()\n# del(sample, train_event_preds, score, sscore, sample_events, sample_event_preds)","metadata":{"_kg_hide-input":true,"papermill":{"duration":122.782803,"end_time":"2023-09-10T17:40:57.981667","exception":false,"start_time":"2023-09-10T17:38:55.198864","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-09-14T16:04:58.454448Z","iopub.execute_input":"2023-09-14T16:04:58.454784Z","iopub.status.idle":"2023-09-14T16:05:40.917867Z","shell.execute_reply.started":"2023-09-14T16:04:58.454754Z","shell.execute_reply":"2023-09-14T16:05:40.916920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Generate Submission","metadata":{}},{"cell_type":"code","source":"test  = pd.read_parquet(\"/kaggle/input/child-mind-institute-detect-sleep-states/test_series.parquet\")\ntest[\"timestamp\"] = pd.to_datetime(test[\"timestamp\"]).apply(lambda t: t.tz_localize(None))\n\ntest = create_features(test, window_size=window_size, verbose=True)\n\ntest.loc[:, ['score', 'event']] = predict_events(\n    classifier,\n    test[['series_id'] + features],\n    smoothing_length=smoothing_length\n).loc[:, ['score', 'event']]\n\nsample_submission = test.loc[test[\"event\"] != 0][[\"series_id\",\"step\",\"event\",\"score\"]].copy().reset_index(drop=True).reset_index(names=\"row_id\")\nsample_submission.to_csv('submission.csv', index=False)\nsample_submission","metadata":{"papermill":{"duration":0.083231,"end_time":"2023-09-10T17:40:58.073470","exception":false,"start_time":"2023-09-10T17:40:57.990239","status":"completed"},"tags":[],"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-09-14T16:05:40.918981Z","iopub.execute_input":"2023-09-14T16:05:40.919638Z","iopub.status.idle":"2023-09-14T16:05:42.452563Z","shell.execute_reply.started":"2023-09-14T16:05:40.919613Z","shell.execute_reply":"2023-09-14T16:05:42.451413Z"},"trusted":true},"execution_count":null,"outputs":[]}],"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}}