{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":53666,"databundleVersionId":6589269,"sourceType":"competition"},{"sourceId":151583600,"sourceType":"kernelVersion"}],"dockerImageVersionId":30558,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Feature Engineering and Random Forest Prediction to Detect Sleep States\nIn this notebook I... \n- Use the Polars library to load and transform the CMI dataset and incorporate features inspired by the work done in my [Sleep Data Exploration](https://www.kaggle.com/code/lccburk/sleep-data-exploration) notebook. \n- Import and implement the [Event Detection AP](https://www.kaggle.com/code/metric/event-detection-ap/notebook) score function to validate results prior to submission.\n- Define helper functions to create training sets and formatted submissions based on classifier results.\n- Train a Random Forest classifier (as well as gradient boost classifier) and use the above work to validate and create a submission.","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport polars as pl\nimport datetime \nfrom tqdm import tqdm\n\nimport plotly.express as px\nfrom plotly.subplots import make_subplots\nimport plotly.graph_objects as go\n\nfrom metric import score # Import event detection ap score function\n\n# These are variables to be used by the score function\ncolumn_names = {\n 'series_id_column_name': 'series_id',\n 'time_column_name': 'step',\n 'event_column_name': 'event',\n 'score_column_name': 'score',\n}\n\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}","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-10-03T14:36:05.062823Z","iopub.execute_input":"2023-10-03T14:36:05.063149Z","iopub.status.idle":"2023-10-03T14:36:06.41494Z","shell.execute_reply.started":"2023-10-03T14:36:05.063124Z","shell.execute_reply":"2023-10-03T14:36:06.414033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Importing data","metadata":{}},{"cell_type":"code","source":"# Importing data \n\n# Column transformations\n\ndt_transforms = [\n pl.col('timestamp').str.to_datetime(), \n (pl.col('timestamp').str.to_datetime().dt.year()-2000).cast(pl.UInt8).alias('year'), \n pl.col('timestamp').str.to_datetime().dt.month().cast(pl.UInt8).alias('month'),\n pl.col('timestamp').str.to_datetime().dt.day().cast(pl.UInt8).alias('day'), \n pl.col('timestamp').str.to_datetime().dt.hour().cast(pl.UInt8).alias('hour')\n]\n\ndata_transforms = [\n pl.col('anglez').cast(pl.Int16), # Casting anglez to 16 bit integer\n (pl.col('enmo')*1000).cast(pl.UInt16), # Convert enmo to 16 bit uint\n]\n\ntrain_series = pl.scan_parquet('/kaggle/input/child-mind-institute-detect-sleep-states/train_series.parquet').with_columns(\n dt_transforms + data_transforms\n )\n\ntrain_events = pl.read_csv('/kaggle/input/child-mind-institute-detect-sleep-states/train_events.csv').with_columns(\n dt_transforms\n ).drop_nulls()\n\ntest_series = pl.scan_parquet('/kaggle/input/child-mind-institute-detect-sleep-states/test_series.parquet').with_columns(\n dt_transforms + data_transforms\n )\n\n# Removing null events and nights with mismatched counts from series_events\nmismatches = train_events.drop_nulls().group_by(['series_id', 'night']).agg([\n ((pl.col('event') == 'onset').sum() == (pl.col('event') == 'wakeup').sum()).alias('balanced')\n ]).sort(by=['series_id', 'night']).filter(~pl.col('balanced'))\n\nfor mm in mismatches.to_numpy(): \n train_events = train_events.filter(~((pl.col('series_id') == mm[0]) & (pl.col('night') == mm[1])))\n\n# Getting series ids as a list for convenience\nseries_ids = train_events['series_id'].unique(maintain_order=True).to_list()\n\n# Updating train_series to only keep these series ids\ntrain_series = train_series.filter(pl.col('series_id').is_in(series_ids))","metadata":{"execution":{"iopub.status.busy":"2023-10-05T03:27:42.390305Z","iopub.execute_input":"2023-10-05T03:27:42.390609Z","iopub.status.idle":"2023-10-05T03:27:42.427921Z","shell.execute_reply.started":"2023-10-05T03:27:42.390587Z","shell.execute_reply":"2023-10-05T03:27:42.427257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Feature Engineering\n\nThe features in this model will consist of\n- The current`hour`.\n- Rolling aggregates (mean, max, std) of `anglez` and `enmo` over a variety of window sizes, from 5 minutes to 8 hours. \n- Rolling aggregates (mean, max, std) of `anglez` and `enmo` [**total variation**](https://en.wikipedia.org/wiki/Total_variation) (or *first variation*, i.e. 1v) over a variety of window sizes, from 5 minutes to 8 hours.\n\n#### Motivation for using total variation \n\nFrom my earlier [Data Exploration](https://www.kaggle.com/code/lccburk/sleep-data-exploration) we observe that during sleeping periods `anglez` resembles a pure [jump process](https://en.wikipedia.org/wiki/Jump_process), while during wakeful periods it resembles a [diffusion process](https://en.wikipedia.org/wiki/Diffusion_process), as shown in the following data sample:\n\n\n\nImportantly, jump and diffusion processes can be distinguished by what is known as their [**total variation**](https://en.wikipedia.org/wiki/Total_variation) - essentially, the sum total of the absolute differences between the points. For diffusion processes, which jiggle around constantly, the total variation is infinite, while for jump processes, which only change by finite amounts a countable number of times, the total variation is finite.\n\n\n\nFor non-continuous, evenly sampled functions such as our time series, the total variation of a function $f(t)$ on an interval $[a,b]$ can be defined simply as\n\n$$V_a^b(f) := \\sum_{j=0}^{n-1} |f(t_{j+1}) - f(t_j)|$$\n\nwhere $t_0=a$, $t_n=b$, and $\\forall j: t_{j+1}-t_j = \\frac{a-b}{n}$. This can be calculated efficiently using Polars' built in `.diff()` and `.abs()` functions. \n\nThese features will give the classification model information which characterizes the stochastic behavior of the variable in the recent past, which can be much more useful for classifying sleep state than the variable value itself.","metadata":{}},{"cell_type":"code","source":"features, feature_cols = [pl.col('hour')], ['hour']\n\nfor mins in [5, 30, 60*2, 60*8] :\n \n for var in ['enmo', 'anglez'] :\n \n features += [\n pl.col(var).rolling_mean(12 * mins, center=True, min_periods=1).abs().cast(pl.UInt16).alias(f'{var}_{mins}m_mean'),\n pl.col(var).rolling_max(12 * mins, center=True, min_periods=1).abs().cast(pl.UInt16).alias(f'{var}_{mins}m_max'),\n pl.col(var).rolling_std(12 * mins, center=True, min_periods=1).abs().cast(pl.UInt16).alias(f'{var}_{mins}m_std')\n ]\n\n feature_cols += [ \n f'{var}_{mins}m_mean', f'{var}_{mins}m_max', f'{var}_{mins}m_std'\n ]\n\n # Getting first variations\n features += [\n (pl.col(var).diff().abs().rolling_mean(12 * mins, center=True, min_periods=1)*10).abs().cast(pl.UInt32).alias(f'{var}_1v_{mins}m_mean'),\n (pl.col(var).diff().abs().rolling_max(12 * mins, center=True, min_periods=1)*10).abs().cast(pl.UInt32).alias(f'{var}_1v_{mins}m_max'),\n (pl.col(var).diff().abs().rolling_std(12 * mins, center=True, min_periods=1)*10).abs().cast(pl.UInt32).alias(f'{var}_1v_{mins}m_std')\n ]\n\n feature_cols += [ \n f'{var}_1v_{mins}m_mean', f'{var}_1v_{mins}m_max', f'{var}_1v_{mins}m_std'\n ]\n\nid_cols = ['series_id', 'step', 'timestamp']\n\ntrain_series = train_series.with_columns(\n features\n).select(id_cols + feature_cols)\n\ntest_series = test_series.with_columns(\n features\n).select(id_cols + feature_cols)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T03:29:10.719516Z","iopub.execute_input":"2023-10-05T03:29:10.719845Z","iopub.status.idle":"2023-10-05T03:29:10.730222Z","shell.execute_reply.started":"2023-10-05T03:29:10.719804Z","shell.execute_reply":"2023-10-05T03:29:10.729218Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_train_dataset(train_data, train_events, drop_nulls=False) :\n \n series_ids = train_data['series_id'].unique(maintain_order=True).to_list()\n X, y = pl.DataFrame(), pl.DataFrame()\n for idx in tqdm(series_ids) : \n \n # Normalizing sample features\n sample = train_data.filter(pl.col('series_id')==idx).with_columns(\n [(pl.col(col) / pl.col(col).std()).cast(pl.Float32) for col in feature_cols if col != 'hour']\n )\n \n events = train_events.filter(pl.col('series_id')==idx)\n \n if drop_nulls : \n # Removing datapoints on dates where no data was recorded\n sample = sample.filter(\n pl.col('timestamp').dt.date().is_in(events['timestamp'].dt.date())\n )\n \n X = X.vstack(sample[id_cols + feature_cols])\n\n onsets = events.filter((pl.col('event') == 'onset') & (pl.col('step') != None))['step'].to_list()\n wakeups = events.filter((pl.col('event') == 'wakeup') & (pl.col('step') != None))['step'].to_list()\n\n # NOTE: This will break if there are event series without any recorded onsets or wakeups\n y = y.vstack(sample.with_columns(\n sum([(onset <= pl.col('step')) & (pl.col('step') <= wakeup) for onset, wakeup in zip(onsets, wakeups)]).cast(pl.Boolean).alias('asleep')\n ).select('asleep')\n )\n \n y = y.to_numpy().ravel()\n \n return X, y","metadata":{"execution":{"iopub.status.busy":"2023-09-29T23:27:53.623791Z","iopub.execute_input":"2023-09-29T23:27:53.624622Z","iopub.status.idle":"2023-09-29T23:27:53.636552Z","shell.execute_reply.started":"2023-09-29T23:27:53.624581Z","shell.execute_reply":"2023-09-29T23:27:53.635486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_events(series, classifier) :\n '''\n Takes a time series and a classifier and returns a formatted submission dataframe.\n '''\n \n series_ids = series['series_id'].unique(maintain_order=True).to_list()\n events = pl.DataFrame(schema={'series_id':str, 'step':int, 'event':str, 'score':float})\n\n for idx in tqdm(series_ids) : \n\n # Collecting sample and normalizing features\n scale_cols = [col for col in feature_cols if (col != 'hour') & (series[col].std() !=0)]\n X = series.filter(pl.col('series_id') == idx).select(id_cols + feature_cols).with_columns(\n [(pl.col(col) / series[col].std()).cast(pl.Float32) for col in scale_cols]\n )\n\n # Applying classifier to get predictions and scores\n preds, probs = classifier.predict(X[feature_cols]), classifier.predict_proba(X[feature_cols])[:, 1]\n\n #NOTE: Considered using rolling max to get sleep periods excluding <30 min interruptions, but ended up decreasing performance\n X = X.with_columns(\n pl.lit(preds).cast(pl.Int8).alias('prediction'), \n pl.lit(probs).alias('probability')\n )\n \n # Getting predicted onset and wakeup time steps\n pred_onsets = X.filter(X['prediction'].diff() > 0)['step'].to_list()\n pred_wakeups = X.filter(X['prediction'].diff() < 0)['step'].to_list()\n \n if len(pred_onsets) > 0 : \n \n # Ensuring all predicted sleep periods begin and end\n if min(pred_wakeups) < min(pred_onsets) : \n pred_wakeups = pred_wakeups[1:]\n\n if max(pred_onsets) > max(pred_wakeups) :\n pred_onsets = pred_onsets[:-1]\n\n # Keeping sleep periods longer than 30 minutes\n sleep_periods = [(onset, wakeup) for onset, wakeup in zip(pred_onsets, pred_wakeups) if wakeup - onset >= 12 * 30]\n\n for onset, wakeup in sleep_periods :\n # Scoring using mean probability over period\n score = X.filter((pl.col('step') >= onset) & (pl.col('step') <= wakeup))['probability'].mean()\n\n # Adding sleep event to dataframe\n events = events.vstack(pl.DataFrame().with_columns(\n pl.Series([idx, idx]).alias('series_id'), \n pl.Series([onset, wakeup]).alias('step'),\n pl.Series(['onset', 'wakeup']).alias('event'),\n pl.Series([score, score]).alias('score')\n ))\n\n # Adding row id column\n events = events.to_pandas().reset_index().rename(columns={'index':'row_id'})\n\n return events","metadata":{"execution":{"iopub.status.busy":"2023-09-29T23:27:54.341892Z","iopub.execute_input":"2023-09-29T23:27:54.342283Z","iopub.status.idle":"2023-09-29T23:27:54.355274Z","shell.execute_reply.started":"2023-09-29T23:27:54.342251Z","shell.execute_reply":"2023-09-29T23:27:54.354367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training Models","metadata":{}},{"cell_type":"code","source":"'''\nfrom sklearn.model_selection import train_test_split\n\ntrain_ids, val_ids = train_test_split(series_ids, train_size=0.7, random_state=42)\n\n# We will collect datapoints at 10 minute intervals for training for validating\ntrain_data = train_series.filter(pl.col('series_id').is_in(train_ids)).take_every(12 * 10).collect()\n\nval_data = train_series.filter(pl.col('series_id').is_in(val_ids)).collect()\nval_solution = train_events.filter(pl.col('series_id').is_in(val_ids)).select(['series_id', 'event', 'step']).to_pandas()\n'''\n# Collecting datapoints at every 5 minutes\ntrain_data = train_series.filter(pl.col('series_id').is_in(series_ids)).take_every(12 * 5).collect()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-29T23:28:14.968045Z","iopub.execute_input":"2023-09-29T23:28:14.968472Z","iopub.status.idle":"2023-09-29T23:33:45.069961Z","shell.execute_reply.started":"2023-09-29T23:28:14.968441Z","shell.execute_reply":"2023-09-29T23:33:45.068781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Creating train dataset\nX_train, y_train = make_train_dataset(train_data, train_events)","metadata":{"execution":{"iopub.status.busy":"2023-09-29T23:33:45.072492Z","iopub.execute_input":"2023-09-29T23:33:45.073202Z","iopub.status.idle":"2023-09-29T23:33:47.553916Z","shell.execute_reply.started":"2023-09-29T23:33:45.073137Z","shell.execute_reply":"2023-09-29T23:33:47.551656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training and validating random forest","metadata":{}},{"cell_type":"code","source":"from sklearn.ensemble import RandomForestClassifier\n\nrf_classifier = RandomForestClassifier(random_state=42)\n\n# Training classifier\nrf_classifier = RandomForestClassifier(n_estimators=500,\n min_samples_leaf=25,\n random_state=42,\n n_jobs=-1)\n\nrf_classifier.fit(X_train[feature_cols], y_train)","metadata":{"execution":{"iopub.status.busy":"2023-09-30T00:07:34.562882Z","iopub.execute_input":"2023-09-30T00:07:34.563276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting feature importances\npx.bar(x=feature_cols, \n y=rf_classifier.feature_importances_,\n title='Random forest feature importances'\n )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Checking performance on validation set\n#rf_submission = get_events(val_data, rf_classifier)\n\n#print(f\"Random forest score: {score(val_solution, rf_submission, tolerances, **column_names)}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Saving classifier \nimport pickle\nwith open('rf_classifier_5m_8h.pkl', 'wb') as f:\n pickle.dump(rf_classifier, f)\n\n#with open('rf_classifier.pkl', 'rb') as f:\n# rf_classifier = pickle.load(f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training and validating gradient boost","metadata":{}},{"cell_type":"code","source":"'''# With SKL\nfrom sklearn.ensemble import GradientBoostingClassifier\n\ngb_classifier = GradientBoostingClassifier(n_estimators=100, learning_rate=0.05, max_depth=10, random_state=42)\ngb_classifier.fit(X_train[feature_cols], y_train)'''","metadata":{"execution":{"iopub.status.busy":"2023-09-20T00:32:29.342993Z","iopub.execute_input":"2023-09-20T00:32:29.343447Z","iopub.status.idle":"2023-09-20T00:51:55.176259Z","shell.execute_reply.started":"2023-09-20T00:32:29.343416Z","shell.execute_reply":"2023-09-20T00:51:55.174488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''# Plotting feature importances\npx.bar(x=feature_cols, \n y=gb_classifier.feature_importances_,\n title='Gradient boosting feature importances'\n )'''","metadata":{"execution":{"iopub.status.busy":"2023-09-20T00:51:55.179395Z","iopub.execute_input":"2023-09-20T00:51:55.180376Z","iopub.status.idle":"2023-09-20T00:51:57.466029Z","shell.execute_reply.started":"2023-09-20T00:51:55.180342Z","shell.execute_reply":"2023-09-20T00:51:57.464688Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n# Checking performance on validation set\ngb_submission = get_events(val_data, gb_classifier)\n\nprint(f\"Gradient boosting score: {score(val_solution, gb_submission, tolerances, **column_names)}\")\n'''","metadata":{"execution":{"iopub.status.busy":"2023-09-20T00:51:57.467479Z","iopub.execute_input":"2023-09-20T00:51:57.468116Z","iopub.status.idle":"2023-09-20T00:56:51.902928Z","shell.execute_reply.started":"2023-09-20T00:51:57.468082Z","shell.execute_reply":"2023-09-20T00:56:51.901915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Applying to test data","metadata":{}},{"cell_type":"code","source":"# Recovering memory\ndel train_data ","metadata":{"execution":{"iopub.status.busy":"2023-09-20T05:22:09.101944Z","iopub.execute_input":"2023-09-20T05:22:09.102333Z","iopub.status.idle":"2023-09-20T05:22:09.107008Z","shell.execute_reply.started":"2023-09-20T05:22:09.102306Z","shell.execute_reply":"2023-09-20T05:22:09.106043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Getting event predictions for test set and saving submission\nsubmission = get_events(test_series.collect(), rf_classifier)\nsubmission.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T05:22:10.01902Z","iopub.execute_input":"2023-09-20T05:22:10.019601Z","iopub.status.idle":"2023-09-20T05:22:10.835397Z","shell.execute_reply.started":"2023-09-20T05:22:10.019563Z","shell.execute_reply":"2023-09-20T05:22:10.83395Z"},"trusted":true},"execution_count":null,"outputs":[]}]}