Exercise4/exercise4_20190226.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# used for manipulating directory paths\n",
    "import os\n",
    "\n",
    "# Scientific and vector computation for python\n",
    "import numpy as np\n",
    "\n",
    "# Plotting library\n",
    "from matplotlib import pyplot\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Optimization module in scipy\n",
    "from scipy import optimize\n",
    "\n",
    "# will be used to load MATLAB mat datafile format\n",
    "from scipy.io import loadmat\n",
    "\n",
    "# library written for this exercise providing additional functions for assignment submission, and others\n",
    "import utils\n",
    "\n",
    "# define the submission/grader object for this exercise\n",
    "grader = utils.Grader()\n",
    "\n",
    "# tells matplotlib to embed plots within the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "\n",
    "Here I am writing what I have learned in great nitty gritty details about neural networks deep learning. The contents are mostly based on Andrew Ng's machine learning course on Coursera (https://www.coursera.org/learn/machine-learning). I strongly believe that you need to implement algorithms such as feedforward, back propagation, gradient descend, and cost function yourself in order to truly understand the inner working of a neural network. I hope through this writing I will deepen my understanding of this topic and others may be able to pick up one insight or two.\n",
    "\n",
    "### Objective\n",
    "Here, I use a small subset of MNIST data set as examples and aim to build a simply neural network to recognize those hand-written digits. To really understand how neural network works, we must understand what forward pass and back propagation are doing. I believe the most effective way to gain the intuition is to implement the key steps such as the cost function, feed forwrad, back propagation process. \n",
    "\n",
    "### neural network model representation\n",
    "\n",
    "The neural network is shown in the following figure. Briefly, This simple neural network has 3 layers: input layer, hidden layer, and output layer\n",
    "* Input layer l1 has 401 units (400 image pixels + 1 bias). Input X dimension is (5000, 401). Each row are pixel gray scale values for a training example. Each column are pixels for a specific square in the 20x20 grid for all training examples.\n",
    "\n",
    "* Hidden layer l2 has 26 units (25 units derived from l1 + 1 bias). \n",
    "* Weight dimension for input layer is (25, 401). Each row are 401 weights for each input unit connecting to a specific unit in the hidden layer. Each column are weights for a specific input unit connecting to all units in hidden layer. \n",
    "\n",
    "* The output layer l3 has 10 units, with each unit representing the probability of being a hand-written number 0 to 9.\n",
    "\n",
    "* Weights dimension for hidden layer is (10, 26). Each row are  26 weights for each of the 26 units in hidden layer transforming into a specific unit in output layer. Each column are weights for a specific unit in hidden layer transforming into 10 units in the output layer.\n",
    "\n",
    "\n",
    "shapes of neural network variables:\n",
    "\n",
    "|Items|dimension_x|dimension_y|\n",
    "|------|------|------|\n",
    "|a1|\t5000|\t401|\n",
    "|theta1|\t26|\t401|\n",
    "|a2|\t5000|26\n",
    "|theta2|\t10|26|\n",
    "|a3\t|5000|10|\n",
    "|J|\t1||\t\t\n",
    "|d1|26||\n",
    "|d2|10||\n",
    "|theta1_grad|\t25|401|\n",
    "|theta2_grad|\t10|\t26|\n",
    "\n",
    "\n",
    "\n",
    "cost is the summation of all costs: costs of each training example to all k labels, and costs of all training examples\n",
    "\n",
    "\n",
    "![](Figures/neural_network.png)\n",
    "\n",
    "### forward pass\n",
    "The forward pass happens when we feed the input signal to the neural network and let it predict the 10 probabilities of the image being each of the 10 numbers. The input signal are values of 400 gray scale intensity of an hand-written image in this case. \n",
    "\n",
    "What really happens in forward pass is just some simple mathatical operations, namely matrix multiplication and activation. This a two-step process. First, we compute the weighted sum of all the units in previous layer for all the units in current layer.  Then we feed the sum into a activation function. In this case, we choose sigmoid function, which squashes the numbers into a value between 0 and 1. Ideally, we would be better off to use softmax function for the last layer so that the sum of all the probabilities is equal to 1. \n",
    "\n",
    "Let's get to the math behind.\n",
    "\n",
    "The implementation of forward pass can be as simple as a few lines of code. In below code snippet,  X is input; a2 is activation in hidden layer; a3 is activation in final output layer.\n",
    "    \n",
    "    a1 = np.concatenate([np.ones((X.shape[0], 1)), X], axis=1) # add bias to input matrix\n",
    "    a2 = sigmoid(a1 @ Theta1.T) # weight multiple input, and then apply activate function    \n",
    "    a2 = np.concatenate([np.ones((a2.shape[0], 1)), a2], axis=1) # add bias to activation in hidden layer\n",
    "    a3 = sigmoid(a2 @ Theta2.T) # weight multiple input, and then apply activate function\n",
    "\n",
    "### back propagation\n",
    "Back propagation is basically the opposite of forward pass. This process propagates the prediction error backward to each of the previous layers except the input layer so that we could compute the gradients for the weights and then finally update the weights in order to minimize the cost function.\n",
    "\n",
    "#### calculate deltas (prediction errors)\n",
    "We could think deltas as prediction errors. It is very straight forward to calculate the delta for the ouput layer (delta2). We have predicted probabilities and the true labels of all input images. Therefore, the error (delta2) is just simply the difference of these two vectors. Let's use one single training image as an example, in our case, predicted vector is the probabilities for each of the 10 numbers and the label is a one-hot encoded vector with 10 elements, which are all zeros except the elements corresponding to the label number. For instance, if the input image is 0, which is one-hot encoded as y, and the predictions is a vector (a3). Detla2 is a vector with shape (10,).\n",
    "\n",
    "$$ y = \n",
    "\\begin{bmatrix} 1 \\\\ 0 \\\\ 0 \\\\\\vdots \\\\ 0 \\end{bmatrix}, \\quad\n",
    "a3 = \n",
    "\\begin{bmatrix} 0.92 \\\\ 0.01 \\\\ 0.02 \\\\ \\vdots \\\\ 0 \\end{bmatrix}, \\quad \n",
    "delta2 = a3 - y =\n",
    "\\begin{bmatrix} -0.08 \\\\ 0.01 \\\\ 0.02 \\\\ \\vdots \\\\ 0 \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "delta2 can be computed as follow:\n",
    "\n",
    "    # a3 is the prediction; y is the label\n",
    "    delta2 = a3 - y # prediction error\n",
    "\n",
    "Now we need to back propagate the error from output layer into the hidden layer.The error of hidden layer (delta1) is the product of 2 components. One is the weighed sum of delta2 backward to hidden layer, and the other is the gradient of each activation in hidden layer. \n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "##### component 1: weighted sum of delta2 (WSD)\n",
    "For each unit in hidden layer, this can be calculated by computing the product of weights(W2) of the unit connecting to each of the 10 output units and the error in output layer (delta2). For example\n",
    "\n",
    "$$\n",
    "W2 = \n",
    "\\begin{bmatrix} -0.82 \\\\ 0.12 \\\\ -0.20 \\\\ \\vdots \\\\ 0 \\end{bmatrix}, \\quad\n",
    "delta2 = \n",
    "\\begin{bmatrix} -0.08 \\\\ 0.01 \\\\ 0.02 \\\\ \\vdots \\\\ 0 \\end{bmatrix}, \\quad\n",
    "weighted\\_sum = 0.0628 \n",
    "$$\n",
    "\n",
    "Weighted sum of delta2 can be computed as follow:\n",
    "    WSD = Theta2.T @ delta2 \n",
    "\n",
    "##### component 2: sigmoid gradient of activation (SGA)\n",
    "This is basically to compute the gradient of the activation for each unit in hidden layer. Let's use one training image to walk through the process. First, we compute the weighted sum by multiplying the input pixel values (400 + 1 bias) and their corresponding weights and summing them for each unit in hidden layer. Then we calculate the activation by applying sigmoid function to the weighted sums for each unit to squash them to values between 0 and 1. Lastly, We compute the derivative of the activation, which is just the product of activation and (1 - activation).\n",
    "\n",
    "Sigmoid function is given by this formula:\n",
    "$$\n",
    "\\sigma(z) = \\frac{1}{1 + e^{-z}} = \\frac{e^{z}}{e^{z} + 1}\n",
    "$$\n",
    "\n",
    "The gradient for the sigmoid function can be computed as following:\n",
    "$$\n",
    "\\sigma'(z) = \\frac{\\partial \\sigma}{\\partial z} = \\sigma(z)(1 - \\sigma(z))\n",
    "$$\n",
    "\n",
    "This code snippet shows how this can be accomplished:\n",
    "\n",
    "    # gradient for each activation in layer 2\n",
    "    SGA = sigmoidGradient(np.hstack((1, X1 @ Theta1.T)))\n",
    "    \n",
    "Then put the two component together:\n",
    "\n",
    "    delta1 = WSD * SGA\n",
    "\n",
    "\n",
    "##### compute the weight gradients\n",
    "To calculate the gradients, you basically just multiply the error in output layer (delta2) by the activations in its previous layer. To be more explicit, the error for the first unit in the output layer multiplying each of the 26 activations in the hidden layer gives the gradient for each weight corresponding to the first unit in the output layer. We repeat the calculations for errors of other remaining units.\n",
    "\n",
    "Similarly, we could multiply the error in hidden layer (delta1) by the activations in its previous layer (they are just input pixel values in this case since the previous layer is input layer) to calculate the gradients for weights connecting input layer and hidden layer (Theta1).\n",
    "\n",
    "We need to do this for all training examples (or a mini batch of training examples) and then calculate the average gradients, which we could use to update the Theta1 and Theta2 (weights).\n",
    "\n",
    "\n",
    "### Intutition of back propagation\n",
    "When a traning image forward passes the network and generates a prediction, the error between the prediction and real label can be calculated. In order for the network to learn and produce better prediction, the network needs to learn from this training process. Depending on how bad the prediction is, the network either has to dramatically adjust the weights or just fine tuning them. This is where the back propagation comes in play.\n",
    "\n",
    "How do the prediction errors affect the weight adjustment?\n",
    "Intuitively, the bigger the error, the deeper the adjustment needs to be made. Recall the weight adjustment (gradient) is the product of error and the gradient of the activation in the previous layer. So indeed the greater the error, the steeper adjustment we need to do for the weights.\n",
    "\n",
    "How do the activation values affect the weight adjustment?\n",
    "By looking at the sigmoid plot, We know that the gradient reaches a maximum of 0.25 when the activation is 0.5. The gradients gradually decline when the activation is moving away from 0.5 in both directions. So this is equivalent to say that we adjust the weights associated with an activation of 0.5 more when others being equal. This makes sense to me because 0.5 basically means that we are absolutely no idea if this unit is supposed to be activated or not. In the final output layer, 0.5 means the network thinks there is 50% chance the input image is belong to that class. Therefore, the network will make relative large adjust to weights associated with this node. On the other hand, if the activation is close to 1 or 0, it is very confident about the activation or prediction, the weights associated with this node only need to be fine tuned.\n",
    "\n",
    "How do the weights affect the adjustment of themselves?\n",
    "The gradient for a weight is the product of the error (delta) and the weight itself. So we do bigger adjustment for large weights if others being equal.   \n",
    "\n",
    "With these intuitions, It is easier to understand that neural network is capable of learning very complex patterns and mapping input to output with universal approximation. \n",
    "\n",
    "If you notice that I misunderstand something here, Please point out and constructively criticize. I appreciate your feedback.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "### some variables \n",
    "\n",
    "k: the number of labels/classes, in this case k=10, representing numbers between 0 and 9\n",
    "\n",
    "m: the number of training examples, in this case m=5000\n",
    "Theta1: weights for input layer\n",
    "Theta2: weights for hidden layer\n",
    "\n",
    "\n",
    "## nnCostFunction implementation notes:\n",
    "\n",
    "### 1) One-hot encoding of y:\n",
    "y_encode is also a 5000x10 matrix, each row is a encoded label for a training iage. If the label is 2, the 3rd element in the row is 1 and the rest 9 elements are all 0s.\n",
    "$$ y = \\begin{bmatrix} 0 & 0 & 1 & 0 \\cdots & 0 \\end{bmatrix}$$\n",
    "\n",
    "### 2) feedforwad outputs (activation of layer 3: a3): \n",
    "a3 is a 5000x10 matrix. Each row is a prediction for a training image. It shows the probability of being each number, 0 to 9. For instance, following y_hat means the predicted digit is 2 with 90% chance and is 0 with 10% probability.\n",
    "$$ y_{hat} = \\begin{bmatrix} 0.1 & 0 & 0.9 & 0 & \\cdots & 0 \\end{bmatrix}$$\n",
    "\n",
    "\n",
    "### 3) calculate the cost (loss),\n",
    "The cost function without regularization is:\n",
    "$$ J(\\theta) = \\frac{1}{m} \\sum_{i=1}^{m}\\sum_{k=1}^{K} \\left[ - y_k^{(i)} \\log \\left( \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) - \\left( 1 - y_k^{(i)} \\right) \\log \\left( 1 - \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) \\right]$$\n",
    "\n",
    "The cost function for neural networks with regularization is given by:\n",
    "$$ J(\\theta) = \\frac{1}{m} \\sum_{i=1}^{m}\\sum_{k=1}^{K} \\left[ - y_k^{(i)} \\log \\left( \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) - \\left( 1 - y_k^{(i)} \\right) \\log \\left( 1 - \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) \\right] + \\frac{\\lambda}{2 m} \\left[ \\sum_{j=1}^{25} \\sum_{k=1}^{400} \\left( \\Theta_{j,k}^{(1)} \\right)^2 + \\sum_{j=1}^{10} \\sum_{k=1}^{25} \\left( \\Theta_{j,k}^{(2)} \\right)^2 \\right] $$\n",
    "use the cost function formula? cost for every class and cost for every training image\n",
    "\n",
    "4) the dimensions do not look right, why does it still work?\n",
    "\n",
    "5) * is element_wise operation, not matrix multiplication @\n",
    "\n",
    "6) for the cost, we only need element-wise multiplication, what is the ituition. when we predicted the exact label at 100% confidence, the cost is 0. when we predict wrong label with 100% confidence, the cost is infinity.  \n",
    "\n",
    "Sigmoid function squashes input (a number) into a value, which is between 0 and 1. We use sigmoid function as activation function because we would like to constrain the activation. When it is close to 1, we are saying the unit is activated. When is close to 0, we are saying the unit is not activated.\n",
    "\n",
    "Sigmoid function is given by this formula:\n",
    "$$\n",
    "\\sigma(z) = \\frac{1}{1 + e^{-z}} = \\frac{e^{z}}{e^{z} + 1}\n",
    "$$\n",
    "\n",
    "The gradient for the sigmoid function can be computed as following:\n",
    "$$\n",
    "\\sigma'(z) = \\frac{\\partial \\sigma}{\\partial z} = \\sigma(z)(1 - \\sigma(z))\n",
    "$$\n",
    "\n",
    "The intuition is that the gradient of a sigmoid function is equal to the sigmoid function times 1 minus the sigmoid function iteself().\n",
    "\n",
    "\n",
    "\n",
    "Now complete the implementation of `sigmoidGradient` in the next cell.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "# neural network structure\n",
    "input_size  = 400  # 20x20 Input Images of Digits\n",
    "hidden_size = 25   # 25 hidden units\n",
    "num_labels = 10    # 10 labels, from 0 to 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5000, 400)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(5000,)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# training images\n",
    "data = np.loadtxt('MNIST_data.csv')\n",
    "X = data[:, :-1]\n",
    "y= data[:,-1]\n",
    "X.shape\n",
    "y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "# each label needs to one-hot encoded into a vector with a shape of (10,)\n",
    "def one_hot(x):\n",
    "    y = np.zeros((x.shape[0], num_labels))\n",
    "    for i, j in enumerate(x.astype(int)):\n",
    "        y[i, j] = 1\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### one-hot encode image labels (y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0.])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_encode = one_hot(y)\n",
    "y[:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0., ..., 0., 0., 0.],\n",
       "       [1., 0., 0., ..., 0., 0., 0.],\n",
       "       [1., 0., 0., ..., 0., 0., 0.],\n",
       "       ...,\n",
       "       [0., 0., 0., ..., 0., 0., 1.],\n",
       "       [0., 0., 0., ..., 0., 0., 1.],\n",
       "       [0., 0., 0., ..., 0., 0., 1.]])"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_encode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10285,)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([-2.25623899e-02, -1.05624163e-08])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load weights previously learned\n",
    "nn_params = np.loadtxt('weights.csv')\n",
    "nn_params.shape\n",
    "nn_params[:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(25, 401)"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(10, 26)"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Reshape nn_params back into the parameters Theta1 and Theta2\n",
    "Theta1 = nn_params[:((input_size+1) * hidden_size)].reshape(hidden_size,input_size+1)\n",
    "Theta2 = nn_params[((input_size +1)* hidden_size ):].reshape(num_labels,hidden_size+1)\n",
    "Theta1.shape\n",
    "Theta2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2MAAAEyCAYAAAB3fAhyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuYnFWZ7/3fnQ7dnRAICBrlNBwmKLy6XznrdrYCggTHAQYHhOy5BASzHbeoIw6COEhEVHTE8ZWwte0GxNkDnkbNRPCAEAZh0LQMZjiTASUN4RxOYncn6fv9oyqxuruqnlVVz2l1fT/XVVeqq+5atfpQv6z1HNZj7i4AAAAAQL5mFd0BAAAAAOhGTMYAAAAAoABMxgAAAACgAEzGAAAAAKAATMYAAAAAoABMxgAAAACgAEzGAAAAAKAATMYAAAAAoABMxgAAAACgALOzfoPx8XHP+j2AmPT29lor9WbW1mfI3Vt6n240NjZGPgE1+vr6Ms8nsinZ6Ogo2QTU6O/vn7FjJ/aMAQAAAEABMt8zBqAzZmxEBlBO5BOAMoopm5iMASUXU6AA6C7kE4AyiimbmIwBJRdToADoLuQTgDKKKZsSJ2Nm9hpJx0raWZJLelTScne/J+O+AZA0axandtZDNgHFI5/qI5+AYsWUTU17amYfk3SNJJP0K0mrqvevNrNzsu8eADNr6zaTkU1AOZBN05FPQPFiGjuZe+OVH83sfkn/j7tvmPJ4r6S73H1hg9ctkbREkpYtW3bAGWeckV6Pgci1urR9f39/W8uzjo6OzthRT7vZVK3Zkk+XXnop+QTUaHVp+3byaSZnk5TO2OnSSy894PTTT8+8r0AsWl3aPqaxU9JhihOSdpL0uymPv6r6XF3uPiBpQOI6Y0CnstpSY2aXS3qHpCfc/bV1njdJX5b0dkkvSTrV3W/PpDOtayubpMn5xHXGgM50w56uNnQ8duI6Y0BnYsqmpMnYhyX93MwekLS2+thukv5U0gey7BiAigwD5UpJl0q6qsHzR0taWL0dIun/VP8tA7IJKIGYBjw5Ip+AgsWUTU0nY+7+YzPbW9LBqpyEapJGJK1y90059A/oelkFirv/m5nt3qTkWElXeeVY5tvMbDsze5W7r8ukQy0gm4ByiGnAkxfyCSheTNmUuJqiu09Iui2HvgCoo8AVgXbWH7fqSpXBxM6SCp+MSWQTUAYxrViWJ/IJKFZM2RRPT4Eu1e6KQGa2xMyGa25LWn3rOo9xHgOALbJarczMLjezJ8zszgbPm5n9f2a2xsxWm9n+qX5jAKKW1WqKWWQTF30GSq7dXe21J4O3aUTSrjVf76LKtXIAQBLntAIop5iyiT1jQMkVeK2M5ZLeXd3K8wZJz5XhfDEA5ZFVNrn7v0l6pknJlnNa3f02SduZ2atS+JYAzABZjZ2yyCb2jAEll9XWHTO7WtKhknY0sxFJn5S0lSS5+1clXavKsvZrVFna/rRMOgIgWu3kk9VcT6tqoLonvxWlPqcVQLEKXMCj5WxiMgaUXFaB4u4nJzzvkv53Jm8OYEZoJ59SOIRa4pxWAE20O3ZKYWNRy9nEZAwouZhWBALQXQrMJ85pBdBQu9lUxPn2jPIAAEBsOKcVQBm1nE1t7xkzs9Pc/Yp2Xw8gTIHHPUeLfALywTmtrSGbgHzElE1WOS2krc487O67NXhuy/GWy5YtO+CMM85o6z2Amai3t7elhHjFK17R1of0iSee6NpZXGg+XXrppeQTUKOvry/zfCKbwrLp9NNPz7VvQJn19/fP2LFT08mYma1u9JSkvd29L+kNxsfHOaEWqNHqZGzBggVtfYYef/zxGT3gSSOfxsbGyCegRquTsXbyiWxKzqbR0VGyCajR6mQsprFT0mGKCyQdJWn9lMdN0q2Z9AjAJBym2BD5BBSMfKqLbAIKFlM2JU3GVkia5+53TH3CzFZm0iMAk8QUKDkjn4CCkU91kU1AwWLKpqaTMXdveMCyuy9OvzsApmJp+/rIJ6B45NN0ZBNQvJiyieuMASUX09YdAN2FfAJQRjFlE5MxoORiChQA3YV8AlBGMWUTkzGg5GIKFCArPT09QXVpf14mJiZSrZtpyCcAZRRTNjEZA0oupkAB0F3IJwBlFFM2MRkDSi6mQAHQXcgnAGUUUzYxGQNKLqYVgQB0F/IJQBnFlE1MxoCSi2nrDoDuQj4BKKOYsilx2mhmrzGzt5rZvCmPL8quWwA2M7O2bjMd2QQUj2yqj3wCihXT2KnpZMzMPijph5LOlHSnmR1b8/RnsuwYgIqYAiUvZBNQDmTTdOQTULyYxk5Je8beK+kAdz9O0qGS/t7MPlR9rmGPzWyJmQ2b2fDg4GA6PQW6VEyBkqO2skkin4A0kU11dTx2GhoayqGbwMwV09gp6ZyxHnd/UZLc/bdmdqik75rZn6hJoLj7gKQBSRofH/eU+goAm7WVTdX6Lfk0NjZGPgFIW8djp9HRUbIJ6BJJe8YeM7PXb/6iGi7vkLSjpNdl2TEAFbNmzWrrNsORTUAJkE11kU9AwWIaOyXtGXu3pI21D7j7RknvNrOvZdYrAFt0yWE9rSKbgBIgn+oin4CCxZRNTSdj7j7S5Llb0u8OgKliCpS8kE1AOZBP05FPQPFiyiauMwaUXJcc1gMgQuQTgDKKKZuYjHWJ0C0EMW1JaGRiYqLoLqRqJvxO0J3S/Nv93e9+F1T3+OOPp/aekrTHHnsE1e24445BdeQTAGQvpmxiMgaUXExbdwB0F/IJQBnFlE1MxoCSi2nrDoDuQj4BKKOYsonJGFByMQUKgO5CPgEoo5iyickYUHIx7WoH0F3IJwBlFFM2MRkDSi6mrTsAugv5BKCMYsqmxMmYmR0syd19lZntK2mRpHvd/drMewcgqq07eSKbgOKRT/WRT0CxYsqmppMxM/ukpKMlzTazn0k6RNJKSeeY2X7uflH2XQS6W0xbd/JCNgHlQD5NRz4BxYspm5KmjX8l6U2S3izpf0s6zt0/JekoSe9q9CIzW2Jmw2Y2PDg4mFpngW40a9astm4zXFvZJJFPQJrIpro6HjsNDQ3l01Nghopp7JR0mOJGd98k6SUz+y93f16S3P0PZtbwypXuPiBpQJLGx8c9td4CXSimrTs5aiubqjVb8mlsbIx8AjpAPtXV8dhpdHSUbAI6EFM2JU3Gxs1srru/JOmAzQ+a2XxJTQc8ANIRU6DkiGwCSoB8qot8AgoWUzYlTcbe7O5jkuTutQGylaRTMusVgC265LCeVpFNQAmQT3WRT0DBYsqmppOxzWFS5/GnJD2VSY8AIAHZBKCsyCcAreA6Y0DJxbSrHUB3IZ8AlFFM2cRkrKRC/4hCd8O++OKLQXUPPfRQUJ172LnFoXUh3+/cuXOD2tptt92C6mbPDvvzD/0eshLTrnZ0h9B8mpgIOz3mzjvvTKw588wzg9pas2ZNUF2oQw45JKjuq1/9alDdggULgupCf3ZFI5+A8Ews8wQh7XFd0WLKJiZjQMmVObwBdDfyCUAZxZRNTMaAkospUAB0F/IJQBnFlE1MxoCSi2lXO4DuQj4BKKOYsimengJdyszaugW2vcjM7jOzNWZ2Tp3ndzOzG83sP8xstZm9PfVvEEC0ssomAOhElmOntLU8GTOzq7LoCID6Zs2a1dYtiZn1SFom6WhJ+0o62cz2nVL2CUnfdvf9JJ0k6bKUv73UkE1A/rLIppmGbALyl9XYKQtND1M0s+VTH5J0mJltJ0nufkxWHQNQkeGWmoMlrXH3B6vvc42kYyXdXVPjkrat3p8v6dGsOtMKsgkoh6zyycwWSfqypB5Jg+7+uSnP7ybpG5K2q9ac4+7XZtKZFpBNQDnElE1JU8BdJD0v6RJJX6zeXqi536ijS8xs2MyGBwcHE94CQDPt7mqv/RxWb0umNL2zpLU1X49UH6t1gaS/NrMRSddKCltfPHttZZNEPgFpyuIwoMj32qeSTUNDQ5l3FJjJsjhMMatsSlrA40BJH5J0nqS/c/c7zOwP7n5Tsxe5+4CkAUkaHx+P44IEQEm1u9u89nPYQL3Umfp5PVnSle7+RTN7o6Rvmtlr3b3oiyC1lU3S5J/L2NgY+QR0IKPDeqLda6+Usml0dJRsAjoQUzY1nYxVB1xfMrPvVP99POk1ANKV4WGKI5J2rfl6F00PjdMlLZIkd/93M+uXtKOkJ7LqVAiyCSiHdvKpupe+dk/9QHUislm9vfZTr759gaSfmtmZkraWdETLHckA2QSUQ7tjp4R8yiSbggLC3UcknWBmf67K7ncAOcnwhNJVkhaa2R6SHlFld/riKTUPS3qrpCvNbB9J/ZKezKpDrSKbgGK1k08zfK+9JLIJKFpGRxVlkk0tba1x9x9J+lErrwHQmaz2jLn7RjP7gKSfqHKS6eXufpeZfUrSsLsvl3SWpK+b2d+qEjinunvpDp8hm4BiZJRP0e61n4psAooRUzax6xwouSyXWq2u8HPtlMfOr7l/t6Q3ZdYBAFHLKJ+i32sPoFgxZROTsQKEztZfeumlxJpbbrklqK0LLrggqO6ee+4JqpszZ05Q3dy5c4Pqnn766dTa+uIXmy5YtcWJJ54YVFe0oi5CiO7T09MTVHf//fcH1Z177rlBdTfdlLi2gSTp+eeTj/baf//9g9p69tlng+p+9KOwnRrnnXdeUN1ll4Ut+tfb25tYMzFR/BF5WeTTTNprj7iFDuhDP4sh47pW2gsR+hkNHdeF/kyK/jjGlE1Mxkoq9AOLmY/JGBA2EUP+MjyMmr32ANoWUzYVc6lpAAAAAOhy7BkDSo49YwDKinwCUEYxZROTMaDkYgoUAN2FfAJQRjFlE5MxoORiChQA3YV8AlBGMWVTS5MxM/szSQdLutPdf5pNlwDUiilQikI2AcUgn5KRT0D+Ysqmpgt4mNmvau6/V9KlkraR9EkzOyfjvgFQJVDauc1kZBNQDmTTdOQTULyYxk5JqyluVXN/iaQj3X2ppLdJ+p+NXmRmS8xs2MyGBwcHU+gm0L1iCpQctZVNEvkEpIlsqqvjsdPQ0FDWfQRmtJjGTkmHKc4ys+1VmbSZuz8pSe7+ezPb2OhF7j4gaUCSxsfHuQgj0IGMriIfu7ayqVqzJZ/GxsbIJ6AD5FNdHY+dRkdHySagAzFlU9JkbL6kX0sySW5mr3T3x8xsXvUxABnrki3JrSKbgBIgn+oin4CCxZRNTSdj7r57g6cmJP1l6r0BME1MgZIXsgkoB/JpOvIJKF5M2dTW0vbu/pKkh1LuC4A6YgqUopFNQL7Ip3DkE5CfmLKJ64wBJRdToADoLuQTgDKKKZuYjKWop6cnqO65555LrLnwwguD2vrKV74SVHfQQQcF1f393/99UN0hhxwSVDd//vyguq9//euJNVdccUVQW5/+9KeD6o4++uigutDvYWJiIqiuVTEFCspp9uywqF+zZk1Q3RlnnBFUd9dddwXVLVu2LLGmv78/qK0DDzwwqG50dDSo7vzzzw+qu+6664Lqvv/97wfVnXDCCYk1odngnt1aEOQT8hC6GEPo32No3fPPPx9UFzo+GR4eDqrbuLHpOlQtCR2bnHjiiUF1xx13XFBd0fkUUzYxGQNKLqZAAdBdyCcAZRRTNjEZA0oupkAB0F3IJwBlFFM2MRkDSi6mQAHQXcgnAGUUUzYxGQNKLqZAAdBdyCcAZRRTNjU9I9LMDjGzbav355jZUjP7VzO72MzCVjUA0BEza+s2k5FNQDmQTdORT0DxYho7JS1Pc7mkl6r3v6zKVeUvrj4WtnQMAKSPbAJQVuQTgGBJk7FZ7r55fc0D3f3D7v4Ld18qac9GLzKzJWY2bGbDg4ODqXUW6EYxbd3JUVvZJJFPQJrIpro6HjsNDQ3l01Nghopp7JR0ztidZnaau18h6TdmdqC7D5vZ3pI2NHqRuw9IGpCk8fHx7C5wAnSBLhm8tKqtbJIm59PY2Bj5BHSAfKqr47HT6Ogo2QR0IKZsStozdoakt5jZf0naV9K/m9mDkr5efQ5AxmLaupMjsgkoAbKpLvIJKFhMY6eme8bc/TlJp5rZNqrsWp8tacTdH8+jcwDi2rqTF7IJKAfyaTryCSheTNkUtLS9u78g6TcZ9wVAHTEFSt7IJqBY5FNj5BNQnJiyieuMASUXU6AA6C7kE4AyiimbmIwBJRdToADoLuQTgDKKKZuYjAVwD1vUaGRkJKju/e9/f2LNqlWrgtoKXZr7+OOPD6rr7+8Pqttqq62C6kI/DB//+McTa37+858HtXXfffcF1Y2OjgbVzZ9f7DU6YwoU5GvWrKQ1mCoeeOCBoLr3vve9QXW//vWvg+rOPffcoLrFixcn1oR+DjZt2hRUN3t22H9/hx9+eFDdddddF1T36U9/Oqhu0aJFiTWh2RT6f1g7yCc0Evq30dPTk1izevXqoLZuueWWoLrQ8dqNN94YVHfvvfcG1YXm09jYWGLNxMREUFtz5swJqrv55puD6p5//vmgulNOOSWoLqt8iimbmIwBJRdToADoLuQTgDKKKZuYjAElF1OgAOgu5BOAMoopm5iMASUXU6AA6C7kE4AyiimbmIwBJRdToADoLuQTgDKKKZuanv1tZh80s13z6gyA6WK6inxeyCagHMim6cgnoHgxjZ2SluK6UNIvzexmM3u/mb08j04B+KOYAiVHZBNQAmRTXeQTULCYxk5Jk7EHJe2iSrAcIOluM/uxmZ1iZts0epGZLTGzYTMbDl16HUB9MQVKjtrKJol8AtJENtXV8dhpaGgor74CM1JMY6ekc8bc3Sck/VTST81sK0lHSzpZ0j9Iqru1x90HJA1I0vj4eHYXOAG6QOi1pLpMW9lUfeGWfBobGyOfgA6QT3V1PHYaHR0lm4AOxJRNSZOxSVNEd98gabmk5WYWdhU5AEgf2QSgrMgnAMGSJmPvavSEu/8h5b4AqKNLDutpFdkElAD5VBf5BBQspmxqOhlz9/vz6giA+mIKlLyQTUA5kE/TkU9A8WLKpngOqAS6VJYnoZrZIjO7z8zWmNk5DWpONLO7zewuM/vnVL85AFGL5QR5AN1lJi3gMaPNnh327d9xxx1Bdccdd1xQ3dq1axNrPvGJTwS1tXjx4qC60BMZJyYmguo2bNgQVJfmH/bGjRuD6vbcc8+guv7+/k66k5uswsHMeiQtk3SkpBFJq8xsubvfXVOzUNK5kt7k7uvN7BWZdAbThHxm169fH9TWWWedFVT361//Oqju+OOPD6o788wzg+o2bdqUWBOaTaHcw9ZHOOSQQ4Lqtt5666C6p556KqguNO+KlmE+LZL0ZUk9kgbd/XN1ak6UdIEkl/Qbdw/7DxEdCf2dh37GHn/88cSa0LHOI488ElSX9uIOoZ//3t7eoLrDDjsssaavry+orR/84AdBdX/4Q9gRtFdddVVQ3amnnhpUl5WYsqmrJ2NADDLcUnOwpDXu/mD1fa6RdKyku2tq3itpmbuvlyR3fyKrzgCITxb5xIYiAJ2KKZs4TBEouXZ3tddes6Z6WzKl6Z0l1e6mHak+VmtvSXub2S1mdlt1ixAASMrsMMUtG4rcfVzS5g1FtdhQBKChjA5TzCSb2DMGlFy7W3dqr1nTqOl6L5vy9WxJCyUdqspFTG82s9e6+7NtdQrAjJLRnvt6G4qmHi+6d/X9b1HlcKEL3P3HWXQGQHzazabqhuvajdcD1fGUlFE2MRkDSi7DwxRHJO1a8/Uukh6tU3Nb9To5D5nZfapMzlZl1SkA8WgnnxIGOxIbigB0KKMN2ZlkU9PJmJn1SjpJ0qPufr2ZLZb03yXdo0p4hq3iAKBtGU7GVklaaGZ7SHpElc/61JNMfyDpZElXmtmOqmzxeTCrDoUim4ByaCefAvbaR72hiHwCipfR2CmTbEraM3ZFtWaumZ0iaZ6kf5H0VlWOmzylle8AQOuymoy5+0Yz+4Ckn6iyK/1yd7/LzD4ladjdl1efe5uZ3S1pk6S/c/enM+lQa8gmoAQyyqdoNxRVkU9AwWLKpqTJ2Ovc/b+Z2ezqm+7k7pvM7J8k/aaNbwJAi9JegreWu18r6dopj51fc98lfaR6KxOyCSiBLPIp8g1FEvkEFC6mbErq6azq7vZtJM2VNL/6eJ+krRq9qHYVt8HBwZDvD0ADGa0IFLu2skkin4A0ZZVN7n6tu+/t7nu5+0XVx86vDnbkFR9x933d/XXufk2G32arOh47DQ0N5dBNYObKauyURTYl7RkbknSvKrO/8yR9x8welPQGVZZzbNTRLceDj4+Ph131D0BdXTCxakdb2SRNzqexsTHyCegA+VRXx2On0dFRsgnoQEzZ1HQy5u5fMrNvVe8/amZXSTpC0tfd/Vd5dBDodjEFSl7IJqAcyKfpyCegeDFlU+LS9u7+aM39ZyV9N9MeAZgkpkDJE9kEFI98qo98AooVUzZxnTGg5GIKFADdhXwCUEYxZVN2y7QBAAAAABqakXvGenp6guqefPLJoLqlS5cG1T322GNBdSeccEJizXHHHRfUVmXl8WQTExNBdWkL/V3ccMMNiTW//e1vg9r67Gc/G1S3zTbbBNUV9bPbLKatO0hPyO/91ltvDWrr5ptvDqp785vfHFS3bNmyoLqtt946qK6Iz1jo5+quu+4KqnvhhReC6t71rncF1c2bNy+xJjT/s0Q+zRyhS4Fv3LgxqO7jH/94UN2vfpV8Gt26deuC2grNkr6+vqC6973vfUF1p5wSdum40M/szjvvnFhz0003BbX1/e9/P6huw4awa5GfdtppQXVFiymbZuRkDJhJYgoUAN2FfAJQRjFlE5MxoORiChQA3YV8AlBGMWUTkzGg5GIKFADdhXwCUEYxZROTMaDkQo/jB4C8kU8AyiimbEqcjJnZXpL+UtKukjZKekDS1e7+XMZ9A6C4tu7kiWwCikc+1Uc+AcWKKZuaThvN7IOSviqpX9JBkuaoEiz/bmaHZt47ADKztm4zGdkElAPZNB35BBQvprFT0j6890pa5O6flnSEpH3d/TxJiyR9qdGLzGyJmQ2b2fDg4GB6vQW6UEyBkqO2skkin4A0kU11dTx2GhoayqmrwMwU09gp5Jyx2ZI2SeqTtI0kufvDZrZVoxe4+4CkAUkaHx8v/kIoQMS6ZPDSjpazqVqzJZ/GxsbIJ6AD5FNDHY2dRkdHySagAzFlU9JkbFDSKjO7TdKbJV0sSWb2cknPZNw3AIorUHJENgElQD7VRT4BBYspm5pOxtz9y2Z2vaR9JF3i7vdWH39SlYABkLGYAiUvZBNQDuTTdOQTULyYsinxMEV3v0vSXTn0BUAdMQVKnsgmoHjkU33kE1CsmLKJ64wBJRdToADoLuQTgDKKKZuimoyF/mDXrVsXVHfmmWcG1a1YsSKo7sgjjwyqu+SSS4LqFixYkFjjXsw5vqG/i40bNwbVXX/99Yk1Bx98cFBbixcvDqoLvSDgxMREUF1WYgoUJAv9u3vxxRcTa5YuXRrUVmhOPPXUU0F1V199dVDdEUccEVS3xx57JNaknXWhv4eRkZFU3/ewww4Lquvr60usKTqbJPIpBj09PUF1oWOniy66KKjuO9/5TlDdhg0bEmte/epXB7V1wAEHBNWFjhMOPPDAoLo5c+YE1YW69dZbE2vOOuusoLbGxsaC6t74xjcG1b3pTW8KqitqfLpZTNkU1WSsm4RMxNAdYrqKPIDuQj4BKKOYsonJGFByMW3dAdBdyCcAZRRTNsUzbQQAAACAGYQ9Y0DJxbR1B0B3IZ8AlFFM2cRkDCi5mAIFQHchnwCUUUzZxGQMKLmYTkIF0F3IJwBlFFM2MRkDSi6mrTsAugv5BKCMYsqmptNGM5tvZp8zs3vN7Onq7Z7qY9s1ed0SMxs2s+HBwcH0ew10ETNr6zbTkU9A8cim6dLIpqGhoTy7DMw4MY2dkvaMfVvSDZIOdffHJMnMXinpFEnfkVT3KsfuPiBpQJLGx8eLveobELluGLy0qeN8GhsbI5+ADpBPdXWcTaOjo2QT0IGYsinpgMrd3f3izWEiSe7+mLtfLGm3bLsGQIpr607OyCegYGRTXWQTULCYxk5Je8Z+Z2ZnS/qGuz8uSWa2QNKpktZm3DcAiusk1JyRT0DByKe6yCagYDFlU1JP3yVpB0k3mdkzZvaMpJWSXibphIz7BkBxbd3JGfkEFIxsqotsAgoW09ip6Z4xd18v6WPV2yRmdpqkKzLqF4CqLhm8tIx8AopHPk1HNgHFiymbOlnafqlSDJSQH9qGDRuC2lq6dGlQ3Y9+9KOgut7e3qC6devWBdXdcccdiTWLFi0Kasu9mHN8Q3f/3nDDDUF1N954Y2LNhRdeGNTWTjvtFFQX+vdUtJgCpURSzacQoZ+J0dHRoLpvfvObiTUjIyNBbfX09ATV3X333UF1l19+eVDdO97xjqC6kL/x0KwL/V4ff/zxoLqbb745qO74448Pqnv7298eVBcL8qlluWdT6GfiF7/4RVDdwMBAUN2OO+4YVLdx48bU3nP//fcPqtu0aVNQXWjuPPfcc0F1IWMdSTr77LMTa9auDTvidZdddgmqu+yyy4LqFi5cGFQ3Pj4eVJeVmLKp6WTMzFY3ekrSgvS7A2CqmI57zhP5BBSPfJqObAKKF1M2Je0ZWyDpKEnrpzxukm7NpEcAJolp607OyCegYORTXWQTULCYsilpMrZC0jx3n3ZcnZmtzKRHACaJKVByRj4BBSOf6iKbgILFlE1JC3ic3uS5xel3B8BUMQVKnsgnoHjk03RkE1C8mLKpkwU8AOQgpkAB0F3IJwBlFFM2xXN2GwAAAAAUxMwWmdl9ZrbGzM5pUvdXZuZmdmBSm5lMxsxsiZkNm9nw4OBgFm8BdI1Zs2a1dQuRRaiUHfkEpIdsSk9tNg0NDRXdHSBqWYydzKxH0jJJR0vaV9LJZrZvnbptJH1Q0i9D+tr2YYpmdp27H13vOXcfkDQgSePj48VcCAuYIbLa1V4TKkdKGpG0ysyWu/vdU+paCpUyCM2nsbEx8gnoQBb5RDZJo6OjZBPQgYzGTgdLWuPuD1bf4xpJx0qaepHOCyV9XtJHQxpNus5Yo6vnmaTXh7wBgM7mHa8xAAAfCklEQVRkeNxzJqGSF/IJKF5MA568kE1A8TLKpp0l1V5te0TSIVPedz9Ju7r7CjPrfDImaZWkm1QJkKm2C3kDAJ1pN1DMbImkJTUPDVS3vG6WSajkiHwCChbTgCdHZBNQsIzGTvUa3bIX28xmSfqSpFNbec+kydg9kv6Xuz9Qp7Nr69QDSFm7V5GvPeSlgUxCJUfkE1CwdvIpYEMR2QSgIxmNnUYk7Vrz9S6SHq35ehtJr5W0sjoZfKWk5WZ2jLsPN3rPpMnYBWq8yMeZCa8FkIIMD1PMJFRydIHIJ6BQ7eRTwIYisglARzIaO62StNDM9pD0iKSTJG25dqC7Pydpx5o+rJT00aRcSrro83ebPL19cp8BdCrDyVgmoZIX8gkoXkwDnryQTUDxssgmd99oZh+Q9BNJPZIud/e7zOxTkobdfXk77XZy0eelkq7o4PWThPzQXnjhhaC2vve97wXVLViwIKjuoIMOCqp76KGHgupe/epXB9UVoaenJ6ju6aefDqq76KKLguoOPfTQxJp3vvOdQW1t3LgxqC4WWU3GsgqVkkgtn0J//qF/d1dcEdatSy65JLFmzpw5QW29/OUvD6rbaaedgupC+iZJr3rVq4LqJiYmEmtCDzlxD1uE7sorrwyqu/3224PqrrvuuqC6+fPnB9Vt2rQpqK5oMQ14SiLVsVOaQj9jW221VVDd+Ph4UN1JJ52UWLNw4cKgtv7whz8E1f3iF78IqgvNidAx0fBw2PaC3//+94k1p556alBbZ599dlBdaP5v2LAhqK5oGY6drpV07ZTHzm9Qe2hIm0mrKa5u9JSksJkMgI5keRX5LEIlL+QTULyYBjx5IZuA4mU5dkpb0p6xBZKOkrR+yuMm6dZMegRgknZPQu0C5BNQMPKpLrIJKFhM2ZQ0GVshaZ673zH1ieox2gAyFtPWnZyRT0DByKe6yCagYDFlU9ICHqc3eW5xo+cApCemQMkT+QQUj3yajmwCihdTNnWygAeAHMQUKAC6C/kEoIxiyiYmY0DJxXTcM4DuQj4BKKOYsqlpT81sWzP7rJl908wWT3nusiavW2Jmw2Y2PDg4mFZfga5kZm3dZjryCSge2TRdGtk0NDSUfUeBGSymsVPSnrErJD0g6XuS3mNm75S02N3HJL2h0YvcfUDSgCSNj4+HXfQFAFrTcT6NjY2RTwDS1nE2jY6Okk1Al0iajO3l7puvtPsDMztP0g1mdkzG/QJQ1Q1bkttEPgEFI5/qIpuAgsWUTUmTsT4zm+XuE5Lk7heZ2Yikf5M0L/PeAUBj5BOAMiKbAARLOrvtXyUdXvuAu39D0lmSxrPqFIA/ium455yRT0DByKa6yCagYDGNnZKuM3Z2g8d/bGafyaZLAGp1yeClZeQTUDzyaTqyCSheTNnUydL2S1U5STUV7snnqm6zzTZBbR1//PFBdd/61reC6vbZZ5+guo985CNBdXvttVdiTegfUU9PT1BdaHvr168PqvvoRz8aVLdw4cKguqVLlybW9PX1BbUV8rcUk5gCpURSy6fQ5XGffPLJoLoPf/jDQXUheXfyyScHtfXud787qK63tzeobs899wyqC82nEC+99FJQ3VVXXRVU9/nPfz6o7qKLLgqqe+1rXxtUt2nTpqC6WJBPLUt17BQi9P/EnXbaKahu/vz5QXXj42E7Ae+5557EmnPOOSeorYmJiaC6n/3sZ0F169atC6oL/Vy/7W1vC6rbb7/9Ems+9KEPBbW17bbbBtWFfg+xjLFiyqamkzEzW93oKUkL0u8OgKliCpQ8kU9A8cin6cgmoHgxZVPSnrEFko6SNHV3iUm6NZMeAZgkpkDJGfkEFIx8qotsAgoWUzYlTcZWSJrn7ndMfcLMVmbSIwCTxBQoOSOfgIKRT3WRTUDBYsqmpAU8Tm/y3OJGzwFA1sgnAGVENgFoRScLeADIQUxbdwB0F/IJQBnFlE1MxoCSiylQAHQX8glAGcWUTWFrNrfIzJaY2bCZDQ8ODmbxFkDXiOnChTEgn4D0kE3pqc2moaGhorsDRC2msVPS0vavlPRJSROSzpd0pqR3SrpH0ofcve4FGNx9QNKAJI2Pj8dxQQKgpBi81JdGPo2NjZFPQAfIp+nSyKbR0VGyCehATNmUtGfsSkl3S1or6UZJf5D055JulvTVTHsGQFJcW3dydqXIJ6BQZFNdV4psAgoV09gp8Tpj7v4VSTKz97v7xdXHv2JmDVcLApCeLhm8tIN8AgpGPtVFNgEFiymbkiZjtXvOrpryXE/KfQFQR0yBkjPyCSgY+VQX2QQULKZsSjpM8YdmNk+S3P0Tmx80sz+VdF+WHQOABOQTgDIimwAES7ro8/kNHl9jZj/KpksAasW0dSdP5BNQPPJpOrIJKF5M2WTu7S3YY2YPu/tuSXVprqYY+oN99tlng+pWrlwZVHfGGWcE1e2+++5Bdfvss09iTW9vb1Bbaf+xPfDAA0F1s2eHXaJuxYoVQXVbb711Ys2mTZuC2iq73t7eln5pt99+e1ufof333z+eJEpZaD6FrKYY+hnbsGFDUN2FF14YVPe1r30tsWbbbbcNautP/uRPgurGxsaC6kKz7nWve11QXYg777wzqO6HP/xhUN1pp50WVPeZz3wmqG7u3LlBde3+n5uXvr6+zPOJbErOpjRXUwzNsImJiaC6yy67LKjus5/9bGrvG5qvs2aFXbGppyfsaNHtt98+qG7x4sVBdWeeeWZQXUi2h/5ey545ofr7+2fs2ClpafvVjZ6StCD97gCYKqatO3kin4DikU/TkU1A8WLKpsTVFCUdJWn9lMdN0q2Z9AjAJDEFSs7IJ6Bg5FNdZBNQsJiyKWkytkLSPHe/Y+oTZrYykx4BmCSmQMkZ+QQUjHyqi2wCChZTNiUt4NHwehjuHnaALICOxBQoeSKfgOKRT9ORTUDxYsqmsFUYABQmpkAB0F3IJwBlFFM2tTwZM7NXuPsTWXQGwHQxBUrRyCcgX+RTGLIJyFdM2dR0DVAze9mU2w6SfmVm25vZy5q8bomZDZvZ8ODgYOqdBrqJmbV1m+nIJ6B4ZNN0aWTT0NBQjj0GZp6Yxk5Je8aekvS7KY/tLOl2SS5pz3ovcvcBSQNSutcZA4AaHedTyHXGAKBFHWdTmtcZA1BuSZOxsyUdIenv3P0/JcnMHnL3PTLvGQBJce1qzxn5BBSMfKqLbAIKFlM2Ja2m+A9mdo2kL5nZWkmfVGWrDoCcxBQoeSKfgOKRT9ORTUDxYsqmxAU83H1E0glm9heSfiZpbua9ArBFTIGSN/IJKBb5VB/ZBBQrpmxquoBHLXf/V0mHqbLrXWZ2WladAvBHMZ2EWhTyCSgG2dQc2QQUI6axk7m3t+fczB52992S6opYwGPWrLA55sTERFDd9ddfH1T34IMPBtXdcMMNqdVt2LAhqK1Qxx13XFDdBRdcEFS35551z1OeZtOmTUF1M0Fvb29Ln/b777+/rc/Q3nvv3V2jnhqh+ZTmAh6hufPcc88F1V100UWJNWvXrg1qKzTnQ/8jKuI/rNCM2GeffYLqPvaxjwXVzZkzJ6gu9P+Tsuvr68s8n8im5GwqYgGP0Ax75plngurWrFnTSXcmSTvDQvNkxx13DKpbuHBhUF2okO+33fF7rPr7+0sxdjKzRZK+LKlH0qC7f27K8x+RdIakjZKelPQed5+6oM8kTQ9TNLPVjZ6StKDZa9GZ0AkbZr5u25IcinwCikc+TUc2AcXLIpvMrEfSMklHShqRtMrMlrv73TVl/yHpQHd/ycz+RtLnJb2rWbtJ54wtkHSUpPVT+yPp1hb6D6CEstjCkyPyCZihyCYAJXSwpDXu/qAkVRfqOVbSlsmYu99YU3+bpL9OajRpMrZC0jx3v2PqE2a2MrnPADqV1ZbnrLbw5Ih8AgoW09bnHJFNQMEyGjvtLKn2/IARSYc0qT9d0nVJjSYtbX96k+cWJzUOoNQy2cKTF/IJmLHIJgCFMLMlkpbUPDRQvSC7VNm7PVXdc9PM7K8lHSjpLUnvmbi0PYBitbt1JyFQpIy28ADoHu3kE9kEIGvtjp2qWTTQ4OkRSbvWfL2LpEfrvPcRks6T9BZ3H0t6TyZjwAyVEChSRlt4AKAZsglApFZJWmhme0h6RNJJkibt7Taz/SR9TdIid38ipFEmY0DJZbhaWSZbeAB0j4zyiWwC0JEsssndN5rZByT9RJXFhS5397vM7FOSht19uaQvSJon6TvVPjzs7sc0a7fpRSWqqxltvj/fzIbMbLWZ/bOZNVye1cyWmNmwmQ0PDg4Gf5MApsvwwoVbtvCYWa8qW3iWT3nvzVt4jgndwpMX8gkoHtk0XRrZNDQ0lE9ngRkqq7GTu1/r7nu7+17uflH1sfOrEzG5+xHuvsDdX1+9NZ2IScl7xj4j6cfV+1+UtE7SX0g6XpUQrHuF4NpDEIq46DMwk2S1ZyyrLTw56jif0rzoM9CNYtr6nKOOs6mIiz4DM0mGRxWlrpXDFA9099dX73/JzE7JokMA8uPu10q6dspj59fcPyL3TrWHfAJmELIJQLdImoy9onphRZO0rZmZu2/eWtP0EEcA6Yhp607OyCegYORTXWQTULCYsilpMvZ1SdtU739D0o6SnjSzV0qadjFDAOmLKVByRj4BBSOf6iKbgILFlE1JF31e2uDxx8zsxnrPAUAeyCcAZUQ2AWiF/XHPeYsvNHvY3XdLqpsJC3j09PQE1YXOwn//+98n1tx3331BbY2Nha3mG/o9vOY1rwmq22abbZKLJG3atCmorpv09va2tLlm7dq1bX2Gdt1113g2C6UsNJ+KWMAjNCcmJiYy7snMFPrzDa1r9//IWPX19WWeT2RTcjaVeQGP0M/OrFnxH5EZ+vknr7PX398/Y8dOTfeMmdnqRk9Jarg8K4D0xLSrPU/kE1A88mk6sgkoXkzZlHTO2AJJR0laP+Vxk3RrJj0CMElMgZIz8gkoGPlUF9kEFCymbEqajK2QNM/dp51wamYrM+kRAIQhnwCUEdkEIFjSAh6nN3lucfrdATBVTFt38kQ+AcUjn6Yjm4DixZRNrVz0GUABYgoUAN2FfAJQRjFlE5MxoORiChQA3YV8AlBGMWVTy+uOmtkOATVLzGzYzIYHBwfb6xkAtIh8AlBGrWbT0NBQHt0CUAJJS9t/TtI/uPtTZnagpG9LmjCzrSS9291vqvc6dx+QNCDNjOuMASifNPKpiOuMAZjZ0simMl9nDEC6kvaM/bm7P1W9/wVJ73L3P5V0pKQvZtozAJIqu9rbuXUB8gkoGNlUF9kEFCymsVPSOWNbmdlsd98oaY67r5Ikd7/fzPqy7x6ALhm8tIN8AgpGPtVFNgEFiymbkvaMLZN0rZkdLunHZvaPZvZmM1sqadr1MwAgR+QTgDIimwAES7rO2FfM7D8l/Y2kvav1e0v6gaQLs+8egJi27uSJfAKKRz5NRzYBxYspmxKXtnf3lZJWTn3czE6TdEX6XQJQK6ZAyRv5BBSLfKqPbAKKFVM2mXt7C/aY2cPuvltSHaspThfyB9LT05NDT6bbtGlTUF27fzeQent7W0qIp59+uq0f9g477BBPEqUsNJ/KvJpiTP+RxIgMq6+vry/zfCKbkrOJ1RSByfr7+2fs2ClpafvVjZ6StCD97gBAGPIJQBmRTQBakXSY4gJJR0laP+Vxk3RrJj0CMAl7SBoin4CCkU91kU1AwWLKpqTJ2ApJ89x92uo/ZrYykx4BmCSmQMkZ+QQUjHyqi2wCChZTNrV9zlgozhmbjnPGulur54ytX7++rR/29ttvH08SFYRzxroXGVZfq+eMtZNPZFMyzhkDJmv1nLGYxk6JqykCKBaDcgBlRT4BKKOYsonJGFByMQUKgO5CPgEoo5iyaVazJ83sdjP7hJnt1UqjZrbEzIbNbHhwcLCzHgJAHeQTgDJKI5uGhoay6h6AkknaM7a9pO0k3Whmj0m6WtK33P3RZi9y9wFJAxLnjAGdimnrTs46zqcynzMGxIB8qqvjbOKcMaAzMWVT0z1jkta7+0erFyg8S9JCSbeb2Y1mtiT77gFAQ+QTgDIimwAES5qMbeHuN7v7+yXtLOliSW/MrFcA0ALyCUAZkU0AkiQdpnj/1AfcfZOkH1dvADIW0672nJFPQMHIp7rIJqBgMWVT0z1j7n5So+fM7LT0uwMAYcgnAGVENgFoRdsXfTazh6vHQzfFAh7AZK1e9PnFF19s6zM0b968eDYLpSw0n1jAA5is1Ys+t5NPZFNyNrGABzBZqxd9jmns1PQwRTNb3egpSQvS7w4AhCGfAJQR2QSgFUnnjC2QdJSk9VMeN0m3ZtIjAJPEdNxzzsgnoGDkU11kE1CwmLIpaTK2QtI8d79j6hNmtjKTHgFAGPIJQBmRTQCCtX3OWCjOGQMma/WcsZdeeqmtz9DcuXPj2SxUEM4ZAyZr9ZyxdvKJbErGOWPAZK2eMxbT2Cn4OmMAAAAAgPQkHaYIoGAxHfcMoLuQTwDKKKZsymTPmJktMbNhMxseHBzM4i0AoC3kE4Ayqs2moaGhorsDoA4zW2Rm95nZGjM7p87zfWb2rerzvzSz3RPbbHbOmJkdKOkLkh6RdK6kyyUdrMrV5Ze4+38kvQHnjAGTtXrOWLvnDrR6fHVs0sgnzhkDJmv1nLF28olsSs4mzhkDJms1N7IYO5lZjyqf4yMljUhaJelkd7+7pub9kv6bu7/PzE6S9Jfu/q5m75m0Z+wySZ+X9CNVlmP9mrvPl3RO9TkAKAr5BMxQWWx9zhHZBMxMB0ta4+4Puvu4pGskHTul5lhJ36je/66kt1rCMZNJk7Gt3P06d79akrv7d1W583NJ/a1+BwDKJfIBD/kEzEDVrc/LJB0taV9JJ5vZvlPKTpe03t3/VNKXJF2cby+bIpuASNUeLly9Lal5emdJa2u+Hqk+pno17r5R0nOSdmj2nkkLeIya2dskzZfkZnacu//AzN4iaVPytwSgU1mdhFoz4Nmyu93MltfublfNgKe6u/1iSU13t+eIfAIKllE+bdn6XH2PzVufa7PpWEkXVO9/V9KlZmae9fV6wpBNQMHazSZ3H5A00KjZei9po2aSpD1j75N0lqT3qHI1+cPM7FlVdrN/MOG1AMotk93tOSKfgAglbHmWMtr6nCOyCZiZRiTtWvP1LpIebVRjZrNV2SjzTLNGm+4Zc/ffqBIkm32oepOZnabKsdAAMpTh3KfegOeQRjXuvtHMNg94nsqqU6HIJ6B47eRTwpZnKaOtz3khm4DiZTR2WiVpoZntocoCPSdJWjylZrmkUyT9u6S/knRD0h77Tpa2X9rBawFkLGDrc9QDngTkExCvTLY+lwTZBESquhf+A5J+IukeSd9297vM7FNmdky1bEjSDma2RtJHVFm4p6mme8bMbHWjpyQtCO08gPwFbH1uZcAzUrYBD/kEzFiZbH3OC9kEzFzufq2ka6c8dn7N/VFJJ7TSZtICHgtU2dW+fsrjJnazA7nI8DDFqAc8Ip+AwmWRT9VDojdvfe6RdPnmrc+Sht19uSpbn79Z3fr8jCr5VRZkE1Cw8pzenixpMrZC0jx3v2PqE2a2MpMeAcjFDBjwkE/ADJXF1ucckU0AglnWG7nHx8fLshUdKIXe3t6WNtds3Lixrc/Q7Nmz49ksVJCxsTHyCajR19eXeT6RTclGR0fJJqBGf3//jB07dbKABwAAAACgTUmHKQIoWEzHPQPoLuQTgDKKKZua7hkzs3nV5RrvMrPnzOxJM7vNzE5NeN2WJbUHBwdT7TAASOQTgHJKI5uGhoZy6i2AojU9Z8zMfijp+5Kul3SipK0lXSPpE5IecfePJ70B54wBk7V6ztjExERbn6FZs2bFs1moDWnkE+eMAZO1es5YO/lENiVnE+eMAZO1es5YTGOnpMnYb9z9/635epW7H2RmsyTd7e6vSXoDJmPAZEzG0pFGPjEZAyZjMta5NLKJyRgw2UyejCUt4PF7M/szSTKzv1D1Yq/uPqHK9TIAZMzM2rp1AfIJKBjZVBfZBBQsprFT0gIe75M0aGZ7S7pT0nskycxeLmlZxn0DgGbIJwBlRDYBCNZ0MubuqyUdXOfxJ83shcx6BWCLLtmS3DLyCSge+TQd2QQUL6Zsavuiz2b2sLvvllTHOWPAZK2eMyap3c9QPEmUstB84pwxYLJWzxlTe/lENiXgnDFgslbPGVNMYyd3b3iTtLrB7T8ljTV7bUK7S9p9bWztlblvfK/laItb27+D0udTmf+G+V7L0VYM7XFr+edf+mxKu70y943vtTztkU31b0mrKT4u6ShJ66c+JelWd9+p4YubMLNhdz+wndfG1l6Z+5Z2e2XuW9rtpd03tC6GfCrz33Da7ZW5b2m3V+a+ZdEeWhNDNqXdXpn7lnZ7Ze5b2dsjm+pLWsBjhaR57n7H1CfMbGUmPQKAMOQTgDIimwAES1rA4/Qmzy1OvzsAEIZ8AlBGZBOAViRdZywrA13UXpn7lnZ7Ze5b2u2l3TeUR5n/TsrcXpn7lnZ7Ze5bFu2hHMr8d1LmvqXdXpn7Vvb2yKY62l5NEQAAAADQvqL2jAEAAABAV2MyBgAAAAAFyH0yZmaLzOw+M1tjZud00M6uZnajmd1jZneZ2YdS6l+Pmf2Hma1Ioa3tzOy7ZnZvtZ9v7KCtv61+n3ea2dVm1t/i6y83syfM7M6ax15mZj8zsweq/27fYXtfqH6vq83s+2a2XSft1Tz3UTNzM9uxk7bM7Mzq395dZvb5TvpmZq83s9vM7A4zGzazg0PbQzmllU3VtlLPp7JmU7W90uRTmbOpWXvkExopezZV2y1lPpUpm5q0V4p8IpsKlOdFzST1SPovSXtK6pX0G0n7ttnWqyTtX72/jaT7221rSrsfkfTPklak0NY3JJ1Rvd8rabs229lZ0kOS5lS//rakU1ts482S9pd0Z81jn5d0TvX+OZIu7rC9t0maXb1/caftVR/fVdJPJP1O0o4d9O0wSddL6qt+/YoOv9efSjq6ev/tklZ2+vfCrbhbmtlUbS/1fCpjNlVfX6p8KnM2Nekf+cSt0e+39NlUbat0+VS2bGrSXinyiWwq7pb3nrGDJa1x9wfdfVzSNZKObachd1/n7rdX778g6R5VPnhtM7NdJP25pMFO2qm2ta0qf4hD1T6Ou/uzHTQ5W9IcM5staa6kR1t5sbv/m6Rnpjx8rCqhp+q/x3XSnrv/1N03Vr+8TdIuHfZPkr4k6WxJwSvNNGjrbyR9zt3HqjVPdNieS9q2en++Wvx9oHRSyyYp/XwqeTZJJcqnMmdTk/bIJzRS6mySSp9PpcmmRu2VJZ/IpuLkPRnbWdLamq9H1GEISJKZ7S5pP0m/7LCpf1Tlj3eiw3akylasJyVdUd11P2hmW7fTkLs/IukfJD0saZ2k59z9pyn0cYG7r6u+xzpJr0ihzc3eI+m6Thows2MkPeLuv0mhP3tL+h9m9kszu8nMDuqwvQ9L+oKZrVXld3Nuxz1EkTLJJim1fCplNklR5lPZskkin9BY2bNJKmk+RZhNUvnyiWzKQd6TMavzWEdr65vZPEnfk/Rhd3++g3beIekJd/91J/2pMVuV3bP/x933k/R7VXZnt9O37VXZErOHpJ0kbW1mf51SP1NnZudJ2ijp/3bQxlxJ50k6P6VuzZa0vaQ3SPo7Sd82s3p/j6H+RtLfuvuukv5W1a14iFbq2SSlk09lzqZq/6LJp5Jmk0Q+obHSZlO1ndLmU0zZJJU2n8imHOQ9GRtR5TjWzXZRB7sozWwrVcLk/7r7v3TYtzdJOsbMfqvKYQCHm9k/ddDeiKQRd9+8xem7qgRMO46Q9JC7P+nuGyT9i6T/3kHfNnvczF4lSdV/g3c/N2Jmp0h6h6T/6e6d/IexlyoB+pvq72QXSbeb2SvbbG9E0r94xa9U2YIXfNJ9Haeo8nuQpO+ocigJ4pVqNkmp5lOZs0mKJJ9KnE0S+YTGypxNUrnzKYpsqrZT1nwim3KQ92RslaSFZraHmfVKOknS8nYaqs7MhyTd4+6XdNoxdz/X3Xdx992r/brB3dveguLuj0laa2avrj70Vkl3t9ncw5LeYGZzq9/3W1U5zrtTy1X5YKj67w87aczMFkn6mKRj3P2lTtpy9/9091e4++7V38mIKicdP9Zmkz+QdHi1n3urclLwUx108VFJb6neP1zSAx20heKllk1SuvlU8mySIsinkmeTRD6hsdJmk1T6fCp9NkmlzyeyKQ+e84ohqqyecr8qqwOd10E7f6bKrvrVku6o3t6eUh8PVTorAr1e0nC1jz+QtH0HbS2VdK+kOyV9U9WVbVp4/dWqHDO9QZUP5+mSdpD0c1U+DD+X9LIO21ujyrHtm38fX+2kvSnP/1bhKwLV61uvpH+q/vxul3R4h9/rn0n6tSorW/1S0gFp/O1xK+6WVjZV28okn8qYTdX2SpNPZc6mJv0jn7g1+x2XPpuqbZcun8qUTU3aK0U+kU3F3az6AwMAAAAA5Cj3iz4DAAAAAJiMAQAAAEAhmIwBAAAAQAGYjAEAAABAAZiMAQAAAEABmIwBAAAAQAGYjAEAAABAAf5/il4TalLH7l4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# look at some examples\n",
    "fig = plt.figure(figsize=(15, 5))\n",
    "axes = fig.subplots(1,3)\n",
    "for i in range(3):\n",
    "    _ = sns.heatmap(X[i,:].reshape(20,20), cmap='Greys', ax=axes[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 524,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "<a id=\"section1\"></a>\n",
    "### 1.3 Feedforward and cost function\n",
    "\n",
    "Now you will implement the cost function and gradient for the neural network. First, complete the code for the function `nnCostFunction` in the next cell to return the cost.\n",
    "\n",
    "Recall that the cost function for the neural network (without regularization) is:\n",
    "\n",
    "$$ J(\\theta) = \\frac{1}{m} \\sum_{i=1}^{m}\\sum_{k=1}^{K} \\left[ - y_k^{(i)} \\log \\left( \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) - \\left( 1 - y_k^{(i)} \\right) \\log \\left( 1 - \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) \\right]$$\n",
    "\n",
    "where $h_\\theta \\left( x^{(i)} \\right)$ is computed as shown in the neural network figure above, and K = 10 is the total number of possible labels. Note that $h_\\theta(x^{(i)})_k = a_k^{(3)}$ is the activation (output\n",
    "value) of the $k^{th}$ output unit. Also, recall that whereas the original labels (in the variable y) were 0, 1, ..., 9, for the purpose of training a neural network, we need to encode the labels as vectors containing only values 0 or 1, so that\n",
    "\n",
    "$$ y = \n",
    "\\begin{bmatrix} 1 \\\\ 0 \\\\ 0 \\\\\\vdots \\\\ 0 \\end{bmatrix}, \\quad\n",
    "\\begin{bmatrix} 0 \\\\ 1 \\\\ 0 \\\\ \\vdots \\\\ 0 \\end{bmatrix}, \\quad \\cdots  \\quad \\text{or} \\qquad\n",
    "\\begin{bmatrix} 0 \\\\ 0 \\\\ 0 \\\\ \\vdots \\\\ 1 \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "For example, if $x^{(i)}$ is an image of the digit 5, then the corresponding $y^{(i)}$ (that you should use with the cost function) should be a 10-dimensional vector with $y_5 = 1$, and the other elements equal to 0.\n",
    "\n",
    "You should implement the feedforward computation that computes $h_\\theta(x^{(i)})$ for every example $i$ and sum the cost over all examples. **Your code should also work for a dataset of any size, with any number of labels** (you can assume that there are always at least $K \\ge 3$ labels).\n",
    "\n",
    "<div class=\"alert alert-box alert-warning\">\n",
    "**Implementation Note:** The matrix $X$ contains the examples in rows (i.e., X[i,:] is the i-th training example $x^{(i)}$, expressed as a $n \\times 1$ vector.) When you complete the code in `nnCostFunction`, you will need to add the column of 1’s to the X matrix. The parameters for each unit in the neural network is represented in Theta1 and Theta2 as one row. Specifically, the first row of Theta1 corresponds to the first hidden unit in the second layer. You can use a for-loop over the examples to compute the cost.\n",
    "</div>\n",
    "<a id=\"nnCostFunction\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## nnCostFunction implementation notes:\n",
    "\n",
    "### 1) One-hot encoding of y:\n",
    "y_encode is also a 5000x10 matrix, each row is a encoded label for a training iage. If the label is 2, the 3rd element in the row is 1 and the rest 9 elements are all 0s.\n",
    "$$ y = \\begin{bmatrix} 0 & 0 & 1 & 0 \\cdots & 0 \\end{bmatrix}$$\n",
    "\n",
    "### 2) feedforwad outputs (activation of layer 3: a3): \n",
    "a3 is a 5000x10 matrix. Each row is a prediction for a training image. It shows the probability of being each number, 0 to 9. For instance, following y_hat means the predicted digit is 2 with 90% chance and is 0 with 10% probability.\n",
    "$$ y_{hat} = \\begin{bmatrix} 0.1 & 0 & 0.9 & 0 & \\cdots & 0 \\end{bmatrix}$$\n",
    "\n",
    "\n",
    "### 3) calculate the cost (loss),\n",
    "The cost function without regularization is:\n",
    "$$ J(\\theta) = \\frac{1}{m} \\sum_{i=1}^{m}\\sum_{k=1}^{K} \\left[ - y_k^{(i)} \\log \\left( \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) - \\left( 1 - y_k^{(i)} \\right) \\log \\left( 1 - \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) \\right]$$\n",
    "\n",
    "The cost function for neural networks with regularization is given by:\n",
    "$$ J(\\theta) = \\frac{1}{m} \\sum_{i=1}^{m}\\sum_{k=1}^{K} \\left[ - y_k^{(i)} \\log \\left( \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) - \\left( 1 - y_k^{(i)} \\right) \\log \\left( 1 - \\left( h_\\theta \\left( x^{(i)} \\right) \\right)_k \\right) \\right] + \\frac{\\lambda}{2 m} \\left[ \\sum_{j=1}^{25} \\sum_{k=1}^{400} \\left( \\Theta_{j,k}^{(1)} \\right)^2 + \\sum_{j=1}^{10} \\sum_{k=1}^{25} \\left( \\Theta_{j,k}^{(2)} \\right)^2 \\right] $$\n",
    "use the cost function formula? cost for every class and cost for every training image\n",
    "\n",
    "4) the dimensions do not look right, why does it still work?\n",
    "\n",
    "5) * is element_wise operation, not matrix multiplication @\n",
    "\n",
    "6) for the cost, we only need element-wise multiplication, what is the ituition. when we predicted the exact label at 100% confidence, the cost is 0. when we predict wrong label with 100% confidence, the cost is infinity.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## modify the implementation myself"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(z):\n",
    "    return 1/(1 + np.exp(-z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is simply applying sigmoid function to each element of matrix, and then for each element, times itself with (1-the element)\n",
    "def sigmoidGradient(z):\n",
    "    sig = sigmoid(z)    \n",
    "    return sig * (1-sig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### sigmoidGradient peaks when activation is 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g is [1.92874985e-22 2.13181162e-22 2.35625206e-22 2.60432193e-22\n",
      " 2.87850898e-22 3.18156288e-22 3.51652275e-22 3.88674772e-22\n",
      " 4.29595055e-22 4.74823489e-22 5.24813642e-22 5.80066837e-22\n",
      " 6.41137173e-22 7.08637090e-22 7.83243502e-22 8.65704593e-22\n",
      " 9.56847315e-22 1.05758568e-21 1.16892995e-21 1.29199670e-21\n",
      " 1.42802012e-21 1.57836430e-21 1.74453694e-21 1.92820449e-21\n",
      " 2.13120885e-21 2.35558583e-21 2.60358556e-21 2.87769508e-21\n",
      " 3.18066329e-21 3.51552845e-21 3.88564874e-21 4.29473587e-21\n",
      " 4.74689232e-21 5.24665251e-21 5.79902823e-21 6.40955891e-21\n",
      " 7.08436722e-21 7.83022039e-21 8.65459814e-21 9.56576766e-21\n",
      " 1.05728665e-20 1.16859943e-20 1.29163139e-20 1.42761635e-20\n",
      " 1.57791801e-20 1.74404367e-20 1.92765929e-20 2.13060625e-20\n",
      " 2.35491978e-20 2.60284939e-20 2.87688141e-20 3.17976395e-20\n",
      " 3.51453443e-20 3.88455007e-20 4.29352153e-20 4.74555014e-20\n",
      " 5.24516902e-20 5.79738855e-20 6.40774661e-20 7.08236411e-20\n",
      " 7.82800639e-20 8.65215105e-20 9.56306293e-20 1.05698770e-19\n",
      " 1.16826901e-19 1.29126618e-19 1.42721269e-19 1.57747185e-19\n",
      " 1.74355054e-19 1.92711424e-19 2.13000382e-19 2.35425393e-19\n",
      " 2.60211343e-19 2.87606797e-19 3.17886487e-19 3.51354070e-19\n",
      " 3.88345171e-19 4.29230753e-19 4.74420833e-19 5.24368594e-19\n",
      " 5.79574933e-19 6.40593481e-19 7.08036157e-19 7.82579302e-19\n",
      " 8.64970465e-19 9.56035897e-19 1.05668884e-18 1.16793868e-18\n",
      " 1.29090108e-18 1.42680914e-18 1.57702582e-18 1.74305755e-18\n",
      " 1.92656935e-18 2.12940156e-18 2.35358826e-18 2.60137768e-18\n",
      " 2.87525476e-18 3.17796604e-18 3.51254724e-18 3.88235366e-18\n",
      " 4.29109388e-18 4.74286690e-18 5.24220328e-18 5.79411058e-18\n",
      " 6.40412353e-18 7.07835959e-18 7.82358027e-18 8.64725894e-18\n",
      " 9.55765577e-18 1.05639006e-17 1.16760844e-17 1.29053607e-17\n",
      " 1.42640571e-17 1.57657992e-17 1.74256470e-17 1.92602461e-17\n",
      " 2.12879947e-17 2.35292278e-17 2.60064214e-17 2.87444178e-17\n",
      " 3.17706747e-17 3.51155406e-17 3.88125592e-17 4.28988057e-17\n",
      " 4.74152585e-17 5.24072105e-17 5.79247229e-17 6.40231276e-17\n",
      " 7.07635818e-17 7.82136814e-17 8.64481392e-17 9.55495334e-17\n",
      " 1.05609136e-16 1.16727830e-16 1.29017117e-16 1.42600239e-16\n",
      " 1.57613414e-16 1.74207199e-16 1.92548003e-16 2.12819755e-16\n",
      " 2.35225749e-16 2.59990681e-16 2.87362903e-16 3.17616915e-16\n",
      " 3.51056117e-16 3.88015850e-16 4.28866760e-16 4.74018518e-16\n",
      " 5.23923923e-16 5.79083446e-16 6.40050250e-16 7.07435733e-16\n",
      " 7.81915665e-16 8.64236959e-16 9.55225167e-16 1.05579275e-15\n",
      " 1.16694825e-15 1.28980638e-15 1.42559919e-15 1.57568849e-15\n",
      " 1.74157941e-15 1.92493560e-15 2.12759580e-15 2.35159239e-15\n",
      " 2.59917168e-15 2.87281651e-15 3.17527109e-15 3.50956856e-15\n",
      " 3.87906138e-15 4.28745498e-15 4.73884489e-15 5.23775783e-15\n",
      " 5.78919710e-15 6.39869275e-15 7.07235705e-15 7.81694578e-15\n",
      " 8.63992595e-15 9.54955076e-15 1.05549423e-14 1.16661830e-14\n",
      " 1.28944168e-14 1.42519610e-14 1.57524296e-14 1.74108698e-14\n",
      " 1.92439132e-14 2.12699422e-14 2.35092748e-14 2.59843677e-14\n",
      " 2.87200422e-14 3.17437328e-14 3.50857622e-14 3.87796457e-14\n",
      " 4.28624270e-14 4.73750498e-14 5.23627685e-14 5.78756020e-14\n",
      " 6.39688352e-14 7.07035734e-14 7.81473553e-14 8.63748301e-14\n",
      " 9.54685062e-14 1.05519579e-13 1.16628844e-13 1.28907709e-13\n",
      " 1.42479313e-13 1.57479756e-13 1.74059469e-13 1.92384720e-13\n",
      " 2.12639281e-13 2.35026275e-13 2.59770206e-13 2.87119216e-13\n",
      " 3.17347572e-13 3.50758417e-13 3.87686807e-13 4.28503076e-13\n",
      " 4.73616545e-13 5.23479629e-13 5.78592377e-13 6.39507480e-13\n",
      " 7.06835819e-13 7.81252591e-13 8.63504075e-13 9.54415124e-13\n",
      " 1.05489743e-12 1.16595867e-12 1.28871260e-12 1.42439027e-12\n",
      " 1.57435228e-12 1.74010253e-12 1.92330323e-12 2.12579158e-12\n",
      " 2.34959821e-12 2.59696755e-12 2.87038033e-12 3.17257842e-12\n",
      " 3.50659240e-12 3.87577189e-12 4.28381916e-12 4.73482629e-12\n",
      " 5.23331615e-12 5.78428779e-12 6.39326658e-12 7.06635961e-12\n",
      " 7.81031691e-12 8.63259919e-12 9.54145262e-12 1.05459916e-11\n",
      " 1.16562899e-11 1.28834822e-11 1.42398752e-11 1.57390713e-11\n",
      " 1.73961052e-11 1.92275941e-11 2.12519051e-11 2.34893386e-11\n",
      " 2.59623326e-11 2.86956873e-11 3.17168137e-11 3.50560091e-11\n",
      " 3.87467601e-11 4.28260791e-11 4.73348752e-11 5.23183642e-11\n",
      " 5.78265228e-11 6.39145888e-11 7.06436159e-11 7.80810854e-11\n",
      " 8.63015831e-11 9.53875477e-11 1.05430097e-10 1.16529941e-10\n",
      " 1.28798394e-10 1.42358488e-10 1.57346211e-10 1.73911864e-10\n",
      " 1.92221575e-10 2.12458961e-10 2.34826970e-10 2.59549917e-10\n",
      " 2.86875735e-10 3.17078458e-10 3.50460969e-10 3.87358044e-10\n",
      " 4.28139700e-10 4.73214911e-10 5.23035711e-10 5.78101722e-10\n",
      " 6.38965168e-10 7.06236413e-10 7.80590078e-10 8.62771812e-10\n",
      " 9.53605766e-10 1.05400286e-09 1.16496992e-09 1.28761976e-09\n",
      " 1.42318236e-09 1.57301721e-09 1.73862690e-09 1.92167224e-09\n",
      " 2.12398887e-09 2.34760571e-09 2.59476528e-09 2.86794619e-09\n",
      " 3.16988802e-09 3.50361874e-09 3.87248515e-09 4.28018640e-09\n",
      " 4.73081106e-09 5.22887818e-09 5.77938258e-09 6.38784493e-09\n",
      " 7.06036715e-09 7.80369354e-09 8.62527849e-09 9.53336117e-09\n",
      " 1.05370482e-08 1.16464050e-08 1.28725565e-08 1.42277992e-08\n",
      " 1.57257239e-08 1.73813525e-08 1.92112882e-08 2.12338823e-08\n",
      " 2.34694182e-08 2.59403149e-08 2.86713513e-08 3.16899155e-08\n",
      " 3.50262787e-08 3.87138994e-08 4.27897584e-08 4.72947301e-08\n",
      " 5.22739922e-08 5.77774785e-08 6.38603803e-08 7.05836993e-08\n",
      " 7.80148595e-08 8.62283835e-08 9.53066397e-08 1.05340669e-07\n",
      " 1.16431095e-07 1.28689138e-07 1.42237726e-07 1.57212730e-07\n",
      " 1.73764324e-07 1.92058495e-07 2.12278703e-07 2.34627723e-07\n",
      " 2.59329681e-07 2.86632297e-07 3.16809371e-07 3.50163529e-07\n",
      " 3.87029260e-07 4.27776266e-07 4.72813173e-07 5.22591625e-07\n",
      " 5.77610818e-07 6.38422504e-07 7.05636521e-07 7.79926912e-07\n",
      " 8.62038686e-07 9.52795283e-07 1.05310684e-06 1.16397930e-06\n",
      " 1.28652453e-06 1.42197145e-06 1.57167833e-06 1.73714649e-06\n",
      " 1.92003527e-06 2.12217870e-06 2.34560392e-06 2.59255146e-06\n",
      " 2.86549773e-06 3.16717987e-06 3.50062315e-06 3.86917133e-06\n",
      " 4.27652020e-06 4.72675463e-06 5.22438949e-06 5.77441497e-06\n",
      " 6.38234657e-06 7.05428044e-06 7.79695445e-06 8.61781576e-06\n",
      " 9.52509549e-06 1.05278912e-05 1.16362581e-05 1.28613099e-05\n",
      " 1.42153301e-05 1.57118951e-05 1.73660103e-05 1.91942606e-05\n",
      " 2.12149764e-05 2.34484173e-05 2.59169751e-05 2.86453981e-05\n",
      " 3.16610390e-05 3.49941290e-05 3.86780802e-05 4.27498203e-05\n",
      " 4.72501624e-05 5.22242132e-05 5.77218246e-05 6.37980924e-05\n",
      " 7.05139072e-05 7.79365635e-05 8.61404321e-05 9.52077034e-05\n",
      " 1.05229208e-04 1.16305325e-04 1.28546982e-04 1.42076761e-04\n",
      " 1.57030125e-04 1.73556760e-04 1.91822074e-04 2.12008835e-04\n",
      " 2.34318991e-04 2.58975677e-04 2.86225426e-04 3.16340612e-04\n",
      " 3.49622146e-04 3.86402452e-04 4.27048741e-04 4.71966637e-04\n",
      " 5.21604154e-04 5.76456101e-04 6.37068919e-04 7.04046017e-04\n",
      " 7.78053645e-04 8.59827359e-04 9.50179134e-04 1.05000518e-03\n",
      " 1.16029454e-03 1.28213852e-03 1.41674105e-03 1.56543006e-03\n",
      " 1.72966989e-03 1.91107492e-03 2.11142450e-03 2.33267915e-03\n",
      " 2.57699834e-03 2.84675980e-03 3.14458048e-03 3.47333932e-03\n",
      " 3.83620191e-03 4.23664706e-03 4.67849547e-03 5.16594036e-03\n",
      " 5.70358038e-03 6.29645443e-03 6.95007865e-03 7.67048521e-03\n",
      " 8.46426287e-03 9.33859892e-03 1.03013220e-02 1.13609455e-02\n",
      " 1.25267104e-02 1.38086267e-02 1.52175126e-02 1.67650289e-02\n",
      " 1.84637079e-02 2.03269736e-02 2.23691501e-02 2.46054561e-02\n",
      " 2.70519797e-02 2.97256309e-02 3.26440640e-02 3.58255659e-02\n",
      " 3.92889021e-02 4.30531137e-02 4.71372567e-02 5.15600757e-02\n",
      " 5.63396037e-02 6.14926796e-02 6.70343772e-02 7.29773399e-02\n",
      " 7.93310199e-02 8.61008233e-02 9.32871678e-02 1.00884468e-01\n",
      " 1.08880071e-01 1.17253169e-01 1.25973746e-01 1.35001593e-01\n",
      " 1.44285469e-01 1.53762489e-01 1.63357801e-01 1.72984657e-01\n",
      " 1.82544942e-01 1.91930240e-01 2.01023485e-01 2.09701205e-01\n",
      " 2.17836353e-01 2.25301655e-01 2.31973363e-01 2.37735251e-01\n",
      " 2.42482670e-01 2.46126420e-01 2.48596211e-01 2.49843502e-01\n",
      " 2.49843502e-01 2.48596211e-01 2.46126420e-01 2.42482670e-01\n",
      " 2.37735251e-01 2.31973363e-01 2.25301655e-01 2.17836353e-01\n",
      " 2.09701205e-01 2.01023485e-01 1.91930240e-01 1.82544942e-01\n",
      " 1.72984657e-01 1.63357801e-01 1.53762489e-01 1.44285469e-01\n",
      " 1.35001593e-01 1.25973746e-01 1.17253169e-01 1.08880071e-01\n",
      " 1.00884468e-01 9.32871678e-02 8.61008233e-02 7.93310199e-02\n",
      " 7.29773399e-02 6.70343772e-02 6.14926796e-02 5.63396037e-02\n",
      " 5.15600757e-02 4.71372567e-02 4.30531137e-02 3.92889021e-02\n",
      " 3.58255659e-02 3.26440640e-02 2.97256309e-02 2.70519797e-02\n",
      " 2.46054561e-02 2.23691501e-02 2.03269736e-02 1.84637079e-02\n",
      " 1.67650289e-02 1.52175126e-02 1.38086267e-02 1.25267104e-02\n",
      " 1.13609455e-02 1.03013220e-02 9.33859892e-03 8.46426287e-03\n",
      " 7.67048521e-03 6.95007865e-03 6.29645443e-03 5.70358038e-03\n",
      " 5.16594036e-03 4.67849547e-03 4.23664706e-03 3.83620191e-03\n",
      " 3.47333932e-03 3.14458048e-03 2.84675980e-03 2.57699834e-03\n",
      " 2.33267915e-03 2.11142450e-03 1.91107492e-03 1.72966989e-03\n",
      " 1.56543006e-03 1.41674105e-03 1.28213852e-03 1.16029454e-03\n",
      " 1.05000518e-03 9.50179134e-04 8.59827359e-04 7.78053645e-04\n",
      " 7.04046017e-04 6.37068919e-04 5.76456101e-04 5.21604154e-04\n",
      " 4.71966637e-04 4.27048741e-04 3.86402452e-04 3.49622146e-04\n",
      " 3.16340612e-04 2.86225426e-04 2.58975677e-04 2.34318991e-04\n",
      " 2.12008835e-04 1.91822074e-04 1.73556760e-04 1.57030125e-04\n",
      " 1.42076761e-04 1.28546982e-04 1.16305325e-04 1.05229208e-04\n",
      " 9.52077034e-05 8.61404321e-05 7.79365635e-05 7.05139072e-05\n",
      " 6.37980924e-05 5.77218246e-05 5.22242132e-05 4.72501624e-05\n",
      " 4.27498203e-05 3.86780802e-05 3.49941290e-05 3.16610390e-05\n",
      " 2.86453981e-05 2.59169751e-05 2.34484173e-05 2.12149764e-05\n",
      " 1.91942606e-05 1.73660103e-05 1.57118951e-05 1.42153301e-05\n",
      " 1.28613099e-05 1.16362581e-05 1.05278912e-05 9.52509549e-06\n",
      " 8.61781576e-06 7.79695445e-06 7.05428044e-06 6.38234657e-06\n",
      " 5.77441497e-06 5.22438949e-06 4.72675463e-06 4.27652020e-06\n",
      " 3.86917133e-06 3.50062315e-06 3.16717987e-06 2.86549773e-06\n",
      " 2.59255146e-06 2.34560392e-06 2.12217870e-06 1.92003527e-06\n",
      " 1.73714649e-06 1.57167833e-06 1.42197145e-06 1.28652453e-06\n",
      " 1.16397930e-06 1.05310684e-06 9.52795283e-07 8.62038686e-07\n",
      " 7.79926912e-07 7.05636521e-07 6.38422504e-07 5.77610818e-07\n",
      " 5.22591625e-07 4.72813173e-07 4.27776266e-07 3.87029260e-07\n",
      " 3.50163529e-07 3.16809371e-07 2.86632297e-07 2.59329681e-07\n",
      " 2.34627723e-07 2.12278703e-07 1.92058495e-07 1.73764324e-07\n",
      " 1.57212730e-07 1.42237726e-07 1.28689138e-07 1.16431095e-07\n",
      " 1.05340668e-07 9.53066397e-08 8.62283835e-08 7.80148595e-08\n",
      " 7.05836993e-08 6.38603803e-08 5.77774786e-08 5.22739921e-08\n",
      " 4.72947301e-08 4.27897586e-08 3.87138994e-08 3.50262788e-08\n",
      " 3.16899156e-08 2.86713515e-08 2.59403148e-08 2.34694182e-08\n",
      " 2.12338822e-08 1.92112881e-08 1.73813524e-08 1.57257240e-08\n",
      " 1.42277992e-08 1.28725565e-08 1.16464049e-08 1.05370481e-08\n",
      " 9.53336123e-09 8.62527853e-09 7.80369352e-09 7.06036713e-09\n",
      " 6.38784488e-09 5.77938272e-09 5.22887819e-09 4.73081105e-09\n",
      " 4.28018641e-09 3.87248521e-09 3.50361872e-09 3.16988812e-09\n",
      " 2.86794609e-09 2.59476528e-09 2.34760566e-09 2.12398898e-09\n",
      " 1.92167215e-09 1.73862680e-09 1.57301727e-09 1.42318246e-09\n",
      " 1.28761979e-09 1.16496990e-09 1.05400288e-09 9.53605860e-10\n",
      " 8.62771853e-10 7.80590036e-10 7.06236402e-10 6.38965103e-10\n",
      " 5.78101788e-10 5.23035615e-10 4.73215022e-10 4.28139746e-10\n",
      " 3.87358145e-10 3.50460993e-10 3.17078363e-10 2.86875634e-10\n",
      " 2.59549937e-10 2.34827047e-10 2.12458939e-10 1.92221572e-10\n",
      " 1.73911774e-10 1.57346136e-10 1.42358569e-10 1.28798305e-10\n",
      " 1.16529897e-10 1.05430109e-10 9.53874757e-11 8.63016325e-11\n",
      " 7.80810971e-11 7.06437131e-11 6.39146513e-11 5.78266324e-11\n",
      " 5.23183719e-11 4.73348027e-11 4.28261870e-11 3.87467836e-11\n",
      " 3.50559581e-11 3.17168514e-11 2.86957125e-11 2.59623434e-11\n",
      " 2.34894326e-11 2.12518891e-11 1.92275085e-11 1.73960846e-11\n",
      " 1.57389657e-11 1.42399426e-11 1.28834721e-11 1.16562315e-11\n",
      " 1.05460085e-11 9.54147872e-12 8.63265015e-12 7.81041898e-12\n",
      " 7.06634751e-12 6.39333031e-12 5.78426196e-12 5.23336929e-12\n",
      " 4.73487916e-12 4.28390656e-12 3.87578858e-12 3.50652840e-12\n",
      " 3.17257332e-12 2.87037061e-12 2.59703370e-12 2.34967601e-12\n",
      " 2.12585505e-12 1.92335037e-12 1.74016357e-12 1.57429625e-12\n",
      " 1.42441614e-12 1.28874689e-12 1.16595622e-12 1.05493392e-12\n",
      " 9.54347712e-13 8.63531469e-13 7.81152920e-13 7.06767977e-13\n",
      " 6.39488462e-13 5.78648240e-13 5.23581178e-13 4.73621142e-13\n",
      " 4.28546088e-13 3.87689880e-13 3.50830476e-13 3.17301740e-13\n",
      " 2.87103674e-13 2.59792188e-13 2.34923192e-13 2.12718732e-13\n",
      " 1.92290628e-13 1.74082970e-13 1.57429625e-13 1.42552636e-13\n",
      " 1.29007915e-13 1.16573418e-13 1.05471187e-13 9.54791801e-14\n",
      " 8.63753513e-14 7.81597009e-14 7.06101844e-14 6.39488462e-14\n",
      " 5.79536419e-14 5.24025268e-14 4.72955008e-14 4.28546088e-14\n",
      " 3.88578059e-14 3.50830476e-14 3.17523785e-14 2.86437540e-14\n",
      " 2.59792188e-14 2.35367281e-14 2.13162821e-14 1.93178806e-14\n",
      " 1.73194792e-14 1.57651669e-14 1.42108547e-14 1.28785871e-14\n",
      " 1.17683641e-14 1.06581410e-14 9.54791801e-15 8.65973959e-15\n",
      " 7.77156117e-15 7.10542736e-15 6.43929354e-15 5.77315973e-15\n",
      " 5.32907052e-15 4.66293670e-15 4.21884749e-15 3.77475828e-15\n",
      " 3.55271368e-15 3.10862447e-15 2.88657986e-15 2.66453526e-15\n",
      " 2.44249065e-15 2.22044605e-15 1.99840144e-15 1.77635684e-15\n",
      " 1.55431223e-15 1.33226763e-15 1.33226763e-15 1.11022302e-15\n",
      " 1.11022302e-15 8.88178420e-16 8.88178420e-16 8.88178420e-16\n",
      " 6.66133815e-16 6.66133815e-16 6.66133815e-16 4.44089210e-16\n",
      " 4.44089210e-16 4.44089210e-16 4.44089210e-16 4.44089210e-16\n",
      " 2.22044605e-16 2.22044605e-16 2.22044605e-16 2.22044605e-16\n",
      " 2.22044605e-16 2.22044605e-16 2.22044605e-16 2.22044605e-16\n",
      " 2.22044605e-16 2.22044605e-16 2.22044605e-16 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7efd06288a58>"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGdVJREFUeJzt3X+UXOV93/H3Z3e1i6UVQkKCGFZC8oloLQzmx0Yoxb9ODVhOXSlN7VrYTnBLoU5RnRynbXA4Ma6cnGBzTmkbEwfFcOq4RjKxa1vtwcUE/OO4sVStBMgWRHgRIC0iYZEE1g/QMjvf/jF3xexodvfuzoxm7tzP65w9mvvc5955rqT97rPf53vvKCIwM7N86Gj2AMzM7PRx0DczyxEHfTOzHHHQNzPLEQd9M7MccdA3M8sRB30zsxxx0DczyxEHfTOzHOlq9gAqLVy4MJYuXdrsYZiZZcqOHTteiohFU/VruaC/dOlSBgYGmj0MM7NMkfRcmn5O75iZ5YiDvplZjjjom5nliIO+mVmOOOibmeWIg76ZWY446JuZ5YiDvplZjqQK+pJWS9ojaVDSLVX2f0rSE5J2SXpY0gVl+0YlPZZ8bann4M3MbHqmvCNXUidwF3ANMARsl7QlIp4o6/Yo0B8RxyX9NvAF4MPJvlcj4tI6j9vMzGYgzUx/JTAYEXsjYgTYDKwt7xAR34+I48nmVqCvvsM0a75CocjfvvALisVis4diNmNpnr1zPrC/bHsIuHKS/jcA3y3bPkPSAFAAbo+Ib097lGZNNjIyyqV/9D2OjxTp7enksT+8lq4uL4lZ9qQJ+qrSFlU7Sh8D+oF3lzUviYgDkt4CPCLppxHxdMVxNwE3ASxZsiTVwM1Ol2Ix+I0//78cHynN8I+eGOWpF4+w4rx5TR6Z2fSlmaoMAYvLtvuAA5WdJF0N3AqsiYgTY+0RcSD5cy/wA+CyymMjYmNE9EdE/6JFUz4Z1Oy0OnhshJ8dODKubcHsWU0ajVlt0gT97cByScskdQPrgHFVOJIuA+6mFPBfLGufL6kneb0QuAooXwA2a3kLZs+it6dzXNsnv/44xWLVX3jNWtqU6Z2IKEhaDzwIdAL3RsRuSRuAgYjYAtwB9AJ/JQlgX0SsAd4K3C2pSOkHzO0VVT9mLe/Q8dd5dWR0XNvO5w5z8NgIi+b2NGlUZjOT6kNUIuIB4IGKts+Uvb56guP+Bri4lgGaNdvC3m6uuGABA88dYnZ3J8dHRrl8yVks7O1u9tDMps3lB2ZTiACIk+ULkkBK2s2yxUHfbAoHj42wc9/LFClV7owW42R6xyxrHPTNprBg9iwu7ptHh6C3p5MO4fSOZZaDvtkkisXgI1/exq79L/P2vnm89c1nOr1jmZZqIdcsrw4eG2HHc4cZDdg19ApI49I7rt6xrPFM32wSpcqd+XR1iCsumD/utdM7lkWe6ZtNQhKbblzFwWMjLOztJoKTr5N7UswyxTN9s2no6BBnz+nmpaMjhJP6lkGe6ZtNolgMrvuLrex47jBXXDCfr91wJR+9Z9vJ7U03rqKjwzN+yw7P9M0mMbaQWygGO547zODw0XHbrtW3rHHQN5tE5ULuhef2ejHXMk2tlpfs7++PgYGBZg/D7KRiMcYt3lZum7UCSTsion+qfp7pm5nliBdyzSbhhVxrN57pm03CC7nWbhz0zSbhhVxrN17INZuCF3ItC7yQa2Zmp/BCrtkkvJBr7cYzfbNJeCHX2o2DvtkkvJBr7cYLuWZT8EKuZUHahVzn9M2m0NGhcZ+QVbltliVO75iZ5YiDvtkkisVg+MiJcR+YUq3NLCuc3jGbQGW55qYbVwGc0uaSTcsSz/TNJlBZrnnw2EjVNrMscdA3m0BluebC3u6qbWZZ4pJNs0lUK890yaa1oro+e0fSakl7JA1KuqXK/k9JekLSLkkPS7qgbN/1kn6efF0/vcswa66x8szy4F6tzSwrpgz6kjqBu4D3AyuA6yStqOj2KNAfEZcA3wC+kBy7ALgNuBJYCdwmaX79hm/WWBNV6riCx7IqzUx/JTAYEXsjYgTYDKwt7xAR34+I48nmVqAvef0+4KGIOBQRh4GHgNX1GbpZY41V7/zqnzzMuo1bKRZj0nazLEgT9M8H9pdtDyVtE7kB+O4MjzVrGRNV6riCx7IsTdCvlrisOrWR9DGgH7hjOsdKuknSgKSB4eHhFEMya7yJKnVcwWNZlubmrCFgcdl2H3CgspOkq4FbgXdHxImyY99TcewPKo+NiI3ARihV76QYk1nDSWLTjatOqdSZqN0sC9LM9LcDyyUtk9QNrAO2lHeQdBlwN7AmIl4s2/UgcK2k+ckC7rVJm1kmTFSp4woey6opZ/oRUZC0nlKw7gTujYjdkjYAAxGxhVI6pxf4q+SbYF9ErImIQ5I+R+kHB8CGiDjUkCsxM7Mp+eYss0lMdCOWb9CyVuPn6ZvVqNoD1zo6NGG7WRb42TtmE3DJprUjB32zCbhk09qRc/pmk3BO37LCOX2zOpjo83D9ObmWVU7vmJnliIO+2QQme5Kmn7JpWeX0jlkVk5VlumTTsswzfbMqJivLdMmmZZmDvlkVk5VlumTTsswlm2YTmKws0yWb1mpcsmlWo8nKMl2yaVnl9I6ZWY446JuZ5YiDvtkEpqrFd62+ZZFz+mZVTFWL71p9yyrP9M2qmKoW37X6llUO+mZVTFWL71p9yyrX6ZtNYKpafNfqWytxnb5ZjaaqxXetvmWR0ztmZjnioG82AZdsWjtyesesCpdsWrvyTN+sCpdsWrty0DerwiWb1q5csmk2AZdsWpa4ZNOsRi7ZtHbk9I6ZWY446JtVkaYc0yWblkWpgr6k1ZL2SBqUdEuV/e+StFNSQdIHK/aNSnos+dpSr4GbNcpYOeav/snDrNu4lWLx1KCepo9ZK5oy6EvqBO4C3g+sAK6TtKKi2z7g48B9VU7xakRcmnytqXG8Zg2XphzTJZuWVWlm+iuBwYjYGxEjwGZgbXmHiHg2InYBxQaM0ey0SlOO6ZJNy6o01TvnA/vLtoeAK6fxHmdIGgAKwO0R8e1pHGt22kli042rJi3HTNPHrBWlCfrV/jdPJ4G5JCIOSHoL8Iikn0bE0+PeQLoJuAlgyZIl0zi1WWOkKcd0yaZlUZr0zhCwuGy7DziQ9g0i4kDy517gB8BlVfpsjIj+iOhftGhR2lObmdk0pQn624HlkpZJ6gbWAamqcCTNl9STvF4IXAU8MdPBmp0uacsxXbZpWTNleiciCpLWAw8CncC9EbFb0gZgICK2SPoV4FvAfOCfSvpPEXER8FbgbklFSj9gbo8IB31raWmfoOknbVoWpXoMQ0Q8ADxQ0faZstfbKaV9Ko/7G+DiGsdodlpVK8eslrtP28+slfiOXLMKacsxXbZpWeSnbJpVkfYJmn7SprUKP2XTrAZpyzFdtmlZ4/SOmVmOOOibVeGSTWtXTu+YVXDJprUzz/TNKqR9gqaftGlZ5KBvVsElm9bOXLJpVoVLNi1rXLJpVgOXbFq7cnrHzCxHHPTNqnDJprUrp3fMKrhk09qZZ/pmFVyyae3MQd+sgks2rZ25ZNOsCpdsWta4ZNOsBi7ZtHbl9I6ZWY446JtVmE4Zpks2LWuc3jErM50yTJdsWhZ5pm9WZjplmC7ZtCxy0DcrM50yTJdsWha5ZNOswnTKMF2yaa3CJZtmMzSdMkyXbFrWOL1jVmG6FTmu4LEs8UzfrMx0K3JcwWNZ45m+WZnpVuS4gseyxkHfrMx0K3JcwWNZ4+odswrTrchxBY+1grTVO6lm+pJWS9ojaVDSLVX2v0vSTkkFSR+s2He9pJ8nX9envwSz5hiryEkbwKfb36yZpgz6kjqBu4D3AyuA6yStqOi2D/g4cF/FsQuA24ArgZXAbZLm1z5ss8Zx9Y61szTVOyuBwYjYCyBpM7AWeGKsQ0Q8m+wrVhz7PuChiDiU7H8IWA1sqnnkZg3g6h1rd2nSO+cD+8u2h5K2NGo51uy0c/WOtbs0Qb/atCXt77GpjpV0k6QBSQPDw8MpT21Wf67esXaXJr0zBCwu2+4DDqQ8/xDwnopjf1DZKSI2AhuhVL2T8txmdSeJTTeuSl2NM93+Zs2WZqa/HVguaZmkbmAdsCXl+R8ErpU0P1nAvTZpM2tZrt6xdjZl0I+IArCeUrB+Erg/InZL2iBpDYCkX5E0BHwIuFvS7uTYQ8DnKP3g2A5sGFvUNWtFM6nEcfWOZYlvzjJLzKQSx9U71irqenOWWR7MpBLH1TuWNQ76ZomZVOK4eseyxukdszIzeY6On71jrcCfnGU2AzP5JCx/epZlidM7ZmVmWonjCh7LCs/0zRIzrcRxBY9liWf6ZomZVuK4gseyxEHfLDHTShxX8FiWuHrHrMxMK3FcwWPN5puzzMzsFF7INUt4IdfywDN9s4QXci0PHPTNEl7ItTzwQq5ZGS/kWlZ5Iddsmhy4LQ+8kGtGbYuxXsi1LPFM34zaFmO9kGtZ4qBvRm2LsV7ItSzxQq5ZolAoMjh8lAvP7aWjY3rzoWIxGD56AoE/JN2aws/TN5uGYjH46D3basrLf3LTo87rW8tzeseM2vPyzutbVjjom1F7Xt55fcsK5/TNErXk9OtxvFktnNM3m4Zac/r1WBMwOx08HTHDOX3LDwd9M5zTt/xwTt+M+tTZO6dvzeScvllK1Z6dM92Y75y+ZYWnI5Z79cjHO6dvWZEq6EtaLWmPpEFJt1TZ3yPp68n+bZKWJu1LJb0q6bHk68/rO3yz2tUjH++cvmXFlOkdSZ3AXcA1wBCwXdKWiHiirNsNwOGI+GVJ64DPAx9O9j0dEZfWedxmdSOJr91w5cl8/Ezy+ZLYdOOqk+sCZq0qzUx/JTAYEXsjYgTYDKyt6LMW+Ery+hvAe+UnTllGjOXjP/CnP+a6v9hGsTjz4oZPbnqUf3T7I6zbuLWm85g1Spqgfz6wv2x7KGmr2iciCsArwNnJvmWSHpX0Q0nvrHG8ZnVXr3y88/qWBWmCfrUZe+UUZqI+LwBLIuIy4FPAfZLOPOUNpJskDUgaGB4eTjEks/pZ2NvN5Uvm09khLq8hH186z1l0Ci5fcpbz+taS0gT9IWBx2XYfcGCiPpK6gHnAoYg4EREHASJiB/A0cGHlG0TExojoj4j+RYsWTf8qzGpQulUlSi8imOmtK6XjBCp9tdgtMGZAuqC/HVguaZmkbmAdsKWizxbg+uT1B4FHIiIkLUoWgpH0FmA5sLc+Qzerj4PHRti572VGA3bue7mm9M7OfYcZLQY7nd6xFjVl0E9y9OuBB4EngfsjYrekDZLWJN3uAc6WNEgpjTNW1vkuYJekxykt8H4iIg7V+yLMauH0juVJqjtyI+IB4IGKts+UvX4N+FCV474JfLPGMZo1VLX0zkxqz6qld1zDZq3Gd+Ra7jm9Y3nioG+5t2D2LC7um0dnjXfTOr1jWeCgb7lWLAYf+fI2du1/mbf3zeO+f33ljO7IBVfvWDY46Fuujd1QNRqwa+gVDh1/vaZzOb1jrc5B33KtXpU7Y+e64oL5dAou6ZvH2XNm1XGkZvXhoG+5Vq8bs+CNB7ddsvgsHh96pebn+Jg1goO+5drBYyPsSCp3dtRQuTPm8Kuvs2v/y4wWg4FnDznFYy3HQd9ybcHsWczu7gRgdncnC2bXlpJZMHsWs3tKt7/M7umq+Xxm9eagb7n20tERjp0oAHD8RKGmhVyAQ8dfP3m+Y3U4n1m9OehbbhWLwb/btJOxtPsVFyyoubZ+wexZzElm+nM807cW5KBvuTWWzwfoEHzxI5fNuEZ/zKHjr3N8ZBQozfSd07dW46BvuVWez5/T08XZc2q/g3ZhbzdXLDkLgGLA+k2PuoLHWoqDvuVWI/LvkvjT6y6nI/mFYYcreKzFOOhbbjUq/76wt/vkeV3BY63GQd9yq96VO2NcwWOtzEHfcqkRlTtjXMFjrcxB33KpEZU7Y1zBY63MQd9yacHsWcyeVd/KnTGVFTw337fTFTzWMhz0LZdeOjrCkSTvfuS1Ai8drd9sXBL/dd1lJ7e3P3uY4aMn6nZ+s1o46FsuFYvFcdtR5088qcwU1fv8ZjPloG+5UywG//a+nePaOjr8CeaWDw76ljsHj43w+NArJ7cvXTyPRXN76voeHadO9et6frOZctC33DkzefQClL4B7r9xVd0qd8YsmttDf7KYC3D1nT+iUChOcoTZ6eGgb7lSLAb//O6fnKzPLwJ7Dx6v+/tI4rNrLzq5ffTEKHv+/hd1fx+z6XLQt1wZPnKCnx14I/jO6e7kwnN7G/Jei3rHp4xu/dbPXLppTeegb7lSKIyO237od99JR0djvg3OOfMMLu078+T2Y0Ov8OKR1xryXmZpOehbboyMjHL1nT8c19bV1TlB79pJ4ksfu2Jc27/5ywHP9q2pHPQtF4rF4Nf/7Me8Wngj4F503ty6V+1U6qz4LeLx53/BgZePNfQ9zSbjoG+5sG/4KE/83dFxbff8Vn/dq3YqLZrbw0Vvnjuu7R1f+CGvveYnb1pzpAr6klZL2iNpUNItVfb3SPp6sn+bpKVl+z6dtO+R9L76Dd1saoVCkYHnXuI9d/5oXPuF58zh3Hlvavj7S+I7N19FT8V32ts++z2eP3jUd+raadc1VQdJncBdwDXAELBd0paIeKKs2w3A4Yj4ZUnrgM8DH5a0AlgHXAScB/y1pAsjYvxqmlkdjYyM8pNnXmKkWODmrz7GSJXy+Ac++Y6Gz/LHdHV18sjvvZur7nhjPaEAXHXHD1k2v5vb1r6Nc858E//wl85s2KKy2Zgpgz6wEhiMiL0AkjYDa4HyoL8W+Gzy+hvAF1X6jloLbI6IE8AzkgaT8/2kPsMfb+ybvRCjHHm1wJlv6iLg5GtJFCOatq/Z75+Hcb9WKPDbX9016f+TpzZcQ1dXmv/69XPegjlcct5cdh04Mq79mcMjfPy/lx4J0QV86WNvp2tWZ0v9nWb1/0LWxn3W7FmcM7fxP/zT/M8/H9hftj0EXDlRn4goSHoFODtp31px7PkzHu0kRkZGufAz/6cRp7Y28tSGa+jurt9jlNOSxLdufgf/7K4fnxL4xxSAG//H46d3YNZyenu6eOwPr6GrqzGBP81Zq/0OXJmInKhPmmORdJOkAUkDw8PDKYZ0qh37Ds/oOMuHMzrg55+7tikBf0xnZwffXv9OfvQf3tW0MVjrO3qiwODw0ak7zlCamf4QsLhsuw84MEGfIUldwDzgUMpjiYiNwEaA/v7+Ga1srVw2fyaHWRv78m++nVmzuljUe0bL5Ms7OsSSs+fy1IbV/OSZlxiNUf74fz3J0wd905aVzO3pathd4gCaqnogCeJPAe8Fnge2Ax+JiN1lfW4GLo6ITyQLub8REf9C0kXAfZTy+OcBDwPLJ1vI7e/vj4GBgRldjHP6HncxghMjwbUrzj3tefuZKhaD5185zs5nDzH3TZ0t+XfaKu/fzuOuNacvaUdE9E/Vb8rviiRHvx54EOgE7o2I3ZI2AAMRsQW4B/hqslB7iFLFDkm/+ykt+haAmxtZudPd3cm7/8G5jTq9WUN0dIjF8+eweP6cZg/FcmDKmf7pVstM38wsr9LO9Juf5DQzs9PGQd/MLEcc9M3McsRB38wsRxz0zcxyxEHfzCxHHPTNzHLEQd/MLEcc9M3McqTl7siVNAw81+xxzMBC4KVmD+I08zXng685Gy6IiEVTdWq5oJ9VkgbS3ALdTnzN+eBrbi9O75iZ5YiDvplZjjjo18/GZg+gCXzN+eBrbiPO6ZuZ5Yhn+mZmOeKgXyeS/r2kkLQw2Zak/yZpUNIuSZc3e4z1IukOSX+bXNe3JJ1Vtu/TyTXvkfS+Zo6zniStTq5pUNItzR5PI0haLOn7kp6UtFvS7yTtCyQ9JOnnyZ9t94HUkjolPSrpfyfbyyRtS67565K6mz3GenHQrwNJi4FrgH1lze8HlidfNwFfasLQGuUh4G0RcQmlz0/+NICkFZQ+KvMiYDXwZ5I6mzbKOkmu4S5K/6YrgOuSa203BeD3IuKtwCrg5uQ6bwEejojllD7nuh1/6P0O8GTZ9ueBO5NrPgzc0JRRNYCDfn3cCfxHoHyBZC3wl1GyFThL0pubMro6i4jvRUQh2dwK9CWv1wKbI+JERDwDDAIrmzHGOlsJDEbE3ogYATZTuta2EhEvRMTO5PURSkHwfErX+pWk21eAX2/OCBtDUh/wT4AvJ9sC/jHwjaRLW12zg36NJK0Bno+Ixyt2nQ/sL9seStrazb8Cvpu8btdrbtfrmpCkpcBlwDbg3Ih4AUo/GIBzmjeyhvgvlCZtxWT7bODlsolNW/17dzV7AFkg6a+BX6qy61bgD4Brqx1WpS0zpVKTXXNEfCfpcyullMDXxg6r0j8z1zyJdr2uqiT1At8EfjciflGa+LYnSR8AXoyIHZLeM9ZcpWvb/Hs76KcQEVdXa5d0MbAMeDz5xugDdkpaSWl2sLisex9woMFDrZuJrnmMpOuBDwDvjTfqfjN9zZNo1+s6haRZlAL+1yLifybNfy/pzRHxQpKifLF5I6y7q4A1kn4NOAM4k9LM/yxJXclsv63+vZ3eqUFE/DQizomIpRGxlFJwuDwi/g7YAvxWUsWzCnhl7FfkrJO0Gvh9YE1EHC/btQVYJ6lH0jJKi9j/rxljrLPtwPKkoqOb0mL1liaPqe6SXPY9wJMR8Z/Ldm0Brk9eXw9853SPrVEi4tMR0Zd8/64DHomIjwLfBz6YdGura/ZMv3EeAH6N0mLmceBfNnc4dfVFoAd4KPkNZ2tEfCIidku6H3iCUtrn5ogYbeI46yIiCpLWAw8CncC9EbG7ycNqhKuA3wR+KumxpO0PgNuB+yXdQKlC7UNNGt/p9PvAZkl/BDxK6YdhW/AduWZmOeL0jplZjjjom5nliIO+mVmOOOibmeWIg76ZWY446JuZ5YiDvplZjjjom5nlyP8HWKAx7R0HmtEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# to compute sigmoidGradient, squash the number to sig first,\n",
    "# and then calculate sig * (1 - sig)\n",
    "# the greater the activation, the higher the gradient\n",
    "a = np.linspace(-50,50,1000)\n",
    "plt.scatter(a, sigmoidGradient(a), s=5)\n",
    "plt.savefig('sigmoidGradient@2x.png', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_pass(Theta1, Theta2, X, y):   \n",
    "#     m = X.shape[0]\n",
    "    # feedforward to produce the prediction output\n",
    "    # a1(X):input layer, a2 hidden layer activation, a3 output layer\n",
    "    bias1 = np.ones((X.shape[0], 1))\n",
    "    a1 = np.concatenate([bias1, X], axis=1) \n",
    "    a2 = sigmoid(a1 @ Theta1.T)\n",
    "    \n",
    "    bias2 = np.ones((a2.shape[0], 1))\n",
    "    a2 = np.concatenate([bias2, a2], axis=1)\n",
    "    a3 = sigmoid(a2 @ Theta2.T) # output layer\n",
    " \n",
    "    return (a1, a2, a3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_costs_loop(Theta1, Theta2, num_labels, X, y, lambda_):     \n",
    "    m = X.shape[0]\n",
    "    J = 0\n",
    "    a1, a2, a3 = forward_pass(Theta1, Theta2, X, y)\n",
    "    # compute cost for each output unit (class)\n",
    "    for j in range(num_labels):\n",
    "        Jj = sum(-y[:,j] * np.log(a3[:,j]) - (1-y[:,j])*np.log(1-a3[:,j]))\n",
    "        print(y[:,j].shape)\n",
    "        J += Jj \n",
    "        print(f'compute cost for {j}', Jj)\n",
    "    \n",
    "    J = 1/m * J\n",
    "    reg_J = J + lambda_/(2*m) * (np.sum(Theta1[:,1:]**2) + np.sum(Theta2[:,1:]**2))\n",
    "    return (J, reg_J)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_costs_vectorized(Theta1, Theta2, num_labels, X, y, lambda_):     \n",
    "    m = X.shape[0]\n",
    "    J = 0\n",
    "    a1, a2, a3 = forward_pass(Theta1, Theta2, X, y)\n",
    "    # below is vectorized J implementation\n",
    "    J = np.sum(-y * np.log(a3) - (1 - y) * np.log(1 - (a3)))\n",
    "  \n",
    "    J = 1/m * J\n",
    "    reg_J = J + lambda_/(2*m) * (np.sum(Theta1[:,1:]**2) + np.sum(Theta2[:,1:]**2))\n",
    "    return (J, reg_J)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2876291651613189"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.2876291651613189"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J, reg_J = compute_costs_vectorized(Theta1, Theta2, num_labels, X, y_encode, lambda_=0.0)\n",
    "J\n",
    "reg_J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5000,)\n",
      "compute cost for 0 71.85847553673608\n",
      "(5000,)\n",
      "compute cost for 1 105.77238482303046\n",
      "(5000,)\n",
      "compute cost for 2 151.13409896391363\n",
      "(5000,)\n",
      "compute cost for 3 197.31569852900847\n",
      "(5000,)\n",
      "compute cost for 4 148.61622965727082\n",
      "(5000,)\n",
      "compute cost for 5 149.92213260735022\n",
      "(5000,)\n",
      "compute cost for 6 101.33146796645212\n",
      "(5000,)\n",
      "compute cost for 7 141.14000796503075\n",
      "(5000,)\n",
      "compute cost for 8 173.73617842845815\n",
      "(5000,)\n",
      "compute cost for 9 197.31915132934316\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.28762916516131876, 0.28762916516131876)"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_costs_loop(Theta1, Theta2, num_labels, X, y_encode, lambda_=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1438.1458258065945, 0.2876291651613189)"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_costs_vectorized(Theta1, Theta2, num_labels, X, y_encode, lambda_=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_gradient(Theta1, Theta2, num_labels,X, y,lambda_): \n",
    "    m = X.shape[0] \n",
    "    Theta1_grad = np.zeros((Theta1.shape))\n",
    "    Theta2_grad = np.zeros((Theta2.shape))\n",
    "    a1, a2, a3 = forward_pass(Theta1, Theta2, X, y)\n",
    "#     for i in range(m):\n",
    "    for i in range(1):\n",
    "        a1i = a1[i,:] # 1 X 401\n",
    "        a2i = a2[i,:] # 1 X 26\n",
    "        a3i = a3[i,:] # 1 x 10\n",
    "        delta2 = a3i - y[i,:]\n",
    "\n",
    "        # Theta2.T @ d2.T is the 10 errors at layer 3 back propagated to layer 2\n",
    "        # sigmoidGradient(np.hstack((1,a1i @ Theta1.T))) is the gradient for each activation in layer 2\n",
    "        delta1 = Theta2.T @ delta2 * sigmoidGradient(np.hstack((1,a1i @ Theta1.T)))\n",
    "        \n",
    "        Theta1_grad= Theta1_grad + delta1[1:][:,np.newaxis] @ a1i[:,np.newaxis].T\n",
    "        Theta2_grad = Theta2_grad + delta2.T[:,np.newaxis] @ a2i[:,np.newaxis].T\n",
    "        \n",
    "    Theta1_grad = 1/m * Theta1_grad\n",
    "    Theta2_grad = 1/m*Theta2_grad\n",
    "    \n",
    "    Theta1_grad_reg = Theta1_grad + (lambda_/m) * np.hstack((np.zeros((Theta1.shape[0],1)),Theta1[:,1:]))\n",
    "    Theta2_grad_reg = Theta2_grad + (lambda_/m) * np.hstack((np.zeros((Theta2.shape[0],1)),Theta2[:,1:]))\n",
    "    grad = np.concatenate([Theta1_grad_reg.ravel(), Theta2_grad_reg.ravel()])\n",
    "\n",
    "    return grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g is [0.19661193 0.04782555 0.07309204 0.00694906 0.24964942 0.20787312\n",
      " 0.00349239 0.01062254 0.0016299  0.00348869 0.01998466 0.06663376\n",
      " 0.02119201 0.04405998 0.04775892 0.00420003 0.10022623 0.00793126\n",
      " 0.08754214 0.12899417 0.16136511 0.02325683 0.01589205 0.02672678\n",
      " 0.05337726 0.09571899]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-5.12078395e-08,  0.00000000e+00,  0.00000000e+00, ...,\n",
       "        1.26052392e-06,  1.22281085e-06,  1.38964573e-07])"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_gradient(Theta1, Theta2, num_labels,X, y_encode,lambda_=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 599,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.01757087838907182"
      ]
     },
     "execution_count": 599,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "# error is the 10 error values for each unit in layer 3\n",
    "error = np.array([-4.26598801e-03,1.12661530e-04,1.74127856e-03,2.52696959e-03\n",
    "                   ,1.84032321e-05,9.36263860e-03,3.99270267e-03,5.51517524e-03\n",
    "                   ,4.01468105e-04,6.48072305e-03])\n",
    "# wts is the 10 weights connecting the first unit in layer 2 to the 10 output units\n",
    "wts = np.array([-0.46089119,-0.76100352,-0.61785176,-0.68934072,-0.67832479,-0.59664339\n",
    "                ,-0.87794907,-0.52746527,-0.7490154,-0.6665468,])\n",
    "# error*wts is the error for first unit in layer2, which is basically a weighted sum of error passed from layer 3 \n",
    "np.sum(error*wts)\n",
    "\n",
    "# component 2:sigmoidGradient(np.hstack((1,a1i @ Theta1.T)))\n",
    "# this is basically the forward pass value z= singal @ weight befor applying sigmoid to get the acitvation value\n",
    "# for each unit in layer 2, use 400 pixel values time their correspoding weights to each unit, then apply sigmoid, use the sigmoid activation to times (1-sigmoid) to get the gradient for each unit.\n",
    "# then use this activation gradient to time the errors for each unit, which is the weighted sum of errors form layer3\n",
    "\n",
    "\n",
    "# compute the theta2_gradient\n",
    "# basically use the error times the activation. the error for the first unit in the output layer times each of th 26 activations to get the gradient for each weight corresponding to the first unit in the layer 3.\n",
    "# similarly idea for theta1 gradient, use delta1 times each activation in the input layer to get gradients.\n",
    "# we need to do this for all training examples. then calculate the average, use this to update the Theta1 and Theta2 (weights)\n"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}