Example_hyperparam_opt.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples of hyperparameter optimization\n",
    "\n",
    "In this notebook, we provide examples for how to optimize hyperparameters of a decoder and then use the decoder with those hyperparameters. We demonstrate how to use 2 different hyperparameter optimization packages, [\"BayesianOptimization\"](https://github.com/fmfn/BayesianOptimization) and [\"hyperopt\"](http://hyperopt.github.io/hyperopt/). Both give similar performance. In the arXiv manuscript, I used \"BayesianOptimization\" (simply because I discovered it first).\n",
    " - The first few sections (1-3) just import packages, load the files, and preprocess them\n",
    " - Section 4 shows examples of [BayesianOptimization](https://github.com/fmfn/BayesianOptimization) for 3 decoders: Wiener Cascade, XGBoost, and Feedforward Neural Net\n",
    " - Section 5 shows examples of [hyperopt](http://hyperopt.github.io/hyperopt/) for 3 decoders: Wiener Cascade, XGBoost, and Feedforward Neural Net \n",
    " - Section 6 shows examples of making test-set predictions using the decoders with the fit hyperparameters\n",
    " \n",
    "Note that the example using the Wiener Cascade is quick, but the examples with XGBoost and the Feedforward Neural Net are slower (depending on your computer, potentially 10's of minutes)."
   ]
  },
  {
   "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": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using Theano backend.\n",
      "WARNING (theano.sandbox.cuda): The cuda backend is deprecated and will be removed in the next release (v0.10).  Please switch to the gpuarray backend. You can get more information about how to switch at this URL:\n",
      " https://github.com/Theano/Theano/wiki/Converting-to-the-new-gpu-back-end%28gpuarray%29\n",
      "\n",
      "Using gpu device 0: GeForce GTX TITAN X (CNMeM is enabled with initial size: 40.0% of memory, cuDNN Mixed dnn version. The header is from one version, but we link with a different version (5103, 5110))\n"
     ]
    }
   ],
   "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",
    "import time\n",
    "# If you would prefer to load the '.h5' example file rather than the '.pickle' example file. You need the deepdish package\n",
    "# import deepdish as dd \n",
    "\n",
    "#Import function to get the covariate matrix that includes spike history from previous bins\n",
    "from Neural_Decoding.preprocessing_funcs import get_spikes_with_history\n",
    "\n",
    "#Import metrics\n",
    "from Neural_Decoding.metrics import get_R2\n",
    "from Neural_Decoding.metrics import get_rho\n",
    "\n",
    "#Import decoder functions\n",
    "from Neural_Decoding.decoders import WienerCascadeDecoder\n",
    "from Neural_Decoding.decoders import WienerFilterDecoder\n",
    "from Neural_Decoding.decoders import DenseNNDecoder\n",
    "from Neural_Decoding.decoders import SimpleRNNDecoder\n",
    "from Neural_Decoding.decoders import GRUDecoder\n",
    "from Neural_Decoding.decoders import LSTMDecoder\n",
    "from Neural_Decoding.decoders import XGBoostDecoder\n",
    "from Neural_Decoding.decoders import SVRDecoder\n",
    "\n",
    "#Import hyperparameter optimization packages\n",
    "#If either are not installed, give a warning\n",
    "try:\n",
    "    from bayes_opt import BayesianOptimization\n",
    "except ImportError:\n",
    "    print(\"\\nWARNING: BayesianOptimization package is not installed. You will be unable to use section 4.\")\n",
    "    pass\n",
    "\n",
    "try:\n",
    "    from hyperopt import fmin, hp, Trials, tpe, STATUS_OK\n",
    "except ImportError:\n",
    "    print(\"\\nWARNING: hyperopt package is not installed. You will be unable to use section 5.\")\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Turn off deprecation warnings\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Load Data\n",
    "The data for this example can be downloaded at this [link](https://www.dropbox.com/sh/n4924ipcfjqc0t6/AACPWjxDKPEzQiXKUUFriFkJa?dl=0&preview=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": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# folder='' #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+'example_data_s1.pickle','rb') as f:\n",
    "#     neural_data,vels_binned=pickle.load(f,encoding='latin1') #If using python 3\n",
    "    neural_data,vels_binned=pickle.load(f) #If using python 2\n",
    "\n",
    "# #If you would prefer to load the '.h5' example file rather than the '.pickle' example file.\n",
    "# data=dd.io.load(folder+'example_data_s1.h5')\n",
    "# neural_data=data['neural_data']\n",
    "# vels_binned=data['vels_binned']"
   ]
  },
  {
   "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).\n",
    "I am using fewer bins in this example than in the manuscript and other examples, to make it run faster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bins_before=3 #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=3 #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": 5,
   "metadata": {
    "collapsed": false
   },
   "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": "markdown",
   "metadata": {},
   "source": [
    "#### Format Output Covariates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Set decoding output\n",
    "y=vels_binned"
   ]
  },
  {
   "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": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Set what part of data should be part of the training/testing/validation sets\n",
    "#I made the ranges smaller for this example so that the hyperparameter optimization runs faster\n",
    "training_range=[0.6, 0.7] \n",
    "testing_range=[0.7, 0.8]\n",
    "valid_range=[0.8,0.9]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Split Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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",
    "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",
    "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",
    "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": 9,
   "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",
    "#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": "markdown",
   "metadata": {},
   "source": [
    "## 4. Optimize Hyperparameters of decoders using \"BayesianOptimization\"\n",
    " - The general idea is that we will try to find the decoder hyperparameters that produce the highest R2 values on the validation set. \n",
    "\n",
    " - We will provide examples for a few decoders (Wiener Cascade, XGBoost, Feedforward Neural Net)\n",
    " \n",
    " A potential downside of BayesianOptimization is that it optimizes over a continuous space. So if a hyperparameter has integer values, the optimizer may unnecessarily test many nearby values (e.g. 2.05, 2.1, and 2.2) which are all treated the same (as 2), when it could just test the one integer value (2)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4A. Wiener Cascade (Linear Nonlinear Model)\n",
    " - The hyperparameter we are trying to optimize is \"degree\" (the degree of the polynomial).\n",
    " - Note that a sophisticated hyperparameter optimization technique is not needed for this decoder with a single hyperparameter - you could easily do a grid search. However, we show the example since it's the simplest and runs quickly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define a function that returns the metric we are trying to optimize (R2 value of the validation set) as a function of the hyperparameter (degree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def wc_evaluate(degree):\n",
    "    model_wc=WienerCascadeDecoder(degree) #Define model\n",
    "    model_wc.fit(X_flat_train,y_train) #Fit model\n",
    "    y_valid_predicted_wc=model_wc.predict(X_flat_valid) #Validation set predictions\n",
    "    return np.mean(get_R2(y_valid,y_valid_predicted_wc)) #R2 value of validation set (mean over x and y position/velocity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set range of hyperparameters, and run optimization\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Define Bayesian optimization, and set limits of hyperparameters \n",
    "#Here, we set the limit of \"degree\" to be [1, 6.99], so we test degrees 1,2,3,4,5,6\n",
    "wcBO = BayesianOptimization(wc_evaluate, {'degree': (1, 6.99)}, verbose=0)\n",
    "#Set number of initial runs (init_points) and subsequent tests (n_iter), and do the optimization\n",
    "#kappa is a parameter that sets exploration vs exploitation in the algorithm\n",
    "#We set kappa=10 (greater than the default) so there is more exploration when there are more hyperparameters\n",
    "wcBO.maximize(init_points=5, n_iter=5, kappa=10) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get best hyperparameters\n",
    "Note that you can also find out more information about each tested hyperparameter in \"wcBO.res\": (each hyperparameter tested and the resulting R2 value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'max_params': {'degree': 3.2266678392926771}, 'max_val': 0.72915601486578996}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Print out the best parameters and associated R2 value (called \"max_val\")\n",
    "wcBO.res['max']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Assign the best hyperparameter to a variable\n",
    "best_params=wcBO.res['max']['max_params'] \n",
    "degree=best_params['degree']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4B. XGBoost\n",
    "The hyperparameters we are trying to optimize are:\n",
    " - \"max_depth\" (maximum depth of the trees)\n",
    " - \"num_round\" (number of trees for fitting)\n",
    " - \"eta\" (learning rate)\n",
    " \n",
    "Note that this example can be somewhat slow (depending on your computer, potentially 10's of minutes)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define a function that returns the metric we are trying to optimize (R2 value of the validation set) as a function of the hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def xgb_evaluate(max_depth,num_round,eta):\n",
    "    #The parameters need to be in the correct format for the decoder, so we do that below\n",
    "    max_depth=int(max_depth) \n",
    "    num_round=int(num_round) \n",
    "    eta=float(eta) \n",
    "    #Define model\n",
    "    model_xgb=XGBoostDecoder(max_depth=max_depth, num_round=num_round, eta=eta) \n",
    "    model_xgb.fit(X_flat_train,y_train) #Fit model\n",
    "    y_valid_predicted_xgb=model_xgb.predict(X_flat_valid) #Get validation set predictions\n",
    "    return np.mean(get_R2(y_valid,y_valid_predicted_xgb)) #Return mean validation set R2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set range of hyperparameters, and run optimization\n",
    "If you want to keep track of progress, set verbose=1 in the cell below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Do bayesian optimization, and set limits of hyperparameters\n",
    "xgbBO = BayesianOptimization(xgb_evaluate, {'max_depth': (2, 6.99), 'num_round': (100,600.99), 'eta': (0.01, 0.8)},verbose=0) #Define Bayesian optimization, and set limits of hyperparameters    \n",
    "#Set number of initial runs and subsequent tests, and do the optimization. Also, we set kappa=10 (greater than the default) so there is more exploration when there are more hyperparameters\n",
    "xgbBO.maximize(init_points=10, n_iter=10, kappa=10) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get best hyperparameters\n",
    "Note that you can also find out more information about each tested hyperparameter in \"xgbBO.res\": (each hyperparameter tested and the resulting R2 value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'max_params': {'eta': 0.08703121928968556,\n",
       "  'max_depth': 5.119221012479418,\n",
       "  'num_round': 302.98790346950642},\n",
       " 'max_val': 0.7746865181664323}"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Print out the best parameters and associated R2 value (called \"max_val\")\n",
    "xgbBO.res['max']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Assign the best hyperparameters to variables, and put them in the correct format\n",
    "best_params=xgbBO.res['max']['max_params'] #Get the hyperparameters that give rise to the best fit\n",
    "num_round=np.int(best_params['num_round']) #We want the integer value associated with the best \"num_round\" parameter (which is what the xgb_evaluate function does above)\n",
    "max_depth=np.int(best_params['max_depth']) #We want the integer value associated with the best \"max_depth\" parameter (which is what the xgb_evaluate function does above)\n",
    "eta=best_params['eta']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4C. Feedforward (Dense) Neural Net\n",
    "The hyperparameters we are trying to optimize are:\n",
    " - \"num_units\" (the number of hidden units in each layer)\n",
    " - \"frac_dropout\" (the proportion of units that are dropped out\"\n",
    " - \"n_epochs\" (the number of epochs used for fitting)\n",
    "\n",
    "Note that this example can be somewhat slow (depending on your computer, potentially 10's of minutes)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define a function that returns the metric we are trying to optimize (R2 value of the validation set) as a function of the hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dnn_evaluate(num_units,frac_dropout,n_epochs):\n",
    "    #The parameters need to be in the correct format for the decoder, so we do that below\n",
    "    num_units=int(num_units)\n",
    "    frac_dropout=float(frac_dropout)\n",
    "    n_epochs=int(n_epochs)\n",
    "    #Declare and fit decoder\n",
    "    model_dnn=DenseNNDecoder(units=[num_units,num_units],dropout=frac_dropout,num_epochs=n_epochs)\n",
    "    model_dnn.fit(X_flat_train,y_train)\n",
    "    #Make predictions and get R2 values on validation set\n",
    "    y_valid_predicted_dnn=model_dnn.predict(X_flat_valid)\n",
    "    return np.mean(get_R2(y_valid,y_valid_predicted_dnn))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set range of hyperparameters, and run optimization\n",
    "If you want to keep track of progress, set verbose=1 in the cell below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda/anaconda2/lib/python2.7/site-packages/keras/models.py:826: UserWarning: The `nb_epoch` argument in `fit` has been renamed `epochs`.\n",
      "  warnings.warn('The `nb_epoch` argument in `fit` '\n"
     ]
    }
   ],
   "source": [
    "#Do bayesian optimization, and set limits of hyperparameters\n",
    "dnnBO = BayesianOptimization(dnn_evaluate, {'num_units': (50, 700.99), 'frac_dropout': (0,.5), 'n_epochs': (2,15.99)},verbose=0)\n",
    "\n",
    "#Set number of initial runs (init_points) and subsequent tests (n_iter), and do the optimization\n",
    "#kappa is a parameter that sets exploration vs exploitation in the algorithm - 10 seems to work pretty welldnnBO = BayesianOptimization(dnn_evaluate, {'num_units': (50, 500), 'frac_dropout': (0.,.5), 'n_epochs': (2,15)})\n",
    "dnnBO.maximize(init_points=10, n_iter=10, kappa=10)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get best hyperparameters\n",
    "Note that you can also find out more information about each tested hyperparameter in \"dnnBO.res\": (each hyperparameter tested and the resulting R2 value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'max_params': {'frac_dropout': 0.26423316372026218,\n",
       "  'n_epochs': 11.021987738367251,\n",
       "  'num_units': 357.79446839049604},\n",
       " 'max_val': 0.80857182708375452}"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Print out the best parameters and associated R2 value\n",
    "dnnBO.res['max']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Assign the best hyperparameters to variables, and put them in the correct format\n",
    "best_params=dnnBO.res['max']['max_params']\n",
    "frac_dropout=float(best_params['frac_dropout'])\n",
    "n_epochs=np.int(best_params['n_epochs'])\n",
    "num_units=np.int(best_params['num_units'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Optimize Hyperparameters of decoders using \"Hyperopt\n",
    "\n",
    " - The general idea is that we will try to find the decoder hyperparameters that produce the highest R2 values on the validation set. \n",
    "\n",
    " - We will provide examples for a few decoders (Wiener Cascade, XGBoost, Feedforward Neural Net)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5A. Wiener Cascade\n",
    " - The hyperparameter we are trying to optimize is \"degree\" (the degree of the polynomial).\n",
    " - Note that a sophisticated hyperparameter optimization technique is not needed for this decoder with a single hyperparameter - you could easily do a grid search. However, we show the example since it's the simplest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define a function that returns the metric we are trying to optimize (R2 value of the validation set) as a function of the hyperparameter (degree)\n",
    " - hyperopt minimizes the parameter, so we will return -R2 (in order to maximize R2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def wc_evaluate2(degree):\n",
    "    \n",
    "    model_wc=WienerCascadeDecoder(degree) #Define model\n",
    "    model_wc.fit(X_flat_train,y_train) #Fit model\n",
    "    y_valid_predicted_wc=model_wc.predict(X_flat_valid) #Validation set predictions\n",
    "    return -np.mean(get_R2(y_valid,y_valid_predicted_wc)) #-R2 value of validation set (mean over x and y position/velocity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set range of hyperparameters, and run optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#The range of values I'll look at for the parameter\n",
    "#\"hp.quniform\" will allow us to look at integer (rather than continuously spaced) values.\n",
    "#Below we consider values of \"degree\" starting at 1, going until 6, and spaced at values of 1 (i.e., 1,2,3,4,5,6)\n",
    "space = hp.quniform('degree', 1, 6, 1)\n",
    "\n",
    "#object that holds iteration results\n",
    "trials = Trials()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Do optimization\n",
    "#Set the number of evaluations below (10 in this example)\n",
    "hyperoptBest = fmin(wc_evaluate2, space, algo=tpe.suggest, max_evals=10, trials=trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get best hyperparameters\n",
    "Note that you can also find out more information about each tested hyperparameter in the \"trials\" object. \"trials.results\" will give the R2 value for each hyperparameters tested, and \"trials.vals\" will give you the values of the hyperparameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('R2_validation', 0.7321137122470418)\n"
     ]
    }
   ],
   "source": [
    "print(\"R2_validation\",-trials.best_trial['result']['loss'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'degree': 5.0}\n"
     ]
    }
   ],
   "source": [
    "print(hyperoptBest)\n",
    "degree=hyperoptBest['degree']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5B. XGBoost\n",
    "The hyperparameters we are trying to optimize are:\n",
    " - \"max_depth\" (maximum depth of the trees)\n",
    " - \"num_round\" (number of trees for fitting)\n",
    " - \"eta\" (learning rate)\n",
    " \n",
    "Note that this example can be somewhat slow (depending on your computer, potentially 10's of minutes)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define a function that returns the metric we are trying to optimize (R2 value of the validation set) as a function of the hyperparameter (degree)\n",
    " - hyperopt minimizes the parameter, so we will return -R2 (in order to maximize R2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def xgb_evaluate2(params):\n",
    "    #Put parameters in correct formats\n",
    "    num_round=int(params['num_round'])\n",
    "    eta=float(params['eta'])\n",
    "    max_depth=int(params['max_depth'])\n",
    "    model_xgb=XGBoostDecoder(max_depth=max_depth, num_round=num_round, eta=eta) #Define model\n",
    "    model_xgb.fit(X_flat_train,y_train) #Fit model\n",
    "    y_valid_predicted_xgb=model_xgb.predict(X_flat_valid) #Get validation set predictions\n",
    "    return -np.mean(get_R2(y_valid,y_valid_predicted_xgb)) #Return mean validation set R2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set range of hyperparameters, and run optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#The range of values I'll look at for the parameter\n",
    "#\"hp.quniform\" will allow us to look at integer (rather than continuously spaced) values.\n",
    "#So for \"num_round\", we are looking at values between 100 and 600 by 50 (100,150,200,...600)\n",
    "#\"hp.uniform\" looks at continuously spaced values\n",
    "space = {\n",
    "    'eta': hp.uniform('eta', 0.01, 0.8),\n",
    "    'num_round': hp.quniform('num_round', 100,600,50),\n",
    "    'max_depth': hp.quniform('max_depth', 2,6,1),\n",
    "}\n",
    "\n",
    "#object that holds iteration results\n",
    "trials = Trials()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Do optimization\n",
    "#Set the number of evaluations below (20 in this example)\n",
    "hyperoptBest = fmin(xgb_evaluate2, space, algo=tpe.suggest, max_evals=20, trials=trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get best hyperparameters\n",
    "Note that you can also find out more information about each tested hyperparameter in the \"trials\" object. \"trials.results\" will give the R2 value for each hyperparameters tested, and \"trials.vals\" will give you the values of the hyperparameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('R2_validation', 0.7716811700071836)\n"
     ]
    }
   ],
   "source": [
    "print(\"R2_validation\",-trials.best_trial['result']['loss'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'num_round': 550.0, 'eta': 0.11823310650011055, 'max_depth': 5.0}\n"
     ]
    }
   ],
   "source": [
    "print(hyperoptBest)\n",
    "\n",
    "best_params=hyperoptBest #Just renamed so it was in the same format as I used with BayesOptimization\n",
    "num_round=np.int(best_params['num_round']) #We want the integer value associated with the best \"num_round\" parameter (which is what the xgb_evaluate function does above)\n",
    "max_depth=np.int(best_params['max_depth']) #We want the integer value associated with the best \"max_depth\" parameter (which is what the xgb_evaluate function does above)\n",
    "eta=best_params['eta']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5C. Feedforward (Dense) Neural Net\n",
    "The hyperparameters we are trying to optimize are:\n",
    " - \"num_units\" (the number of hidden units in each layer)\n",
    " - \"frac_dropout\" (the proportion of units that are dropped out\"\n",
    " - \"n_epochs\" (the number of epochs used for fitting)\n",
    "\n",
    "Note that this example can be somewhat slow (depending on your computer, potentially 10's of minutes)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define a function that returns the metric we are trying to optimize (R2 value of the validation set) as a function of the hyperparameter (degree)\n",
    " - hyperopt minimizes the parameter, so we will return -R2 (in order to maximize R2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dnn_evaluate2(params):\n",
    "    #Put parameters in proper format\n",
    "    num_units=int(params['num_units'])\n",
    "    frac_dropout=float(params['frac_dropout'])\n",
    "    n_epochs=int(params['n_epochs'])\n",
    "    model_dnn=DenseNNDecoder(units=[num_units,num_units],dropout=frac_dropout,num_epochs=n_epochs) #Define model\n",
    "    model_dnn.fit(X_flat_train,y_train) #Fit model\n",
    "    y_valid_predicted_dnn=model_dnn.predict(X_flat_valid) #Get validation set predictions\n",
    "    return -np.mean(get_R2(y_valid,y_valid_predicted_dnn)) #Return -R2 value of validation set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set range of hyperparameters, and run optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#The range of values I'll look at for the parameter\n",
    "#\"hp.quniform\" will allow us to look at integer (rather than continuously spaced) values.\n",
    "#So for \"num_units\", we are looking at values between 50 and 700 by 10 (50,60,70,...700)\n",
    "#\"hp.uniform\" looks at continuously spaced values\n",
    "space = {\n",
    "    'frac_dropout': hp.uniform('frac_dropout', 0., 0.5),\n",
    "    'num_units': hp.quniform('num_units', 50,700,10),\n",
    "    'n_epochs': hp.quniform('n_epochs', 2,15,1),\n",
    "}\n",
    "\n",
    "#object that holds iteration results\n",
    "trials = Trials()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Do optimization\n",
    "#Set the number of evaluations below (20 in this example)\n",
    "hyperoptBest = fmin(dnn_evaluate2, space, algo=tpe.suggest, max_evals=20, trials=trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get best hyperparameters\n",
    "Note that you can also find out more information about each tested hyperparameter in the \"trials\" object. \"trials.results\" will give the R2 value for each hyperparameters tested, and \"trials.vals\" will give you the values of the hyperparameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('R2_validation', 0.8082196200199168)\n"
     ]
    }
   ],
   "source": [
    "print(\"R2_validation\",-trials.best_trial['result']['loss'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'frac_dropout': 0.35466311053401717, 'num_units': 550.0, 'n_epochs': 11.0}\n"
     ]
    }
   ],
   "source": [
    "print(hyperoptBest)\n",
    "\n",
    "best_params=hyperoptBest #Just renamed so it was in the same format as I used with BayesOptimization\n",
    "frac_dropout=float(best_params['frac_dropout'])\n",
    "n_epochs=np.int(best_params['n_epochs'])\n",
    "num_units=np.int(best_params['num_units'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Use the optimal hyperparameters to fit the decoder on the test set\n",
    "This can be run after running either section 4 or section 5 (both don't need to be run).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6A. Wiener Cascade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('R2s_wc:', array([ 0.65910948,  0.71765739]))\n"
     ]
    }
   ],
   "source": [
    "#\"degree\" was determined during hyperparameter optimization\n",
    "model_wc=WienerCascadeDecoder(degree) #Declare model w/ fit hyperparameter\n",
    "model_wc.fit(X_flat_train,y_train) #Fit model on training data\n",
    "y_test_predicted_wc=model_wc.predict(X_flat_test) #Get test set predictions\n",
    "#Print R2 values on test set\n",
    "R2s_wc=get_R2(y_test,y_test_predicted_wc)\n",
    "print('R2s_wc:', R2s_wc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 6B. XGBoost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('R2s:', array([ 0.71761818,  0.7576485 ]))\n"
     ]
    }
   ],
   "source": [
    "# Run model w/ above hyperparameters\n",
    "model_xgb=XGBoostDecoder(max_depth=max_depth, num_round=num_round, eta=eta) #Declare model w/ fit hyperparameters\n",
    "model_xgb.fit(X_flat_train,y_train) #Fit model\n",
    "y_test_predicted_xgb=model_xgb.predict(X_flat_test) #Get test set predictions\n",
    "#Print R2 values on test set\n",
    "R2s_xgb=get_R2(y_test,y_test_predicted_xgb)\n",
    "print('R2s:', R2s_xgb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6C. Feedforward Neural Net"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('R2s:', array([ 0.76140749,  0.80849555]))\n"
     ]
    }
   ],
   "source": [
    "# Run model w/ above hyperparameters\n",
    "model_dnn=DenseNNDecoder(units=[num_units,num_units],dropout=frac_dropout,num_epochs=n_epochs) #Declare model w/ fit hyperparameters\n",
    "model_dnn.fit(X_flat_train,y_train) #Fit model\n",
    "y_test_predicted_dnn=model_dnn.predict(X_flat_test) #Get test set predictions\n",
    "#Print R2 values on test set\n",
    "R2s_dnn=get_R2(y_test,y_test_predicted_dnn)\n",
    "print('R2s:', R2s_dnn)    "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [py35]",
   "language": "python",
   "name": "Python [py35]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}