Examples_all_decoders-grat30:636380500192558468.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples of all decoders (except Kalman Filter)\n",
    "\n",
    "In this example notebook, we:\n",
    "1. Import the necessary packages\n",
    "2. Load a data file (spike trains and outputs we are predicting)\n",
    "3. Preprocess the data for use in all decoders\n",
    "4. Run all decoders and print the goodness of fit\n",
    "5. Plot example decoded outputs\n",
    "\n",
    "See \"Examples_kf_decoder\" for a Kalman filter example. <br>\n",
    "Because the Kalman filter utilizes different preprocessing, we don't include an example here. to keep this notebook more understandable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Import Packages\n",
    "\n",
    "Below, we import both standard packages, and functions from the accompanying .py files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Import standard packages\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy import io\n",
    "from scipy import stats\n",
    "import pickle\n",
    "\n",
    "#Import function to get the covariate matrix that includes spike history from previous bins\n",
    "from preprocessing_funcs import get_spikes_with_history\n",
    "\n",
    "#Import metrics\n",
    "from metrics import get_R2\n",
    "from metrics import get_rho\n",
    "\n",
    "#Import decoder functions\n",
    "from decoders import WienerCascadeDecoder\n",
    "from decoders import WienerFilterDecoder\n",
    "from decoders import DenseNNDecoder\n",
    "from decoders import SimpleRNNDecoder\n",
    "from decoders import GRUDecoder\n",
    "from decoders import LSTMDecoder\n",
    "from decoders import XGBoostDecoder\n",
    "from decoders import SVRDecoder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Load Data\n",
    "The data file for this example can be downloaded at https://dl.dropboxusercontent.com/u/2944301/Decoding_Data/example_data_s1.pickle. It was recorded by Raeed Chowdhury from Lee Miller's lab at Northwestern.\n",
    "\n",
    "\n",
    "The data that we load is in the format described below. We have another example notebook, \"Example_format_data\", that may be helpful towards putting the data in this format.\n",
    "\n",
    "Neural data should be a matrix of size \"number of time bins\" x \"number of neurons\", where each entry is the firing rate of a given neuron in a given time bin\n",
    "\n",
    "The output you are decoding should be a matrix of size \"number of time bins\" x \"number of features you are decoding\"\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "folder='/Users/guitchounts/Dropbox (coxlab)/Ephys/Data/GRat30/636380500192558468' #ENTER THE FOLDER THAT YOUR DATA IS IN\n",
    "# folder='/home/jglaser/Data/DecData/' \n",
    "# folder='/Users/jig289/Dropbox/Public/Decoding_Data/'\n",
    "\n",
    "with open(folder+'/jerk_10hz.pickle','rb') as f:\n",
    "#     neural_data,vels_binned=pickle.load(f,encoding='latin1') #If using python 3\n",
    "    jerk=pickle.load(f) #If using python 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "jerk = jerk.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'jerk' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-6761a849301d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mjerk\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'jerk' is not defined"
     ]
    }
   ],
   "source": [
    "jerk.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import h5py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = h5py.File(folder+'/lfp_spec_win1step01.mat','r')\n",
    "freqs = f['f'][:]\n",
    "time = f['t'][:]\n",
    "lfp_spec = f['lfp_spec'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lfp_spec.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_freq_idx(freqs,desired_freq): # make desired_freq a tuple, e.g. (0,4)\n",
    "    idx = []\n",
    "    for counter,value in enumerate(freqs):\n",
    "        if  desired_freq[0] <= value <= desired_freq[1]:\n",
    "            #yield counter\n",
    "            idx.append(counter)\n",
    "    return idx\n",
    "\n",
    "\n",
    "def get_power(spec,freq_range,freqs):\n",
    "\n",
    "   \n",
    "\n",
    "    idx = get_freq_idx(freqs,freq_range)\n",
    "        #idx_0_4 = get_freq_idx(freqs,(0,4))\n",
    "        #idx_5_30 = get_freq_idx(freqs,(5,30))\n",
    "        #idx_30_100 = get_freq_idx(freqs,(30,100))\n",
    "\n",
    "    power = np.mean(spec[idx,:],0)\n",
    "    #print 'shape of Pxx[idx,:] = ', Pxx[idx,:].shape\n",
    "    #power = np.trapz(spec[idx,:],axis=0)\n",
    "        #power_0_4 = np.mean(Pxx[idx_0_4,:],0)\n",
    "        #power_5_30 = np.mean(Pxx[idx_5_30,:],0) # take mean on 0-th dim of Pxx -> this gives power over time. \n",
    "        #power_30_100 = np.mean(Pxx[idx_30_100,:],0)\n",
    "\n",
    "\n",
    "    #print 'shape of power = ', power.shape\n",
    "\n",
    "    return power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "####### get LFP power for each channel in the four bands:\n",
    "\n",
    "lfp_power = np.zeros([64*4,time.shape[0]])  ## 64 channels x 4 bands\n",
    "\n",
    "for ch in range(64):\n",
    "    power_0_4 = get_power(lfp_spec[ch,:,:],[0,4],freqs)\n",
    "    power_5_15 = get_power(lfp_spec[ch,:,:],[5,15],freqs)\n",
    "    power_15_40 = get_power(lfp_spec[ch,:,:],[15,40],freqs)\n",
    "    power_40_100 = get_power(lfp_spec[ch,:,:],[40,100],freqs)\n",
    "    lfp_power[ch*4:(ch+1)*4,:] = power_0_4,power_5_15,power_15_40,power_40_100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "neural_data = lfp_power.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "neural_data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "jerk_power = np.mean(jerk_spec,axis=1).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "jerk_power.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(jerk_power[:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for ch in range(64):\n",
    "    power_0_4 = get_power(lfp_spec[ch,:,:],[0,4],freqs)\n",
    "    power_5_15 = get_power(lfp_spec[ch,:,:],[5,15],freqs)\n",
    "    power_15_40 = get_power(lfp_spec[ch,:,:],[15,40],freqs)\n",
    "    power_40_100 = get_power(lfp_spec[ch,:,:],[40,100],freqs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f, axarr = plt.subplots(3, sharex=True,dpi=600)\n",
    "axarr[0].plot(spike_time_vec,neural_data[:,0],linewidth=.25)\n",
    "axarr[1].plot(jerk_time,jerk_power[:,0],linewidth=.25)\n",
    "axarr[2].plot(spike_time_vec,raw_jerk[:,0],linewidth=.25)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lfp_spec.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(np.mean(lfp_spec,axis=(0,1)),jerk[:,0],alpha=.1,marker='o')\n",
    "#plt.ylim([-.002,.002])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "neural_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Preprocess Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3A. User Inputs\n",
    "The user can define what time period to use spikes from (with respect to the output)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bins_before=15 #How many bins of neural data prior to the output are used for decoding\n",
    "bins_current=1 #Whether to use concurrent time bin of neural data\n",
    "bins_after=15 #How many bins of neural data after the output are used for decoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3B. Format Covariates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Format Input Covariates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Format for recurrent neural networks (SimpleRNN, GRU, LSTM)\n",
    "# Function to get the covariate matrix that includes spike history from previous bins\n",
    "X=get_spikes_with_history(neural_data,bins_before,bins_after,bins_current)\n",
    "\n",
    "# Format for Wiener Filter, Wiener Cascade, XGBoost, and Dense Neural Network\n",
    "#Put in \"flat\" format, so each \"neuron / time\" is a single feature\n",
    "X_flat=X.reshape(X.shape[0],(X.shape[1]*X.shape[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_flat.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Format Output Covariates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Set decoding output\n",
    "#y=jerk_power\n",
    "y=jerk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3C. Split into training / testing / validation sets\n",
    "Note that hyperparameters should be determined using a separate validation set. \n",
    "Then, the goodness of fit should be be tested on a testing set (separate from the training and validation sets)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### User Options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Set what part of data should be part of the training/testing/validation sets\n",
    "training_range=[0, 0.5]\n",
    "testing_range=[0.7, 0.85]\n",
    "valid_range=[0.5,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Split Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "num_examples=X.shape[0]\n",
    "\n",
    "#Note that each range has a buffer of\"bins_before\" bins at the beginning, and \"bins_after\" bins at the end\n",
    "#This makes it so that the different sets don't include overlapping neural data\n",
    "training_set=np.arange(np.int(np.round(training_range[0]*num_examples))+bins_before,np.int(np.round(training_range[1]*num_examples))-bins_after)\n",
    "testing_set=np.arange(np.int(np.round(testing_range[0]*num_examples))+bins_before,np.int(np.round(testing_range[1]*num_examples))-bins_after)\n",
    "valid_set=np.arange(np.int(np.round(valid_range[0]*num_examples))+bins_before,np.int(np.round(valid_range[1]*num_examples))-bins_after)\n",
    "\n",
    "#Get training data\n",
    "X_train=X[training_set,:,:]\n",
    "X_flat_train=X_flat[training_set,:]\n",
    "\n",
    "\n",
    "y_train=y[training_set,:]\n",
    "\n",
    "#Get testing data\n",
    "X_test=X[testing_set,:,:]\n",
    "X_flat_test=X_flat[testing_set,:]\n",
    "\n",
    "y_test=y[testing_set,:]\n",
    "\n",
    "#Get validation data\n",
    "X_valid=X[valid_set,:,:]\n",
    "X_flat_valid=X_flat[valid_set,:]\n",
    "\n",
    "\n",
    "y_valid=y[valid_set,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3D. Process Covariates\n",
    "We normalize (z_score) the inputs and zero-center the outputs.\n",
    "Parameters for z-scoring (mean/std.) should be determined on the training set only, and then these z-scoring parameters are also used on the testing and validation sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Z-score \"X\" inputs. \n",
    "X_train_mean=np.nanmean(X_train,axis=0)\n",
    "X_train_std=np.nanstd(X_train,axis=0)\n",
    "X_train=(X_train-X_train_mean)/X_train_std\n",
    "X_test=(X_test-X_train_mean)/X_train_std\n",
    "X_valid=(X_valid-X_train_mean)/X_train_std\n",
    "\n",
    "\n",
    "#Z-score \"X_flat\" inputs. \n",
    "X_flat_train_mean=np.nanmean(X_flat_train,axis=0)\n",
    "X_flat_train_std=np.nanstd(X_flat_train,axis=0)\n",
    "X_flat_train=(X_flat_train-X_flat_train_mean)/X_flat_train_std\n",
    "X_flat_test=(X_flat_test-X_flat_train_mean)/X_flat_train_std\n",
    "X_flat_valid=(X_flat_valid-X_flat_train_mean)/X_flat_train_std\n",
    "\n",
    "#Zero-center outputs\n",
    "y_train_mean=np.mean(y_train,axis=0)\n",
    "y_train=y_train-y_train_mean\n",
    "y_test=y_test-y_train_mean\n",
    "y_valid=y_valid-y_train_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f, axarr = plt.subplots(2, sharex=True,dpi=600)\n",
    "axarr[0].plot(X_flat_train[0:1000,100],linewidth=.25)\n",
    "axarr[1].plot(y_train[0:1000,0],linewidth=.25)\n",
    "#axarr[2].plot(spike_time_vec,raw_jerk[:,0],linewidth=.25)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Run Decoders\n",
    "Note that in this example, we are evaluating the model fit on the validation set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4A. Wiener Filter (Linear Regression)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn import linear_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ridgeCV_model(train_x,train_y,test_x,test_y):\n",
    "    \n",
    "    model = linear_model.RidgeCV(alphas=[0.1,1.0,10.],normalize=True,fit_intercept=True)\n",
    "    model.fit(train_x,train_y)\n",
    "    prediction = model.predict(test_x)\n",
    "    score = model.score(test_x,test_y)\n",
    "    print 'Model score R^2 = ', score\n",
    "    plt.scatter(test_y,prediction,alpha=0.1,marker='o')\n",
    "    plt.axis('equal')\n",
    "    return model\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_flat_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_flat_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ridge_model = ridgeCV_model(X_flat_train,y_train,X_flat_test,y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Declare model\n",
    "model_wf=WienerFilterDecoder()\n",
    "\n",
    "#Fit model\n",
    "model_wf.fit(X_flat_train,y_train)\n",
    "\n",
    "#Get predictions\n",
    "y_valid_predicted_wf=model_wf.predict(X_flat_valid)\n",
    "\n",
    "#Get metric of fit\n",
    "R2s_wf=get_R2(y_valid,y_valid_predicted_wf)\n",
    "print('R2s:', R2s_wf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f, axarr = plt.subplots(2, sharex=True,dpi=600)\n",
    "axarr[0].plot(y_valid,linewidth=0.1)\n",
    "\n",
    "axarr[1].plot(y_valid_predicted_wf,linewidth=0.1,color='red')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(y_valid,y_valid_predicted_wf,alpha=0.1,marker='o')\n",
    "plt.axis('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4B. Wiener Cascade (Linear Nonlinear Model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Declare model\n",
    "model_wc=WienerCascadeDecoder(degree=3)\n",
    "\n",
    "#Fit model\n",
    "model_wc.fit(X_flat_train,y_train)\n",
    "\n",
    "#Get predictions\n",
    "y_valid_predicted_wc=model_wc.predict(X_flat_valid)\n",
    "\n",
    "#Get metric of fit\n",
    "R2s_wc=get_R2(y_valid,y_valid_predicted_wc)\n",
    "print('R2s:', R2s_wc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4C. XGBoost (Extreme Gradient Boosting)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Declare model\n",
    "model_xgb=XGBoostDecoder(max_depth=3,num_round=200,eta=0.3,gpu=-1) \n",
    "\n",
    "#Fit model\n",
    "model_xgb.fit(X_flat_train, y_train)\n",
    "\n",
    "#Get predictions\n",
    "y_valid_predicted_xgb=model_xgb.predict(X_flat_valid)\n",
    "\n",
    "#Get metric of fit\n",
    "R2s_xgb=get_R2(y_valid,y_valid_predicted_xgb)\n",
    "print('R2s:', R2s_xgb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4D. SVR (Support Vector Regression)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#The SVR works much better when the y values are normalized, so we first z-score the y values\n",
    "#They have previously been zero-centered, so we will just divide by the stdev (of the training set)\n",
    "y_train_std=np.nanstd(y_train,axis=0)\n",
    "y_zscore_train=y_train/y_train_std\n",
    "y_zscore_test=y_test/y_train_std\n",
    "y_zscore_valid=y_valid/y_train_std\n",
    "\n",
    "#Declare model\n",
    "model_svr=SVRDecoder(C=5, max_iter=4000)\n",
    "\n",
    "#Fit model\n",
    "model_svr.fit(X_flat_train,y_zscore_train)\n",
    "\n",
    "#Get predictions\n",
    "y_zscore_valid_predicted_svr=model_svr.predict(X_flat_valid)\n",
    "\n",
    "#Get metric of fit\n",
    "R2s_svr=get_R2(y_zscore_valid,y_zscore_valid_predicted_svr)\n",
    "print('R2s:', R2s_svr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4E. Dense Neural Network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Declare model\n",
    "model_dnn=DenseNNDecoder(units=400,dropout=0.25,num_epochs=10)\n",
    "\n",
    "#Fit model\n",
    "model_dnn.fit(X_flat_train,y_train)\n",
    "\n",
    "#Get predictions\n",
    "y_valid_predicted_dnn=model_dnn.predict(X_flat_valid)\n",
    "\n",
    "#Get metric of fit\n",
    "R2s_dnn=get_R2(y_valid,y_valid_predicted_dnn)\n",
    "print('R2s:', R2s_dnn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4F. Simple RNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Declare model\n",
    "model_rnn=SimpleRNNDecoder(units=400,dropout=0,num_epochs=5)\n",
    "\n",
    "#Fit model\n",
    "model_rnn.fit(X_train,y_train)\n",
    "\n",
    "#Get predictions\n",
    "y_valid_predicted_rnn=model_rnn.predict(X_valid)\n",
    "\n",
    "#Get metric of fit\n",
    "R2s_rnn=get_R2(y_valid,y_valid_predicted_rnn)\n",
    "print('R2s:', R2s_rnn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4G. GRU (Gated Recurrent Unit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Declare model\n",
    "model_gru=GRUDecoder(units=400,dropout=0,num_epochs=5)\n",
    "\n",
    "#Fit model\n",
    "model_gru.fit(X_train,y_train)\n",
    "\n",
    "#Get predictions\n",
    "y_valid_predicted_gru=model_gru.predict(X_valid)\n",
    "\n",
    "#Get metric of fit\n",
    "R2s_gru=get_R2(y_valid,y_valid_predicted_gru)\n",
    "print('R2s:', R2s_gru)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4H. LSTM (Long Short Term Memory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Declare model\n",
    "model_lstm=LSTMDecoder(units=400,dropout=0,num_epochs=5)\n",
    "\n",
    "#Fit model\n",
    "model_lstm.fit(X_train,y_train)\n",
    "\n",
    "#Get predictions\n",
    "y_valid_predicted_lstm=model_lstm.predict(X_valid)\n",
    "\n",
    "#Get metric of fit\n",
    "R2s_lstm=get_R2(y_valid,y_valid_predicted_lstm)\n",
    "print('R2s:', R2s_lstm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Make Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#As an example, I plot an example 1000 values of the x velocity (column index 0), both true and predicted with the Wiener filter\n",
    "#Note that I add back in the mean value, so that both true and predicted values are in the original coordinates\n",
    "fig_x_wf=plt.figure()\n",
    "plt.plot(y_valid[1000:2000,0]+y_train_mean[0],'b')\n",
    "plt.plot(y_valid_predicted_wf[1000:2000,0]+y_train_mean[0],'r')\n",
    "\n",
    "#Save figure\n",
    "# fig_x_wf.savefig('x_velocity_decoding.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}