{"cells":[
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# -*- coding: utf-8 -*-\n\"\"\"\nCreated on Mon Jun 29 14:00:37 2015\n\n@author: alexandrebarachant\n\nBeat the benchmark with CSP and Logisitic regression.\n\nGeneral Idea :\n\nThe goal of this challenge is to detect events related to hand movements. Hand \nmovements are caracterized by change in signal power in the mu (~10Hz) and beta\n(~20Hz) frequency band over the sensorimotor cortex. CSP spatial filters are\ntrained to enhance signal comming from this brain area, instantaneous power is\nextracted and smoothed, and then feeded into a logisitic regression.\n\nPreprocessing :\n\nSignal are bandpass-filtered between 7 and 30 Hz to catch most of the signal of\ninterest. 4 CSP spatial filters are then applied to the signal, resutlting to\n4 new time series.  In order to train CSP spatial filters, EEG are epoched \nusing a window of 2 second before and after the event 'HandStart'. CSP training\nneeds two classes. the epochs after Replace event are assumed to contain \npatterns corresponding to hand movement, and epochs before are assumed to \ncontain resting state.\n\nFeature extraction :\n\nPreprocessing is applied, spatialy filtered signal are the rectified and \nconvolved with a 0.5 second rectangular window for smoothing. Then a logarithm\nis applied. the resutl is a vector of dimention 4 for each time sample.\n\nClassification :\n\nFor each of the 6 event type, a logistic regression is trained. For training \nonly, features are downsampled in oder to speed up the process. Prediction are\nthe probailities of the logistic regression.\n\n\"\"\"\n\nprint(__doc__)\n\nimport numpy as np\nimport pandas as pd\nfrom mne.io import RawArray\nfrom mne.channels import read_montage\nfrom mne.epochs import concatenate_epochs\nfrom mne import create_info, find_events, Epochs, concatenate_raws, pick_types\nfrom mne.decoding import CSP\n\nfrom sklearn.linear_model import LogisticRegression\nfrom glob import glob\n\nfrom scipy.signal import butter, lfilter, convolve, boxcar\nfrom joblib import Parallel, delayed\n\ndef 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\n\nsubjects = range(1,13)\nids_tot = []\npred_tot = []\n\n# design a butterworth bandpass filter \nfreqs = [7, 30]\nb,a = butter(5,np.array(freqs)/250.0,btype='bandpass')\n\n# CSP parameters\n# Number of spatial filter to use\nnfilters = 4\n\n# convolution\n# window for smoothing features\nnwin = 250\n\n# training subsample\nsubsample = 10\n\n# submission file\nsubmission_file = 'beat_the_benchmark.csv'\ncols = ['HandStart','FirstDigitTouch',\n        'BothStartLoadPhase','LiftOff',\n        'Replace','BothReleased']\n\nfor subject in subjects:\n    epochs_tot = []\n    y = []\n\n    ################ READ DATA ################################################\n    fnames =  glob('../input/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    # The function parallelized for speeding up the script\n    raw._data[picks] = np.array(Parallel(n_jobs=-1)(delayed(lfilter)(b,a,raw._data[i]) for i in picks))\n    \n    ################ CSP Filters training #####################################\n    # get event posision corresponding to HandStart\n    events = find_events(raw,stim_channel='HandStart', verbose=False)\n    # epochs signal for 2 second after the event\n    epochs = Epochs(raw, events, {'during' : 1}, 0, 2, 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 2 second before the event, this correspond to the \n    # rest period.\n    epochs_rest = Epochs(raw, events, {'before' : 1}, -2, 0, 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 #################################\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.array(Parallel(n_jobs=-1)(delayed(convolve)(feat[i],boxcar(nwin),'full') for i in range(nfilters)))\n    feattr = np.log(feattr[:,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 \n    fnames =  glob('../input/test/subj%d_series*_data.csv' % (subject))\n    raw = concatenate_raws([creat_mne_raw_object(fname, read_events=False) for fname in fnames])\n    raw._data[picks] = np.array(Parallel(n_jobs=-1)(delayed(lfilter)(b,a,raw._data[i]) for i in picks))\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.array(Parallel(n_jobs=-1)(delayed(convolve)(feat[i],boxcar(nwin),'full') for i in range(nfilters)))\n    featte = np.log(featte[:,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)\n\n# create pandas object for sbmission\nsubmission = pd.DataFrame(index=np.concatenate(ids_tot),\n                          columns=cols,\n                          data=np.concatenate(pred_tot))\n\n# write file\nsubmission.to_csv(submission_file,index_label='id',float_format='%.5f')        \n"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": ""
 }
],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}}, "nbformat": 4, "nbformat_minor": 0}