filterbank_gradient_boosting_standard_scaler.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "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.ensemble import GradientBoostingClassifier\n",
    "from sklearn.preprocessing import StandardScaler\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": 15,
   "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-6*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": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "scaler= StandardScaler()\n",
    "\n",
    "def data_preprocess_train(X):\n",
    "    X_prep = scaler.fit_transform(X)\n",
    "    return X_prep\n",
    "\n",
    "def data_preprocess_test(X):\n",
    "    X_prep = scaler.transform(X)\n",
    "    return X_prep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "subjects = range(1,13)\n",
    "ids_tot = []\n",
    "pred_tot = []\n",
    "y_tot = []\n",
    "\n",
    "# design a butterworth bandpass filter bank\n",
    "freqs = []\n",
    "for i in range(4, 28, 2):\n",
    "    freqs.append([i, i+4])\n",
    "    \n",
    "filter_bank_length = len(freqs)\n",
    "\n",
    "b = []\n",
    "a = []\n",
    "for i in range(filter_bank_length):\n",
    "    numer, denom = butter(5,np.array(freqs[i])/250.0,btype='bandpass')\n",
    "    b.append(numer)\n",
    "    a.append(denom)\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 = 1\n",
    "\n",
    "# submission file\n",
    "submission_file = 'filterbank_gradient_boosting_standard_scaler.csv'\n",
    "cols = ['HandStart','FirstDigitTouch',\n",
    "        'BothStartLoadPhase','LiftOff',\n",
    "        'Replace','BothReleased']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "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",
    "    feattr_tot = []\n",
    "    featte_tot = []\n",
    "    \n",
    "    ################ READ DATA ################################################\n",
    "    fnames =  glob('train/subj%d_series*_data.csv' % (subject))\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",
    "    for i in range(filter_bank_length):\n",
    "        # y labels will be rebuilt during each iteration, last will be final\n",
    "        y = []\n",
    "        \n",
    "        raw._data[picks] = lfilter(b[i],a[i],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",
    "        feattr_tot.append(feattr)\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",
    "    # vertically stack filter bank features together\n",
    "    feattr = np.concatenate(feattr_tot)\n",
    "    \n",
    "    ################ Create test Features #####################################\n",
    "    # read test data \n",
    "    fnames =  glob('test/subj%d_series*_data.csv' % (subject))\n",
    "    raw = concatenate_raws([creat_mne_raw_object(fname, read_events=False) for fname in fnames])\n",
    "    \n",
    "    for i in range(filter_bank_length):\n",
    "        raw._data[picks] = lfilter(b[i],a[i],raw._data[picks])\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",
    "        featte_tot.append(featte)\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",
    "    # vertically stack filter bank test features together\n",
    "    featte = np.concatenate(featte_tot)\n",
    "    \n",
    "    ################ Train classifiers ########################################\n",
    "    clf = GradientBoostingClassifier()\n",
    "    pred = np.empty((len(ids),6))\n",
    "    for i in range(6):\n",
    "        print('Train subject %d, class %s' % (subject, cols[i]))\n",
    "        \n",
    "        # Use standard scaler from train data\n",
    "        X_train=data_preprocess_train(feattr[:,::subsample].T)\n",
    "        X_test=data_preprocess_test(featte.T)\n",
    "        \n",
    "        clf.fit(X_train[::subsample,:],labels[i,::subsample])\n",
    "        pred[:,i] = clf.predict_proba(X_test)[:,1]\n",
    "    \n",
    "    pred_tot.append(pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create pandas object for sbmission\n",
    "submission = pd.DataFrame(index=np.concatenate(ids_tot),\n",
    "                          columns=cols,\n",
    "                          data=np.concatenate(pred_tot))\n",
    "\n",
    "# write file\n",
    "submission.to_csv(submission_file,index_label='id',float_format='%.3f')"
   ]
  },
  {
   "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
}