{"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"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Introduction and Objective\n\nIn research of freezing of gait for patients with parkinson's disease, I came across a research paper that explained the phenomena and its challenges quite well. The following information is from [Freezing of gait: moving forward on a mysterious clinical phenomenon (Nutt et. al.)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7293393/)\n<br></br>\nFreezing of gait (FoG) is a disabling clinical phenomenon common in advanced Parkinson’s disease (PD). Although it is easy to recognize, the phenomenon is pretty difficult to define. This is due to the definition including a variety of episodes:\n\n1. inability to initiate gait (“start hesitation”) \n2. arrests in forward progression during walking (“turn” hesitation)\n3. as well as episodes of shuffling forward with steps that are millimetres to a couple of centimetres in length\n\nFoG can last anywhere from a few seconds, which is most common, to over 30 seconds. FoG also seems to be nearly continuous meaning that the patient is inable to generate any steps long enough to provide useful assistance (or not) in moving their bodies forward. Several important features have been found to accompany FoG:\n<br>\n(1) the foot or toe does not leave the ground or only barely clears the support surface; \n<br>\n(2) alternate trembling of the legs occurs at a frequency of 3–8 Hz;\n<br>\n(3) hastening, or an increase in cadence with a decrease in step length, often precedes FoG; \n<br>\n(4) a subjective feeling of the feet being glued to the floor accompanies episodes of freezing; \n<br>\n(5) FoG is commonly precipitated or relieved by various cues; and (6) FoG can be asymmetrical, affecting mainly one foot or being elicited more easily by turning in one direction.\n<br></br>\nDue to this varied definition with a collective of features, FoG appears to be able to represent several different syndromes with different underlying mechanisms. Three different patterns have been suggested for FoG: \n<br>\n(1) trembling in place: alternating tremor of the legs (knees); \n<br>\n(2) shuffling forward: very short, shuffling steps; \n<br>\n(3) complete akinesia: no movement of the limbs or trunk, but this pattern is uncommon.\n<br>\n\nThe objective of this notebook is to model 3D lower-back sensor data to be able to model event types that indicate freezing of gait, ultimately helping researchers better understand when and why FOG episodes occur.","metadata":{}},{"cell_type":"markdown","source":"## Dataset Curation and Modeling Overview\nHere's a quick recap of the 3 datasets that were initially provided in the competition:\n<br>\n1. The tDCS FOG (tdcsfog) dataset, comprising data series collected in the lab, as subjects completed a FOG-provoking protocol.\n2. The DeFOG (defog) dataset, comprising data series collected in the subject's home, as subjects completed a FOG-provoking protocol.\n3. The Daily Living (daily) dataset, comprising one week of continuous 24/7 recordings from sixty-five subjects. Forty-five subjects exhibit FOG symptoms and also have series in the defog dataset, while the other twenty subjects do not exhibit FOG symptoms and do not have series elsewhere in the data.\n<br>\nTrials from the tdcsfog and defog datasets were videotaped and annotated by expert reviewers documented the freezing of gait episodes. That is, the start, end and type of each episode were marked by the experts. Series in the daily dataset are unannotated. You will be detecting FOG episodes for the tdcsfog and defog series. You may wish to apply unsupervised or semi-supervised methods to the series in the daily dataset to support your detection modelling.\n\nSince the evaluation criteria will be only on valid examples, we will only be using valid, annotated records in our supervised learning models. Some key notes on the datasets:\n- the tDCS FOG dataset is all valid examples with annotated labels of event types\n- the tDCS FOG dataset recorded acceleration at 128 Hz in units of m/s^2\n- the DeFOG dataset has a combination of valid and invalid examples where invalid resembles examples that experts were inable to distinguish the appropriate event type label\n- the DeFOG dataset recorded acceleration at 100 Hz in units of g.\n- the Daily Living dataset is unannotated, but recorded acceleration the same as the DeFOG dataset.\n<br>\n\nThe DeFOG dataset will first be filtered to only valid examples; according to the data dictionary this corresponds to where \"Valid\" = 1 and \"Task\" = 1. This isn't relevant for the tDCS FOG dataset. Additionally, each dataset will get joined to its corresponding metadata file to incorporate medication data.\n\nOur baseline will be leveraging only the acceleration features to predict for our target variable (EventType, which will need to be created) and over time feature engineering will be employed to incorporate time-based features as well. The dataset being at the sub-second level makes it difficult to incorporate features from the subject, tasks, events, and metadata files.\n\nDue to the differences in data collection, units, and environment between the two main datasets, I will create a model for each and ensemble the two together via a VotingClassifier which will be leveraged on the test set. This will allow for each model to pick up on the subtle patterns in each environment's accelerometer data while also covering for the potential biases of the other.","metadata":{}},{"cell_type":"markdown","source":"## Data Preprocessing and EDA","metadata":{}},{"cell_type":"markdown","source":"### Setup and Helper Functions","metadata":{}},{"cell_type":"code","source":"!pip install tsflex\n!pip install seglearn","metadata":{"execution":{"iopub.status.busy":"2023-04-08T16:53:19.642894Z","iopub.execute_input":"2023-04-08T16:53:19.643606Z","iopub.status.idle":"2023-04-08T16:53:41.081382Z","shell.execute_reply.started":"2023-04-08T16:53:19.643566Z","shell.execute_reply":"2023-04-08T16:53:41.080137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Imports and constant variables\n\nimport pandas as pd\nimport numpy as np\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy import stats\nfrom scipy.stats import pearsonr\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import train_test_split\nimport xgboost as xgb\nimport hyperopt\nfrom hyperopt import fmin, hp, tpe, rand, Trials, STATUS_OK\nfrom hyperopt.early_stop import no_progress_loss\nfrom hyperopt.pyll.base import scope\nfrom sklearn.calibration import CalibratedClassifierCV\nfrom sklearn.metrics import (\n    accuracy_score,\n    precision_score,\n    recall_score,\n    f1_score\n)\nfrom seglearn.feature_functions import base_features, emg_features\nfrom tsflex.features import FeatureCollection, MultipleFeatureDescriptors\nfrom tsflex.features.integrations import seglearn_feature_dict_wrapper\n\nMAIN_ROOT = '/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction'\nDEFOG_PATH = f'{MAIN_ROOT}/train/defog/'\nDEFOG_METADATA_PATH = f'{MAIN_ROOT}/defog_metadata.csv'\nTDCSFOG_PATH = f'{MAIN_ROOT}/train/tdcsfog/'\nTDCSFOG_METADATA_PATH = f'{MAIN_ROOT}/tdcsfog_metadata.csv'\nDEFOG_TEST_PATH = f'{MAIN_ROOT}/test/defog/'\nTDCSFOG_TEST_PATH = f'{MAIN_ROOT}/test/tdcsfog/'\n\nRANDOM_SEED = 0","metadata":{"execution":{"iopub.status.busy":"2023-04-08T16:53:41.083922Z","iopub.execute_input":"2023-04-08T16:53:41.085482Z","iopub.status.idle":"2023-04-08T16:53:42.125258Z","shell.execute_reply.started":"2023-04-08T16:53:41.085433Z","shell.execute_reply":"2023-04-08T16:53:42.124241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Reduce Memory Usage\n# reference : https://www.kaggle.com/code/arjanso/reducing-dataframe-memory-size-by-65 @ARJANGROEN\n\ndef reduce_memory_usage(df):\n    \n    start_mem = df.memory_usage().sum() / 1024**2\n    print('Memory usage of dataframe is {:.2f} MB'.format(start_mem))\n    \n    for col in df.columns:\n        col_type = df[col].dtype.name\n        if ((col_type != 'datetime64[ns]') & (col_type != 'category')):\n            if (col_type != 'object'):\n                c_min = df[col].min()\n                c_max = df[col].max()\n\n                if str(col_type)[:3] == 'int':\n                    if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                        df[col] = df[col].astype(np.int8)\n                    elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                        df[col] = df[col].astype(np.int16)\n                    elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                        df[col] = df[col].astype(np.int32)\n                    elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                        df[col] = df[col].astype(np.int64)\n\n                else:\n                    if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n                        df[col] = df[col].astype(np.float16)\n                    elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                        df[col] = df[col].astype(np.float32)\n                    else:\n                        pass\n            else:\n                df[col] = df[col].astype('category')\n    mem_usg = df.memory_usage().sum() / 1024**2 \n    print(\"Memory usage became: \",mem_usg,\" MB\")\n    \n    return df","metadata":{"execution":{"iopub.status.busy":"2023-04-08T16:53:42.126899Z","iopub.execute_input":"2023-04-08T16:53:42.127703Z","iopub.status.idle":"2023-04-08T16:53:42.157486Z","shell.execute_reply.started":"2023-04-08T16:53:42.127663Z","shell.execute_reply":"2023-04-08T16:53:42.155834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Helper functions for visualizations\n\ndef plot_value_counts_bar_chart(df, column):\n    \"\"\"\n    Creates bar chart for one column based on the value counts for its unique values. Best for small number of unique integer values.\n    Args:\n        df: dataframe to visualize\n        column: column in cleaned dataframe that has numeric values.\n    Returns:\n        bar chart showing value counts of each integer values in the column. Max value is colored.\n    \"\"\"\n\n    df_viz = df[column].value_counts().reset_index()\n    values = df[column].value_counts().reset_index().sort_values('index')[column].values\n    colors = ['slategrey' if (x < max(values)) else 'deeppink' for x in values]\n\n    fig, ax = plt.subplots(figsize=(25,5))\n\n    ax = sns.barplot(\n        data=df_viz, \n        x=df_viz['index'], \n        y=df_viz[column], \n        palette=colors\n    )\n\n    ax.set(xlabel=column)\n    ax.set(ylabel='Counts')\n    ax.set_title(f'{column} Value Counts')\n    custom_style = {\n        'axes.labelcolor': 'black',\n        'xtick.color': 'black',\n        'ytick.color': 'black',\n        'font.size': 12\n    }\n    ax = sns.set_style('darkgrid', rc=custom_style)\n    plt.show()\n    \ndef plot_numerical_distribution(df, column):\n    \"\"\"\n    Creates histogram for one numerical column. Best for columns with float dtype.\n    Args:\n        df: dataframe to visualize\n        column: column in cleaned dataframe that has numeric values.\n    Returns:\n        histogram chart showing the column's distribution.\n    \"\"\"\n    k2, p = stats.normaltest(df[column])\n    alpha = 1e-3\n    if p < alpha:\n        color = 'slategrey'\n    else:\n        color = 'deeppink'\n\n    fig, ax = plt.subplots(figsize=(25,5))\n\n    ax = sns.histplot(\n        data=df, \n        x=df[column], \n        stat='density', \n        kde=True, \n        color=color\n    )\n\n    ax.set(xlabel=column)\n    ax.set_title(f'{column} Distribution, p = {p}')\n    custom_style = {\n        'axes.labelcolor': 'black',\n        'xtick.color': 'black',\n        'ytick.color': 'black',\n        'font.size': 12\n    }\n    ax = sns.set_style('darkgrid', rc=custom_style)\n    plt.show()\n    \ndef plot_line_chart(df, x, y, hue=None):\n    \"\"\"\n    Creates line chart for one column with optional hue for multiple lines.\n    Args:\n        df: dataframe to visualize\n        x: column in cleaned dataframe that is the time value\n        y: values to plot over time\n    Returns:\n        line chart of selected values over time\n    \"\"\"\n    palette = sns.color_palette('flare', n_colors=3)\n\n    fig, ax = plt.subplots(figsize=(25,5))\n\n    ax = sns.lineplot(\n        data=df, \n        x=df[x], \n        y=df[y], \n        hue=hue,\n        palette=palette\n    )\n\n    ax.set(xlabel=x)\n    ax.set(ylabel=y)\n    ax.set_title(f'{y} Over Time')\n    custom_style = {\n        'axes.labelcolor': 'black',\n        'xtick.color': 'black',\n        'ytick.color': 'black',\n        'font.size': 12\n    }\n    ax = sns.set_style('darkgrid', rc=custom_style)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-08T16:53:42.164899Z","iopub.execute_input":"2023-04-08T16:53:42.165799Z","iopub.status.idle":"2023-04-08T16:53:42.189503Z","shell.execute_reply.started":"2023-04-08T16:53:42.165749Z","shell.execute_reply":"2023-04-08T16:53:42.188415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Exploring Subjects","metadata":{}},{"cell_type":"code","source":"# Read in subjects file\n\nsubjects = pd.read_csv(\n    f'{MAIN_ROOT}/subjects.csv',\n    index_col='Subject')\nsubjects = reduce_memory_usage(subjects)\n\nsubjects","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:53:53.535311Z","iopub.execute_input":"2023-04-07T17:53:53.536905Z","iopub.status.idle":"2023-04-07T17:53:53.612970Z","shell.execute_reply.started":"2023-04-07T17:53:53.536867Z","shell.execute_reply":"2023-04-07T17:53:53.611989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read in defog metadata file\n\nmetadata_defog = pd.read_csv(\n    DEFOG_METADATA_PATH\n)\nmetadata_defog = reduce_memory_usage(metadata_defog)\n\n# Read in tdcsfog metadata file\n\nmetadata_tdcsfog = pd.read_csv(\n    TDCSFOG_METADATA_PATH\n)\nmetadata_tdcsfog = reduce_memory_usage(metadata_tdcsfog)","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:53:53.617318Z","iopub.execute_input":"2023-04-07T17:53:53.619657Z","iopub.status.idle":"2023-04-07T17:53:53.655921Z","shell.execute_reply.started":"2023-04-07T17:53:53.619619Z","shell.execute_reply":"2023-04-07T17:53:53.654821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(subjects['Sex'].value_counts())\nplot_value_counts_bar_chart(subjects, 'Sex')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:53:53.660405Z","iopub.execute_input":"2023-04-07T17:53:53.662892Z","iopub.status.idle":"2023-04-07T17:53:53.952162Z","shell.execute_reply.started":"2023-04-07T17:53:53.662855Z","shell.execute_reply":"2023-04-07T17:53:53.951002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_numerical_distribution(subjects, 'Age')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:53:53.956389Z","iopub.execute_input":"2023-04-07T17:53:53.958760Z","iopub.status.idle":"2023-04-07T17:53:54.418342Z","shell.execute_reply.started":"2023-04-07T17:53:53.958720Z","shell.execute_reply":"2023-04-07T17:53:54.417194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Preprocessing and Exploring DeFOG Dataset","metadata":{}},{"cell_type":"code","source":"# Read in all defog train files. Creating time based features for window of 100 due to DeFOG files read in at 100Hz. Stride is at 50 so model can pick up time continuity in the overlap.\n\ndefog_files = Path(DEFOG_PATH).glob('*.csv')\ntrain_defog = pd.concat(\n    (pd.read_csv(f).assign(FileName=f.stem) for f in defog_files), \n    ignore_index=True\n)\ntrain_defog = reduce_memory_usage(train_defog)\ntrain_defog.set_index(['FileName'], inplace=True)\ntrain_defog","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:53:54.422953Z","iopub.execute_input":"2023-04-07T17:53:54.425440Z","iopub.status.idle":"2023-04-07T17:54:20.920711Z","shell.execute_reply.started":"2023-04-07T17:53:54.425396Z","shell.execute_reply":"2023-04-07T17:54:20.919622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train_defog['Valid'].value_counts())\nplot_value_counts_bar_chart(train_defog, 'Valid')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:20.925856Z","iopub.execute_input":"2023-04-07T17:54:20.927031Z","iopub.status.idle":"2023-04-07T17:54:21.962190Z","shell.execute_reply.started":"2023-04-07T17:54:20.926991Z","shell.execute_reply":"2023-04-07T17:54:21.960865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train_defog['Task'].value_counts())\nplot_value_counts_bar_chart(train_defog, 'Task')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:21.964015Z","iopub.execute_input":"2023-04-07T17:54:21.964775Z","iopub.status.idle":"2023-04-07T17:54:22.633694Z","shell.execute_reply.started":"2023-04-07T17:54:21.964728Z","shell.execute_reply":"2023-04-07T17:54:22.632525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Filtering to only valid records\n\ntrain_defog_valid = train_defog[(train_defog['Valid'] == 1) & (train_defog['Task'] == 1)].drop(columns=['Valid', 'Task'])\ntrain_defog_valid.shape","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:22.635310Z","iopub.execute_input":"2023-04-07T17:54:22.635703Z","iopub.status.idle":"2023-04-07T17:54:23.055679Z","shell.execute_reply.started":"2023-04-07T17:54:22.635663Z","shell.execute_reply":"2023-04-07T17:54:23.054622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This takes forever to run, so uncommenting for now\n\n# defog_acceleration_timesteps = train_defog_valid[['Time', 'AccV', 'AccML', 'AccAP']].reset_index(drop=True).set_index('Time').stack().reset_index().rename(columns={'level_1': 'Acc', 0: 'Value'})\n# plot_line_chart(defog_acceleration_timesteps, 'Time', 'Value', hue='Acc')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:23.057347Z","iopub.execute_input":"2023-04-07T17:54:23.057753Z","iopub.status.idle":"2023-04-07T17:54:23.063309Z","shell.execute_reply.started":"2023-04-07T17:54:23.057714Z","shell.execute_reply":"2023-04-07T17:54:23.062118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Creating EventType target column, as a stacking of StartHesitation, Turn, and Walking columns\nevent_type_conditions = [\n    (train_defog_valid['StartHesitation'] == 1),\n    (train_defog_valid['Turn'] == 1),\n    (train_defog_valid['Walking'] == 1)\n]\nevent_types = ['StartHesitation', 'Turn', 'Walking']\ntrain_defog_valid['EventType'] = np.select(\n    event_type_conditions, \n    event_types, \n    default='Normal'\n)\nle = LabelEncoder()\ntrain_defog_valid['Target'] = le.fit_transform(train_defog_valid['EventType'])\n\ntrain_defog_valid.drop(columns=event_types, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:23.074286Z","iopub.execute_input":"2023-04-07T17:54:23.075083Z","iopub.status.idle":"2023-04-07T17:54:24.979068Z","shell.execute_reply.started":"2023-04-07T17:54:23.075046Z","shell.execute_reply":"2023-04-07T17:54:24.977997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train_defog_valid['EventType'].value_counts())\nplot_value_counts_bar_chart(train_defog_valid, 'EventType')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:24.980764Z","iopub.execute_input":"2023-04-07T17:54:24.981155Z","iopub.status.idle":"2023-04-07T17:54:25.836409Z","shell.execute_reply.started":"2023-04-07T17:54:24.981113Z","shell.execute_reply":"2023-04-07T17:54:25.835202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Preprocessing and Exploring tDCSFOG Dataset","metadata":{}},{"cell_type":"code","source":"# Read in all tdcsfog train files. Creating time based features for window of 128 due to tDCSFOG files read in at 100Hz. Stride is at 64 so model can pick up time continuity in the overlap.\n\ntdcsfog_files = Path(TDCSFOG_PATH).glob('*.csv')\ntrain_tdcsfog = pd.concat(\n    (pd.read_csv(f).assign(FileName=f.stem) for f in tdcsfog_files), \n    ignore_index=True\n)\ntrain_tdcsfog = reduce_memory_usage(train_tdcsfog)\ntrain_tdcsfog.set_index(['FileName'], inplace=True)\ntrain_tdcsfog","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:25.838247Z","iopub.execute_input":"2023-04-07T17:54:25.838684Z","iopub.status.idle":"2023-04-07T17:54:40.726333Z","shell.execute_reply.started":"2023-04-07T17:54:25.838641Z","shell.execute_reply":"2023-04-07T17:54:40.725197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This takes forever to run, so uncommenting for now\n\n# tdcsfog_acceleration_timesteps = train_tdcsfog[['Time', 'AccV', 'AccML', 'AccAP']].reset_index(drop=True).set_index('Time').stack().reset_index().rename(columns={'level_1': 'Acc', 0: 'Value'})\n# plot_line_chart(tdcsfog_acceleration_timesteps, 'Time', 'Value', hue='Acc')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:40.727971Z","iopub.execute_input":"2023-04-07T17:54:40.728360Z","iopub.status.idle":"2023-04-07T17:54:40.733013Z","shell.execute_reply.started":"2023-04-07T17:54:40.728312Z","shell.execute_reply":"2023-04-07T17:54:40.731843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Creating EventType target column, as a stacking of StartHesitation, Turn, and Walking columns\nevent_type_conditions = [\n    (train_tdcsfog['StartHesitation'] == 1),\n    (train_tdcsfog['Turn'] == 1),\n    (train_tdcsfog['Walking'] == 1)\n]\nevent_types = ['StartHesitation', 'Turn', 'Walking']\ntrain_tdcsfog['EventType'] = np.select(\n    event_type_conditions, \n    event_types, \n    default='Normal'\n)\nle = LabelEncoder()\ntrain_tdcsfog['Target'] = le.fit_transform(train_tdcsfog['EventType'])\n\ntrain_tdcsfog.drop(columns=event_types, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:40.734736Z","iopub.execute_input":"2023-04-07T17:54:40.735082Z","iopub.status.idle":"2023-04-07T17:54:44.153515Z","shell.execute_reply.started":"2023-04-07T17:54:40.735048Z","shell.execute_reply":"2023-04-07T17:54:44.152164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train_defog_valid['EventType'].value_counts())\nplot_value_counts_bar_chart(train_defog_valid, 'EventType')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:44.155182Z","iopub.execute_input":"2023-04-07T17:54:44.155592Z","iopub.status.idle":"2023-04-07T17:54:45.028298Z","shell.execute_reply.started":"2023-04-07T17:54:44.155533Z","shell.execute_reply":"2023-04-07T17:54:45.026570Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Modeling XGBoost with Hyperopt and Evaluation","metadata":{}},{"cell_type":"markdown","source":"### Baseline Models\n\nCreating a class makes it really easy to make the model for each dataset modular and consistent. Key things to note:\n1. Using 'hist' for tree_method significantly speeds up training, but 'gpu_hist' is far better\n2. 'multi:softprob' objective allows us to output probabilities\n3. 'map' eval_metric ensures model is evaluated with mean average precision on validation data\n4. In HyperOpt, our loss function is the weighted precision on our validation set due to high class imbalance","metadata":{}},{"cell_type":"code","source":"class XGBHyperopt:\n    '''\n    This class fits a XGBoost model with hyperparameters tuned with HyperOpt, given a dataset.\n    '''\n    def __init__(self, dataset: pd.DataFrame) -> None:\n        self.dataset = dataset\n        self.target = 'Target'\n        self.search_space = {\n            'tree_method': 'gpu_hist',\n            'objective': 'multi:softprob',\n            'eval_metric': 'map',\n            'max_depth': scope.int(hp.quniform('max_depth', 4, 10, 1)),\n            'subsample': hp.uniform('subsample', 0.6, 0.9),\n            'learning_rate': hp.uniform('learning_rate', .001, 0.5),\n            'min_child_weight': hp.quniform('min_child_weight', 1, 20, 1),\n            'reg_alpha': hp.quniform('reg_alpha', 0, 1, .001),\n            'reg_lambda': hp.quniform('reg_lambda', 0, 100, 1),\n            'gamma': hp.uniform ('gamma', 1, 9),\n            'verbosity': 1,\n            'random_state': RANDOM_SEED\n        }        \n        \n    def preprocess(self, test_size=0.25):\n        '''\n        Splits the provided dataset into train and validation sets.\n        '''\n        X = self.dataset.drop(columns=[self.target])\n        y = self.dataset[self.target].values\n\n        # Splitting training/validation sets\n        print('Splitting dataset into training and validation sets...')\n        self.X_train, self.X_val, self.y_train, self.y_val = train_test_split(\n            X, \n            y, \n            stratify=y,\n            test_size=test_size,\n            random_state=RANDOM_SEED\n        )\n\n    def evaluate_multiclassif_model(self, y_test, y_pred):\n\n        \"\"\"\n        Creating evaluation metrics for classification predictions.\n\n        Args:\n            y_test: representing the true values of the target.\n            y_pred: representing the predicted values of the target.\n            y_pred_proba: representing the predicted probability values of the target.\n        Returns:\n            Dictionary of classification evaluation metrics\n        \"\"\"\n\n        accuracy = accuracy_score(y_test, y_pred).round(3)\n        # Getting weighted and macro averages of each below metric since task is multiclass classification\n        weighted_precision = precision_score(y_test, y_pred, average='weighted', labels=np.unique(y_pred), zero_division=1).round(3)\n        weighted_recall = recall_score(y_test, y_pred, average='weighted', labels=np.unique(y_pred), zero_division=1).round(3)\n        weighted_f1 = f1_score(y_test, y_pred, average='weighted', labels=np.unique(y_pred), zero_division=1).round(3)\n        \n        macro_precision = precision_score(y_test, y_pred, average='macro', labels=np.unique(y_pred), zero_division=1).round(3)\n        macro_recall = recall_score(y_test, y_pred, average='macro', labels=np.unique(y_pred), zero_division=1).round(3)\n        macro_f1 = f1_score(y_test, y_pred, average='macro', labels=np.unique(y_pred), zero_division=1).round(3)\n\n    #     print(\"Weighted Precision Score : \", weighted_precision)\n    #     print(\"Macro Precision Score : \", macro_precision)\n\n        return {\n            'Accuracy': accuracy, \n            'Weighted Precision': weighted_precision, \n            'Weighted Recall': weighted_recall, \n            'Weighted F1': weighted_f1, \n            'Macro Precision': macro_precision, \n            'Macro Recall': macro_recall, \n            'Macro F1': macro_f1, \n        }\n\n    def calibrate_best_model(self):\n        self.best_calibrated_model = CalibratedClassifierCV(\n            self.best_model, \n            method='isotonic', \n            cv='prefit'\n        )\n        self.best_calibrated_model = self.best_calibrated_model.fit(self.X_val, self.y_val)\n        \n    def hpopt_fit(self, calibrate=False):\n        '''\n        Sweeps the search space with a XGBoost model to minimize the loss function. Can be adapted to provide alternative loss functions to measure the validation set with.\n        '''\n        print('Modeling XGBoost with HyperOpt optimization...')\n        def train_model(params):\n            \n            model = xgb.XGBClassifier(**params)\n            model = model.fit(self.X_train, self.y_train)\n            y_train_pred = model.predict(self.X_train)\n            y_train_pred_proba = model.predict_proba(self.X_train)[:, 1]\n            y_pred = model.predict(self.X_val)\n            y_pred_proba = model.predict_proba(self.X_val)[:, 1]\n            \n            # Evaluating classification model evaluation metrics for testing set predictions and validation set predictions\n            training_metrics = self.evaluate_multiclassif_model(self.y_train, y_train_pred)\n            validation_metrics = self.evaluate_multiclassif_model(self.y_val, y_pred)\n\n            # Set the loss to validation Weighted Precision so fmin maximizes the Weighted Precision\n            return {\n                'status': STATUS_OK, \n                'loss': -1*validation_metrics['Weighted Precision'],\n                'model': model,\n                'training_metrics': training_metrics,\n                'validation_metrics': validation_metrics,\n                'y_pred': y_pred,\n                'y_pred_proba': y_pred_proba\n            }\n    \n        self.trials = Trials()\n        self.best_params = fmin(\n            fn=train_model, \n            space=self.search_space, \n            algo=tpe.suggest, \n            max_evals=200,\n            early_stop_fn=no_progress_loss(25),\n            trials=self.trials\n        )\n        \n        self.best_model = self.trials.results[np.argmin([r['loss'] for r in self.trials.results])]['model']\n        if calibrate:\n            self.best_model = self.calibrate_best_model()\n            \n    \n    def validation_predictions(self):\n        '''\n        Returns the predictions and probabilities for the validation set\n        '''\n        self.y_pred = self.trials.results[np.argmin([r['loss'] for r in self.trials.results])]['y_pred']\n        self.y_pred_proba = self.trials.results[np.argmin([r['loss'] for r in self.trials.results])]['y_pred_proba']\n        predictions_df = pd.DataFrame([self.y_pred, self.y_pred_proba], columns=['predictions', 'probabilities'])\n        return predictions_df\n    \n    def training_metrics(self):\n        '''\n        Returns the training metrics for the best model\n        '''\n        return self.trials.results[np.argmin([r['loss'] for r in self.trials.results])]['training_metrics']\n    \n    def validation_metrics(self):\n        '''\n        Returns the validation metrics for the best model\n        '''\n        return self.trials.results[np.argmin([r['loss'] for r in self.trials.results])]['validation_metrics']\n\n    def predict(self, X_test):\n        '''\n        Returns the predictions for the test set\n        '''\n        return self.best_model.predict(X_test)\n\n    def feature_importances(self):\n        '''\n        Returns the feature importances for the best model\n        '''\n        return self.best_model.feature_importances_","metadata":{"execution":{"iopub.status.busy":"2023-04-08T16:53:42.190875Z","iopub.execute_input":"2023-04-08T16:53:42.191499Z","iopub.status.idle":"2023-04-08T16:53:42.236076Z","shell.execute_reply.started":"2023-04-08T16:53:42.191457Z","shell.execute_reply":"2023-04-08T16:53:42.234975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"defog_baseline_dataset = train_defog_valid.reset_index(drop=True).drop(columns=['EventType']).set_index('Time')\ntdcsfog_baseline_dataset = train_tdcsfog.reset_index(drop=True).drop(columns=['EventType']).set_index('Time')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T17:54:45.059819Z","iopub.execute_input":"2023-04-07T17:54:45.060273Z","iopub.status.idle":"2023-04-07T17:54:45.509213Z","shell.execute_reply.started":"2023-04-07T17:54:45.060235Z","shell.execute_reply":"2023-04-07T17:54:45.508106Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"defog_model = XGBHyperopt(defog_baseline_dataset)\ndefog_model.preprocess()\ndefog_model.hpopt_fit()","metadata":{"execution":{"iopub.status.busy":"2023-04-06T19:17:09.976699Z","iopub.execute_input":"2023-04-06T19:17:09.977883Z","iopub.status.idle":"2023-04-06T19:36:31.629319Z","shell.execute_reply.started":"2023-04-06T19:17:09.977834Z","shell.execute_reply":"2023-04-06T19:36:31.627919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"defog_model.validation_metrics()","metadata":{"execution":{"iopub.status.busy":"2023-04-06T19:36:31.631612Z","iopub.execute_input":"2023-04-06T19:36:31.632887Z","iopub.status.idle":"2023-04-06T19:36:31.643325Z","shell.execute_reply.started":"2023-04-06T19:36:31.632835Z","shell.execute_reply":"2023-04-06T19:36:31.641615Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tdcsfog_model = XGBHyperopt(tdcsfog_baseline_dataset)\ntdcsfog_model.preprocess()\ntdcsfog_model.hpopt_fit()","metadata":{"execution":{"iopub.status.busy":"2023-04-06T19:38:21.184224Z","iopub.execute_input":"2023-04-06T19:38:21.185428Z","iopub.status.idle":"2023-04-06T20:45:32.478137Z","shell.execute_reply.started":"2023-04-06T19:38:21.185384Z","shell.execute_reply":"2023-04-06T20:45:32.476747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tdcsfog_model.validation_metrics()","metadata":{"execution":{"iopub.status.busy":"2023-04-06T20:45:32.480782Z","iopub.execute_input":"2023-04-06T20:45:32.481337Z","iopub.status.idle":"2023-04-06T20:45:32.490761Z","shell.execute_reply.started":"2023-04-06T20:45:32.481291Z","shell.execute_reply":"2023-04-06T20:45:32.489466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Models with Feature Engineering and Calibration","metadata":{}},{"cell_type":"code","source":"def read_data_and_engineer_features(filename, w, s):\n    '''\n    Function to read in training file and create time series features using tsflex.\n    Args:\n        filename: filename of training data to read in\n        w: window size\n        s: stride size\n    Returns:\n        df: dataframe with time series features joined on\n    '''\n    \n    df = pd.read_csv(filename, index_col=['Time']).assign(FileName=filename.stem)\n    \n    # reference: https://www.kaggle.com/code/jeroenvdd/time-series-tsflex\n    basic_feats = MultipleFeatureDescriptors(\n        functions=seglearn_feature_dict_wrapper(base_features()),\n        series_names=['AccV', 'AccML', 'AccAP'],\n        windows=[w],\n        strides=[s],\n    )\n\n    emg_feats = emg_features()\n    del emg_feats['simple square integral']\n\n    emg_feats = MultipleFeatureDescriptors(\n        functions=seglearn_feature_dict_wrapper(emg_feats),\n        series_names=['AccV', 'AccML', 'AccAP'],\n        windows=[w],\n        strides=[s],\n    )\n\n    fc = FeatureCollection([basic_feats, emg_feats])\n    df_ts = fc.calculate(\n        df, \n        return_df=True, \n        include_final_window=True, \n        approve_sparsity=True, \n        window_idx='begin'\n    ).astype(np.float32)\n    \n    df = df.join(df_ts, how='left')\n    df.fillna(method='ffill', inplace=True)\n    \n    df = reduce_memory_usage(df)\n    \n    return df","metadata":{"execution":{"iopub.status.busy":"2023-04-08T16:53:42.238327Z","iopub.execute_input":"2023-04-08T16:53:42.238649Z","iopub.status.idle":"2023-04-08T16:53:42.251888Z","shell.execute_reply.started":"2023-04-08T16:53:42.238610Z","shell.execute_reply":"2023-04-08T16:53:42.250468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Reading in defog files and creating features for each file - 100 window size due to 100Hz and stride of 50 to ensure time continuity is captured via overlap\ndefog_files = Path(DEFOG_PATH).glob('*.csv')\ntrain_defog_ts = pd.concat(\n    (read_data_and_engineer_features(f, 100, 50) for f in defog_files), \n    ignore_index=True\n)\ntrain_defog_ts.set_index(['FileName'], inplace=True)\n\ntrain_defog_ts = train_defog_ts[(train_defog_ts['Valid'] == 1) & (train_defog_ts['Task'] == 1)].drop(columns=['Valid', 'Task'])\n\n# Creating EventType target column, as a stacking of StartHesitation, Turn, and Walking columns\nevent_type_conditions = [\n    (train_defog_ts['StartHesitation'] == 1),\n    (train_defog_ts['Turn'] == 1),\n    (train_defog_ts['Walking'] == 1)\n]\nevent_types = ['StartHesitation', 'Turn', 'Walking']\ntrain_defog_ts['EventType'] = np.select(\n    event_type_conditions, \n    event_types, \n    default='Normal'\n)\nle = LabelEncoder()\ntrain_defog_ts['Target'] = le.fit_transform(train_defog_ts['EventType'])\n\ntrain_defog_ts.drop(columns=event_types+['EventType'], inplace=True)","metadata":{"execution":{"iopub.status.busy":"2023-04-08T16:53:42.253298Z","iopub.execute_input":"2023-04-08T16:53:42.253704Z","iopub.status.idle":"2023-04-08T17:02:36.309528Z","shell.execute_reply.started":"2023-04-08T16:53:42.253665Z","shell.execute_reply":"2023-04-08T17:02:36.308203Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_defog_ts","metadata":{"execution":{"iopub.status.busy":"2023-04-08T17:02:36.311617Z","iopub.execute_input":"2023-04-08T17:02:36.311973Z","iopub.status.idle":"2023-04-08T17:02:37.605909Z","shell.execute_reply.started":"2023-04-08T17:02:36.311937Z","shell.execute_reply":"2023-04-08T17:02:37.604898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Reading in tdcsfog files and creating features for each file - 128 window size due to 128Hz and stride of 64 to ensure time continuity is captured via overlap\ntdcsfog_files = Path(TDCSFOG_PATH).glob('*.csv')\ntrain_tdcsfog_ts = pd.concat(\n    (read_data_and_engineer_features(f, 128, 64) for f in tdcsfog_files), \n    ignore_index=True\n)\ntrain_tdcsfog_ts.set_index(['FileName'], inplace=True)\n\n# Creating EventType target column, as a stacking of StartHesitation, Turn, and Walking columns\nevent_type_conditions = [\n    (train_tdcsfog_ts['StartHesitation'] == 1),\n    (train_tdcsfog_ts['Turn'] == 1),\n    (train_tdcsfog_ts['Walking'] == 1)\n]\nevent_types = ['StartHesitation', 'Turn', 'Walking']\ntrain_tdcsfog_ts['EventType'] = np.select(\n    event_type_conditions, \n    event_types, \n    default='Normal'\n)\nle = LabelEncoder()\ntrain_tdcsfog_ts['Target'] = le.fit_transform(train_tdcsfog_ts['EventType'])\n\ntrain_tdcsfog_ts.drop(columns=event_types+['EventType'], inplace=True)","metadata":{"execution":{"iopub.status.busy":"2023-04-08T17:02:37.609243Z","iopub.execute_input":"2023-04-08T17:02:37.610082Z","iopub.status.idle":"2023-04-08T17:11:27.491170Z","shell.execute_reply.started":"2023-04-08T17:02:37.610049Z","shell.execute_reply":"2023-04-08T17:11:27.489848Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_tdcsfog_ts","metadata":{"execution":{"iopub.status.busy":"2023-04-08T17:11:27.493943Z","iopub.execute_input":"2023-04-08T17:11:27.494720Z","iopub.status.idle":"2023-04-08T17:11:29.716552Z","shell.execute_reply.started":"2023-04-08T17:11:27.494667Z","shell.execute_reply":"2023-04-08T17:11:29.715598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"defog_model_ts = XGBHyperopt(train_defog_ts)\ndefog_model_ts.preprocess()\ndefog_model_ts.hpopt_fit()","metadata":{"execution":{"iopub.status.busy":"2023-04-07T18:14:45.800156Z","iopub.execute_input":"2023-04-07T18:14:45.800887Z","iopub.status.idle":"2023-04-07T19:47:23.940842Z","shell.execute_reply.started":"2023-04-07T18:14:45.800846Z","shell.execute_reply":"2023-04-07T19:47:23.939543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"defog_model_ts.validation_metrics()","metadata":{"execution":{"iopub.status.busy":"2023-04-07T19:47:23.943752Z","iopub.execute_input":"2023-04-07T19:47:23.946195Z","iopub.status.idle":"2023-04-07T19:47:23.955005Z","shell.execute_reply.started":"2023-04-07T19:47:23.946163Z","shell.execute_reply":"2023-04-07T19:47:23.953813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tdcsfog_model_ts = XGBHyperopt(train_tdcsfog_ts)\ntdcsfog_model_ts.preprocess()\ntdcsfog_model_ts.hpopt_fit()","metadata":{"execution":{"iopub.status.busy":"2023-04-08T17:11:29.718069Z","iopub.execute_input":"2023-04-08T17:11:29.718443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tdcsfog_model_ts.validation_metrics()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Calibrate Models","metadata":{}},{"cell_type":"code","source":"defog_model_ts_calibrated = XGBHyperopt(train_defog_ts)\ndefog_model_ts_calibrated.preprocess()\ndefog_model_ts_calibrated.hpopt_fit(calibrate=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"defog_model_ts_calibrated.validation_metrics()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tdcsfog_model_ts_calibrated = XGBHyperopt(train_tdcsfog_ts)\ntdcsfog_model_ts_calibrated.preprocess()\ntdcsfog_model_ts_calibrated.hpopt_fit(calibrate=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tdcsfog_model_ts_calibrated.validation_metrics()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Combining Models with a VotingClassifier","metadata":{}},{"cell_type":"code","source":"estimators = [\n    ('defog-xgb', defog_model.best_model),\n    ('tdcsfog-xgb', tdcsfog_model.best_model)\n]\nvoting_xgb = VotingClassifier(estimators=estimators, voting='soft')\nvoting_xgb.estimators_ = [defog_model.best_model, tdcsfog_model.best_model]\nvoting_xgb.le_ = le\nvoting_xgb.classes_ = voting_xgb.le_.classes_","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Predict on Test Set and Submit","metadata":{}},{"cell_type":"code","source":"# Reading in defog test file and creating features for each file - 100 window size due to 100Hz and stride of 50 to ensure time continuity is captured via overlap\ndefog_test_file = Path(DEFOG_TEST_PATH).glob('*.csv')\ntest_defog_ts = pd.concat(\n    (read_data_and_engineer_features(f, 100, 50) for f in defog_test_file), \n    ignore_index=True\n)\n\n# Reading in tdcsfog test file and creating features for each file - 128 window size due to 128Hz and stride of 64 to ensure time continuity is captured via overlap\ntdcsfog_test_file = Path(TDCSFOG_TEST_PATH).glob('*.csv')\ntest_tdcsfog_ts = pd.concat(\n    (read_data_and_engineer_features(f, 128, 64) for f in tdcsfog_test_file), \n    ignore_index=True\n)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_defog_predictions = best_defog_model.predict(test_defog_ts)\ntest_tdcsfog_predictions = best_tdcsfog_model.predict(test_tdcsfog_ts)","metadata":{},"execution_count":null,"outputs":[]}]}