split_train_test_pywavelets.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: MacOSX\n"
     ]
    }
   ],
   "source": [
    "%matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pywt\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from mne.io import RawArray\n",
    "from mne.channels import read_montage\n",
    "from mne.epochs import concatenate_epochs\n",
    "from mne import create_info, find_events, Epochs, concatenate_raws, pick_types\n",
    "from mne.decoding import CSP\n",
    "\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from glob import glob\n",
    "\n",
    "from scipy.signal import butter, lfilter, convolve, boxcar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def creat_mne_raw_object(fname,read_events=True):\n",
    "    \"\"\"Create a mne raw instance from csv file\"\"\"\n",
    "    # Read EEG file\n",
    "    data = pd.read_csv(fname)\n",
    "    \n",
    "    # get chanel names\n",
    "    ch_names = list(data.columns[1:])\n",
    "    \n",
    "    # read EEG standard montage from mne\n",
    "    montage = read_montage('standard_1005',ch_names)\n",
    "\n",
    "    ch_type = ['eeg']*len(ch_names)\n",
    "    data = 1e-3*np.array(data[ch_names]).T\n",
    "    \n",
    "    if read_events:\n",
    "        # events file\n",
    "        ev_fname = fname.replace('_data','_events')\n",
    "        # read event file\n",
    "        events = pd.read_csv(ev_fname)\n",
    "        events_names = events.columns[1:]\n",
    "        events_data = np.array(events[events_names]).T\n",
    "        \n",
    "        # define channel type, the first is EEG, the last 6 are stimulations\n",
    "        ch_type.extend(['stim']*6)\n",
    "        ch_names.extend(events_names)\n",
    "        # concatenate event file and data\n",
    "        data = np.concatenate((data,events_data))\n",
    "        \n",
    "    # create and populate MNE info structure\n",
    "    info = create_info(ch_names,sfreq=500.0, ch_types=ch_type, montage=montage)\n",
    "    info['filename'] = fname\n",
    "    \n",
    "    # create raw object \n",
    "    raw = RawArray(data,info,verbose=False)\n",
    "    \n",
    "    return raw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "subjects = range(1,13)\n",
    "ids_tot = []\n",
    "pred_tot = []\n",
    "y_tot = []\n",
    "\n",
    "# design a butterworth bandpass filter \n",
    "freqs = [7, 30]\n",
    "b,a = butter(5,np.array(freqs)/250.0,btype='bandpass')\n",
    "\n",
    "# CSP parameters\n",
    "# Number of spatial filter to use\n",
    "nfilters = 4\n",
    "\n",
    "# convolution\n",
    "# window for smoothing features\n",
    "nwin = 250\n",
    "\n",
    "# training subsample\n",
    "subsample = 10\n",
    "\n",
    "# submission file\n",
    "submission_file = 'split_train_test_pywavelets.csv'\n",
    "cols = ['HandStart','FirstDigitTouch',\n",
    "        'BothStartLoadPhase','LiftOff',\n",
    "        'Replace','BothReleased']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(192, 32, 501)\n",
      "(192, 32, 501)\n",
      "(192, 32, 501)\n",
      "(192, 32, 501)\n",
      "(191, 32, 501)\n",
      "(191, 32, 501)\n"
     ]
    }
   ],
   "source": [
    "subject = 1\n",
    "epochs_tot = []\n",
    "y = []\n",
    "\n",
    "################ READ DATA (FIRST SIX SERIES AS TRAINING DATA) ###########\n",
    "fnames = []\n",
    "for i in range(1, 7):\n",
    "    fnames.append('train/subj%d_series%d_data.csv' % (subject, i))\n",
    "\n",
    "# read and concatenate all the files\n",
    "raw = concatenate_raws([creat_mne_raw_object(fname) for fname in fnames])\n",
    "\n",
    "# pick eeg signal\n",
    "picks = pick_types(raw.info,eeg=True)\n",
    "\n",
    "# Filter data for alpha frequency and beta band\n",
    "# Note that MNE implement a zero phase (filtfilt) filtering not compatible\n",
    "# with the rule of future data.\n",
    "# Here we use left filter compatible with this constraint\n",
    "raw._data[picks] = lfilter(b,a,raw._data[picks])\n",
    "\n",
    "# determine start time/position of each event\n",
    "hand_start_events = find_events(raw,stim_channel='HandStart', verbose=False)\n",
    "first_digit_touch_events = find_events(raw,stim_channel='FirstDigitTouch', verbose=False)\n",
    "both_start_load_phase_events = find_events(raw,stim_channel='BothStartLoadPhase', verbose=False)\n",
    "lift_off_events = find_events(raw,stim_channel='LiftOff', verbose=False)\n",
    "replace_events = find_events(raw,stim_channel='Replace', verbose=False)\n",
    "both_released_events = find_events(raw,stim_channel='BothReleased', verbose=False)\n",
    "\n",
    "# epochs signal for 1.0 second after start of event\n",
    "hand_start_epochs = Epochs(raw, hand_start_events, {'HandStart' : 1}, tmin=0, tmax=1, proj=False, picks=picks, baseline=None, preload=True, add_eeg_ref=False, verbose=False)\n",
    "first_digit_touch_epochs = Epochs(raw, first_digit_touch_events, {'FirstDigitTouch' : 1}, tmin=0, tmax=1, proj=False, picks=picks, baseline=None, preload=True, add_eeg_ref=False, verbose=False)\n",
    "both_start_load_phase_epochs = Epochs(raw, both_start_load_phase_events, {'BothStartLoadPhase' : 1}, tmin=0, tmax=1, proj=False, picks=picks, baseline=None, preload=True, add_eeg_ref=False, verbose=False)\n",
    "lift_off_epochs = Epochs(raw, lift_off_events, {'LiftOff' : 1}, tmin=0, tmax=1, proj=False, picks=picks, baseline=None, preload=True, add_eeg_ref=False, verbose=False)\n",
    "replace_epochs = Epochs(raw, replace_events, {'Replace' : 1}, tmin=0, tmax=1, proj=False, picks=picks, baseline=None, preload=True, add_eeg_ref=False, verbose=False)\n",
    "both_released_epochs = Epochs(raw, both_released_events, {'BothReleased' : 1}, tmin=0, tmax=1, proj=False, picks=picks, baseline=None, preload=True, add_eeg_ref=False, verbose=False)\n",
    "\n",
    "hand_start_tensor = hand_start_epochs.get_data()\n",
    "first_digit_touch_tensor = first_digit_touch_epochs.get_data()\n",
    "both_start_load_phase_tensor = both_start_load_phase_epochs.get_data()\n",
    "lift_off_tensor = lift_off_epochs.get_data()\n",
    "replace_tensor = replace_epochs.get_data()\n",
    "both_released_tensor = both_released_epochs.get_data()\n",
    "\n",
    "print hand_start_tensor.shape\n",
    "print first_digit_touch_tensor.shape\n",
    "print both_start_load_phase_tensor.shape\n",
    "print lift_off_tensor.shape\n",
    "print replace_tensor.shape\n",
    "print both_released_tensor.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in range(32):\n",
    "    plt.plot(hand_start_tensor[0][i,:])\n",
    "plt.savefig('hand_start_plot_0.png', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(32):\n",
    "    plt.plot(first_digit_touch_tensor[0][i,:])\n",
    "plt.savefig('first_digit_touch_plot_0.png', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1150"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_instances = len(hand_start_tensor) + len(first_digit_touch_tensor) + len(both_start_load_phase_tensor) + len(lift_off_tensor) + len(replace_tensor) + len(both_released_tensor)\n",
    "total_instances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create feature vector by appending wavelet features (NOT COMPLETE)\n",
    "#feattr = np.empty((total_instances, ))\n",
    "curr = 0\n",
    "for m in hand_start_tensor:\n",
    "    for i in range(32):\n",
    "        wd = pywt.wavedec(m[i,:], 'db6', level=6)\n",
    "        wave_feats = map(lambda x: np.dot(x.T, x), wd)\n",
    "        #feattr[curr] = wave_feats\n",
    "        curr += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.84978878525210078,\n",
       " 0.17711777008117557,\n",
       " 0.69056067333006921,\n",
       " 1.1033921615109803,\n",
       " 0.14139449457167164,\n",
       " 0.0010081872178661863,\n",
       " 3.4548664661936582e-05]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wave_feats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.84978878525210089,\n",
       " 0.17711777008117557,\n",
       " 0.69056067333006921,\n",
       " 1.1033921615109803,\n",
       " 0.14139449457167164,\n",
       " 0.0010081872178661863,\n",
       " 3.4548664661936582e-05]"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "map(lambda x: np.dot(x.T, x) ,hand_start_wd1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = [np.array([1, 2, 3])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[14]"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "map(lambda x: np.dot(x.T, x) , a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train subject 1, class HandStart\n",
      "Train subject 1, class FirstDigitTouch\n",
      "Train subject 1, class BothStartLoadPhase\n",
      "Train subject 1, class LiftOff\n",
      "Train subject 1, class Replace\n",
      "Train subject 1, class BothReleased\n",
      "Train subject 2, class HandStart\n",
      "Train subject 2, class FirstDigitTouch\n",
      "Train subject 2, class BothStartLoadPhase\n",
      "Train subject 2, class LiftOff\n",
      "Train subject 2, class Replace\n",
      "Train subject 2, class BothReleased\n",
      "Train subject 3, class HandStart\n",
      "Train subject 3, class FirstDigitTouch\n",
      "Train subject 3, class BothStartLoadPhase\n",
      "Train subject 3, class LiftOff\n",
      "Train subject 3, class Replace\n",
      "Train subject 3, class BothReleased\n",
      "Train subject 4, class HandStart\n",
      "Train subject 4, class FirstDigitTouch\n",
      "Train subject 4, class BothStartLoadPhase\n",
      "Train subject 4, class LiftOff\n",
      "Train subject 4, class Replace\n",
      "Train subject 4, class BothReleased\n",
      "Train subject 5, class HandStart\n",
      "Train subject 5, class FirstDigitTouch\n",
      "Train subject 5, class BothStartLoadPhase\n",
      "Train subject 5, class LiftOff\n",
      "Train subject 5, class Replace\n",
      "Train subject 5, class BothReleased\n",
      "Train subject 6, class HandStart\n",
      "Train subject 6, class FirstDigitTouch\n",
      "Train subject 6, class BothStartLoadPhase\n",
      "Train subject 6, class LiftOff\n",
      "Train subject 6, class Replace\n",
      "Train subject 6, class BothReleased\n",
      "Train subject 7, class HandStart\n",
      "Train subject 7, class FirstDigitTouch\n",
      "Train subject 7, class BothStartLoadPhase\n",
      "Train subject 7, class LiftOff\n",
      "Train subject 7, class Replace\n",
      "Train subject 7, class BothReleased\n",
      "Train subject 8, class HandStart\n",
      "Train subject 8, class FirstDigitTouch\n",
      "Train subject 8, class BothStartLoadPhase\n",
      "Train subject 8, class LiftOff\n",
      "Train subject 8, class Replace\n",
      "Train subject 8, class BothReleased\n",
      "Train subject 9, class HandStart\n",
      "Train subject 9, class FirstDigitTouch\n",
      "Train subject 9, class BothStartLoadPhase\n",
      "Train subject 9, class LiftOff\n",
      "Train subject 9, class Replace\n",
      "Train subject 9, class BothReleased\n",
      "Train subject 10, class HandStart\n",
      "Train subject 10, class FirstDigitTouch\n",
      "Train subject 10, class BothStartLoadPhase\n",
      "Train subject 10, class LiftOff\n",
      "Train subject 10, class Replace\n",
      "Train subject 10, class BothReleased\n",
      "Train subject 11, class HandStart\n",
      "Train subject 11, class FirstDigitTouch\n",
      "Train subject 11, class BothStartLoadPhase\n",
      "Train subject 11, class LiftOff\n",
      "Train subject 11, class Replace\n",
      "Train subject 11, class BothReleased\n",
      "Train subject 12, class HandStart\n",
      "Train subject 12, class FirstDigitTouch\n",
      "Train subject 12, class BothStartLoadPhase\n",
      "Train subject 12, class LiftOff\n",
      "Train subject 12, class Replace\n",
      "Train subject 12, class BothReleased\n"
     ]
    }
   ],
   "source": [
    "for subject in subjects:\n",
    "    epochs_tot = []\n",
    "    y = []\n",
    "    \n",
    "    ################ READ DATA (FIRST SIX SERIES AS TRAINING DATA) ###########\n",
    "    fnames = []\n",
    "    for i in range(1, 7):\n",
    "        fnames.append('train/subj%d_series%d_data.csv' % (1, i))\n",
    "    \n",
    "    # read and concatenate all the files\n",
    "    raw = concatenate_raws([creat_mne_raw_object(fname) for fname in fnames])\n",
    "       \n",
    "    # pick eeg signal\n",
    "    picks = pick_types(raw.info,eeg=True)\n",
    "    \n",
    "    # Filter data for alpha frequency and beta band\n",
    "    # Note that MNE implement a zero phase (filtfilt) filtering not compatible\n",
    "    # with the rule of future data.\n",
    "    # Here we use left filter compatible with this constraint\n",
    "    raw._data[picks] = lfilter(b,a,raw._data[picks])\n",
    "    \n",
    "    ################ CSP Filters training #####################################\n",
    "    # get event posision corresponding to Replace\n",
    "    events = find_events(raw,stim_channel='Replace', verbose=False)\n",
    "    # epochs signal for 1.5 second before the movement\n",
    "    epochs = Epochs(raw, events, {'during' : 1}, -2, -0.5, proj=False,\n",
    "                    picks=picks, baseline=None, preload=True,\n",
    "                    add_eeg_ref=False, verbose=False)\n",
    "    \n",
    "    epochs_tot.append(epochs)\n",
    "    y.extend([1]*len(epochs))\n",
    "    \n",
    "    # epochs signal for 1.5 second after the movement, this correspond to the \n",
    "    # rest period.\n",
    "    epochs_rest = Epochs(raw, events, {'after' : 1}, 0.5, 2, proj=False,\n",
    "                    picks=picks, baseline=None, preload=True,\n",
    "                    add_eeg_ref=False, verbose=False)\n",
    "    \n",
    "    # Workaround to be able to concatenate epochs with MNE\n",
    "    epochs_rest.times = epochs.times\n",
    "    \n",
    "    y.extend([-1]*len(epochs_rest))\n",
    "    epochs_tot.append(epochs_rest)\n",
    "        \n",
    "    # Concatenate all epochs\n",
    "    epochs = concatenate_epochs(epochs_tot)\n",
    "    \n",
    "    # get data \n",
    "    X = epochs.get_data()\n",
    "    y = np.array(y)\n",
    "    \n",
    "    # train CSP\n",
    "    csp = CSP(n_components=nfilters, reg='lws')\n",
    "    csp.fit(X,y)\n",
    "    \n",
    "    ################ Create Training Features #################################x`\n",
    "    # apply csp filters and rectify signal\n",
    "    feat = np.dot(csp.filters_[0:nfilters],raw._data[picks])**2\n",
    "    \n",
    "    # smoothing by convolution with a rectangle window    \n",
    "    feattr = np.empty(feat.shape)\n",
    "    for i in range(nfilters):\n",
    "        feattr[i] = np.log(convolve(feat[i],boxcar(nwin),'full'))[0:feat.shape[1]]\n",
    "    \n",
    "    # training labels\n",
    "    # they are stored in the 6 last channels of the MNE raw object\n",
    "    labels = raw._data[32:]\n",
    "    \n",
    "    ################ Create test Features #####################################\n",
    "    # read test data (remaining series of train data)\n",
    "    fnames = []\n",
    "    for i in range(7, 9):\n",
    "        fnames.append('train/subj%d_series%d_data.csv' % (1, i))\n",
    "        \n",
    "    raw = concatenate_raws([creat_mne_raw_object(fname) for fname in fnames])\n",
    "    raw._data[picks] = lfilter(b,a,raw._data[picks])\n",
    "    \n",
    "    # get actual y values\n",
    "    y = raw._data[32:]\n",
    "    y_tot.append(y.T)\n",
    "    \n",
    "    # read ids\n",
    "    ids = np.concatenate([np.array(pd.read_csv(fname)['id']) for fname in fnames])\n",
    "    ids_tot.append(ids)\n",
    "    \n",
    "    # apply preprocessing on test data\n",
    "    feat = np.dot(csp.filters_[0:nfilters],raw._data[picks])**2\n",
    "    featte = np.empty(feat.shape)\n",
    "    for i in range(nfilters):\n",
    "        featte[i] = np.log(convolve(feat[i],boxcar(nwin),'full'))[0:feat.shape[1]]\n",
    "    \n",
    "    ################ Train classifiers ########################################\n",
    "    lr = LogisticRegression()\n",
    "    pred = np.empty((len(ids),6))\n",
    "    for i in range(6):\n",
    "        print('Train subject %d, class %s' % (subject, cols[i]))\n",
    "        lr.fit(feattr[:,::subsample].T,labels[i,::subsample])\n",
    "        pred[:,i] = lr.predict_proba(featte.T)[:,1]\n",
    "    \n",
    "    pred_tot.append(pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final AUC: 0.633816663314\n"
     ]
    }
   ],
   "source": [
    "y_score = np.concatenate(pred_tot)\n",
    "y_true = np.concatenate(y_tot)\n",
    "\n",
    "# Convert y_true to int matrix, auc calc will not work with floats\n",
    "y_true = y_true.astype(int)\n",
    "auc = roc_auc_score(y_true, y_score, average='macro')\n",
    "print 'Final AUC:', auc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}