1_Baseline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**BEFORE WE START**: Make sure that any open notebooks that are no longer needed are shut down to free up compute resources (note that closing the tab is not enough!). This will reduce the risk of running into memory and/or performance issues."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Random forest + temporal models\n",
    "\n",
    "In the previous section, we loaded and processed the data for one participant,\n",
    "used it to train a classifier, and tested the model on another participant.\n",
    "\n",
    "In this section, let's make use of all the data from the 151 participants.\n",
    "We will train the classifier on data from 100 participants, and set aside 51\n",
    "participants for testing.\n",
    "\n",
    "Finally, we will explore ways to account for the temporal dependency such as\n",
    "mode smoothing and hidden Markov models.\n",
    "\n",
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from glob import glob\n",
    "import scipy.stats as stats\n",
    "from imblearn.ensemble import BalancedRandomForestClassifier\n",
    "from sklearn import metrics\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from joblib import Parallel, delayed\n",
    "import urllib\n",
    "import shutil\n",
    "from tqdm.auto import tqdm\n",
    "import utils  # helper functions -- check out utils.py\n",
    "import zipfile\n",
    "\n",
    "# For reproducibility\n",
    "np.random.seed(42)\n",
    "N_JOBS = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Process dataset\n",
    "\n",
    "First of all, we need to process the Capture-24 dataset so that it is ready for training machine learning models.\n",
    "\n",
    "The following code reads each accelerometer data file and slices it into 30-second windows. It creates four arrays: X, Y, T, and pid.\n",
    "\n",
    "**Note**: The following cell downloads the already processed files from https://wearables-files.ndph.ox.ac.uk/files/processed_data.zip instead of running the actual data processing (which can take a long time). We recommend this option for devices with less than 16 GB of RAM. The actual data processing, which is worth understanding, is detailed at the bottom of the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: The required files are already available at /srv.\n",
    "# The following creates a symbolic link (shortcut in Windows) to it in the current directory.\n",
    "if not os.path.exists(\"capture24\"):\n",
    "    os.symlink(\"/srv/capture24\", \"capture24\")\n",
    "if not os.path.exists(\"processed_data\"):\n",
    "    os.symlink(\"/srv/processed_data\", \"processed_data\")\n",
    "\n",
    "# -----------------------------------------------------------------------\n",
    "#\n",
    "# You could also download the files if you uncomment and run the code below\n",
    "# (note that this will take a while).\n",
    "#\n",
    "# -----------------------------------------------------------------------\n",
    "# print(f\"Downloading processed files...\")\n",
    "# os.makedirs(\"processed_data/\", exist_ok=True)\n",
    "# url = \"https://wearables-files.ndph.ox.ac.uk/files/processed_data.zip\"\n",
    "# with urllib.request.urlopen(url) as f_src, open(f\"processed_data.zip\", \"wb\") as f_dst:\n",
    "#     shutil.copyfileobj(f_src, f_dst)\n",
    "# print(f\"Finished downloading\")\n",
    "\n",
    "# print(\"Unzipping...\")\n",
    "# with zipfile.ZipFile(\"processed_data.zip\", \"r\") as f:\n",
    "#     f.extractall(\".\")\n",
    "# print(f\"Finished unzipping\")\n",
    "\n",
    "# -----------------------------------------------------------------------\n",
    "#\n",
    "# If you would like to run the actual data processing, uncomment and run the\n",
    "# code below (note that this will take a while).\n",
    "#\n",
    "# -----------------------------------------------------------------------\n",
    "# def load_all_and_make_windows(datafiles, N=999):\n",
    "\n",
    "#     def worker(datafile):\n",
    "#         X, Y, T = utils.make_windows(utils.load_data(datafile), winsec=30)\n",
    "#         pid = os.path.basename(datafile).split(\".\")[0]  # participant ID\n",
    "#         pid = np.asarray([pid] * len(X))\n",
    "#         return X, Y, T, pid\n",
    "\n",
    "#     results = Parallel(n_jobs=N_JOBS)(\n",
    "#         delayed(worker)(datafile) for datafile in tqdm(datafiles[:N])\n",
    "#     )\n",
    "\n",
    "#     X = np.concatenate([result[0] for result in results])\n",
    "#     Y = np.concatenate([result[1] for result in results])\n",
    "#     T = np.concatenate([result[2] for result in results])\n",
    "#     pid = np.concatenate([result[3] for result in results])\n",
    "\n",
    "#     return X, Y, T, pid\n",
    "\n",
    "# DATAFILES = 'capture24/P[0-9][0-9][0-9].csv.gz'\n",
    "# X, Y, T, pid = load_all_and_make_windows(glob(DATAFILES))\n",
    "# # Save arrays for future use\n",
    "# os.makedirs(\"processed_data/\", exist_ok=True)\n",
    "# np.save(\"processed_data/X.npy\", X)\n",
    "# np.save(\"processed_data/Y.npy\", Y)\n",
    "# np.save(\"processed_data/T.npy\", T)\n",
    "# np.save(\"processed_data/pid.npy\", pid)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load processed files\n",
    "X = np.load('processed_data/X.npy', mmap_mode='r')\n",
    "Y = np.load('processed_data/Y.npy')\n",
    "T = np.load('processed_data/T.npy')\n",
    "pid = np.load('processed_data/pid.npy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1**: Can you tell what each of the arrays (`X`, `Y`, `T` and `pid`) represent? Hint: check their content as well as their shape (`.shape`).\n",
    "\n",
    "**Exercise 2**: How many unique annotations are found in `Y`?\n",
    "\n",
    "**Exercise 3**: What is the shape of `X`? What does each dimension represent?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extract features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_features(xyz):\n",
    "    ''' Extract features. xyz is an array of shape (N,3) '''\n",
    "\n",
    "    feats = {}\n",
    "    feats['xMean'], feats['yMean'], feats['zMean'] = np.mean(xyz, axis=0)\n",
    "    feats['xStd'], feats['yStd'], feats['zStd'] = np.std(xyz, axis=0)\n",
    "    v = np.linalg.norm(xyz, axis=1)  # magnitude stream\n",
    "    feats['mean'], feats['std'] = np.mean(v), np.std(v)\n",
    "\n",
    "    return feats\n",
    "\n",
    "# Extract features\n",
    "X_feats = pd.DataFrame(Parallel(n_jobs=N_JOBS)(delayed(extract_features)(x) for x in tqdm(X)))\n",
    "# The features array will be saved to disk in the current directory for future use\n",
    "X_feats.to_pickle('X_feats.pkl')\n",
    "print(X_feats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3**: Feel free to edit the function above to implement your own features!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Map annotations\n",
    "\n",
    "As in the previous chapter, summarise the annotations by mapping to a reduced set of classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# As before, let's map the text annotations to simplified labels\n",
    "anno_label_dict = pd.read_csv(\n",
    "    \"capture24/annotation-label-dictionary.csv\", \n",
    "    index_col='annotation', \n",
    "    dtype='string'\n",
    ")\n",
    "Y = anno_label_dict.loc[Y, 'label:Willetts2018'].to_numpy()\n",
    "\n",
    "print(\"Label districution (# windows)\")\n",
    "print(pd.Series(Y).value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4**: How many hours of \"bicycling\" in total do we have?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train/test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hold out participants P101-P151 for testing (51 participants)\n",
    "test_ids = [f'P{i}' for i in range(101,152)]\n",
    "mask_test = np.isin(pid, test_ids)\n",
    "mask_train = ~mask_test\n",
    "X_train, Y_train, T_train, pid_train = \\\n",
    "    X_feats[mask_train], Y[mask_train], T[mask_train], pid[mask_train]\n",
    "X_test, Y_test, T_test, pid_test = \\\n",
    "    X_feats[mask_test], Y[mask_test], T[mask_test], pid[mask_test]\n",
    "print(\"Shape of X_train:\", X_train.shape)\n",
    "print(\"Shape of X_test:\", X_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train a random forest classifier\n",
    "\n",
    "*This may take a while*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Argument oob_score=True to be used for HMM smoothing (see later below)\n",
    "clf = BalancedRandomForestClassifier(\n",
    "    n_estimators=100,\n",
    "    replacement=True,\n",
    "    sampling_strategy='not minority',\n",
    "    oob_score=True,\n",
    "    n_jobs=N_JOBS,\n",
    "    random_state=42,\n",
    "    verbose=1\n",
    ")\n",
    "clf.fit(X_train, Y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 5**: Tweak the parameters of the classifier. See the [API reference](https://imbalanced-learn.org/stable/references/generated/imblearn.ensemble.BalancedRandomForestClassifier.html#imblearn.ensemble.BalancedRandomForestClassifier) for parameter descriptions. Here are some ideas:\n",
    "\n",
    "- Increase the number of `n_estimators`\n",
    "- Set `max_features` (default is `sqrt` i.e. square root of the total number of features)\n",
    "- Set `max_depth` (default is `None`, i.e. as deep as needed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_test_pred = clf.predict(X_test)\n",
    "print('\\nClassifier performance')\n",
    "print('Out of sample:\\n', metrics.classification_report(Y_test, Y_test_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overall, the model seems to do well in distinguishing between very inactive\n",
    "periods (\"sit-stand\" and \"sleep\") and very active ones (\"bicycling\"), but there\n",
    "seems to be confusion between the remaining activities.\n",
    "\n",
    "#### Plot predicted vs. true activity profiles\n",
    "\n",
    "Using our utility function, let's plot the activity profile for one participant, e.g.\n",
    "`P101`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = pid_test == 'P101'\n",
    "fig, axs = utils.plot_compare(T_test[mask],\n",
    "                              Y_test[mask],\n",
    "                              Y_test_pred[mask],\n",
    "                              trace=X_test.loc[mask, 'std'])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The profile plots look good at first glance. After all, the majority of\n",
    "activities happen to be of the sedentary type for which the model performs\n",
    "well &mdash; this is reflected by the relatively high `weighted avg` scores in\n",
    "the table report.\n",
    "However, the `macro avg` scores are still low, and we see that the model\n",
    "struggles to classify relevant activities such as bicycling and walking.\n",
    "Moreover, we often find discontinuities during the sleep periods. This is\n",
    "because the model is only trained to classify each window instance independently\n",
    "and does not account for temporal dependencies.\n",
    "\n",
    "## Accounting for temporal dependency\n",
    "\n",
    "### Rolling mode smoothing\n",
    "Let's use rolling mode smoothing to smooth the model predictions: Pick the\n",
    "most popular label within a rolling time window.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mode(alist):\n",
    "    ''' Mode of a list, but return middle element if ambiguous '''\n",
    "    m, c = stats.mode(alist, keepdims=False)\n",
    "    m, c = m.item(), c.item()\n",
    "    if c==1:\n",
    "        return alist[len(alist)//2]\n",
    "    return m\n",
    "\n",
    "def rolling_mode(t, y, window_size='100s'):\n",
    "    y_dtype_orig = y.dtype\n",
    "    # Hack to make it work with pandas.Series.rolling()\n",
    "    y = pd.Series(y, index=t, dtype='category')\n",
    "    y_code_smooth = y.cat.codes.rolling(window_size).apply(mode, raw=True).astype('int')\n",
    "    y_smooth = pd.Categorical.from_codes(y_code_smooth, dtype=y.dtype)\n",
    "    y_smooth = np.asarray(y_smooth, dtype=y_dtype_orig)\n",
    "    return y_smooth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Smooth the predictions of each participant\n",
    "Y_test_pred_smooth = []\n",
    "unqP, indP = np.unique(pid_test, return_index=True)\n",
    "unqP = unqP[np.argsort(indP)]  # keep the order or else we'll scramble our arrays\n",
    "for p in unqP:\n",
    "    mask = pid_test == p\n",
    "    Y_test_pred_smooth.append(rolling_mode(T_test[mask], Y_test_pred[mask]))\n",
    "Y_test_pred_smooth = np.concatenate(Y_test_pred_smooth)\n",
    "\n",
    "print('\\nClassifier performance -- mode smoothing')\n",
    "print('Out of sample:\\n', metrics.classification_report(Y_test, Y_test_pred_smooth))\n",
    "\n",
    "# Check again participant\n",
    "mask = pid_test == 'P101'\n",
    "fig, axs = utils.plot_compare(T_test[mask],\n",
    "                              Y_test[mask],\n",
    "                              Y_test_pred_smooth[mask],\n",
    "                              trace=X_test.loc[mask, 'std'])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simple mode smoothing already improved performance slightly.\n",
    "\n",
    "### Hidden Markov Model\n",
    "\n",
    "A more principled approch is to use a Hidden Markov Model (HMM), which can help\n",
    "smooth the sequence of predictions and avoid weird sequences such as\n",
    "\"sleep\"-\"bicycling\"-\"sleep\". \n",
    "Here the random forest predictions are interpreted as \"observations\" of the\n",
    "\"hidden ground truth\". The emission matrix can be estimated from probabilistic\n",
    "predictions of model, and the transition matrix can be estimated from the ground\n",
    "truth sequence of activities. The prior probabilities can be set as the rates\n",
    "observed in the dataset, or a uniform (uninformative) prior.\n",
    "\n",
    "Check `utils.train_hmm` and `utils.viterbi` for implementationd details.\n",
    "\n",
    "**Exercise 5**: How does random forest provide out-of-bag estimations for free?\n",
    "\n",
    "**Exercise 6**: Why do we need out-of-bag estimations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the conveniently provided out-of-bag probability predictions from the\n",
    "# random forest training process.\n",
    "Y_train_prob = clf.oob_decision_function_  # out-of-bag probability predictions\n",
    "labels = clf.classes_  # need this to know the label order of cols of Y_train_prob\n",
    "hmm_params = utils.train_hmm(Y_train_prob, Y_train, labels)  # obtain HMM matrices/params\n",
    "Y_test_pred_hmm = utils.viterbi(Y_test_pred, hmm_params)  # smoothing\n",
    "print('\\nClassifier performance -- HMM smoothing')\n",
    "print('Out of sample:\\n', metrics.classification_report(Y_test, Y_test_pred_hmm))\n",
    "\n",
    "# Check again participant\n",
    "mask = pid_test == 'P101'\n",
    "fig, ax = utils.plot_compare(T_test[mask],\n",
    "                             Y_test[mask],\n",
    "                             Y_test_pred_hmm[mask],\n",
    "                             trace=X_test.loc[mask, 'std'])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HMM further improves the performance scores.\n",
    "\n",
    "## Is a simple logistic regression enough?\n",
    "\n",
    "*This takes a while*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf_LR = LogisticRegression(\n",
    "    penalty=None,\n",
    "    max_iter=1000,\n",
    "    random_state=42)\n",
    "scaler = StandardScaler()\n",
    "pipe = make_pipeline(scaler, clf_LR)\n",
    "pipe.fit(X_train, Y_train)\n",
    "\n",
    "Y_test_pred_LR = pipe.predict(X_test)\n",
    "\n",
    "# HMM smoothing\n",
    "Y_train_LR_prob = pipe.predict_proba(X_train)  # sorry! LR doesn't provide OOB estimates for free\n",
    "labels = pipe.classes_\n",
    "hmm_params_LR = utils.train_hmm(Y_train_LR_prob, Y_train, labels)\n",
    "Y_test_pred_LR_hmm = utils.viterbi(Y_test_pred_LR, hmm_params_LR)  # smoothing\n",
    "\n",
    "print('\\nClassifier performance -- Logistic regression')\n",
    "print('Out of sample:\\n', metrics.classification_report(Y_test, Y_test_pred_LR_hmm))\n",
    "\n",
    "# Check again participant\n",
    "mask = pid_test == 'P101'\n",
    "fig, axs = utils.plot_compare(T_test[mask],\n",
    "                      Y_test[mask],\n",
    "                      Y_test_pred_LR_hmm[mask],\n",
    "                      trace=X_test.loc[mask, 'std'])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The LR model performed well on the easier classes \"sleep\" and \"sit-stand\",\n",
    "but was much worse on all the other classes."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}