lab-1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import json\n",
    "from xgboost.sklearn import XGBClassifier\n",
    "import xgboost as xgb\n",
    "import sklearn\n",
    "import pandas as pd\n",
    "from typing import Tuple\n",
    "from typing import List\n",
    "\n",
    "dataset = pd.read_csv('data/synchronized-eeg.csv',\n",
    "                      parse_dates=['indra_time'],\n",
    "                      index_col='indra_time')\n",
    "# convert to arrays from strings\n",
    "dataset.raw_values = dataset.raw_values.map(json.loads)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "# 1. EEG\n",
    "\n",
    "Neurons produce electrical activity in the brain.\n",
    "\n",
    "If we zoom out far enough (and, for our purposes, we must), these neuronal firings appear most clearly as patterns of constructive or destructive interference. They are typically low-frequency waves, 0.5Hz-70Hz.\n",
    "\n",
    "But, don't take my word for it. Let's open up a signal from our corpus.\n",
    "\n",
    "Here's 512 readings, or one second, of the 30,000 seconds in our dataset. Let's take the 100th reading in our dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "pd.Series(dataset.raw_values[100].values).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "This is *time-domain* data. It represents voltages over time.\n",
    "\n",
    "We can turn this into *frequency-domain* data using an [FFT](https://www.youtube.com/watch?v=iWZNTM139xQ). (You don't need to understand how this works right now). We'll produce what's called a *power spectrum*.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_power_spectrum (\n",
    "    raw_readings: np.array,\n",
    "    sampling_rate: int = 512,\n",
    ") -> pd.Series:\n",
    "    '''\n",
    "    Take raw voltages,\n",
    "    transform into frequency domain,\n",
    "    return a power spectrum.\n",
    "    '''\n",
    "    # FFT the raw readings\n",
    "    fftd = np.fft.fft(raw_readings)\n",
    "    # take absolute value\n",
    "    # producing a symmetrical power spectrum\n",
    "    ps = np.abs(fftd)**2\n",
    "    # since the power spectrum is symmetrical, \n",
    "    # take half\n",
    "    half_len = len(ps)//2\n",
    "    ps = ps[:half_len]\n",
    "#     # we'll calculate the frequencies\n",
    "    window_size = len(raw_readings)\n",
    "    freqs = numpy.fft.fftfreq(window_size, d=1/sampling_rate)\n",
    "#     # splitting that in half to match\n",
    "    freqs = freqs[:half_len]\n",
    "    power_spectrum = pd.DataFrame({\n",
    "        'frequency (Hz)': freqs,\n",
    "        'magnitude': ps,\n",
    "    })\n",
    "    return power_spectrum\n",
    "\n",
    "ps = to_power_spectrum(dataset.raw_values[1000])\n",
    "ps.plot(x='frequency (Hz)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can see some low-frequency data, toward the left side of the graph\n",
    "\n",
    "Let's do a \"band-pass\" between 0.5Hz and 20Hz, effectively \"zooming in\" or \"cropping\" to these frequencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bandpass (\n",
    "    power_spec: pd.Series,\n",
    "    low_cutoff: float,\n",
    "    high_cutoff: float,\n",
    ") -> pd.Series:\n",
    "    return power_spec[(power_spec['frequency (Hz)'] > low_cutoff) & \n",
    "                      (power_spec['frequency (Hz)'] < high_cutoff)]\n",
    "\n",
    "bandpass(ps, 0.5, 20).plot(x='frequency (Hz)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, here's the problem: There's a bunch of other crap in the signal, too!\n",
    "\n",
    "For example, along with the EEG signal we just saw, there's also a bunch of EMG (electromyographic) data, which is produced by the muscles in your face, or from moving your eyes!\n",
    "\n",
    "Most of this signal will be in the *50Hz to 150Hz* range. Eyeblinks will cause big spikes in this frequency band.\n",
    "\n",
    "### TODO\n",
    "Bandpass to the EMG frequency bands, and plot them. If you want, see if you can find any readings with big spikes in this band!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We might also see ECG (electrocardiographic) signals, from the heart.\n",
    "\n",
    "We'll even see background radiation from the electrical system! Here in North America, we'll see a 60Hz spike no matter what we do.\n",
    "\n",
    "# 2. Classification\n",
    "\n",
    "We can go crazy trying to filter all of these artifacts (and many have). \n",
    "\n",
    "However, we are at Berkeley, so we will take a lazy (computational) view.\n",
    "\n",
    "We will shift our perspective in two ways, and in the end, we will not feel obligated to filter artifacts:\n",
    "\n",
    "1. *The electromagnet is a general-purpose **biosensor***. We are sensing a multitude of human activity, including, but not limited to, EEG.\n",
    "\n",
    "2. *Our machine learning will handle the artifacts*. If some artifacts really are artifacts (i.e., they are not informative with respect to our problem domain), then our ML will simply learn to ignore them.\n",
    "\n",
    "`<OPTIONAL NOTE>`\n",
    "\n",
    "These two \"perspective shifts\" are tightly intertwined. If you are up for it, try using the second as a lens for viewing the first. Our assumptions about EEG melt away, and we are left only with a classification problem and a signal (and the assumptions [built into those infrastructures](https://mitpress.mit.edu/books/raw-data-oxymoron) ;) ). \n",
    "\n",
    "We do not need to make specific commitments about what mechanisms produce the effects we are interested in. Thus, we avoid all hard epistemolgical questions (and complaints)! We are also free to discover signals that might come from unexpected places.  \n",
    "\n",
    "If this is over your head, just remember this: If we're trying to classify interest (say), perhaps an unconscious eyebrow raise is our key signal. It would be a shame if we filtered it out, in our myopic quest for a \"clean\" EEG signal!\n",
    "\n",
    "`</OPTIONAL NOTE>`\n",
    "\n",
    "Fine, so what kind of classifier will just handle the artifacts for us?\n",
    "I know, I know, you're probably thinking of neural networks. [That's an option](https://arxiv.org/pdf/1611.08024.pdf).\n",
    "\n",
    "But, for now, we will use something much simpler, and much easier to run on your laptop: [XGBoost](https://github.com/dmlc/xgboost). You can [read up](https://xgboost.readthedocs.io/en/latest/model.html) on XGBoost on your own time, but like the FFT, you don't need to know how it works right now.\n",
    "\n",
    "In fact, you don't even need to know what machine learning is. I will just show you an example of how to use XGBoost, and we can go from there.\n",
    "\n",
    "First, let's make a plain, *untrained* classifier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fresh_clf () -> XGBClassifier:\n",
    "    return XGBClassifier(\n",
    "        # Don't worry about those parameters for now,\n",
    "        # though feel free to look them up if you're interested.\n",
    "        objective= 'binary:logistic',\n",
    "        seed=27)\n",
    "\n",
    "clf = fresh_clf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Right now, this classifier is fresh to the world. It has seen no evil, heard no evil...\n",
    "\n",
    "The \"learning\" part happens when our  classifier is *trained*, or *fit*, using examples. These examples have *labels*, and *features*. Our classifier will study the relationship between features and labels, so that in the future, it can map features to labels (more on this in a second).\n",
    "\n",
    "Our **features** here are pretty clear: electromagnetic signals, our power spectra.\n",
    "\n",
    "Our **labels** can really be anything, but for now, let's pick stuff that's easy from the existing dataset: People who are doing math, and people who are doing nothing. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All the readings during the \"relax\" task\n",
    "relax = dataset[(dataset.label == 'relax') &\n",
    "                (dataset.id == 1) ]\n",
    "# All the readings taken during \"math\" tasks\n",
    "# (Math tasks were labeled math1...math12 -\n",
    "#  `math\\d` is a \"regular expression,\" or RegEx)\n",
    "math = dataset[(dataset.label.str.match('math\\d')) &\n",
    "               (dataset.id == 1) ]\n",
    "\n",
    "len(relax), len(math)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll go about preparing our data for the classifier. \n",
    "\n",
    "Since our features are the power spectra, we should assemble a big collection of power spectra - and nothing else."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def to_features (\n",
    "    df: pd.DataFrame\n",
    ") -> np.array:\n",
    "    power_specs = df.raw_values.apply(to_power_spectrum)\n",
    "    return np.array([row.magnitude.values for row in power_specs])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can make a big list of all features, with their labels removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "relax_features = to_features(relax)\n",
    "math_features = to_features(math)\n",
    "features = np.concatenate([relax_features, math_features])\n",
    "    \n",
    "assert np.all( [ len(feat) == 256 for feat in features ] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, for each of these features, we will need corresponding labels.\n",
    "\n",
    "Since our machine has no semantic knowledge of what these labels \"mean,\" we will simply refer to them as 0 and 1, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = np.array([ 0 for feature in relax_features  ] \\\n",
    "                + [ 1 for feature in math_features ])\n",
    "\n",
    "# list of labels should be the same \n",
    "# as the number of features\n",
    "assert len(labels) == len(features)\n",
    "# first label in the list should be 0\n",
    "assert labels[0] == 0\n",
    "# last label in the list should be 1\n",
    "assert labels[-1] == 1\n",
    "\n",
    "labels[:5], labels[-5:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we are left with are two matrices (numpy arrays): one of the features, and one of the labels. The labels matrix is one-dimensional. \n",
    "\n",
    "The features matrix can be have another dimension of any length. However, along its first dimension, it must be the same shape as the labels, so that the two align (every feature has a label).\n",
    "\n",
    "So, here, our labels are `n x 1`, and our features are `256 x n`. (Our sampling rate is 512, producing power spectra of size 256.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert features.shape[0] == labels.shape[0]\n",
    "\n",
    "features.shape, labels.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By convention, we refer to the features matrix as `X` and the labels matrix as `y`. \n",
    "\n",
    "(I know, I know. Mathematicians have no idea how to name a variable. But, the variable names `X` and `y` are so common in literature and documentation, I feel you should be familiar with them. And thus I perpetuate the problem...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = features\n",
    "y = labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can finally train our classifier!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.fit(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, we trained it! So, how does it do?\n",
    "\n",
    "Aha, but if we've used all of our data to trian the classifier, what can we use to *test* the classifier?\n",
    "\n",
    "We could test the classifier on the exact same stuff we trained it on, but *WE SHOULD NEVER DO THIS!!! NEVER TEST ON THE TRAIN SET!!* \n",
    "\n",
    "Why? The answer has to do with a concept called *overfitting*. \n",
    "\n",
    "> A model *overfits* to its data when it explains the data it was trained on *so exactly* that it no longer generalizes to new, unseen data. \n",
    "\n",
    "For example, a model may memorize the data, including noise in the data, to produce perfect accuracy. Effectively, it will memorize the data! How useful will this memorization be in classifying data it's never seen before? Not very!\n",
    "\n",
    "To get around this, we *split* the data into train and test sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(\n",
    "    X, y, \n",
    "    # We will withold 33% of the data for the test split.\n",
    "    test_size=0.33, \n",
    "    # We'll fix this parameter,\n",
    "    # so that we'll always produce the same\n",
    "    # \"random\" split on each run.* \n",
    "    random_state=42)\n",
    "\n",
    "# (Remember, there's nothing \"random\" on a computer,\n",
    "#  but, we do have *pseudo-*random number generators\n",
    "#  (PRNGs) which are almost as good. If we fix the seed,\n",
    "#  they will always generate the same sequence of numbers.\n",
    "#  However, to someone who doesn't know the seeds, the\n",
    "#  numbers look like white noise. Magic! You ever wonder\n",
    "#  how garage door openers work? Research more\n",
    "#  on your own time.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can train a *fresh* classifier on the train set,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = fresh_clf()\n",
    "clf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and test on the test set!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But XGBoost is even smarter than that. \n",
    "\n",
    "Oh yes, XGBoost is very smart. It can take some data, randomly split it into train and test sets, then see how well it does on the test set. It will use these results to tweak its parameters. Then, it can repeat this process a number of times, until it starts doing very well!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def xgb_cross_validate (\n",
    "    X: np.array,\n",
    "    y: np.array,\n",
    "    nfold: int=7\n",
    ") -> Tuple[XGBClassifier, pd.DataFrame]:\n",
    "    # eval_metrics:\n",
    "    # http://xgboost.readthedocs.io/en/latest//parameter.html\n",
    "    metrics = ['error@0.1', 'auc']\n",
    "#     metrics = [ 'auc' ]\n",
    "    # we use the @ syntax to override the default of 0.5 as the threshold for 0 / 1 classification\n",
    "    # the intent here to to minimize FAR at the expense of FRR\n",
    "    alg = fresh_clf()\n",
    "    xgtrain = xgb.DMatrix(X,y)\n",
    "    param = alg.get_xgb_params()\n",
    "    cvresults = xgb.cv(param,\n",
    "                      xgtrain,\n",
    "                      num_boost_round=alg.get_params()['n_estimators'],\n",
    "                      nfold=nfold,\n",
    "                      metrics=metrics,\n",
    "                      early_stopping_rounds=100\n",
    "                      )\n",
    "    alg.set_params(n_estimators=cvresults.shape[0])\n",
    "    alg.fit(X,y,eval_metric=metrics)\n",
    "    return alg, cvresults"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll now split the data into a *train* set (for which XGBoost will manage the train/test splitting), and a *validation* set, which we will use to see the results of XGBoost's efforts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train, X_validate, y_train, y_validate = sklearn.model_selection.train_test_split(\n",
    "    X, y, \n",
    "    test_size=0.33, \n",
    "    random_state=42)\n",
    "\n",
    "clf, cvres = xgb_cross_validate(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This `cvres` value is actually a DataFrame detailing each \"tweaking\" step XGBoost did in training the parameters. We can inspect the last few rows to see what it did right before it delievered the classifier to us.\n",
    "\n",
    "(If you notice in the `xgb_cross_validate` method above, we specified an `early_stopping_rounds`. We keep tweaking parametmers until the results are pretty high, or until we have tried for 100 rounds - whichever comes first)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvres.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.score(X_validate, y_validate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TODO\n",
    "\n",
    "For how many rounds did XGBoost tweak our classifier, above, before stopping?\n",
    "\n",
    "*Your answer here...*\n",
    "\n",
    "What is the estimated accuracy of our classifier at distinguishing between \"math\" and \"relax\" tasks, across all subjects?\n",
    "\n",
    "*Your answer here...*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TODO\n",
    "\n",
    "In this notebook, why did our classifier do as well as it did?\n",
    "\n",
    "In answering this question, do research, and make educated guesses based on what you learn. Optionally make arguments based on any any evidence you can collect. \n",
    "\n",
    "There are no right or wrong answers here. Nor is there only one answer - there are potentially a number of possible explanations! \n",
    "\n",
    "Some questions I always ask myself whenever I come up with explanations:\n",
    "\n",
    "- Could I test this theory? How?\n",
    "- Are there any *other* explanations? Can I test those?\n",
    "- Why might the current results be misleading?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here...*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# optionally, your code here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extra credit\n",
    "\n",
    "Discover something new, and describe (or show) it below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# your code, or words, here"
   ]
  }
 ],
 "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.6.1"
  },
  "name": "lab-1.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}