IntroToDLwithKeras-SOS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An Introduction to Deep Learning with Keras\n",
    "\n",
    "This tutorial is meant to teach a beginning HEP under-graduate or graduate student who may be unfamiliar with python or data science in python to train Deep Learning models using Keras. Usually, tutorials attempt to familiarize you with a software package by leading you through the steps of a few tasks. Like most tutorials, there are many sections where you can simply follow the instructions and execute example code like a robot (usualy via copy/paste into a terminal, but in this case using Jupyter notebooks). But this tutorial also aims to teach you key concepts in scientific computing, Machine Learning, and Deep Learning in python through exercises that require you to slow down, think critically, and apply what you read. The exercises were derived from labs for a Deep Learning in HEP course for undergrads (taught at University of Texas at Arlington by Amir Farbin).  \n",
    "\n",
    "The tutorial is divided into three sections:\n",
    "\n",
    "A. <a href='#Basics'>Basics</a>\n",
    "   0. <a href='#Jupyter'>Jupyter</a>\n",
    "   1. <a href='#Python'>Python</a>\n",
    "   2. <a href='#Numpy'>Numpy</a>\n",
    "   3. <a href='#HDF5'>HDF5</a>\n",
    "   \n",
    "B. <a href='#MachineLearning'>MachineLearning</a>\n",
    "   1. <a href='#Dataset'>Dataset</a>\n",
    "   2. <a href='#Pandas'>Pandas</a>\n",
    "   3. <a href='#Scikit-learn'>Scikit-learn</a>\n",
    "\n",
    "C. <a href='#DeepLearning'>DeepLearning</a>\n",
    "   1. <a href='#Keras'>Keras</a>\n",
    "\n",
    "You are very likely to find that the first sections of this tutorial are very basic. If you have some familiarity with data science in python, we suggest you skip what you know and only follow sections B.1, B.2, and C. \n",
    "\n",
    "Please quickly skim the beginning sections to find the appropriate starting point for you. If there isn't sufficient time for you to finish the exercises during the tutorial session please read through the explanations, execute the cells containing the examples, and think about how you would solve the exercises. You can go back and try the exercises at home."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before you skip, please select the location where you are running this tutorial by uncommenting the appropriate line and run the following cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Please Select location of data by [un]commenting appropriate line\n",
    "\n",
    "# Running at CERN\n",
    "# basedir = \"/afs/cern.ch/user/a/afarbin/public/DLClass/\" \n",
    "\n",
    "# Running at UTA DL Cluster\n",
    "# basedir = \"/data/afarbin/DLClass/\" \n",
    "\n",
    "# Running local\n",
    "basedir = \"/Users/afarbin/Tutorial/DLClass/\"\n",
    "\n",
    "SUSY_filename= basedir + \"SUSY/SUSY.csv\"\n",
    "\n",
    "HiggsML_filename_train = basedir + \"HiggsML/training.csv\"\n",
    "HiggsML_filename_test = basedir + \"HiggsML/test.csv\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Basics'></a>\n",
    "# A. Basics\n",
    "\n",
    "Data Science in python usually starts with loading hdf5 files into numpy tensors for manipulation in an interactive python session. While you can run the session in a terminal, Jupyter provides a nice web-based alternative environment for data analysis. As you can see (since you should be running this in Jupyter Notebbok), it allows you to combine text, code, and results all in one interactive document in your browser. There are many excellent primers already out there for python, numpy, h5py, and jupyter. You are encouraged to study them on your own as needed.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Jupyter'></a>\n",
    "## 0. Jupyter\n",
    "If you are seeing this page, you have successfully connected to a python server via ssh tunnel and navigated to this notebook. Jupyter notebooks consist of cells that can hold text or code (usually python). This text that you are reading, was written into a text cell as simple text \"coding\" language known as mark-down. When this cell is run (either automatically at start of the notebook or manually by pressing shift-enter), the mark-down text is interpreted into nice looking text. Running a code cell will execute the code in that cell and give you the results. If you make a mistake, you can usually simply change the cell and re-run. But be aware that since you ran the mistaken cell already, whatever code was properly executed before your mistake/error, was already executed and has therefore changed your current python environment accordingly. In some cases this situation will be problematic, and you will need to rerun the notebook from the start by pressing the \"reload\" botton (next to the \"stop\" button) above.\n",
    "\n",
    "You are encouraged to add cells to this notebook (using the \"+\" button on the tool bar) and play around a bit. If you don't want to mess up this notebook, you can work in a copy of this notebook by selecting Make Copy from the File menu."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Python'></a>\n",
    "## 1. Python\n",
    "\n",
    "Here we are assuming you have some basic level of python knowledge, such as the syntax. There are many great python tutorials available. For an introductory level interactive tutorial you can try this one: http://www.learnpython.org/\n",
    "\n",
    "\n",
    "We will lead you through exercises that show you common fundamental problems you might face when doing a deep learning problem. Lets start with generating some fake random data. You can get a random number between 0 and 1 using the python random module as follow: (remember to execute this in Jupyter notebooks click in the cell and hit shift-enter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "x=random.random()\n",
    "print \"The Value of x is\", x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise A.1.1\n",
    "Using the random method (shown above), write a function GenerateData(N, mymin, mymax), that returns a python list containing N random numbers between a specified minimum and maximum value. Note that you may want to quickly work out on paper how to turn numbers between 0 and 1 to other values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Skeleton\n",
    "def GenerateData(N,min,max):\n",
    "    out = []\n",
    "    ### BEGIN SOLUTION\n",
    "\n",
    "    # Fill in your solution here \n",
    "    \n",
    "    ### END SOLUTION\n",
    "    return out\n",
    "\n",
    "Data=GenerateData(1000,-10,10)\n",
    "print \"Data Type:\", type(Data)\n",
    "print \"Data Length:\", len(Data)\n",
    "if len(Data)>0: \n",
    "    print \"Type of Data Contents:\", type(Data[0])\n",
    "    print \"Data Minimum:\", min(Data)\n",
    "    print \"Data Maximum:\", max(Data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise A.1.2\n",
    "\n",
    "Write a function that computes the mean of values in a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Skeleton\n",
    "def mean(Data):\n",
    "    m=0\n",
    "    ### BEGIN SOLUTION\n",
    "\n",
    "    # Fill in your solution here        \n",
    "    \n",
    "    ### END SOLUTION   \n",
    "    return m\n",
    "\n",
    "print \"Mean of Data:\", mean(Data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise A.1.3\n",
    "\n",
    "Write a function the applies a booling function (that returns true/false) to every element in data, and returns a list of indices of elements where the result was true. Use this function to find the indices of positive entries. (This might be something you want to do if you are applying some selection criteria to your dataset and only want to keep events/entries/examples that pass the criteria.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def where(mylist,myfunc):\n",
    "    out= []\n",
    "    ### BEGIN SOLUTION\n",
    "\n",
    "    # Fill in your solution here        \n",
    "    \n",
    "    ### END SOLUTION    \n",
    "    return out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise A.1.4\n",
    "\n",
    "The inrange(mymin,mymax) function below returns a function that tests if it's input is between the specified values. Use this function, in conjunction to your solution to Exercise 1.3, to demonstrate that your data is \"flat\". Hint: pick several sub-ranges and show that the number of data point divided by the size of the range is roughly constant. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inrange(mymin,mymax):\n",
    "    def testrange(x):\n",
    "        return x<mymax and x>=mymin\n",
    "    return testrange\n",
    "\n",
    "# Examples:\n",
    "F1=inrange(0,10)\n",
    "F2=inrange(10,20)\n",
    "\n",
    "print F1(0), F1(1), F1(10), F1(15), F1(20)\n",
    "print F2(0), F2(1), F2(10), F2(15), F2(20)\n",
    "\n",
    "print \"Number of Entries passing F1:\", len(where(Data,F1))\n",
    "print \"Number of Entries passing F2:\", len(where(Data,F2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise A.1.5\n",
    "\n",
    "Repeat Exercise 1.4 using the built in python functions sum and map instead of your solution to 1.3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### BEGIN SOLUTION\n",
    "\n",
    "# Fill in your solution here        \n",
    "    \n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise A.1.6\n",
    "\n",
    "Write a new function called GenerateDataFromFunction(N,mymin,mymax,myfunc), that instead of generating a flat distribution, generates a distribution with a functional form coded in myfunc. Note that for this exercise myfunc should always be > 0.  \n",
    "\n",
    "For this exercise, let us make myfunc a Gaussian distribution, which is given below. Generate 1000 numbers that follow this Gaussian distribution. Confirm that the mean of the generated data is close to mean you specified when building the Gaussian. \n",
    "\n",
    "Hint: A simple, but slow, solution to generate data with a given distribution is to a draw random number, let's call it test_x within the specified range (mymin, mymax) and another number p between the minimum and maximum values of the function myfunc (which you will have to determine). If p<=function(test_x), then place test_x on the output. If not, repeat the process, drawing two new numbers. Repeat until you have the specified number of generated numbers, N, in the output. This method is often called \"the accept and reject sampling method\". For this problem, it's OK to determine the min and max by numerically sampling the function.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def GenerateDataFromFunction(N,mymin,mymax,myfunc):\n",
    "    out = []\n",
    "    ### BEGIN SOLUTION\n",
    "\n",
    "    # Fill in your solution here        \n",
    "    \n",
    "    ### END SOLUTION   \n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def gaussian(mean, sigma):\n",
    "    def f(x):\n",
    "        return (1/math.sqrt(2*math.pi*sigma**2))*math.exp(-( (x-mean)**2)/(2*(sigma**2)   ))\n",
    "    return f\n",
    "\n",
    "# Example Instantiation\n",
    "g1=gaussian(0,1)\n",
    "g2=gaussian(10,3)\n",
    "\n",
    "### BEGIN SOLUTION\n",
    "\n",
    "# Fill in your solution here        \n",
    "    \n",
    "### END SOLUTION\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Numpy'></a>\n",
    "## 2. Numpy\n",
    "\n",
    "[Numpy](http://www.numpy.org) is the tensor manipulation package most commonly used in python-based scientific computing. Numpy tensor interface is also adopted by all packages that provide tensors (e.g. h5py, theano, TensorFlow, ...). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise A.2.1\n",
    "\n",
    "Let start with some basic reshape manipulations. Consider a classification task. We can imagine the training data X consisting of N examples each with M inputs, so the shape of X is (M,N). The output of the Neural Network for the training sample encodes the true class of each of the N examples in X, in a \"one-hot\" matrix of shape (N,C), where C is the number of classes and each row corresponds to the true class for the corresponding example in X. So for a given row Y[i], all elements are 0 except for the column corresponding to the true class.\n",
    "\n",
    "For example consider a classification task of separating between 4 classes. We'll call them A, B, C, and D.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "Y=np.array( [ [0, 1, 0, 0], # Class B\n",
    "              [1, 0, 0, 0], # Class A\n",
    "              [0, 0, 0, 1], # Class C\n",
    "              [0, 0, 1, 0]  # Class D\n",
    "            ])\n",
    "\n",
    "print \"Shape of Y:\", Y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets imagine that we want to change to 2 classes instead by combining classes A with B and C with D. Use np.reshape and np.sum to create a new vector Y1. Hint: change the shape of Y into (8,2), sum along the correct axes, and change shape to (4,2). LH: solution given"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print \"Transpose:\", np.transpose(Y)\n",
    "print \"Reshape 8,2:\", np.transpose(Y).reshape((8,2))\n",
    "print \"Sum:\", np.sum(np.transpose(Y).reshape((8,2)),axis=1)\n",
    "\n",
    "Y1= np.sum(np.transpose(Y)\n",
    "           .reshape((8,2)),axis=1).reshape(4,2)\n",
    "print \"Answer: \",Y1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise A.2.2\n",
    "\n",
    "Oftentimes we find that neutral networks work best when their input is mostly between 0,1. Below, we create a random dataset that is normal distributed (mean of 4, sigma of 10). Shift the data so that the mean is 0.5 and 68% of the data lies between 0 and 1. LH: solution given."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X=np.random.normal(4,10,1000)\n",
    "print np.mean(X)\n",
    "print np.min(X)\n",
    "print np.max(X)\n",
    "print np.var(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "X1=(X-np.mean(X))/math.sqrt(np.var(X)) # Replace X with your answer\n",
    "\n",
    "print np.mean(X1)\n",
    "print np.var(X1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise A.2.3\n",
    "\n",
    "Using np.random.random and np.random.normal to generate two datasets. Then use np.where to repeat exercise 1.4 showing that one creates a flat distribution and the other does not by binning the data by hand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X0=np.random.random(1000)\n",
    "\n",
    "def CheckFlatness(D,steps=10):\n",
    "    maxD=np.max(D)\n",
    "    minD=np.min(D)\n",
    "    i=minD\n",
    "    stepsize=(maxD-minD)/steps\n",
    "    while i<maxD:\n",
    "        print i,i+stepsize,\":\",np.shape(np.where((D<=(i+stepsize)) & (D>i) ))\n",
    "        i+=stepsize\n",
    "        \n",
    "CheckFlatness(X0)\n",
    "CheckFlatness(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='h5py'></a>\n",
    "## 3. h5py\n",
    "\n",
    "[HDF5](https://support.hdfgroup.org/HDF5/) is a \"data model, library, and file format for storing and managing data.\" It is also the most common storage format in data science. [h5py](http://www.h5py.org) provides a python API for HDF5. In most cases, you do not need to know very much about HDF5 or h5py, just how to read/write tensors into/from files, which you can easily pick up from the [h5py Quick Start](http://docs.h5py.org/en/latest/quick.html#quick). We won't be using HDF5 for this tutorial. This section is here as reference for when you do encounter an HDF5 file.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='MachineLearning'></a>\n",
    "# B. Machine Learning\n",
    "\n",
    "For the remainder of this tutorial, we will attempt to follow the first paper on Deep Learning in High Energy physics [P. Baldi, et al](https://arxiv.org/pdf/1402.4735.pdf). This paper demonstrates that Deep Neural Networks can learn from raw data the features that are typically used for searches for exotics particles. The authors publically provide the two benchmark scenarios considered in the paper. We will focus on the SUSY benchmark. \n",
    "\n",
    "<a id='Dataset'></a>\n",
    "## 1. The Dataset\n",
    "\n",
    "The data is distributed as a comma separated values (CSV) file. If you are running on lxplus, you can find a local copy to use as input. Otherwise, download the ~ GB compressed file from [UCI's ML Archive](http://archive.ics.uci.edu/ml/datasets/SUSY), use `gunzip` to decompress it, and change the path to the file below accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note filenames are defined in cell at top of this tutorial\n",
    "\n",
    "# print out the first 5 lines using unix head command (note in jupyter ! => shell command)\n",
    "!head -5 $SUSY_filename"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each row represents a LHC collision event. Each column contains some observable from that event. The variable names are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VarNames=[\"signal\", \"l_1_pT\", \"l_1_eta\",\"l_1_phi\", \"l_2_pT\", \"l_2_eta\", \"l_2_phi\", \"MET\", \"MET_phi\", \"MET_rel\", \"axial_MET\", \"M_R\", \"M_TR_2\", \"R\", \"MT2\", \"S_R\", \"M_Delta_R\", \"dPhi_r_b\", \"cos_theta_r1\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of these variables represent the \"raw\" kinematics of the observed final state particles, while others are \"features\" that are derived from these raw quantities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RawNames=[\"l_1_pT\", \"l_1_eta\",\"l_1_phi\", \"l_2_pT\", \"l_2_eta\", \"l_2_phi\"]\n",
    "FeatureNames=[ \"MET\", \"MET_phi\", \"MET_rel\", \"axial_MET\", \"M_R\", \"M_TR_2\", \"R\", \"MT2\", \"S_R\", \"M_Delta_R\", \"dPhi_r_b\", \"cos_theta_r1\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Pandas'></a>\n",
    "## 2. Pandas\n",
    "\n",
    "We will use [pandas](http://pandas.pydata.org) to read in the file, and [matplotlib](https://matplotlib.org) to make plots. Pandas provides \"data structures and data analysis tools for the Python Programming Language\". Many machine learning tasks can be accomplished with [numpy](http://www.numpy.org) tensors and [h5py](http://www.h5py.org) files. In this case, pandas just makes it very easy to read a CSV file.\n",
    "\n",
    "The following ensures pandas is installed and sets everything up:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can read the data into a pandas dataframe. It's a ~GB file, so be patient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(SUSY_filename, dtype='float64', names=VarNames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another nice feature of pandas is that you can see the data in Jupyter by just evaluating the dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first column stores the \"truth\" label of whether an event was signal or background. Pandas makes it easy to create dataframes that store only the signal or background events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_sig=df[df.signal==1]\n",
    "df_bkg=df[df.signal==0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following example plots the signal and background distributions of every variable. Note that we use VarNames[1:] to skip the first variable, which was the true label. \n",
    "\n",
    "We will use matplotlib for plotting. There are lots of tutorials and primers out there that you can find searching the web. A good tutorial can be found in the [Scipy Lectures](http://www.scipy-lectures.org/intro/matplotlib/matplotlib.html). Look through these on your own time, it is not necessary for doing these exercise. The code below is all you need to know for making histograms with matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for var in VarNames[1:]:\n",
    "    print var\n",
    "    plt.figure()\n",
    "    plt.hist(np.array(df_sig[var]),bins=100,histtype=\"step\", color=\"red\",label=\"background\",stacked=True)\n",
    "    plt.hist(np.array(df_bkg[var]),bins=100,histtype=\"step\", color=\"blue\", label=\"signal\",stacked=True)\n",
    "    plt.legend(loc='upper right')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Scikit-learn'></a>\n",
    "## 3. Scikit-learn\n",
    "\n",
    "[Scikit-learn](http://scikit-learn.org) is a rich python library for data science, including machine learning. As an example, we can easily build a Fisher Discriminant (aka Linear Discriminant Analysis, or LDA). The [LDA Documentation](http://scikit-learn.org/stable/modules/lda_qda.html#dimensionality-reduction-using-linear-discriminant-analysis) does as great job explaining this classifier. Here's how we instanciate the classifier: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn.discriminant_analysis as DA\n",
    "Fisher=DA.LinearDiscriminantAnalysis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets separate the data into inputs (X) vs outputs (Y) and training vs testing samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_Train=4000000\n",
    "\n",
    "Train_Sample=df[:N_Train]\n",
    "Test_Sample=df[N_Train:]\n",
    "\n",
    "X_Train=Train_Sample[VarNames[1:]]\n",
    "y_Train=Train_Sample[\"signal\"]\n",
    "\n",
    "X_Test=Test_Sample[VarNames[1:]]\n",
    "y_Test=Test_Sample[\"signal\"]\n",
    "\n",
    "Test_sig=Test_Sample[Test_Sample.signal==1]\n",
    "Test_bkg=Test_Sample[Test_Sample.signal==0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can train the classifier as follow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fisher.fit(X_Train,y_Train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the output, comparing signal and background:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(Fisher.decision_function(Test_sig[VarNames[1:]]),bins=100,histtype=\"step\", color=\"blue\", label=\"signal\",stacked=True)\n",
    "plt.hist(Fisher.decision_function(Test_bkg[VarNames[1:]]),bins=100,histtype=\"step\", color=\"red\", label=\"background\",stacked=True)\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can make a ROC curve and evaluate the AUC:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_curve, auc\n",
    "fpr, tpr, _ = roc_curve(y_Test, Fisher.decision_function(X_Test))\n",
    "\n",
    "roc_auc = auc(fpr, tpr)\n",
    "\n",
    "plt.plot(fpr,tpr,color='darkorange',label='ROC curve (area = %0.2f)' % roc_auc)\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise B.3.1\n",
    "\n",
    "Train the Fisher performance using the raw, features, and raw+features as input. Compare the performance one a single plot. Add cells to this notebook as needed. Or start new notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_Train_Raw=Train_Sample[RawNames]\n",
    "X_Test_Raw=Test_Sample[RawNames]\n",
    "\n",
    "X_Train_Features=Train_Sample[FeatureNames]\n",
    "X_Test_Features=Test_Sample[FeatureNames]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TrainFisher(X_Train,X_Test,y_Train):\n",
    "    Fisher=DA.LinearDiscriminantAnalysis()\n",
    "    Fisher.fit(X_Train,y_Train)\n",
    "\n",
    "    fpr, tpr, _ = roc_curve(y_Test, Fisher.decision_function(X_Test))\n",
    "    roc_auc = auc(fpr, tpr)\n",
    "\n",
    "    plt.plot(fpr,tpr,color='darkorange',label='ROC curve (area = %0.2f)' % roc_auc)\n",
    "    plt.legend(loc=\"lower right\")\n",
    "    plt.xlabel('False Positive Rate')\n",
    "    plt.ylabel('True Positive Rate')\n",
    "    plt.show()\n",
    "    \n",
    "    return Fisher\n",
    "\n",
    "RawFisher=TrainFisher(X_Train_Raw,X_Test_Raw,y_Train)\n",
    "FeatureFisher=TrainFisher(X_Train_Features,X_Test_Features,y_Train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise B.3.2\n",
    "\n",
    "Select 3 different classifiers from the techniques listed [here](http://scikit-learn.org/stable/supervised_learning.html#supervised-learning). Note that you can use the multi-layer peceptron to build a deep network, though training may be prohibitively slow, so avoid this technique. Perform the comparison in exercise 1 for each classifier. Compare your conclusions for your selected techniques to the paper.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise B.3.3\n",
    "\n",
    "The following function calculates the significance of the observation of the signal given the number of expected Signal and Background events, using the simple formula $\\sigma_S= \\frac{N_S}{\\sqrt{N_S+N_B}}$. Read through the code carefully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PlotSignificance(N_S,N_B, N_S_min=1):\n",
    "    plt.figure()\n",
    "    eff_sig,bins_sig,p_sig=plt.hist(Fisher.decision_function(Test_sig[VarNames[1:]]),bins=100,histtype=\"step\", color=\"blue\", label=\"signal\",cumulative=-1,stacked=True,normed=True)\n",
    "    eff_bkg,bins_bkg,p_bkg=plt.hist(Fisher.decision_function(Test_bkg[VarNames[1:]]),bins=100,histtype=\"step\", color=\"red\", label=\"background\",cumulative=-1,stacked=True,normed=True)\n",
    "    plt.legend(loc='upper right')\n",
    "    plt.show()\n",
    "    \n",
    "    good_bins = np.where(eff_sig*N_S>=N_S_min)\n",
    "\n",
    "    print len(good_bins[0])\n",
    "    if len(good_bins[0])<1:\n",
    "        print \"Insufficient Signal.\"\n",
    "        return 0,0,0\n",
    "    \n",
    "    significance=(N_S*eff_sig)/np.sqrt((N_B*eff_bkg)+(N_S*eff_sig))\n",
    "\n",
    "    plt.figure()\n",
    "    plt.plot(bins_sig[:-1],significance)\n",
    "    \n",
    "    max_sign=np.max(significance[good_bins])\n",
    "    max_signI=np.argmax(significance[good_bins])\n",
    "    \n",
    "    plt.show()\n",
    "    print \"Max significance at \", bins_sig[max_signI], \" of\", max_sign\n",
    "    return bins_sig[max_signI],max_sign, max_signI\n",
    "    \n",
    "PlotSignificance(1000000,1e11)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer the following questions:\n",
    "   * What are we computing when making a normalized cummulative plot? \n",
    "   * Assume that the experiment produces 1 signal event for every $10^{11}$ background events. For each of your classifiers, how many signal events need to be produced to be able to make a $5\\sigma$ discovery claim?\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise B.3.4\n",
    "\n",
    "Read the Baldi, et al. paper and attempt to reproduce the results, as closely as possible, using scikit-learn. \n",
    "Try using the [multi-layer peceptron](http://scikitlearn.org/stable/modules/neural_networks_supervised.html#multi-layer-perceptron) to build a deep network. Or if you are capable, try it using [Keras Scikit-learn interface](https://keras.io/scikit-learn-api/). \n",
    "                                                "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='DeepLearning'></a>\n",
    "# C. Deep Learning\n",
    "\n",
    "This section is meant to get you started in using Keras to design Deep Neural Networks. The goal here is to simply repeat section B with Deep Learning.\n",
    "\n",
    "If you are starting here and have not run the cells above that load the data, you will need to run the following cell: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "VarNames=[\"signal\", \"l_1_pT\", \"l_1_eta\",\"l_1_phi\", \"l_2_pT\", \"l_2_eta\", \"l_2_phi\", \"MET\", \"MET_phi\", \"MET_rel\", \"axial_MET\", \"M_R\", \"M_TR_2\", \"R\", \"MT2\", \"S_R\", \"M_Delta_R\", \"dPhi_r_b\", \"cos_theta_r1\"]\n",
    "RawNames=[\"l_1_pT\", \"l_1_eta\",\"l_1_phi\", \"l_2_pT\", \"l_2_eta\", \"l_2_phi\"]\n",
    "FeatureNames=[ \"MET\", \"MET_phi\", \"MET_rel\", \"axial_MET\", \"M_R\", \"M_TR_2\", \"R\", \"MT2\", \"S_R\", \"M_Delta_R\", \"dPhi_r_b\", \"cos_theta_r1\"]\n",
    "\n",
    "df = pd.read_csv(SUSY_filename, dtype='float64', names=VarNames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets define training and test samples. Note that DNNs take very long to train, so for testing purposes we will use only about 10% of the 5 million events in the training/validation sample. Once you get everything working, you can go back and make the final version of your plots with the full sample. \n",
    "\n",
    "Also note that Keras had trouble with the Pandas tensors, so after doing all of the nice manipulation that Pandas enables, we convert the Tensor to a regular numpy tensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_Max=550000\n",
    "N_Train=500000\n",
    "\n",
    "Train_Sample=df[:N_Train]\n",
    "Test_Sample=df[N_Train:N_Max]\n",
    "\n",
    "X_Train=np.array(Train_Sample[VarNames[1:]])\n",
    "y_Train=np.array(Train_Sample[\"signal\"])\n",
    "\n",
    "X_Test=np.array(Test_Sample[VarNames[1:]])\n",
    "y_Test=np.array(Test_Sample[\"signal\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Keras'></a>\n",
    "## 1. Keras\n",
    "\n",
    "Training Deep Learning models can take a very long time. If you have access to a GPU, training with the GPU will be about 2 orders of magnitude faster that training with just the CPU. Unforunately, there are no GPUs on lxplus. But, if you are running this notebook on a system with NVidia GPU(s) properly setup, you can tell Keras to use a specific GPU:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Since lxplus does not have any GPUs, please DO NOT RUN THIS CELL ON LXPLUS.\n",
    "# Selecting First GPU in the system\n",
    "\n",
    "import os\n",
    "\n",
    "# Selecting GPU manually:\n",
    "gpuid= 0\n",
    "\n",
    "# If running on UTA DL Cluster and using GPU:\n",
    "# print \"Using Queue:\", os.environ[\"PBS_QUEUE\"]\n",
    "# gpuid=int(os.environ[\"PBS_QUEUE\"][3:4])\n",
    "\n",
    "print \"Selected GPU:\", gpuid\n",
    "\n",
    "# If Using TensorFlow:\n",
    "import tensorflow as tf\n",
    "# Comment if using GPU\n",
    "os.environ['CUDA_VISIBLE_DEVICES'] = '-1'\n",
    "config = tf.ConfigProto(intra_op_parallelism_threads=0,\n",
    "                        inter_op_parallelism_threads=0,\n",
    "                        allow_soft_placement=True,\n",
    "                        # Comment if using GPU\n",
    "                        device_count={'CPU' : 1, 'GPU' : 0},                                                \n",
    "                        # Uncomment to use GPU\n",
    "                        # gpu_options=tf.GPUOptions(visible_device_list=\\\"{}\\\".format(gpu_id),\n",
    "                        #                          force_gpu_compatible=True,\n",
    "                        #                          allow_growth=True)\n",
    "                       )\n",
    "\n",
    "# If using Theano and GPU:\n",
    "# os.environ['THEANO_FLAGS'] = \"\" # \"mode=FAST_RUN,device=gpu%s,floatX=float32,force_device=True\" % (gpuid)\n",
    "# import theano\n",
    "# theano.config.profile=False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will build a simple model. Note that this is a very small model, so things run fast. You should attempt more ambitious models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.models import Sequential\n",
    "from keras.layers import Dense\n",
    "\n",
    "model = Sequential()\n",
    "model.add(Dense(12, input_dim=X_Train.shape[1], init='uniform', activation='relu'))\n",
    "model.add(Dense(8, init='uniform', activation='relu'))\n",
    "model.add(Dense(1, init='uniform', activation='sigmoid'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model has to be compiled. At this time we set the loss function and the optimizer too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we train. We are running only 10 epochs in this example. Models may need hundreds of epochs before they stop improving."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "history=model.fit(X_Train, y_Train, validation_data=(X_Test,y_Test), nb_epoch=10, batch_size=2048)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model history keeps track of the loss and accuracy for each epoch. Note that the training above was setup to run on the validation sample at the end of each epoch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print history.history"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can plot the loss versus epoch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_history=history.history[\"loss\"]\n",
    "plt.plot(range(len(loss_history)),loss_history)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise C.1.1\n",
    "\n",
    "You will need to create several models and make sure they are properly trained. Write a function that takes this history and plots the values versus epoch. For every model that you train in the remainder of this lab, assess:\n",
    "\n",
    "    * Has you model's performance plateaued? If not train for more epochs. \n",
    "\n",
    "    * Compare the performance on training versus test sample. Are you over training?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Your Solution Here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise C.1.2\n",
    "\n",
    "Explicitly overtrain your model by running for a large number of epochs. Demonstrate that you are over-fitting. Then add dropout layers to your model to alleviate the problem.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise C.1.3\n",
    "\n",
    "Following section B, make a comparison of the performance between models trained with raw, features, and raw+features data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can evaluate how the trained model does on the test sample as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scores = model.evaluate(X_Test, y_Test)\n",
    "print scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can make ROC curves as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_curve, auc\n",
    "fpr, tpr, _ = roc_curve(y_Test, model.predict(X_Test))\n",
    "                        \n",
    "roc_auc = auc(fpr, tpr)\n",
    "\n",
    "plt.plot(fpr,tpr,color='darkorange',label='ROC curve (area = %0.2f)' % roc_auc)\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise C.1.4\n",
    "\n",
    "Again, following section B, design and implement at least 3 different DNN models. Train them and compare performance. You may try different architectures, loss functions, and optimizers to see if there is an effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Your solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise C.1.5\n",
    "\n",
    "Write a function that evaluates the performance (AUC) as a function of a given input variable. You will need to bin the test data in the variable (i.e. make sub-samples for events which have the particular variable in a range), evaluate the performance in each bin, and plot the results.\n",
    "\n",
    "Apply your function to each input variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Your solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HiggsML Dataset\n",
    "\n",
    "The HiggsML challenge was organized by ATLAS colleagues and ran on Kaggle's platform. The challenge provided a training and test set of Higgs and background events in CSV format. Using this data, participants were tasked with creating a classifier, which they submitted to Kaggle. Kaggle evaluated the classifier against another test set. At the end of the competition, the best performing classifiers were awarded a cash prize.\n",
    "\n",
    "The challenge is described in https://higgsml.lal.in2p3.fr\n",
    "\n",
    "The Kaggle site is https://www.kaggle.com/c/higgs-boson\n",
    "\n",
    "Detail description of the data and challenge: https://higgsml.lal.in2p3.fr/files/2014/04/documentation_v1.8.pdf\n",
    "\n",
    "You are welcome to use Scikit or any other package you like.\n",
    "\n",
    "Please separate different steps into different Jupyter Notebooks. For example:\n",
    "\n",
    "A copy of the data CSV files are on the cluster at: /data/afarbin/DLClass/HiggsML\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}