exercise_1-Submission.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import _pickle as cPickle\n",
    "import os\n",
    "import gzip\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First exercise: Classifying MNIST with MLPs\n",
    "In this exercise you will implement a Neural Network (or MLP) and classify the MNIST digits with it.\n",
    "MNIST is a \"well hung\" dataset that has been used a lot over the years to benchmark different classification algorithms. \n",
    "To learn more about it have a look here: http://yann.lecun.com/exdb/mnist/ ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Loading\n",
    "We first define a function for downloading and loading MNIST.\n",
    "**WARNING**: Executing it will obviously use up some space on your machine ;). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mnist(datasets_dir='./data'):\n",
    "    if not os.path.exists(datasets_dir):\n",
    "        os.mkdir(datasets_dir)\n",
    "    data_file = os.path.join(datasets_dir, 'mnist.pkl.gz')\n",
    "    if not os.path.exists(data_file):\n",
    "        print('... downloading MNIST from the web')\n",
    "        try:\n",
    "            import urllib\n",
    "            urllib.urlretrieve('http://google.com')\n",
    "        except AttributeError:\n",
    "            import urllib.request as urllib\n",
    "        url = 'http://www.iro.umontreal.ca/~lisa/deep/data/mnist/mnist.pkl.gz'\n",
    "        urllib.urlretrieve(url, data_file)\n",
    "\n",
    "    print('... loading data')\n",
    "    # Load the dataset\n",
    "    f = gzip.open(data_file, 'rb')\n",
    "    try:\n",
    "        train_set, valid_set, test_set = cPickle.load(f, encoding=\"latin1\")\n",
    "    except TypeError:\n",
    "        train_set, valid_set, test_set = cPickle.load(f)\n",
    "    f.close()\n",
    "\n",
    "    test_x, test_y = test_set\n",
    "    test_x = test_x.astype('float32')\n",
    "    test_x = test_x.astype('float32').reshape(test_x.shape[0], 1, 28, 28)\n",
    "    test_y = test_y.astype('int32')\n",
    "    valid_x, valid_y = valid_set\n",
    "    valid_x = valid_x.astype('float32')\n",
    "    valid_x = valid_x.astype('float32').reshape(valid_x.shape[0], 1, 28, 28)\n",
    "    valid_y = valid_y.astype('int32')\n",
    "    train_x, train_y = train_set\n",
    "    train_x = train_x.astype('float32').reshape(train_x.shape[0], 1, 28, 28)\n",
    "    train_y = train_y.astype('int32')\n",
    "    rval = [(train_x, train_y), (valid_x, valid_y), (test_x, test_y)]\n",
    "    print('... done loading data')\n",
    "    return rval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural Network Layers\n",
    "We now define \"bare bone\" neural network layers.\n",
    "The parts marked with **TODO** are where you should finish the implementation!\n",
    "Conceptually we will implement the layers as follows:\n",
    "\n",
    "Each layer has a constructor that takes an input layer plus some additional arguments such as layer size and the activation function name. The layer then uses the provided input layer to compute the layer dimensions, weight shapes, etc. and setup all auxilliary variables.\n",
    "\n",
    "Each layer then has to provide three functions (as defined in the Layer class below): *output_shape()*, *fprop()* and *brop()*. The output_shape function is used to figure out the shape for the next layer and the *fprop()/bprop()* functions are used to compute forward and backward passes through the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# start by defining simple helpers\n",
    "def sigmoid(x):\n",
    "    return 1.0/(1.0+np.exp(-x))\n",
    "\n",
    "def sigmoid_d(x):\n",
    "    # TODO ################\n",
    "    # TODO: implement this!\n",
    "    # TODO ################\n",
    "    inter=sigmoid(x)\n",
    "    return (inter*(1.0-inter))\n",
    "    #return x\n",
    "\n",
    "def tanh(x):\n",
    "    return np.tanh(x)\n",
    "\n",
    "def tanh_d(x):\n",
    "    # TODO ################\n",
    "    # TODO: implement this!\n",
    "    # TODO ################\n",
    "    inter=np.tanh(x)\n",
    "    return( 1.0 - np.square(inter))\n",
    "\n",
    "\n",
    "def relu(x):\n",
    "    #print(\"input--main\",x)\n",
    "    #dx=np.zeros_like(x)\n",
    "    #dx[x > 0] = x\n",
    "    return np.maximum(x,0.0)\n",
    "\n",
    "def relu_d(x):\n",
    "    # TODO ################\n",
    "    # TODO: implement this!\n",
    "    # TODO ################\n",
    "    dx=np.zeros_like(x)\n",
    "    dx[x >= 0] = 1. \n",
    "    return dx\n",
    "\n",
    "def softmax(x, axis=1):\n",
    "    # to make the softmax a \"safe\" operation we will \n",
    "    # first subtract the maximum along the specified axis\n",
    "    # so that np.exp(x) does not blow up!\n",
    "    # Note that this does not change the output.\n",
    "    x_max = np.max(x, axis=axis, keepdims=True)\n",
    "    x_safe = x - x_max\n",
    "    e_x = np.exp(x_safe)\n",
    "    return e_x / np.sum(e_x, axis=axis, keepdims=True)\n",
    "\n",
    "def one_hot(labels):\n",
    "    \"\"\"this creates a one hot encoding from a flat vector:\n",
    "    i.e. given y = [0,2,1]\n",
    "     it creates y_one_hot = [[1,0,0], [0,0,1], [0,1,0]]\n",
    "    \"\"\"\n",
    "    classes = np.unique(labels)\n",
    "    n_classes = classes.size\n",
    "    one_hot_labels = np.zeros(labels.shape + (n_classes,))\n",
    "    for c in classes:\n",
    "        one_hot_labels[labels == c, c] = 1\n",
    "    return one_hot_labels\n",
    "\n",
    "def unhot(one_hot_labels):\n",
    "    \"\"\" Invert a one hot encoding, creating a flat vector \"\"\"\n",
    "    return np.argmax(one_hot_labels, axis=-1)\n",
    "\n",
    "# then define an activation function class\n",
    "class Activation(object):\n",
    "    \n",
    "    def __init__(self, tname):\n",
    "        if tname == 'sigmoid':\n",
    "            self.act = sigmoid\n",
    "            self.act_d = sigmoid_d\n",
    "        elif tname == 'tanh':\n",
    "            self.act = tanh\n",
    "            self.act_d = tanh_d\n",
    "        elif tname == 'relu':\n",
    "            self.act = relu\n",
    "            self.act_d = relu_d\n",
    "        else:\n",
    "            raise ValueError('Invalid activation function.')\n",
    "            \n",
    "    def fprop(self, input):\n",
    "        # we need to remember the last input\n",
    "        # so that we can calculate the derivative with respect\n",
    "        # to it later on\n",
    "        self.last_input = input #z0\n",
    "        return self.act(input) #h0=g0(z0)\n",
    "    \n",
    "    def bprop(self, output_grad):\n",
    "        #return np.dot (self.act_d(self.last_input).T, output_grad ) #dl/dz\n",
    "        return  output_grad * self.act_d(self.last_input)  #dl/dz # hadamard product\n",
    "\n",
    "# define a base class for layers\n",
    "class Layer(object):\n",
    "    \n",
    "    def fprop(self, input):\n",
    "        \"\"\" Calculate layer output for given input \n",
    "            (forward propagation). \n",
    "        \"\"\"        \n",
    "        raise NotImplementedError('This is an interface class, please use a derived instance')\n",
    "\n",
    "    def bprop(self, output_grad):\n",
    "        \"\"\" Calculate input gradient and gradient \n",
    "            with respect to weights and bias (backpropagation). \n",
    "        \"\"\"\n",
    "        raise NotImplementedError('This is an interface class, please use a derived instance')\n",
    "\n",
    "    def output_size(self):\n",
    "        \"\"\" Calculate size of this layer's output.\n",
    "        input_shape[0] is the number of samples in the input.\n",
    "        input_shape[1:] is the shape of the feature.\n",
    "        \"\"\"\n",
    "        raise NotImplementedError('This is an interface class, please use a derived instance')\n",
    "\n",
    "# define a base class for loss outputs\n",
    "# an output layer can then simply be derived\n",
    "# from both Layer and Loss \n",
    "class Loss(object):\n",
    "\n",
    "    def loss(self, output, output_net):\n",
    "        \"\"\" Calculate mean loss given real output and network output. \"\"\"\n",
    "        raise NotImplementedError('This is an interface class, please use a derived instance')\n",
    "\n",
    "    def input_grad(self, output, output_net):\n",
    "        \"\"\" Calculate input gradient real output and network output. \"\"\"\n",
    "        raise NotImplementedError('This is an interface class, please use a derived instance')\n",
    "\n",
    "# define a base class for parameterized things        \n",
    "class Parameterized(object):\n",
    "    \n",
    "    def params(self):\n",
    "        \"\"\" Return parameters (by reference) \"\"\"\n",
    "        raise NotImplementedError('This is an interface class, please use a derived instance')\n",
    "    \n",
    "    def grad_params(self):\n",
    "        \"\"\" Return accumulated gradient with respect to params. \"\"\"\n",
    "        raise NotImplementedError('This is an interface class, please use a derived instance')\n",
    "\n",
    "# define a container for providing input to the network\n",
    "class InputLayer(Layer):\n",
    "    \n",
    "    def __init__(self, input_shape):\n",
    "        if not isinstance(input_shape, tuple):\n",
    "            raise ValueError(\"InputLayer requires input_shape as a tuple\")\n",
    "        self.input_shape = input_shape\n",
    "\n",
    "    def output_size(self):\n",
    "        return self.input_shape\n",
    "    \n",
    "    def fprop(self, input):\n",
    "        return input\n",
    "    \n",
    "    def bprop(self, output_grad):\n",
    "        return output_grad\n",
    "        \n",
    "class FullyConnectedLayer(Layer, Parameterized):\n",
    "    \"\"\" A standard fully connected hidden layer, as discussed in the lecture.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, input_layer, num_units, \n",
    "                 init_stddev, activation_fun=Activation('relu')):\n",
    "        self.num_units = num_units\n",
    "        self.activation_fun = activation_fun\n",
    "        # the input shape will be of size (batch_size, num_units_prev) \n",
    "        # where num_units_prev is the number of units in the input \n",
    "        # (previous) layer\n",
    "        self.input_shape = input_layer.output_size()\n",
    "        \n",
    "        #print(\"printing details of each layer\\n\")\n",
    "        #print(\"input shape\",self.input_shape)\n",
    "        #print(\"number of units\", self.num_units)\n",
    "        \n",
    "        \n",
    "       # self.last_input = input\n",
    "        # TODO ################################\n",
    "        # TODO: implement weight initialization\n",
    "        # TODO ################################\n",
    "        # this is the weight matrix it should have shape: (num_units_prev, num_units)\n",
    "        #print (self.input_shape[1])\n",
    "        #print (self.num_units)\n",
    "        self.W = np.random.randn(self.input_shape[1],self.num_units) * init_stddev #FIXME\n",
    "        #self.W = self.W /float(self.W.shape[0])\n",
    "        #print(\"shape of weight \",self.W.shape)\n",
    "        # and this is the bias vector of shape: (num_units)\n",
    "\n",
    "        self.b = np.random.randn(self.num_units)* init_stddev  #FIXME\n",
    "        #self.b= self.b/float(self.b.shape[0])\n",
    "        #print(\"shape of bias \",self.b.shape)\n",
    "        # create dummy variables for parameter gradients\n",
    "        # no need to change these here!\n",
    "        self.dW = None\n",
    "        self.db = None\n",
    "    \n",
    "    def output_size(self):\n",
    "        return (self.input_shape[0], self.num_units)\n",
    "    \n",
    "    def fprop(self, input):\n",
    "        # TODO ################################################\n",
    "        # TODO: implement forward propagation\n",
    "        # NOTE: you should also handle the case were \n",
    "        #       activation_fun is None (meaning no activation)\n",
    "        #       then this is simply a linear layer\n",
    "        # TODO ################################################\n",
    "        # you again want to cache the lasfpropt_input for the bprop\n",
    "        # implementation below!\n",
    "        self.last_input = input\n",
    "        \n",
    "        ##print(\"shape of X\",input.shape)\n",
    "        #print(\"shape of weight\",self.W.shape)\n",
    "        #print(\"shape of bias\",self.b.shape)\n",
    "        \n",
    "        \n",
    "        #print(\"printing input in fprop---------\", self.last_input)\n",
    "        \n",
    "        #print(\"shape of last input\",self.last_input.shape)\n",
    "        # NOTE: Use numpy dot product\n",
    "        # h = t(a), a = W * x + b \n",
    "        activation=np.dot(self.last_input,self.W)+self.b #a = W * x + b : in matrix notation - a= XW+b\n",
    "        if self.activation_fun is None :\n",
    "            hidden_layer=activation  # h = t(a)\n",
    "        else  :  \n",
    "            #a=Activation(activation_fun)\n",
    "            hidden_layer= self.activation_fun.fprop(activation)  # h = t(a)\n",
    "        return  hidden_layer\n",
    "        #FIXME\n",
    "        raise NotImplementedError(\"You should implement this!\")\n",
    "        return np.zeros(self.output_size)\n",
    "        \n",
    "    def bprop(self, output_grad):\n",
    "        \"\"\" Calculate input gradient (backpropagation). \"\"\"\n",
    "        #import pdb; pdb.set_trace()\n",
    "        # TODO ################################\n",
    "        # TODO: implement backward propagation\n",
    "        # TODO ###############################\n",
    "        \n",
    "        # HINT: you may have to divide dW and db by n\n",
    "        #       to make gradient checking work \n",
    "        #       OR you divide the gradient in the output layer by n\n",
    "        n = output_grad.shape[0]\n",
    "        # accumulate gradient wrt. the parameters first\n",
    "        # we will need to store these to later update\n",
    "        # the network after a few forward backward passes\n",
    "        # NOTE: you should also handle the case were \n",
    "        #       activation_fun is None (meaning no activation)\n",
    "        # the gradient wrt. W should be stored as self.dW\n",
    "        # the gradient wrt. b should be stored as selfdb\n",
    "       \n",
    "        # NOTE: self.dW is also a numpy dot product\n",
    "        \n",
    "        #self.dW = np.dot(np.dot(output_grad,self.act_d(self.last_input)),self.fprop(input))/n #FIXME\n",
    "        #self.db = np.dot(output_grad,self.act_d(self.last_input))/n #FIXME\n",
    "        #a=Activation(activation_fun)\n",
    "        \n",
    "        #self.dW = np.dot(self.activation_fun.fprop(self.activation_fun.last_input).T,\n",
    "        #                 self.activation_fun.bprop(output_grad)  )#h.T*dl/dz\n",
    "        \n",
    "        \n",
    "        #print(\"last input in bprop------\",self.last_input)\n",
    "        \n",
    "       # print(\"activation function in bprop--\",self.activation_fun)\n",
    "        \n",
    "        if self.activation_fun is None:\n",
    "            inter = output_grad\n",
    "        else:\n",
    "            inter = self.activation_fun.bprop(output_grad) # hadamard product, dl/dz\n",
    "        \n",
    "        grad_input = np.dot(inter,self.W.T)\n",
    "        \n",
    "        self.dW = np.dot(self.last_input.T,inter)\n",
    "        \n",
    "        #self.db = inter/n #dl/dz\n",
    "        self.db = np.dot(np.ones((1, n)), inter)\n",
    "        \n",
    "        \n",
    "        #print(\"dw:\",self.dW)\n",
    "        #print(\"db:\",self.db)\n",
    "        \n",
    "        # the gradient wrt. the input should be calculated here\n",
    "        \n",
    "        #print(\"printing shapes for b prop\")\n",
    "        #print(output_grad.shape)\n",
    "        #print(self.W.T.shape)\n",
    "        #print(self.activation_fun.act_d(self.activation_fun.last_input).shape)\n",
    "        \n",
    "        #maybe we have to consider the weights and activation gradients of prev layer\n",
    "        ##############check this!!!!!!!!!!!!!!!!!!!!!!!!!!!###########\n",
    "        \n",
    "        #grad_input=  np.dot(output_grad ,self.W.T ,\n",
    "        #self.activation_fun.act_d(self.activation_fun.last_input))\n",
    "        #grad_input = np.dot(self.activation_fun.act_d(self.activation_fun.last_input),self.W)\n",
    "        \n",
    "        #activation=np.dot(self.last_input,self.W)+self.b\n",
    "        \n",
    "        #grad_input=np.dot(self.activation_fun.act_d(activation),self.W.T)\n",
    "\n",
    "        return grad_input\n",
    "      \n",
    "        raise NotImplementedError(\"you should implement this\")\n",
    "        \n",
    "    def params(self):\n",
    "        return self.W, self.b\n",
    "\n",
    "    def grad_params(self):\n",
    "        \n",
    "        #print(\"inside grad_params\")\n",
    "        #print(\"dw:\",self.dW)\n",
    "        #print(\"db:\",self.db)\n",
    "        \n",
    "        return self.dW, self.db\n",
    "\n",
    "# finally we specify the interface for output layers \n",
    "# which are layers that also have a loss function\n",
    "# we will implement two output layers:\n",
    "#  a Linear, and Softmax (Logistic Regression) layer\n",
    "# The difference between output layers and and normal \n",
    "# layers is that they will be called to compute the gradient\n",
    "# of the loss through input_grad(). bprop will never \n",
    "# be called on them!\n",
    "class LinearOutput(Layer, Loss):\n",
    "    \"\"\" A simple linear output layer that  \n",
    "        uses a squared loss (e.g. should be used for regression)\n",
    "    \"\"\"\n",
    "    def __init__(self, input_layer):\n",
    "        self.input_size = input_layer.output_size()\n",
    "        \n",
    "    def output_size(self):\n",
    "        return (1,)\n",
    "        \n",
    "    def fprop(self, input):\n",
    "        return input\n",
    "\n",
    "    def bprop(self, output_grad):\n",
    "        raise NotImplementedError(\n",
    "            'LinearOutput should only be used as the last layer of a Network'\n",
    "            + ' bprop() should thus never be called on it!'\n",
    "        )\n",
    "    \n",
    "    def input_grad(self, Y, Y_pred):\n",
    "        # TODO ####################################\n",
    "        # TODO: implement gradient of squared loss\n",
    "        # TODO ####################################\n",
    "        #return np.zeros_like(Y_pred) #FIXME\n",
    "        n_samples = Y.shape[0]\n",
    "        return (Y_pred-Y)*1.0/n_samples\n",
    "\n",
    "    def loss(self, Y, Y_pred):\n",
    "        loss = 0.5 * np.square(Y_pred - Y)\n",
    "        return np.mean(np.sum(loss, axis=1))\n",
    "\n",
    "class SoftmaxOutput(Layer, Loss):\n",
    "    \"\"\" A softmax output layer that calculates \n",
    "        the negative log likelihood as loss\n",
    "        and should be used for classification.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, input_layer):\n",
    "        self.input_size = input_layer.output_size()\n",
    "        \n",
    "    def output_size(self):\n",
    "        return (1,)\n",
    "    \n",
    "    def fprop(self, input):\n",
    "        return softmax(input)\n",
    "    \n",
    "    def bprop(self, output_grad):\n",
    "        raise NotImplementedError(\n",
    "            'SoftmaxOutput should only be used as the last layer of a Network'\n",
    "            + ' bprop() should thus never be called on it!'\n",
    "        )\n",
    "    \n",
    "    def input_grad(self, Y, Y_pred):\n",
    "        # TODO #######################################################\n",
    "        # TODO: implement gradient of the negative log likelihood loss\n",
    "        # TODO #######################################################\n",
    "        # HINT: since this would involve taking the log \n",
    "        #       of the softmax (which is np.exp(x)/np.sum(x, axis=1))\n",
    "        #       this gradient computation can be simplified a lot,\n",
    "        #       you may find a connection to the LinearOutput layer!\n",
    "        # return np.zeros_like(Y_pred)\n",
    "        #return(-np.dot(Y,1/(Y_pred)) + np.dot((1-Y),1/(1-Y_pred)))\n",
    "        n_samples = Y.shape[0]\n",
    "        return (Y_pred-Y)*1.0/n_samples\n",
    "             \n",
    "\n",
    "    def loss(self, Y, Y_pred):\n",
    "        # Assume one-hot encoding of Y\n",
    "        out = Y_pred\n",
    "        # to make the loss numerically stable \n",
    "        # you should add an epsilon in the log ;)\n",
    "        eps = 1e-10\n",
    "        # TODO ####################################\n",
    "        # calculate negative log likelihood\n",
    "        # TODO ####################################\n",
    "        #loss= (Y)* np.log(out+eps) # + np.dot((1-Y),np.log(1-Y_pred+eps))\n",
    "        #loss = np.zeros_like(out) #FIXME\n",
    "        \n",
    "        \n",
    "        #print(\"loss\", -1*np.mean(np.sum(loss)))\n",
    "        loss = -np.sum(np.multiply(Y, np.log(out + eps)), axis=1)\n",
    "        #return -1*np.mean(np.sum(loss))\n",
    "        return np.mean(loss)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural Network class\n",
    "With all layers in place (and properly implemented by you) we can finally define a neural network.\n",
    "For our purposes a neural network is simply a collection of layers which we will cycle through and on which we will call fprop and bprop to compute partial derivatives with respect to the input and the parameters.\n",
    "\n",
    "Pay special attention to the *check_gradients()* function in which you should implement automatic differentiation. This function will become your best friend when checking the correctness of your implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "class NeuralNetwork:\n",
    "    \"\"\" Our Neural Network container class.\n",
    "    \"\"\"\n",
    "    def __init__(self, layers):\n",
    "        self.layers = layers\n",
    "        \n",
    "    def _loss(self, X, Y):\n",
    "        Y_pred = self.predict(X)\n",
    "        return self.layers[-1].loss(Y, Y_pred)\n",
    "\n",
    "    def predict(self, X):\n",
    "        \"\"\" Calculate an output Y for the given input X. \"\"\"\n",
    "        # TODO ##########################################\n",
    "        # TODO: implement forward pass through all layers\n",
    "        # TODO ##########################################\n",
    "        input_here = X\n",
    "        #loop through all layers in self.layers and calculate the complete chain\n",
    "        n = len(self.layers)\n",
    "        for i in range(1, n):\n",
    "            input_here = self.layers[i].fprop(input_here)\n",
    "\n",
    "        return input_here\n",
    "    \n",
    "    def backpropagate(self, Y, Y_pred, upto=0):\n",
    "        \"\"\" Backpropagation of partial derivatives through \n",
    "            the complete network up to layer 'upto'\n",
    "        \"\"\"\n",
    "        next_grad = self.layers[-1].input_grad(Y, Y_pred)\n",
    "        # TODO ##########################################\n",
    "        # TODO: implement backward pass through all layers\n",
    "        # TODO ##########################################\n",
    "        #i = np.shape(self.layers)[0]\n",
    "        #print(\"shape of i is \",i)\n",
    "        #i=i-1\n",
    "        #while i>=upto:\n",
    "            #print(\"input to gradient is \", self.layers[i-1].last_input)\n",
    "         #   a=self.layers[i-1].activation_fun\n",
    "          #  print(next_grad.shape)\n",
    "           # print(self.layers[i-1].W.shape)\n",
    "            #print(a.act_d(self.layers[i-1].last_input).shape)\n",
    "            #next_grad=np.dot(next_grad,self.layers[i-1].W,a.act_d(self.layers[i-1].last_input))\n",
    "            #i=i-1\n",
    "            \n",
    "            \n",
    "        #iterate in reverse order for backpropagation and leave out the last entry\n",
    "        n = len(self.layers)\n",
    "        for i in range(n-2, upto-1, -1):\n",
    "            next_grad = self.layers[i].bprop(next_grad)\n",
    "        \n",
    "        return next_grad\n",
    "    \n",
    "    def classification_error(self, X, Y):\n",
    "        \"\"\" Calculate error on the given data \n",
    "            assuming they are classes that should be predicted. \n",
    "        \"\"\"\n",
    "        Y_pred = unhot(self.predict(X))\n",
    "        error = Y_pred != unhot(Y)\n",
    "        return np.mean(error)\n",
    "    \n",
    "    def sgd_epoch(self, X, Y, learning_rate, batch_size):\n",
    "        n_samples = X.shape[0]\n",
    "        n_batches = n_samples // batch_size\n",
    "        \n",
    "        for b in range(n_batches):\n",
    "            # TODO #####################################\n",
    "            # Implement stochastic gradient descent here\n",
    "            # TODO #####################################\n",
    "            # start by extracting a batch from X and Y\n",
    "            # (you can assume the inputs are already shuffled)\n",
    "\n",
    "            # TODO: then forward and backward propagation + updates\n",
    "            # HINT: layer.params() returns parameters *by reference*\n",
    "            #       so you can easily update in-place\n",
    "            # Determine the sample slice for this minibatch\n",
    "            batch_start = b * batch_size\n",
    "            batch_end = (b + 1) * batch_size\n",
    "        \n",
    "            # Extract the batch from the training set\n",
    "            X_batch = X[batch_start:batch_end]\n",
    "            #print(\"X_batch size\",X_batch.shape)\n",
    "            Y_batch = Y[batch_start:batch_end]\n",
    "            #print(\"Y_batch size\",Y_batch.shape)\n",
    "         \n",
    "            Y_pred = self.predict(X_batch)\n",
    "            self.backpropagate(Y_batch, Y_pred, upto=0)\n",
    "            num_lay = len(self.layers)\n",
    "            for i in range(1, num_lay - 1):\n",
    "                w, b = self.layers[i].params()\n",
    "                dw, db = self.layers[i].grad_params()\n",
    "                w -= learning_rate * dw\n",
    "                b -= learning_rate * np.transpose(np.squeeze(db))\n",
    "                self.layers[i].W = w\n",
    "                self.layers[i].b = b\n",
    "\n",
    "                \n",
    "                \n",
    "    def adam(self, X, Y, learning_rate, batch_size,decay_rate1= 0.01,decay_rate2=0.02):\n",
    "        n_samples = X.shape[0]\n",
    "        n_batches = n_samples // batch_size\n",
    "        moment1=0\n",
    "        moment2=0\n",
    "        \n",
    "        for b in range(n_batches):\n",
    "            # TODO #####################################\n",
    "            # Implement stochastic gradient descent here\n",
    "            # TODO #####################################\n",
    "            # start by extracting a batch from X and Y\n",
    "            # (you can assume the inputs are already shuffled)\n",
    "\n",
    "            # TODO: then forward and backward propagation + updates\n",
    "            # HINT: layer.params() returns parameters *by reference*\n",
    "            #       so you can easily update in-place\n",
    "            # Determine the sample slice for this minibatch\n",
    "            batch_start = b * batch_size\n",
    "            batch_end = (b + 1) * batch_size\n",
    "        \n",
    "            # Extract the batch from the training set\n",
    "            X_batch = X[batch_start:batch_end]\n",
    "            Y_batch = Y[batch_start:batch_end]\n",
    "            Y_pred = self.predict(X_batch)\n",
    "            self.backpropagate(Y_batch, Y_pred, upto=0)\n",
    "            num_lay = len(self.layers)\n",
    "            for i in range(1, num_lay - 1):\n",
    "                w, b = self.layers[i].params()\n",
    "                dw, db = self.layers[i].grad_params()\n",
    "                \n",
    "                moment1_= decay_rate1 * moment1 + (1-decay_rate1)*dw\n",
    "                moment2_= decay_rate2 * moment2 + (1-decay_rate2)*np.multiply(dw,dw)\n",
    "                delta=moment1_*1/(np.sqrt(0.001+moment2_))     \n",
    "                w -= learning_rate*delta \n",
    "                \n",
    "                moment1= decay_rate1 * moment1 + (1-decay_rate1)*db\n",
    "                moment2= decay_rate2 * moment2 + (1-decay_rate2)*db*db\n",
    "                delta=moment1*1/(np.sqrt(0.001+moment2))                \n",
    "                \n",
    "                print(\"shape of delta\",delta.shape)\n",
    "                print(\"shape of b\",b.shape)  \n",
    "                \n",
    "                b -= learning_rate * np.squeeze(delta)\n",
    "                self.layers[i].W = w\n",
    "                self.layers[i].b = b\n",
    "                   \n",
    "\n",
    "                    \n",
    "    def gd_epoch(self, X, Y,learning_rate):\n",
    "        # TODO ##################################################\n",
    "        # Implement batch gradient descent here\n",
    "        # A few hints:\n",
    "        #   There are two strategies you can follow:\n",
    "        #   Either shove the whole dataset throught the network\n",
    "        #   at once (which can be problematic for large datasets)\n",
    "        #   or run through it batch wise as in the sgd approach\n",
    "        #   and accumulate the gradients for all parameters as\n",
    "        #   you go through the data. Either way you should then\n",
    "        #   do one gradient step after you went through the\n",
    "        #   complete dataset!\n",
    "        # TODO ##################################################\n",
    "        Y_pred = self.predict(X)\n",
    "        self.backpropagate(Y, Y_pred, upto=0)\n",
    "        num_lay = len(self.layers)\n",
    "        for i in range(1, num_lay - 1):\n",
    "                w, b = self.layers[i].params()\n",
    "                dw, db = self.layers[i].grad_params()\n",
    "                w -= learning_rate * dw\n",
    "                b -= learning_rate * np.transpose(np.squeeze(db))\n",
    "                self.layers[i].W = w\n",
    "                self.layers[i].b = b\n",
    "                   \n",
    "       \n",
    "    \n",
    "    def train(self, X, Y, X_valid, Y_valid, X_test,Y_test ,learning_rate, max_epochs, batch_size,\n",
    "              descent_type, y_one_hot=True):\n",
    "\n",
    "        \"\"\" Train network on the given data. \"\"\"\n",
    "        n_samples = X.shape[0]\n",
    "        n_batches = n_samples // batch_size\n",
    "        if y_one_hot:\n",
    "            Y_train = one_hot(Y)\n",
    "            Y_valid = one_hot(Y_valid)\n",
    "            Y_test = one_hot(Y_test)\n",
    "        else:\n",
    "            Y_train = Y\n",
    "            \n",
    "        training_loss =np.ones(max_epochs+1) * np.nan \n",
    "        validation_loss = np.ones(max_epochs+1) * np.nan \n",
    "        train_error =np.ones(max_epochs+1) * np.nan \n",
    "        valid_error = np.ones(max_epochs+1) * np.nan\n",
    "       \n",
    "        \n",
    "        print(\"... starting training\")\n",
    "        for e in range(max_epochs+1):\n",
    "            if descent_type == \"sgd\":\n",
    "                self.sgd_epoch(X, Y_train, learning_rate, batch_size)\n",
    "            elif descent_type == \"gd\":\n",
    "                self.gd_epoch(X, Y_train, learning_rate)\n",
    "            elif descent_type == \"adam\":\n",
    "                self.adam(X, Y_train, learning_rate,batch_size)\n",
    "            else:\n",
    "                raise NotImplementedError(\"Unknown gradient descent type {}\".format(descent_type))\n",
    "\n",
    "            # Output error on the training data\n",
    "            train_loss = self._loss(X, Y_train)\n",
    "            training_loss[e]=train_loss\n",
    "            train_error[e] = self.classification_error(X, Y_train)\n",
    "            te = train_error[e]\n",
    "            validation_loss[e]=self._loss(X_valid,  Y_valid)\n",
    "            v=validation_loss[e]\n",
    "            validation_error=self.classification_error(X_valid, Y_valid)\n",
    "            valid_error[e]=validation_error\n",
    "            \n",
    "            if(e%5==0):\n",
    "                print('epoch {:.1f}, train loss {:.4f}, train error {:.4f}, validation loss {:.4f}, validation error {:.4f}'.format(\n",
    "                    e, train_loss, te,v,validation_error))\n",
    "            # TODO ##################################################\n",
    "            # compute error on validation data:\n",
    "            # simply make the function take validation data as input\n",
    "            # and then compute errors here and print them\n",
    "            # TODO ##################################################\n",
    "            #validation_loss=self._loss(X_valid, Y_valid)  \n",
    "        print(\"learning rate: \", learning_rate)\n",
    "        print(\"descent type\", descent_type)\n",
    "        print(\"validation error\", valid_error[len(valid_error)-1])\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        print(\"test error\",  self.classification_error(X_test, Y_test))\n",
    " \n",
    "        #plt.plot(training_loss, c='b', label='Training loss')\n",
    "        #plt.plot(validation_loss, c='m', label='Validation loss')\n",
    "        #plt.plot(train_error, c='g', label='Training error')\n",
    "        #plt.plot(valid_error, c='r', label='Validation error')\n",
    "        #plt.legend(loc='lower left')\n",
    "        \n",
    "        return training_loss,train_error,validation_loss,valid_error\n",
    "        \n",
    "\n",
    "        \n",
    "    def check_gradients(self, X, Y):\n",
    "        \"\"\" Helper function to test the parameter gradients for\n",
    "        correctness. \"\"\"\n",
    "        for l, layer in enumerate(self.layers):\n",
    "            if isinstance(layer, Parameterized):\n",
    "                print('checking gradient for layer {}'.format(l))\n",
    "                for p, param in enumerate(layer.params()):\n",
    "                    # we iterate through all parameters\n",
    "                    #print(\"parameters of layer--\",param)\n",
    "                    param_shape = param.shape\n",
    "                    # define functions for conveniently swapping\n",
    "                    # out parameters of this specific layer and \n",
    "                    # computing loss and gradient with these \n",
    "                    # changed parametrs\n",
    "                    def output_given_params(param_new):\n",
    "                        \"\"\" A function that will compute the output \n",
    "                            of the network given a set of parameters\n",
    "                        \"\"\"\n",
    "                        # copy provided parameters\n",
    "                        param[:] = np.reshape(param_new, param_shape)\n",
    "                        # return computed loss\n",
    "                        return self._loss(X, Y)\n",
    "\n",
    "                    def grad_given_params(param_new):\n",
    "                        \"\"\"A function that will compute the gradient \n",
    "                           of the network given a set of parameters\n",
    "                        \"\"\"\n",
    "                        # copy provided parameters\n",
    "                        param[:] = np.reshape(param_new, param_shape)\n",
    "                        # Forward propagation through the net\n",
    "                        #pdb.set_trace()\n",
    "                        Y_pred = self.predict(X)\n",
    "                        # Backpropagation of partial derivatives\n",
    "                        #pdb.set_trace()\n",
    "                        self.backpropagate(Y, Y_pred, upto=l)\n",
    "                        # return the computed gradient \n",
    "                        return np.ravel(self.layers[l].grad_params()[p])\n",
    "\n",
    "                    # let the initial parameters be the ones that\n",
    "                    # are currently placed in the network and flatten them\n",
    "                    # to a vector for convenient comparisons, printing etc.\n",
    "                    param_init = np.ravel(np.copy(param))\n",
    "                    \n",
    "                    # To debug you network's gradients use scipys\n",
    "                    # gradient checking!\n",
    "                    epsilon = 1e-5\n",
    "                    import scipy.optimize\n",
    "                    err = scipy.optimize.check_grad(output_given_params, \n",
    "                          grad_given_params,\n",
    "                          param_init) \n",
    "                    print('diff scipy {:.2e}'.format(err))\n",
    "                    assert(err < epsilon)\n",
    "                    \n",
    "                    # reset the parameters to their initial values\n",
    "                    param[:] = np.reshape(param_init, param_shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient Checking\n",
    "After implementing everything it is always a good idea to setup some layers and perform gradient\n",
    "checking on random data. **Note** that this is only an example! It is not a useful network architecture ;). We also expect you to play around with this to test all your implemented components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_shape = (5, 10)\n",
    "n_labels = 6\n",
    "layers = [InputLayer(input_shape)]\n",
    "#print(InputLayer(input_shape))\n",
    "\n",
    "layers.append(FullyConnectedLayer(\n",
    "                layers[-1],\n",
    "                num_units=15,\n",
    "                init_stddev=0.1,\n",
    "                activation_fun=Activation('relu')\n",
    "))\n",
    "layers.append(FullyConnectedLayer(\n",
    "                layers[-1],\n",
    "                num_units=6,\n",
    "                init_stddev=0.1,\n",
    "                activation_fun=Activation('tanh')\n",
    "))\n",
    "layers.append(FullyConnectedLayer(\n",
    "                layers[-1],\n",
    "                num_units=n_labels,\n",
    "                init_stddev=0.1,\n",
    "                activation_fun=Activation('relu')\n",
    "))\n",
    "layers.append(SoftmaxOutput(layers[-1]))\n",
    "nn = NeuralNetwork(layers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create random data\n",
    "X = np.random.normal(size=input_shape)\n",
    "# and random labels\n",
    "Y = np.zeros((input_shape[0], n_labels))\n",
    "for i in range(Y.shape[0]):\n",
    "    idx = np.random.randint(n_labels)\n",
    "    Y[i, idx] = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import pdb; pdb.set_trace()\n",
    "\n",
    "#nn.check_gradients(X, Y)\n",
    "\n",
    "#nn.train(X,Y,X,Y,learning_rate=0.6, \n",
    "         #max_epochs=30, batch_size=64, y_one_hot=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training on MNIST\n",
    "Finally we can let our network run on the MNIST dataset!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First load the data and reshape it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "... loading data\n",
      "... done loading data\n"
     ]
    }
   ],
   "source": [
    "# load\n",
    "Dtrain, Dval, Dtest = mnist()\n",
    "X_train, y_train = Dtrain\n",
    "X_valid, y_valid = Dval\n",
    "X_test, y_test = Dtest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Dtrain* contains 50k images which are of size 28 x 28 pixels. Hence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_train shape: (50000, 1, 28, 28)\n",
      "y_train shape: (50000,)\n"
     ]
    }
   ],
   "source": [
    "print(\"X_train shape: {}\".format(np.shape(X_train)))\n",
    "print(\"y_train shape: {}\".format(np.shape(y_train)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "y_train will automatically be converted in the *train()* function to one_hot encoding.\n",
    "\n",
    "\n",
    "But we need to reshape X_train, as our Network expects flat vectors of size 28*28 as input!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reshaped X_train size: (50000, 784)\n",
      "Reshaped X_valid size: (10000, 784)\n",
      "Reshaped X_test size: (10000, 784)\n"
     ]
    }
   ],
   "source": [
    "X_train = X_train.reshape(X_train.shape[0], -1)\n",
    "print(\"Reshaped X_train size: {}\".format(X_train.shape))\n",
    "X_valid = X_valid.reshape((X_valid.shape[0], -1))\n",
    "print(\"Reshaped X_valid size: {}\".format(X_valid.shape))\n",
    "X_test = X_test.reshape((X_test.shape[0], -1))\n",
    "print(\"Reshaped X_test size: {}\".format(X_test.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ah, much better ;-)! \n",
    "\n",
    "Now we can finally really start training a Network!\n",
    "\n",
    "\n",
    "I pre-defined a small Network for you below. Again This is not really a good default and will not produce state of the art results. Please play around with this a bit. See how different activation functions and training procedures (gd / sgd) affect the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot(tl,vl,max_epochs):\n",
    "    fig, ax = plt.subplots(1, 1)\n",
    "    ax.set_title('Training progress')\n",
    "    ax.set_xlim([0, max_epochs])\n",
    "    ax.set_ylim([0, 1])\n",
    "    ax.set_xlabel('Epoch')\n",
    "    ax.set_ylabel('Error/loss')\n",
    "\n",
    "    if ax.lines:\n",
    "        ax.lines[0].set_ydata(tl)\n",
    "        ax.lines[1].set_ydata(vl)\n",
    "        ax.set_ylim([0,  1])\n",
    "    else:\n",
    "        ax.plot(tl, c='b', label='Training loss')\n",
    "        ax.plot(vl, c='m', label='Validation loss')\n",
    "        ax.legend(loc='upper right')\n",
    "    fig.canvas.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_graphs(tl, te, vl, ve):    \n",
    "    epochs = np.arange(tl.shape[0])\n",
    "    fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(12, 12))\n",
    "    y = [tl, vl, te, ve]\n",
    "    lab = [\"Training Loss\", \"Validation Loss\", \"Training Error\", \"Validation Error\"]\n",
    "    #y=[tl,vl]\n",
    "    #lab = [\"Training Loss\", \"Validation Loss\"]\n",
    "    i = 0\n",
    "    for row in ax:\n",
    "        for col in row:\n",
    "            col.plot(epochs, y[i])\n",
    "            col.grid(True)\n",
    "            col.set_xlabel(\"EPOCHS\")\n",
    "            col.set_ylabel(lab[i])\n",
    "            # col.annotate(s=y[i], xy=zip(epochs, y[i]))\n",
    "            i += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "# Setup a small MLP / Neural Network\n",
    "# we can set the first shape to None here to indicate that\n",
    "# we will input a variable number inputs to the network\n",
    "input_shape = (None, 28*28)\n",
    "layers = [InputLayer(input_shape)]\n",
    "\n",
    "layers.append(FullyConnectedLayer(\n",
    "                layers[-1],\n",
    "                num_units=500,\n",
    "                init_stddev=0.01,\n",
    "                activation_fun=Activation('relu')\n",
    "))\n",
    "\n",
    "layers.append(FullyConnectedLayer(\n",
    "                layers[-1],\n",
    "                num_units=300,\n",
    "                init_stddev=0.01,\n",
    "                activation_fun=Activation('relu')\n",
    "))\n",
    "layers.append(FullyConnectedLayer(\n",
    "                layers[-1],\n",
    "                num_units=10,\n",
    "                init_stddev=0.01,\n",
    "                # last layer has no nonlinearity \n",
    "                # (softmax will be applied in the output layer)\n",
    "                activation_fun=None\n",
    "))\n",
    "layers.append(SoftmaxOutput(layers[-1]))\n",
    "\n",
    "nn = NeuralNetwork(layers)\n",
    "# Train neural network\n",
    "t0 = time.time()\n",
    "tl, te, vl, ve=nn.train(X_train, y_train, X_valid, y_valid,X_test, y_test, learning_rate=0.13, \n",
    "         max_epochs=250, batch_size=512, y_one_hot=True,descent_type=\"sgd\")\n",
    "t1 = time.time()\n",
    "print('Duration: {:.1f}s'.format(t1-t0))\n",
    "print(\"2-Layer network with 500+300 HU\")\n",
    "draw_graphs(tl, te, vl, ve)\n",
    "plot(tl,vl,max_epochs=250)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Figure out a reasonable Network that achieves good performance\n",
    "As the last part of this task, setup a network that works well and gets reasonable accuracy, say ~ 1-3 percent error on the **validation set**. \n",
    "Train this network on the complete data and compute the **test error**. \n",
    "\n",
    "Visualize the validation loss and training loss for each iteration in a plot, e.g. using matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}