Ex6-Solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 6: Neural Networks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are using again the wine-data from Exercise 4, but now building a neural network for prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data: (4898, 12)\n",
      "First example:\n",
      "Features: [6.600e+00 1.900e-01 4.100e-01 8.900e+00 4.600e-02 5.100e+01 1.690e+02\n",
      " 9.954e-01 3.140e+00 5.700e-01 9.800e+00]\n",
      "Quality: 6.0\n"
     ]
    }
   ],
   "source": [
    "# load all examples from the file\n",
    "data = np.genfromtxt('winequality-white.csv',delimiter=\";\",skip_header=1)\n",
    "\n",
    "print(\"data:\", data.shape)\n",
    "\n",
    "# Prepare for proper training\n",
    "np.random.shuffle(data) # randomly sort examples\n",
    "\n",
    "# take the first 3000 examples for training\n",
    "X_train = data[:3000,:11] # all features except last column\n",
    "y_train = data[:3000,11]  # quality column\n",
    "\n",
    "# and the remaining examples for testing\n",
    "X_test = data[3000:,:11] # all features except last column\n",
    "y_test = data[3000:,11] # quality column\n",
    "\n",
    "print(\"First example:\")\n",
    "print(\"Features:\", X_train[0])\n",
    "print(\"Quality:\", y_train[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework\n",
    "\n",
    "The goal is to mplement the training of a neural network with one input layer, one hidden layer and one output layer using gradient descent. \n",
    "   * Define the matrices and initialise with random values. We need W,b,W' and b'. The shapes should be:\n",
    "     * W: (number of hidden nodes, number of inputs) named `W`\n",
    "     * b: (number of hidden nodes) named `b`\n",
    "     * W': (number of hidden nodes) named `Wp`\n",
    "     * b': (one) named `bp`\n",
    "   * Implement a forward pass of the network as `dnn` (see below)\n",
    "   * Implement a function that uses one example to update the weights using gradient descent. You can follow the `update_weights` skeleton below\n",
    "   * Now you can use the code below (training loop and evaluation) to train the network for multiple epochs. Try to find a set of hyperparameters (number of nodes in the hidden layer, learning rate, number of training epochs) that gives stable results. What is the best result (as measured by the loss on the training sample) you can get?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can use this implementation of the ReLu activation function\n",
    "def relu(x):\n",
    "    return np.maximum(x, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "hidden_nodes = 50\n",
    "\n",
    "W = np.random.randn(hidden_nodes,11)*np.sqrt(2./11)\n",
    "b = np.random.randn(hidden_nodes)*np.sqrt(2./11)\n",
    "Wp = np.random.randn(hidden_nodes)*np.sqrt(2./hidden_nodes)\n",
    "bp = np.random.randn((1))\n",
    "\n",
    "lr = 0.000005\n",
    "\n",
    "\n",
    "train_losses = []\n",
    "test_losses = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dnn(x,W,b,Wp,bp):\n",
    "    # Calculate and return network output of forward pass\n",
    "    return np.dot(Wp, relu(np.dot(W,x) + b)) + bp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_weights(x,y, W, b, Wp, bp):\n",
    "    \n",
    "    # Calculate the network output\n",
    "    \n",
    "    # Use the formulas derived to calculate the gradient for each of W,b,Wp,bp\n",
    "    \n",
    "    # Update the weights/bias following the rule:  X_new = X_old - learning_rate * gradient\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_weights(x,y, W, b, Wp, bp):\n",
    "    \n",
    "    # Calculate the network output\n",
    "    phi = dnn(x,W,b,Wp,bp)\n",
    "\n",
    "    # Use the formulas derived to calculate the gradient for each of W,b,Wp,bp\n",
    "    delta_bp = 2 * (phi - y)\n",
    "    delta_Wp = 2 * (phi - y) * relu(np.dot(W,x) + b)\n",
    "    delta_b  = 2 * (phi - y) * Wp * np.heaviside(np.dot(W,x) + b, 0.5)\n",
    "    delta_W  = 2 * (phi - y) * np.outer(Wp * np.heaviside(np.dot(W,x) + b, 0.5), x)\n",
    "          \n",
    "    # Update the weights/bias following the rule:  X_new = X_old - learning_rate * gradient    \n",
    "    bp -= lr * delta_bp\n",
    "    Wp -= lr * delta_Wp\n",
    "    b  -= lr * delta_b\n",
    "    W  -= lr * delta_W\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training loop and evaluation below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch: 0 Train Loss: 0.6276763491725015 Test Loss: 0.6462291935683923\n",
      "Best loss: 0.5848485857988583 Final loss: 0.6462291935683923\n",
      "Correlation coefficient: 0.5037907844764149\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHiFJREFUeJzt3X+QHGd95/H3d0YtmBU2I8HaQWstwsZZDiMkmYm1RHcUxoAgMc5aJyOUU6VCCnThchwcd5uKgiu2KRMltcHlkKvjYkjlSAEKYOwt8+NYnAsXIIWVW3mxZQfvgcG2vEvwYnsNWIM9mv3eHzOznpmdHz2z2/Nj+/OqUmnV83Q/32e6n8/29vS2zN0REZH1L9HtAkREpDMU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmNnS7gHIvfvGLffv27d0uQ0Skb5w8efIn7j4Ypm1PBf727duZnp7udhkiIn3DzB4O21aXdEREYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYiDXwze5+Z3Wdm95vZ+6PsS0REGovstkwzexXwbuAy4Fngq2b2ZXf/XlR9ikh3TM7MMTE1y/xilq3pFOP7RhjbPdRym16ocz2L8gz/XwF3ufsZdz8L/ANwdYT9iUgXTM7McfS2U8wtZnFgbjHL0dtOMTkz11KbXqhzvYsy8O8DXmdmLzKzAeDXgG0R9iciXTAxNUs2l69Yls3lmZiabalN1Hqhhm6L7JKOu3/XzP4UuBP4OXAPcLa6nZkdAY4ADA8PR1WOiERkfjHbdHmYNlHrhRq6LdIPbd39r9z9Und/HfAEsOL6vbvf4u4Zd88MDoZ6HISI9JCt6VTT5WHaRK0Xaui2qO/SOa/49zCwHzgeZX8i0nnj+0ZIBcmKZakgyfi+kZbaRK0Xaui2qB+e9gUzexGQA37P3Z+MuD8R6bDSXS6N7n4J06YX6lzvzN27XcOyTCbjelqmiEh4ZnbS3TNh2uo3bUVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISExsiHLjZvafgXcBDpwC3unuv4iyT1m9PR++kx//7NmKZec+L8m9N7yl47VMzswxMTXL/GKWrekU4/tGGNs9VLftDV+8nyfP5ABIpwKuv+qSuu1b6RNgYmqWucXscjsDEgZ5f+7fqSBBNrdEeiDgF7k82dwSAEECzi4VJkL5+l71dWn9RsrXK0kY/OaeYQCOnzhN3p2kGYf2bOPGsR0rxvfB20/x9LP5iuV7L9rCNZnhFeMsqfV+lt6rucUsSTPy7mweCHgml+dMnXEkDJaqBwBsTBrP5mu8UEet96FeveX7tHrfQOF9f36QZPFMrulxVq6V47MXmHv4N7ilDZsNAd8CXunuWTP7HPAVd/+f9dbJZDI+PT0dST0STq2wL+l06E/OzHH0tlNkc88FUypIcmz/jhWTanJmjvFb7yFXFRhBwpi4ZmfoSVirzyBp4JCrlVJ94PDo8HLoT87M8V8+fw/5NsdS/n7Weq96TZAwDl62jS+cnGupznrHWblWjs8omdlJd8+EaRv1JZ0NQMrMNgADwHzE/ckq1Qt7gJ8+09mJPTE1u2KSZnN5JqZma7atDnsohHSt9q30mct734Y9FM74SyamZtsOe6h8P2u9V70mt+QcP3G65TrrHWflWjk+e0Vkge/uc8CfAY8APwKecvevVbczsyNmNm1m0wsLC1GVI31ovsZlhXrL67Vt9tpq2vaLfNlP8WsxvtI2+uW9yrd5FaPZ+Fo5PntFZIFvZpuB3wBeBmwFNpnZ4ep27n6Lu2fcPTM4OBhVOdKHtqZToZfXa9vstdW07RdJs+Wv12J8pW30y3tVPv5WNBtfK8dnr4jyks4bgR+6+4K754DbgF+NsD9ZA+efs7Hua+c+L9nBSmB83wipoLLPVJBc/hC1um2QXDmxg4TVbN9Kn0HSCBLthUYvOLRn2/LX4/tGSK5iLOXvZ633qtcEicIH163WWe84K9fK8dkrogz8R4BRMxswMwOuAL4bYX+yBk588E01Q78bd+mM7R7i2P4dDKVTGDCUTtX9QGxs9xATB3ayeSBYXpZOBS19YFuvz4kDO5m4ZidDVWduBpR/jzFgIEhgwOaBgFTw3PQKEoXXq9ev/rp8nXpqxXXCCh/OHh4dXj6jTZpVfGBbGt9HrtnJpo0rA3DvRVu4+eCuFeMsqX4/y9+rUn9QGPtAg3HU+36zscY37EaatS7Ve+PYjop9Wr1voPC+bx4Imh5n5Vo5PntFZHfpAJjZDcBB4CwwA7zL3Z+p11536YiItKaVu3QivQ/f3a8DrouyDxERCUe/aSsiEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhMKPBFRGJiQ1QbNrMR4LNliy4E/sjdb46qzzibnJljYmqWucXs8rKkGaMXbuahx7PML2Z5fpDgmbNLLHnluheft4k7P/B6AF5+9MucrXq92lA6xcDGBN977OmW60waLAHuhfoO7dkGwPETp8l7k45F1tjmgYDr3nYJY7uHlufQ/GKWrekU4/tGAFbMq3KluVNv/h3as40bx3bUXLdWf2O7h9Z+kGXMOzDJzCwJzAF73P3heu0ymYxPT09HXs96Mzkzx9HbTpHN5dvexsXnbeKHC083DXuR9SZIGgd/ZRtfODlXMYeChIFBLt94Upx/zkZ++ot83fl3eHR4RejXmrOpIMmx/TtaDn0zO+numTBtO3VJ5wrgwUZhL+2bmJpdVdgDfO8xhb3EUy7vHD9xesUcyi1507AH+PHPnm04/46fOL1iWa05m83lmZiaDVl1ezoV+O8Ajtd6wcyOmNm0mU0vLCx0qJz1Zb7Oj5siEk6UlxNrbbvenI16Lkce+Ga2EbgK+Hyt1939FnfPuHtmcHAw6nLWpa3pVLdLEOlrSbOObrvenI16LnfiDP+twN3u/uMO9BVL4/tGSAXJVW3j4vM2sSG6Y16kZwXJwoer1XMoSBhBsvmkOP+cjQ3nX+nGhHK15mwqSC5/UByVTgT+IepczpG1MbZ7iGP7dzBUdXaQNGPvRVsYSqcwIBUkSNQ4fkt3Gnz/2K+HCv2hdIqLz9vUVq1Jg9IJT9KMw6PDHB4djvQMS6SezQMBEwd2cuPYjuU5ZBSO8YlrdjJxYOeKeVXu4vM2ceKDb6o7/2p9YAuVc7bUXzsf2LYq0rt0zGwAOA1c6O5PNWuvu3RERFrTyl06kd2HD+DuZ4AXRdmHiIiEo9+0FRGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYmJ0IFvZs+LshAREYlW08A3s8vM7BTwveK/d5rZX0RemYiIrKkwZ/gfBa4EHgdw93uAy6MsSkRE1l6YwE+4+8NVy/JRFCMiItHZEKLNaTO7DHAzSwLvBf5ftGWJiMhaC3OG/x7gA8Aw8GNgtLhMRET6SNMzfHd/DHhHOxs3szTwCeBVgAO/4+7fbmdbIiKyOk0D38w+TiGsK7j7kRDb/3Pgq+5+wMw2AgOtl9ieyZk5JqZmmV/MsjWdYnzfCGO7h9a0j2snT/Hpux5ZfnM2bUzy4at3rOgnbC3XTp7i+InT5H3F282mjUl2bXsh3/7BEyytfFlEymzamOTqS4f4+gMLzC9m2ZCA3NLKdulUgBksnsm1lBPlczVpxqE927hxbEcEI1lb5jXCpaKB2cGyfz4fuBo47e7vbbLeucA9wIXerJOiTCbj09PTYZo2NDkzx9HbTpHNPffZcipIcmz/yjBu17WTp/jUXY+sWJ5MGB+5ZudyP2Frqbc9EemcMDlRb64eHh3uSuib2Ul3z4Rp2/Qavrt/tuzPJ4H9wCtDbPtCYAH4azObMbNPmNmmMEWt1sTUbEXAAmRzeSamZtesj+MnTtdcnl/yin7C1lJveyLSOWFyot5c7Yc53M6jFV4GvDREuw3ApcDH3H038DTwB9WNzOyImU2b2fTCwkIb5aw0v5htaXk7al12qdVP2FoabU9EOqdZTtSbq/0wh8P8pu2TZvZE8c8icCfwhyG2/SjwqLufKP77VgrfACq4+y3unnH3zODgYCu117U1nWppeTuSZqH6D1tLo+2JSOc0y4l6c7Uf5nDDwDczA3YCg8U/m939Qnf/XLMNu/u/ULiHf6S46Argn1dZbyjj+0ZIBcmKZakgyfi+kTprtO7Qnm01lycTVtFP2FrqbU9EOidMTtSbq/0whxsGfvHD1tvdPV/80+rPLO8FPm1m9wK7gD9us86WjO0e4tj+HQylUxgwlE6t6Qe2ADeO7eDw6DDl39M3bUxWfGDbSi2l7dU7S9i0Mcnei7aQ6P2TCJGu27QxyeHR4eV5F9RJunQqYPNA0FJOVM/VpFnXPrBtVZi7dD4GfNzd7466mLW6S0dEJC5auUun7n34ZrbB3c8C/xp4t5k9SOGDV6Nw8r/ieryIiPSuRr949U8UPmQd61AtIiISoUaBbwDu/mCHahERkQg1CvxBM/tAvRfd/aYI6hERkYg0Cvwk8AJA94WIiKwDjQL/R+7+oY5VIiIikWp0H77O7EVE1pFGgX9Fx6oQEZHI1Q18d3+ik4WIiEi02nlapoiI9CEFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYa/ReHq2ZmDwE/A/LAWXfPRNmfiIjUF2ngF13u7j+JsoPJmTkmpmaZW8ySNCPvzlA6xfi+EcZ2D1W0mV/MsrX4GsANX7yfJ8/klrdlgAPpVMCZZ8/ybN4Lyw3+3Z5hMi/dwvV33M9itrBOwmDJqej38lcM8vUHFphfzJJMGGeXPMrhS4wNBAn+eP+rmX74CY6fOE3ewx1rmzYmufrSIb5874+Wj/90KuD6qy4BqDjGzcCdunNqbjG7PG8ANg8EXPe2SypqSppxaM82Mi/dsmKupoIEz5xdWp5H5e3K52upX2mfecgDpK2NF87wM2EDP5PJ+PT0dEt9TM7McfS2U2Rz+RWvpYIkx/bvAFjRJkgYS0BeYSyyLAFQPImppdGcqthOnW0kgKWQdZS3K/Wr0F/JzE6GvXoS9Rm+A18zMwf+0t1vWesOJqZm6x502VyeianZ5a/L5RT0IisswXOn6jU0mlMV26mzjTBhX6tdqV8F/upEHfh73X3ezM4D7jSzB9z9G+UNzOwIcARgeHi45Q7mF7Orel1EWtOtOaW5vHqR3qXj7vPFvx8Dbgcuq9HmFnfPuHtmcHCw5T62plNNX2/WRkTC69ac0jxevcgC38w2mdk5pa+BNwP3rXU/4/tGSAXJmq+lgiTj+0ZqtgkSRjJha12OSF9LULj+Xk+jOVWxnTrbCBs41e1K/crqRHmGfz7wLTO7B/gn4Mvu/tW17mRs9xDH9u9gqPjdP2mFI20onVr+kKe8jRVfm7hmJx+5ZiebB4KK7ZWO03QqYGPyuaPWDA6PDnPzwV2kU8+tUzqwy/s9PDq83NcGfVORCA0ECW4+uIvDo8PLx2AYmzYmOTw6XHH8p1MBNx3cxU1vrzzGS5utN6fguXkDhbt0bnp7ZU1JMw6PDnPTwV0r5moqSFTMo/J2VtWvrE6kd+m0qp27dERE4qyVu3T0m7YiIjGhwBcRiQkFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmIg98M0ua2YyZfSnqvkREpL4NHejjfcB3gXOj7GRyZo4bvng/T57JAZBOBVy58yV86Z4fsZgtLNs8EPDKl5zDXT94krw7STMO7dnGjWM7VmxrYmqW+cUsW9Mptr8otbyOGZjDUrFtwuA39wwDcPzE6YrtZl66ZXk76YGAn/8iR26pois2DwS4w1PZHFvTKS5/xSBff2Bhue/xfSMAXH/H/cvjaGRDwji75E3bpVMB1191CWO7h5icmeMPb7uXM9XFrUE/JQnghQMBi2dyy+Ma2z20/Pq1k6f49IlH8PCbXJa0wv5wh6QZoxdu5qHHs8wtZkmakXcnYVBdbvn+r97nteor7d+S8n33wlSAGRXjAyq2WWvflvfRbc3eg6jWlc4xb2eGhd242QXAJ4EPAx9w9ysbtc9kMj49Pd1yP5Mzc4zfeg+5fHtjOTw6vBz6kzNzHL3tFNlcvq1tlasVMq0KkkY+74SP4ha2nTAOXraNz5x4ZNV1tioVJDm2fwdju4e4dvIUn7rrkc4WUGbvRVu4+5GnKvb5ausLkgYOuQZvbHkf3VbruA9b32rWldUzs5PungnTNupLOjcDvw+R5NWyianZtsMeCmfm5dtai7CH1Yc9QC6isIdCGB0/cbrjYQ+QzeWZmJoFKt//bvjHB59Ysc9XW18u7w3DvrqPbqt13IetbzXrSmdFFvhmdiXwmLufbNLuiJlNm9n0wsJCW33NL2bbWq+k/Mf01W6r3+Qj/AmvmdJ73c0aGulEfb1yvNWrI0x9q1lXOivKM/y9wFVm9hDwt8AbzOxT1Y3c/RZ3z7h7ZnBwsK2OtqZTqyo0abZm2+o35WPvtNJ73c0aGulEfb1yvNWrI0x9q1lXOiuywHf3o+5+gbtvB94B/L27H46ir/F9I4Vrpm06tGdbxbZSQXItyiKxBjkRJC2ynRQkCh9arkWdrUoFyeUPNsvf/27Ye9GWFft8tfUFSSNo8saW99FttY77sPWtZl3prHVxH/7Y7iEmDuxk80CwvCydCjg8Okw69dyyzQMBey/asnzGljSr+MC2tK1j+3cwlE5hwFA6VbGOWeWblrDCh76HR4dXbPemt+9a3s7mgYCgxru9eSAgnQqW+zo8OlzR98SBndx0cFfFOBrZEDK906mAiWt2cuPYDm56+y4GahW3Bv2UJCiMtTSu8g/0bhzbweHRYdo9kU4ay+smzdh70RaGqs7Oa5Vb2k+ffvdrV+zzWvVVn+mX77t0KqgY38SBnUxcs7Nim9X7tpc+1Kx13IetbzXrSmdFepdOq9q9S0dEJK566S4dERHpEQp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYmJDVBs2s+cD3wCeV+znVne/Lqr+SiZn5piYmmV+McvWdIrxfSOM7R6KutuW+y61nVvMYoAXl28eCLjubZesWK+b42qmUW3Vr13+ikG+/sBC3X/XG1cvj78d62080h/M3Zu3amfDZgZscvefm1kAfAt4n7vfVW+dTCbj09PTbfc5OTPH0dtOkc3ll5elgiTH9u+IfDK10nettuWCpDFxYGdFaHZrXM00qg1oOM5aao2rl8ffjvU2HukuMzvp7pkwbSO7pOMFPy/+Myj+iea7S9HE1OyKcMnm8kxMzUbZbct912pbLpf3ivW6Oa5mGtXWbJy11BpXL4+/HettPNI/Ir2Gb2ZJM/sO8Bhwp7ufqNHmiJlNm9n0wsLCqvqbX8y2tHwttdJ3mHrK23RzXM00qq3d+qrX6+Xxt2O9jUf6R6SB7+55d98FXABcZmavqtHmFnfPuHtmcHBwVf1tTadaWr6WWuk7TD3lbbo5rmYa1dZufdXr9fL427HexiP9oyN36bj7IvB/gLdE2c/4vhFSQbJiWSpIMr5vJMpuW+67VttyQdIq1uvmuJppVFuzcdZSa1y9PP52rLfxSP+I8i6dQSDn7otmlgLeCPxpVP0Byx94dePuh1b6Lm8b5i6dbo6rmTC1rfYunV4efzvW23ikf0R5l86rgU8CSQo/SXzO3T/UaJ3V3qUjIhI3rdylE9kZvrvfC+yOavsiItIa/aatiEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jERGS3ZbbDzBaAh7tdRwteDPyk20WsIY2nt2k8va1b43mpu4d6TEFPBX6/MbPpsPe/9gONp7dpPL2tH8ajSzoiIjGhwBcRiQkF/urc0u0C1pjG09s0nt7W8+PRNXwRkZjQGb6ISEwo8EMq/u9dM2b2pRqv/baZLZjZd4p/3tWNGsMys4fM7FSx1hWPJ7WCj5rZ983sXjO7tBt1hhViPK83s6fK9s8fdaPOsMwsbWa3mtkDZvZdM3tt1ev9tn+ajadv9o+ZjZTV+R0z+6mZvb+qTc/un8ielrkOvQ/4LnBundc/6+7/sYP1rNbl7l7vnuG3AhcX/+wBPlb8u5c1Gg/AN939yo5Vszp/DnzV3Q+Y2UZgoOr1fts/zcYDfbJ/3H0W2AWFk0BgDri9qlnP7h+d4YdgZhcAvw58otu1dMhvAH9T/I/o7wLSZvaSbhcVB2Z2LvA64K8A3P3Z4v8YV65v9k/I8fSrK4AH3b36l0V7dv8o8MO5Gfh9YKlBm39b/PHtVjPb1qG62uXA18zspJkdqfH6EHC67N+PFpf1qmbjAXitmd1jZv/LzC7pZHEtuhBYAP66eAnxE2a2qapNP+2fMOOB/tk/5d4BHK+xvGf3jwK/CTO7EnjM3U82aPZFYLu7vxr4Owr/01cv2+vul1L40fP3zOx1Va9bjXV6+XauZuO5m8Kvn+8E/gKY7HSBLdgAXAp8zN13A08Df1DVpp/2T5jx9NP+AaB4aeoq4PO1Xq6xrCf2jwK/ub3AVWb2EPC3wBvM7FPlDdz9cXd/pvjPjwOv6WyJrXH3+eLfj1G4/nhZVZNHgfKfUi4A5jtTXeuajcfdf+ruPy9+/RUgMLMXd7zQcB4FHnX3E8V/30ohMKvb9Mv+aTqePts/JW8F7nb3H9d4rWf3jwK/CXc/6u4XuPt2Cj/C/b27Hy5vU3V97ioKH+72JDPbZGbnlL4G3gzcV9XsDuC3incbjAJPufuPOlxqKGHGY2a/ZGZW/PoyCsf9452uNQx3/xfgtJmNFBddAfxzVbO+2T9hxtNP+6fMIWpfzoEe3j+6S6dNZvYhYNrd7wD+k5ldBZwFngB+u5u1NXE+cHtxfm0APuPuXzWz3wVw9/8BfAX4NeD7wBngnV2qNYww4zkAvMfMzgJZ4B3e279x+F7g08XLBj8A3tnH+weaj6ev9o+ZDQBvAv592bK+2D/6TVsRkZjQJR0RkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6sG2aWLz7B8D4z+3zx9rl2t/V6Kz4Z1cyuMrPq3w4tb5s2s//QRh/Xm9l/bbdGkVYp8GU9ybr7Lnd/FfAs8LvlLxZ/EablY97d73D3P2nQJA20HPginabAl/Xqm8DLzWx78Rns/53CM1u2mdmbzezbZnZ38SeBFwCY2Vus8Mz2bwH7Sxuywv938N+KX59vZrcXH/R1j5n9KvAnwEXFny4miu3Gzez/Fh+od0PZtj5oZrNm9nfACCIdpMCXdcfMNlB41smp4qIRCo+rLT2861rgjcUHrk0DHzCz51N4DtLbgH8D/FKdzX8U+Ifig74uBe6n8DCwB4s/XYyb2ZspPAv9MgrPTn+Nmb3OzF5D4fEcuyl8Q/mVNR66SEN6tIKsJykz+07x629SeAb7VuDh4nPJAUaBVwL/WHwcw0bg28ArgB+6+/cAig/Iq/Wo5TcAvwXg7nngKTPbXNXmzcU/M8V/v4DCN4BzgNvd/UyxjztWNVqRFinwZT3Juvuu8gXFUH+6fBFwp7sfqmq3i7V7hK0Bx9z9L6v6eP8a9iHSMl3Skbi5C9hrZi+HwoOwzOyXgQeAl5nZRcV2h+qs/7+B9xTXTVrhf3T6GYWz95Ip4HfKPhsYMrPzgG8AV5tZqviEz7et8dhEGlLgS6y4+wKFp5keN7N7KXwDeIW7/4LCJZwvFz+0rf5v60reB1xuZqeAk8Al7v44hUtE95nZhLt/DfgM8O1iu1uBc9z9buCzwHeAL1C47CTSMXpapohITOgMX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMTE/wcsh3S0++jk1wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xmc3HWd5/HXp+6u7k53ku4c5CAJxHCGq+UQRwgjR3AElRkGFMVr4ujMqruDK7o7sOrOrjPuOiw7CoNjZFglKpeiohIwiMplggECSchhIE1CunP33V1Vn/3j9+tOdac6qSRdXZ3U+/l4/B71+31/v6r69A/ye9fv+v7M3RERETmYSLkLEBGRo4MCQ0REiqLAEBGRoigwRESkKAoMEREpigJDRESKosAQEZGiKDBERKQoCgwRESlKrNwFjKSGhgafNWtWucsQETlqrFixYru7Nxaz7DEVGLNmzWL58uXlLkNE5KhhZq8Vu6wOSYmISFEUGCIiUhQFhoiIFOWYOochInKo+vr6aG5upru7u9yllFQqlWL69OnE4/HD/gwFhohUtObmZmpra5k1axZmVu5ySsLd2bFjB83NzcyePfuwP0eHpESkonV3dzNx4sRjNiwAzIyJEyce8V6UAkNEKt6xHBb9RuJvVGAA/PqfYP1j5a5CRGRMU2AA/O7/wPrHy12FiFSg3bt3881vfvOQ33fllVeye/fuElQ0PAUGQKIGetvLXYWIVKDhAiObzR7wfY888gj19fWlKqsgXSUF9ESqyLTvpbrchYhIxbn55pvZsGEDZ555JvF4nJqaGqZOncrKlSt55ZVXeM973sPmzZvp7u7mM5/5DIsWLQL2dYXU3t7OwoULefvb385TTz3FtGnT+PGPf0xVVdWI16rAANbvcZLWyonlLkREyupLP3mZV7bsHdHPPOW4cdz67lOHnf/Vr36VVatWsXLlSp544gne9a53sWrVqoHLXxcvXsyECRPo6urirW99K9dccw0TJ04c9Bnr1q1jyZIlfOtb3+Laa6/lgQce4IYbbhjRvwMUGAB0W5p0pqPcZYiIcO655w66V+L222/noYceAmDz5s2sW7duv8CYPXs2Z555JgDnnHMOmzZtKkltCgyCQ1KxrM5hiFS6A+0JjJbq6n0Hx5944gkee+wxnn76adLpNBdffHHBeymSyeTAeDQapaurqyS16aQ30BetIp7tLHcZIlKBamtraWtrKzhvz549jB8/nnQ6zZo1a3jmmWdGubrBtIcB9EbTJPoUGCIy+iZOnMiFF17IaaedRlVVFZMnTx6Yd8UVV3DnnXcyf/585s2bx/nnn1/GShUYAGRiaZI9pdmFExE5mHvvvbdgezKZ5Oc//3nBef3nKRoaGli1atVA+0033TTi9fXTISkgF6sm5d3gXu5SRETGrJIFhpnNMLNlZrbazF42s88UWOYDZvZiODxlZmfkzdtkZi+Z2UozK+lzV7PxaqLkoE97GSIiwynlIakM8Hfu/ryZ1QIrzGypu7+St8wfgYvcfZeZLQTuAs7Lm7/A3beXsEYAPB5eldDbAYl0qb9OROSoVLI9DHff6u7Ph+NtwGpg2pBlnnL3XeHkM8D0UtVzIJboDwxdWisiMpxROYdhZrOAs4BnD7DYx4D8szsOPGpmK8xsUemqA5K1wRf2FL60TURERuEqKTOrAR4APuvuBe+5N7MFBIHx9rzmC919i5lNApaa2Rp3f7LAexcBiwBmzpx5WDVGUjUA9Ha2kTzIsiIilaqkexhmFicIi++5+4PDLDMf+Dfganff0d/u7lvC1xbgIeDcQu9397vcvcndmxobGw+rzmgyCIzuzpHtQ0ZE5GAOt3tzgNtuu43OztG7h6yUV0kZ8G1gtbt/fZhlZgIPAh9091fz2qvDE+WYWTVwGbCq0GeMhGjVOCDYwxARGU1HU2CU8pDUhcAHgZfMbGXY9kVgJoC73wncAkwEvhk+PjDj7k3AZOChsC0G3OvuvyhVoYmq4BxGpkt7GCIyuvK7N7/00kuZNGkSP/zhD+np6eG9730vX/rSl+jo6ODaa6+lubmZbDbL3//937Nt2za2bNnCggULaGhoYNmyZSWvtWSB4e6/BQ74EFl3/zjw8QLtG4Ez9n9HacSrgkNSfV26Skqkov38ZnjzpZH9zCmnw8KvDjs7v3vzRx99lPvvv5/nnnsOd+eqq67iySefpLW1leOOO46f/exnQNDHVF1dHV//+tdZtmwZDQ0NI1vzMHSnN5CsDp5ale3WHoaIlM+jjz7Ko48+yllnncXZZ5/NmjVrWLduHaeffjqPPfYYn//85/nNb35DXV1dWepTX1JAVVUVGY+Q69EzMUQq2gH2BEaDu/OFL3yBT3ziE/vNW7FiBY888ghf+MIXuOyyy7jllltGvT7tYQDpZJxOUrhu3BORUZbfvfnll1/O4sWLaW8PtkVvvPEGLS0tbNmyhXQ6zQ033MBNN93E888/v997R4P2MIDqRJR2UtCjwBCR0ZXfvfnChQt5//vfzwUXXABATU0N3/3ud1m/fj2f+9zniEQixONx7rjjDgAWLVrEwoULmTp16qic9DY/hnpobWpq8uXLD72fwrbuPrb9j/lEppzKnE/dX4LKRGSsWr16NSeffHK5yxgVhf5WM1sRXp16UDokBaQTMdpJYX06hyEiMhwFBhCNGN1WRUyBISIyLAVGqMeqiOq53iIV6Vg6ND+ckfgbFRihnmiaeEaBIVJpUqkUO3bsOKZDw93ZsWMHqVTqiD5HV0mF+qJVJLSHIVJxpk+fTnNzM62treUupaRSqRTTpx/ZI4cUGKFMNE1Sj2gVqTjxeJzZs2eXu4yjgg5JhTKxalLeDblcuUsRERmTFBihbP9zvXWllIhIQQqMfvF08NqrwBARKUSBEfJE8EwMBYaISGEKjFAkGR6S6tFT90REClFghCKp4CFKmW4FhohIIQqMUCwVPNe7R8/1FhEpqGSBYWYzzGyZma02s5fN7DMFljEzu93M1pvZi2Z2dt68G81sXTjcWKo6+8XC53r3dOwp9VeJiByVSnnjXgb4O3d/3sxqgRVmttTdX8lbZiEwNxzOA+4AzjOzCcCtQBPg4XsfdvddpSo2kQ4Co69LexgiIoWUbA/D3be6+/PheBuwGpg2ZLGrgXs88AxQb2ZTgcuBpe6+MwyJpcAVpaoVIFUdPCNXgSEiUtionMMws1nAWcCzQ2ZNAzbnTTeHbcO1l0wq3MPIKDBERAoqeWCYWQ3wAPBZd987dHaBt/gB2gt9/iIzW25my4+k87B0uooej5HTY1pFRAoqaWCYWZwgLL7n7g8WWKQZmJE3PR3YcoD2/bj7Xe7e5O5NjY2Nh11rTTJGBykFhojIMEp5lZQB3wZWu/vXh1nsYeBD4dVS5wN73H0r8EvgMjMbb2bjgcvCtpKpScboJAUKDBGRgkp5ldSFwAeBl8xsZdj2RWAmgLvfCTwCXAmsBzqBj4TzdprZV4Dfh+/7srvvLGGtVCdjbPeUHtMqIjKMkgWGu/+Wwuci8pdx4G+GmbcYWFyC0gpKxiJ0UkW9AkNEpCDd6R0yM7ojKaIZBYaISCEKjDy9kTQxPddbRKQgBUaevmiaeFaPaRURKUSBkScTS5PMaQ9DRKQQBUaebCxNMqc9DBGRQhQYebLxGpL0QjZT7lJERMYcBUa+ePjUvV7dvCciMpQCI1//Y1r1XG8Rkf0oMPJYMuix1rWHISKyHwVGnmgyeK53T8fQTnVFRESBkSeaDp7r3dWxu8yViIiMPQqMPLF0PQC9bSV7EqyIyFFLgZEnVj0BgExHSTvGFRE5Kikw8iRrJwIKDBGRQhQYeWpqxtHjMbIKDBGR/Sgw8tSlE+ylGu/SSW8RkaEUGHnq03F2ew3WpZPeIiJDKTDy1CRj7KGGaI/2MEREhlJg5DEzOiM1xPt0456IyFAlCwwzW2xmLWa2apj5nzOzleGwysyyZjYhnLfJzF4K5y0vVY2FdMXGkezbM5pfKSJyVCjlHsbdwBXDzXT3r7n7me5+JvAF4Nfunn950oJwflMJa9xPb3wc6WzbaH6liMhRoWSB4e5PAsVen3o9sKRUtRyKvkQdVd4J2b5ylyIiMqaU/RyGmaUJ9kQeyGt24FEzW2Fmiw7y/kVmttzMlre2th5xPdlk0D0I3TosJSKSr+yBAbwb+N2Qw1EXuvvZwELgb8zsHcO92d3vcvcmd29qbGw84mKsKgwMXVorIjLIWAiM6xhyOMrdt4SvLcBDwLmjVYylg/6kcp0KDBGRfGUNDDOrAy4CfpzXVm1mtf3jwGVAwSutSqG/A8KuvdtH6ytFRI4KsVJ9sJktAS4GGsysGbgViAO4+53hYu8FHnX3/GeiTgYeMrP++u5191+Uqs6hEjVBYHTv3UH1aH2piMhRoGSB4e7XF7HM3QSX3+a3bQTOKE1VB1dV1wBAT9uOcpUgIjImjYVzGGNKelzQxXlfu3qsFRHJp8AYoq66ir1eRbZTgSEikk+BMUR9Os4e9VgrIrIfBcYQdVVxdlFDpFuBISKST4ExRCoeZQ+1xHsUGCIi+RQYBXRGx5FQj7UiIoMoMAroiteTzughSiIi+RQYBfQl6knnOiCbKXcpIiJjhgKjgExqfDCiK6VERAYoMAqpCroHoUv3YoiI9FNgFNDfY613qnsQEZF+CowC4rVB9yDd6rFWRGSAAqOAZG3wIKau3S1lrkREZOxQYBRQVR8ERo/2MEREBigwChg3rp4ej5Fp1zkMEZF+CowCxlcn2E0NOZ30FhEZUFRgmNkJZpYMxy82s0+bWX1pSyufuqoEu7wW02W1IiIDit3DeADImtmJwLeB2cC9JauqzOrTcXZTQ1Q91oqIDCg2MHLuniF4Bvdt7v4fgakHeoOZLTazFjNbNcz8i81sj5mtDIdb8uZdYWZrzWy9md1c7B8zUuLRCG1WS6JX/UmJiPQrNjD6zOx64Ebgp2Fb/CDvuRu44iDL/MbdzwyHLwOYWRT4BrAQOAW43sxOKbLOEdMVqyOlHmtFRAYUGxgfAS4A/sHd/2hms4HvHugN7v4kcDgnAc4F1rv7RnfvBb4PXH0Yn3NEuhPjSWf3QC472l8tIjImFRUY7v6Ku3/a3ZeY2Xig1t2/OgLff4GZvWBmPzezU8O2acDmvGWaw7ZR1ZWaRJQcdOheDBERKP4qqSfMbJyZTQBeAL5jZl8/wu9+Hjje3c8A/i/wo/6vK7CsH6C2RWa23MyWt7a2HmFJ+/SlJwcjbVtG7DNFRI5mxR6SqnP3vcD7gO+4+znAO4/ki919r7u3h+OPAHEzayDYo5iRt+h0YNittrvf5e5N7t7U2Nh4JCUNkqsJz+m3vTlinykicjQrNjBiZjYVuJZ9J72PiJlNMTMLx88Na9kB/B6Ya2azzSwBXAc8PBLfeShidccB0LfrjdH+ahGRMSlW5HJfBn4J/M7df29mc4B1B3qDmS0BLgYazKwZuJXwyip3vxP4c+CTZpYBuoDr3N2BjJn9bfh9UWCxu798yH/ZEaptOI6sG507NlM32l8uIjIGFRUY7n4fcF/e9EbgmoO85/qDzP8X4F+GmfcI8EgxtZXKpPoatlNHZLfOYYiIQPEnvaeb2UPhjXjbzOwBM5te6uLKafK4JG/6BHyPAkNEBIo/h/EdgvMIxxFc4vqTsO2YNWVcim0+nljntnKXIiIyJhQbGI3u/h13z4TD3cDIXZI0BtVVxWm1CVR1KTBERKD4wNhuZjeYWTQcbiC4oumYZWZ0JBupyu6Fvu5ylyMiUnbFBsZHCS6pfRPYSnCF00dKVdRY0VfVf/Pe1vIWIiIyBhTbNcjr7n6Vuze6+yR3fw/BTXzHtFytbt4TEel3JE/c+08jVsUYFa0LLwTbq5v3RESOJDAK9fl0TEk2HA9Az/aNZa5ERKT8jiQwhu0Q8FjRMGE8rV5Hd8sfy12KiEjZHfBObzNro3AwGFBVkorGkEm1KV73SZywS4EhInLAwHD32tEqZCyaWpdipTcyb+/r5S5FRKTsjuSQ1DHvuPoqNvsk0p1bIZspdzkiImWlwDiARCxCW2oaEbKwt7nc5YiIlJUC4yAydTODkV2bylqHiEi5KTAOIjZxdjCy67XyFiIiUmYKjIOom3w8fR4ls0NXSolIZVNgHMT0ibU0ewNdLRvKXYqISFkpMA5ixoQ0G/042L623KWIiJRVyQLDzBaHT+hbNcz8D5jZi+HwlJmdkTdvk5m9ZGYrzWx5qWosxswJadb5dNJ7NurSWhGpaKXcw7gbuOIA8/8IXOTu84GvAHcNmb/A3c9096YS1VeUidUJNkVmEPU+XSklIhWtZIHh7k8COw8w/yl33xVOPgOMyWeEmxkd404MJlpXl7cYEZEyGivnMD4G/Dxv2oFHzWyFmS0qU00D4lNOCkZa1pS3EBGRMjpgX1KjwcwWEATG2/OaL3T3LWY2CVhqZmvCPZZC718ELAKYOXNmSWo8fsokNq9rZOq21eVfYSIiZVLWPQwzmw/8G3C1uw88I9zdt4SvLcBDwLnDfYa73+XuTe7e1NjYWJI63zK5hldz0+l785WSfL6IyNGgbIFhZjOBB4EPuvuree3VZlbbPw5cBhS80mq0zJ1cwzqfRmL3el0pJSIVq2RHWMxsCXAx0GBmzcCtQBzA3e8EbgEmAt80M4BMeEXUZOChsC0G3OvuvyhVncU4fmI1rzKLaK4PWtfAlNPKWY6ISFmULDDc/fqDzP848PEC7RuBM/Z/R/nEoxF2158G7cCWPygwRKQijZWrpMa86qlzaSMdBIaISAVSYBRp7pQ6XszNJtu8otyliIiUhQKjSKdNG8dLuTlYy8uQ6Sl3OSIio06BUaTTp9XzYm42kVwfbHu53OWIiIw6BUaRGmuTvFlzSjDxhg5LiUjlUWAcgsbpJ9JqE+G1p8pdiojIqFNgHIL5M8bzVGYeuU2/A/dylyMiMqoUGIdg/vQ6nsudRKRjG+zcWO5yRERGlQLjEMyfXs+zfnIwocNSIlJhFBiHoK4qTnzSSeyN1MFrvyt3OSIio0qBcYjOmzOR32VOxjc+ofMYIlJRFBiH6Pw5E3g8ewbWthW2lbUTXRGRUaXAOETnzp7Ir7Nh34jrHi1vMSIio0iBcYgmVCcYP3k6f4yfCOuWlrscEZFRo8A4DO+Y28gj3afjm5+Fzp3lLkdEZFQoMA7DJSdN4ueZczDPwZqflrscEZFRocA4DE2zJrApMZedianw8o/KXY6IyKhQYByGRCzCO97SyE8y5+F//LUOS4lIRVBgHKZLTprMfV1NWC4Dq39S7nJEREqupIFhZovNrMXMCt6wYIHbzWy9mb1oZmfnzbvRzNaFw42lrPNwXHryZNZG5rA9NRNeWFLuckRESq7Uexh3A1ccYP5CYG44LALuADCzCcCtwHnAucCtZja+pJUeorp0nHfMncQP+i6C15+G7evLXZKISEmVNDDc/UngQAf4rwbu8cAzQL2ZTQUuB5a6+0533wUs5cDBUxbvmj+VuzvOxy0KK79b7nJEREqq3OcwpgGb86abw7bh2vdjZovMbLmZLW9tbS1ZoYVcespk9sYmsqb2AvjDd6Gve1S/X0RkNJU7MKxAmx+gff9G97vcvcndmxobG0e0uIOpTcVZeNoUvrb3EuhohRd/MKrfLyIymsodGM3AjLzp6cCWA7SPOX/51pn8qnseu+pOgae/AblcuUsSESmJcgfGw8CHwqulzgf2uPtW4JfAZWY2PjzZfVnYNuacP2cCsyZWcw9/BtvXwnr1LyUix6ZSX1a7BHgamGdmzWb2MTP7azP763CRR4CNwHrgW8CnANx9J/AV4Pfh8OWwbcwxM244/3j+77bT6K0+Dn57m56TISLHJPNjaOPW1NTky5cvH/Xvbe/J8Lb/+ThfbPgt122/Hd5/H7zlslGvQ0TkUJnZCndvKmbZch+SOibUJGN88ILj+fs3muirmwNLb4FcttxliYiMKAXGCPnw22Zj0QQ/qP8otK6GlfeWuyQRkRGlwBghjbVJrm2azpfWz6Fnyjnwq69A165ylyUiMmIUGCPo05fMJRaJclvyr6FjO/zii+UuSURkxCgwRtCkcSk+/iezuWNtNVvnfwpeuBdeHZNXA4uIHDIFxgj7xEUnMGVcik+8dgneeBI8/Olgb0NE5CinwBhhNckYt7z7FF58s4sHZv234DzGg4t0B7iIHPUUGCWw8LQpvPPkyXzxKdj6tlthw+Pw26+XuywRkSOiwCgBM+MfrzmdcVVxbnzhVLKnvBeW/QOs/UW5SxMROWwKjBKZWJPkf/3FfF5t6eBLkU/iU+bD/R+FrS+UuzQRkcOiwCihi+dN4pMXn8A9y7dz/1u+BlX18L1rYceGcpcmInLIFBgldtNl83jnyZO4eel2nv+Tb0G2F/79Kti1qdyliYgcEgVGiUUjxm3XncWJjTV8+GftvPZnS6C3He5+N+x+vdzliYgUTYExCmqSMf7txiaS8SjXPNTO5nd/H3r2wOIrYNvL5S5PRKQoCoxRMmNCmiV/dR4A1/yog+ar7wPPBaGxYVmZqxMROTgFxig6cVItS/7qPLI5530PtrHu3Q/BuGnw3ffBb/63bu4TkTFNgTHK5k6uZcmi84mY8b7vvc4zf/p9OPW98PiXYcl10DkmHywoIqLAKIe3TK7lwU+9jan1KW645xW+M+W/4gu/Bht+Bf96EfzxN+UuUURkP6V+pvcVZrbWzNab2c0F5v+zma0Mh1fNbHfevGzevIdLWWc5HFdfxX1//TYunjeJL/10NX+7vomOG34G0Rj8+5/Bz26CnvZylykiMqBkz/Q2syjwKnAp0Az8Hrje3V8ZZvn/AJzl7h8Np9vdveZQvrNcz/Q+Eu7OXU9u5J9+uZaZE9L8y1+cxKmv3AbP3gl1M+Dyf4CT3w1m5S5VRI5BY+WZ3ucC6919o7v3At8Hrj7A8tcDS0pYz5hkZnziohNY8lfn09mb4ap/fZ6vRT5M74d+Cska+OEH4f+9R5ffikjZlTIwpgGb86abw7b9mNnxwGzgV3nNKTNbbmbPmNl7Slfm2HDu7Ak8+tmLeN9Z0/jGsg1c+aMsv/3Th2DhP8Ebf4A73gY/vBFaVpe7VBGpUKUMjELHUIY7/nUdcL+7Z/PaZoa7Se8HbjOzEwp+idmiMFiWt7a2HlnFZVaXjvO1vziDuz/yVnoyWW74zgo+vuYcXvvgU/AnN8H6x+CbFwSdGL65qtzlikiFKWVgNAMz8qanA1uGWfY6hhyOcvct4etG4AngrEJvdPe73L3J3ZsaGxuPtOYx4eJ5k1j6Hy/i81ecxNMbtvPOO17klvb38uZHn4W3fzboJv3OC+Hbl8OLP4RMT7lLFpEKUMqT3jGCk95/CrxBcNL7/e7+8pDl5gG/BGZ7WIyZjQc63b3HzBqAp4Grhzth3u9oPOl9MC1t3fzz0nXct3wzETOuO3cGnzpvPFM2PgjLF8POjZBugLNugKaPwPhZ5S5ZRI4ih3LSu2SBERZyJXAbEAUWu/s/mNmXgeXu/nC4zH8DUu5+c9773gb8K5Aj2Au6zd2/fbDvOxYDo9/mnZ1884n13Le8mYgZf/nWGXzyotkct+OZIDjWPgLucMICmH8dnPSu4KS5iMgBjJnAGG3HcmD0C4JjA/ev2EzO4fJTJ3PD+cdzwcQu7Pl74IXvw57XIZ4OQuPU98IJl0C8qtyli8gYpMCoAM27Ornn6df44fLN7O7s48RJNXzgvJlcNX8KE3euhBd/AK/8CLp2QaIG5l4K866EE98J6QnlLl9ExggFRgXp7svy0xe38v+e3sQLzXuIRox3zG3gPWdN49J5E0hveRpe+TGseQQ6WsCiMPN8eMvlMGcBTD4NIuohRqRSKTAq1Jo39/KjP2zhxyvfYOueblLxCBee0MCCkyaxYF4D0zrWBOc61v4cWsJrD6obYc7FwWGrOQtg3NRy/gkiMsoUGBUul3Oe27STR17ayq/WtNC8qwuAeZNrufikRhbMm8Q547uIb/o1bFwGG5+AjvAelsaTguCY9fZgT6S6oXx/iIiUnAJDBrg7G1o7WLamhWVrW3jujzvJ5JzaVIx3zG1kwUmTuGjuRBo71wcPctrwK3j9ach0Bx8w8cQgOGZeEAwT5qhfK5FjiAJDhtXW3cfv1m9n2ZpWlq1toaUtuOlv/vQ6FsybxIKTJjF/corImy8EwfH6M8Frd9iRcKoOpp4RDmcGrxNO0HkQkcPhDn1d0NcJvR3B0D++X1s79HYWXjZRDe//wWGVoMCQorg7L2/ZyxNrW1i2tpU/vL6LnMO4VIxzjh9P06wJvHXWBOZPqyW1e0MQHFtfgK0rg84Qs73BByVqYMrpeUFyBjTMC7pqFzlWuAcb5562YOPdszd4BEFP2762QRv14caHbOyH7TGpgEgsCId4NSTS+8ZrJ8Nf3H1Yf5YCQw7Lro5enlzXyjMbd7J8007WtQTP40hEI5w+vY6zZ9Zz+vR65k+r4/j6OLZ9bRgg4fDmS8E/AoBYKjgf0jgPGuYGAdLwluCQVixRxr9SKkr/L/je/g173kZ+YKPfVqBtSBD0v3oxj1G2cEMebtAHxtPBj6sC47lYmt5oFT2WoosUnaTo8ASdnqQtl6Qtl2BvNkF7JkJHTyYYerN09mbo6MlSFY9y5wfPOaxVpMCQEbGro5flr+1i+aadPLdpJy9v2UtvJvgHMy4V4/TpdZw+rZ750+s4Zeo4ZtYniezasC9Atq2C7etg7xv7PtSiMGF2EB4Nc4PX8bOhfibUTtVeiey/ke9tDzbgh7TRz1tuUJ+mw7Bo0DNCclywIU/WhkM4nqgdaMslauiJpum2NJ2RNJ2kaSdFu6doyyZpy8bo6M3tt1Hv7A2nezK092TozJvX1VdEjaF41KhOxqhOxEgnoqSTMSbVJvnWh4ra5u//pyswpBR6Mzle3dbGqjf28OIbe3ipeQ9r3txLXzb4f6gqHmXu5BrmTa5l3pRwmFxLY6IX27E+CI/tr+4bdmyAXN++L7Ao1E2D+uODAOkf6mZA7RSomRT8g9ZJ97HFPdiz7N+oD93A93YMbhtuPL+t2F+8RfS9AAALlklEQVTy/Rv2gY18/+u4QW2eqKUvVkOnpei0NB1U0ZZLsTtXxZ5skj19Udp7s7T1ZGjvDjbo7d2ZQdNBAGTo7iumtkAiFqE6ESWdiFGdHPxakww2+NX9r4kY6WTwGgRCEAaDXhMxErGRPV+owJBR05PJsvbNNtZsbWPttjbWvhm8trbt60G3rirOjAlVzJyQZsb4NDMmhENdnOnWQmLvZtizGXa/Hg7heNtW9ju+G01CzeQgPAaGycH9JFXjoao+fB0PqfrgJH0kOrorZazLZvKOt3cM3lAP2rh3DNmoD53Oey32OHwsFWzIE9V5G/pwOhFu8AfaagZ+zXdRFWzkww39nlyK3Zk47T25YOPek6Gtu4+2vI190JYZmJ/NHbzGZCxCbSrYmNf0vybj1CSj1KT6f9UP3vgP3tCH7WFbPDr2LwZRYEjZ7Wjv4dVt7ax9cy/rW9vZvLOLzbs6ad7ZRW923y80M5gyLpUXJGGwTEgzY1yMSbntRPZuhvYWaN8WDi3BXev9bR3bGX6DZUFoJMcN3hgNHGaogWgCYskgjGKJYDq/LRoPxxP7XoeOR+PBL23YV0vB6SHzPAe5DGT7gr2tbCa4mCDXF7aF09lwOtMFfd3Ba6YnOHST6d73muneN7+ve3BbXxgO2UPoDj+eztvAh+tv0HRt3ga/ZiAE+qJpuqyKDlK0eRVtuSR7c3H29trABr2jZ/Av+EEb+HB+e2+GYjZRVfFgg14bbuhr8zb2QwNg37xw2WR8YN5I/3o/GigwZMzK5ZyWth427+rk9R2dwevOIEhe39nJtrbuQRuIRCzC5HFJJlYnaahJ0liboKEmmTckaKiO0hjpoDrXRrRnT9B/Vteu4FLg/vH+E5j7HRMPN6CZHg7papWxIpqAWBXEU0GADYxXBdPxquBXfbxq30nYZO2+DXyiGk/U0BOposuqwl/ySdpJ0ZZJ0Jlh0LH3gePuvdmBk6/5v+L7D+P0n+s6EDOCDXj+r/lU+Gu+/5d9XghUhxv7/OVrk3Gqk1FiR8Ev+bHqUAJDZxhlVEUixpS6FFPqUrx11v6dIPZksryxKwiPzbu62Lyzk5a93Wxv76V5VycrN+9iZ0cvwx1dCH5BxqlNHce41PHBBiYVZ1xVjNr64NfmuNS+X53JWJRkPEIyaiQiOVKWJRnJkCBDkj4SliHufVi2N/iln+kJf+2HIdPfnu0Lz62E51cGzrMUmB46HokHJ/sjcXKROBmiZC1Ghih9FiPjMbLE6CVKD0m6iNOVi9GVi9Odha6+LN19Obr7snRngvGevmww3ZcL2tqzdPZm6erNDjoR29WbobMvi3s70F7Uf8N03jH5/l/qU+tSA7/Yq/tDYCAAYgUO8wTH7U3no44qCgwZU5KxKHMaa5jTOPyzPLI5Z2dHL9vbewaGHe29tHUHv3b3dvcNHM9ube9h4/aOoL2rj0wRx7GHMgsuLU7EIkQjRtQMsyTRSDIct6A9YvuywCHnTs7BcXK54L4XJ6/dIZPLkck6vdkcmWyOnGeADHBkT1FMxiKk4lFS8QhV8SipeJSqRJR0Isr46sTAidR0fPCJ1XQ8Ouj4fFV88PH6VCxKJKKNfKVSYMhRJxoxGmuTNNYmD+l97k53X4627j72dvfR3pOlN5OjJ9P/OmS8L0dvNvi13hO2BRt7J5sLDq/l3Mm6h+OQDY+nRcyIWLB/EQlDxQwi1j8NYMSjRjwaIRY14pHIwHiivy0aIR41YpFgOhkLQiAVhkAqHiEVyxuPR0nGIvrlLiWhwJCKYWZUJYJf2pPGpcpdjshRR2eKRESkKCUNDDO7wszWmtl6M7u5wPwPm1mrma0Mh4/nzbvRzNaFw42lrFNERA6uZIekzCwKfAO4FGgGfm9mD7v7K0MW/YG7/+2Q904AbgWaCK51XBG+d1ep6hURkQMr5R7GucB6d9/o7r3A94Gri3zv5cBSd98ZhsRS4IoS1SkiIkUoZWBMAzbnTTeHbUNdY2Yvmtn9ZjbjEN8rIiKjpJSBUei6vqEXwf8EmOXu84HHgH8/hPcGC5otMrPlZra8tbX1sIsVEZEDK2VgNAMz8qanA1vyF3D3He7ef4fSt4Bzin1v3mfc5e5N7t7U2Ng4IoWLiMj+ShkYvwfmmtlsM0sA1wEP5y9gZlPzJq8CVofjvwQuM7PxZjYeuCxsExGRMinZVVLunjGzvyXY0EeBxe7+spl9GVju7g8Dnzazqwj6QtgJfDh8704z+wpB6AB82d13Huw7V6xYsd3MXjvMkhuA7Yf53mOR1sc+WheDaX0MdrSvj+OLXfCY6q32SJjZ8mJ7bKwEWh/7aF0MpvUxWCWtD93pLSIiRVFgiIhIURQY+9xV7gLGGK2PfbQuBtP6GKxi1ofOYYiISFG0hyEiIkWp+MA4WI+6xyIzW2xmLWa2Kq9tgpktDXsHXhre/4IFbg/Xz4tmdnb5Ki8NM5thZsvMbLWZvWxmnwnbK3KdmFnKzJ4zsxfC9fGlsH22mT0bro8fhPdXYWbJcHp9OH9WOesvBTOLmtkfzOyn4XRFrouKDoy8HnUXAqcA15vZKeWtalTczf6dOd4MPO7uc4HHw2kI1s3ccFgE3DFKNY6mDPB37n4ycD7wN+H/B5W6TnqAS9z9DOBM4AozOx/4R+Cfw/WxC/hYuPzHgF3ufiLwz+Fyx5rPsO/GYqjQdVHRgcGR9ah71HL3JwlulMx3Nfv68vp34D157fd44Bmgfsgd+kc9d9/q7s+H420EG4ZpVOg6Cf+u9nAyHg4OXALcH7YPXR/96+l+4E/tGHpGrJlNB94F/Fs4bVTouqj0wFCvuPtMdvetEGxAgUlhe0Wto/AQwlnAs1TwOgkPwawEWggeL7AB2O3umXCR/L95YH2E8/cAE0e34pK6DfjPQC6cnkiFrotKD4yie8WtYBWzjsysBngA+Ky77z3QogXajql14u5Zdz+ToOPPc4GTCy0Wvh6z68PM/gxocfcV+c0FFj3m1wUoMIruFbcCbOs/rBK+toTtFbGOzCxOEBbfc/cHw+aKXicA7r4beILg3E69mfX3P5f/Nw+sj3B+Hfsf8jxaXQhcZWabCA5ZX0Kwx1GJ66LiA+OgPepWkIeB/men3wj8OK/9Q+GVQecDe/oP0xwrwmPM3wZWu/vX82ZV5Doxs0Yzqw/Hq4B3EpzXWQb8ebjY0PXRv57+HPiVHyM3eLn7F9x9urvPItg+/MrdP0AFrgsA3L2iB+BK4FWCY7T/pdz1jNLfvATYCvQR/CL6GMFx1seBdeHrhHBZI7iSbAPwEtBU7vpLsD7eTnDY4EVgZThcWanrBJgP/CFcH6uAW8L2OcBzwHrgPiAZtqfC6fXh/Dnl/htKtF4uBn5ayetCd3qLiEhRKv2QlIiIFEmBISIiRVFgiIhIURQYIiJSFAWGiIgURYEhcgjMLGtmK/OGEevh2Mxm5fcgLDLWxA6+iIjk6fKgywyRiqM9DJERYGabzOwfw+dIPGdmJ4btx5vZ4+FzMx43s5lh+2Qzeyh85sQLZva28KOiZvat8DkUj4Z3WouMCQoMkUNTNeSQ1F/mzdvr7ucC/0LQ3xDh+D3uPh/4HnB72H478GsPnjlxNvBy2D4X+Ia7nwrsBq4p8d8jUjTd6S1yCMys3d1rCrRvInjo0MawI8M33X2imW0Hprp7X9i+1d0bzKwVmO7uPXmfMQtY6sFDeTCzzwNxd//vpf/LRA5OexgiI8eHGR9umUJ68saz6DyjjCEKDJGR85d5r0+H408R9HIK8AHgt+H448AnYeBhReNGq0iRw6VfLyKHpip8El2/X7h7/6W1STN7luCH2PVh26eBxWb2OaAV+EjY/hngLjP7GMGexCcJehAWGbN0DkNkBITnMJrcfXu5axEpFR2SEhGRomgPQ0REiqI9DBERKYoCQ0REiqLAEBGRoigwRESkKAoMEREpigJDRESK8v8Br8/CCT5j4/gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# how many epochs to train\n",
    "n_epochs = 1 # This will just train for one epoch. You will want a higher number\n",
    "\n",
    "# Loop over the epochs\n",
    "for ep in range(n_epochs):\n",
    "        \n",
    "    # Each epoch is a complete over the training data\n",
    "    for i in range(X_train.shape[0]):\n",
    "        \n",
    "        # pick one example\n",
    "        x = X_train[i]\n",
    "        y = y_train[i]\n",
    "\n",
    "        # use it to update the weights\n",
    "        update_weights(x,y,W,b,Wp,bp)\n",
    "    \n",
    "    # Calculate predictions for the full training and testing sample\n",
    "    y_pred_train = [dnn(x,W,b,Wp,bp)[0] for x in X_train]\n",
    "    y_pred = [dnn(x,W,b,Wp,bp)[0] for x in X_test]\n",
    "\n",
    "    # Calculate aver loss / example over the epoch\n",
    "    train_loss = sum((y_pred_train-y_train)**2) / y_train.shape[0]\n",
    "    test_loss = sum((y_pred-y_test)**2) / y_test.shape[0] \n",
    "    \n",
    "    # print some information\n",
    "    print(\"Epoch:\",ep, \"Train Loss:\", train_loss, \"Test Loss:\", test_loss)\n",
    "    \n",
    "    # and store the losses for later use\n",
    "    train_losses.append(train_loss)\n",
    "    test_losses.append(test_loss)\n",
    "    \n",
    "    \n",
    "# After the training:\n",
    "    \n",
    "# Prepare scatter plot\n",
    "\n",
    "\n",
    "y_pred = [dnn(x,W,b,Wp,bp)[0] for x in X_test]\n",
    "\n",
    "print(\"Best loss:\", min(test_losses), \"Final loss:\", test_losses[-1])\n",
    "\n",
    "print(\"Correlation coefficient:\", np.corrcoef(y_pred,y_test)[0,1])\n",
    "plt.scatter(y_pred_train,y_train)\n",
    "plt.xlabel(\"Predicted\")\n",
    "plt.ylabel(\"True\")\n",
    "plt.show()\n",
    "\n",
    "# Prepare and loss over time\n",
    "plt.plot(train_losses,label=\"train\")\n",
    "plt.plot(test_losses,label=\"test\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Epoch\")\n",
    "plt.ylabel(\"Loss\")\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}