Data Preparation.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:2929e46ac1f87fde21eadf599b8af51b0dd6fe3bcdd7a0c316cf5ac728b82504"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os, sys, pickle, itertools\n",
      "\n",
      "import joblib\n",
      "\n",
      "from scipy.io import loadmat\n",
      "from scipy.signal import resample, butter, lfilter\n",
      "\n",
      "import sklearn.linear_model\n",
      "import sklearn.cross_validation\n",
      "import sklearn.metrics"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def butter_bandpass(lowcut, highcut, fs, order=5):\n",
      "    nyq = 0.5 * fs\n",
      "    low = lowcut / nyq\n",
      "    high = highcut / nyq\n",
      "    b, a = butter(order, [low, high], btype='band')\n",
      "    return b, a\n",
      "\n",
      "def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):\n",
      "    b, a = butter_bandpass(lowcut, highcut, fs, order=order)\n",
      "    y = lfilter(b, a, data)\n",
      "    return y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Prepare 8 overlapping log-spaced bandpass filters between 5 and 200Hz\n",
      "filters = []\n",
      "for i in range(8):\n",
      "    x = logspace(log10(5),log10(200),10)\n",
      "    filters.append(butter_bandpass(x[i], floor(x[i+2]),400,3))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# List of all possible combinations of 4 or fewer filters\n",
      "wis = list(itertools.chain(itertools.combinations(range(8),1),itertools.combinations(range(8),2),itertools.combinations(range(8),3),itertools.combinations(range(8),4)))\n",
      "\n",
      "# Labels to read and write data for each subject\n",
      "inlabels = ['Dog_1','Dog_2','Dog_3','Dog_4','Patient_1','Patient_2','Patient_3','Patient_4','Patient_5','Patient_6','Patient_7','Patient_8']\n",
      "outlabels = ['dog_0', 'dog_1', 'dog_2', 'dog_3', 'patient_0', 'patient_1', 'patient_2', 'patient_3', 'patient_4',\n",
      "      'patient_5', 'patient_6', 'patient_7']\n",
      "\n",
      "def process(subject):\n",
      "    \"\"\"\n",
      "        Given a subject index, filter each clip, calculate covariance matrices, \n",
      "        calculate top-3 filter sets, and save processed data to pickle.\n",
      "    \"\"\"\n",
      "    dn = './data/clips/%s/'%inlabels[subject]\n",
      "    fns = [fn for fn in os.listdir(dn) if '.mat' in fn]\n",
      "\n",
      "    allcovs = []\n",
      "    labels = []\n",
      "\n",
      "    print dn\n",
      "    # For each clip, resample to 400Hz, apply each filter, calculate and normalize covariance.\n",
      "    for fn in fns:\n",
      "        covs = []   \n",
      "        m = loadmat(dn+fn)\n",
      "        d = m['data']\n",
      "        d = resample(d, 400, axis=1)\n",
      "        if 'inter' in fn:\n",
      "            l = 0\n",
      "        elif '_ictal' in fn:\n",
      "            l = 1\n",
      "        else:\n",
      "            l = -1\n",
      "\n",
      "        labels.append(l)\n",
      "\n",
      "        for b, a in filters:\n",
      "            f = lfilter(b,a,d)\n",
      "            c = cov(f)\n",
      "            c = (c-c.mean())/c.std()\n",
      "            covs.append(c)\n",
      "        allcovs.append(covs)\n",
      "    allcovs = array(allcovs)\n",
      "    labels = array(labels)\n",
      "    \n",
      "    # For each filter combination, test prediction quality by CV of logistic regression.\n",
      "    scores = []\n",
      "    for w in wis:\n",
      "        y = labels[labels != -1]\n",
      "        X = allcovs[labels != -1]\n",
      "        X = X[:,w,::2,::2]\n",
      "        X = X.reshape((X.shape[0],-1))\n",
      "\n",
      "        ps = []\n",
      "        test_size = 0.25\n",
      "        for tri, tei in sklearn.cross_validation.ShuffleSplit(X.shape[0], n_iter=15, test_size=test_size, random_state=42):\n",
      "            X_train = X[tri]\n",
      "            X_test = X[tei]\n",
      "            y_train = y[tri]\n",
      "            y_test = y[tei]\n",
      "\n",
      "            clf = sklearn.linear_model.SGDClassifier(loss='log', penalty='l1', alpha=0.0001)\n",
      "            clf.fit(X_train, y_train)\n",
      "            p = clf.predict_proba(X_test)[:,1]\n",
      "            cv = sklearn.metrics.roc_auc_score(y_test, p)\n",
      "            ps.append(cv)\n",
      "        ps = array(ps)\n",
      "        scores.append(ps.mean())    \n",
      "\n",
      "    # Select 3 best filter sets and save processed features and labels to pickle.\n",
      "    best = sorted(zip(scores,wis))[-3:]\n",
      "    sets = 'ABC'\n",
      "    i = 0\n",
      "    for cv, w in best:\n",
      "        print outlabels[subject], cv, w\n",
      "        y = labels\n",
      "        X = allcovs\n",
      "        X = X[:,w,:,:]\n",
      "        d = {'y':y, 'covs':X, 'w':w, 'cv':cv, 'fns':fns}\n",
      "        pickle.dump(d, open('./data/cov_opt_%s_%s.pickle'%(outlabels[subject],sets[i]), 'w'))\n",
      "        i += 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Process all subjects in parallel\n",
      "# reduce n_jobs if out of memory\n",
      "r = joblib.Parallel(n_jobs=-1)(joblib.delayed(process)(n) for n in range(12))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}