3-backpropagation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Backpropagation\n",
    "\n",
    "There are many ways to assign weight to a network. The most popular being backpropagation.\n",
    "\n",
    "Similar other machine learning algorithm, backpropagation fit a model to a training data.\n",
    "\n",
    "With a single case of training data inputted into the network, the output is compared to to correct value (target).\n",
    "\n",
    "If the output of the network is correct then nothing needs to be done, but if it is wrong than backpropagation looks at the difference between the correct value and network output to see how to change the weights so the network would output the correct value.\n",
    "\n",
    "The difference between network output and target is calculated by an loss function. There are many types of loss functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One loss function is the mean square error, which is the mean of square error function.\n",
    "\n",
    "Square error is for compare a single value with a target, and mean square error for a vector of values.\n",
    "\n",
    "Square error (SE):\n",
    "\n",
    "$$E_\\text{SE}(\\text{target}, \\text{output}) = \\frac{1}{2} (\\text{target} - \\text{output})^2$$\n",
    "\n",
    "Note, the $\\frac{1}{2}$ is there so the derivative of square error becomes: $\\frac{dE_\\text{SE}}{d\\text{output}} = \\text{output} - \\text{target}$\n",
    "\n",
    "Mean square error (MSE):\n",
    "\n",
    "$$E_\\text{MSE}(\\vec{\\text{target}}, \\vec{\\text{output}}) = \\frac{1}{2n} \\sum_{i=0}^n(\\text{target}_i - \\text{output}_i)^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the loss function ($E(\\cdot)$) is called the loss ($E$).\n",
    "\n",
    "$$E = E(\\vec t, \\vec {o_n})$$\n",
    "\n",
    "Where $\\vec t$ is the vector of targets, $\\vec {o_n}$ is the vector of network output, and $n$ is the size of these vectors.\n",
    "\n",
    "Backpropagation alter the weights by performing some form of gradient descent on the loss."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient descent\n",
    "\n",
    "Gradient descent is an optimization algorithm. It is for finding the minimum of a function. \n",
    "\n",
    "By the way, gradient is just the multi-variable version of derivative.\n",
    "\n",
    "![gradient_descent.png](img/gradient_descent.png)\n",
    "\n",
    "In backpropagation, we are finding a set of weights the yields the minimum loss (difference between target and output). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea is to use gradient descent to find the weights that gives the minimum loss.\n",
    "\n",
    "Since the slope/derivative/gradient always points away from the local minimum, by subtracting the gradient of the loss to the weights, the weights will move toward a lower loss. This is repeatedly done until the minimum is reached.\n",
    "\n",
    "There will be visualization of how this work in practice next lesson. Just know that gradient descent changes the input variables (e.g., weights) to minimize the output (e.g., loss)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The gradient descent update rule\n",
    "\n",
    "$\\vec w_{n+1} = \\vec w_n - \\eta \\nabla E $\n",
    "\n",
    "Where $\\vec w_n$ is a vector of all the weights at the $n$ iteration, $\\eta$ is a learning rate, $E$ is the loss. $\\nabla E$ means the gradient of loss.\n",
    " \n",
    "Usually the learning rate $\\eta$ is positive and set far below 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or the change in $\\vec w$ equals negative learning rate times gradient of loss:\n",
    "\n",
    "$\\Delta \\vec w = - \\eta \\nabla E$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Backpropagation is just repeatedly applying some form of gradient descent update rule to each weight. The most popular update rule right now is the [Adam optimizer](https://arxiv.org/abs/1412.6980)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic gradient descent\n",
    "Generally, batch gradient descent or stochastic gradient descent is used. The difference with gradient descent is how much data is used to calculate the loss.\n",
    "\n",
    "- gradient descent: The loss of all the training data is added up.\n",
    "\n",
    "- stochastic gradient descent: Only use a single case of training data.\n",
    "\n",
    "- batch gradient descent: The loss of a handful of training data is added up.\n",
    "\n",
    "The original backpropagation uses stochastic gradient descent, so that is the one we will be using for now.\n",
    "\n",
    "We will be considering the loss of a single case of training data at each iteration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Breaking down the update rule\n",
    "\n",
    "The gradient is a vector consisting of all the partial derivative with respect to the weights:\n",
    "\n",
    "$$\\nabla E = (\\frac{\\partial E}{\\partial w_{ij}}, \\dots)$$ $$ i, j \\in N $$\n",
    "\n",
    "Where $E$ is the loss, $w_{ij}$ is the weight between neuron $i$ and $j$, and $N$ denotes the set of all the neurons in the network.\n",
    "\n",
    "Note $\\partial$ is the partial derivative symbol. Partial derivative calculated just like derivatives. The difference is just the multi-variable setting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The update rule describes how to change all the weight:\n",
    "\n",
    "$$\\Delta \\vec w = - \\eta \\nabla E$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a single weight, this formula becomes:\n",
    "\n",
    "$$\\Delta w_{ij} = - \\eta \\frac{\\partial E}{\\partial w_{ij}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So in backpropagation we need to calculate a bunch of partial derivative to update the weights."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating partial derivatives\n",
    "\n",
    "This section is going to be very long. Calculating the partial derivatives is not difficult at all, it only requires the basics learned in a calculus class. But due to the complex structure of the neuron network, there is a lot stuff to keep track of. This section might require several rereading to fully comprehend.\n",
    "\n",
    "In order to calculate a partial derivative of the loss with respect to a weight, we need to apply the chain rule for derivative.\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial w_{ij}} = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial o_j}{\\partial w_{ij}} $$\n",
    "\n",
    "Where $w_{ij}$ is the weight of $i$ input of $j$ neuron, $o_j$ is the output of $j$ neuron."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this formula is saying is that the partial derivative of the loss with respect to a weight ($\\frac{\\partial E}{\\partial w_{ij}}$) can be split into: partial derivative of the loss with respect to the neuron output ($\\frac{\\partial E}{\\partial o_j}$) and partial derivative of the neuron output with respect to the weight ($\\frac{\\partial o_j}{\\partial w_{ij}}$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![backpropagation_chain_rule.png](img/backpropagation_chain_rule.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The $\\frac{\\partial o_j}{\\partial w_{ij}}$ part\n",
    "\n",
    "Recall the neuron function: $n(\\vec x) = h(\\vec w \\cdot \\vec x)$\n",
    "\n",
    "So the output of neuron $j$ ($o_j$) is equal to:\n",
    "\n",
    "$o_j = h_j(\\vec w \\cdot \\vec x)$\n",
    "\n",
    "Where $h_j(\\cdot)$ is the activation function of neuron $j$.\n",
    "\n",
    "Or breaking it down to two parts, the weighted sum ($s_j$) and the activation function of neuron $j$ ($h_j(\\cdot)$):\n",
    "\n",
    "$s_j = \\vec w \\cdot \\vec x$ \n",
    "\n",
    "$o_j = h_j(s_j)$\n",
    "\n",
    "Using the chain rule:\n",
    "\n",
    "$$ \\frac{\\partial o_j}{\\partial w_{ij}} = \\frac{\\partial o_j}{\\partial s_j} \\cdot \\frac{\\partial s_j}{\\partial w_{ij}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The partial derivative of first term is:\n",
    "\n",
    "$$\\frac{\\partial o_j}{\\partial s_j} = \\frac{\\partial}{\\partial s_j} h_j(s_j) = \\frac{\\partial h_j(s_j)}{\\partial s_j}$$\n",
    "\n",
    "Which is the derivative of the activation function with respect to the weighted sum.\n",
    "\n",
    "If sigmoid is the activation function.\n",
    "\n",
    "$$h(x) = \\sigma(x) = \\frac{1}{1+e^{\\text{-}x}}$$\n",
    "\n",
    "Then the derivative would be:\n",
    "\n",
    "$$ \\frac{\\partial h_j(s_j)}{\\partial s_j} = \\sigma^\\prime(s_j) = \\frac{e^{\\text{-}s_j}}{(1+e^{\\text{-}s_j})^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the second term $\\frac{\\partial s_j}{\\partial w_{ij}}$:\n",
    "\n",
    "Recall the weighted sum function: \n",
    "\n",
    "$s_j = \\vec w \\cdot \\vec x$ where $\\vec x$ is a vector of the inputs and $\\vec w$ is the weights associated with them.\n",
    "\n",
    "$$s_j = \\sum_{i \\in I} w_{ij} \\cdot x_i$$\n",
    "\n",
    "where $I$ is the set of all inputs connect to neuron $j$, and $w_{ij}$ is the weight associated with the connection from input $i$ to neuron $j$. \n",
    "\n",
    "Or in the case where the input is the output of the neuron that comes before it.\n",
    "\n",
    "$$s_j = \\sum_{l \\in L} w_{lj} \\cdot o_l$$\n",
    "\n",
    "Where $L$ is the set of all the neurons whose output is the input for neuron $j$, the weight $w_{lj}$ is associated with the connection from neuron $l$ to neuron $j$, and $o_l$ is the output of neuron $l$ which is also the input of neuron $j$.\n",
    "\n",
    "For a fully connected network, the set $L$ is the layer of neurons before neuron $j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{\\partial s_j}{\\partial w_{ij}} = \\frac{\\partial}{\\partial w_i} (\\vec w \\cdot \\vec x) = \\frac{\\partial}{\\partial w_{ij}} \\sum_{l \\in L} w_{lj} \\cdot o_l = \\frac{\\partial}{\\partial w_{ij}} (w_{ij} \\cdot o_i) = o_i$$\n",
    "\n",
    "So this term is just the input of neuron $j$ associated with weight $w_{ij}$, aka output of neuron $i$ ($o_i$). \n",
    "\n",
    "In the case, which neuron $j$ is in the next layer after input layer, then: \n",
    "\n",
    "$$o_i = x_i$$\n",
    "\n",
    "Since the output of input layer is the input of the network ($\\vec x = (x_1, \\dots , x_n)$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting it together:\n",
    "\n",
    "$$ \\frac{\\partial o_j}{\\partial w_{ij}} = \\frac{\\partial h_j(s_j)}{\\partial s_j} \\cdot o_i$$\n",
    "\n",
    "Which is the derivative of the activation function $h(\\cdot)$ with respect to the weighted sum ($s_j$) multiplied by the output of neuron 𝑖."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The $\\frac{\\partial E}{\\partial o_j}$ part\n",
    "\n",
    "If neuron $j$ is in the output layer then it would be simple, $\\frac{\\partial E}{\\partial o_j}$ is the derivative of the loss function with respect to one of the output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the error function is mean square error (MSE):\n",
    "\n",
    "$$E = E_\\text{MSE}(\\vec{\\text{target}}, \\vec{\\text{output}}) = \\frac{1}{2n} \\sum_{i=0}^n(\\text{target}_i - \\text{output}_i)^2$$\n",
    "\n",
    "Where $\\vec t$ is the vector containing the outputs, $\\vec {o_n}$ is the vector containing network output, and n is the size of these vectors.\n",
    "\n",
    "Then the partial derivative will be:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial o_j} = \\frac{\\partial E_\\text{MSE}(\\vec t, \\vec {o_n})}{\\partial o_j} = \\frac{\\partial}{\\partial o_j} \\biggl( \\frac{1}{2n} \\sum_{i=0}^n(\\text{target}_i - \\text{output}_i)^2 \\biggr)$$\n",
    "\n",
    "Since neuron $j$ is in output layer, its output $o_j$ is in the vector $\\vec {o_n}$.\n",
    "\n",
    "$$ \\frac{\\partial}{\\partial o_j} \\biggl( \\frac{1}{2n} \\sum_{i=0}^n(\\text{target}_i - \\text{output}_i)^2 \\biggr) = \\frac{\\partial}{\\partial o_j} \\frac{1}{2n}\\cdot (t_j - o_j)^2 = \\frac{1}{n} \\cdot (o_j - t_j)$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general form for $\\frac{\\partial E}{\\partial o_j}$ is:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial o_j} = \\frac{\\partial E(t_j, o_j)}{\\partial o_j}$$\n",
    "\n",
    "Where $o_j$ is the output of neuron $j$, and $t_j$ is the target value for neuron $j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, if neuron $j$ is not in the output layer then $\\frac{\\partial E}{\\partial o_j}$ needs to be broken down with the chain rule.\n",
    "\n",
    "$$ \\frac{\\partial E}{\\partial o_j} = \\sum_{l \\in L} \\frac{\\partial E}{\\partial s_l} \\cdot \\frac{\\partial s_l}{\\partial o_j}$$\n",
    "\n",
    "Where $L$ denotes the set that contains all neurons which takes the output of neuron $j$ ($o_j$) as input.\n",
    "\n",
    "For fully connected network, the set $L$ is consist of the next layer after neuron $j$.\n",
    "\n",
    "$s_l$ are the weighted sum component of these neurons.  \n",
    "\n",
    "$$\\sum_{l \\in L} \\frac{\\partial E}{\\partial s_l} \\cdot \\frac{\\partial s_l}{\\partial o_j} = \\sum_{l \\in L} \\frac{\\partial E}{\\partial o_l} \\cdot \\frac{\\partial o_l}{\\partial s_l} \\cdot \\frac{\\partial s_l}{\\partial o_j}$$\n",
    "\n",
    "Let us break down what each term means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last term $\\frac{\\partial s_l}{\\partial o_j}$ is quite simple.\n",
    "\n",
    "Recall the weighted sum function: \n",
    "\n",
    "$s_l = \\vec w \\cdot \\vec x$ where $\\vec x$ is a vector of the inputs and $\\vec w$ is the weights associated with them.\n",
    "\n",
    "In this case the input is the output of the neuron that comes before it.\n",
    "\n",
    "$$s_l = \\sum_{i \\in I} w_{il} \\cdot o_i$$\n",
    "\n",
    "Where $I$ is the set of neuron that which neuron $l$ takes as input. \n",
    "\n",
    "By the way, the neuron $j$ is in set $I$. Or in mathematics expression: $I = \\{j,  \\dots\\}$.\n",
    "\n",
    "Taking the partial derivative of $s_l$ with respect to $o_j$ will be:\n",
    "\n",
    "$$\\frac{\\partial s_l}{\\partial o_j} = \\frac{\\partial}{\\partial o_j} \\sum_{i \\in I} w_{il} \\cdot o_i = \\frac{\\partial}{\\partial o_j} w_{jl} \\cdot o_j = w_{jl} $$\n",
    "\n",
    "This is just the weight between neuron $j$ and neuron $l$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have already seen the $\\frac{\\partial o_l}{\\partial s_l}$. \n",
    "\n",
    "It is just the derivative of the activation function. Recall we have concluded that: \n",
    "\n",
    "$$\\frac{\\partial o_j}{\\partial s_j} = \\frac{\\partial h_j(s_j)}{\\partial s_j}$$\n",
    "\n",
    "Same here:\n",
    "\n",
    "$$\\frac{\\partial o_l}{\\partial s_l} = \\frac{\\partial h_l(s_l)}{\\partial s_l}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also have seen the $\\frac{\\partial E}{\\partial o_l}$ term. Recall at the start of this section:\n",
    "\n",
    "$$ \\frac{\\partial E}{\\partial o_j} = \\sum_{l \\in L} \\frac{\\partial E}{\\partial o_l} \\cdot \\frac{\\partial o_l}{\\partial s_l} \\cdot \\frac{\\partial s_l}{\\partial o_j}$$\n",
    "\n",
    "Note the $\\frac{\\partial E}{\\partial o_j}$ and $\\frac{\\partial E}{\\partial o_l}$ just differ on which neuron the output is from.\n",
    "\n",
    "We can substitute $o_j$ for $o_l$ recursively in this formula.\n",
    "\n",
    "Therefore, the $\\frac{\\partial E}{\\partial o_l}$ term can be defined recursively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Example:\n",
    "\n",
    "Say neuron $j$ is in the second last layer, we can expend out the recursion to:\n",
    "\n",
    "$$ \\frac{\\partial E}{\\partial o_j} = \\sum_{l_{n-1} \\in L_{n-1}} \\Bigl( \\sum_{l_{n} \\in L_n} \\frac{\\partial E}{\\partial o_{l_n}} \\cdot \\frac{\\partial o_{l_n}}{\\partial s_{l_n}} \\cdot \\frac{\\partial s_{l_n}}{\\partial o_j} \\Bigr) \\cdot \\frac{\\partial o_{l_{n-1}}}{\\partial s_{l_{n-1}}} \\cdot \\frac{\\partial s_{l_{n-1}}}{\\partial o_j}$$\n",
    "\n",
    "if we substitute the out the term we define in this section then:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial o_j} = \\sum_{l_{n-1} \\in L_{n-1}} \\Bigl( \\sum_{l_{n} \\in L_n} \\frac{\\partial E}{\\partial o_{l_n}} \\cdot \\frac{\\partial h_{l_n}(s_{l_n})}{\\partial s_{l_n}} \\cdot w_{jl_n} \\Bigr) \\cdot \\frac{\\partial h_{l_{n-1}}(s_{l_{n-1}})}{\\partial s_{l_{n-1}}} \\cdot w_{jl_{n-1}}$$\n",
    "\n",
    "Where $L_n$ is the set of neurons in the last layer, and $L_{n-1}$ is the set of neurons in the second last layer.\n",
    "\n",
    "As it mention at the start of this section, if neuron $j$ in $\\frac{\\partial E}{\\partial o_j}$ is the output layer then $\\frac{\\partial E}{\\partial o_j}$ is the partial derivative of the error function ($\\frac{\\partial E(t_j, o_j)}{\\partial o_j}$).\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial o_{l_n}} = \\frac{\\partial E(t_j, o_j)}{\\partial o_j} \\quad \\text{if neuron $j$ is an output neuron}$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Summing it up:\n",
    "\n",
    "$$ \\frac{\\partial E}{\\partial o_j} = \\begin{cases} \\frac{\\partial E(t_j, o_j)}{\\partial o_j} & \\text{if $j$ is an output neuron} \\\\ \\sum_{l \\in L} \\frac{\\partial E}{\\partial o_l} \\cdot \\frac{\\partial h_l(s_l)}{\\partial s_l} \\cdot w_{jl} & \\text{otherwise}\\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Putting it together\n",
    "\n",
    "Now each term of the of the partial derivative is defined we can put it together. \n",
    "\n",
    "$$\\frac{\\partial E}{\\partial w_{ij}} = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial o_j}{\\partial w_{ij}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In \"the $\\frac{\\partial o_j}{\\partial w_{ij}}$ part\", we learn to expand out $\\frac{\\partial o_j}{\\partial w_{ij}}$:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial o_j}{\\partial w_{ij}} = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial o_j}{\\partial s_j} \\cdot \\frac{\\partial s_j}{\\partial w_i}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also learned that $\\frac{\\partial s_j}{\\partial w_i} = o_i$. Therefore:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial o_j}{\\partial s_j} \\cdot \\frac{\\partial s_j}{\\partial w_i} = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial o_j}{\\partial s_j} \\cdot o_i$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\frac{\\partial o_j}{\\partial s_j}$ is the partial derivative of the activation function. So:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial o_j}{\\partial s_j} \\cdot o_i = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} \\cdot o_i$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we learned the recursive definition of $\\frac{\\partial E}{\\partial o_j}$:\n",
    "\n",
    "$$ \\frac{\\partial E}{\\partial o_j} = \\begin{cases} \\frac{\\partial E(t_j, o_j)}{\\partial o_j} & \\text{if $j$ is an output neuron} \\\\ \\sum_{l \\in L} \\frac{\\partial E}{\\partial o_l} \\cdot \\frac{\\partial h_l(s_l)}{\\partial s_l} \\cdot w_{jl} & \\text{otherwise}\\end{cases}$$\n",
    "\n",
    "Together with:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial w_{ij}} = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} \\cdot o_i$$\n",
    "\n",
    "The is complete formula for partial derivative of loss with respect to any weight in the network."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using the shorthand $\\delta$:\n",
    "\n",
    "However, it is more common to rewrite it as with a shorthand variable $\\delta$:\n",
    "\n",
    "$$\\delta_j = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial o_j}{\\partial s_j} = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j}$$\n",
    "\n",
    "Expanding out the $\\frac{\\partial E}{\\partial o_j}$ gives:\n",
    "\n",
    "$$ \\delta_j = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} = \\begin{cases} \\frac{\\partial E(t_j, o_j)}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{if $j$ is an output neuron} \\\\ \\Bigl( \\sum_{l \\in L} \\frac{\\partial E}{\\partial o_l} \\cdot \\frac{\\partial h_l(s_l)}{\\partial s_l} \\cdot w_{jl} \\Bigr) \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{otherwise}\\end{cases}$$\n",
    "\n",
    "By substitute of $\\delta_l = \\frac{\\partial E}{\\partial o_l} \\cdot \\frac{\\partial h_l(s_l)}{\\partial s_l}$, gives:\n",
    "\n",
    "$$\\delta_j = \\begin{cases} \\frac{\\partial E(t_j, o_j)}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{if $j$ is an output neuron} \\\\ \\bigl( \\sum_{l \\in L} \\delta_l \\cdot w_{jl} \\bigr) \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{otherwise}\\end{cases}$$\n",
    "\n",
    "Now the formula can be rewritten with $\\delta$:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial w_{ij}} = \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} \\cdot o_i = \\delta_j \\cdot o_i$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The backpropagation partial derivative formula\n",
    "\\begin{align}\n",
    "& \\frac{\\partial E}{\\partial w_{ij}} = \\delta_j \\cdot o_i \\\\\n",
    "& \\delta_j = \\begin{cases} \\frac{\\partial E(t_j, o_j)}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{if $j$ is an output neuron} \\\\ \\bigl( \\sum_{l \\in L} \\delta_l \\cdot w_{jl} \\bigr) \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{otherwise}\\end{cases}\n",
    "\\end{align}\n",
    "\n",
    "Where:\n",
    "\n",
    "$o_i$ is the output of neuron $i$\n",
    "\n",
    "$E(\\cdot)$ is the loss function\n",
    "\n",
    "$o_j$ is the output of neuron $j$\n",
    "\n",
    "$t_j$ is the target for output neuron $j$\n",
    "\n",
    "$h_j(\\cdot)$ is the activation function of neuron $j$\n",
    "\n",
    "$s_j$ is the weight sum part of neuron $j$\n",
    "\n",
    "$L$ is the set of neurons that is connected to the output of neuron $j$\n",
    "\n",
    "$w_{jl}$ is the weight of the connection from neuron $j$ to neuron $l$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example:\n",
    "\n",
    "To give a feel of what the backpropagation partial derivative formula looks like when expended out, we will use an example.\n",
    "\n",
    "A small network with 1 hidden layer and 2 neurons in all the layers:\n",
    "\n",
    "![labeled_neural_network.png](img/labeled_neural_network.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us start with the partial derivative of $w_{h2o2}$:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial w_{h2o2}} = \\delta_{o2} \\cdot o_{h2}$$\n",
    "\n",
    "with $\\delta_{o2}$ being:\n",
    "\n",
    "$$\\delta_{o2} = \\frac{\\partial E(t_{o2}, o_{o2})}{\\partial o_{o2}} \\cdot \\frac{\\partial h_{o2}(s_{o2})}{\\partial s_{o2}}$$\n",
    "\n",
    "Together:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial w_{h2o2}} = \\frac{\\partial E(t_{o2}, o_{o2})}{\\partial o_{o2}} \\cdot \\frac{\\partial h_{o2}(s_{o2})}{\\partial s_{o2}} \\cdot o_{h2}$$\n",
    "\n",
    "Each of the term:\n",
    "\n",
    "$\\frac{\\partial E(t_{o2}, o_{o2})}{\\partial o_{o2}}$ = partial derivative loss function with respect to output of neuron $o2$.\n",
    "\n",
    "$\\frac{\\partial h_{o2}(s_{o2})}{\\partial s_{o2}}$ = derivative of activation function of neuron $o2$.\n",
    "\n",
    "$o_{h2}$ = output of the neuron $h2$. Which is the input corresponding to the weight $w_{h2o2}$.\n",
    "\n",
    "One thing to note here is that we would need the output of neurons $o2$ and $h2$ to calculate $\\frac{\\partial E}{\\partial w_{h2o2}}$. To calculate $\\frac{\\partial E}{\\partial w_{h1o1}}$, we will need output of neuron $h1$ and $01$, etc.\n",
    "\n",
    "So it would be time saving if we just save the input/output of all the neurons. Rather than recalculating it every time we need it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $w_{i2h2}$ the partial derivative would be more complicated:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial w_{i2h2}} = \\delta_{h2} \\cdot o_{i2}$$\n",
    "\n",
    "with $\\delta_{h2}$ being:\n",
    "\n",
    "$$\\delta_{h2} = \\bigl( \\sum_{l \\in L} \\delta_l \\cdot w_{h2l} \\bigr) \\cdot \\frac{\\partial h_{h2}(s_{h2})}{\\partial s_{h2}}$$\n",
    "\n",
    "The set $L$ contains all the neurons connected to the output of neuron $h2$, therefore: $L=\\{o1, o2 \\}$\n",
    "\n",
    "Expanding it out:\n",
    "\n",
    "$$\\delta_{h2} = \\bigl( \\delta_{o1} \\cdot w_{h2o1} +  \\delta_{o2} \\cdot w_{h2o2} \\bigr) \\cdot \\frac{\\partial h_{h2}(s_{h2})}{\\partial s_{h2}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we need $\\delta_{o1}$ and $\\delta_{o2}$, which is also used to calculate partial derivative of weights for neuron $o1$ and $o2$.\n",
    "\n",
    "In other words, we could calculate the partial derivative of neuron $o1$ and $o2$ first, saving the value of $\\delta_{o1}$ and $\\delta_{o2}$ then used it when calculating neuron $h1$ and $h2$.\n",
    "\n",
    "Calculating the partial derivatives from the output layer first then moving toward the input layer is what give backpropagation its name. Since we are moving backwards from the output layer first toward the input layer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuing on, if we expand out the $\\delta$ then the formula becomes:\n",
    "\n",
    "\n",
    "$$\\delta_{h2} = \\Biggl( \\frac{\\partial E(t_{o1}, o_{o1})}{\\partial o_{o1}} \\cdot \\frac{\\partial h_{o1}(s_{o1})}{\\partial s_{o1}} \\cdot w_{h2o1} +  \\frac{\\partial E(t_{o2}, o_{o2})}{\\partial o_{o2}} \\cdot \\frac{\\partial h_{o2}(s_{o2})}{\\partial s_{o2}} \\cdot w_{h2o2} \\Biggr) \\cdot \\frac{\\partial h_{h2}(s_{h2})}{\\partial s_{h2}}$$\n",
    "\n",
    "And the whole thing becomes:\n",
    "\\begin{align}\n",
    "\\frac{\\partial E}{\\partial w_{i2h2}} & = \\Biggl( \\frac{\\partial E(t_{o1}, o_{o1})}{\\partial o_{o1}} \\cdot \\frac{\\partial h_{o1}(s_{o1})}{\\partial s_{o1}} \\cdot w_{h2o1} +  \\frac{\\partial E(t_{o2}, o_{o2})}{\\partial o_{o2}} \\cdot \\frac{\\partial h_{o2}(s_{o2})}{\\partial s_{o2}} \\cdot w_{h2o2} \\Biggr) \\cdot \\frac{\\partial h_{h2}(s_{h2})}{\\partial s_{h2}} \\cdot o_{i2} \\\\\n",
    "& = \\frac{\\partial E(t_{o1}, o_{o1})}{\\partial o_{o1}} \\cdot \\frac{\\partial h_{o1}(s_{o1})}{\\partial s_{o1}} \\cdot w_{h2o1} \\cdot \\frac{\\partial h_{h2}(s_{h2})}{\\partial s_{h2}} \\cdot o_{i2} +  \\frac{\\partial E(t_{o2}, o_{o2})}{\\partial o_{o2}} \\cdot \\frac{\\partial h_{o2}(s_{o2})}{\\partial s_{o2}} \\cdot w_{h2o2} \\cdot \\frac{\\partial h_{h2}(s_{h2})}{\\partial s_{h2}} \\cdot o_{i2} \n",
    "\\end{align}\n",
    "\n",
    "There are two large terms in the formula separated by a $+$ sign.\n",
    "\n",
    "Note each of the large term contains two derivatives for activation functions. \n",
    "\n",
    "This is because the neuron $h2$ is at the second last layer. The partial derivative for weights of neurons $n$ layers away from output layer would have $n + 1$ derivatives of activation function in each term in a fully connected network.\n",
    "\n",
    "This is important to remember when we learn about the vanishing gradient problem and exploding gradient problem next lesson. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recap:\n",
    "\n",
    "How Backpropagation works:\n",
    "\n",
    "With a single case of training data calculate the network output.\n",
    "\n",
    "Calculate the difference of network output and target with loss function.\n",
    "\n",
    "Calculate the partial derivatives of the loss with respect to all the weights.\n",
    "\n",
    "Update each weight according to the update rule:\n",
    "\n",
    "$$\\Delta w_{ij} = - \\eta \\frac{\\partial E}{\\partial w_{ij}}$$\n",
    "\n",
    "Repeat."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Backpropagation algorithm\n",
    "\n",
    "Backpropagation gradient calculation is consist of two phases: the forward pass, and the backward pass.\n",
    "\n",
    "During the forward pass, the input/output of the neurons are saved in a cache. It will be used when calculating the gradients during the backward pass.\n",
    "\n",
    "During the backward pass, the partial derivatives of each weight are calculated starting from neurons in the output layer and working backwards. While computing the partial derivatives, the $\\delta$ are saved. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Full algorithm:\n",
    "\n",
    "With a neural network, a loss function, an input vector, and a target vector:\n",
    "\n",
    "1. Compute the network output from the neural network and input vector. Store the output of each neuron in a cache.\n",
    "\n",
    "2. Calculate the partial derivatives of loss with respect to each weight connected to the output neurons, using variables in the cache. Store the $\\delta$ of output neurons in the cache.\n",
    "\n",
    "3. Calculate the partial derivatives of loss with respect to each weight connected to the neurons of the next layer closer to input layer, using variables in the cache. Store the $\\delta$ of these neurons in the cache.\n",
    "\n",
    "4. Repeat step 3 until all partial derivatives are calculated.\n",
    "\n",
    "5. Use the update rule to update all the weights."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "We will use the same dataset as the last lesson.\n",
    "\n",
    "Dataset: Iris flower classification dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = os.path.join('data', 'iris.data')\n",
    "dataset = pd.read_csv(data_path, names=['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)', 'class'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is what some of the data look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sepal length (cm)</th>\n",
       "      <th>sepal width (cm)</th>\n",
       "      <th>petal length (cm)</th>\n",
       "      <th>petal width (cm)</th>\n",
       "      <th>class</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>99</th>\n",
       "      <td>5.7</td>\n",
       "      <td>2.8</td>\n",
       "      <td>4.1</td>\n",
       "      <td>1.3</td>\n",
       "      <td>Iris-versicolor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>100</th>\n",
       "      <td>6.3</td>\n",
       "      <td>3.3</td>\n",
       "      <td>6.0</td>\n",
       "      <td>2.5</td>\n",
       "      <td>Iris-virginica</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     sepal length (cm)  sepal width (cm)  petal length (cm)  petal width (cm)  \\\n",
       "99                 5.7               2.8                4.1               1.3   \n",
       "100                6.3               3.3                6.0               2.5   \n",
       "\n",
       "               class  \n",
       "99   Iris-versicolor  \n",
       "100   Iris-virginica  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset[99:101]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"class\" column is the type of flower.\n",
    "\n",
    "The other four column are the length and width of sepal and petal, these are the four features given in the dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Out of the 4 features given we are only going to use the petal length and petal width.\n",
    "\n",
    "And we will only be classifying two type of flowers (Iris-setosa and Iris-versicolor)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of our network is to determine which class a data point belongs to.\n",
    "\n",
    "Given set of petal length and petal width, the network would output the class of the flower."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So our training data should be input and target pairs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "training = {}\n",
    "training['input'] = dataset.loc[:99, ['petal length (cm)', 'petal width (cm)']].to_numpy(dtype=np.float_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A input vector would look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.7, 1.4])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "training['input'][50]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to convert the class to an array.\n",
    "\n",
    "Each element of array represents a type of flower. The element corresponding to the correct class would be 1, and other will have a value of 0.\n",
    "\n",
    "For example, \"Iris-setosa\" be may represented with $(1, 0, 0)$.\n",
    "\n",
    "However, we are only using two classes so it would be $(1, 0)$.\n",
    "\n",
    "This type of representation is called one-hot vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_class = {0: \"Iris-setosa\", 1: \"Iris-versicolor\"}\n",
    "to_index = {\"Iris-setosa\": 0, \"Iris-versicolor\": 1}\n",
    "\n",
    "def to_one_hot(num_class, indexes):\n",
    "    \"\"\"Outputs a matrix of one hot vectors.\"\"\"\n",
    "    # Create vectors of zeros.\n",
    "    n = len(indexes)\n",
    "    one_hot = np.zeros((n, num_class), dtype=np.float_)\n",
    "    # Set correct class to 1.\n",
    "    one_hot[np.arange(n), indexes] = 1\n",
    "    return one_hot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "training['target'] = to_one_hot(2, [to_index[t] for t in dataset.loc[:99, 'class']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A target vector would look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 1.])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "training['target'][50]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we are hoping is the network would learn a decision boundary like the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAFgCAYAAADEuRy2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd5hU9fX48ffZXmEXFkSaCIgFUUCkioAYYkGxo9gbaoIlmhhDEmP52kI0MTYsP40dVERsgBVQsQAqIqKIgEivy/Z+fn98Zt2ZrcOys3dm57yeZ57d+cy9M2cROXvvPfccUVWMMcaYaBHjdQDGGGNMc7LEZ4wxJqpY4jPGGBNVLPEZY4yJKpb4jDHGRJU4rwNoSscdd5zOmTPH6zCMMdFDvA7A7LkWdcS3fft2r0MwxhgT5lpU4jPGGGMaYonPGGNMVLHEZ4wxJqpY4jPGGBNVLPEZY4yJKpb4jDHGRBVLfMYYY6KKJT5jjDFRxRKfMcaYqGKJzxhjTFSxxGeMaRkqKryOwESIFtWk2hgThfK3wY/vwZr50Ps06DwAUtp4HZUJY5b4jDGRq2AXvH41/DDbPV/6Ihz9Jxh+A8QnexubCVt2qtMYE7lK86uSXqVPH4Ki3d7EYyKCJT5jTMsi9s+aqZ/9DTHGRK6EVHddz9/wGyApw5t4TESwa3zGmMiVnAknTIE+Z8Daj+GgsdD+YIhP8joyE8Ys8RljIltqFhx0onsYE4SQJT4ReRIYC2xV1UNref1PwLl+cRwMtFPVnSKyFsgFyoEyVR0QqjiNMcZEl1Be4/sfcFxdL6rqFFXtq6p9gb8A81V1p98mo3yvW9IzxhjTZEKW+FR1AbCzwQ2dc4AXQxWLMcYYU8nzqk4RScEdGc7wW1bgHRFZIiITvYnMGGNMSxQOxS0nAZ9UO805TFU3ikh74F0R+d53BFmDLzFOBOjatWvoozXGGBPRPD/iA86m2mlOVd3o+7oVmAkMrGtnVX1MVQeo6oB27dqFNFBjjDGRz9PEJyKtgRHALL+1VBFJr/weGAN8602ExhhjWppQ3s7wIjASyBKR9cA/gHgAVZ3q2+xU4B1VzffbdR9gpohUxveCqs4JVZzGGGOii6iq1zE0mQEDBujixYu9DsMYE47Ky6Bgh2tsHZ8CyRkQt9cdXqQpQjPNKxyKW4wxJvS2rYBnT3Xz+xLS4IwnofuIpkh+JsKEQ3GLMcaEVt42mHGpS3oAJXnueeEub+MynrDEZ4xp+bQMtv0QuFacC6WF3sRjPGWJzxjT8sUmQtfBgWtp+0B8qjfxGE9Z4jPGtHwpbeC0J2C/oe55+4Ph/JlusoOJOlbcYoyJDhldYPzzUF4CMbGQag0vopUlPmNM9Ehp43UEJgzYqU5jjDFRxRKfMcaYqGKJzxhjTFSxxGeMMSaqWOIzxhgTVSzxGWOMiSqW+IwxxkQVS3zGmPBQUghFOQ1vV1YMhdlQUdE0n1tR4d6vrLhp3s+EPbuB3RjjrfJy2L0OPrwL8rfAoCuh6xA3L6+6nI3w8X9g63I49Aw4+GRIbdv4z87fASteh29fgX0OhWHXQquOjX8/ExEs8RljvFWwFR49Gop9R3ur57nWYgePDdwubys8PRZ2/OSer/3YJcKj/whxiXv+uWXF8PkjsGBK1futehcungNp7Rv945jwZ6c6jTHeWr+kKulV+uyhmrPyinZXJb1Ki/+fO03ZGIXZsOSpwLUdP0FRI9/PRAxLfMYYb9XWPzOlLcTEB67VdlSXnAkijftcEUiq5XSqTWRv8SzxGWO81fYA6Ni/6nl8Chzzd0hMC9wuMR0OO7vqucTAcXc3fspCajs4/h73PpUOGw8JaXXvY1oEUVWvY2gyAwYM0MWLF3sdhjFmT+Vtha0r3Nf9hrqkFJdQc7v8HbBrDWxfCV2Hunl61RPknijOg/ztsG4hZPWCzP33tFimkYebxktW3GKM8V5a++AKSlLbukfnAU3zuYlp7tGmW9O8n4kIdqrTGGNMVLHEZ4wxJqpY4jPGGBNVLPEZY4yJKpb4jDHGRBVLfMYYY6KKJT5jjDFRxRKfMcaYqBKyG9hF5ElgLLBVVQ+t5fWRwCxgjW/pVVW9zffaccD9QCzwhKreHao4jTFNJG8r5G6G2HjXeSU1K/h9d2+A4lzXrLpVR0jtACW5ULQLCnZC686Q0s63thvytkBGF0huW3uHl4KdUJIHu36GzP1cG7LaeoKaqBTKzi3/Ax4Enqlnm49UNWD2iIjEAg8BvwHWA4tE5HVV/S5UgRpj9lLuZvjfibBjlXveZZAbLZQWRB/N3Rtg9o3w/ZvueXoHNxpo5Tsw50a3lpwJE+fD0mkw7063lpAKF70FHfsFvl9JAfz0Psy8AirKISYWTn0UDjwRElKa5uc1ES1kpzpVdQGwsxG7DgRWqepqVS0BpgHjmjQ4Y0zTqaiAL5+tSnoAv3wOv3wa3P45G6qSHrgkOu9u6HlM1VrhLnc0OP+uqrWSfHjjOtdr01/hDnjrjy7pgfv61h+hsDH/HJmWyOtrfENEZKmIzBaR3r61TsAvftus963VSkQmishiEVm8bdu2UMZqjKlNRRlsW1Fzfev3we2/c03NtV2rA6cmxMRBUQ5Ub6q/a637/IB4ymvO1CvKrkqEJup5mfi+BPZT1cOBB4DXfOu1dTuvc4SEqj6mqgNUdUC7do0cT2KMaby4BOh/QeCaCBwS5ImaroPddUF/h54B2RuqnleUQXJGzUbWvU+BxFaBa7EJNZtYdx7g1o3Bw8Snqjmqmuf7/m0gXkSycEd4Xfw27Qxs9CBEY0yw9j0cxj0EWQfAPofChJcgvWNw+8anwAVvuOSU2Q1G/Q0OORnS2rgRRa27wPAbIH1fuGg29DjGFbsMusptW/26XauOcMb/XFJs1dEl4DP+B632beIf2kSqkM7jE5FuwJt1VHV2ALaoqorIQOAVYD9cJedKYDSwAVgETFDV5Q19ns3jM8ZDFeVQsMMd7TVmOOzuX9z1wtQsV7gCULALykvc0V7lBPbCbCgrgqTWEJ9c9/sV7HDXAUNb0Wnz+CJQKG9neBEYCWSJyHrgH0A8gKpOBc4ArhKRMqAQOFtdFi4TkUnAXFwSfDKYpGeM8VhMbHAz9erSukvNtZTMmmvJGcG9X0pb9zCmGpvAbowxjWdHfBHI66pOY4wxpllZ4jPGGBNVLPEZY4yJKpb4jDHGRBVLfMYYY6KKJT5jjDFRJZTTGYwxkaas2PW1jE0M/n65xti9AVBI26dmuzJ/BbugogSSMqvGD+VtdXEmZUJSmlsryoHSQkhMtwkMpkGW+IwxTv42WPggLH8VMveHE/4FbXpAbGzTfUbBTtj4Fbx/m5uXN+hKX3uyfQK3Ky+F7T/C7D9B9i9w2HgYNBEKsmHOn2H7Suh1PIz8M5QWwOybYPNS6DEajvlrzfczxk/03sBeVgYPPAC/+x0kJoY2MGPCXWkhfHA7fPpQ1VpyJvzuc0hvwiSyfSU8NAi0omrtrGdd8vOXuxkeOMIlx0oj/uxGC33xeNXaJXNh1u9gx09VaweeAKc8Etoj1ip2A3sEit5rfJMnw/XXw7HHwvbtDW9vTEtWtBu+nRG4VrjLJaCm9P3bgUkPYOmLrv+mv52rA5MewLKXocPhNd/TP+kBrJztenkaU4foTXwTJkDnzvDxxzBoEHwf5OwwY1qi2HjI2K/mem29MvdGVq+aa2161Gw2XVuT68xukLclcC02oea4oVadXKNsY+oQvYmvb1/4/HM44ghYvRqGDIH33/c6KmO8kdIWTrzPFYdUGnpNzVl3e6tjP9hvWNXzjP1g8JVVkxf84zni4qrnSRnw27tc0U2l+BQ3xeG3d1YlutgENx4pxWZzmrpF7zW+Svn5cN558NprEBcHjzwCl10WmgCNCWflpW6UT/Y6d8SVnBma62Q5m9znlBZC647uCK02BTvd6daCHZDR1SWzwh1QnAu5m1wBTnKmG1JblA05G912SZmQUM+4oqZlh5YRyBIfuBlgf/kL/POf7vmf/gR33w0x0XtAbIwJiiW+CGT/soNLcPfcA4895o76pkyBM85wR4PGGGNaFEt8/i6/HObMgdatYeZMGDECNm70OipjjDFNyBJfdaNHw6efQvfusGSJq/hcutTrqIwxxjQRS3y1Ofhg+OwzGDYM1q93X9980+uojDHGNAFLfHVp1w7ee8/d75efD+PGwf33QwsqBjLGmGhkia8+SUnw3HNw662u8vO662DSJNfuzBhjTESyxNcQEbj5Znj+eUhIgIcfhpNOgpwcryMzxhjTCJb4gjVhAnzwAWRlucrPYcPg55+9jsoYY8wesrFEe2LYMNfmbOxY+PZbGDgQXn/dVX4aEw4KdkHuRtj0DXQ+wo3nSWrd+Pcr3O26oqz9CFp1hqwDICUL8rfCz59AShvY51BIbe/6aK77FOKSoFN/aNURdm+EDYtcY+rOgyC9AxRsh83LoDgPugyEtPauS8u2FW7W3n5D3PsVZcOO1ZD9M3QdAqlZNXt6GtMIlvj2VPfusHAhnHmmK34ZORKefhrOOsvryEy0K86Dzx+F+XdVrZ30Xzj87Jq9MIO1aw08OcYNfgU46EQ49lZ4fJRrHQbQoQ+cMx0eHe7ai4FrRn3+a27f3evdWnoHuOx9eP4s2LrcrSVnwsR58Po1sGa+W4tPgd9/Bm9eD6vec2ux8W4EUacjGvdzGOPHTnU2RkYGvP02TJwIRUUwfjzceadVfBpvFefAx/8KXHv3767fZWPkbob3bqlKeuB6YS6YUpX0wB29bVgS2NezdRc35qgy6VW+31fPQabfFIjCXfDJA4FrMbHuSLEy6YHrIzp3sjsyNGYvWeJrrPh4mDoV7r3XFcD89a9w8cVQXNzwvsaEQkW5SxD+inMb/wtZeUnNpJmQBvm1zK/M3xY4ySExzZ0OrS5vS+AECHDb+a/FJbr5gNUV7HQ/ozF7yRLf3hBxw2xffRVSUtwpzzFjYMcOryMz0SghJXDkD0DvU93onsZI7wgDLg5cW/cZDLqy2uemQo9j3JHfr9t9Cn3Pc0dvlURgwKWwbmHg/oOvgp8/rXqev91dS0xrH7jdwCsguU3jfhZj/Nh0hqby5ZfuNoeNG6FnT3jrLehVy9BNY0Ipbwt8NtUlngPGQP/zax/qGqycTe7a21fPufFBI250R2frF8MXj7o5eSP/4gpetn0Hn/zHFbeM+DOkd4Lste7UqJbD8D9C254usX14h5uwPuw62Pdwd21w3l2uuGXQFdBtuCtumT8Fdq2GfhfAAb9xxTThxaYzRCBLfE1p/XqX/L7+GjIz3ZHgyJHexWOiU1kJlOa7BBXTRPVrOZtcQvOfyJ63zRWd+F/by98OEhu4XcFOQN1w2Uol+e60ZZLf6dHSAneq1r8KtbTQXWMMxVzApmGJLwJZ4mtqeXnunr833nDXAR991F37M8a0RJb4IlDIrvGJyJMislVEvq3j9XNF5BvfY6GIHO732loRWSYiX4uIx5lsD6WluZFGf/gDlJbCJZfA5Mmu5ZkxxhjPhbK45X/AcfW8vgYYoaqHAbcDj1V7fZSq9lXVASGKL3RiY+G+++CRR9z3d93lbnkoLPQ6MmOMiXohS3yqugCo86YbVV2oqpW10p8BnUMVi2euvNLd79eqFbzyirvet3mz11EZY0xUC5fbGS4FZvs9V+AdEVkiIhPr21FEJorIYhFZvG3btpAG2ShjxrhOL926wRdfuPZmy5Y1uJsxxpjQ8DzxicgoXOL7s9/yMFXtDxwP/F5Ejq5rf1V9TFUHqOqAdu32omw7lHr3dj0+Bw+Gdetcz885c7yOyhhjopKniU9EDgOeAMap6q93favqRt/XrcBMYKA3ETah9u3ddIfx4yE3F048ER56yOuojDEm6niW+ESkK/AqcL6qrvRbTxWR9MrvgTFArZWhESc5GV54Af7+d1flOWkSXHstlFsbJmOMaS4hm84gIi8CI4EsEVkP/AOIB1DVqcDNQFvgYREBKPNVcO4DzPStxQEvqGrLOS8YEwO33QYHHACXXQb//S/89BO8+CKkpze8vzGhVLgbyorc3Wmp7V2bsaJcd3O54Dq0+Lch81daCEW+Ac3JmRCXUPt25aVVzaYT012rtWCpur6gWuFuqA/fG9tNGLMb2L300Udw6qmut+dhh7mb3rt29ToqE61yt8Cbf4CVsyFjPzj1UWizP8z9Gyyf4Wb7jf2P6weaWK3/Z8EOWPggfD4VYhPgmL9CnzNdAvRXmA3fvQ7v/wNKCuDIS2H4DYFdXepSVgKbv4GZE2HnauhxLIx70I078o7dwB6BPC9uiWrDh8Nnn7ment984yo+Fy3yOioTjUry4b1/wA9vuaOpXWtg6Yvw6UOwbDpUlEHOBph2NhTVMuZozcfw8X3uyLAoG97+E+xcW3O7nA3wxtXuiK+syL3/yrnBTZAo3AnPngI7fnLbr3oXZv85cESSMUGwxOe1nj3h009h1Ch3j9+IEa7HpzHNqTg3cP4dQIdD4ce5gWsV5bDlu8C1smJYXsvf2er7Avz0Yc2172a5hNmQwp01k9xPH7ikbcwesMQXDtq0cbc3XHKJ6+5y+ulwzz022NY0n7gkN0nd3861sG/fmtu27R74PDYB9htac7vOtRRjd+pfc63LIIhNajjGpIyaTbf36e0+35g9YIkvXCQkwBNPwN13u+c33eSKX0pKvI3LRIfkDDjxPndtD1wy2edgGPVXaHegW4uJdc9TsgL3FYHep8H+frfbHjYe9j2s5udkHQj9L3T7AHQZDP3Oc639GpLYCk5+wA2qBTcm6eQHw3FUkQlzVtwSjmbMgPPPd0d/o0a555mZDe9nzN6orJgsyXfJJam1GzKbt9UVosQluCrM6hPUK+XvcOOQJMbtV72wpVJhtpvFV1HutkvNqn272pQUuOnsZYVu35R2rlLaO1bcEoEs8YWrRYvg5JPddb8DD4Q333TXA40x4cQSXwSyU53h6sgjXZuzPn3ghx9cu7OPPvI6KmOMiXiW+MJZ167wySdwwgnuXr9jj4XnnvM6KmOMiWhBJT4RyRSR3iLSXUQsWTan9HSYNQuuucYVupx/Ptx8s1V8GmNMI9WZxESktYhMFpFluHl5jwIvAT+LyMu+qQqmOcTFwf33wwMPuAv5t98OEyZAUZHXkRljTMSp7+jtFeAXYLiqHqiqR/nG/3QB7gbGicilzRKlcSZNckUu6ekwbRoccwxs3ep1VMYYE1GsqjMSLVsGY8e62X7durlk2Lu311EZE42sqjMCBXuN7zAROVlETqt8hDowU48+fVzF55FHwtq1MHQovPOO11EZY0xEaDDxiciTwJPA6cBJvsfYEMdlGtKhA8ybB2ecATk5rvJz6lSvozLGmLAXzDy+wap6SMgjMXsuJQWmT4e//Q3uuguuugpWroQpU4JrAWVCo7wU8jbD0pcgNh76nA6p+0BsyMZf7rnCXbD9R1jxBnQ5EroOdbHu+hm+neF6YHYfBWntvI7UmCYXzP+Jn4rIIar6XcObmmYXEwN33ukG215xBfz737BqlZv0npbmdXTRKWcjPDKkamrAR/fC7xa63pLhoKwEvnkZZv+pau3AE2DUZHh0eNVa5yPhnGl71lLMmAgQzDW+p3HJ7wcR+UZElonIN6EOzOyhiy921/kyM91A2+HDYf16r6OKPqqw6InAUTlF2S7RhIvCXTD/rsC1H952R6r+1i+qmpRuTAsSzBHfk8D5wDKgIrThmL0ycqQbbHviifD1126w7RtvQP9aRsGY0KmeQMANcg0XgmsQXZ3W8r93bWvGRLhgjvjWqerrqrpGVX+ufIQ8MtM4vXq55Hf00bBxozvymzXL66iihwgMuqJqdA64KQKHjfcupuqSMmDYtYFr+x/trvH52+dQSG3bfHEZ00wavI9PRB4GMoA3gOLKdVUNuzHhUXMfXzCKi901v6efdv8YT5kC119fNQfNhE5pMeRugM8fg5h4GHgZtOpYM7F4qWAnbPgSvn3ZFbYcdIKLdet38OUzsO/hcOjpkNbe60jDnf0PFYGCSXxP1bKsqnpJaEJqPEt81ai6as+//tU9nzgRHnwQ4sPoH+CWrKLC/aIRzr9sVFTUnGdX25qpSxj/xzV1sc4t0eCll+DCC11vz2OPhZdfhowMr6MypiWwxBeBgrmB/WkRyfB7num7qd1EirPOcje7t28P773nOr2sXu11VMYY44lgzmccpqrZlU9UdRfQL3QhmZAYNMi1OevdG1ascM8XLvQ6KmOMaXbBJL4YEcmsfCIibQjuNggTbrp1c4Ntx4yB7dvddIcXX/Q6KmOMaVbBJL57gYUicruI3AYsBP4Z2rBMyLRuDW+95dqbFRe7uX633mqDbY0xUaPBxKeqz+AaVG8BtgGnqeqzoQ7MhFBcHDz0EPznP67i8JZb3GR3G2xrjIkCdVZ1ikiaqubVu3MQ2zQnq+pshDfegHPOgfx8GDYMZs6EdtaY2JggWVVnBKrviG+WiNwrIkeLSGrlooh0F5FLRWQucFzoQzQhddJJ8PHH0Lmzu/43eDB8/73XURljTMjUmfhUdTTwPnAFsFxEckRkB/Ac0AG4UFVfqe/NReRJEdkqIt/W8bqIyH9FZJWvAXZ/v9cuFJEffY8LG/PDmSD17esqPo84wt3mMGQIvP++11FFr8JsyP4Fdq5xkx7qUlIAu9fDjp/c17pUlEPuFjdyKHczlO9h39D87S6enE3uM8F1ftn9i4uvOL/+/Y0JMyG9gV1EjgbygGdU9dBaXj8BuBo4ARgE3K+qg3yVo4uBAYACS4AjfLdS1MlOde6l/Hw47zx47TV3HfCRR+Cyy7yOKrrkbYWl0+HD/4OyIjca6MynoXW1kUZFufDzRzDzKjf9oW0PN0Ioq1fgdqqwZTm8cBbkbHAtyMY/Dx37BzcfMHczTD8X1i+G+GT4zf/BwWPh1YmwZr5rw3bUH2HQREhp03R/DpHDTnVGoJD2JVLVBUB9c03G4ZKiqupnQIaI7Av8FnhXVXf6kt272GnV0EtNhRkz4MYboawMLr/cfV9hHfqbTVE2vPs3l/TAjQZa8E8oygncrjgHXrnUbQ/uqG/WJHdU5i9/G7x0vkt64BLrtAlQsKPhWEoKYN7dLukBlBbCuoWw6EmX9MBNoph/V/1HpsaEGa8b8nUCfvF7vt63Vtd6DSIyUUQWi8jibdu2hSzQqBETA/fcA48/7o76pkyBM85wR4Mm9LbUMu95/eJaEl8ulBYErm1YUnOMUHkp7KzWpSd/W1VirU9JPvzyWeBaVi9Y90nNbTd+1fD7GRMmvE58tZ0m0HrWay6qPqaqA1R1QDurRmw6l10Gc+a4+/5mzoQRI9yYIxNaHfrUXOt2FCRX662amO4e/vYbCjGxgWtxCdDuoMC1Vh0hPqnhWBLTYf+RgWtblkOP0TW37Xxkw+9nTJgIKvGJSKyIdBSRrpWPJvr89UAXv+edgY31rJvmNHq0m+3XvTssWeLanC1d6nVULVtSBox7GJIz3T2WB/wGhl0HiWmB2yVnwtkvQivfiZCO/eHkByC9Q+B2qe1g/HNVya9Nd5jwMqQE8UtifBIMvx56/qbqMw85BfqdC71PA4lxyfGEf9X8XGPCWDBjia4G/oG7gb3yPIqq6mFBfYBIN+DNOopbTgQmUVXc8l9VHegrblkCVFZ5fokrbqnveqEVt4TK9u1wyinudofUVJg2DcaO9TqqlqukAAp3usKUuMS6Z+KVl0LeFle1GRvvjuTqkrcNykvcdqnt9mxUUmG2O60qMZDcxh1FFu12cYq4ZB3MEWTLZMUtESiYxLcKGKSqQVwNr7Hvi8BIIAuXOP8BxAOo6lQREeBBXOFKAXCxqi727XsJMNn3Vneoam1zAQNY4guhoiK49FJ44QV3HfC+++Caa8J71pwxoWf/A0SgYBLfh8BvVHUPb/5pfpb4QkwVbr8d/vEP9/x3v4P773dFMMZEJ0t8EajOf7FE5Hrft6uBeSLyFlBc+bqq3hfi2Ey4EYGbb4aePeHii+Hhh+Gnn2D6dFcEY4wxEaC+4pZ032Md7j66BL+1tHr2My3dhAnwwQeQlQVz57oen2vXeh2VMcYEpc4jPlW9FUBEzlTVl/1fE5EzQx2YCXPDhrk2Z2PHwvLlruLz9dfdV2OMCWPB3M7wlyDXTLTp3t1NcT/2WNi6FUaOhJde8joqY4ypV33X+I7H3WbQSUT+6/dSKyDsC11MM8nIgLffhkmT4LHHYPx4WLUK/vIXq/g0xoSl+o74NuLupSvyfa18vI7rpWmMEx8PU6fCvfe6ZPfXv7ril+Lihvc1xoSMiNwiIn/0Oo5wU981vqXAUhF5XlVLmzEmE4lE4PrroUcPV/zy9NOwZg28+iq0bet1dMYY86s6j/hEZJmIfAMs8c3KC3g0Y4wmkowbBx99BB07woIFbrDtypVeR2VMVBCRC3z/Ri8VkWervXa5iCzyvTZDRFJ862eKyLe+9QW+td4i8oWIfO17vwO8+HlCpb5TnWOBk4A5vse5vsfbQL0DaE2U698fvvgC+vVz1/sGD4Z587yOypgWTUR6A38FjlHVw4Frq23yqqoe6XttBXCpb/1m4Le+9ZN9a1fi5qP2xc1FrWfSceSpbwL7z6r6MzBMVW9U1WW+x03YNT7TkE6d3BHfSSfBrl0wZgw81WDXOWNM4x0DvKKq2wFq6W18qIh8JCLLcAcxvX3rnwD/E5HLgcrxHp8Ck0Xkz8B+qloY+vCbTzC3M6SKyFGVT0RkKJAaupBMi5GW5kYa/eEPUFoKl1wCkyfbYFtjQkOoY3ybz/+ASaraB7gVSAJQ1SuBv+Em4nwtIm1V9QXc0V8hMFdEjgll4M0tmMR3KfCQiKwVkbXAw8AlIY3KtByxsa6h9SOPuO/vusvd8lBQ0PC+xpg98T5wloi0BfBNufGXDmwSkXjcER++7Xqo6ueqejOwHegiIt2B1ar6X1wlf1DTeCJFg92FVXUJcLiItMI1td4d+rBMi3Plle6G9zPPhFdegZ9/dp1eOtgcN2OagqouF5E7gPkiUg58Baz12+TvwOfAz8AyXCIEmOIrXoaVPbUAACAASURBVBFc8lwK3AScJyKlwGbgtmb5IZpJndMZROQ8VX3Or1l1gHBsUm3TGSLA8uWuzdnatdC1K7z5JvSpZeq4MZHBujREoPpOdVZex0uv42HMnuvd2/X4HDwY1q1zPT9nz/Y6KmNMFAlmHl+SqhY1Uzx7xY74IkhhoevuMn26G2z73//C73/vdVTG7Ck74otAwRS3fCsin4jI3SJygojY4DWz95KT3TT3v//dVXlOmgTXXgvl5V5HZoxp4RpMfKraEzgHdzF0LK6N2dehDsxEgZgYuO02eOYZSEhwR33jxkFurteRGWNasAYTn4h0BoYBw4F+wHJgeojjMtHk/PPhvfdcT8+33oKjjnLX/4wxJgSCOdW5DrgOmK2qQ1T1RFW9K8RxmWgzfDh89hn06gXffOMG2i5a5HVUxpgWKJjE1w94BpggIp+KyDMicmlDOxmzx3r2dMlv1CjYvBlGjHDTHYwxDRKRvHpeWxjCz50cqvcOlWCu8S0FngaeAj4ARuBuhDSm6WVmwpw5rr1ZYSGcfjrccw80UH1sjKlJRGIBVHVoCD+m5SU+EVmMa1h6KvA9cLSqdgtxXCaaJSTAE0/A3Xe75zfdBJddBiUl3sZlTBPpdtNbE7rd9Nbabje9VeH7OqGp3ltERorIhyLyAq4o8dejQRHZV0QW+MYNfSsiw2vZv9aRRCJynt/6oyISKyJ3A8m+ted9213ve+9vReQ631qqiLzlG330rYiM963f7BuV9K2IPCYizXJ7SDD38bVT1W3NEczesvv4WqAZM1zxS2GhOwU6Y4Y7KjQmPOzxP9S+JPc4kOK3XABcvvbuE19odCAieaqaJiIjgbeAQ1V1TbXXbgCSVPUO39FgiqrmVnufB4DPVPV5EUnATWzoBvwTOE1VS0XkYd82z1S+t2/fI3DNsAfj/mw+B84DugPHqerlvu1aq+puEWlTOUXCNz/wJVV9o7F/BsEK5lRnRCQ900KdfjrMn+96en74oev4smqV11EZszfuJDDp4Xt+ZxN+xheVSa+aRcDFInIL0Kd60vOpbSTRaOAIYJHvdrbRuGRW3VHATFXNV9U84FXcHQHLgGNF5B4RGe7X83mUiHzuG5V0DFWjkkIqmOIWY7x15JGuzVmfPm6a+6BBbsq7MZGp6x6uN0Z+bYuqugA4GtgAPOub2H6q71Tl1yIyoI6RRAI8rap9fY8DVfWWWj6i1iNgVV2JS5zLgLt8pziTcNN+zvCNSnoc36ikULPEZyJD167wySdwwgmwcyeMHg3PPut1VMY0Rl03qYb85lUR2Q/YqqqPA/8P6K+qM/0S2uI6RhK9D5whIu1979PG914Apb5RRwALgFNEJEVEUnG1IR+JSEegQFWfA/4F9KcqyW0XkTTgjFD//JXqHEskIqfVt6OqWp25aV7p6TBrFtxwg+vycsEF8OOPcOut0DzXxI1pCpOp/Rpfc1RHjgT+5Bs3lAdcUMs246k2kkhVd4rI34B3RCQGKAV+jxtx9BjwjYh8qarnisj/gC987/WEqn4lIr/FjT+q8O17lapmi8jjuKPAtbjTsM2ivrFET9Wzn6pq2A2jteKWKPLgg663Z0UFnH02PPUUJDXLWRJj/DXqNy5fgcuduNOb64DJe1PYYvZMg1WdkcQSX5SZPdtNc8/NhSFD4LXXoH17r6My0cVONUSgBiewA4jIibhqm19/pVbVBifyishxwP24ctgnVPXuaq//Gxjle5oCtFfVDN9r5fjuQQHWqerJwcRqosjxx7vrfmPHwqefuqKXN990M/+MMaYOwdzAPhV3zvdq3G83ZwL71bsTv3YMeAg4HjgEOEdEDvHfRlX/UHlRFXgAV/paqdDvgqslPVO7Pn1cxefAgW6q+9Ch8M47XkdljAljwVR1DlXVC4BdqnorMAToEsR+A4FVqrpaVUuAacC4erY/B3gxiPc1JlCHDjBvHpx5JuTkuMrPqVO9jsoYE6aCSXyFvq8FvpLUUmD/IPbrBPzi93y9b60GX1ns/rheoJWSRGSxiHwmIqfU9SEiMtG33eJt2+xe+6iVnAzTpsFf/uKG2V51FVx/vQ22NcbUEEzie1NEMoApwJe4stNpQexX20XfuippzgZeUVX/f6W6quoAYALwHxHpUduOqvqYqg5Q1QHt2rULIizTYsXEwJ13wpNPQnw8/PvfcOqpkFdn03pjTBQKJvH9U1WzVXUG7treQcD/BbHfegJPiXYGNtax7dlUO82pqht9X1cD83DjkYxp2MUXu+t8mZnwxhtu1t/69V5HZUxIeTWWKFgi8rbvIGpP97tFRP7YlLEEk/g+rfxGVYt9PdY+rWf7SouAA0Rkf1+j07NxXQACiMiBQKb/e4pIpogk+r7Pwk2A/y6IzzTGGTnSzfbr2RO+/tpVfC5Z4nVUxjSrZhpL5P95dd4poKonqGq2lzFUqjPxiUgHX6ftZBHpJyL9fY+R1GywWoOqlgGTgLnAClzX7eUicpuI+FdpngNM08AbCg8GFovIUuBD4G5VtcRn9kyvXi75HX00bNzovr72mtdRGQO3tJ7ALa3XckvrCt/XcBpL9LmI9PZ7Pk9EjvCNFnrSN0boKxEZ53v9IhF5WUTewHV2qfUzRGSt70AGX4/Qb8SNKXrWt7afiLzvW39fRGr0LhWRvr66j29EZKaIZPrFeKeIzAeubfDPqJ7OLRcCFwEDAP+7wnNwzUrDrmWZ3cBualVcDFdcAU8/7VqbTZniCl+szZnZe3v+l8gluVrHEnHL7nAYS/QHIENV/yEi+wLzVbWXiNwJfKeqz/lOWX6BuwR1Ju7y12G+1ma1foaIrMXlk31wt64NU9Xt4htN5Eucr6jq0yJyCXCyqp4ibpJEnqr+S0S+Aa5W1fkichvQSlWvE5F5vth+F8yfVZ1HfKr6tKqOAi5S1VF+j3HhmPSMqVNiomtpdscdbpL7H/8IV14JpaVeR2aiU7iPJXoJl8wAzgJe9n0/BrhJ3FiiebiGJpVHZe9WztUL4jOOwSW47QB++w0BKhP/s7gRR78Skda4hDzft/Q0btJEpem1/Cy1CuYa3yci8v9EZLbvww8RkUuD/QBjwoIITJ4M06e7np6PPebu98sO+SUHY6oL97FEG4AdInIYrnlJZRW/AKf7NRbpqqorqn9ebZ9RLQyh7gr/gHCD/Fkr1foz1yaYxPcU7jpdR9/zlcB1exiQMeHhrLPcze7t28N777lOL6tXex2ViS5hPZbIt+k04EagtapWto6cC1wt4q4RiEitlfa1fUa1Td4HzhKRtr7t2/jWF+KKIAHOBT7238lXWLnL77rk+cB8GiGYxJelqi8BFb4PLwPsrmATuQYNcm3OeveGFSvc84WeV3ub6DEZd03PX3OOJfpaRL4CTsf1Uq7NK7gk9JLf2u1APG4E0be+53v8Gaq6HLgDmO8rYLzP99I1uFOk3+CSWm1FKhfixht9A/QFGuwZXZsGpzP4LhqejjuH219EBgP3qOqIxnxgKFlxi9kju3e76Q5z51ZdBzznHK+jMpGlcRVSrsAlYCzR3hS2mD0TTOLrj2sgfSjwLdAONyr+m9CHt2cs8Zk9Vlbm5vo9/LB7fsstcPPNVvFpgmV/USJQUPP4fDcEHoj7j/yDqoZlOZwlPtMoqm6i+x/+4L4/91x44gkbbGuCYYkvAgUzligJd+71duBW4Pe+NWNaBhF31DdrFqSmwvPPw7HHgjU9N6ZFCqa45RncENoHgAdxs/WeDWVQxnjipJPg44+hc2c34HbwYPj+e6+jMsY0sWAS34Gqeqmqfuh7TAR6hTowYzzRt6+r+DziCHebw+DB8P77XkdljGlCwSS+r3yVnACIyCDgk9CFZIzHOnaE+fPhlFNc5edxx8Hjj3sdlTGmiQST+AYBC30NRtfipiiMEJFlvnspjGl5UlNhxgy48UZX+Tlxovu+osLryIwxe6nB8Q3AcSGPwphwFBMD99wDBxzgJrpPmQKrVsGzz7rEaIyJSEHdzhAp7HYGEzLvvw+nn+5OfR5xBLz+ujslaqKd3c4QgYI51WlMDeUVytacIuav3MbSX7LZkVfsdUihNXq0m+3XvbsbaDtoECxd6nVUxphGCOZUpzE1bMgu5KQHPmZ3oetlcGS3TB457wiy0hI9jiyEDjrIVXyecoq73WHYMJg2DcaO9ToyY8wesCM+s8cKS8q5/72VvyY9gEVrd/HT1jwPo2omWVnutOe550J+PowbB/ff7zq+GGMigiU+s8dKyivYuLuoxvrmnJprLVJioitwufVWV+V53XUwaZKr/jTGhD1LfGaPtU6O54LB+wWsJcbFcGS3NnXs0QKJuGbWL7wACQmuyfXYsa74xRgT1izxmUYZ0qMt/zm7L4d3bs3IA9sx6/fDyEpL8Dqs5nfOOfDBB+4U6Ny57rrf2rVeR2WMqYfdzmAaTVXZVVBKfKyQnhTvdTjeWr3aHfGtWOGmu8+a5dqdmZbObmeIQHbEZxpNRGiTmmBJD9xtDgsXuqkOW7fCyJEwfbrXURljamGJz5imkpEBb7/t2psVF8PZZ8Mdd1jFpzFhxhKfMU0pPh6mToV773UFMH/7G1x0kUuExpiwYInPmKYmAtdfDzNnQkoKPPMMjBkDO3Z4HZkxBkt8xoTOuHHw0Ueup+eCBa7YZeVKr6MyJupZ4jMmlPr3hy++gH793GSHwYNh3jyvozImqlniMybUOnVyR3wnnwy7drnTnk895XVUxkQtS3zGNIe0NHj1VXftr7QULrkEJk+2wbbGeCCkiU9EjhORH0RklYjcVMvrF4nINhH52ve4zO+1C0XkR9/jwlDGaZpOaVkFW3OL2JZbRHm5/aMeIDbWVXtOneq+v+suGD8eCgq8jsyYqBKyzi0iEgusBH4DrAcWAeeo6nd+21wEDFDVSdX2bQMsBgYACiwBjlDVXfV9pnVu8dau/BJeXLSOJz5aQ1JcDH8+/iBGHdSeVnaDe03vvANnngk5OXDkkW6wbYcOXkdl9px1bolAoTziGwisUtXVqloCTAPGBbnvb4F3VXWnL9m9CxwXojhNE/l8zQ7+OecHduaXsHF3EddO+5qN2YVehxWexoxxnV66dYNFi9xg22XLvI7KmKgQysTXCfjF7/l631p1p4vINyLyioh02cN9TZgoKi3nta831lh/77utHkQTIXr3doNtBw+Gdetcg+vZs72OypgWL5SJr7ZTANXPq74BdFPVw4D3gKf3YF+3ochEEVksIou3bdvW6GDN3omPjeHwzhk11vt0auVBNBGkfXs33WH8eMjNdY2uH3zQ66iMadFCmfjWA138nncGAg4JVHWHqlb2cnocOCLYff3e4zFVHaCqA9q1a9ckgZs9FxsjnHFEZw71S3RjDtmHQzu19jCqCJGc7Ob6/f3vrsrz6qvhmmtssK0xIRLK4pY4XHHLaGADrrhlgqou99tmX1Xd5Pv+VODPqjrYV9yyBOjv2/RLXHHLzvo+04pbvLc9r5i8ojJiY4TUxDjapEbhjL698eyzcNllUFICJ5wA06ZBerrXUZm6WXFLBIoL1RurapmITALmArHAk6q6XERuAxar6uvANSJyMlAG7AQu8u27U0RuxyVLgNsaSnomPGSlJZKVluh1GJHr/PNdwcupp7pJD0cdBW+8AV27eh2ZMS2GDaI1JhytWgUnnuh6e3bo4G53OPJIr6MyNdkRXwSyzi3GhKOePeGzz2DUKNi8GUaMcJ1fjDF7zRKfMeEqMxPmzIFLL4XCQjj9dLjnHhtsa8xessRnTDhLSIDHH3cJD+Cmm6qKX4wxjWKJz5hwJwI33ggzZrhbH558Eo47zk16MMbsMUt8xkSK006D+fNdscuHH7qOL6tWeR2VMRHHEp8xkeTII12bsz59XMXnoEFuyrsxJmiW+IyJNF27wiefuBvcd+6E0aPdje/GmKCE7AZ2E/52F5ayYVcB81du4/AuGRzUoVWtnVZ25BWTXVDK7G8307VtMgO7tSEzNYFtucXMXb6ZrLREhvbIIiMlnm25xby3Ygsp8bEc3asd7Vsl1frZW3OK+OjH7eQVlzHmkH3ISk8gPjY21D9yy5GeDrNmwQ03wH//Cxdc4I4Ab70VYuz3WWPqYzewR6nS8gpmLFnPTa9WjcI5rX8nbh57CBkpgcnv89U7mPDE55RXuL8rfTq15pFz+3PMvfMp8Q2b7Z6VynOXDeI3980nv6QcgA6tknh90rAayW9rbhHjHvyETbuLAEhJiGXOdUfTtU1KyH7eFu3BB+Haa12fz/Hj4amnXBGMaQ52A3sEsl8No1R2QQn3zPk+YO3VLzdQ4EtalTbvLuRf7/zwa9IDWLZhN2t35JOaWHWEtnp7Pl+t20WH1lVJbnNOEfNX1pyY8dHK7b8mPYCCknIenf8TpWU2sb1RJk2CN990R4HTp8Mxx8BWGwdlTF0s8UUpBYpKayYa/wRX+bywtLzGdkWlFcTHBv71yS8pJyGu2lpxzQkD+SW1rBWXUdGCzj40u+OPd9f9unZ1HV8GDYLlyxvez5goZIkvSqUnxXPR0G4BawO6ZZKaGHjZt0OrRC49qnvA2r6tk+i1Txpbc4t/XctIiWdoj7as3JL761pKQixjeneo8dnHHrwPqQlVR4sxApcf3Z3EeLvGt1f69HEVnwMHwtq1MHQovPOO11EZE3bsGl8U25lfwsc/buOtZZs4Yr9MTuvXmaz0mpMVtuQUsWJTDtMX/ULnzBQuGtaNjOQ4Vm7J4/99vIb26S45ZqbE8/POAh5b8BMpCXFccXQP9s1IrFG0UlpewabdRTy24CfyisqYeHQP9mubUiPpmkYqLIQLL4SXX4bYWHcN8MorvY6qpbJrfBHIEp+h0HeKMjam/v+HcwpLSYgTkuKrElRxWTmxIsT5nfYsKStHEOLj6j+hUFpWgaIkxNmRXpOrqHCDbe+80z3/wx9gyhSXCE1TssQXgexUpyE5IbbBpAfQKjk+IOkBJMbFBiQ9gIS42AaTHkB8XIwlvVCJiYE77nAVnvHx8O9/uxl/eXleR2aM5yzxGdOSXXQRvPuum/TwxhswfDisX+91VMZ4yhKfMS3diBGu0rNnT/j6a1fxuWSJ11EZ4xlLfMZEg169XPI7+mjYuNF9fe01r6MyxhOW+IyJFm3butsbLrwQCgrctId777XBtibqWOIzJpokJrqClzvucAnvj3+EK66A0lKvIzOm2VjiMybaiMDkya69WVKSm/B+/PGQne11ZMY0C0t8xkSrs86CefOgfXt4/30YMgRWr/Y6KmNCzhJflNu0u5BfdhawMbuQUt+khY3Zbm3DrgIKfL02swtK2LS7kK05RZTU00w6t6iUzbuL2JJTREEtPTlNmBk0yLU5690bvv/ePV+40OuojAkp6xEVxdZsz+d3zy9hxaZcOrRK4qFz+7FPqySufuErvvolm7apCdx9+mEc1rk117/0NZ+s2kGr5DhuH3coow/ah7SkwL8+O/KLuevt75n51QbiYoSrRvbgwqHdyEypOePPhJFu3VyD6/HjYe5cN93hqafgnHO8jsyYkLAjvii1eXcRN7y0lBWbXFPpzTlFfLFmJ7e/8R1f/eKu9ezIL+H3z39JdkEpn6zaAUBOYRnXTf+a7MKSgPdTVd5fsZVXlqynvEIpLqvgP+/9yJpt+c37g5nGad3ajTb63e+guBgmTHBDba3i07RAlviilKJ8uW5XwNr+WWl8sXZnwFpJeQVbcopI9pucoAo/bQtsfVVcVsH7K2rOgPvkp+1NGLUJqbg419D6/vtdy7NbboHzz4eiogZ3NSaSWOKLUgIc3rl1wNraHfn075oZsBYfK7RvlUhRWeBMvu5ZaQHPE+NiGNErq8bnDO7etmkCNs1DBK65BmbNgtRUeP55OPZY2FZzoLAxkcoSX5Tq0DqZe8/qy/5ZqQC0TU2gb5cMbjm5N4fs2wqAVslx/PusvmQkx9Ovi0uIKQmx3H1aHzJS4gPeT0T4be8OjD1sX0Rcwrzi6O70bBeYIE2EGDvWXffr3Nl9HTzYFb8Y0wLYWKIotym7kJJyN029TWoCSfGxbMguoKxciYsRWifHk5YUz868EorKyon1rSXVMTR2d2EJBSVuLFFaUhxpNmMvsm3aBCed5Hp7tm4NM2bA6NFeRxVObCxRBLLEZ4ypX34+nHee6+0ZFwcPPwyXX+51VOHCEl8ECumpThE5TkR+EJFVInJTLa9fLyLficg3IvK+iOzn91q5iHzte7weyjiNMfVITXVHejfeCGVlMHGi+76i7vs5jQlnIUt8IhILPAQcDxwCnCMih1Tb7CtggKoeBrwC/NPvtUJV7et7nByqOI0xQYiJgXvuce3N4uLcNPfTT3dHg8ZEmFAe8Q0EVqnqalUtAaYB4/w3UNUPVbXA9/QzoHMI4zHG7K3LLoM5c9z1vtdeqxpzZEwECWXi6wT84vd8vW+tLpcCs/2eJ4nIYhH5TEROqWsnEZno227xNiu5Nib0Ro92s/26d4cvv4SBA92AW2MiRCgTX20XfWutpBGR84ABwBS/5a6qOgCYAPxHRHrUtq+qPqaqA1R1QLt27fY2ZmNMMA46yPX4HDYMNmyAo45ynV+MiQChTHzrgS5+zzsDNc6JiMixwF+Bk1W1uHJdVTf6vq4G5gH9QhirMWZPZWW5qQ7nnuuu9Y0b57q+tKBKcdMyhTLxLQIOEJH9RSQBOBsIqM4UkX7Ao7ikt9VvPVNEEn3fZwHDgO9CGKsxpjESE+HZZ11fz4oKuO46mDTJVX8aE6ZClvhUtQyYBMwFVgAvqepyEblNRCqrNKcAacDL1W5bOBhYLCJLgQ+Bu1XVEl89ysqDKy0vKi2vsVZYEtxaUWkZ5eXlNdbKqo0pKq+ooKKi4d/6KyqUciuJj3wicPPN8MILLhE+/LDr/LJ7t9eRGVMru4E9wm3PLebNbzbx5bpdnNq/E307Z5CZWnMM0NacIhb8uI0FK7czpEdbjjmoParKkp+zmbt8M4d3ac2JfToCyvebc5n51QYO3Ced0/p3JjYGVm/LZ9qiX+jSJoWzj+xCYpywflcRz332M21SEzh/yH5kpSawPruQJz9ZQ2pCHBcO7UaHVknExQb+flVeXsHmnGKeXriWvOJSLh62P50yk0lJsC4vEW/hQjjlFNfbs3dvd92vWzevowolu4E9Alnii2A78ou5/JnFfPlz9q9rfzn+IC4etj8JcVXJZld+MXfN/oGXFlcV2c6+5ije/GYTD8376de1/4zvS3ZhKbe8vvzXtYP3TeeBc/px7H0Lfl3rnJnMC5cNZsS/Pvz1ck5WWgKvTzqKUfd9SHGJW0xLjOPd649m39bJAXFv3l3Eb+6bT65vyG1sjPD2NUdxYIdWTfCnYjy3Zg2ceCKsWOGmu8+a5Xp9tkyW+CKQNamOYPlF5QFJD+DRBavJLgiclVdQUsGrX64PWGuVnMD/Fq4NWOu1TxpPfLQ6YG3FplyyC0pJ9Euk63cV8uPWXLLSEn9d255XwuK1O7lh9IG/ruUVl/Hu8i014n5vxZZfkx5AeYXy6ILVlJTVPL1qItD++7sjv2OPha1bYeRImD7d66iM+ZUlvggWU8t/vbgYQar9Dirijqr8qWqNU5AVCvGxNd80JkaoqHZmID42psZ1vIS4GPKrXRv0P/KslFjLWkJsDFI9cBO5MjLg7bdde7PiYjj7bLjjDqv4NGHBEl8ES02I45gDA+9dvGFMLzJTEqptF8uFQ7sFrOUUlTJpVM+AtVVbcrl29AEBa4P2b0NGcjxlfknuoA7p7J+Vyu7C0l/XOmcmc3iXDB5dUHXqNCstgVEHta8R94he7WiXXnW0mBgXw5UjetSadE0Ei4+HqVPh3nvdb19/+xtcdJFLhMZ4yK7xRbgdecUs+XkXS9dnc1zvfenaJpnWKTWLW7bkFPHjltxfi1sO6diKClXW7yrk3eVb6LdfBn27ZCAI23KLeGvZJg7cJ50hPdoSFyNszy/h9a830qVNCqMObE9SnLCzoJSZX20gKy2RMYfsQ9u0BLbmFvP61xtJSYjj+D4daJ+eWONITlXZmlvM3G83k1dcxsl9O9I+PZGEuNpHHZkWYNYsmDABCgpg+HCYORPatoghxXaaIgJZ4jPGNI8vv3Sz/TZuhJ49XcXngQc2vF94s8QXgezckjGmefTvD198Af36wapVMGQIzJvndVQmClniM8Y0n06dYMECOPlk2LULxoyBp57yOioTZSzxGWOaV1oavPoqXH89lJbCJZfA5Mk22NY0G0t8xpjmFxvrqj2nTnXf33UXjB/vil+MCTFLfMYY71xxBcyeDa1awSuvuJvdN2/2OirTwlniM8Z46ze/gU8/dT09Fy2CQYNg2TKvozItmCU+Y4z3DjnEDbYdMgTWrXMDbmfP9joq00JZ4jPGhIf27eGDD1x7s9xcN9rowQe9jsq0QDYHpgUqKi1nd2EpBSVlpCbEkZmaEHQ7MFVlW14xBcXlJMXHkp4UR2piHBuzCykqLSchLoaE2Bjat0oK8U9holJSEjz/PBxwANx+O1x9NaxcCffdB3H2z5VpGvY3qYUpLi1n4U87+P3zX1JYWk7r5HieuWQgh3VuHVQT6LU7Cpjw+Gds2l1EfKzwwDn9OKRjay74f5+zdkcBsTHCVSN7cN6grnSoNm7ImCYREwO33eaS32WXwQMPwE8/wbRpkJ7udXSmBbBTnS1MdmEp1774FYW+Seu7C0u5+sWv2J7XcGPg7IISJs9cxqbdRQCUliudMpL555zvWbvDlZmXVygPfrCKvGIbIWRC7Pzz4b33XE/Pt9921/3WrfM6KtMCWOJrYYpLywNm3QGs21lAeRD3BpeUVfD9ppyAtZSEWFZsyq2x7drt+XsVpzFBGT4cPvsMevVylZ4DB7rKT2P2giW+FiY5IZZOGYGnIPt3zah1Ll51KYlxjDww2f3+HAAACkdJREFUcIzQ9rwSRlYbfRQXIxzYwU45mWbSs6dLfqNGwZYtMGIEzJjhdVQmglnia2HapibyzCUD6dOpNSIwuHsbHpjQnzapNUcVVZeWGMdfTjiI43p3IDZG6NEuleTEGCYe3Z1T+3UkPlbonJnMExcOIDne/uqYZpSZCXPmwKWXQmEhnHEG3HOPDbY1jWJjiVqonXkllGkFCbExZNQyn68+OUWlFJWWE4OQ5RsYuz23iJJyRVVpn5ZIfLzNzjMeUIUpU+DPf3bPL7kEHnkEEvbs73gTsrFEEcgSnzEm8rz6Kpx3njv6GzXKnfrMzPQiEkt8EcjOVxljIs9pp7nxRh06wIcfwuDBbsafMUGwxGeMiUwDBrjBtocd5m5yHzQIPvrI66hMBLDEZ4yJXF26wMcfw4knws6dMHo0PPus11GZMGeJzxgT2dLTYdYsuOYaN9j2ggvg73+3wbamTpb4jDGRLzYW7r/ftTeLiYH/+z+YMMEVvxhTjSU+Y0zLMWkSvPmmOwqcPh2OOcbd9G6MH0t8xpiW5fjj4ZNPoGtX1/Fl0CBYvtzrqEwYscRnjGl5+vRxg20HDoSff4ahQ2HuXK+jMmEipIlPRI4TkR9EZJWI3FTL64kiMt33+uci0s3vtb/41n8Qkd+GMk5jTAvUoQPMmwdnngk5Oa7yc+pUr6MyYSBkiU9EYoGHgOOBQ4BzROSQaptdCuxS1Z7Av4F7fPseApwN9AaOAx72vZ8xxgQvOdnN8Zs8GcrL4aqr4Prr3fcmaoXyiG8gsEpVV6tqCTANGFdtm3HA077vXwFGi5uWOg6YpqrFqroGWOV7P2OM2TMxMXDHHfDUUxAfD//+N5x7rtdRGQ+FMvF1gv/f3t3H2FHVYRz/PtkW2yJiEKK1rWkDBSMoILWJVkSIGlSiJDamBCtVExIBgfAHVhLfYwL/GI2KpkCDShUbEFMx0hYBTX2hb5ZCKSBoTTclaRvkLQK19PGPOZteym579967TPfO80k2OzM7c+aZNNvfzsy557C9ZX2wbBt2H9t7gWeAN7V5LACSLpa0XtL6Xbt29Sh6RPSdRYtg9epqYtv58+tOEzWaMIZtDzd464EjYo+0TzvHVhvtJcASqAapHk3AiGiYs86CJ56Ao4+uO0nUaCzv+AaBGS3r04EdI+0jaQJwNPBUm8dGRIxeil7jjWXhWwfMljRL0hFUnVVWHLDPCuCisjwfuMfVPEkrgAWl1+csYDawdgyzRkREQ4zZo07beyVdBqwEBoCltrdI+haw3vYK4Cbg55Iep7rTW1CO3SJpOfAwsBe41Ha6YUVERNcyEW1EROcyEe04lJFbIiKiUVL4IiKiUVL4IiKiUVL4IiKiUVL4IiKiUVL4IiKiUVL4IiKiUVL4IiKiUfrqA+ySdgH/HuVhxwK7xyDOa61frgNyLYerXMur7bZ9bg/aiddQXxW+Tkhab3tO3Tm61S/XAbmWw1WuJfpFHnVGRESjpPBFRESjpPCVSWz7QL9cB+RaDle5lugLjX/HFxERzZI7voiIaJQUvoiIaJTGFj5J50p6VNLjkhbXnadTkpZK2inpobqzdEvSDEn3StoqaYukK+rO1ClJkyStlfRAuZZv1p2pG5IGJP1d0p11Z+mGpG2SHpS0SVJmrW6oRr7jkzQAPAZ8GBgE1gEX2H641mAdkPQB4HngZ7ZPqTtPNyRNBaba3ijpKGADcP44/XcRcKTt5yVNBNYAV9j+W83ROiLpKmAO8Abb59Wdp1OStgFzbPfLB/GjA02945sLPG77n7b3ALcCn6w5U0ds/wl4qu4cvWD7Sdsby/JzwFZgWr2pOuPK82V1Yvkal39lSpoOfBy4se4sEb3Q1MI3Ddjesj7IOP0Ptl9JmgmcDtxfb5LOlceDm4CdwGrb4/VavgdcDeyrO0gPGFglaYOki+sOE/VoauHTMNvG5V/j/UjS64HbgSttP1t3nk7Zftn2acB0YK6kcfcoWtJ5wE7bG+rO0iPzbL8b+ChwaXlVEA3T1MI3CMxoWZ8O7KgpS7Qo78NuB5bZ/nXdeXrB9tPAfcB4HMx4HvCJ8m7sVuAcSbfUG6lztneU7zuBO6hee0TDNLXwrQNmS5ol6QhgAbCi5kyNVzqE3ARstf3duvN0Q9Jxkt5YlicDHwIeqTfV6Nn+iu3ptmdS/Z7cY/szNcfqiKQjS6cpJB0JfAQY972hY/QaWfhs7wUuA1ZSdaBYbntLvak6I+mXwF+BkyQNSvpC3Zm6MA9YSHVXsal8fazuUB2aCtwraTPVH1qrbY/rjwL0gTcDayQ9AKwFfmf7rpozRQ0a+XGGiIhorkbe8UVERHOl8EVERKOk8EVERKOk8EVERKOk8EVERKOk8MVhQ9IiSW9tY7+bJc1vd3sPcl3Tsjyz3ZkwJF0p6bM9OP9lkj7XbTsRUUnhi8PJIuCQha8G1xx6l1eSNAH4PPCLHpx/KXB5D9qJCFL4YoyUO6NHJP1U0mZJt0maUn52hqQ/loGCV0qaWu7U5gDLygfXJ0v6mqR1kh6StKSM7NLu+V91jrL9PknXlbnyHpN0Ztk+RdLykvVXku6XNEfStcDkkmlZaX5A0g1lnr1VZWSWA50DbCyDJSDpBEl3l/n5Nko6XtIHS8blJcu1ki4s2R6UdDyA7f8C2yRleK2IHkjhi7F0ErDE9ruAZ4FLylicPwDm2z6D6m7mO7ZvA9YDF9o+zfYLwA9tv6fMMzgZaGseuJHO0bLLBNtzgSuBr5dtlwD/KVm/DZwBYHsx8ELJdGHZdzbwI9snA08Dnxomxjyq+QSHLCvHnAq8D3iybD8VuAJ4J9WoNSeWbDcCX2o5fj1wZjvXHxEHN6HuANHXttv+c1m+hepx3V3AKcDqcgM3wP4icKCzJV0NTAGOAbYAv23jvCcd4hxDg19vAGaW5fcD3wew/VAZamwk/7K9aZg2Wk2lGg6PMj7kNNt3lPZfLNsB1tl+sqw/Aawqxz8InN3S3k7g7QfJFBFtSuGLsXTgeHimmhJqi+33HuxASZOA66lmy94u6RvApDbPe6hzvFS+v8z+34G2H6O2HD/UxnCPOl9gf96Dtd3a1r6W9X288vdzUmkzIrqUR50xlt4maaj4XACsAR4FjhvaLmmipJPLPs8BR5XloaKxu8zPN5remgc7x0jWAJ8u+7+D6tHjkP+Vx6ejsRU4AaDMKTgo6fzS/uuG3neOwolkJoGInkjhi7G0FbioPDY8Bvix7T1URey6Mkr+Jqp3XgA3Az9RNWv5S8ANVI/8fkM1w0FbDnGOkVxPVSw3A18GNgPPlJ8tATa3dG5px++B1klOFwKXl/b/ArxlFG1B9c7w7lEeExHDyOwMMSYkzQTuLB1TDnuSBoCJtl8svSn/QNXRZE8Xbd4BXG37H11mOx24yvbCbtqJiEre8UVUplDNnzeR6p3cF7spesViqk4uXRU+4Fjgq122ERFF7vgiIqJR8o4vIiIaJYUvIiIaJYUvIiIaJYUvIiIaJYUvIiIa5f8gA+Z/MPnBhQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 463.5x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.relplot(x='petal length (cm)', y='petal width (cm)', hue='class', data=dataset.loc[:99])\n",
    "plt.plot([5, 0], [0, 1.4], 'r', linewidth=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is best to shuffle the data before training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.utils import shuffle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "training['input'], training['target'] = shuffle(training['input'], training['target'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing the network\n",
    "\n",
    "We will be implementing a fully connected feed forward network.\n",
    "\n",
    "Unlike the formula above, the implementation will be processing a layer at once.\n",
    "\n",
    "E.g,. calculating neurons output of a layer at once in matrix dot product, calculating loss and partial derivatives of loss at once using vectors operations, etc.\n",
    "\n",
    "In other words, to increase performance the implementation will take a vector of inputs, process them in vector form and outputs vectors of results.\n",
    "\n",
    "This is done with numpy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us start with implementing an activation function and loss function:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implementation of sigmoid activation function, and its derivative.\n",
    "\n",
    "$$\\sigma(x) = \\frac{1}{1+e^{\\text{-}x}}$$\n",
    "\n",
    "$$\\sigma^\\prime(x) = \\frac{e^{\\text{-}x}}{(1+e^{\\text{-}x})^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Sigmoid:\n",
    "    \n",
    "    @staticmethod\n",
    "    def forward(x):\n",
    "        return 1 / (1 + np.exp(-x))\n",
    "    \n",
    "    @staticmethod\n",
    "    def backward(x):\n",
    "        return np.exp(-x) / (1 + np.exp(-x)) ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `forward` function is the vector version of sigmoid activation function.\n",
    "\n",
    "The `backward` function is the vector version of sigmoid derivative function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.21416502, 0.42555748, 0.95689275, 0.9168273 ])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Sigmoid.forward(np.array([-1.3, -0.3, 3.1, 2.4]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.16829836, 0.24445831, 0.04124902, 0.076255  ])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Sigmoid.backward(np.array([-1.3, -0.3, 3.1, 2.4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normally, sigmoid takes in one input scalar and output one scalar. Here the implementation works on vectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loss function for classification\n",
    "\n",
    "Cross entropy loss is the most popular classification loss.\n",
    "\n",
    "$$E_{CEL}(\\vec x, \\vec t) = - \\sum_i t_i \\cdot \\text{log}(x_i)$$\n",
    "\n",
    "Where $\\vec x$ is a vector of output probability and $\\vec t$ is the vector of target probability.\n",
    "\n",
    "For performance, many implementations require the target to be in one-hot vectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the output need to be in probability form, we can use the softmax function to convert network output to probability.\n",
    "\n",
    "$$\\text{Softmax}(\\vec x) = \\frac{e^{\\vec x}}{\\sum_{x_i \\in \\vec x}e^{x_i}}$$\n",
    "\n",
    "There the vector $e^{\\vec x}$ is being element-wise divided by the scalar $\\sum_{x_i \\in \\vec x}e^{x_i}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def softmax(x):\n",
    "    return np.exp(x) / np.sum(np.exp(x))\n",
    "\n",
    "\n",
    "class CrossEntropyLoss:\n",
    "    \n",
    "    @staticmethod\n",
    "    def forward(x, t):\n",
    "        x = softmax(x)\n",
    "        return -np.log(x) @ t\n",
    "        \n",
    "    @staticmethod\n",
    "    def backward(x, t):\n",
    "        x = softmax(x)\n",
    "        return x - t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `forward` function is the cross entropy loss with softmax assuming one-hot targets.\n",
    "\n",
    "The `backward` function is the derivative of cross entropy loss with softmax assuming one-hot targets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple network implementation\n",
    "\n",
    "This is a simple fully connected feed forward neural network implementation.\n",
    "\n",
    "The same activation function is used in all the neurons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def network(inputs, weights, activation_fn):\n",
    "    x = inputs\n",
    "    for w in weights:\n",
    "        x = activation_fn.forward(x @ w)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the implementation of backpropagation for the network.\n",
    "\n",
    "Recall the formula:\n",
    "\n",
    "\\begin{align}\n",
    "& \\frac{\\partial E}{\\partial w_{ij}} = \\delta_j \\cdot o_i \\\\\n",
    "& \\delta_j = \\begin{cases} \\frac{\\partial E(t_j, o_j)}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{if $j$ is an output neuron} \\\\ \\bigl( \\sum_{l \\in L} \\delta_l \\cdot w_{jl} \\bigr) \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{otherwise}\\end{cases}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backpropagation(inputs, targets, weights, activation_fn, loss_fn, learning_rate):\n",
    "    \"\"\"\n",
    "    Backpropagation for a fully connected feed forward network.\n",
    "    \n",
    "    :param inputs: A vector of network input\n",
    "    :param targets: A vector of the target for network output\n",
    "    :param weights: A list of layer weight matrix\n",
    "    :param activation_fn: The activation function to use\n",
    "    :param loss_fn: The loss function to use\n",
    "    :param learning_rate: The learning rate for update rule\n",
    "    :returns: A list of updated layer weight matrix\n",
    "    \"\"\"\n",
    "    \n",
    "    # Rename variable to better match the formula\n",
    "    x = inputs\n",
    "    t = targets\n",
    "    h = activation_fn\n",
    "    loss = loss_fn\n",
    "    \n",
    "    n = len(weights)  # Number of hidden and ouput layer(s)\n",
    "    params = {'w': weights}\n",
    "    grads = {'w': [None] * n}\n",
    "    cache = {'input': [None] * n, 'delta': [None] * n}\n",
    "    \n",
    "    # Forward pass\n",
    "    for i in range(n):\n",
    "        w = params['w'][i]\n",
    "        cache['input'][i] = x  # Save layer inputs\n",
    "        s = x @ w\n",
    "        x = h.forward(s)\n",
    "\n",
    "    output = x  # Network output\n",
    "    \n",
    "    # Backward pass\n",
    "    for i in range(n):\n",
    "        # Start from last layer\n",
    "        layer = (n - 1) - i  # Index of current layer\n",
    "        o_i = cache['input'][layer]  # Input of current layer / Output of previous layer\n",
    "        w_j = params['w'][layer]  # Weight of current layer\n",
    "        s_j_der = h.backward(o_i @ w_j)\n",
    "        \n",
    "        if layer == (n - 1):  # If last layer\n",
    "            delta_j = loss.backward(output, t) * s_j_der\n",
    "        else:\n",
    "            delta_l = cache['delta'][layer + 1]  # Delta of next layer\n",
    "            w_l = params['w'][layer + 1]  # Weight of next layer\n",
    "            delta_j = w_l @ delta_l * s_j_der\n",
    "            \n",
    "        grads['w'][layer] = np.outer(o_i, delta_j)\n",
    "        cache['delta'][layer] = delta_j\n",
    "    \n",
    "    # Update parameters\n",
    "    for i in range(n):\n",
    "        params['w'][i] -= learning_rate * grads['w'][i]\n",
    "    \n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the above implementation, we will create a network with: input layer of 2 units, a single hidden layer of 4 units, output layer of 2 units.\n",
    "\n",
    "The activation function for all the neurons is sigmoid and the loss function is cross entropy loss."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let see how well a network will random weight can predict the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Network accuracy: 50.0%\n"
     ]
    }
   ],
   "source": [
    "# Weights for the input layer, hidden layers, and output layers.\n",
    "weights = [np.random.randn(2, 4), np.random.randn(4, 2)]\n",
    "# Test the network:\n",
    "correct_count = 0\n",
    "total_count = 0\n",
    "for inputs, targets in zip(training['input'], training['target']):\n",
    "    outputs = network(inputs, weights, Sigmoid)\n",
    "    \n",
    "    total_count += 1\n",
    "    if np.argmax(outputs) == np.argmax(targets):\n",
    "        correct_count += 1\n",
    "\n",
    "accu = correct_count / total_count\n",
    "\n",
    "print(f\"Network accuracy: {accu * 100}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected the accuracy is not very high. At around 50%, the network performs as well as random guesses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will train the network with the training data.\n",
    "\n",
    "We did not learn about when to stop the training yet. We will stop the training after a fix number of iterations for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At    10 epoch, mean loss of 10 epoch: 69.02026085753617\n",
      "At    20 epoch, mean loss of 10 epoch: 66.50513172885509\n",
      "At    30 epoch, mean loss of 10 epoch: 61.922925339851645\n",
      "At    40 epoch, mean loss of 10 epoch: 56.445128531079305\n",
      "At    50 epoch, mean loss of 10 epoch: 50.87498130849378\n",
      "At    60 epoch, mean loss of 10 epoch: 46.382373547445304\n",
      "At    70 epoch, mean loss of 10 epoch: 43.33148947752215\n",
      "At    80 epoch, mean loss of 10 epoch: 41.36343649376842\n",
      "At    90 epoch, mean loss of 10 epoch: 40.06364224312124\n",
      "At   100 epoch, mean loss of 10 epoch: 39.163992864283706\n"
     ]
    }
   ],
   "source": [
    "# Weights for the input layer, hidden layers, and output layers.\n",
    "weights = [np.random.rand(2, 4), np.random.rand(4, 2)]\n",
    "\n",
    "learning_rate = 0.1\n",
    "\n",
    "epoch = 100  # How many time to train the network with the whole training data.\n",
    "num_e = 10  # Print info every num_e epoch\n",
    "training_loss = 0\n",
    "\n",
    "for e in range(epoch):\n",
    "    for inputs, targets in zip(training['input'], training['target']):        \n",
    "        outputs = backpropagation(inputs, targets, weights, Sigmoid, CrossEntropyLoss, learning_rate)\n",
    "\n",
    "        training_loss += CrossEntropyLoss.forward(outputs, targets) / num_e\n",
    "\n",
    "    # Print training info every set amount of epoch\n",
    "    if not (e + 1) % num_e:\n",
    "        print(f\"At {e + 1: 5} epoch, mean loss of {num_e} epoch: {training_loss}\")\n",
    "        training_loss = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We often print out the loss during training to see if the training is going well.\n",
    "The loss should be decreasing during training. Once the loss stops decreasing, the training is complete."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us test the trained network on the training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Network accuracy: 98.0%\n"
     ]
    }
   ],
   "source": [
    "# Test the network:\n",
    "correct_count = 0\n",
    "total_count = 0\n",
    "for inputs, targets in zip(training['input'], training['target']):\n",
    "    outputs = network(inputs, weights, Sigmoid)\n",
    "\n",
    "    total_count += 1\n",
    "    if np.argmax(outputs) == np.argmax(targets):\n",
    "        correct_count += 1\n",
    "\n",
    "accu = correct_count / total_count\n",
    "\n",
    "print(f\"Network accuracy: {accu * 100}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we are not performing model selection, there are some chance that the accuracy is low. If encounter retrain the model.\n",
    "\n",
    "Else, the accuracy should be very high. In fact, it is probably overfitting, since we are testing on the training data.\n",
    "\n",
    "Regardless, this shows the network is learning from the training data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To examine the test of a single case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The network prediction:\tIris-versicolor\n",
      "The correct answer:\tIris-versicolor\n"
     ]
    }
   ],
   "source": [
    "# Test a single case (0 to 99)\n",
    "i = 99\n",
    "\n",
    "inputs = training['input'][i]\n",
    "outputs = network(inputs, weights, Sigmoid)\n",
    "output_label = to_class[np.argmax(outputs)]\n",
    "\n",
    "targets = training['target'][i]\n",
    "target_label = to_class[np.argmax(targets)]\n",
    "\n",
    "print(f\"The network prediction:\\t{output_label}\\nThe correct answer:\\t{target_label}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding  bias terms\n",
    "\n",
    "While in the deriving the backpropagation formula, we focus on the tuning the weights. However, backpropagation can used to tune any parameter in the network.\n",
    "\n",
    "A common tunable parameter is the bias.\n",
    "\n",
    "Recall the bias term is added to the weighted sum of a neuron.\n",
    "\n",
    "$s_j = \\bigl( \\sum_{i \\in I} w_{ij} \\cdot x_i \\bigr) + b_j$\n",
    "\n",
    "where the set $I$ is all of the inputs for neuron $j$, and $w_{ij}$ is the weight associated with the connection from input $i$ to neuron $j$. \n",
    "\n",
    "The full neuron formula becomes:\n",
    "\n",
    "$o_j = h_j\\Bigl(\\bigl( \\sum_{i \\in I} w_{ij} \\cdot x_i \\bigr) + b_j \\Bigr)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like tuning weights, the update rule for bias is the same:\n",
    "\n",
    "$$\\Delta \\vec b = - \\eta \\nabla E$$\n",
    "\n",
    "Where $\\vec b$ is the vector of all the bias in the network, $\\eta$ is a learning rate, $E$ is the loss.\n",
    "\n",
    "or for a single bias:\n",
    "\n",
    "$$\\Delta b_{i} = - \\eta \\frac{\\partial E}{\\partial b_{i}}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like for partial derivatives for weights, the chain rule is used.\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial E}{\\partial b_j} &= \\frac{\\partial E}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} \\cdot \\frac{\\partial s_j}{\\partial b_j} \\\\\n",
    "&=\\delta_j \\cdot \\frac{\\partial s_j}{\\partial b_j}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\frac{\\partial s_j}{\\partial b_j}$ term is:\n",
    "$$ \\frac{\\partial s_j}{\\partial b_j} = \\frac{\\partial}{\\partial b_j} \\bigl( \\sum_{i \\in I} w_{ij} \\cdot x_i \\bigr) + b_j = \\frac{\\partial}{\\partial b_j} b_j = 1$$\n",
    "\n",
    "Therefore:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial b_j} = \\delta_j \\cdot 1 = \\delta_j$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The whole formula becomes:\n",
    "\n",
    "$$\\frac{\\partial E}{\\partial b_j} = \\delta_j = \\begin{cases} \\frac{\\partial E(t_j, o_j)}{\\partial o_j} \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{if $j$ is an output neuron} \\\\ \\bigl( \\sum_{l \\in L} \\delta_l \\cdot w_{jl} \\bigr) \\cdot \\frac{\\partial h_j(s_j)}{\\partial s_j} & \\text{otherwise}\\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practice "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Modify the network function to include bias term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "Modify the backpropagation function to include bias term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the same network with bias term with Sigmoid activation function may be more unstable and requires good learning rate adjustment for each parameter to produce consistent good result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3\n",
    "\n",
    "Train the network with bias term.\n",
    "\n",
    "Same network structure and training data.\n",
    "\n",
    "Network: input layer 2 units, one hidden layer 10 units, output layer 2 units.\n",
    "\n",
    "Initialize biases to 0.\n",
    "\n",
    "Learning rate: 0.1\n",
    "\n",
    "Epoch: 100\n",
    "\n",
    "Use the ReLU activation function as implemented below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ReLU:\n",
    "    \n",
    "    @staticmethod\n",
    "    def forward(x):\n",
    "        return np.maximum(0, x)\n",
    "    \n",
    "    @staticmethod\n",
    "    def backward(x):\n",
    "        x[x <= 0] = 0\n",
    "        x[x > 0] = 1\n",
    "        return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4\n",
    "\n",
    "Test the network. The accuracy may be lower than without bias. That is fine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The functions used in the network must be differentiable\n",
    "\n",
    "Note that the partial derivative formula for backpropagation implies that both the loss function $E(\\cdot)$ and all the activation functions $h_j(\\cdot)$ need to be differentiable. Meaning they must have a derivative.\n",
    "\n",
    "Technically, it is possible to just have a defined derivative for each value in input space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, a version of ReLU (Rectified Linear Unit):\n",
    "\n",
    "$$\\text{ReLU}(x) = \\text{max}(0, x)$$\n",
    "\n",
    "Which looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def relu(x):\n",
    "    return np.maximum(0, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAYwElEQVR4nO3dd3hUdfoF8PMSCB0CJHRC6EUhhUi1ggVB0NUVEeKuFQ2gYF2UVVd31bUsYgGUdfnpbkJHBBVRLIhYgHRIQugQCJAJkEp63t8fCS6LgdyEuXPvzJzP8/gIJmZOHuDk8s2dM6KqICIi+6pndQAiIrowFjURkc2xqImIbI5FTURkcyxqIiKbq2/GB/X399egoCAzPjQRkUeKjY3NUtWA6t5mSlEHBQUhJibGjA9NROSRROTg+d7Gow8iIptjURMR2RyLmojI5ljUREQ2x6ImIrI5Q3d9iMgBAHkAygGUqWq4maGIiOi/anN73jWqmmVaEiIiqhaPPoiInGDr/pP44Id9MGM62mhRK4CvRCRWRKZU9w4iMkVEYkQkxuFwOC8hEZHNZeYVYdriOERvOYTC0nKnf3yjRT1CVcMA3Ahgmohcee47qOpCVQ1X1fCAgGqfBUlE5HHKyivw8OJ45BWVYkFEGJr4Ov8J34aKWlUzqv6dCWA1gMFOT0JE5Ibe+GoXtuw/iZd/NwB927cw5TFqLGoRaSoizc/8GMD1AHaYkoaIyI1sSDmO977fi0lDAnFrWGfTHsfINXo7AKtF5Mz7L1bV9aYlIiJyAwdPFOCx5QkY0Kklnrupv6mPVWNRq+o+AMGmpiAiciNFpeWIjIpDPRHMnxyGRg18TH08U2ZOiYg82fNrkpFyNBeL7g5Hl9ZNTH883kdNRFQLy7elY1lMOqZf0xMj+7ZzyWOyqImIDErOyMGza3ZgRM82ePS63i57XBY1EZEBOYWliIyKQ6smvnhrYih86onLHptn1ERENVBVPLEiERnZhVj24DD4N2vo0sfnFTURUQ3e37QPG1KO45kx/TCoayuXPz6LmojoAn7ZdwKvrd+JsQM74J4RQZZkYFETEZ1HZm4Rpi+OR5B/U7x620BUPfHP5XhGTURUjbLyCkxfEo+C4jIsfmAImjW0ri5Z1ERE1Xj9yzRs3X8Sc+8IQe92zS3NwqMPIqJzfJl8DO9v2oeIoYG4JbST1XFY1EREZzuQVYAnliciuHNLPGvy2JJRLGoioipFpeWIjI6Dj49g3uQwNKxv7tiSUTyjJiJC5ZNa/vzJDuw8lotFd1+Gzq3MH1syilfUREQAlm1Lx8rYw3j4mp64pk9bq+P8DxY1EXm9HUdy8NzaZFzRyx8zrnXd2JJRLGoi8mo5p0sRGR2LNk19MfeOEJeOLRnFM2oi8loVFYrHVyTgWE4Rlj04DG1cPLZkFK+oichrvbdpL75OzcTsMf0QFuj6sSWjWNRE5JV+2puFN75Mw7jgjvjj8CCr41wQi5qIvM7x3CI8siQe3fyb4u+3DrBsbMkonlETkVcpLa/A9MVxOF1SjiUPDEVTC8eWjLJ/QiIiJ3pt/U5sO3AKb00MQS+Lx5aM4tEHEXmN9TuO4p8/7McfhnXFzSHWjy0ZxaImIq+wP6sAT65IQnAXP8we28/qOLXCoiYij1dYUo7IqFjU9xHMt9HYklE8oyYij6aqmP3JdqQdz8OH9wxGJ7/GVkeqNV5RE5FHW7I1HR/HHcEjI3vhqt4BVsepExY1EXms7Ydz8Je1ybiydwAeGdXL6jh1xqImIo+UfboEkdGx8G9m37Elo3hGTUQep6JC8djyRBzPLcKKh4ajdVNfqyNdFMNX1CLiIyLxIvKZmYGIiC7Wgu/34tudmXj2pv4I6eJndZyLVpujjxkAUs0KQkTkDD/uycI/vkrD+OCOuGtoV6vjOIWhohaRzgDGAvjA3DhERHV3LKdybKl7QDO84gZjS0YZvaKeC+ApABXnewcRmSIiMSIS43A4nBKOiMioM2NLhaXleC8izC3GloyqsahF5CYAmaoae6H3U9WFqhququEBAe55ryIRua+/f7ETMQdP4dXbBqJnW/cYWzLKyBX1CADjReQAgKUARopIlKmpiIhqYd32o/jX5v24e3gQxgV3tDqO09VY1Kr6tKp2VtUgABMBfKuqEaYnIyIyYK8jH0+tTEJooB+eGeNeY0tG8QkvROS2TpeUITIqFr7162HepDD41vfMSqvVabuqbgSw0ZQkRES1oKqYvXoHdmfm49/3DkZHNxxbMsozv/wQkceL3nIIq+OPYOao3riil2ffwMCiJiK3k3Q4Gy9+moKr+wTg4ZE9rY5jOhY1EbmVUwUliIyKQ0DzhnhzQgjqufHYklGec0c4EXm8igrFo8sT4MgrxoqHhqGVm48tGcUraiJyG/O+24ONaQ48O64/gj1gbMkoFjURuYXNu7Mw5+tduCWkIyKGBFodx6VY1ERke0dzCvHI0nj0atsML3vQ2JJRLGoisrWSsgpMi45DcWk5FkQMQhNf7/vWmvd9xkTkVl75IhVxh7Ixb1IYegQ0szqOJXhFTUS29VlSBv7vxwO4Z0QQxg7sYHUcy7CoiciW9mTm408rkxAW6Ienb/TMsSWjWNREZDsFxZVjSw0b+GDeZM8dWzKKZ9REZCuqimdWb8ceRz7+c+8QdGjpuWNLRnn3lykisp2oXw5iTUIGHr+uNy7v5W91HFtgURORbSSkZ+PFz1Iwsm9bTL3a88eWjGJRE5EtnCoowbToOLRr0QhzJgR7xdiSUTyjJiLLVVQoZi6rHFtaGTkMfk28Y2zJKF5RE5Hl3vl2D77f5cDz4/tjYGfvGVsyikVNRJbatMuBud/swq2hnTBpsHeNLRnFoiYiy2RkF2LG0nj0btscL/3O+8aWjGJRE5ElSsoqMDU6DqXligURYWjs62N1JNviNxOJyBIvr0tFQno25k8OQ3cvHVsyilfURORyaxMz8OFPB3Df5d0wZoD3ji0ZxaImIpfafTwPs1YlIbxrK8y6sa/VcdwCi5qIXKaguAyR0XFo4uuDdyeFoYEPK8gInlETkUuoKmZ9vB37HPmIum8I2rdsZHUkt8EvZ0TkEv/++SA+TczA49f3wfCeHFuqDRY1EZku7tAp/O3zFIzq2xaRV/WwOo7bYVETkalOFpRgenQc2rdshDkTQji2VAc8oyYi05RXKGYsjUdWQQk+jhyOlk0aWB3JLdV4RS0ijURkq4gkikiyiLzgimBE5P7e/mY3ftidhRfGX4JLO7W0Oo7bMnJFXQxgpKrmi0gDAJtF5AtV/cXkbETkxjamZeLtb3fjtrDOmHhZF6vjuLUai1pVFUB+1U8bVP2jZoYiIvd2JLsQM5cloE+75vjbLZdybOkiGfpmooj4iEgCgEwAG1R1SzXvM0VEYkQkxuFwODsnEbmJ4rJyTI2OQ3m5YkHEII4tOYGholbVclUNAdAZwGARubSa91moquGqGh4QEODsnETkJv72WSoS07Px+u0D0c2/qdVxPEKtbs9T1WwAGwGMNiUNEbm1NQlH8J9fDuKBK7ph9KUcW3IWI3d9BIiIX9WPGwO4FsBOs4MRkXvZdTwPs1Ztx2VBrfDUaI4tOZORuz46APhIRHxQWezLVfUzc2MRkTvJLy7DQ1GxaNqwPseWTGDkro8kAKEuyEJEbkhV8adVSTiQVYDo+4eiXQuOLTkbv+wR0UX58KcD+DzpKJ68oS+G9WhjdRyPxKImojqLPXgKL32eimv7tcNDV3W3Oo7HYlETUZ2cyC/G9MVx6OjXGP+YEMwntZiIo0xEVGuVY0sJOHFmbKkxx5bMxCtqIqq1t77ehc17svDXmzm25AosaiKqle/SMvH2t3tw+6DOuOOyQKvjeAUWNREZln7yNB5dloB+HVrgr7f8ZkmCTMKiJiJDisvKMW1x1djS5DA0asCxJVfhNxOJyJAXP01B0uEcvH/XIARxbMmleEVNRDVaHX8Y0VsO4cEru+OGS9pbHcfrsKiJ6ILSjuXh6Y+3Y3C31njyhj5Wx/FKLGoiOq+8olJERsWiWcMGePfOUNTn2JIleEZNRNU6M7Z08ORpLL5/CNpybMky/PJIRNVa9OMBrNt+DE/d0AdDunNsyUosaiL6jZgDJ/HKulRc378dplzJsSWrsaiJ6H9k5Rdj2uI4dGrVGK/fzrElO+AZNRH9qnJsKR7Zp0uxeupgji3ZBIuaiH715oZd+HHPCbz2+4Ho37GF1XGoCo8+iAgA8E3qcbz73R7cEd4FE8K7WB2HzsKiJqJfx5b6d2iBF26+xOo4dA4WNZGXKyotR2R0LBTAexGDOLZkQzyjJvJyL3yagh1HcvHPP4QjsE0Tq+NQNXhFTeTFVsUexpKth/DQVT1wXf92Vseh82BRE3mpncdyMfuT7RjavTWeuL631XHoAljURF4ot6gUkVFxaNGoAd7m2JLt8YyayMuoKp5akYRDJ09jyQND0bY5x5bsjl9GibzMvzbvx/rkY5g1ui8Gd2ttdRwygEVN5EW2HTiJV77YidGXtMf9V3SzOg4ZxKIm8hKOvGJMi45Dl1aN8drtAzm25EZ4Rk3kBcrKK/DIknjkFpXio3sHo0Ujji25ExY1kRf4x4Zd+HnfCbxxezD6deDYkrup8ehDRLqIyHcikioiySIywxXBiMg5NqQcx4KNe3Hn4C74/aDOVsehOjByRV0G4HFVjROR5gBiRWSDqqaYnI2ILtKhE6fx2PIEXNqpBZ4fx7Eld1XjFbWqHlXVuKof5wFIBdDJ7GBEdHHOjC0JgAWTObbkzmp114eIBAEIBbClmrdNEZEYEYlxOBzOSUdEdfaXtclIzsjFm3eEoEtrji25M8NFLSLNAKwCMFNVc899u6ouVNVwVQ0PCAhwZkYiqqUVMelYui0dU6/ugVH9OLbk7gwVtYg0QGVJR6vqx+ZGIqKLkZKRiz9/sgPDurfBY9dxbMkTGLnrQwD8C0Cqqs4xPxIR1VVuUSmmRsfCrwnHljyJkV/FEQDuAjBSRBKq/hljci4iqiVVxRPLE3H4VCHmTQpDQPOGVkciJ6nx9jxV3QyAzzUlsrl//rAPX6Ucx5/H9kN4EMeWPAn/XkTkAbbsO4FX16dhzID2uO9yji15GhY1kZvLzCvC9CXx6Nq6CV69jWNLnohbH0RurKy8AtMXxyOvqBT/uW8wmnNsySOxqInc2OtfpWHr/pOYMyEYfdtzbMlT8eiDyE19lXwM73+/D5OGBOLWMI4teTIWNZEbOniiAI+vSMSATi3x3E39rY5DJmNRE7mZotJyPBQVh3oimD85jGNLXoBn1ERu5rk1O5B6NBeL7g7n2JKX4BU1kRtZvi0dy2MOY/o1PTGyL8eWvAWLmshNJGfk4Nk1OzCiZxs8yrElr8KiJnIDOYWliIyKQ6smvnh7Yih86vFJLd6EZ9RENqeqeGJFIjKyC7HswWFo04xjS96GV9RENvf+pn3YkHIcz4zph0FdW1kdhyzAoiaysZ/3nsBr63di7MAOuGdEkNVxyCIsaiKbyswtwsNL4hHk35RjS16OZ9RENlRaNbZUUFyGxQ8MQbOG/KPqzfirT2RDr3+Zhq0HTmLuHSHo3a651XHIYjz6ILKZ9TuOYeGmfYgYGohbQjtZHYdsgEVNZCP7swrw5IpEBHduiWc5tkRVWNRENlFYUo7IqFj4+AjmTQ5Dw/ocW6JKPKMmsgFVxbNrdiDteB4W3X0ZOrfi2BL9F6+oiWxg2bZ0rIw9jIdH9sI1fdpaHYdshkVNZLEdR3Lw3NpkXNHLHzNG9bI6DtkQi5rIQjmnSxEZHYs2TX3xFseW6Dx4Rk1kkYoKxeMrEnAspwjLHhyG1k19rY5ENsUraiKLvLdpL75OzcTsMf0QFsixJTo/FjWRBX7am4U3vkzDuOCO+OPwIKvjkM2xqIlc7FhOER5ZEo9u/k3x91sHcGyJasQzaiIXqhxbisPpknIseWAomnJsiQzg7xIiF3r1i52IOXgKb00MQS+OLZFBNR59iMgiEckUkR2uCETkqb7YfhQfbN6PPwzriptDOLZExhk5o/4QwGiTcxB5tH2OfDy5MgnBXfwwe2w/q+OQm6mxqFV1E4CTLshC5JEKS8oxNToODXwE8zm2RHXgtLs+RGSKiMSISIzD4XDWhyVya6qK2Z9sR9rxPMydGIpOfo2tjkRuyGlFraoLVTVcVcMDAgKc9WGJ3NqSren4OO4IZozqhat6888F1Q3voyYyyfbDOfjL2mRc2TsAj4zk2BLVHYuayATZp0sQGR0L/2a+mHtHCOpxbIkugpHb85YA+BlAHxE5LCL3mR+LyH1VVCgeW56I47lFmB8xiGNLdNFqfMKLqt7piiBEnmL+xj34dmcmXrz5EoR08bM6DnkAHn0QOdGPe7IwZ8MujA/uiLuGdrU6DnkIFjWRk5wZW+oe0AyvcGyJnIhbH0ROUFpegWmL41BYWo5lEWEcWyKn4u8mIid4Zd1OxB48hXfuDEXPthxbIufi0QfRRfo86SgW/bgfdw8PwrjgjlbHIQ/Eoia6CHsd+XhqZSJCA/3wzBiOLZE5WNREdXS6pAyRUbFo2MAH8yaFwbc+/ziROXhGTVQHqorZq3dgd2Y+/n3vYHTk2BKZiJcARHUQveUQVscfwaPX9sYVvTi2ROZiURPVUtLhbLz4aQqu7hOA6df0tDoOeQEWNVEtnCooQWRUHAKaN8SbEzi2RK7BM2oigyoqFI8uT4AjrxgrHhqGVhxbIhfhFTWRQe9+twcb0xx4dlx/BHNsiVyIRU1kwA+7HXjz6124JaQjIoYEWh2HvAyLmqgGGdmFmLE0Ab3aNsPLHFsiC7CoiS6gpKxybKm4tBwLIgahiS+/rUOux991RBfw8rpUxB/KxrxJYegR0MzqOOSleEVNdB6fJmbgw58O4J4RQRg7sIPVcciLsaiJqrEnMx+zViUhLNAPT9/IsSWyFoua6BwFxWeNLU3m2BJZj78Dic6iqnhm9XbsdeTjnTtD0aElx5bIeixqorNE/XIQaxIy8Nh1vTGip7/VcYgAsKiJfpWQno0XP0vByL5tMfVqji2RfbCoiQCcLCjB1KhYtGvRCHMmBHNsiWyF91GT1yuvUMxcloCs/BKsjBwGvyYcWyJ74RU1eb13vt2NTbsceH58fwzszLElsh8WNXm173c58NY3u3FraCdMGsyxJbInFjV5rSPZhZi5NB692zbHS7/j2BLZF4uavFJJWQWmRcehtFyxICIMjX19rI5EdF78ZiJ5pZc+T0FCejbmTw5Dd44tkc3xipq8ztrEDHz080Hcd3k3jBnAsSWyP0NFLSKjRSRNRPaIyCyzQxGZZffxPMxalYTwrq0w68a+VschMqTGohYRHwDzANwIoD+AO0Wkv9nBiJytoLgMkdFxaOLrg3cnhaGBD/9CSe7ByBn1YAB7VHUfAIjIUgA3A0hxdphx72xGUWm5sz8sEQAgr6gMmXlFiLp/CNq3bGR1HCLDjBR1JwDpZ/38MIAh576TiEwBMAUAAgPrdj9qj4CmKCmvqNP/S2TEDZe0x/AeHFsi92KkqKu7uVR/8x9UFwJYCADh4eG/ebsRcyeG1uV/IyLyaEYO6Q4D6HLWzzsDyDAnDhERnctIUW8D0EtEuomIL4CJANaaG4uIiM6o8ehDVctEZDqALwH4AFikqsmmJyMiIgAGn5moqusArDM5CxERVYM3khIR2RyLmojI5ljUREQ2x6ImIrI5Ua3Tc1Mu/EFFHAAOOv0Dm8sfQJbVIVyMn7N34OfsHrqqakB1bzClqN2RiMSoarjVOVyJn7N34Ofs/nj0QURkcyxqIiKbY1H/10KrA1iAn7N34Ofs5nhGTURkc7yiJiKyORY1EZHNsairISJPiIiKiMe/FIiIvC4iO0UkSURWi4if1ZnM4G0v0CwiXUTkOxFJFZFkEZlhdSZXEREfEYkXkc+szuIsLOpziEgXANcBOGR1FhfZAOBSVR0IYBeApy3O43Re+gLNZQAeV9V+AIYCmOYFn/MZMwCkWh3CmVjUv/UmgKdQzcuNeSJV/UpVy6p++gsqX8HH0/z6As2qWgLgzAs0eyxVPaqqcVU/zkNlcXWyNpX5RKQzgLEAPrA6izOxqM8iIuMBHFHVRKuzWOReAF9YHcIE1b1As8eX1hkiEgQgFMAWa5O4xFxUXmh51KtkG3rhAE8iIl8DaF/Nm2YDeAbA9a5NZL4Lfc6quqbqfWaj8q/L0a7M5iKGXqDZE4lIMwCrAMxU1Vyr85hJRG4CkKmqsSJytdV5nMnrilpVr63uv4vIAADdACSKCFB5BBAnIoNV9ZgLIzrd+T7nM0TkjwBuAjBKPfPGeq98gWYRaYDKko5W1Y+tzuMCIwCMF5ExABoBaCEiUaoaYXGui8YnvJyHiBwAEK6q7rbAVSsiMhrAHABXqarD6jxmEJH6qPxG6SgAR1D5gs2TPPm1P6XyauMjACdVdabVeVyt6or6CVW9yeoszsAzanoXQHMAG0QkQUTeszqQs1V9s/TMCzSnAljuySVdZQSAuwCMrPp1Tai60iQ3xCtqIiKb4xU1EZHNsaiJiGyORU1EZHMsaiIim2NRExHZHIuaiMjmWNRERDb3/71ISaBdfUfsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-5, 5, 3)  # We only need 3 points to draw ReLU\n",
    "y = relu(x)\n",
    "\n",
    "plt.plot(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For ReLU function, there is no derivative at $x=0$ defined by calculus.\n",
    "\n",
    "\n",
    "However, we can just define a value for it. Usually the derivative of ReLU is defined as:\n",
    "\n",
    "$$\\text{ReLU}^ \\prime (x) = \\begin{cases}0 & x \\le 0 \\\\1 & x \\gt 0\\end{cases}$$\n",
    "\n",
    "So at $x=0$ the derivative is defined to be $0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training neural networks\n",
    "\n",
    "Training complex networks can be difficult. There are massive improvements to be made to our implementation.\n",
    "\n",
    "The update rule we implement is the original gradient descent, which can be massively improved. No learning rate schedule was used either.\n",
    "\n",
    "The implementation of our function does not consider numerical stability. For example, our softmax function is not  numerical stable.\n",
    "\n",
    "There can be better choice of activation functions. Since sigmoid activation have vanishing gradient problem in large networks.\n",
    "\n",
    "We did not perform model selection using validation data.\n",
    "\n",
    "In the next lesson, we will about about common problems and best practice for training neural networks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra:\n",
    "\n",
    "Here is a website that lets you play around with the training of a neural network. To see how different setting affects the training and result.\n",
    "\n",
    "http://playground.tensorflow.org"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution to exercise 1\n",
    "\n",
    "def network_bias(inputs, weights, biases, activation_fn):\n",
    "    x = inputs\n",
    "    for w, b in zip(weights, biases):\n",
    "        x = activation_fn.forward((x @ w) + b)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution to exercise 2\n",
    "\n",
    "def backpropagation_bias(inputs, targets, weights, biases, activation_fn, loss_fn, learning_rate):\n",
    "    \"\"\"\n",
    "    Backpropagation for a fully connected feed forward network with bias.\n",
    "    \n",
    "    :param inputs: A vector of network input\n",
    "    :param targets: A vector of the target for network output\n",
    "    :param weights: A list of layer weight matrix\n",
    "    :param biases: A list of layer bias vector\n",
    "    :param activation_fn: The activation function to use\n",
    "    :param loss_fn: The loss function to use\n",
    "    :param learning_rate: The learning rate for update rule\n",
    "    :returns: A list of updated layer weight matrix and bias vector\n",
    "    \"\"\"\n",
    "    \n",
    "    # Rename variable to better match the formula\n",
    "    x = inputs\n",
    "    t = targets\n",
    "    h = activation_fn\n",
    "    loss = loss_fn\n",
    "    \n",
    "    n = len(weights)  # Number of hidden and ouput layer(s)\n",
    "    params = {'w': weights, 'b': biases}\n",
    "    grads = {'w': [None] * n, 'b': biases}\n",
    "    cache = {'input': [None] * n, 'delta': [None] * n}\n",
    "    \n",
    "    # Forward pass\n",
    "    for i in range(n):\n",
    "        w = params['w'][i]\n",
    "        b = params['b'][i]\n",
    "        cache['input'][i] = x  # Save layer inputs\n",
    "        s = x @ w + b\n",
    "        x = h.forward(s)\n",
    "\n",
    "    output = x  # Network output\n",
    "    \n",
    "    # Backward pass\n",
    "    for i in range(n):\n",
    "        # Start from last layer\n",
    "        layer = (n - 1) - i  # Index of current layer\n",
    "        o_i = cache['input'][layer]  # Input of current layer / Output of previous layer\n",
    "        w_j = params['w'][layer]  # Weight of current layer\n",
    "        s_j_der = h.backward(o_i @ w_j)\n",
    "        \n",
    "        if layer == (n - 1):  # If last layer\n",
    "            delta_j = loss.backward(output, t) * s_j_der\n",
    "        else:\n",
    "            delta_l = cache['delta'][layer + 1]  # Delta of next layer\n",
    "            w_l = params['w'][layer + 1]  # Weight of next layer\n",
    "            delta_j = w_l @ delta_l * s_j_der\n",
    "            \n",
    "        grads['w'][layer] = np.outer(o_i, delta_j)\n",
    "        grads['b'][layer] = delta_j\n",
    "        cache['delta'][layer] = delta_j\n",
    "    \n",
    "    # Update parameters\n",
    "    for i in range(n):\n",
    "        params['w'][i] -= learning_rate * grads['w'][i]\n",
    "        params['b'][i] -= learning_rate * grads['b'][i]\n",
    "    \n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At    10 epoch, mean loss of 10 epoch: 69.9127474776679\n",
      "At    20 epoch, mean loss of 10 epoch: 58.11183851043386\n",
      "At    30 epoch, mean loss of 10 epoch: 54.60694051082548\n",
      "At    40 epoch, mean loss of 10 epoch: 53.516423076351465\n",
      "At    50 epoch, mean loss of 10 epoch: 52.569727966237814\n",
      "At    60 epoch, mean loss of 10 epoch: 52.03470276247441\n",
      "At    70 epoch, mean loss of 10 epoch: 51.74193517988944\n",
      "At    80 epoch, mean loss of 10 epoch: 51.47856187251277\n",
      "At    90 epoch, mean loss of 10 epoch: 50.14140025878953\n",
      "At   100 epoch, mean loss of 10 epoch: 49.29044547863606\n",
      "At   110 epoch, mean loss of 10 epoch: 49.1259387078269\n",
      "At   120 epoch, mean loss of 10 epoch: 49.22909755635836\n",
      "At   130 epoch, mean loss of 10 epoch: 49.3607949340306\n",
      "At   140 epoch, mean loss of 10 epoch: 49.49863885067637\n",
      "At   150 epoch, mean loss of 10 epoch: 49.46716564586022\n",
      "At   160 epoch, mean loss of 10 epoch: 49.25086335696626\n",
      "At   170 epoch, mean loss of 10 epoch: 49.07025250473208\n",
      "At   180 epoch, mean loss of 10 epoch: 49.01377746276818\n",
      "At   190 epoch, mean loss of 10 epoch: 48.960626721092474\n",
      "At   200 epoch, mean loss of 10 epoch: 48.9099264281752\n"
     ]
    }
   ],
   "source": [
    "# Solution to exercise 3\n",
    "\n",
    "# Weights for the input layer, hidden layers, and output layers.\n",
    "weights_b = [np.random.rand(2, 10), np.random.rand(10, 2)]\n",
    "biases_b = [np.zeros(10), np.zeros(2)]\n",
    "\n",
    "learning_rate = 0.1\n",
    "\n",
    "epoch = 200  # How many time to train the network with the whole training data.\n",
    "num_e = 10  # Print info every num_e epoch\n",
    "training_loss = 0\n",
    "losses = []\n",
    "\n",
    "for e in range(epoch):\n",
    "    loss = 0\n",
    "    for inputs, targets in zip(training['input'], training['target']):        \n",
    "        outputs = backpropagation_bias(inputs, targets, weights_b, biases_b, ReLU, CrossEntropyLoss, learning_rate)\n",
    "\n",
    "        training_loss += CrossEntropyLoss.forward(outputs, targets) / num_e\n",
    "        loss += CrossEntropyLoss.forward(outputs, targets)\n",
    "    losses.append(loss)\n",
    "\n",
    "    # Print training info every set amount of epoch\n",
    "    if not (e + 1) % num_e:\n",
    "        print(f\"At {e + 1: 5} epoch, mean loss of {num_e} epoch: {training_loss}\")\n",
    "        training_loss = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Network accuracy: 90.0%\n"
     ]
    }
   ],
   "source": [
    "# Solution to exercise 4\n",
    "\n",
    "# Test the network:\n",
    "correct_count = 0\n",
    "total_count = 0\n",
    "for inputs, targets in zip(training['input'], training['target']):\n",
    "    outputs = network_bias(inputs, weights_b, biases_b, ReLU)\n",
    "\n",
    "    total_count += 1\n",
    "    if np.argmax(outputs) == np.argmax(targets):\n",
    "        correct_count += 1\n",
    "\n",
    "accu = correct_count / total_count\n",
    "\n",
    "print(f\"Network accuracy: {accu * 100}%\")"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}