StanfordMLPython/ex2/ex2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Programming Exercise 2: Logistic Regression</h1>\n",
    "\n",
    "<h3>Introduction</h3>\n",
    "\n",
    "In this exercise we will implement logistic regression and apply it to two different datases.\n",
    "\n",
    "<h3>Files Included in this exercise</h3>\n",
    "\n",
    "- ex2data1.txt\n",
    "- ex2data2.txt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>1 Logistic Regression</h3>\n",
    "Here we will build a logistic regression model to predict whether a student gets admitted into a university given the results of two exams. Our training set consists of samples of applicants' scores on two exams and an admissions decision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# used for manipulating directory paths\n",
    "import os\n",
    "\n",
    "# Scientific and vector computation for python\n",
    "import numpy as np\n",
    "\n",
    "# Plotting library\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Optimization module in scipy\n",
    "from scipy import optimize\n",
    "\n",
    "# tells matplotlib to embed plots within the notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Visualizing the data</h4>\n",
    "\n",
    "Before we begin on the algorithm we load and visualize the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "# The first two columns contains the exam scores and the third column\n",
    "# contains the label.\n",
    "data = np.loadtxt(os.path.join('Data', 'ex2data1.txt'), delimiter=',')\n",
    "X, y = data[:, 0:2], data[:, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotData(X,y):\n",
    "    # New figure\n",
    "    fig = plt.figure()\n",
    "\n",
    "    # Find indeces of positive and negative examples \n",
    "    # Then plot them seperately (Don't try to plot then label after)\n",
    "    pos = y == 1\n",
    "    neg = y == 0\n",
    "\n",
    "    plt.plot(X[pos,0],X[pos,1],'k*', lw=2, ms=7)\n",
    "    plt.plot(X[neg,0],X[neg,1],'yo',mec='k',ms=7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3hU5bX48e8iUZKANYKWg0VLQESPoihgsSImolhaoS31AqgEpVUopHJsvURttdCfKPo8WqEinoOCnIBWlFOvPajgwXgpF0XEK5AEi6QaUVRkuLp+f8yeEMJkkpnMvs2sz/PMk8yeZPbKJJm139t6RVUxxhhjANr4HYAxxpjgsKRgjDGmniUFY4wx9SwpGGOMqWdJwRhjTL1cvwNojcMPP1y7du3qdxjGGBMqq1at+kxVj4j3WKiTQteuXVm5cqXfYRhjTKiIyMamHrPuI2OMMfUsKRhjjKnnWlIQkQdF5FMRWdvgWAcReV5E1jkfD3OOi4jcKyLrRWSNiJzqVlzGGGOa5mZLYQ7wo0bHbgBeVNUewIvOfYAhQA/ndiUw08W4jDHGNMG1gWZVXSYiXRsd/ilQ7Hw+F3gJuN45/rBGCzG9LiKFItJZVWvdiq81IpEIixYtYsOGDXTv3p3hw4eTl5fnd1jGGNNqXo8pdIq90Tsfv+sc/x7wzwZft8k5dgARuVJEVorIyrq6OleDjWfFihUUFR3J9OlX8eGHtzB9+lV07dqZFStWeB6LMcakW1CmpEqcY3HLt6rqA8ADAH379vW0xGskEmHo0MFMnLiVAQNiR7dRWQlDhw6mpqbWWgzGmFDzuqXwiYh0BnA+fuoc3wQc1eDrugCbPY6tWYsWLaKoaE+DhBA1YAAUFe3hiSee8CcwY4xJE6+TwpNAqfN5KfC3BsdHO7OQ+gNfBnE8YcOGDRxzzDdxHzvmmG+oqqryOCJjjEkvN6ekLgBeA3qKyCYRGQvcDpwrIuuAc537AM8CVcB64D+BX7sVV2t0796d9evbxX1s/fp2dOvWzdXz19bWctZZZ/Gvf/3L1fOYYLHfu/GSa0lBVUeqamdVPUhVu6jqbFXdoqqDVLWH8/Fz52tVVSeoandV7aWqgaxdMXz4cKqrc6ms3P94ZSVUV+cyfPhwV88/ZcoUKisrmTx5sqvnMcFiv3fjJQnzdpx9+/ZVr2sfrVixgqFDB1NUtIdjjvmG9evbUV2dy1NPLaZfv36unDM/P58dO3YccDwvL49IJOLKOc0+tbW1jBgxgkcffZR/+7d/8+y89ns3bhGRVaraN95jVuYiSf369aO6ejNlZbPo2XMyZWWzqKmpdS0hAFRVVVFaWkp+fj4QfbMYM2YM1dXVrp0zmzXurvHrSj0sv3fr3sowqhraW58+fTRbPPzww5qbm6sFBQWak5Oj8+bN8zukjDV+/Hht06aN5uTkKNGp0fvd8vLyPIsl1d/75s2bdeDAgVpbW+tyhPter/Hjx7t+LpMewEpt4n3V9zf21tzcTgrbt2/XiooKnTx5slZUVGgkEnH1fIkUFxdrbm6uTps2TXNzc7WkpMS3WDJVXl5e3CQQu+Xn5+uYMWM8eaONSfX37sUbdVOvl5dJ06QmUVKw7qMmBG3l8tFHH83zzz/Ptddey+LFizn66KN9iSOTxeuuOfPMM8nJyaGgoIBdu3YxaNAgT8cVkv295+fnIyLMnDmTb7/9lpkzZyIi9T9TOoWle8skqalsEYabWy2F7du3a6dOhTplCrp06b7blClop06FvrYYjLsad9ccf/zxoWqhbd68WUtLSzU/P9+T1o11a4YT1lJIjq1czl4PPvggALfeeisiwueffx6qFlrnzp0ZNGgQu3fv9qR10/j1it034WVJIQ5buZy9GnfX/OhHP6K4uBiAkpIS5syZ42t8LeHlG7V1a2aeoBTEC5Tu3bvz7LPtgG0HPLZ+fTuGDHF35bLxz9y5c+s/LykpoaSkxMdoUhN7oy4uLqZv3777/Uzplgmvl9mfLV6LY8eOHXTt2rlRNdToyuUZMwqtGqoxJtQSLV6zlkIceXl5PPXUYoYOHcwzzxy4ctkSgjEmU9mYQhP8WLlsTEvZKmLjFmspJJCfn8+oUaP8DsOYAzQsvXHffff5HY7JIDamYEyIWJE8kw5WEM+YDBG2VcTWzRU+lhSMCRGvF6e1lu0FET6WFIwJmTCsIvayBpNJL0sKxoRMGFYRN+7mEhEuuuiiwHZzmX1s9pExIROGVcSxbq6KigpycnLYu3cvn332WWC7ucw+1lIwxgDpHxQeM2YMe/bsYe/evQAsWbLEupBCwJekICJXi8haEXlHRCY5xzqIyPMiss75eJgfsRmTrdI9KDx8+HDOO++8+iSQk5MT6JlSJsrzpCAiJwK/Ak4DTgbOF5EewA3Ai6raA3jRuW+McVFtbS1t2rRxZVD4scce45JLLqmfKQUEeqaUifKjpXA88LqqblfVPcD/AT8HfgrEOkvnAj/zITYTEDa/3RtTpkwBoGfPnq6sfQjDTCmzPz+SwlpgoIh0FJEC4MfAUUAnVa0FcD5+N943i8iVIrJSRFbW1dV5FrTxls1vd1fDKaOqygcffFC/Ijqdax/CMFPKNNLUlmxu3oCxwBvAMuB+4G5ga6Ov+aK553FrO06zz+bNm3XgwIGebVZvm8F7o/G2nW3atFER0d///veh2HbUa17/H7iNoG3HqaqzVfVUVR0IfA6sAz4Rkc4AzsdP/YjN7M/rK/awlXEIq8Yro1WV8vJyJk+ebFf0cWRVy7WpbOHmDfiu8/Fo4H3gMOBO4Abn+A3AtOaex1oK7vHzit02g/dGcXGx5ubm6rRp06x10IRMbbkStJYC8LiIvAs8BUxQ1S+A24FzRWQdcK5z3/jEzyt2G5z0RqL+/mwe6G/4s2dly7WpbBGGm7UU9pfufk+/rthHjx6tS5cuVVXVJUuWaGlpqSfnNfuMHz9e27Rpo+PHj/c7FM81/tkzseVKgpaC72/srblZUthfuv+RrXsh+2Rqd0lLNPWzt2nTxtP/Ay8GtRMlBStzkQHcqkhp0wmzT1Z2lzia+tl/8YtfePp/4PugdlPZIgw3aylENZ5emJ+fr2PGjMmY6XNuy7Tphq2Vid0lLeXnz+5lKw1rKWS2sG28EjS+X5kFRGyAdebMmUB2DvT7OckhKK00SwoZwmbsJC+oG8H4NfMnlhy3bNmStd2GfnaZBubirqkmRBhu1n20j83YSV5Qu928nvmTzYPLQePV5A5s9pEx8QWp/9zrN+fYWMrq1asDmRyzkVcXd4mSgnUfmawWpG43r/uUY91Fs2bNCka3hWHu3LkUFxcD0V315syZ43kMlhRMVgvStFuv+pTjjaWMHj2aPXv2BCI5Gn9ZUjBZLQhXZg150XKJ1yLp3r07jz/+eCCSo/GXJQVjAsSLlku8Fsmtt97K8OHDgWAkR+MfiY45hFPfvn115cqVfodhTOiUlJRQWVnJbbfdxo033siZZ57JkiVL/A7LeEREVqlq33iP5XodjDHGf7EWSXFxMX379mXu3LnNf5PJCtZSMMaYLJOopWBjCsYYY+pZUjDGGFPPkoIxxph6lhSMSYNs3r7SZBZLCsakgZXfNpnCZh8Z0wr5+fns2LHjgON5eXlEIhEfIjKmeYGbfSQi/yEi74jIWhFZICJ5IlIkIv8QkXUi8qiIHOxHbMYkIygboxiTLp4nBRH5HvAboK+qngjkACOAO4C7VbUH8AUw1uvYjElWYDZGMSZN/BpTyAXyRSQXKABqgbOBhc7jc4Gf+RRb6EQiEebPn8+UKVOYP39+3O4M454gld82prU8Twqq+jFwF/AR0WTwJbAK2Kqqe5wv2wR8L973i8iVIrJSRFbW1dV5EXKgrVixgqKiI5k+/So+/PAWpk+/iq5dO7NixQq/Q8saQSq/bUxreT7QLCKHAY8DFwNbgcec+7eo6jHO1xwFPKuqvRI9V7YPNEciEYqKjmTixK0MGLDveGUlzJhRSE1NLXl5ef4FaDxVW1vLiBEjePTRR637yiQUtIHmc4BqVa1T1d3AE8APgUKnOwmgC7DZh9iS5mfXzaJFiygq2rNfQgAYMACKivbwxBNPeBaL8Z9NizXp4EdS+AjoLyIFIiLAIOBdYClwgfM1pcDffIgtKX533WzYsIFjjvkm7mPHHPMNVVVVnsRh/F28Fm8nNRGpnxFlTDL8GFP4B9EB5TeAt50YHgCuB64RkfVAR2C217ElIxKJMHToYCZO3MrUqdsYO1aZOnUbEyduZejQwZ60GLp378769e3iPrZ+fTu6devmegwmys+rdJsWa9LJl9lHqnqLqh6nqieq6mWqulNVq1T1NFU9RlUvVNWdfsTWUkHouhk+fDjV1blUVu5/vLISqqtz63fSMu4JwlW6TYsNlrCXPLEyFykKQtdNXl4eTz21mBkzCikvb8/s2UJ5eXtmzCjkqacW2yCzB4JylW7TYoMj7GM7lhRS5HbXTUsHsPv160d19WbKymbRs+dkyspmUVNTS79+/Vp1ftMyQblKD/O02LBfWccEodWYFqoa2lufPn3UL5FIRDt1KtQpU9ClS/fdpkxBO3Uq1EgkkvJzL1++XDt1KtT+/dvrpZeK9u/fXjt1KtTly5en8Scw6VJcXKy5ubk6bdo0zc3N1ZKSEr9DCpXx48drmzZtdPz48X6H0iqbN2/W0tJSzc/PV0Dz8/N1zJgxWltb63doBwBWahPvq76/sbfm5mdSUHXnzXv79u2uJRsvbd68WQcOHBiofwi3Yho9erQuXbpUVVWXLFmipaWlaX3+TJWXl6fAAbe8vDy/Q0vZww8/rLm5uVpQUKA5OTk6b948v0OKK1FSsO6jFEUiEdatW8eVV5Zx6qmj6dbt92npugnCAHY6BLFf1a2Y5s6dS3FxMQAlJSXMmTMnrc8fFsl2AwVlPCadMmFsx5JCChquT6iu/hNvvPEws2bdS48ePVo9uBuEAezWCGK/ahBjykTJJt2gjMekU5jHdmIsKSTJ7fUJYV97EMSrvyDGlElak3Qz4cq6oUxoNTabFETkWBF5UUTWOvdPEpGb3Q8tmNzu3gn72oMgXv0FMaZM0pqkmwlX1pmmJS2F/wTKgd0AqrqG6P4HWcnt7p1MWHsQxKu/IMaUKZpKuqra7BhDJlxZZ5qWJIUCVV3e6NieuF+ZBbzo3gn72oMgXv0FMaZMEi/pBnGygWles6WzReQ5YCLwmKqeKiIXAGNVdYgXASbiR+nsHTt20LVr52bLVUciERYtWsSGDRvo3r07w4cPD8VVvjGpKC0t5fLLL6e4uJi2bduya9euA77G9q0OjkSls5tdCwB0A14AtgMfA5XA95v7Pi9ufq1TaG59gi0+C58grqsIqzAt4spWJFinkBs3U+zLJm2I7qV8joi0A9qo6tfpy1fhFOveWbRoEVVVVQwZ0q2+JdBwdtK+lsQ2Kith6NDBtvFNQDXs6rjvvvtsw5pWiI0xVFRUUFBQwM6dO21gP0yayhaxG7Csua/x6+b3iuZ4KioqtH//9vutRo7d+vdvrxUVFX6HaBpoalVtTk5ORpRe8IuV/nBXa1u2tHJF8/Mi8jsROUpEOsRurmWpkAv74rNs03g6ZczevXttoVsr2MC+u9wcxG9JUrgCmAAsA1Y5t+zdGLkZbsxO8nPLz0zXeDplTk4OAwYMsIVurWRTTd3hxer8ZpOCqhbFuQV7Wa2P0r34zO8tP7NB4+mUW7ZssYVuJpC8WJ2fcKAZQEQOAsYDA51DLwGzVHV32qLIILHFZ0OHDuaZZ/ZwzDHfsH59O6qrc5NefGaD1t6IdXUUFxfTt29fRo4cCUSTxI033siDDz7IpZde6nOUxngziN+S7qOZQB/gPufWxzlmmpCuxWeZUjE16Bp3dZx33nnWH24Cy+3V+c22FIB+qnpyg/tLROSttEaRgfLz8xk1alSrnsMGrf0xd+7c+s9LSkooKSnxMZrgsGm6wdC4Zdvw7zUdWtJS2Csi3WN3RKQbsDfVE4pITxFZ3eD2lYhMcmY1PS8i65yPh6V6jkwR9oqpmSRTtoyMSeXnyaSyFWH+fbo+iN/UXNXYDRgEfER0LOH/gBqgpLnva8kNyAH+BXwfmAbc4By/Abijue8P4jqFdHJzy08/hXH1cKZsGRmTzM+TiTukZdrvM1m0djtOoC1wEnAy0LYl39PC5x0MvOJ8/gHQ2fm8M/BBc9+f6UlBNTNLZoTpHzLT3hBT+XkyqWxFpv0+U5UoKbRkP4UJQL6qrlHVt4ACEfl1Eo2RREYAC5zPO6lqLYDz8btNxHOliKwUkZV1dXVpCiO4wl4xtaEw7oCWaRv0pPLzxCuN3adPHy6++OLQdb8E/feZqFvLsy6vprKF7ruaXx3n2JvNfV8Lnvdg4DOiyQBga6PHv2juObKhpZBJwnrFGZbN2FsqlZ+ncdmKI488MjStvcaC/PtM1IpOZwub1nQfAWtwSmzrvnGAd5r7vhY870+BxQ3uW/dRFgjyP2RTMq2OTyo/z+jRo3Xp0qUZ0f0SxN9notfVjdc8UVJoyeyj/wX+KiKDRORsot09f0+6SXKgkezrOgJ4Eih1Pi8F/paGc5iACeMOaJlWxyeVnyc24yXo3S8tEcTfZ6LX1fPXvKlsofuu4NsA44CFwOPAVUBOc9/XzHMWAFuAQxsc6wi8CKxzPnZo7nmspRA+sStOVdUlS5ZoaWmpr/GY5IWxtRcGiV7XdL/mtHb2ke574+4AnJTM97h5s6RgjPeC2P2SCRK9rul+zRMlhZbUPnoJGEZ09fNqoE5E/k9Vr0lbc8V4xrYJNa3l9orabJXodfXyNW/JHs1vquopIvJL4ChVvUVE1qjqSa5F1UJ+7NEcZitWrGDo0MEUFR1YqC+MU1wzgZWOMH5ItEdzS2of5YpIZ+Ai4Ka0RmY8YxVXg6nxNqDG+K0ls48mE52BtF5VVzi1j9a5G5ZJN6u4GixhXMiXqcJcB8kNLdlk5zFVPUlVf+3cr1LVX7gfWnZxe3c1q7gaLJkwtTNTuFHoL8yJpiUtBeMyL3ZXs4qrwRKvdITt8OYtN1trYa4oa0nBZw37+qdO3cbYscrUqduYOHErQ4cOTluLId3bhPolzFdgjYVxIV8mcaO1lgndgpYUfOZVX39sm9AZMwopL2/P7NlCeXl7ZswoTHqbUD+F+QqssSCurM0mbrTWMqFbMGFSEJHjnPIW7Rsd/5G7YWUPL/v6w1xxNROuwBpzfbMU06x0t9YyoVuwyaQgIr8hWn+oDFgrIj9t8PBtbgeWLbzu649tE3rzzTczatSo0LQQMuEKzASPG621sHcLJmop/Aroo6o/A4qB34vI1c5j4nZg2SJT+vpbKtVZVplwBWaCx43WWti7BZtc0Swi76rqvze4355oUbx3gbNVtbc3ITYtU1Y0Z8tK49b+nCUlJVRWVnLbbbdx4403cuaZZ7JkyRIPIjcms6S6ovlfItJbVVcDqOo2ETkfeBDo5UKcWSvW179o0SKqqqoYMqRbxtUkSseKaqu5Y4z7ErUUugB7VPWAuX8icoaqvuJ2cM3JlJZCNpg/fz7Tp1/F1KnbDnisvLw9ZWWzGDVqlA+RGZN9UmopqOqmBI/5nhBMuNiKamPCwdYpGE/YimpjwsGSgvFEts2yMiasWpwUROQ7ItIhdnMzKJN5Gq+onjULfvWrttx+ex4TJmTXfk2ZVKrDeMerv5tmk4KIXCUinwBrgFXOzUZ3HW5XN80ksVlWP/7xdTz3XB6FhcLPf76TZ5+dlvYCgEGWSaU6jHe8+rtpyc5r64DTVfUzVyNJgd+zj7JlfUFzktniMxKJUFR0ZKOpqdFupBkzCjN6s5/8/Py4Fw15eXlEIhEfIjJh4MbfTaLZRy3pPtoAbE/pzE0HVCgiC0XkfRF5T0ROd7qlnheRdc7Hw9J5znTzqrpp0CVb9jubN/uxUh0mFV7/3bQkKZQDr4rILBG5N3Zr5Xn/DPxdVY8DTgbeA24AXlTVHsCLzv3AyuY3t5hUEqPbU1OD3F9vpTpMKrz+u2lJUpgFLAFeZ9+YwqpUTygi3wEGArMBVHWXqm4FfgrElqjOBX6W6jm8YPPuU0uMbk9NDXp/fdiLpRl/ePl305KksEdVr1HVh1R1buzWinN2A+qAh0TkTRH5LxFpB3RS1VoA5+N3432ziFwpIitFZGVdXV0rwmgdm3efWmJ0a2pqWEprh71YmvGHl383LRlo/n/ARuApYGfsuKp+ntIJRfoSbXWcoar/EJE/A18BZapa2ODrvlDVhOMKfg4079ixg65dO2flgGlMqqUr3Bigr62tpby8nL/+9a9EIhHy8/O5+OKLmTp1qnXPGNNIqgXxYmL/1eUNjinRK/5UbAI2qeo/nPsLiY4ffCIinVW1VkQ6A5+m+PyeiM27Hzp0MM88c+CbW6YnBIhe9V9zzQQqKzkgMSa66nejAGCs37WiooKCggJ27txp/fXGpKDZloIrJxV5Gfilqn4gIrcCsX6YLap6u4jcAHRQ1esSPY/fU1Jh33TMqqoqunXLvOqmzQnStFwrrW1My7S2pYCInAj8O1D/bqeqD7cipjKgQkQOBqqAy4mOb/xVRMYCHwEXtuL5PRPbySxbBanst5XWNqb1WjKmcAvRndf+HXgWGAJUquoFrkfXjCC0FIwx/qitrWXEiBE8+uij1k2YpNYuXrsAGAT8S1UvJ7quoG0a4zM+sPIcJuyCPv04rFqSFCKq+i2wx1lj8CmpDzKbAEh2FbIxfom3GDEs04/DqiVJYaWIFAL/SXTR2hvAclejMq6x8hwmTOK1BqxciLuaTQqq+mtV3aqq9wPnAqVON5IJISvPYcIgUWvAyoW4qyWls8fGPlfVGuAdZ/DZhJCV5zBh0FxrwMqFuKcl3UeDRORZEensTE19HTjE5biMS6w8hwmD5loDVi7EPS1avCYiFwN/IVpCe6SqvuJ2YC1hU1KTZ+U5TFjYYkT3tGrxmoj0AK4GHgeOBy4TkTdVNa17LBhvWHkOExa2GNEfLVm89j4wQVVfFBEBrgGuUNUTvAgwEWsppC7by3MYk80StRRakhS+o6pfNTrWQ1XXpTHGlFhSMMaY5KW0ollErgNQ1a9EpHEdIpuSaowxGSjR7KMRDT4vb/TYj1yIxRhjjM8SJQVp4vN4941xRTbXaAryftMmcyWafaRNfB7vvslCscHqDRs20L1797QPVjfeq+HZZ9txzTUTPN+rwe2fsykNSzzcd999rp/PGEgw0Cwie4FviLYK8omuUcC5n6eqB3kSYQI20OwftzfXiUQiFBUd6ft6Cj82EcrPz4/bIsrLyyMSibhyTpNdUlqnoKo57oVkwqxhUb19b9jbqKyEoUMHp+UNO1GNpmeeidZocntzIy9+zniqqqqa3G/aGLe1pMyFMfvxoqheEGo0+VU80Aq+GT9ZUjBJ8+INOwg1mvxMTFbwzfjFkoJJmhdv2MOHD6e6OpfKyv2PV1ZCdXUuw4cPb/U5mnPUUUfx/vvxNxl0OzFZwTfjlxYVxAsqG2j2R3NF9d57r4rnnnuu1bN1/BjkbXju888/l23bvqS8HCseaDJKq8pcuEFEaoCvgb3AHlXtKyIdgEeBrkANcJGqfpHoeVJJCn5NL8w0Tb1h33HHPVx//aS0vZH7UaOp4cynww+Hm2+Gbt2gRw9YswY+/bSQp5/2dlqsMekU1KTQV1U/a3BsGvC5qt4uIjcAh6nq9YmeJ9mk4OeVZyZq/IY9ZMgQjj++m+/TSFtr/vz5TJ9+FVOnbgNg5054+WWorYVly9py0033M2bMGH+DNKYVWlU620M/BYqdz+cCLwEJk0Iy/JpemMny8/P3mxY6f/5836eRpkPjAea2beGcc6Kf79q1i02bNvkUmTHu82ugWYHFIrJKRK50jnVS1VoA5+N3432jiFwpIitFZGVdXV2LT2h7E7sv0WydoqJtLFu2zOOIUhOEmU/G+MWvpHCGqp4KDAEmiMjAln6jqj6gqn1Vte8RRxzR4hMGYd57pkv0ZvrBB7BgwbxQ1C4KwsynbGd1n/zjS1JQ1c3Ox0+BRcBpwCci0hnA+fhpOs9pV3/uGz58OB988G3cN9ONG+G44yQULbLY7nQzZhRSXt6e2bOF8vL2zJhRaLvTeaRh3SfjLc8HmkWkHdBGVb92Pn8emAwMArY0GGjuoKrXJXquZAaabW9ib4wbN46Kilkcfzwceyx8+CFUVcGf/gSvvCL07DmZm2++2e8wW8R2p/Oe1X3yRtAGmjsBi6I7e5ILzFfVv4vICuCvIjIW+AhovLFPq3ixN7FNd4WBAwfy5pv/zbnnfkNtLQweDAMHwsEHw9y57RgyJDwtssYD6dmstraWESNG8Oijj7pabsPqPvkv6xavuXX1Z9Ndo5pqkS1dCvfcU8CkSb+jZ8+eWZkwEwn6BcWvf/1rZs2axVVXXeV6Ge958+ZxxRVXcPDBB7Nz507mzJnDpZde6uo5s03g1imkS1BWNAelzHNQNE6Qq1fnsX79Dnr1yqdnz0jWJsymBPmCwo/unJKSEiorK7ntttu48cYbOfPMM1myZIkr58pWKe3RbFrOprvur1+/flRXb6asbBbduv2ejz8WbrpJmTZtO2PHKlOnbmPixK0MHTo4FLOR3NRw/czUqdt8eX0SzfSpqqqitLSU/Px8IJokxowZQ3V1tWvxWN0nf1lSSAOb7nqgWH98UVERHTvupaoKXngBdu2KPp6tCbOxIFxQJJrp40cZ77lz51JcXAxEWw1z5sxx7VzmQJYU0sCmu8a3YsUKJk0aT/v2O9m9GxYvhlGj4P33o49na8JsyM8Livz8fESEmTNn8u233zJz5kxEpL5VEGNlvFMXxvUWlhTSwBY7HSjWLfK73+3g7rth7FiYNg0mTYoWmNu1K7sTZoyfFxQt7RrKhu4ct968Q7neQlVDe+vTp48GxfLly7VTp0Lt37+9XnqpaP/+7bVTp0Jdvny536H5oqKiQvv3b69Ll3LArV8/9IIL0E6dCjUSifgdqq8ikYh26lSoU6bs/xpNmeLN6/Pwww9rbm6uFhQUaE5Ojs6bN8/V8wXV+PHjtU2bNjdRqL4AABafSURBVDp+/Pi0PF9eXp4SLeez3y0vLy8tz99awEpt4n3VWgpp0nBwtWfPyZSVzaKmptb32SN+SdQt0qMHPPdcnq0Oxv/V0253DQW9+6SlXWjJ8mOAPl2CVCXVM+mYEx7vOWyx0z7du3fn2WfbAdsOeOzDD/OYMWNm1ibMxmIXFLH1M0OGeLd6OtY1VFxcTN++fZk7d25an79h94nb6xtS4dZiudgAfUVFBQUFBezcuTM8+2w31YQIwy2V7qN0dPNYV1Hz/O4WMf4KevdJQ251oRUXF2tubq5OmzZNc3NztaSkJC3Pmw4k6D7y/Y29Nbdkk8L27duTfqPavn27VlRU6OTJk7WiokI///xze7NrIUue2Wvz5s1aWlqq+fn5Cmh+fr6OGTNGa2trU36+gQMHpvz9ibj15j169GhdunSpqqouWbJES0tL0/K86ZAoKWTVmEKyc8JXrFhBUdGRTJ9+FR9+eAvTp19F9+5dOOKInbZQrQVsnCWxSCTC/PnzmTJlCvPnz8+ohXzpXt/g5iyexrOrjjjiiLSMg4R1vUVWjSkkMyc80U5td9wRnVJ58MGJn8OEo6icH3WHGpe2ePbZdlxzzYRAlLZIl4aD2DfeeCMPPvjgATWMmnvtG5fZmDlzJjNnzkxrmY2G4yglJSU89thjLFy4MLDjIG7LqpZCMnPCE7UqunWDeJuI2bz78InXGuzatTMrVqxw7ZxBKG3hhebWN7TktfdyFo9bM5HCJquSQjKLzBK1Kk44AV5/vfnnMMHm15tzEEpbeCFR90lLX3svy2yEeRppOmVVUkhmTniiVsWGDQWsWFFgu3KFnF9vzkGpleXnmEYyr71XZTb8qPMURFmVFKDlg5+JWhUbNx5MTc0mG0ANOb/enINQK8uPbrOGknntvSyzYXWesmygOaYlg5/N7dR22GGHBX4A1SSWaIHd+vXu7BIXiUTYtWsXb7+9h7/8BX71q30TFrzqgkw0iWLo0MGe7P+RzGvfeCC4pKTEtbjcXswXBrbJTjNsn97M5fW+3QduPqRUVUW3K926tb1nG+vMnz+f6dOvYurUA9+Qy8vbU1Y2y/ULHtsz3V9B26M5VMIwpdKkxot9u2PiX51H3wTvuiuPe+6ZzogRIzx5IwzCmEZeXh533HEP48ZdzrHHKieeCGvXwocfCvfff48lBB9ZUjBZzau6Q4kGVp95JpeDDz7YszdCP7rNGotEIlx//SSuu07Zuxdqa2HoUMjJUa6/fhIXX3yxJQaf+JYURCQHWAl8rKrni0gR8AjQAXgDuExVd/kVn8keXrQGg3B1HjN8+HCuuWYClZUc0Grxalp1LEnGGx5YvDg6+8ha6P7ws6VwNfAe8B3n/h3A3ar6iIjcD4wFZvoVnAk+P1YipyoIV+cxXnabNSVISdLsz5ekICJdgJ8A/w+4RkQEOBuIXRrMBW4lhaSwe/duNm3alDGrQjNBXl4eXbp04aCDDkrbc4atTEQQrs4b8rNcNwQrSaYqTBclyfBl9pGILASmAocAvwPGAK+r6jHO40cBz6nqiXG+90rgSoCjjz66z8aNG/d7vLq6mkMOOYSOHTsSzTXGT6rKli1b+PrrrykqKkrLc0YiEYqKjgzdzJXGiazh1XkQE5mbwj77KOy/y0DNPhKR84FPVXWViBTHDsf50rjZSlUfAB6A6JTUxo9H/9i6WkIICBGhY8eO1NXVpe05Ew/aBrc/2u+r8yAJQhdWqlJd5xGWloUf3UdnAMNE5MdAHtExhXuAQhHJVdU9QBdgc6onsIQQLOn+fYS5P9qmOO8T1iSZykVJmLo7PU8KqloOlAM4LYXfqeolIvIYcAHRGUilwN+8js2EQyb0R5uoMCbJZC9KgrCCPBlBqn10PdFB5/VAR2C2Vyd2Y3PxRYsWISK8//77cR8fM2YMCxcubPHzbd68mQsuuACA1atX8+yzz9Y/9tJLL/Hqq68mHWPXrl357LPPkv4+vyVT7daYdEu2dlXYquL6mhRU9SVVPd/5vEpVT1PVY1T1QlXd6VUcbuzqtGDBAgYMGMAjjzySluc78sgj65NIupJCWCVT7daYdEv2oiRs3Z1ZvaLZrV2dtm3bxiuvvMLSpUsZNmwYt956K6pKWVkZS5YsoaioiIazvrp27cqoUaNYunQpu3fv5oEHHqC8vJz169dz7bXXMm7cOGpqajj//PN54403+MMf/kAkEqGyspKRI0dy//33k5OTw3//938zffp0jjvuOMaNG8dHH30EwD333MMZZ5zBli1bGDlyJHV1dZx22mmEue5VWPujTfglO0geuu7OpjZvDsOtT58+B2xI/e6777Zo42rV9G8uHjNv3jy94oorVFX19NNP11WrVunjjz+u55xzju7Zs0c//vhjPfTQQ/Wxxx5TVdXvf//7et9996mq6qRJk7RXr1761Vdf6aeffqpHHHGEqqpWV1frCSecoKqqDz30kE6YMKH+fLfccoveeeed9fdHjhypL7/8sqqqbty4UY877jhVVS0rK9M//vGPqqr69NNPK6B1dXWt+llbKpnfizFhsH37dq2oqNApU6ZoRUWFRiKRuF8XiUS0U6dCnTIFXbp0323KFLRTp8Imv89NwEpt4n01q1sKsU01KioqKCgoYOfOnWnZVGPBggVMmjQJgBEjRrBgwQJ2797NyJEjycnJ4cgjj+Tss8/e73uGDRsGQK9evdi2bRuHHHIIhxxyCHl5eWzdujWp87/wwgu8++679fe/+uorvv76a5YtW1bff/mTn/yEww47rDU/pjFZraWD5GGbfpvVSQFatrl4MrZs2cKSJUtYu3YtIsLevXsREX7+858nnJrZtm1bANq0aVP/eez+nj17korh22+/5bXXXou7t2yYp+uGZZ63MY2FqbszSLOPfJHuXZ0WLlzI6NGj2bhxIzU1Nfzzn/+kqKiIDh068Mgjj7B3715qa2tZunRpyuc45JBD+Prrr5u8P3jwYGbMmFF/f/Xq1QAMHDiQiooKAJ577jm++OKLlGPwmt87hRnTWrGWxc0338yoUaMCmRDAWgpp39VpwYIF3HDDDfsd+8UvfsF7771Hjx496NWrF8ceeyxnnXVWyucoKSnh9ttvp3fv3pSXlzN06FAuuOAC/va3vzF9+nTuvfdeJkyYwEknncSePXsYOHAg999/P7fccgsjR47k1FNP5ayzznJ1W8N0Cts8b2PCLON2Xnvvvfc4/vjjfYrINKU1v5cg7BRmTCZJVPso67uPTPCFbZ63MWFmScEEXrIrSI0xqbOkYALPyloY452sH2g2wRe2ed7GhJklBRMKYZrnbUyYZX1SsAVR4RHGMsvGhE1Wjym4tSBKRPjtb39bf/+uu+7i1ltvTfg9//M//7NfaYpUJFsK+8knn+T222+Pe/45c+aweXNy+xzV1NRw4okH7KBqjAmRrE0KDRdETZ26jbFjlalTtzFx4laGDh28X/XUZLVt25YnnngiqTfodCSFZA0bNqx+oV06koIxJvyyNim4ufFFbm4uV155JXffffcBj23cuJFBgwZx0kknMWjQID766CNeffVVnnzySa699lp69+7Nhg0b9vuep556ih/84AeccsopnHPOOXzyySdAtM7S4MGDOeWUU7jqqqvqS2HX1NRw3HHH8ctf/pITTzyRSy65hBdeeIEzzjiDHj16sHz5ciD6xj9x4sQDzn/HHXewcuVKLrnkEnr37k0kEmHVqlWcddZZ9OnTh/POO4/a2loAVq1axcknn8zpp5/OX/7yl5RfM2NMMGRtUnB7QdSECROoqKjgyy+/3O/4xIkTGT16NGvWrOGSSy7hN7/5DT/84Q8ZNmwYd955J6tXr6Z79+77fc+AAQN4/fXXefPNNxkxYgTTpk0D4I9//CMDBgzgzTffZNiwYfX7JwCsX7+eq6++mjVr1vD+++8zf/58Kisrueuuu7jtttv2e/7G57/++uvp27cvFRUVrF69mtzcXMrKyli4cCGrVq3iiiuu4KabbgLg8ssv59577+W1115r1etljAmGrB1odnvji+985zuMHj2ae++9d79qpa+99lp9K+Syyy7juuuua/a5Nm3axMUXX0xtbS27du2iqKgIIGEp7KKiInr16gXACSecwKBBgxARevXqRU1NTVI/ywcffMDatWs599xzAdi7dy+dO3fmyy+/ZOvWrfV1nC677DKee+65pJ7b+MsmWpjGsral4MWCqEmTJjF79my++SZ+iwRaVsq6rKyMiRMn8vbbbzNr1qz9xjua+v7G5bcbluZOthS3qnLCCSewevVqVq9ezdtvv83ixYtR1VCX4s52VnnWxJO1ScGLfX47dOjARRddxOzZs+uP/fCHP6zft7miooIBzqBG4/LXDX355Zd873vfA/av6prOUtiJynH37NmTurq6+i6i3bt3884771BYWMihhx5KpZNZY7GY4HNzooUJN8+TgojkichyEXlLRN4RkT86x4tE5B8isk5EHhWRg92OJbYgqqxsFj17TqasbBY1NbX069cvbef47W9/u98spHvvvZeHHnqIk046iXnz5vHnP/8ZiO7Qduedd3LKKaccMNB86623cuGFF3LmmWdy+OGH1x+/5ZZbWLZsGaeeemqr94JofP4xY8Ywbtw4evfuzd69e1m4cCHXX389J598Mr179+bVV18F4KGHHmLChAmcfvrpcTf1McHk5kQLE26el86WaH9DO1XdJiIHAZXA1cA1wBOq+oiI3A+8paozEz2Xlc4OD/u9BMuUKVP48MNbGDv2wP//2bOFnj0nc/PNN/sQmfFCoEpnO/tGx0Z3D3JuCpwNLHSOzwV+5nVsxmQLqzxrmuLLmIKI5IjIauBT4HlgA7BVVWMjoJuA7zXxvVeKyEoRWVlXV+dNwMZkGKs8a5riy5RUVd0L9BaRQmAREK9fIW6/lqo+ADwA0e6jJr7GZsUESJh398tUVnnWNMXXdQqqulVEXgL6A4Uikuu0FroAKdVYyMvLY8uWLXTs2NESQwCoKlu2bLE3mQCyyrMmHs+TgogcAex2EkI+cA5wB7AUuAB4BCgF/pbK83fp0oVNmzZhXUvBkZeXR5cuXfwOw8RhlWdNY360FDoDc0Ukh+iYxl9V9WkReRd4RET+BLwJzE70JE056KCD6lf8GmOMSY7nSUFV1wCnxDleBZzmdTzGGGP2ydoVzcYYYw5kScEYY0w9z1c0p5OI1AEbU/z2w4GW74LjvzDFG6ZYweJ1U5hihXDF25pYv6+qR8R7INRJoTVEZGVTy7yDKEzxhilWsHjdFKZYIVzxuhWrdR8ZY4ypZ0nBGGNMvWxOCg/4HUCSwhRvmGIFi9dNYYoVwhWvK7Fm7ZiCMcaYA2VzS8EYY0wjlhSMMcbUy4qkEKQtQFvK2XPiTRF52rkf5FhrRORtEVktIiudYx1E5Hkn3udF5DC/4wQQkUIRWSgi74vIeyJyeoBj7em8prHbVyIyKajxAojIfzj/Y2tFZIHzvxfIv10RudqJ8x0RmeQcC8xrKyIPisinIrK2wbG48UnUvSKyXkTWiMipqZ43K5ICsBM4W1VPBnoDPxKR/kSrs96tqj2AL4CxPsbY2NXAew3uBzlWgBJV7d1g3vQNwItOvC8694Pgz8DfVfU44GSir3EgY1XVD5zXtDfQB9hOdP+RQMYrIt8DfgP0VdUTgRxgBAH82xWRE4FfEa23djJwvoj0IFiv7RzgR42ONRXfEKCHc7sSSLiVcUKqmlU3oAB4A/gB0dWAuc7x04H/9Ts+J5Yuzi/8bOBpQIIaqxNPDXB4o2MfAJ2dzzsDHwQgzu8A1TgTLIIca5zYBwOvBDleorsl/hPoQLTY5tPAeUH82wUuBP6rwf3fA9cF7bUFugJrG9yPGx8wCxgZ7+uSvWVLS6FVW4D64B6if6DfOvc7EtxYIbpL3mIRWSUiVzrHOqlqLYDz8bu+RbdPN6AOeMjpmvsvEWlHMGNtbASwwPk8kPGq6sfAXcBHQC3wJbCKYP7trgUGikhHESkAfgwcRUBf2waaii+WkGNSfp2zJimo6l6NNsO7EG0ytngLUC+JyPnAp6q6quHhOF/qe6wNnKGqpxJtwk4QkYF+B9SEXOBUYKaqngJ8Q0C6XhJx+uCHAY/5HUsiTv/2T4Ei4EigHdG/icZ8/9tV1feIdms9D/wdeAvYk/Cbgi1t7xFZkxRiVHUr8BINtgB1Hkp5C9A0OwMYJiI1RHehO5toyyGIsQKgqpudj58S7fM+DfhERDoDOB8/9S/CepuATar6D+f+QqJJIoixNjQEeENVP3HuBzXec4BqVa1T1d3AE8APCejfrqrOVtVTVXUg8DmwjuC+tjFNxbeJaEsnJuXXOSuSgogcISKFzuexLUDfY98WoNCKLUDTSVXLVbWLqnYl2mWwRFUvIYCxAohIOxE5JPY50b7vtcCTROOEgMSrqv8C/ikiPZ1Dg4B3CWCsjYxkX9cRBDfej4D+IlIgIsK+1zeof7vfdT4eDQwn+hoH9bWNaSq+J4HRziyk/sCXsW6mpPk94OPRYM1JRLf4XEP0DesPzvFuwHJgPdGmeVu/Y20UdzHwdJBjdeJ6y7m9A9zkHO9IdLB8nfOxg9+xOnH1BlY6fwv/AxwW1FideAuALcChDY4FOd4/Au87/2fzgLYB/tt9mWjSegsYFLTXlmiSqgV2E20JjG0qPqLdR38hOlb6NtEZYCmd18pcGGOMqZcV3UfGGGNaxpKCMcaYepYUjDHG1LOkYIwxpp4lBWOMMfUsKZiMJCJ7G1UY9WzlcrzqlsaEhU1JNRlJRLapanufzj0Q2AY8rNFqoV6cM0dV93pxLpPZrKVgsoaIHCoiH8RWNDv1/n/lfD5TRFZKg/02nOM1InKbiLzmPH6qiPyviGwQkXHxzqOqy4iWTUgUy4VOLf+3RGSZcyxHRO6S6N4Ua0SkzDk+yCng97bTCmnbILY/iEglcKGIdBeRvzuFCV8WkePS8bqZ7JLb/JcYE0r5TlXcmKmq+qiITATmiMifgcNU9T+dx29S1c9FJAd4UUROUtU1zmP/VNXTReRuojXuzwDyiK7gvj/F+P4AnKeqH8dKsBCtg18EnKKqe5wNVfKccw5S1Q9F5GFgPNF6WAA7VHUAgIi8CIxT1XUi8gPgPqK1s4xpMUsKJlNFNFoVdz+q+ryIXEi0JMDJDR66yCn7nUu0Tv2/Ey2FAdG6MhAtH9BeVb8GvhaRHSJSqNEii8l6hWhy+ivRwnEQrcl1vzplpp0kdTLRInMfOl8zF5jAvqTwKICItCdafO6xaNkhIFpiwpikWFIwWUVE2hAtmx4huhnMJhEpAn4H9FPVL0RkDtGWQMxO5+O3DT6P3U/pf0hVxzlX8z8BVotIb6L1axoP8sUridzQN87HNkT3LTggERqTDBtTMNnmP4hWyB0JPCgiBxHdke0b4EsR6UT8PQDSSkS6q+o/VPUPRHcmOwpYDIyLlZkWkQ5Ei8t1FZFjnG+9DPi/xs+nql8B1U4rKLZn78mNv86Y5lhSMJkqv9GU1NtF5Fjgl8BvVfVlYBlws6q+RbSK7jvAg0S7dlImIguA14CeIrJJROLtSXynM3C81onjLeC/iJafXiMibwGjVHUHcDnRbqG3ibZOmhrHuAQY63zvO0Q3vDEmKTYl1RhjTD1rKRhjjKlnScEYY0w9SwrGGGPqWVIwxhhTz5KCMcaYepYUjDHG1LOkYIwxpt7/ByXhUjvRUQGzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotData(X,y)\n",
    "plt.xlabel(\"Exam 1 score\")\n",
    "plt.ylabel(\"Exam 2 score\")\n",
    "plt.legend([\"Admitted\", \"Not admitted\"])\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>1.2 Implementation</h4>\n",
    "First we construct the sigmoid function defined as:\n",
    "$$h_\\theta(x) = g(\\theta^Tx)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(z):\n",
    "    \"\"\"\n",
    "    Compute sigmoid function given the input z.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    z : array_like\n",
    "        The input to the sigmoid function. This can be a 1-D vector \n",
    "        or a 2-D matrix. \n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    g : array_like\n",
    "        The computed sigmoid function. g has the same shape as z, since\n",
    "        the sigmoid is computed element-wise on z.\n",
    "    \"\"\"\n",
    "    # convert input to a numpy array\n",
    "    z = np.array(z)\n",
    "\n",
    "    g = 1 + np.exp(-1*z)\n",
    "    g = np.reciprocal(g)\n",
    "\n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sigmoid of 0 is  0.5\n",
      "Sigmoid of 100 is  1.0\n",
      "Sigmoid of -100 is  3.7200759760208356e-44\n"
     ]
    }
   ],
   "source": [
    "# Check a few values of sigmoid\n",
    "print(\"Sigmoid of 0 is \",sigmoid(0))\n",
    "print(\"Sigmoid of 100 is \",sigmoid(100))\n",
    "print(\"Sigmoid of -100 is \",sigmoid(-100))\n",
    "\n",
    "# sigmoid of 0 should be exactly 0.5\n",
    "# sigmoid of large positive numbers should be close to 1\n",
    "# sigmoid of large negative numbers should be close to 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a working sigmoid function, we can now implement a cost function which returns the cost and gradient for cost defined as:\n",
    "\n",
    "$$\\begin{align}\n",
    "J(\\theta) & = \\dfrac{1}{m} \\sum_{i=1}^m \\mathrm{Cost}(h_\\theta(x^{(i)}),y^{(i)}) \\\\\n",
    "& = - \\dfrac{1}{m} [\\sum_{i=1}^{m} y^{(i)} \\log(h_\\theta(x^{(i)})) + (1 - y^{(i)}) \\log(1-h_\\theta(x^{(i)}))] \\\\\n",
    "\\end{align}$$\n",
    "\n",
    "and derivative:\n",
    "\n",
    "$$\\frac{\\partial}{\\partial \\theta_j} J(\\theta) = \\dfrac{1}{m} \\sum_{i=1}^{m} (h_\\theta(x^{(i)}) - y^{(i)}) x_j^{(i)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the data matrix appropriately, and add ones for the intercept term\n",
    "m, n = X.shape\n",
    "\n",
    "# Add intercept term to X\n",
    "X = np.concatenate([np.ones((m, 1)), X], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def costFunction(theta,X,y):\n",
    "    \"\"\"\n",
    "    Compute cost and gradient for logistic regression. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    theta : array_like\n",
    "        The parameters for logistic regression. This a vector\n",
    "        of shape (n+1, ).\n",
    "    \n",
    "    X : array_like\n",
    "        The input dataset of shape (m x n+1) where m is the total number\n",
    "        of data points and n is the number of features. We assume the \n",
    "        intercept has already been added to the input.\n",
    "    \n",
    "    y : arra_like\n",
    "        Labels for the input. This is a vector of shape (m, ).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    J : float\n",
    "        The computed value for the cost function. \n",
    "    \n",
    "    grad : array_like\n",
    "        A vector of shape (n+1, ) which is the gradient of the cost\n",
    "        function with respect to theta, at the current values of theta.\n",
    "    \"\"\"\n",
    "    ## Initialize some useful values\n",
    "    m = y.size  # number of training examples\n",
    "    J = 0\n",
    "    grad = np.zeros(theta.shape)\n",
    "    h = sigmoid(X.dot(theta))\n",
    "    logh = np.log(h)\n",
    "    tempLog = np.log(1-h)\n",
    "    yTrans = y.transpose()\n",
    "    Xtrans = X.transpose()\n",
    "    tempTrans = (1-y).transpose()\n",
    "    \n",
    "    \n",
    "    J = ((-yTrans).dot(logh))\n",
    "    J = J - tempTrans.dot(tempLog)\n",
    "    J = J * (1/m)\n",
    "    \n",
    "    diff = np.subtract(sigmoid(X.dot(theta)),y)\n",
    "    grad = Xtrans.dot(diff)\n",
    "    grad = grad * (1/m)\n",
    "    \n",
    "    # =============================================================\n",
    "    return J, grad\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now test our cost function with varying initial thetas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cost at initial theta (zeros): 0.693\n",
      "Expected cost (approx): 0.693\n",
      "\n",
      "Gradient at initial theta (zeros):\n",
      "\t[-0.1000, -12.0092, -11.2628]\n",
      "Expected gradients (approx):\n",
      "\t[-0.1000, -12.0092, -11.2628]\n",
      "\n",
      "Cost at test theta: 0.218\n",
      "Expected cost (approx): 0.218\n",
      "\n",
      "Gradient at test theta:\n",
      "\t[0.043, 2.566, 2.647]\n",
      "Expected gradients (approx):\n",
      "\t[0.043, 2.566, 2.647]\n"
     ]
    }
   ],
   "source": [
    "# Initialize fitting parameters\n",
    "initial_theta = np.zeros(n+1)\n",
    "\n",
    "cost, grad = costFunction(initial_theta, X, y)\n",
    "\n",
    "print('Cost at initial theta (zeros): {:.3f}'.format(cost))\n",
    "print('Expected cost (approx): 0.693\\n')\n",
    "\n",
    "print('Gradient at initial theta (zeros):')\n",
    "print('\\t[{:.4f}, {:.4f}, {:.4f}]'.format(*grad))\n",
    "print('Expected gradients (approx):\\n\\t[-0.1000, -12.0092, -11.2628]\\n')\n",
    "\n",
    "# Compute and display cost and gradient with non-zero theta\n",
    "test_theta = np.array([-24, 0.2, 0.2])\n",
    "cost, grad = costFunction(test_theta, X, y)\n",
    "\n",
    "print('Cost at test theta: {:.3f}'.format(cost))\n",
    "print('Expected cost (approx): 0.218\\n')\n",
    "\n",
    "print('Gradient at test theta:')\n",
    "print('\\t[{:.3f}, {:.3f}, {:.3f}]'.format(*grad))\n",
    "print('Expected gradients (approx):\\n\\t[0.043, 2.566, 2.647]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a working cost function, we can implement gradient descent using a built in optimization function scipy.optimize.\n",
    "\n",
    "To use this function we need to pass in:\n",
    "\n",
    "- The initial values of the parameters we are trying to optimize\n",
    "- A function that, when given training set and theta, computes the logistic regression cost and gradient with respect to theta for (X,y)\n",
    "- jac: which is an indication if we would like the function to return the jacobian (gradient) as well\n",
    "- method: which is the method/algorithm we would like to implement\n",
    "- options: options specific to our chosen algorithm (chosen iterations in our case)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cost at theta found by optimize.minimize: 0.203\n",
      "Expected cost (approx): 0.203\n",
      "\n",
      "theta:\n",
      "\t[-25.161, 0.206, 0.201]\n",
      "Expected theta (approx):\n",
      "\t[-25.161, 0.206, 0.201]\n"
     ]
    }
   ],
   "source": [
    "# set options for optimize.minimize\n",
    "options= {'maxiter': 400}\n",
    "\n",
    "# see documention for scipy's optimize.minimize  for description about\n",
    "# the different parameters\n",
    "# The function returns an object `OptimizeResult`\n",
    "# We use truncated Newton algorithm for optimization which is \n",
    "# equivalent to MATLAB's fminunc\n",
    "# See https://stackoverflow.com/questions/18801002/fminunc-alternate-in-numpy\n",
    "res = optimize.minimize(costFunction,\n",
    "                        initial_theta,\n",
    "                        (X, y),\n",
    "                        jac=True,\n",
    "                        method='TNC',\n",
    "                        options=options)\n",
    "\n",
    "# the fun property of `OptimizeResult` object returns\n",
    "# the value of costFunction at optimized theta\n",
    "cost = res.fun\n",
    "\n",
    "# the optimized theta is in the x property\n",
    "theta = res.x\n",
    "\n",
    "# Print theta to screen\n",
    "print('Cost at theta found by optimize.minimize: {:.3f}'.format(cost))\n",
    "print('Expected cost (approx): 0.203\\n');\n",
    "\n",
    "print('theta:')\n",
    "print('\\t[{:.3f}, {:.3f}, {:.3f}]'.format(*theta))\n",
    "print('Expected theta (approx):\\n\\t[-25.161, 0.206, 0.201]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have an optimal theta, we can use it to get a decision boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mapFeature(X1, X2, degree=6):\n",
    "    \"\"\"\n",
    "    Maps the two input features to quadratic features used in the regularization exercise.\n",
    "\n",
    "    Returns a new feature array with more features, comprising of\n",
    "    X1, X2, X1.^2, X2.^2, X1*X2, X1*X2.^2, etc..\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    X1 : array_like\n",
    "        A vector of shape (m, 1), containing one feature for all examples.\n",
    "\n",
    "    X2 : array_like\n",
    "        A vector of shape (m, 1), containing a second feature for all examples.\n",
    "        Inputs X1, X2 must be the same size.\n",
    "\n",
    "    degree: int, optional\n",
    "        The polynomial degree.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    : array_like\n",
    "        A matrix of of m rows, and columns depend on the degree of polynomial.\n",
    "    \"\"\"\n",
    "    if X1.ndim > 0:\n",
    "        out = [np.ones(X1.shape[0])]\n",
    "    else:\n",
    "        out = [np.ones(1)]\n",
    "\n",
    "    for i in range(1, degree + 1):\n",
    "        for j in range(i + 1):\n",
    "            out.append((X1 ** (i - j)) * (X2 ** j))\n",
    "\n",
    "    if X1.ndim > 0:\n",
    "        return np.stack(out, axis=1)\n",
    "    else:\n",
    "        return np.array(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotDecisionBoundary(plotData, theta, X, y):\n",
    "    \"\"\"\n",
    "    Plots the data points X and y into a new figure with the decision boundary defined by theta.\n",
    "    Plots the data points with * for the positive examples and o for  the negative examples.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    plotData : func\n",
    "        A function reference for plotting the X, y data.\n",
    "\n",
    "    theta : array_like\n",
    "        Parameters for logistic regression. A vector of shape (n+1, ).\n",
    "\n",
    "    X : array_like\n",
    "        The input dataset. X is assumed to be  a either:\n",
    "            1) Mx3 matrix, where the first column is an all ones column for the intercept.\n",
    "            2) MxN, N>3 matrix, where the first column is all ones.\n",
    "\n",
    "    y : array_like\n",
    "        Vector of data labels of shape (m, ).\n",
    "    \"\"\"\n",
    "    # make sure theta is a numpy array\n",
    "    theta = np.array(theta)\n",
    "    \n",
    "    # Plot the data (note: first collumn is x-intercepts so we can ignore it)\n",
    "    plotData(X[:,1:3],y)\n",
    "    \n",
    "    if X.shape[1] <= 3:\n",
    "        # Only need 2 points to define line, so we choose the two endpoints\n",
    "        plot_x = np.array([np.min(X[:, 1]) - 2, np.max(X[:, 1]) + 2])\n",
    "        \n",
    "        # Calculate the decision boundary line ( given form y = theta0*x0 + \n",
    "        #   theta1*x1 + theta2*x2, we just solve for y)\n",
    "        plot_y = (-1. / theta[2]) * (theta[1] * plot_x + theta[0])\n",
    "        \n",
    "        # Plot and adjust axes\n",
    "        plt.plot(plot_x, plot_y)\n",
    "        \n",
    "        # Setup legend\n",
    "        plt.legend(['Admitted', 'Not admitted', 'Decision Boundary'])\n",
    "        plt.xlim([30, 100])\n",
    "        plt.ylim([30, 100])\n",
    "        \n",
    "    else:\n",
    "        # Setup grid range\n",
    "        u = np.linspace(-1, 1.5,50)\n",
    "        v = np.linspace(-1,1.5,50)\n",
    "        \n",
    "        z = np.zeros((u.size, v.size))\n",
    "        # Evaluate z = theta*x over the grid\n",
    "        for i, ui in enumerate(u):\n",
    "            for j, vj in enumerate(v):\n",
    "                z[i, j] = np.dot(mapFeature(ui, vj), theta)\n",
    "                \n",
    "        z = z.T  # important to transpose z before calling contour\n",
    "        # Plot z = 0\n",
    "        plt.contour(u, v, z, levels=[0], linewidths=2, colors='g')\n",
    "        plt.contourf(u, v, z, levels=[np.min(z), 0, np.max(z)], cmap='Greens', alpha=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeVhVVdfAf5vxAs6mKJqKijOIiqU5kqlNWtlkVg7Vq2na9GZlk6WlZX1N+mpW5hSWZVlqWmZSSpNA4ZA5MuRAiTghXkFgfX9cIEDmO5172b/nOc+955x7zl73cNlrr7XXWluJCBqNRqPRVBYPZwug0Wg0GtdCKw6NRqPRVAmtODQajUZTJbTi0Gg0Gk2V0IpDo9FoNFVCKw6NRqPRVAm7KQ6l1AdKqWNKqV1FjjVQSn2rlNqf/1o//7hSSr2tlDqglNqhlOpuL7k0Go1GYx32tDiWAFeXOPYk8J2IhADf5e8DXAOE5G/jgQV2lEuj0Wg0VmA3xSEiW4ATJQ7fACzNf78UuLHI8WVi4RegnlKqqb1k02g0Gk318XJwe4EikgogIqlKqcb5x5sBh4p87nD+sdSSN1BKjcdilRAQENCjQ4cOVgmUl5fHqVOnyMrKwtfXl/r166OUsuqeGo1GY2Ti4+OPi0ij6l7vaMVRFqX11KXWQhGRd4F3ASIiIiQuLq7ajcbGxjJs2BCCg3No2zaTAwcCSEo6ydq1G+nZs2e176vRaDRGRimVYs31jlYc/yilmuZbG02BY/nHDwOXFvlcc+CoPQUxm80MGzaEyZNP0bdvwdGzxMTAsGFDSE5OxWQy2VMEjUajcUkcHY67BhiT/34M8GWR46Pzo6t6AacLXFr2YvXq1QQH5xRRGhb69oXg4Bw+//xzezav0Wg0Los9w3E/An4G2iulDiul7gVeBgYrpfYDg/P3AdYDicAB4D1gkr3kKuDgwYO0bZtZ6rm2bTNJTEy0twgaTaVJTU1lwIAB/P33384WRaNBuXJZ9dLmOC5cuMDhw4c5f/58uddmZmaSkZFOvXoXf/9TpxS1azckICDApvKWRm5uLmlpaTRq1AhPT0+7t1cTMJlMNG/eHG9vb2eLYjMmTZrEwoULmTBhAvPnz3e2OBoXRykVLyIR1b7e3RRHUlIStWvXpmHDhuVGR+Xl5bFz53YaNcqldu1/j2dkQFqaJ6GhXfHwsL8nLyUlpVBxtGzZ0u7tuTsiQnp6OhkZGQQHBztbHKvx8/MrdRBkMpkwm81OkEjjDlirONyu5Mj58+crVBoAHh4etG3bjrQ0Tw4f9iAtjfxXT9q2bWd3pREfH09cXBxpaWkApKWlERcXR3x8vF3bdXeUUjRs2LBCi7MyGME9lJiYyJgxY/Dz8wMsimTs2LEkJSU5TSaNxu0UB1DpPIyAgABCQ7vSuHFLTKYgGjduSWhoV4e4qEJDQ2nYsGGhgvLw8KBhw4aEhYXZvW13p7p5OCUVxcyZM4mJiWHGjBm2FK9KNG3alEGDBnHhwgX8/f3Jzs5m0KBBNGnSxGkyVRYjKF6NfXBLxVEVCjrsoKCgYh25vfHx8aFOnTqICB4eHuTl5VGnTh238su7GgWKonnz5iilWLBgAXl5eSxYsAClVOGo39F88MEHADz//PMopQr3K4uzOnAjKF6NfagRiiMvL4/09HSOHj1Keno6eXl5dm9z9erVKKXYs2dPqefHjh3Lhx9+CEBQUBBKKY4fP17m/Y4ePcott9wCQEJCAuvXry889/333/PTTz9VWcZWrVqV22ZNwc/Pr5iiyM3Nvei8M91DLVq04Ntvv2Xq1Kls3LiRFi1aVOl6R3fgJZ+nsxWvxva4veLIzMxk587tHDuWwvnzRzl2LIWdO7eTmflvKK49RmQfffQRffv25eOPPy7zM15eXrRr144mTZoQEhKCj49PmZ8NCgpi1apVgO0Uh8ZCafMI/fr1w9PT0xDuoaVLlzJw4EAAIiMjWbJkSaWuc1YHrudl3B+3Vhx5eXkcOLCPRo1yad48j0aNyH/N5cCBfYWWh61HZGfPnuXHH39k0aJFhYpDRJg8eTKdOnXiuuuu49ixYwQGBlK7dm1atWrFyy+/zKhRo4iIiOC3335j6NChtGnThnfeeQeA5ORkunTpQnZ2Ns899xwrV64kPDycV155hXfeeYc33niD8PBwtm7dSlpaGjfffDM9e/akZ8+e/PjjjwCkp6czZMgQunXrxoQJE3DliDpbUto8wvHjx1FKVds9ZASc1YG78ryMppKIiMtuPXr0kJLs3r278P3x48dl9+54OXMm9qJt9+548fX1FSw1sYptJpPpovtWheXLl8s999wjIiK9e/eW+Ph4+eyzz+Sqq66SnJwcOXLkiNStW1c+/fRTERFp2bKlzJ8/X0REHn74YQkNDZUzZ87IsWPHpFGjRiIikpSUJJ07dxYRkcWLF8sDDzxQ2N706dPl1VdfLdy/4447ZOvWrSIikpKSIh06dBARkSlTpsgLL7wgIiLr1q0TQNLS0qz6rkal6O+gMgwcOFC8vLxkzpw54uXlJYGBgRIdHS0iIps3b5YxY8bYXkgHsGzZMvHy8hJ/f3/x9PSU5cuXO6Tdks8zMjLSIe1qKgcQJ1b0vW5tcVgq3pY+n+Hrm8fPP/9slxHZRx99xMiRIwEYOXIkH330EVu2bOGOO+7A09OToKAgrrzyymLXDB8+HLBEW11++eXUrl2bRo0aYTKZOHXqVJXa37RpE5MnTyY8PJzhw4dz5swZMjIy2LJlC3fddRcA1113HfXr17fqe7oTJecRrr766mq5h4yGtRPr1cXaeRmNsTFKdVy74Ovry+nTHsDFyiMry4MWLVowaNAgoqKi8Pf3Jysry2qTOj09nc2bN7Nr1y6UUuTm5qKU4qabbio3TNTX1xewRHkVvC/Yz8nJqZIMeXkWpViaL1uXjC+dpUuXFr6PjIwkMjLSidLYjoIOfODAgURERBT7nvbEXZ+nxoJbWxz169cnO1uRkVH8eEYGZGcr6tevb/MR2apVqxg9ejQpKSkkJydz6NAhgoODadCgAR9//DG5ubmkpqYSHR1d7TZq165NRpEvVXJ/yJAhzJs3r3A/ISEBgP79+xMVFQXAhg0bOHnyZLVl0LgG1Z1Y12jKw6UVh/lCbrnnK5MdbmuT+qOPPuKmm24qduzmm2/m77//JiQkhNDQUCZOnMiAAQOq3UZkZCS7d+8mPDyclStXMmzYMFavXl04Of72228TFxdHWFgYnTp1Kpxgnz59Olu2bKF79+7afaDRaKqNS9eq8m0aIrOWrmPylW3x9bIUCPzzzz/p2LFjsc/l5eVx8uTJYqv8OSrRT+McSvsdaDQaCzW6VlV9fx/mbj7A8Lk/svPw6TI/56zscI3GEejSHhpH49I9aPP6fnwwNoJT5mxunP8jr32zV+cmaGocurSHxtG4tOIAuLJDIBsfGcBN3ZoxL/oAxzKyOJddtSgkjcYV0aU9NM7C5RUHQF0/b167tSuLx/YkT+DgsUz+Pm0mT1sfGjfG1Ut7aBeb6+IWiqOAyA6NCazjS31/b45lZHHgn7Pa+tC4La5e2kO72FwXt1IcAB5K0byBP8GXBJArwsFjZ0k9bSYvT1sfGvfDWZnh1qBdbK6P2ymOAmqbvGkXWIv6/j6kZWSx/9hZzmVdbH2YzWZWrFjBzJkzWbFihU1WjlNK8d///rdw/7XXXuP5558v95ovvviC3bt3W9VuVcukr1mzhpdffrnU9pcsWcLRo0er1H5BIUaN43DF0h4FLjaTyQRYlsF1JRebxo0VB4Cnh0eh9ZEnwsG04tZHbGwswcFBzJ07gX37pjN37gRatWpKbGysVe36+vry+eefV6kTt4XiqCrDhw/nySefLLX96igOjeNxxczwAhdbVlYWYFnu2ZVcbBo3VxwFFFofAf9aH8dPZTBs2BAmTz7F7NlnufdeYfbss0yefIphw4ZYZXl4eXkxfvx43njjjYvOpaSkMGjQIMLCwhg0aBB//fUXP/30E2vWrGHq1KmEh4dz8ODBYtesXbuWyy+/nG7dunHVVVfxzz//AGWXSU9OTqZDhw7cd999dOnShTvvvJNNmzbRp08fQkJC2LZtG2BRDpMnT76o/VdeeYW4uDjuvPNOwsPDMZvNxMfHM2DAAHr06MHQoUNJTU0FLGund+3ald69e/O///2v2s9MY0zsMYHt5+fH6NGji4XO33333dpV5UpYU1q3uhvwELAL+AN4OP9YA+BbYH/+a/2K7lNRWfXSOGPOlt1HT8vst9+Vyy+vJdHRXLT16lVLoqKiyq9LXA4BAQFy+vRpadmypZw6dUpeffVVmT59uoiIXH/99bJkyRIREVm0aJHccMMNIiIyZsyYwjLrJTlx4oTk5eWJiMh7770njz76qIiUXSY9KSlJPD09ZceOHZKbmyvdu3eXcePGSV5ennzxxReFbRYtz16y/QEDBkhsbKyIiGRnZ0vv3r3l2LFjIiLy8ccfy7hx40REJDQ0VL7//nsREXnssccKS787m6qWVdeUzsSJE8XDw0MmTpxos3sePXpU2rRpU7isgY+Pj7Rt21ZSU1Nt1oamfHC1supKqS7Af4DLgK7A9UqpEOBJ4DsRCQG+y9+3OQXWx/HUQ4SEZJb6mbZtM0lMTLSqnTp16jB69GjefvvtYsd//vlnRo0aBVhGWTExMRXe6/DhwwwdOpTQ0FBeffVV/vjjD4Byy6QHBwcTGhqKh4cHnTt3ZtCgQSilCA0NJTk5uUrfZe/evezatYvBgwcTHh7Oiy++yOHDhzl9+jSnTp0qrLt19913V+m+GuNizwnspk2bMn36dHJzc/H39yc3N5fp06drV5UL4QxXVUfgFxE5JyI5wA/ATcANQEEt5qXAjfYSwNPDg+5dOnDgQECp5w8cCKB169ZWt/Pwww+zaNGiYsvUlqQyZc6nTJnC5MmT2blzJwsXLizmRivr+pKl2YuWba9qmXYRoXPnziQkJJCQkMDOnTvZuHEjIqLLtLsZBa6pX375xa45Iq4YDab5F2cojl1Af6VUQ6WUP3AtcCkQKCKpAPmvjUu7WCk1XikVp5SKS0tLq7YQI0aMICnJi5ID/pgYSEryYsSIEdW+dwENGjTgtttuY9GiRYXHrrjiisLlZKOioujbty9wcWn0opw+fZpmzZoBxdc5sGWZ9PJKtbdv3560tDR+/vlnAC5cuMAff/xBvXr1qFu3bqHVVCBLVcnOzmbPnj1cuHCh2vJrbENBbsXChQvtmiPiitFgmiJY4+eq7gbcC/wGbAHeAd4ATpX4zMmK7lOdOY6ibNu2TQID60mvXrXkrruUXHZ5LbmkUV356rstkpubV+n7lCQgIKDw/d9//y1+fn6FcxxJSUkSGRkpoaGhcuWVV0pKSoqIiMTExEjHjh0lPDxcDhw4UOx+X3zxhQQHB0vfvn3lsccekwEDBoiIZWncwYMHS7du3eThhx+WFi1aFM5xFJ1rKDp/UdYStCXbX7VqlbRr1066du0q586dk99//1369esnYWFh0qlTJ3n33XdFRCQuLk7CwsKkV69eMn369GrNcSQnJ0tsbKwkJydX+dqy0HMcVcNkMpW6jDKgl3+tIkePHpX+/fsbes4GK+c4nF5WXSk1CziMZcJ8oIikKqWaAt+LSPvyro2IiJC4uLhix6paTttsNrN69WoSExNp1SqY3oOuITPXA18vT5rX9yPA160XSQQsI/7ExETatGmDt7e3w9qNj48vtSilUooePXpYdW9dVr1qpKamMm3aND755BPMZjN+fn4EBQUxZ84cRowYQXR0NEuXLnWJcF9nM2nSJBYuXMiECROYP3++s8UpFWvLqjtFcSilGovIMaVUC2Aj0Bt4CkgXkZeVUk8CDUTk8fLuYwvFURpnz1/g8Ekz2bl5XFLLlyZ1THh4uK8vPyUlhbS0NBo1akTLli0d1m52djZHjhzh5MmT5OXl4eHhQf369WnevLnVCkwrjqqzfPly7rnnHnx8fMjKymLJkiWFwReaivHz8ys1jN9kMmE2m50gUdm46nocnymldgNrgQdE5CTwMjBYKbUfGJy/7xRqmbwJCaxNwwBfjp/NYv+xDDJLyTp3deLj44mLi6NgrigtLY24uDji4+Md0r6Pjw916tSxhPd5eJCXl0edOnUcavVo/kVPWFuHqxedrApOURwi0k9EOolIVxH5Lv9YuogMEpGQ/NcTzpCtAE8PRbP6frS+JAABDqad5egpM7luVPMqNDS02MJWBQtehYWFOUyGguz6oKAglFJVyrbX2JaKJqx1NdvilHwerl50sirUiMxxa6hl8iakcW0a1vrX+jjrYOvDXlFHRhjx+/j40K5dO5o0aUJISAg+Pj4Oa1tTnIrKl+hqtsUp7XnUFKtNK45K4OmhaFbPj9aX1AIg0cHWR2pqKmfPnrVL7Shnj/iDg4OpXbs2YEmaDA4Odmj7morR1WyLU97zqClhxlpxVIFaJi9CGtfmEgdZH46Yg9Ajfk1F1CTffWUo73k4q+iko92IWnFUEU8PRVA9P1o3+tf6OFLC+vD09CQ8PJzOnTvTtWtXXn/9dfLy8qrcVmhoKMuWLSus1lvaHMQ777zDsmXLqv19Ckb8ycnJBAYGctNNN9G1a1euuOIK9u7dW+37VpZatWrZvY3S0P76ylOTfPeVwYjPw+FuRGuSQJy9WZsAaC05uXly5OQ52X7opPyZeloyzBdEpHgC4D///CODBg2S5557rlptHD9+XOLi4iQ+Pl5iY2Pl+PHjNpG9JCWTBt955x0ZPXq0XdoqStFnVRF5eXmSm5tbqc9W9DuwR/E+d2bgwIHi5eWlkwHzMcrzKCtx02QylXsdrlbk0J24yPo4brE+itK4cWPeffdd5s2bh4iQm5vL1KlT6dmzJ2FhYSxcuLDws3PmzCE0NJSuXbsWrpNx3333sWnTJoKCgpg3bx6XX345YWFhPPbYY4BlEu61114DICEhgV69ehEWFsZNN91UWIJk4MCBPPHEE1x22WW0a9eOrVu3Vvjdzpw5U1g08fz584wbN47Q0FC6detGdHQ08G9Z9gKuv/56vv/+e8BiSTz99NN07dqVXr16FZaCT0pKonfv3vTs2ZNnn3228NqzZ88yaNAgunfvTmhoKF9++SVgKRHfsWNHJk2aRPfu3Zk5cyaPPPJI4XXvvfcejz76aIXfpwBX8dcbwSIqKkNN8d1XFqM8D2e5Ed06LfqFtX+w++gZm96zU1Adpg/rXOxYLV8v2jWuzd9nznP8bBYiliTCWiZLdFLr1q3Jy8vj2LFjfPnll9StW5fY2FiysrLo06cPQ4YMYc+ePXzxxRf8+uuv+Pv7c+KEJRrZ09OToKAgfHx8+PHHH9m4cSOtW7fm1KlTF8k2evRo5s6dy4ABA3juued44YUXePPNNwHIyclh27ZtrF+/nhdeeIFNmzZddP3BgwcJDw8nIyODc+fO8euvvwIUrrOxc+dO9uzZw5AhQ9i3b1+5zykzM5NevXrx0ksv8fjjj/Pee+/xzDPP8NBDDzFx4kRGjx5dbP0Ok8nE6tWrqVOnDsePH6dXr14MHz4csFTnXbx4MfPnzyczM5OwsDDmzJmDt7c3ixcvLqZ8KyIxMfGiDOnbb7+d2bNnV/oejqCo68FZ2cdFZShaIy0yMpLIyEinyGQUjPI8CtxmUVFR+Pv7k5WV5RC3mbY4bIRHvvXRptD6yOTIyXOFcx+Sn6G/ceNGli1bRnh4OJdffjnp6ens37+fTZs2MW7cOPz9/QFLgUSwjNz9/PyoU6cO/v7+vPTSS3z++eeFnyugZInzMWPGsGXLlsLzBUUbe/ToUWZZ9TZt2pCQkMDBgwd58803GT9+PAAxMTGFJdM7dOhAy5YtK1QcPj4+XH/99Re1+eOPP3LHHXcAxcuwiwhPPfUUYWFhXHXVVRw5cqTQSmnZsiW9evUCICAggCuvvJJ169YVhiiHhoaWK0tRjOifLoozLaICC8NkMrmEVaax4IwQYLe2OEpaBo4gwNcLpSiMvMo4n8OFU6l4enrSuHFjRIS5c+cydOjQYtd9/fXX5ZYo9/LyYtu2bXz33Xd8/PHHzJs3j82bN1daroKy6p6enpUqqz58+HDGjRsH/Kv0SpOp6KR/0XIL3t7ehd+nZJulfc+oqCjS0tKIj4/H29ubVq1aFd4vIKB4+fv77ruPWbNm0aFDh0IZq0LRf7SnnnqKDz74wDClNZxpERVYGHfddRdKKcNbZRoLBW6zgQMHEhERUcwashfa4rATBdbHyfR0xk+YyOh7J5AnMHToUBYsWFCYzLdv3z4yMzMZMmQIH3zwAefOnQModFUVcPbsWU6fPs21117Lm2++SUJCQrHzdevWpX79+oXzF8uXLy+0PqpDTEwMbdq0AYqXb9+3bx9//fUX7du3p1WrViQkJJCXl8ehQ4cKl6Qtjz59+hQrK1/A6dOnady4Md7e3kRHR5OSklLmPS6//HIOHTrEihUrCq2XqmAU/3RpOMMiKmnlLFu2jKVLl2I2mw1plWmK44wQYLe2OJyF2WwmPDycCxcu4OXlxQ23jGTEmAns/yeDkXeNITk5me7duyMiNGrUiC+++IKrr76ahIQEIiIi8PHx4dprr2XWrFmF98zIyOCGG27g/PnziEip65kvXbqU+++/n3PnztG6dWsWL15cJbkL5jhEBB8fH95//33AUu3z/vvvJzQ0FC8vL5YsWYKvry99+vQpXGmwS5cudO/evcI23nrrLUaNGsVbb73FzTffXHj8zjvvZNiwYURERBAeHk6HDh3Kvc9tt91GQkJCsVUPK4tR/NNl4WiLqDQrp2CuyYhWmcYAWBOS5ezN2eG4VeHs+QuyJ/WMbD90Ug6dyJScSoaVakrnuuuuk02bNpV53qi/g8owevRoiY6OFhGRzZs3y5gxY+ze5rJly8TLy0v8/f3F09NT+vbt63AZNI4DHY7rGgT4ehHSuBaNavtyMjOb/f+cJeO8XvGuqpw6dYp27drh5+fHoEGDnC2OXXCG66HkBKu3t7dTMqA1roF2VTkQDw9F07p+1DF5c/ikmaTjmTQI8KFpXROeHlqHV4Z69epVGNGlqTrOmGDVuC5uqThEpNwIJWdTYH38k3Ge4xmWyKvm9f2obdLrUNgCKSMKTFM2Rp/30RgLtxvmmkwm0tPTDd95FFgfbRrVwkMpko5ncvjEOXKrUdNK8y8iQnp6OiaTydmiaDRui9tZHM2bN+fw4cOFFWVdARHh3Pkcjp3PYb+Hop6/NyZvT2eL5bKYTCaaN2/ubDE0GrfF7RSHt7e3y67pkHDoFI99up0Dx85yW0Rznr6uE3X9tPtKo9EYC7dzVbky4ZfWY92Uvkwc2IZV8YcZ+sYWovcec7ZYGo1GUwytOAyGyduTJ67uwOpJfajj58W4xbFM/XQ7p806dFej0RgDrTgMStdL67F2Sl8eiGzD578fYcgbPxC9R1sfRsUIZdA1GkehFYeB8fXyZOrQDqyedAX1/HwYtySW/36yndPntPVhNBy+AptG40SUM8JWlVKPAPdhWa1qJzAOaAp8DDQAfgPuFpHs8u4TEREhcXFxdpbWGGTl5DJv8wHmf3+QS2r5MHtEKFd2CHS2WDUePz+/YlWBCzCZTJjN5lKu0Gicj1IqXkQiqnu9wy0OpVQz4EEgQkS6AJ7ASOAV4A0RCQFOAvc6WjYj4+vlyX+HtOeLSX2o7+/DPUviePSTBG19OBlnrcCm0TgTZ7mqvAA/pZQX4A+kAlcCq/LPLwVudJJshia0eV3WTO7Lg1e25cuEowx+4wc27f7H2WLVWIy+MJRGYw8crjhE5AjwGvAXFoVxGogHTolIwWo/h4FmpV2vlBqvlIpTSsW5UpJfdTCbzaxYsYKZM2eyYsWKQpeIj5cHjw5pz5cP9KFBgA/3LYvj0ZUJnDpXrmdPYyecsQKbRuNMnOGqqg/cAAQDQUAAcE0pHy118kVE3hWRCBGJaNSokf0EdTKxsbEEBwcxd+4E9u2bzty5E2jVqimxsbGFn+nSLN/6GBTCmu1HGfzGFr7V1ofDMfLCUBqNPXD45LhS6lbgahG5N39/NNAbuBVoIiI5SqnewPMiMrScW7nt5LjZbCY4OIjJk0/Rt++/x2NiYN68eiQnp15Ui2nXkdNMXbWDP1PPcGN4ENOHdaZ+gI+DJdcYhdTUVEaOHMnKlSu120xzES43OY7FRdVLKeWvLCVsBwG7gWjglvzPjAG+dLRgZbmGHM3q1asJDs4ppjQA+vaF4OAcPv/884uu6dKsLl8+0IeHrwph3Y5UBr+xhW/+0DkFNRUdHqyxJ86Y4/gVyyT4b1hCcT2Ad4EngEeVUgeAhsAiR8pVGdeQozh48CBt22aWeq5t20wSExNLPefj5cHDV7Xjy8l9aFTblwnL43nwo985mannPuyJkZL/Sq4fvmDBApRShVFfGo0tcEpUlYhMF5EOItJFRO4WkSwRSRSRy0SkrYjcKiJZjpLHbDYzbNgQJk8+xezZZ7n3XmH27LNMnnyKYcOGONzyaNOmDQcOBJR67sCBAFq3bl3u9Z2D6rJmch8euaod63emMviNH/h6l/M7NXfFSKN7HR6scQQ6c5zquYbsyYgRI0hK8iImpvjxmBhISvJixIgRFd7D29ODh64KYc3kvgTWMXH/h/FM+eh3Tmjrw2YYcXSvw4NdDyNZrJVFKw6q7xqyFyaTibVrNzJvXj2mTavFokWKadNqMW9ePdau3VilRYo6BdXhiwf68Ojgdny9K5Uhb/zA17tS7Sh9zcGoo3sdHuxaGMlirSxacWC9a8ge9OzZk6Sko0yZspD27WcwZcpCkpNT6dmzZ5Xv5e3pYQnZndyXJnVN3P/hb0xe8RvpZx3mDXRLjDq6d6fwYFccjVcWI1qslUZEXHbr0aOH2AKz2SyBgfVk5kwkOvrfbeZMJDCwnpjNZpu0c+7cOYmKipIZM2ZIVFSUze5bFbJzcmXud/uk7VNfSfcZG+WrHUcdLoM7MXDgQPHy8pI5c+aIl5eXREZGOlskt2LixIni4eEhEydOdLYoNufo0aMyZswY8fPzE0D8/Pxk7NixkhHnufAAACAASURBVJqaave2gTixou91eudvzWYrxSEism3bNgkMrCe9etWSu+5S0qtXLQkMrCfbtm1ziftXlT2pZ+T6t7dKyyfWyaQP4+V4xnmnyGFvjh49Kv3797fbP+Po0aMlOjpaREQ2b94sY8aMsUs7NQ2TySRYkoCLbSaTydmi2ZRly5aJl5eX+Pv7i6enpyxfvtwh7WrFYQMKLIFnn31WJk2aJM8995xNLYJz5845xKKpKhdycmXe5v0S8tR66TZjo6zbXnXrw94ds7W484jVnXHmaNyROMtitVZx1Pg5jqL5G0lJL/Lbb8tYuPBtQkJCqjQJXR5Gi9oqwMvTgwci27J2Sl+a1/fjgRW/MSkqnuNVmPsw6sSeS/uP3YzqzFMYdf7I1rjqfFSNVhyOyt8wWtRWSdo3qc3nE6/g8avbs2n3MYa8sYW1249aTNIyMHrHbNSIp5pIdQcXNSE6bOnSpQwcOBCAyMhIlixZ4lR5KkuNVhyOsgSMGLVVEi9PDyYNbMtXD/bl0vp+TPnodyZ++BtpGaVbH0bvmGvKiNXIWDu4cNXReE2gRisOR1kCtkjocxQhgbX5bOIVPHF1BzbvOcaQN35gTSnWhyt0zDVhxGpkrB1cuOpovCZQoeJQSrVTSn2nlNqVvx+mlHrG/qLZH0dZArZM6HMEXp4eTBzYhq8e7EuLhgE8+NHv3P9hPMcyirvujN4x6xGrcylrcCEibpubUWOoaPYc+AG4DPi9yLFd1szI22qzNqrKUfkbBRREb82cOdNpeRxV5UJOrrzz/QEJeXq9dH3hG/ni98OSl5cnIjoUVVMxpUUN6Ug354OVUVUVrsehlIoVkZ5Kqd9FpFv+sQQRCbevSqsYW6zHERsby7BhQwgOzqFt20wOHAggKcmLtWs3VitL2105cCyDxz7dQcKhUwzpFMiLN3WhcW1jWUoa4zFmzBjGjRvHwIED8fX1JTv74lppJpMJs9nsBOlqLo5Yj+O4UqoN+SvyKaVuwbLkq1tQmdIeRlmnw5m0bWyZ+3jq2g58vy+Nwa9v4Yvfj5QbeeUKuHNJCyNQdJ4iOTnZ0AEVmipQkUkCtAY2AeeAI0AM0NIaM8dWmy0zx8vCaBnfRmD/Pxly0/9ipOUT6+TeJbHyz2nju9zKoqTbxOgJja6OszKlNcXBnpnjWCyS2/LfBwC1rWnM1pu9FYdRM76NQE5unry35aC0e3q9hD3/jXz+26HCuQ9XoKySFp6entr/bkd0bS/nUjAwAhLEir63XFeViOQBk/PfZ4pIhs1NHgNj1IxvI+DpobivX2s2PNSPto1r8cjK7fxnWRz/nHENN17JUNECcnNzDZnQ6C7oSDfnUpCMCQRZc5/KzHF8q5R6TCl1qVKqQcFmTaOugtEzvo1A60a1+GRCb565riNb9x9n8Os/8Fn8YcPPfZQMFfX09KRv377a/25ndG6GcyiZjAk0suZ+lVEc9wAPAFuA+PzNulAmF8GeeR7uNOFeYH18/XB/2gXW5r+fbue+pca3PkrmoaSnpxs6oVGjqS6lWNh51tyvQsUhIsGlbM6vkeEA7JXxXbSw4r5905k7dwKtWjUlNjbWBlI7j+BLAlg5oTfPXd+JHw9arI9VBrY+SrpNTpw4ARg3oVGjqS4lLWxAWXO/yuRxeAMTgf75h74HForIBWsatgW2yOOoCFvneZjNZoKDg5g8+VSxuZOYGJg3rx7JyamGyySvDknHM3l81XZik08S2b4Rs0eE0aSusb9X0ZyD6Oholi5dql0pGrchMjKSmJgYZs2axeOPPy4iUu2SU5VRHO8D3sDS/EN3A7kicl91G7UVjlAcYOnsV69eTWJiIq1bt2bEiBHV7txXrFjB3LkTmD377EXnpk2rxZQpCxk1apS1IhuCvDxh6c/JvPL1Hrw9PXj2+k7c2qM5Slk12NEYiNTUVEaOHMnKlSu1W8/gFB0YKaX2iUj76t7LqxKf6SkiXYvsb1ZKba9ug0qp9sDKIodaA88By/KPtwKSsYQBn6xuO7bEz8/PZp15TZpw9/BQjOsTTGT7xjz+2Q4eX7WDr3ak8vLNoTStq6OV3IGiJdPnz5/vbHE05bB06dKiu1ZFyFbGVMnNzxwHQCnVGsitboMisldEwsVSsqQHlsTC1cCTwHciEgJ8l7/vdrhCiXVb0+qSAD7+Ty+eH9aJbUknGPL6Fj6JPWTYuY+i1JTM8qp+T6Ovx2Irasrfv6pURnFMBaKVUt8rpX4ANgP/tVH7g4CDIpIC3MC/7rClwI02asNQuFKJdVtQ8I937Ng/jO0TzNcP96NTUB0e/2wHYxbHcvSUsWsUGXWFQ1tT1e9p9PVYbEVN+ftXmcpkCQK+QBjQFfC1JuOwxH0/ACbnvz9V4tzJiq53RMkRe1CTypiUVgk1NzdPlv6UJB2e2SCdn/taPvo1xXBZ52VllptMJmeLZlOs+Z7uXD7E3f/+2LPkiOX+PADUK7JfH5hkTaP59/EBjgOBUgXFAYzHkkcS16JFC1s/T4fhiiXWq0Jl/vFSjmfK7Qt/kpZPrJO73v9FDp8850SJi3P06FEZM2aM+Pn5CSB+fn4yduxYt6thZc33LFk+pE+fPm5T58vV/v4V1Vgred4RiuOimiYUWZuj2g1bXFMbi+zvBZrmv28K7K3oHq5qcdQEKvuPl5ubJ8t+SpKOz1qsjxUGsj7ceURdlOp+z5LrsbRv396t6ny50t+/ojVOSp53hOLYQX7Ybv6+J/CHNY3m3+djYFyR/VeBJ/PfPwnMqegeWnEYm6r84/2VnikjF/5caH0cOpHpQElLp6YU5LP2e7qrW8cV/v4VPfuyzgN5YkX/XZnJ8W+AT5RSg5RSVwIfAV9X4royUUr5A4OBolUCXwYGK6X255972Zo2NM6nKkvLXtrAn6j7LufFG7vwW8pJrn5zKyt+/atgUOEUakpBPmu/p7tOlLvC37+iZ1/WeSwGQfWpSLNgiby6H1gFfAZMADyt0Va22rTFYWyqu7TsX+mZMuo9i/Vx53vGsD405eNKbh13o6JnX9p57O2qKvZhaACEWdOgLTetONyXvLw8+fCXZOn07Abp9OwGWf5zsmHmPjQX4wpuHXelomdf2nlrFUeFrqr8/I06+aXUE4DFSqnXrTJzNJoKUEpx5+Ut+eaR/nRrUZ9nvtjFne//yqET55wtmqYUXMGt465U9Ozt8bepTK2q30Wkm1LqPuBSEZmulNohImFWt24ljqpV5SoU1NQ6ePAgbdq0saqmlpEQET6OPcRLX/1JngjTru3InZe1wMPDPWpe6XpPGkejlIoXkYjqXl+ZyXEvpVRT4DZgXXUb0tgXdy3VDhbr447LWvDNI/3p0bI+z7qZ9aGzkzWuRmUsjluBZ4EYEZmUX6vqVRG52RECloe2OCzUlFLtYLE+VsYe4sV86+PJazpw1+UtXdL68PPzK3UBL5PJhNls7FIsGtfG7haHiHwqImEiMil/P9EISkPzLzVpbXSlFCPzrY+IVg147ss/uOO9X/gr3fWsD3cNY3Vn7FX00NWKKVZ7IQ9NxThqediaVKq9gGb1/Fg6riev3BzK7qNnGPrmFpb+lExenvEr7hZQclU2vVyt8bGXW9HV3JVacdgJR845uGup9opGYUopbu9psT4uC27A9DV/MPK9X0hJL12JGpGqJElqnIe9ysi7anl6rTjsgNlsZtiwIUyefIrZs89y773C7NlnmTz5FMOGDbG55eGupdorOwoLqufHknE9mXNLGH8ePcPVb25l8Y9JLmF96DBW18BebkWXdVeWl+QBdMCyZkatEsevtiZ5xFabURMAo6KipFevWhIdzUVbr161JCoqyuZtulOpdmtqHx09dU7GfPCrtHxindy64CdJSjvrAIk1NQF7Zcc7I+seeyUAKqUeBL4EpgC7lFI3FDk9yx5KzF1wxpxDz549SUo6ypQpC2nffgZTpiwkOTmVnj172rwte2PNKKxpXT8Wj+3Jq7eE8effZ7j6rS18EOMa1ofG2NjLreiK7sry1hz/D9BDRM4qpVoBq5RSrUTkLcD1Yh8dSJs2bVi/PgA4e9G5AwcCuOYa+8w52HJtdGdSMGkcFRWFv78/WVlZVZo0Vkpxa8Sl9AtpxFOrdzJj3W427Eplzi1dCb6k9LkgjaYiCtyKAwcOJCIiouQa3oa7rz0pM49DKbVbRDoV2a+FpdDhbuBKsawZ7lSMmsdx/vx5WrVqWiPyKsrDmkz2yMhIYmJimDVrFk899RT9+vVj8+bNVZZBRPj8tyO8sPYPsnPzmDq0A2OvaIWnC+Z9aDS2wto8jvIsjr+VUuEikgCQb3lcj2W519DqNlgTMJlMrF27kWHDhvDVVzm0bZvJgQMBJCV5sXbtxhqhNGJjYxk2bAjBwZbvv359AI8++gBr126slPvMVqMwpRQ392hO35BLeOrzncxct5sNO1OZc0sYrRvVqtY9NZqaTnkWR3MgR0QuioVUSvURkR/tLVxFGNXiKKBgxJ2YmEjr1q3dpnZURRg1k11EWP37EZ5f8wdZOXlMHdqecX2CtfWhqXHYzeIQkcPlnHO60nAF3GXOoaqUl8n+1VeWTHZnPBelFCO6N6dv20t4avVOXvzqT9bvTOXVW7vSRlsfGk2l0XkcGptj9Ez2xnVMvDc6gjdvD+dgWibXvrWVd7ccJFdHXmk0lUIrDo3NcYVMdqUUN3ZrxreP9KdfSCNmrd/DLe/8xIFjF0fCaTSa4lRacRQs5lSw2VMojWvjSpnsFuujB2+NDCfpeCbXvr2VhT9o60OjKY/yoqoAUEpNAGYAZiwZvOS/On/YqDEkRaPK1q69gLd3JikpPpw548VXX601XICAUoobwpvRu01Dnlm9i9kb9rBh19+8dmsYbRvXdrZ4Go3hqIzF8RjQWURaiUhw/uaWSsNR1WxrAj179uTTT9ewf79gNvswYEA2nTp5cMstwwy7uFTj2iYW3m2xPpLTM7n27Rje+eEgObl5Nm/L1cpoa1wXu/zWKqpJAnwN+FtT18Remy1rVblTrScjcO7cOQkMrCczZxav1TVzJhIYWE/MZrOzRSyXY2fOy/hlsdLyiXUyfF6M7Pv7jE3vP3HiRPHw8JCJEyfa9L4aTUlK+61hZa2qyqwA2A1YDPwKZBVROA9WV1kppeoB7wNdsLi97gH2AiuBVkAycJuInCzvPrbK4zBq3oHRqEom+IoVK5g7dwKzZ1882TxtWi2mTFlo+FBlEWHdjlSe+3IXmVm5PDw4hPH9WuPlWf2YEr3qn8ZRlPdbO3/+vN3XHF8IbAZ+AeKLbNbwFvC1iHQAugJ/Ak8C34lICPBd/r5DqEkr6FWXqq4vYvSQ3MqglGJY1yA2PjKAKzs0Zs7Xe7l5wU/s+yej2vd02TLaGpfDnr+1yiiOHBF5VEQWi8jSgq26DSql6gD9gUUAIpItIqeAG4CC+y4FbqxuG1XFHTo5e1Kd9UUcFZLriLmCRrV9WXBXd+aN6sahk2aufzuG/0UfqNbch171T+Mo7Plbq4ziiFZKjVdKNbVROG5rIA1YrJT6XSn1vlIqAAgUkVSA/NfGVrRRJVwh78CZVMcic1RIrqOW3FRKcX1YEBsf6c/gToG8+s1eRiz4ib1/V936cMUy2hrXxF6/tcoojlHANOAn/nVTWTOx4AV0BxaISDcgkyq4pfKVWJxSKi4tLc0KMf7FlfIOnEF1LLKCkNx58+oxbVotFi1STJtWi3nz6tmk0KOzlty8pJYv/7uzO/8b1Z0jJ80Mm1t160Ov+qdxFPb6rVU4OW5rlFJNgF9EpFX+fj8siqMtMFBEUpVSTYHvRaR9efeyZZHDktVci1azdcXFkGyJNRPd9ir0mJqayrRp0/jkk08wm834+flx++23M3v2bIe5fdLPZvHcmj/4akcqoc3q8uqtYXRoUschbWs01mBtkcNKKQ6lVBegE1D4Hy8iy6rdqFJbgftEZK9S6nmgwE+ULiIvK6WeBBqIyOPl3cfW1XFrajXbijDq+iLLly/nnnvuwcfHh6ysLJYsWcJdd93lcDnW70zl2S92ceb8BR68MoT7B7bB24rIK43G3thdcSilpgMDsSiO9cA1QIyI3FLtRpUKxxKO6wMkAuOwuM0+AVoAfwG3isiJ8u5j9LLq7oQRLTJbLfZkC05kZjN9zR+s3X6UzkF1eO3WrnRsqq0PjTGx50JOBdyCJWT2dxEZp5QKxNLpVxuxLA5VmtCDrLmvxn4UrGleYJFdc43zLTIjLbnZIMCHuXd047rQJjzzxS6Gz4thcmQIkyK19aFxPypjcWwTkcuUUvFAJJAB7BKRzo4QsDy0xaExIicys3l+zR+s2X6UTk0t1kenIG19aIyDtRZHZYZCcfmZ3u9hiaj6DdhW3QY1GnenQYAPb9/RjXfu6sGxjCyGz4vhzU37yM6xfc0rTdnoemD2o0pRVUqpVkAdEdlhL4GqgrY4SqcqpUE09uVkZjbPr/2DLxOO0rFpHV67NYzOQXWdLVaNYNKkSSxcuJAJEyYwf/58Z4tjKBwxOX6viCwqsu8JPCMiL1S3UVuhFcfFGHESWwPf/PE3T6/exalz2TwQ2ZYHItvi46XnPqwlNTWVkSNHsnLlysIwbF0PrGIc4aoapJRan5853gVLzSq9SIEBqU5pEI1jGNq5CZse7c+wrkG89d1+hs+LYdeR084Wy+UprXKArgdmfypUHCIyCkvtqJ1YwnEfFpHH7C2YpuroYo3Gpp6/D2/cHs57oyNIz8zmxv/9yOsb9+q5j2pQXuUAXQ/M/lSoOJRSIcBDwGdYyp3frZTyt7NcmmqgizW6BoM7BfLtI/0Z3jWItzcf0NZHNajIqtD1wOxLZVxVa4FnRWQCMADYDxhzCbcaji7W6DrU8/fh9dvDeX90BCcys7nhfz/yfxv3kpWT62zRXIKKrApdD8y+VGZyvI6InClxLERE9ttVskqgJ8eLY9TSIJryOX3uAjPW7eaz3w7TPrA2r93aldDmOvKqIoxUOcDVsNvkuFLqcQAROaOUurXE6XHVbVBjP+xdkVZjH+r6e/N/t3Xlg7ERnDJnc+P8H3ntG219VIS2KpxHmRaHUuo3Eele8n1p+85CWxylo4s1ui6nzReYuW43q+IP0y6wFq/d2pWw5vWcLZbGzbBbHodS6vf89TKKvS9t31loxaFxV6L3HGPa5ztJO5vFhP6teeiqEHy9PJ0tlsZNsGceh5TxvrR9jUZjQyI7NOabR/pzc/dmzP/+INe/HUPCoVPOFkujAcpXHF2VUmeUUhlAWP77gv1QB8mn0dRY6vp5M+eWriwZ15OzWTmMmP8jL2/Yw/kLeu5D41zKVBwi4ikidUSktoh45b8v2Pd2pJAaTU2mfZ1cvL59mes6NeSdHw5y/dwYfv/rpLPF0tRgdLEcjWEwm82sWLGCmTNnsmLFCqeXSDGKPDNnzuSnHzaT+/Mylt5zGZlZOdy84Cdmb/hTWx8ap+DwNcdtiZ4ctz+OqrRrtOKMRpCnrGJ9fnUa8OAH0Xwce4g2jQJ49daudG9R3yEyadwDh6w5blS04rAvjuo8zWYzwcFBhklcNIo8qampTJs2jU8++QSz2Yyfnx+33347s2fPpkmTJmzZl8aTn+3g7zPn+U+/1jwyuB0mbx15pakYR1TH1dRAHFlp12jFGY0iT0VlNfq3a8Q3j/Tn9p4tWLglkWvf3kp8ip770NgfrTg0peLIztNoxRmNJE9Fxfpqm7yZPSKU5fdeRtaFPG555yde+mq3nvvQ2BWtODSl4sjO02jFGY0kT2XLavQLacTXD/dj1GUteG9rEte+tZX4lBMOk1NTs9CKQ1Mqjuw8R4wYQVKSFzExxY/HxEBSkhcjRoywWVuV4ZprrmHv3jxDyLN06VIGDhwIWIr6LVmypMzP1jZ589JNoUTddzlZOXnc8s7PvLhuN+ZsbX1obIueHNeUSkWVdv/8M5ENGzbYLNrKCFFMReVo1CiL5GQzrVtD585w8KA/KSk+LrME79msHF7e8Ccf/vIXwZcE8OotYUS0auBssTQGwSWjqpRSyUAGkAvkiEiEUqoBsBJohWXBqNtEpNyZvmbNmsmrr76qi/jZibI681deeZMnnnjY5p28s4szloymysqCrVvhl18gNtaflJQj1KvnWgUHfzpwnMc/28GRU2bGXRHM1KHt8fPRkVc1HVdWHBEicrzIsTnACRF5WSn1JFBfRJ4o7z4NGypp166WU+P93Z2Snfk111xDx46tnR6qag9WrFjB3LkTmD377EXnpk2rxZQpCxk1apQTJLOOzKwcXt6wh+W/pNCqoT9zbunKZcHa+qjJuFM47g1Y1jYn//XGii645BLsFiKqseDn58eoUaN45plnGDVqFBs2bDBEqKo9MFI0lS0J8PVi5o1dWPGfy8kV4fZ3f+aFtX9wLjvH2aJpXBRnKQ4BNiql4pVS4/OPBYpIKkD+a+PSLlRKjVdKxSmlCic33KHTchXK61yDg8+yZcsWB0tkO4wUTWUPrmhzCV8/1J/RvVqy+MdkrnlrK78mpjtbLI0L4izF0Sd/IahrgAeUUv0re6GIvCsiESXNLFceEboS5XWue/fCRx8td1nLz2jRXfYgwNeLF27owkf/6YUI3P7uLzy/RlsfmqrhFMUhIkfzX48Bq4HLgH+UUk0B8l+PVeWe7jAidAVGjBhRZqhqSgp06KBc1vKrSUvv9m7TkK8f7sfYK1qx5Kdkrn5zK7+4kPWRmprKgAED+Pvvv50tSo3E4YpDKRWglKpd8B4YAuwC1gBj8j82Bviysvd0pxGh0TGZTNx229288go8/ji8/77l9c034cUXoV27cy5t+fXs2ZOkpKNMmbKQ9u1nMGXKQpKTU90y8MLfx4vnh3fm4/G9ABj57i9M/3IXmVnGtz5mzpxJTEwMM2bMcLYoNRKHR1UppVpjsTIAvIAVIvKSUqoh8AnQAvgLuFVEyk191VFVzmHFihW89dZ4Bg/OJDUVmjaF/v3Bx8e1o49qMueyc5jz9V6W/JTMpQ38mHNzV3q3aVjuNampqYwcOZKVK1cW1s+yN2VVDDaZTJjNZofI4A64XFSViCSKSNf8rbOIvJR/PF1EBolISP5rhfUSTKYgu40IjbIWgxEZMWIEKSnemExw991w1VUWpRETA4mJnmRnZ+vn5mIUWB8rx/fCQynueO8XnqvA+nDGqD8xMZExY8bg5+cHWBTJ2LFjSUpKcpgMGp05XipGyWI2MqU9o337wMND0bat6OdWDRy19kmFcmTn8uo3e1n8UxLN6vkx55YwrmhzSeF5Z4/6ly9fzj333IOPjw9ZWVksWbKEu+66y+7tuhPWWhyIiMtuPXr0EFtz7tw5CQysJzNnItHR/24zZyKBgfXEbDbbvE1X5dy5cxIVFSUzZ86UxYsXS+PG+rlVl23btklgYD3p1auW3HWXkl69aklgYD3Ztm2bw2U5evSo9O/fXzbE7ZcBczZLyyfWydOrd8jZ8xcKz48ZM0b8/PwEED8/Pxk7dqykpqY6RL6BAweKl5eXzJkzR7y8vCQyMtIh7boTQJxY0fcaKQHQEBhlLQZXoGhyoI+PD61b6+dWHRy59kllKHBBrVn0Ohse6s+9fYOJ+vUvhr65hZ8OHK9wnRB7U9mKwRr7oRVHCdw1e9je7NmzB0/PsyxbBps2QXb2v+f0cysfowxW/Pz8UEqxYMEC8vLyWLBgAf6+Xsy6tQer7u+Nj6cHo97/ladX7+T9JcuBstcJsSdVqRissQ9acZTA3bOH7UFsbCzz5v0fGRlw4QJs3AijRsGePZbz+rmVj1EGK+VNPPdo2YD1D/XjP/2CWbHtL/7udh9vfvS1HvXbGFfJT9GKowQ1IXvYlhS4WR599BxvvQX33gtz5sDDD8Mzz0B0tH5uFWGUwUpFLiiTtydPX9eJVff3plmTxrwad55pn+8konffGjPqt3fH7jL5KdZMkDh7s8fkuIixJiqNTlRUlPTqVavYhHjBFhaG1Knjr59bBZjNZsMEZFR24tmcnSMvfbVbgp9cJ1fM/k627DvmMBmdycSJE8XDw0MmTpxo0/uaTCbBUsOv2GYymWzaTgHoyXHbU5Oyh62lPDdLaKjikUem6udWAUYqdVLZiWeTtydPXduRT++/Al9vD+5etI1pn+8g4/yFSrXjKi6ZAkqb/1FKFbr1rMXl8lOs0TrO3mxlcRSElc6YMUOioqJ06GgVKM/i6NWrlkRFRTlbRJehaHizK/0Ozdk5Mmu9xfroPWuTfL+3YuvDXiN3e+GIEORly5aJl5eX+Pv7i6enpyxfvtxm9y4JVlocTu/8rdmCg4Ot7uy1W8o6jORm0TiX31JOyJWvRUvLJ9bJ459ul9Pm7Is+42iXjC2xd8fuyPyUGq04AgI8Kt3Zl2ZV6GQ/26CVr6YAc3aOzF7/pwQ/uU56zdok0Xv+KXbe2cmD1mDvjn306NESHR0tIiKbN2+WMWPG2PT+RanRiqNdu8p19mV1bDNmzNBuFhvhqm4WjX34/a+TMuj/vpeWT6yTqZ8myKlz/1of9hi5F2S721MBldaxO6Jde2Ct4nDpWlXt2ytZuPDf/dIqs5rNZoKDg0pdI/u110xce+15xo/nIhYtUrRvP4NnnnnGjt9AU9MxSn0qe3D+Qi5vfbefhT8cpHFtE7NvDiWyfWMiIyOJiYlh1qxZPPXUU/Tr14/NmzdfdH1Vns2kSZNYuHAhEyZMYP78+fb+ak5v11pcrjquPSktWaq8rNyQEIiN9S31Xjppzf1xdgXk2NhYgoODmDt3Avv2TWfu3Am0atWU2NhYh8phL0zenjxxdQdWT+pDbZMX4xbHMvXT7QS1bFNh5FZln429o53KwlntGgW3UhyldfblhYt26JBFaqrSyX41EGd32karT2VPul5aj3UP9uWByDZ8/vsRDrQbiTTpBJReMqQqj5w8/wAAEPRJREFUz8ZZYawuFz5rY9xGcZTV2VeUlfvEE08ZIn5e4ziM0GkbpT6Vo/D18mTq0A6snnQF9fx8GLcklv9+sp3T5y7O+6jKs3FWwUVnF3p0Ni6tOI4c8aiws6+ohMjUqVN1sl8NwwidtlHqUxXFEa67sOb1WDOlD1OubMsXCUcY8uYPbN7zT7HPVPXZFBRYdHTBRWe1awS8nC2ANTRu3JL27e/hmmtalzlxVpCVO2zYEL766uKFmQqu0Uud1hyc2WkXTPju2LGDhAQf7r47Cx+f4p85cCCAa65x7PxayYW51q8P4NFHH7DLIly+Xp78d0h7hnRqwtRV27lnSRwjujdj+vWdqevvTZs2bVi/PgA4e9G1pT2bgmz3gQMHEhERwdKlS20qb1k4q10j4NJRVVVZAbDgHzYxMZHWrctWNBr3Z8WKFcydO4HZsy/umOy5ZnrJzjkhQTh0CF5+GTp0sHwmJgbmzatHcnKqw36f5UUe2luW7Jw85m3ez/++P0jDAB9m3RRK39Z1adWqqVPkqSlYG1VVYxSHRlPA+fPnHd4xldc5z54NN9wASUm1nLLUrrMUaVF2HTnNY59uZ8/fGdzUrRnDmmcx6uarCQ7OoU2bTP7805ukJC/WrPmGviV9jJoqY63icGlXlUZTHSrrvrQl5c2rdOpk4sKF65ky5SanWMJGmG/p0qwuayb3ZV70AeZHHyDmgA/T5n3Gc/cO4/Bhb1q2zCYkxJtbbhmm17A3AFpxaGokBRWQC9yX5c2T2YKKwsLbt+/qtHm2qs4p2AsfLw8eHdyOIZ0CeeyTBN6IO0vElIk80v9davlkA5nExMCwYUO0u8rJOC2qSinlqZT6XSm1Ln8/WCn1q1Jqv1JqpVLKp6J7aDTWUHTN9FGjRtm1IzLKYk2lYbTFy7o0q8uoRkeolfwpybn9eCpmPvH/9ALcN1zZ1XBmOO5DwJ9F9l8B3hCREOAkcK9TpNIYFmdneluD0TrnohhpPZACUpIO0s17Gc/1foS6vieZ+/szvLP9Mc5m19Zr2BsAp7iqlFLNgeuAl4BHlVIKuBIosNWXAs8DC5whn8Z4ODJc1B44Y16lKjjadVcRBe6ze+skMb33I3yVeCtrDo5kd3pXPE9+wDUuWA7IneqSOSWqSim1CpgN1AYeA8YCv4hI2/zzlwIbRKRLKdeOB8YDtGjRokdKSoqjxNY4CWeGi9oaHRZeOUqLfPvrTDBv/fIQ6XltubZLIC/eFEaDANfwaJcc+BQdNDhj4ONyUVVKqeuBYyISr5QaWHC4lI+WqtFE5F3gXbCE49pFSI2hKC8i6auvLP5uV0ngLJhX0ZRP6RZaGkeSZzLmlVV8/ucxtiX/wIs3duHqLk2dLW65FC1x8+9v+GylJvqNaqU4Y46jDzBcKZUMfIzFRfUmUE8pVaDImgNHnSCbxoAYIVxU43gK3GfFygElHeaV0QNZM7kvTeqauP/D35i84jfSz2Y5W9wyqW6JG2cX4iwPh1scIjINmAaQb3E8JiJ3KqU+BW7BokzGAF86WjaNMTFKuKjG8ZRloXVsWofVk/qw8IeDvPXdfn4+mM7MG7twbajxrI/qDHyssVIcgZGKHD6BZaL8ANAQWORkeTQGwcgRSRrn4e3pweQrQ1g7pS9B9fyYFPUbD0QZz/qoTii2EQpxlodTEwBF5Hvg+/z3icBlzpRHY0yMHpGkcS4dmtRh9aQrWLglkbc27efnxHRm3tCF68KMYX2MGDGCRx99gJgYLgruKGvgY3T3rM4c17gERgsX1RgLL08PHohsy1UdA5m6ajsPrPiNr3Y2YcYNXbikVumrfDqK6gx8jO6e1UUONRqNW5GTm8e7WxN589v9BPh6MuOGLlwf1hRLupjzqEootr0LcerquFpxaDSaUtj/TwaPfbqd7YdPc3XnJsy8sQuNajvX+qgK9sz90IpDKw63xagx7BrXISc3j/e2JvHGt/sI8PXkhRu6MMwA1kdlsVfCqFYcWnG4JUbLtNW4Nvv/yeCxVTvYfugUQzsHMvPGLjSuXXMHIVpxaMXhdrhTiRGNccjJzWNRTBL/9+0+/H08eWF4Z4Z3DXIZ68OWWKs4jJTHodEAxo9h17gmXp4eTBjQhvUP9qVVwwAe+jiBCcvjOZbhOlWWjYJWHBrDYfQYdo1r07ZxbT6beAVPXduB7/elMfj1LXzx+xFc2fviaLTi0BgOIy96pHEPPD0U4/u3Yf2D/WjTKICHVybwn2XxHDujrY/KoBWHxnDoEiMaR9G2cS0+vf8Knr62I1v3pzH4jS2s/v2wtj4qQE+OawyJjqrSOJqDaWd5fNUO4lNOclXHxrx0UyiBddwzCENHVWnF4bboRY80jiY3T1j8YxKvfrMXXy8Ppg/rzIjuzdwu8korDq04NBqNjUk6nsnUT7cTl3KSKzs0ZtZNoTSp6z6DFh2Oq9FoNDYm+JIAVk7ozXPXd+Kng8cZ/MYPfBp3SM995KMVh0aj0ZSCp4finr7BbHioPx2a1Gbqqh3csySWv0/ryCutODQajaYcgi8JYOX43kwf1omfE9MZ/MYPfFLDrQ+tODQajaYCPDwU4/oE8/VD/enYtA6Pr9rB2MWxpJ42O1s0p6AVh0aj0VSSVpcE8PF/evH8sE5sSzrBkNe38ElszbM+tOLQaDSaKuDhoRjbJ5ivH+5Hp6A6PP7ZDsYsjuXoqZpjfehwXI1GU4heA6Vq5OUJH/6awuz1e/D0UDxzXUdu73mp4fM+dB6HVhwajU3Q2frV56/0czz+2XZ+STxBv5BLePnmMJrV83O2WGWiFYdWHBqN1eg1UKwnL0+I+jWF2Rv24KEUT1/XkZEGtT5cLgHw/9u7+xipqjOO49/fLpa3tcBiNVugBYViG8qu21aW0mxaMBWtsaYpjdRg06CExrZi2jS1tiSmTV+iidjEYilQsKmIEm3NJgUJfYEaCwFl2V0R0UBhlYotiIJiZH36xzkTptOZXQbGvXt2n08ymbl37sz8uHNnD/fce58jaYikbZJaJXVIujPOnyBpq6S9ktZKel9vZ3NuoPIxUM5dVZWYN308GxY18/ExI7j90TZuXLmNzqNvZh2t4rI4OP42MNPM6oEGYLakJuAXwD1mNgk4CszPIJtzA5KPgVI542qH8fubpvHj66aw459Hmb1kCw9uPdCvzrzq9YbDguNx8rx4M2AmsC7OXw1c19vZnBuofAyUyqqqEvOaPsyGRc1MHTuCHzzWxrwV2zh4pH/sfWRyjENSNbADmAjcB9wF/MPMJsbnxwF/MrMpRV67AFgQJ6cA7b0S+r1xAfDvrEOcA8+fnUpnl0R9XR3VNTWnZx4/DocO0WVGK+E/eJWS8rqH9PNPNrPzz/bFgyqZ5EyZWRfQIGkk8Bjw0WKLlXjtMmAZgKTt53KAJ2ueP1sp5085O3j+rEk6p7OKMr0A0MxeA/4KNAEjJeUasrHAy1nlcs45V1oWZ1V9IO5pIGkocAWwG/gL8OW42NeAP/Z2Nueccz3LoquqDlgdj3NUAQ+bWYukZ4GHJP0EeAZYcQbvtew9zNkbPH+2Us6fcnbw/Fk7p/xJXwDonHOu93mRQ+ecc2XxhsM551xZkmk4+kOpEknVkp6R1BKnU8q+X1KbpJ25U/kk1UraGPNvlDQq65ylSBopaZ2k5yTtljQ9lfySJsf1nru9LmlRKvkBJN0Wf7ftktbE33MS27+kW2PuDkmL4rw+ve4lrZR0WFJ73ryimRX8UtILknZJauzp/ZNpOOgfpUpuJZxBlpNSdoDPmVlD3vnr3wc2xfyb4nRfdS+w3swuBeoJ30MS+c1sT1zvDcAngDcJ1z8lkV/SGODbwCfjRb3VwPUksP1LmgLcDFxO2G6ukTSJvr/uVwGzC+aVynwVMCneFgBLe3x3M0vuBgwDngamEa7eHBTnTwc2ZJ2vROax8cuaCbQASiV7zLcfuKBg3h6gLj6uA/ZknbNE9vcD+4gng6SWvyDz54EnU8oPjAEOArWEMzlbgCtT2P6BOcDyvOkfAd9LYd0D44H2vOmimYFfA3OLLVfqltIeR66rZydwGNgIvAi8Zman4iKdhI20L1pC2ODejdOjSSc7hCv5n5C0I5Z9AbjIzA4BxPsLM0vXvYuBV4Hfxq7C5ZKGk07+fNcDa+LjJPKb2UvA3cAB4BBwjFByKIXtvx1oljRa0jDgamAciaz7AqUy5xr2nB6/i6QaDjPrsrC7Ppaw63jGpUqyJOka4LCZ7cifXWTRPpc9zwwzayTs1t4iqTnrQGUYBDQCS83sMuAEfa9roUfxGMC1wCNZZylH7Ev/IjAB+CAwnLAdFepz27+Z7SZ0qW0E1gOtwKluX5Sesv8WJdVw5Fh6pUpmANdK2g88ROiuWkIa2QEws5fj/WFC//rlwCuS6gDi/eHsEnarE+g0s61xeh2hIUklf85VwNNm9kqcTiX/FcA+M3vVzN4BHgU+TSLbv5mtMLNGM2sGjgB7SWfd5yuVuZOwF5XT43eRTMORcqkSM7vdzMaa2XhCV8OfzewGEsgOIGm4pPNzjwn97O3A44Tc0Ifzm9m/gIOSJsdZs4BnSSR/nrmc7qaCdPIfAJokDZMkTq//VLb/C+P9h4AvEb6DVNZ9vlKZHwdujGdXNQHHcl1aJWV9AKeMAz1TCaVIdhH+aC2O8y8GtgEvEHbhB2edtYd/x2eBlpSyx5yt8dYB3BHnjyYc8N8b72uzztrNv6EB2B63nz8AoxLLPwz4DzAib15K+e8Enou/3d8BgxPa/rcQGrpWYFYK657QuB0C3iHsUcwvlZnQVXUf4ZhxG+Hst27f30uOOOecK0syXVXOOef6Bm84nHPOlcUbDuecc2XxhsM551xZvOFwzjlXFm84XL8kqaugomyvXSlerDKpc/2Jn47r+iVJx82sJqPPbgaOAw9YqAbbG59ZbWZdvfFZzvkehxswJI2QtCd3BXkcF+Lm+HippO3KG+slzt8v6aeSnorPN0raIOlFSQuLfY6ZbSaUpuguy5w4xkOrpM1xXrWkuxXGPdkl6Vtx/qxYnLEt7s0Mzsu2WNLfgTmSLpG0Phai3CLp0kqsN+cKDep5EeeSNDRWUs75mZmtlfRNYJWke4FRZvab+PwdZnZEUjWwSdJUM9sVnztoZtMl3UMY52AGMIRwFf39Z5lvMXClmb2UK6VDGAthAnCZmZ2KA+8MiZ85y8yel/QA8A1CrTOAk2b2GQBJm4CFZrZX0jTgV4S6aM5VlDccrr96y0Il5f9hZhslzSGUWKjPe+orsVz8IMJYBR8jlCeBUMsHQjmGGjN7A3hD0klJIy0U3SzXk4QG7GFC0T8I9dfut1hqPDZk9YQCgc/HZVYDt3C64VgLIKmGUDjwkVAOCghlPZyrOG843IAiqYpQjv8twsBCnZImAN8FPmVmRyWtIuxR5Lwd79/Ne5ybPqvfkJktjHsFXwB2Smog1AwqPOhYrOR1vhPxvoowvsX/NZbOVZof43ADzW2EqspzgZWSziOMEHgCOCbpIoqPFVFRki4xs61mtpgwEt444AlgYa7UuKRaQmHA8ZImxpfOA/5W+H5m9jqwL+5N5caRri9czrlK8IbD9VdDC07H/bmkjwA3Ad8xsy3AZuCHZtZKqLzcAawkdCOdNUlrgKeAyZI6JRUbS/uueLC7PeZoBZYTSpDvktQKfNXMTgJfJ3RBtRH2ckodV7kBmB9f20EYPMm5ivPTcZ1zzpXF9zicc86VxRsO55xzZfGGwznnXFm84XDOOVcWbzicc86VxRsO55xzZfGGwznnXFn+C+Ye1hsAxktxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotDecisionBoundary(plotData, theta, X, y)\n",
    "plt.xlabel(\"Exam 1 score\")\n",
    "plt.ylabel(\"Exam 2 score\")\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a decision boundary we can create a function to predict whether a given student (a single sample of two exam scores) will be admitted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(theta, X):\n",
    "    \"\"\"\n",
    "    Predict whether the label is 0 or 1 using learned logistic regression.\n",
    "    Computes the predictions for X using a threshold at 0.5 \n",
    "    (i.e., if sigmoid(theta.T*x) >= 0.5, predict 1)\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    theta : array_like\n",
    "        Parameters for logistic regression. A vecotor of shape (n+1, ).\n",
    "    \n",
    "    X : array_like\n",
    "        The data to use for computing predictions. The rows is the number \n",
    "        of points to compute predictions, and columns is the number of\n",
    "        features.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    p : array_like\n",
    "        Predictions and 0 or 1 for each row in X.\n",
    "    \"\"\"\n",
    "    # Number of training samples\n",
    "    m = X.shape[0]\n",
    "    \n",
    "    # initialize p\n",
    "    p = np.zeros(m)\n",
    "    \n",
    "    temp = sigmoid(X.dot(theta))\n",
    "    for i in range(m):\n",
    "        if temp[i] >= 0.5:\n",
    "            p[i] = 1\n",
    "            \n",
    "    return p\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For a student with scores 45 and 85,we predict an admission probability of 0.776\n",
      "Train Accuracy: 89.00 %\n"
     ]
    }
   ],
   "source": [
    "#  Predict probability for a student with score 45 on exam 1 \n",
    "#  and score 85 on exam 2 \n",
    "prob = sigmoid(np.dot([1, 45, 85], theta))\n",
    "print('For a student with scores 45 and 85,'\n",
    "      'we predict an admission probability of {:.3f}'.format(prob))\n",
    "\n",
    "# Compute accuracy on our training set\n",
    "p = predict(theta, X)\n",
    "print('Train Accuracy: {:.2f} %'.format(np.mean(p == y) * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>2 Regularized Logistic Regression</h3>\n",
    "\n",
    "In this part of the exercise, we will implement regularized logistic regression to predict whether microchips from a fabrication plant pass quality assurance. \n",
    "\n",
    "First we visualize the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "# The first two columns contains the exam scores and the third column\n",
    "# contains the label.\n",
    "data = np.loadtxt(os.path.join('Data', 'ex2data2.txt'), delimiter=',')\n",
    "X, y = data[:, 0:2], data[:, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO29e5gU5bXo/Vsw4syIZhSRgKgzgDGSmONlcHNy0ZmQgJigkW0U+UgGxUdEJTmbLzlx1J24nZM9yvlOEoU4G48iSGYEdYdEokm8gObjJIbBBG+41WEgCZlRUWO2SDOCrPNH1UDRdPf0pW7dvX7PU093vfVW1erq7lr1vusmqophGIZh5MqgqAUwDMMwihNTIIZhGEZemAIxDMMw8sIUiGEYhpEXpkAMwzCMvKiIWoAwOfbYY7W2tjZqMQzDMIqKZ5999i1VHZ7cXlYKpLa2lo0bN0YthmEYRlEhIn9K1W5TWIZhGEZemAIxDMMw8sIUiGEYhpEXZWUDMQzDKIQ9e/awfft2du/eHbUogVBZWcno0aM57LDDsupvCsQoaxKJBKtXr2bLli2MHTuW6dOnU1lZGbVYRkzZvn07Rx55JLW1tYhI1OL4iqry9ttvs337durq6rLax6awjLKls7OTurpRLFo0l1df/R6LFs2ltnYknZ2dUYtmxJTdu3czbNiwklMeACLCsGHDchpd2QjEKEsSiQTTpk3muuve5bOf7W/dyfr1MG3aZLZt67WRiJGSUlQe/eT62WwEYpQlq1evpq5ur0d5OHz2s1BXt5ef/vSn0QgWIIlEgo6ODlpaWujo6CjZeXwjPEyBGGXJli1bGDfu/ZTbxo17n+7u7pAlChabritvFi9ezLhx4xAR3nrrLd+OawrEKEvGjh1LV9cRKbd1dR3BmDFjQpYoOLzTda2tO5kzR2lt3cl1173LtGmTbSQSML29vZx77rm8/vrrkcnwmc98hieeeIKTTjrJ1+OaAjHKkunTp7N1awXr1x/cvn49bN1awfTp06MRLADKcbouTrS0tLB+/XpuueWWgo/1z//8z9x+++3712+88UbuuOOOAfc744wzCCIPoBnRjbKksrKSNWseY9q0yTzyyF7GjXufrq4j2Lq1gjVrHispA3q5TdfFhaqqqoNGd21tbbS1tVFZWUkikcjrmHPmzGH69Ol885vfZN++faxcuZK1a9dy+umnp+zf0dHB+PHj8zpXNpgCMcqWCRMmsHVrD6tXr6a7u5upU8eUZBzI2LFjefTRI4Cdh2zr6jqCqVNLZ7ouTnR3d9Pc3MwDDzxAIpGgqqqKSy+9lNbW1ryPWVtby7Bhw/jjH//IG2+8wRlnnMFJJ53Epk2bfJQ8e0yBGGVNVVUVM2fOjFqMQJk+fToLFlzL+vUcNI1VitN1cWLkyJFMmjSJ9vZ2qqur6evrY9KkSXz0ox8t6LhXXnkly5Yt4/XXX+eKK67gvffe43Of+1zKvjYCMQyjIMppui5uLF26FICbb76ZG264gaVLlzJr1qyCjnnRRRfx3e9+lz179tDR0cHgwYNtBGIYRnCUy3Rd3DjxxBN5/PHHaWhooL6+nuXLlxd8zCFDhtDY2EhNTQ2DBw/Oap877riDhQsX8vrrr/OpT32K888/n7vvvrtgWURVCz5IsVBfX69WUOoAlgcqd+yalTcvv/wyp556aqQy7Nu3jzPPPJMHH3yQk08+2ffjp/qMIvKsqtYn943UjVdElorImyLyYprtIiJ3iEiXiDwvImd6tjWJyGvu0hSe1KVBuQaWFRKNXa7XzIgPmzdvZty4cUyaNCkQ5ZErUU9hLQMWA/el2T4VONld/gFoA/5BRI4BvgfUAwo8KyIPq+rfApe4BCjXPFCdnZ1MmzaZujrHDvDoo0ewYMG1rFnzGBMmTMi4b7leMyNejB8/PlZu15GOQFT1N8A7GbpcCNynDs8ANSIyEpgCPK6q77hK43HgvOAlLg3KMbCs0GjscrxmhjEQcY9EPx74i2d9u9uWrv0QROQqEdkoIht37NgRmKDFRDkGlhWqAMrxmhnGQMRdgaTKLawZ2g9tVL1LVetVtX748OG+ClesxDUPVJA5gwpVAHG9ZoYRJXFXINuBEzzro4GeDO1GFsQ1D5SfOYOSKVQBxPWaGUaURO7GKyK1wC9U9ZMptn0JuA44H8eIfoeqnu0a0Z8F+r2y/gCcpaqZ7Cnmxush2aDsDSwbyKDsN8k5g/opJGdQMrt376a2dmSSEdxRAIsX12RlBI/ympn7cDzI1Y03Lt/b1q1bmTFjBu+88w5nnnkmK1asYMiQISn75uLGi6pGtgD3A73AHpxRxRzgauBqd7sAPwa2AC8A9Z59rwC63OXybM531llnqXGAXbt2aXt7u7a0tGh7e7smEom8j3HLLbfkfYyenh5tamrSqqoqBbSqqkpnz56tvb29OR8rExs2bNARI2p04sShOmuW6MSJQ3XEiBrdsGFD1sfI9pr5cV38lNvwh82bN2fdN07f21e/+lW9//77VVV17ty5euedd6btm+ozAhs11T08VWOpLqZA/MXPP8h9992nFRUVWl1drYMHD9YVK1YEILE/SnMg/Lwuu3bt0hEjarSlBV237sDS0oKOGFETiPxGerJVIEF9bzfddJP+6Ec/2r9+ww036O23355xn3379umwYcN0z549qqr629/+VidPnpy2fy4KJOo4ECNk/BpS+x0XEUTOoFQEnTzR7+uSyXvskUcc77FSTwZZjAT1veWTzv24446jpqaGigrndj969Gj++te/5nzuVJgCKSMKCaRLxu8/SBA5g6LA7+ti7sPZExd7AwT3veWTzj1V+IJIKkfW3DEFUib4/WTs9x/EqzAaGxtpbGzMaf+44Pd1sVoe2eHnw5EfBPm95ZrO/dRTT+Xdd99l7969VFRUsH37dkaNGpX3+b3E3Y3X8Am/I6ktLiI1fl8Xcx8emDjWfA/ye7vooov41a9+RWdnJ1OmTOHII49k06ZNKZfx48cjIjQ2NvLQQw8BzsPahRdeWMjH248pkDLB7ydju7Glxu/r0l/LY/HiGpqbh3LPPUJz81AWL66xWh4ucUwzE+T31p/O/ZJLLsk6nfttt93GD37wA8aNG8fbb7/NnDlz8j6/F5vCKhP8HlJbkaLUBHFdrJZHZuJqJwrqe9u3bx/PPPMMDz74YNb7jBkzhg0bNhR03lSYAikTgihraje21ARxXcqh9G6+xNlO5Pf3tnnzZr785S9z0UUXxSKde+SR6GFS7pHocYo+Nwy/8CPLQLbEoaBU0OQSiW4jkDLCRgxGocTJVbafsKdTVdU3N9i4keuAwkYghmFkRdxHsP3Krbu7mzFjgnk42rp1K0ceeSTDhg0rOSWiqrz99tu899571NXVHbQt3QjEFIhhpKC3t5cZM2awatUqPvrRj0YtTuQkEgnq6kaFMk0UZ/bs2cP27dsjcQ0Og8rKSkaPHs1hhx12ULtNYRlGDnhTy995551RixM5llLF4bDDDjvk6bycsTgQo6QotChVVVUVIkJbWxv79u2jra0NEaGqqspnSYuLuLrKGtFiCsQoKQotStXd3U1TU9N+hVFVVcXs2bPZunWrbzImEgk6OjpoaWmho6Oj4OkQv4+XCss8YKTCFIgRGX6WsPVr5DBy5EgmTZrEnj17qK6u5oMPPmDSpEm+2UE6OzupqxvFokVzefXV77Fo0Vxqa0fS2dkZi+OlwzIPGClJleO9VBerBxIv5s2bp4MGDdJ58+YVfCw/i1I1NDRoRUWFLly4UCsqKrSxsbFg+VT9rxERdq0Qv+uc+FVwywge0tQDsRGIETpB2Bn8HDn0p5b/9re/zWOPPcaJJ56Yt1xe/M7ZFHYOqP44ovnzl3DKKbcwf/4Stm3rzdmFN6xRkxE8kXphich5wO3AYOBuVb01afsPgf683tXAcapa4277EKfMLcCfVfWCcKQ2CqW7u5vm5mYeeOABEokEVVVVXHrppbS2thZ0XL+KUgWVWt5vQ3QUhu1CU3P4XVbAiJbIRiAiMhin3vlUYDxwmYiM9/ZR1X9S1dNV9XRgEeB9pEr0bzPlUVwEZWcIauTgF34boovRsO33qCkMBwIjPVFOYZ0NdKlqt6p+AKwEMiWpvwy4PxTJjMDxjhZEZP96ISxfvpyGhgbAGTksW7as4GP6id+G6GI0bPs5arKpsOiJcgrreOAvnvXtwD+k6igiJwF1wFpPc6WIbAT2Areq6s/S7HsVcBUQuyfScqZUStjmgt85m4oxpb5fmXNtKiweRJbKRES+CkxR1Svd9a8BZ6vq/BR9vwOM9m4TkVGq2iMiY3AUyyRV3ZLpnOWSyiTIhHdxTKZXbPidsymMHFB+4Vfm3I6ODhYtmktr66GKqLl5KPPnLymLyPiwiGMqk+3ACZ710UBPmr4zgGu9Dara4752i8hTwBlARgVSDgRZGzpudae9xC13VSZ5/K4RUUy1QvwaNVlkfDyIUoF0AieLSB3wVxwlcci/QEROAY4GfudpOxrYpap9InIs8BlgYShS+0BQT/FBDuvjMmWQ7sYct9xVcZNnIMJUwH6UFYhzEamyIlVwSFgLcD7wKs7I4Ua37RbgAk+fm3FsHN79Po3jwvuc+zonm/PFIZDQz2CsZNrb23XixKEHBZX1LxMnDtX29vZYHjsXkoMPKysrFThkqaysHPBYPT09es455+QVbJiOQuSJEj+DOsMgkUiEGkRZ7pAmkDDSOBBVfRR4NKntu0nrN6fY77fAaYEKFwBBP8UHOayPesqgqqrqIBfNtrY22traOPzww2lqasorpiSIUUJQMS5Bke66VlZWkkgk8jpmGKOZODkQlLNd0CLRQyToyOEg4wKijjlIl+Rw27ZtOceUBJlxN+hcWn4TRPLIQhNaZotfkfGFUO6uxKZAQiTop/gg4wKijjnIdGPONaYk6Iy7QcS4BIWfCi+KVPj9DgQ33XQTM2fOTPnkH1SwoXdGobV1J3PmKK2tO7nuuneZNm1yWQQ1mgIJkaCf4vuH9YsX19DcPJR77hGam4eyeHFNwcP6II+dLeluzCeeeCIrV67kF7/4BStXrhww3ifoUULcI+KT8UvhhZEKP1eCHCGEnYssjpgCCZEwnuLzGdZnm1Y96imDdDfm5cuX8+STT7J+/XqefPLJrCLQgxwlxD0iPhm/FF7cpu+CHiFEbReMBaks66W6lLoXVr4UmweOl3y9nr7+9a/runXrVFV17dq12tTUFLywMSHIVOpBpcLPh6A9B+PimRgGWDr3eBD1U7yXUijfmu+0SbGNEvwiaKNvnKbvitnmWCxE6sZbrsQlcrjYXE5T0T9t0t7eTnV1NX19fbH2eoqSMIJBg0qFn0w2rrNBBxvGyZU4KmwEUoT4VQo2bnPW+VJMXk9RUipG32xHUXG1OZYSNgIpQvwMgPOrCFOUlGNm33woBaNvLqOosEYIcZlRiAIbgRQRQdgs4jRnnS+lZM/wa3SZikxu5K+9Vk17e3sg5/WTXEdR5T5CCBpTIEVEEH72pXTzLQWCjOLONKWzefMeXnnllcCjxwsln1FUNsGGRn6YAikiSsVmYRxKGB5xqYJB58+HlhZ4770PUNXYe+JFnVLHOBhTICHhVzoFMxjnRqYpoSCni3IlrCju5Cmdq69ezCWXzIpV9HgmzHU2ZqQKDinVJapAQj+DB8s5AC4fMgVJxi2A8r777tOKigqtrq7WwYMH64oVK1Q1mLTz2Zw3rsQxGLfUIU0gYeQ39TCXKBTIrl27rG5BBGSKUI9rzY50UdxBK7o4RY9nS380fUtLi+/R9MahpFMgNoUVMKXie19sZJoSimPSPzjUI+7pp5/O2i5SyHRcMXrimWE8HpgCCZhS8L0vRjI5HMTVGSHZI2779u1ZK7pCvLfMEy89QaWCLxXSKhAR+YSIrBeRrSJyp4h8xLPtd+n2ywUROU9EXhGRLhG5PsX22SKyQ0Q2ucuVnm1NIvKauzT5IU8QBOU1ErcfdpwM0v1kcjgoBmeEbBRdKeQziyvlXiwqK1LNazlTXvz/wJeBY4HrcWqP17nb/phuv2wXYDBOLfQxwBCc+ubjk/rMBhan2PcYoNt9Pdp9f/RA54zCBhJE7eY4GhHjZpBWzexwUCzOCAPZJ3p6erSpqUmrqqoU0KqqKp09e3ZgBvdywWyXB0OuRnRgU9L6F4DXgAnAH9Ltl+0C/Ffg1571ZqA5qU86BXIZsMSzvgS4bKBzloIXVtx+2HE1SJcK2Si6YvOiGoigvc6yoZxStWdDOgWSyQYySESO8oxUngC+CnQAfljZjgf+4lnf7rYl848i8ryIPCQiJ+S4LyJylYhsFJGNO3bs8EHs3PEznULcjPJ+GKTjOP0VF7KxTxTDdFwuhFVTPRNmu8yOTArkfwKf8Dao6ibgi8AaH84tKdo0aX0NUKuqnwKeAPqz5GWzr9Ooepeq1qtq/fDhw/MWtlD88hop9Ift983aD4N0HG4YxUwxelGlIk72HIt4z460CkRVV6jqIcZyVd2mqpf7cO7twAme9dFAT9K53lbVPnf1fwNnZbtvqVLoDzuIm3W+T8BxumEUM6XiRRUn92qLeM+SVPNaYSw4qeS7gToOGNE/kdRnpOf9RcAzesCIvhXHgH60+/6Ygc4Zh5K2hZKvUT5IW0W+BmkzABvJRGXPSWV3iaOzSlQQx0h04HzgVRxvrBvdtluAC9z3rcBLrnJZB3zcs+8VQJe7XJ7N+UpBgajm98OO68261AzARmFEFRWfzovQIt4d8lYgwMRs2ophKRUFoprfDzuON+tiTKNhBEfY7tXmRZgdhSiQQ1x2gWcH2i+OSykpkHyI4826WOIxjNIkriPzuJFOgaQtaSsiZ+PEagwXkW94Nh0FHDagccWIHXEs/eqVobGxkcbGxgilMcqNfi/C9vZ2qqur6evri0Vam2IhkxvvEThR6BXAcM/yAU48iFFklIq3TtRY3EppUWpxNGGSdgSiquuAdSJyr6p2A4iIANWqmjoQwTDKAK8r9J133hm1OEaBxHFkXiyIM72VoYPIfcB1wF5gI86o5FZV/UHw4vlLfX29bty4MWoxjCKlqqoqZdLKyspKEolEBBIZRjiIyLOqWp/cnk0699NU9T+BrwCP4QTtzfZXPMOIP3EKdDNKh7hl1s6FbBTIEBGpAC4EfqaqHwD7ghXLMNITlQ0irnVEjOKl2FPGZ6NA7gb+jBPx/bSInAjsDFQqw8hAlLmzzOBq+EUikWDatMlcd927tLbuZM4cpbV1J9dd9y7Tpk0uipHIgApEVX+oqqNUdbLrD7wd+HzwohnGwcQhd1apJC40oidumbXzYUAFIiLDRWSJiPzCbfo4MDNYsQwv5jbqEAcbhLlCG35RCinjs5nCWgY8zYHst68B/29QAsWNOBi4LN25g9kgjFKiFFLGZ6NAjlPVDlzDuaruAT4MVKqYELWBKw5TNrkS9GjJbBBGqVAKKePTBhJ6eF9EjsEt2CQiE4D3ApUqBngNXAfmKHeyfj1MmzaZbdt6UxaFSiQSrF69mi1btjB27FimT5+ed/Go7u5umpubeeCBB0gkElRVVXHppZfS2tqa/wcLmKCD7CzoyygVKisrWbPmMaZNm8wjj+xl3Lj36eo6gq1bK1iz5rG87xthkk0gYT1wO051wudwSsderE51wqIil0DCjo4OFi2aS2vroQ5nzc1DmT9/CTNnHmwK6uzsZNq0ydTVHfpjyKd8LcCKFSu44oorGDJkCH19fSxbtoxZs2bldawgsSA7w8iP/ofO7u5uxowZU9BDZ1CkCyTMlExxoqo+o6obRaQROBWnlOxmNxakpMnVwJXviGUgvFM2N9xwA0uXLo2lAinG0ZJh+EGhsw795a6LkUw2kP3zD6r6gao+p6qbykF5QO4GrqBc8orFbdQM3EY5ErWdNGqyMaKXJbkauIJyySsmt1EzcBtBEUdX9lIIBCyUTApkjIg8nG7x4+Qicp6IvCIiXSJyfYrtC0Rks4g8LyJPishJnm0fisgmd/FFHi/9Bq7Fi2tobh7KPfcIzc1DWby4JqWBqxRc8gqlWEZLRvERR1f2UggELJS0RnQReQ24Mt2Oqvp0QScWGYxTD/2LONHtncBlqrrZ06cR+L2q7hKReUCDql7qbtupqkNzOWc+2XizNXDt3r2b2tqRSTYQZ8SyeHFN3jYQwyhn4uyc0dLSwquvfo85cw69h95zj3DKKbdw0003RSCZ/+RsRAfeK1RJDMDZQJen1shKnISN+xWIW5Okn2eA0K3H2Rq4SsElzzDiRpydM8aOHcujjx5BqtSAXV1HMHVq6c86ZJrC2hbwuY8H/uJZ3+62pWMO8EvPeqWIbBSRZ0TkK0EImCsTJkxg69Ye5s9fwimn3ML8+UvYtq03bxdewyh34uycUQqBgAWTqlB6GAtOWdy7PetfAxal6TsLZwRyuKdtlPs6BkfZjU2z71U4hbA2nnjiibnUkTdUtaenR8855xzt7e2NWhSjTGloaNCKigpduHChVlRUaGNjY9Qi7WfDhg06YkSNTpw4VGfNEp04caiOGFGjGzZsiFo0XwE2aor7a5ReWNs5kF8LnEJVPcmdROQLwI3ABara19+uqj3uazfwFHBGqpOo6l2qWq+q9cOHD/dP+iIjXy+WOBovjfIizs4Z5T7rMGAkemAndopUvQpMAv6KY0SfqaovefqcATwEnKeqr3najwZ2qWqfiBwL/A64UD0G+FSUc0nba665hiVLljB37tysUozE2XhpGEa4FFLSFhGZLiI/EJH/JSIX+SGQqu7FqbX+a+Bl4AFVfUlEbhGRC9xu/xMYCjyY5K57KrBRRJ4D1uHUaM+oPMqVfBMyxiF1umEY8SabeiB3AlcDLwAvAnNF5Md+nFxVH1XVj6nqWFX9vtv2XVV92H3/BVUdoaqnu8sFbvtvVfU0Vf0v7us9fshTiuSrCOJsvDQMIx5kMwI5F5iiqveq6r3A+UBDoFIZvlGIIohrZHkco5INoxzJRoG8AnitVicAzwcjjhEE+SqCuBovzbBvJGMPFdGQTTr3p4EJwAa3aQKO0XoXQP+0UjFQrkb0pqYmLr/8choaGli3bh3Lly+PdU6tdJhh30hHrk4iRm6kM6Jno0DOzbRdg41W95VyVSClQm9vb9qoZLPNlCf2UBEOeXthqerTmZZgxDWMQzHDvpGMeQtGS1oFIiLr3df3ROQ/Pct7IvKf4YloGAeIq2HfiAZ7qIiWtApEVT/rvh6pqkd5liNV9ajwRDSMA8TVsG9Ehz1UpCeRSNDR0UFLSwsdHR2+1yjJlI13P27q9RHe/qr6Z18lMYwsWL58+f73jY2NNDY2RiiNEQf6HyoaGhqor68/6DdSznR2djJt2mTq6pzs4I8+egQLFlzLmjWP+ZZqJRsj+nzge8AbwD63WVX1U75IECJmRDcMo1gopNZ6IpGgrm6Ub/WJCkll8k3gFFX9hBv1fVoxKg/DMIxiodBa62FVS8xmCusvwN99OZthGIaREW+t9QMKYCfr18O0aZOzGj1s2bKFcePeT7lt3Lj36e7u9kXWTF5YC0RkAdANPCUizf1tbrthGIbhM36MHsaOHUtX1xEpt3V1HcGYMf5US8w0hXWku/wZeBwY4mk70pezG4ZhGAfhx+ghrGqJaaewVPVffDmDcQiFGMcMwyht/Ki1XllZyZo1jzFt2mQeecTxwurqOoKtWytYs+Yx3+432aRzf1xEajzrR4vIr305exlSqHEsLljyOsMIBr9GD2FUS8zGiD5cVd/tX1HVv4nIcb5JUEb4YRyLC96MuPkkr+vt7WXGjBmsWrXKooYNw4Ofo4eqqipmzpwZmKzZuPF+KCL7w31F5CQgmjq4RU5YrnVBkm+Fw2QsJbthHIw3avy1117j5Ze7Y19rPRsFciOwXkRWiMgK4DdAc7BilSZ+u9ZFMY1UaPI6vxSQURzYVGd2pJraPvXUMZx88sncdNNNzJw5M5azE9lk4/0VcCawyl3OUlVfbCAicp6IvCIiXSJyfYrth4vIKnf770Wk1rOt2W1/RUSm+CFP0PjtWhfFU3yhyesse2p5YSPNgfFObbe27mTOHKW1dSfXXfcu06ZN9j1/lZ9kMwIB+DROGdsGYKIfJ3bza/0YmAqMBy4TkfFJ3eYAf1PVccAPgdvcfccDM4BPAOcBd7rHizV+GceifoovJHmdZU8tD6L+jRYTxTy1nY0X1q046Uw2u8s3RaTVh3OfDXSpareqfgCsBC5M6nMh0J8Z7SFgkoiI275SVftUdSvQ5R4v1vQbxxYvrqG5eSj33CM0Nw9l8eKanIxjUT/FF5oRN18FZNMhxUPYv9Ggs84GSVhR40GQzQjkfOCLqrpUVZfiPPF/yYdzH4+TJqWf7W5byj6quhcnpcqwLPcFQESuEpGNIrJxx44dPohdGH641kX9FL98+XIaGhoAJyNuruVx81VANh1SPIT5Gy121/iwosYDQVUzLsDzwDGe9WOA5wfaL4vjfhW427P+NWBRUp+XgNGe9S04CuTHwCxP+z3APw50zrPOOktLhYaGBq2oqNCFCxdqRUWFNjY2Ri1SYFRWViqO599BS2VlZdSiGRkI4ze6a9cuHTGiRlta0HXrDiwtLeiIETWaSCR8P6ffJBKJ2H8GYKOmuKdmMwJpBf4oIstEZDnwLPCvBWsuZ9Rwgmd9NNCTro+IVAAfAd7Jct+SphQKK2U77RD1lJ2RH2H8RovZftCPX1PbUZAxkNC1N6zHMZxPAAT4jqr6MQndCZwsInXAX3GM4skRLw8DTcDvgIuBtaqqIvIw0CEiPwBGAScDG3yQKdZ4U6BMmTKFiRMdf4ZiLKyUS7Gb/umQ9vZ2qqur6evrM8N7ERBG8a9ith946Z/aXr16Nd3d3UydOqYoUhxlVCDuzfpnqnoWzs3cN1R1r4hcB/waGAwsVdWXROQWnOHSwzhTUytEpAtn5DHD3fclEXkAx6i/F7hWVT/0U764EUZ1sbDIJyLfa3i/4YYbWLp0KbNmzQpXcCN2+JE3Ki4EHTUeBNmkMnlGRCaoqu8WKVV9FHg0qe27nve7cWwlqfb9PvB9v2UKgmddmTYAABnKSURBVEKTJ5ZSChTIPO3wyCPOtEPyH8nKlhqpmD59OgsWXMv69RxSec/PrLP5UA5JU7NRII3AXBH5E/A+zjSWqlUlzAo/Rg753HDjTD7TDlYL3UhFWFlnc6WUZgwykY0CmRq4FCWKXyOHUpnn7aeUph2M6Imb/aDUZgwykY0X1kjgHVX9k6r+CccWYdbLLPDLQ6So/cRTEFaxG6N86LcfxCFvVCl4hmVLNiOQNpxcWP28n6LNSIFfI4c4z/PmQ1ynHQzDD0ptxiAT2SgQcQNJAFDVfW5MhjEAfk3VhHHDDdvgF7dpB8Pwi3KaohWPbkjdQeSnwFM4ow6Aa4BGVf1KsKL5T319vW7cuDG08+3evZva2pFJc6HOyGHx4pqc50L7b/Ld3d2MGePfDTfZ4OdVTqVk8DOMMPD7fx8HRORZVa0/pD0LBXIccAfweZwUEk8C/01V3wxC0CAJW4FA/G/OiUSCurpRJfVjN8qbOFS7jPv/PlfSKZABp6JcRTEjEKnKgLhP1ZSai7BhFFpu2Q/i/r/3i7QKRET+u6ouFJFFpChhq6rfCFSyEiLOEablZPCLgjg8DZcLVVVVB+VTa2tro62tjcrKShKJRCTyxPV/7xeZ3Hhfdl834iRQTF6MEqDUXITjhqWgDw9Luhk+aRWIqq5xX5enWsIT0QgSi8k4GL+KVllFvkMJuiCYHzVIrGhZjqTK8e4a1h/OtKTbL85LKdUD8ZMNGzboiBE1OnHiUJ01S3TixKE6YkSNbtiwIZDz9fT06DnnnKO9vb2BHL8Q5s2bp4MGDdJ58+YVdJyenh5tamrSqqoqBbSqqkpnz54dy88cFn5d20wUWoMkDBmLEdLUA0nrhSUiO3Cq/t0P/B4nB5ZX8TwdjEoLjii8sHIhyuRrQbkIp+Kaa65hyZIlzJ07NzIjZzLJ8+f9FDJ/vmLFCq644gqGDBlCX18fy5YtK8sMwkFc23Q0NTVx+eWX09DQwLp161i+fHlWFTPDlLEYkyym88LKNAIZjFO+djnwR+B/AJ9I178YljiPQMIeBURBPpUFwxqtBDFiKKeqkZkohtFYWDIW6/+cNCOQrG68wOHAbGAHMD+bfeK4xFWBlEJZzmzI508a5pTCfffdpxUVFVpdXa2DBw/WFStWFHS8r3/967pu3TpVVV27dq02NTUVLmSR4ve1DYKgZSzm/3k6BZIxmaKIHC4i04GfANfiBBSWTiawmFAuyddyMXJGYYT2Fq0Skf3r+bJ8+XIaGhoAJwV9NlMpYRCFodjvaxsEQctYiv/ztArErX/+W5ykif+iqhNUtUVV/xqadGVCOcViZPsnjcIlsxTqzGdDFK7FxXBtg5axFP/nmSLRv4aTefdjwDec8ujAgYJSR+V7UhE5BlgF1ALbgEtU9W9JfU7Hyb91FPAh8H1VXeVuWwacC/zd7T5bVTflK08Q5GIoK8bka/kaArOtLBhFHfRSL1oVZaBdWNe2kMDNoGUsxv/5gKSa1wp6ARYC17vvrwduS9HnY8DJ7vtRQC9Q464vAy7O9bxh2UByNZQlEomimhsNyxBoRmh/KQZjdqHE2Q232P7nXkhjA4kqLfuFQIP7fjlOtt/veDuo6que9z0i8iYwHHg3HBHzI59qZMVUHyPMamtWB91fohjVhUXc0pikopj+59mSTUXCIBihqr0A7utxmTqLyNnAEGCLp/n7IvK8iPxQRA7PsO9VIrJRRDbu2LHDD9kzkq+hrD/52vz5SzjllFuYP38J27b1xipzZ29vL2eeeSa1tXtCMQT6YYS2yOKDKQZjdj4USxqTYvif50JgIxAReYLUpW9vzPE4I4EVQJOq7nObm4HXcZTKXTijl5QWQVW9y+1DfX195tz1PlCIoSzuyddaWlp45ZVXmDkz9WWMoyEwDplZ40SpjuqKaXQV9/95LgSmQFT1C+m2icgbIjJSVXtdBZGytoiIHAU8Atykqs94jt3rvu0TkXuBb/koekEUi6EsFyN48vTApjTuCnH6fMUwpREFpewo4B1d3XDDDSxdurQsI//DJKoprIeBJvd9E/Dz5A4iMgRYDdynqg8mbRvpvgrwFeDFQKXNgWJITtjZ2Uld3SgWLZrLq69+j0WL5lJbO5LOzs6U/ZOnB7q6iPXng+KZ0jD8oxhchUuOVJb1oBdgGE5lw9fc12Pc9nrgbvf9LGAPsMmznO5uWwu8gKM4fgIMzea8cfXCCpN8o2G9UbqDBg3So4+ujuXn81IM0c+GUQyQTyR6UKjq26o6SVVPdl/fcds3quqV7vufqOphqnq6Z9nkbvu8qp6mqp9U1Vmqeuh8UYTE2VCWr5HfOz0waNAgTjutPpafz0upGowNIy5E5cZb8sTVUJavkT+V8TWOn89LqRqMDSMumAIpM/I18hej8bUYZTaMYiIqI7oREcVg5DcMoziwEUiZUYrRsMkUY8Eeo3Qop99f2oqEpUjcKxKGSZgVCMOks7OTadMmU1d3qHKMm5HfKD1K9feXriKhKRCjZEgkEtTVjUrK0+VMzy1eXONrni7DSKaUf3/pFIjZQIySoRQL9hjFQzn+/kyBGFmTSCTo6OigpaWFjo6Og1KFxIGBXJSfe+65QBIrWsJGA0qzYNRAmAIpQqK4keea/iQKxo4dS1fXESm3dXUdwbPPPhtIJb4oKvylIo6KLO4PHX4y0O9vzJh45InzE7OBFBlRGOmKZW539+7d1NaOTClnSwt88MHB/QtNrJicsNGv4+bLNddcw5IlS5g7d24sMg+XqkE5HZl+f3H6n+RDOhtIJLmwolrCyoUVFPnmsSqU9vZ2nThx6EHn7F8mThyq7e3tgZw3H1LlIRs+/CP6pS99yfdKfHGp8FdZWanAIUtlZWWocniJ6rcaNXHOg1cIxCkXlpEfURnpimluNzkP2dy5i1i48Eccfvjh9PX1UVVVxQcffOBLrYj+GhR79uyhurrat+PmShwzD5ejQRninQcvCEyBFBFR3ciLbW63Pw/ZlClTuP76f2LJkvlUVf2U8eP3MWSI08evxIpxSNgYF0XmpZgeOvym//d30003MXPmzKKdtsoGUyBFRFQ38mJMf+Kt3d7aupMrr4RFi+Bb30pQU1PF8ccf78t54lKDIg6KzEuxPXQY+WGpTIqI6dOns2DBtaxfzyFGuiBv5MWY/iTTFMojj8DUqVN9OU9cEjbGLfNwVL9VI1xMgRQRUd7I++d2+9OfTJ0a7/Qn5TaFEhdF1k8xPnQYuWMKpMiI8kYe1xonqSiW2vSlTLE9dBi5E0kciIgcA6wCaoFtwCWq+rcU/T7EKV0L8GdVvcBtrwNWAscAfwC+pqofJO+fTCnEgRjZUco++YYRNuniQKIagVwPPKmqt4rI9e76d1L0S6jq6SnabwN+qKorReTfgDlAW3DiGsWGTaEYRvBENQJ5BWhQ1V4RGQk8paqnpOi3U1WHJrUJsAP4qKruFZH/CtysqlMGOq+NQMqPYkpb39vby4wZM1i1alWkLrhG7pR6DZBYRaID7yat/y1Nv73ARuAZ4Ctu27FAl6fPCcCL2Zy32CPRjdJm3rx5OmjQIJ03b17UokRCT0+PnnPOOaFH8hdKqUafeyFNJHqQSuIJ4MUUy4U5KJBR7usYHFvJWGB4CgXyQgY5rnKV0MYTTzwxmKtrGAUQx1QkfpKtYghSgQalnMolZUs6BRJYIKGqfkFVP5li+Tnwhjt1hfv6Zppj9Liv3cBTwBnAW0CNiPTbb0YDPRnkuEtV61W1fvjw4b59PsPwizimIvGTgbIVV1VVISK0tbWxb98+2traEJH91yMMGfKlXFO29BNVJPrDQJP7vgn4eXIHETlaRA533x8LfAbY7GrDdcDFmfY3jGIhjqlI/CBbxRCkAvVbOSWnp/+P//iPsoo3SiYqBXIr8EUReQ34oruOiNSLyN1un1OBjSLyHI7CuFVVN7vbvgMsEJEuYBhwT6jSG4bPxC0VSTpyqTmSrWIIUoH6qZxS1cRZvPh/sWlTamN5OaRsiUSBqOrbqjpJVU92X99x2zeq6pXu+9+q6mmq+l/c13s8+3er6tmqOk5Vv6qqfVF8DqM0iaIIkjen1po1a9i7d29kRZgyff5cpoJyUQxBKVC/lFNybrU5c5TW1p0sWLCLrq7drFt3cP+ySdmSyjBSqot5YRkDEbVHTVzPP2TIkLwM/Q0NDVpRUaELFy7UiooKbWxsPGj7rl27tL29XU8//XS98cYbNZFI6Nq1a7Wpqcm3zzSQDNmQqSbO2WdX61FHVZelF5alMjEMF+9T5gGj6E7Wr4dp0yYHHr0e5/PfccdHmDLlEv793/+dRCJBVVUVl156Ka2trRmPmSnJo7di4Sc/+T5PPtnF3Xf/mDVrHmPZsmW+fS4/Ek1myq32sY8lOO+8f+aUU04pu5QtpkAMwyVzBl/HoybIXGC5nD+IoMPM5/+QoUOH7p8K6uvry2oqKF2SxzCVpR+JJgfOrXZK0eSJ8xOrB2IYLlFn8M3l/EG4pQ50/l/96leAP3aKYnN/LcaaOGFgIxDDcIk6g28256+qqjrIqN3W1kZbWxuVlZUkEolAzz969GjuvfdeX2qORK2sc8Vyq6XGRiBGbAnbGyrqp8xszh9kzMRA5++3I4AzFVSInaIYKxaWW73zrEhlWS/VxbywioeovJHi6gXlPf99992nFRUVWl1drYMHD9YVK1aEen4/SCQSZZECpFQgjRdWJNl4o8Ky8RYHiUSCurpRkdXyiDqD70Dnb2xsZP369fzrv/4rN9xwA5/73OdYu3ZtaOf3C68XVvKUUFk/1ceQuNUDMYy0RO0NFXXlxYHOH3T987A+v1UsLH5MgRixo9gMrGETt/rnhRC1sjYKw4zoRuwoRgOrYZQjpkCM2BG1N5RhGNlhU1hG7DCfe8MoDkyBGLHEDKyGEX9MgRixxQyshhFvzAZiGIZh5IWNQAzDKBr6gxy3bNnC2LFjbVozYkyBGEZI2M2vMJIj1x999AgWLLjWItcjJJJUJiJyDLAKqAW2AZeo6t+S+jQCP/Q0fRyYoao/E5FlwLnA391ts1V100DntVQmRlRY2o7CiDq9TbkTt1Qm1wNPquqtInK9u/4dbwdVXQecDvsVThfwmKfLt1X1oZDkNYy8ibrSoFeOYh0BRZ3exkhNVEb0C4H+fAzLga8M0P9i4JequitQqQwjAOJQPKmzs5O6ulEsWjSXV1/9HosWzaW2diSdnZ2Bn9sPLL1NPIlqBDJCVXsBVLVXRI4boP8M4AdJbd8Xke8CTwLXq2pfqh1F5CrgKnCS0BlG2ER984vLCKgQoi72ZaQmsBGIiDwhIi+mWC7M8TgjgdOAX3uam3FsIhOAY0ia/vKiqnepar2q1g8fPjyPT2IYhRF1bq84jIAKxdLbxJPARiCq+oV020TkDREZ6Y4+RgJvZjjUJcBqVd3jOXav+7ZPRO4FvuWL0IYRANOnT2fBgmtZv55DDMBh3PyiHgH5gaW3iSdRTWE9DDQBt7qvP8/Q9zKcEcd+PMpHcOwnLwYlqGEUStQ3v1KZ/rH0NvEjKjfeYcADwInAn4Gvquo7IlIPXK2qV7r9aoH/A5ygqvs8+68FhgMCbHL3OfTfkYS58RpRElWlw927d1NbO9JcYI28SefGayVtDaMMsDgUoxDiFgdiGEaI2PSPEQSmQAyjTLDsxobfWDZewzAMIy9MgRiGYRh5YQrEMAzDyAtTIIZhGEZelJUbr4jsAP4U4imPBd4K8Xy5YvIVhslXGCZfYYQp30mqekguqLJSIGEjIhtT+U7HBZOvMEy+wjD5CiMO8tkUlmEYhpEXpkAMwzCMvDAFEix3RS3AAJh8hWHyFYbJVxiRy2c2EMMwDCMvbARiGIZh5IUpEMMwDCMvTIEUiIgcIyKPi8hr7uvRKfo0isgmz7JbRL7iblsmIls9204PWz6334ceGR72tNeJyO/d/VeJyJCw5ROR00XkdyLykog8LyKXerYFcv1E5DwReUVEukTk+hTbD3evR5d7fWo925rd9ldEZIof8uQo2wIR2exeqydF5CTPtpTfcwQyzhaRHR5ZrvRsa3J/D6+JSFNE8v3QI9urIvKuZ1ug11BElorImyKSslCeONzhyv68iJzp2Rb4tTsIVbWlgAVYCFzvvr8euG2A/scA7wDV7voy4OKo5QN2pml/AJjhvv83YF7Y8gEfA052348CeoGaoK4fMBjYAowBhgDPAeOT+lwD/Jv7fgawyn0/3u1/OFDnHmdwyLI1en5f8/ply/Q9R3D9ZgOLU+x7DNDtvh7tvj86bPmS+s8HloZ1DYFzgDOBF9NsPx/4JU5BvYnA78O6dsmLjUAK50Jguft+OU6J3UxcDPxSVXcFKtUBcpVvPyIiwOeBh/LZP0sGlE9VX1XV19z3PcCbOBUpg+JsoEtVu1X1A2ClK6cXr9wPAZPc63UhsFJV+1R1K9DlHi802VR1nef39Qww2sfz+yJjBqYAj6vqO6r6N+Bx4LyI5bsMuN9nGdKiqr/BechMx4XAferwDFAjIiMJ59odhCmQwhmhqr0A7utxA/SfwaE/xu+7Q9EfisjhEclXKSIbReSZ/uk1YBjwrqrudde3A8dHJB8AInI2zlPjFk+z39fveOAvnvVUn3t/H/f6/B3nemWzb9CyeZmD87TaT6rv2W+ylfEf3e/tIRE5Icd9w5APd/qvDljraQ7jGmYinfxhXLuDsIJSWSAiTwAfTbHpxhyPMxI4Dfi1p7kZeB3npngX8B3glgjkO1FVe0RkDLBWRF4A/jNFv5z9vn2+fiuAJlXd5zYXfP1SnSpFW/LnTtcnm30LIevji8gsoB4419N8yPesqltS7R+wjGuA+1W1T0SuxhnNfT7LfcOQr58ZwEOq+qGnLYxrmImofnuHYAokC1T1C+m2icgbIjJSVXvdG9ybGQ51CbBaVfd4jt3rvu0TkXuBb0Uhnzs1hKp2i8hTwBnAv+MMjyvcp+zRQE8U8onIUcAjwE3usL3/2AVfvxRsB07wrKf63P19totIBfARnGmHbPYNWjZE5As4CvpcVe3rb0/zPft98xtQRlV927P6v4HbPPs2JO37VNjyeZgBXOttCOkaZiKd/GFcu4OwKazCeRjo93ZoAn6eoe8hc6nuTbPf3vAVIKXnRZDyicjR/VM/InIs8BlgszqWuXU4dpu0+4cg3xBgNc6874NJ24K4fp3AyeJ4oA3BuYkke9t45b4YWOter4eBGeJ4adUBJwMbfJApa9lE5AxgCXCBqr7paU/5PfsoWy4yjvSsXgC87L7/NTDZlfVoYDIHj9hDkc+V8RQcY/TvPG1hXcNMPAx83fXGmgj83X2QCuPaHUyQFvpyWHDmvZ8EXnNfj3Hb64G7Pf1qgb8Cg5L2Xwu8gHPj+wkwNGz5gE+7Mjznvs7x7D8G5wbYBTwIHB6BfLOAPcAmz3J6kNcPx9PlVZwnyxvdtltwbsoAle716HKvzxjPvje6+70CTA3gNzeQbE8Ab3iu1cMDfc8RyNgKvOTKsg74uGffK9zr2gVcHoV87vrNwK1J+wV+DXEeMnvd3/x2HDvW1cDV7nYBfuzK/gJQH+a18y6WysQwDMPIC5vCMgzDMPLCFIhhGIaRF6ZADMMwjLwwBWIYhmHkhSkQwzAMIy9MgRhlgYioiKzwrFeIkw32F+76BZIiK6uP579ZRFIGOYrIb3M4zmo3C2yXiPxdDmSF/XSO8nzejSFIte0T4mQ/7hOR/5bLcY3ywiLRjXLhfeCTIlKlqgngizhxOQCo6sOkCCZLhRu0KHognUpBqGrWN39VvciVoQH4lqp+Oc/Tfh54CyfZYjJv4WSgvTjFNsPYj41AjHLil8CX3PcHZQUQpz7FYvf9CPdJ/zl3+bSI1IrIyyJyJ/AH4AQRuUxEXhCRF0XkNs+xzhORP7j7Puk5/3gReUpEukXkG57+O93XBhH5jXvuzSLybyKS9X9URCaIyNMi8qyI/FJERrjt/+Qe7zkR+YmIjAWuBL6davSiqm+o6kZgb4rTGMZ+bARilBMrge+601afApYCn0vR7w7gaVW9SEQGA0NxUlqcghPde42IjMLJ33QW8DfgMXEys/4fnNxO56jqVhE5xnPcj+PU6jgSeEVE2tSTF83lbJyaIn8CfgVM50A6/bS46TVux4mkfktE/h+gBbgK+O/ASar6gYjUqOq7InI38Jaq/migYxtGOkyBGGWDqj4vTuXAy4BHM3T9PPB1d58Pgb+7uYX+pAcSOU4AnlLVHQAi0o5TCOhD4Dfq1AJBVb11HR5RJ7Fhn4i8CYzASVXhZYOqdrvHvB/4LFkoEOBU4BPAE84MG4M9x34J+ImI/Bz4WRbHMoysMAVilBsPA/8fTtbSYTnu+77nfarU2f3t6fID9Xnef0jq/1/yvtnmGhLgeVVNNaKagpPS/ULgJhH5ZJbHNIyMmA3EKDeWAreo6gsZ+jyJUwoWERksTir5ZH4PnCsix7rTXJcBT+Nkbj3XzcRL0hRWNpztZokdBFwKrM9yv83A8eIU3EJEhrjeVIOB0aq6Fvg2TiXHauA9nKk0w8gbUyBGWaGq21X19gG6fRNoFKeo1rM4U0PJx+nFKWa1Dicz6x9U9efulNZVwE9F5DlgVY4i/g64FSe78FacNPYD4k6NXQz8wD3vH4F/wBnldIjI8zjG/9tU9T2ctPmXiMgfk43oIjJaRLYD3wBuFpHtIlKd4+cwygDLxmsYMcEH11zDCBUbgRiGYRh5YSMQwzAMIy9sBGIYhmHkhSkQwzAMIy9MgRiGYRh5YQrEMAzDyAtTIIZhGEZe/F8jb+Itqk+dTgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotData(X,y)\n",
    "plt.xlabel(\"Microchip Test 1\")\n",
    "plt.ylabel(\"Microchip Test 2\")\n",
    "plt.legend([\"y=1\", \"y=0\"])\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to create a more complex boundary, we will now map the features into all polynomial terms of x1 and x2 up to the sixth power. This results in a conversion of our vector of two features becoming a vector of 28 features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that mapFeature also adds a column of ones for us, so the intercept\n",
    "# term is handled\n",
    "X = mapFeature(X[:, 0], X[:, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compute the cost function and gradient for our newly mapped features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def costFunctionReg(theta, X, y, lambda_):\n",
    "    \"\"\"\n",
    "    Compute cost and gradient for logistic regression with regularization.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    theta : array_like\n",
    "        Logistic regression parameters. A vector with shape (n, ). n is \n",
    "        the number of features including any intercept. If we have mapped\n",
    "        our initial features into polynomial features, then n is the total \n",
    "        number of polynomial features. \n",
    "    \n",
    "    X : array_like\n",
    "        The data set with shape (m x n). m is the number of examples, and\n",
    "        n is the number of features (after feature mapping).\n",
    "    \n",
    "    y : array_like\n",
    "        The data labels. A vector with shape (m, ).\n",
    "    \n",
    "    lambda_ : float\n",
    "        The regularization parameter. \n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    J : float\n",
    "        The computed value for the regularized cost function. \n",
    "    \n",
    "    grad : array_like\n",
    "        A vector of shape (n, ) which is the gradient of the cost\n",
    "        function with respect to theta, at the current values of theta.\n",
    "    \"\"\"\n",
    "    # Initialize some useful values\n",
    "    m = y.size  # number of training examples\n",
    "    temp, n = X.shape\n",
    "    J = 0\n",
    "    grad = np.zeros(theta.shape)\n",
    "\n",
    "    h = sigmoid(X.dot(theta))\n",
    "    logh = np.log(h)\n",
    "    tempLog = np.log(1-h)\n",
    "    yTrans = y.transpose()\n",
    "    Xtrans = X.transpose()\n",
    "    tempTrans = (1-y).transpose()\n",
    "    \n",
    "    tempTheta = theta[0]\n",
    "    theta[0] = 0\n",
    "    J = ((-yTrans).dot(logh))\n",
    "    J = J - tempTrans.dot(tempLog)\n",
    "    J = J * (1/m)\n",
    "    J = J + (lambda_ / (2*m)) * np.sum(np.square(theta))\n",
    "    theta[0] = tempTheta\n",
    "    \n",
    "    diff = np.subtract(sigmoid(X.dot(theta)),y)\n",
    "    grad = Xtrans.dot(diff)\n",
    "    grad = grad * (1/m)\n",
    "    for i in range(1,n):\n",
    "        grad[i] = grad[i] + (lambda_ / m)*theta[i]\n",
    "        \n",
    "    \n",
    "    # =============================================================\n",
    "    return J, grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cost at initial theta (zeros): 0.693\n",
      "Gradient at initial theta (zeros) - first two values only:\n",
      "\t[0.0085, 0.0188]\n",
      "------------------------------\n",
      "\n",
      "Cost at test theta    : 3.16\n",
      "Gradient at test theta - first two values only:\n",
      "\t[0.3460, 0.1614]\n"
     ]
    }
   ],
   "source": [
    "# Initialize fitting parameters\n",
    "initial_theta = np.zeros(X.shape[1])\n",
    "\n",
    "# Set regularization parameter lambda to 1\n",
    "# DO NOT use `lambda` as a variable name in python\n",
    "# because it is a python keyword\n",
    "lambda_ = 1\n",
    "\n",
    "# Compute and display initial cost and gradient for regularized logistic\n",
    "# regression\n",
    "cost, grad = costFunctionReg(initial_theta, X, y, lambda_)\n",
    "\n",
    "print('Cost at initial theta (zeros): {:.3f}'.format(cost))\n",
    "print('Gradient at initial theta (zeros) - first two values only:')\n",
    "print('\\t[{:.4f}, {:.4f}]'.format(*grad[:5]))\n",
    "\n",
    "\n",
    "# Compute and display cost and gradient\n",
    "# with all-ones theta and lambda = 10\n",
    "test_theta = np.ones(X.shape[1])\n",
    "cost, grad = costFunctionReg(test_theta, X, y, 10)\n",
    "\n",
    "print('------------------------------\\n')\n",
    "print('Cost at test theta    : {:.2f}'.format(cost))\n",
    "print('Gradient at test theta - first two values only:')\n",
    "print('\\t[{:.4f}, {:.4f}]'.format(*grad[:4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a working cost function, we can now apply linear regression to fit our parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train Accuracy: 83.1 %\n",
      "Expected accuracy (with lambda = 1): 83.1 % (approx)\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeVyU1f7H32eAYUAQBBFQVBTFrVwSNU2LLK+VWmrrrV/ZYnvX6uZNb2VWZqllt8WyRUuzsjKXNC21BZXUzHJHRcF9wR1FGGCY8/tjBkSWYZZnNjjv12teMM92zjPKfJ/zXT5fIaVEoVAoFApn0Xl7AgqFQqHwb5QhUSgUCoVLKEOiUCgUCpdQhkShUCgULqEMiUKhUChcQhkShUKhULiEMiQKv0YIsVcIca0brpsqhDjowPH3CiHStZ6HQuEPKEOiUPgpQojbhBCrhRD5Qog0O46/UwixTwhxXgixQAgRVW5flBBivnXfPiHEnW6dvKJWoQyJQuG/nALeBibUdKAQogPwEXA3EAvkAx+UO+R9oMi67y5gqvUchaJGlCFR1BqEEN2FEGuEEGeEEEeEEFOEEPpy+6UQ4jEhxC4hxDkhxDghRJL1nLNCiG/LH2895zkhxAmrC+2uctujhRALreetA5IqnPeOEOKAdf9fQog+Wt+vlPJnKeW3wGE7Dr8LWCSlXCmlzAPGAEOFEOFCiHrAzcAYKWWelDIdWIjF6CgUNaIMiaI2UQI8DTQEegLXAI9VOOY6oCtwOfAs8DGWL9mmwCXAP8sdG2e9VhNgGPCxEKKNdd/7gBGIB+63vsrzJ9AZiAK+AuYIIQxVTVoIMdpq/Kp8OfgZVEcHYFPpGyllFpYVSLL1VSKlzCx3/CbrOQpFjShDoqg1SCn/klKulVKapJR7sbhyrqpw2EQp5Vkp5TZgK7BMSpktpcwFfgS6VDh+jJSyUEq5AlgM3CaECMDyBP+ilPK8lHIrMLPCXL6QUp60zmUyEAy0oQqklBOklJHVvVz7VMoIA3IrbMsFwmvYp1DUiDIkilqDECJZCPGDEOKoEOIs8BqWFUV5csr9XlDF+7By709LKc+Xe78PaAzEAIHAgQr7ys/lGSHEdiFErnVVEVHFXDxJHlC/wrb6wLka9ikUNaIMiaI2MRXYAbSWUtYHngOEC9drYI0flNIMSzziOGDC4g4rvw8AazxkFHAb0MC6qsitbi7WOExedS8X5l+ebUCncmO2xLJKyrS+AoUQrcsd38l6jkJRI8qQKGoT4cBZIE8I0RZ4VINrviyE0FuNw0BgjpSyBJgHvCSECBVCtMcSQyk/DxMWgxMohHiRyk/8ZUgpX5NShlX3qu48IUSANe4SCOiEEAYhRFA1h38JDBJC9LEax1eAeVLKc9ZV1zzgFSFEPSHEFcBNwCz7PiJFXUcZEkVtYiRwJxaXzCfANy5e7yhwGssq5EvgESnlDuu+J7C4wY4CM4DPyp23FEu8JROLy8vIxW4wrbgbiztuKtDH+vsnpTutK5o+ANaY0CPW+ziGxdiVT0R4DAix7psNPGo9R6GoEaEaWykUCoXCFdSKRKFQKBQu4VVDIoT4VAhxTAixtZr9qdasl43W14uenqNCoVAobBPo5fFnAFOAz20cs0pKOdAz01EoFAqFo3h1RSKlXIlFL0ihUCgUfoq3VyT20FMIsQlL5szIqjJJhBAPAQ8BhISGdG3eqrmHp6hQKBT+zY7NO05IKWOcOdfXDcnfQHMpZZ4Q4gZgAdC64kFSyo+xaCbRrlM7OXPZzIqHKBQKhcIGPeJ67Kv5qKrx6awtqyZSnvX3JUCQEMKbMhMKhUKhqIBPGxIhRJwQQlh/745lvie9OyuFQqFQlMerri0hxGwgFWhobWs6FggCkFJ+CNwCPCqEMGGp2r1DqgpKhUKh8Cm8akiklP+sYf8ULOnBCoVCoS0lEHQ2CGFyRdfT/5CBkuL6xRCg3TV9PdiuUCgUbiHobBAxkTFERkVi9aDXeqSUnDl1huNnjlPcoFiz6/p0jEShUCjchTCJOmVEAISw3rPGqzBlSBQKRZ2lLhmRUtxxz8qQKBQKhcIllCFRKBSKWkL6qnSu6H4FESERzJ8732PjKkOiUCgUdnL0yFH6X9OfnKM53p5KlTRt2pSPpn3EbXfc5tFxlSFRKBQKO5nw2gTW/L6G18e/7vK1Xhn7Cu+/937Z+5fGvMQHUz5w6ZrNE5tzScdL0Ok8+9Wu0n8VCoWiBqLrR1NoLCx7P+2jaUz7aBrBhmBOnnVObGPYfcO487Y7efxfj2M2m5k7Zy5pv6dVOq7f1f3IO5dXaftrE1/j6muudmpsrVGGRKFQKGpg285tjH1hLHO/m4uxwIghxMAtt97Cy6++7PQ1myc2Jyo6ik0bNnHs2DE6dupIdHR0peOW/7bclal7BGVIFAqFogbi4uNIvTqVr7/6mpDQEIoKi0i9OpXYuFiXrjvsvmF8MesLco7mcM+991R5jFqRKBQKRS3h8xmfI4Tg+THP89KYl/h85ufccdcdLl3zxsE3Mv7l8RSbivls1mdVHuMPKxIVbFcoFAo7aNqsKQt/XMhTzzzF90u+p2nTpi5fU6/X0ye1D0NvHkpAgOviV3+t/4vkFsnMnzufJx9/kpROKS5f0x7UikShUCjs4ONPPy77/arUq7gq9SqXr2k2m/nzjz+ZNXuWy9cC6JrSlcw9mZpcyxHUikShUCi8wPaM7XRs15HUvqm0at3K29NxCbUiUSgUCi/Qrn07tu7c6u1paIJakSgUCoXCJZQhUSgUCoVLKEOiUCgUCpdQhkShUCgULqGC7QqFQmEHBQUFLFqwiOzsbFq2bMmNQ27EYDB4e1oXUVhYyIP3PcjGDRuJiopi5pczaZ7Y3O3jqhWJQqFQ1MBf6/+iQ3JrPp46gr1Z4/l46gjat27FX+v/8vbULmLmZzOJbBDJ5u2beXzE44x5boxHxlUrEoVCobBBQUEBtw65iX89cYbevUu35pGeDrcOuYmMXbudWpm8MvYVohtG8/i/HgcsMvKNYhvx2BOPOT3XxYsW89yY5wAYcvMQnnnqGaSUbm8prFYkCoVCYYNFCxbRsoWpnBGx0Ls3tGxhYuH8hU5dd9h9w/hq1lcAZTLyt//z9krH9bu6Hz1TelZ6/fbLb5WOPXzoMAkJCQAEBgYSERHByZPOydw7glqRKBQKhQ2ys7Np1ep8lftatTrPnj17nLquO2TkpZSVtrl7NQLKkCj8EGOBkRU/ruDg3oMkJCaQekMqwYZgb09LUUtp2bIlPy+tB1SWct+9ux59+7Vw+tpay8g3SWjCwYMHaZLQBJPJRG5uLlFRUU7Pz15EVRbMn2nXqZ2cuWymt6ehcBMZGzIYec8IWrYsoU3rAnbuCiE7O4A3P3+X9l3ae3t6Cj9Cf1xP6zatazzOaDTSvnWrCjESSE+H96ZEOh0jASgqKqJHlx4Um4rZlLHJZQXgj6Z+xLat23j3/XeZ880cFi5YWKUg5K6duyiKKbpoW4+4Hn9JKZ2SC1YrEoXfYCwwMvKeETw14ly5P+h80tNh5D0jmP/nYr9amaiVlX9gMBiYM/97bh1yE4sXm2jV6jy7d9cje08gc+Z/71IKcKmMfGREpCYy8sPuG8bwe4fTsV1HGjRowIwvZrh8TXtQhkThN6z4cQUtW5ZUGfRc+EMJaUvS6D+0v3cm5yAVV1aLVoTwztg31MrKR+ma0pVtmbtYtGARe/bsoW+/FprUkWgtI28wGPji6y80uZYjKEOi8BsO7j1Im9YFVe5LblXAoX2HPDwj56htK6u6QkhICLf98zbNrrc9Yzu3DrmVQTcN8nsZeZX+q/AbEhIT2LkrpMp9mbtDaNK8iYdn5By2VlYtW1pWVoraT6mM/OuTXvf2VFxGGRKF35B6QyrZ2QGkp1+8PT0dsrMDSL0h1SvzcpTasrKqDdS2ZCN7cMc9K9eWwm8INgTz5ufvMvKeESz8oYTkVgVk7r6QteUv7qCExAQWrQgB8ivty9wdwqAr/GNl5e/IQMmZU2eIjIr0SK2FLyCl5Z5loLbGRBkShV/Rvkt75q37gRU/ruDQvkMMuqKJ32U7pd6Qyjtj3yA9nUrppP60svJ3iusXc/zMcU6cOOHtqXgUGSgprl+s6TWVIVH4HYYQg99kZ1VFbVlZ+T0BUNxA2y/UuooyJLUYVafgu59BbVhZKRSleLWyXQjxKTAQOCalvKSK/QJ4B7gBi0P5Xinl37auqSrbLdS2CnBnDEJt+wwUCnfiSmW7tw3JlVgEbD6vxpDcAPwLiyHpAbwjpexh65rKkFi+dId2H1ihTsHig3/73XC/q1NwxiDUts9AoXA3rhgSr6b/SilXAqdsHHITFiMjpZRrgUghRLxnZue/1KY6hfLFe5Mm5PPAA5JJE/J5asQ5Rt4zgkJjYZXn1abPQKHwdXy9jqQJcKDc+4PWbRchhHhICLFeCLH+zKkzHpucr+ILdQonck7wyOBHOHnMtV4IzhoEX/gMFIq6gq8bkqqSuyv54qSUH0spU6SUKZFRkR6Ylm/jCxXg09+azsY/NjJt8jSXruOsQfCFz0ChqCv4uiE5CDQt9z4BOOylufgN3qwA79O8Dz3iejBv5jyklMybOY8ecT3o07yPU9dz1iB48jMwFhhZOm8p09+aztJ5S6t1tykUtRVfT/9dCDwhhPgaS7A9V0p5xMtz8nlcrVNwJWV2/rr5fPDaB/z8/c8UGgsJNgTT76Z+PPa8c32onS3es/czcDU9WKn4KhTez9qaDaQCDYEcYCwQBCCl/NCa/jsFuA5L+u99Usr1tq5Z17O2yn8xxjaOBQHHDh+jSXP76hS0SJldMmcJrz71KkH6IIqLihnzzhiuv+V6p++p/JwqGoSa5lT6eRzad6jSZ+DqvarMMEVtwm/Tf91BXTYkvvLF+OjQR9m0bhOP/vdRpr4+lc49OvPB3A9cuTWbBsHZ67l6r0vnLWXRzAlMmlBZM+vZ0aEMGjbaryvwbeGrhZ4K51EdEhWa9LjQqnFUXJM4hn8znK5XdKVdp3Ys+XaJczdVDq1lUbS417qaGabceYqKKENSS/ClL8ax740t+z2ldwopvZ16yHErWtxrXVTxVU25FFXh61lbCjvR6ouxrqTManGvtaU/iiOoQk9FVShDUktQX4yOocW9lmaGvf1uOM+ODmXaNMGzo0N5+93wWqviW1fdeQrbKNdWLUGLHhd1Sd5cq3utayq+ddGdp6gZlbVVi3AlTbY8WmdI1YRZmikuKaZElmAymygxl1Bitv4uSxAIdEKHEJaf5V/BAcHoA/XohHOLa0/fq1Z4K2uq0FjIkG4DVMpzLUSl/5ajLhsS8N4Xo1maOV1wmpP5JzljPGN5FVz8M68oj4LiAvJN+ZafxfnkF+dTaHK9ElwfoEcfoCc4MJjggGAMgQbq6esRpg+jnr5e2StMH0a4PpwIQwSRhsiyn5GGSAxBBg0+CffjbXl8rR5YFL6FMiTlqOuGxN2cKThD1qks9ufu50DugbKfh88epqikyOnrBopAAkQAOhFg/akjAMt7kJiRSGnGjBkpJWbMmKUZk7mYYmntcncO+A64BQh3fA7BgcHEhMYQFx5HbFgs8WHxxIXHERcWV/YzKCDI6XvUAl8pgvTXlZyielQdicItnDh/gh0ndrDzxE52Ht/JzhM7OZp3tNrj6wWEEaGPJDywPmGBYYQF1ScsMJzwwHDCAsMJCQzFoAuxrBh0BoIDDATrggnSXXBN5Z46zScT3+fBUU8QYacAp1maMUkTs9+fwZr9q+i2tSf9hw/EWGKkoCQfY0kBBSX5FJQUUFBSQL4pjzxTHudL8sgrPkee6Rx5pjwKTYUcPHuQg2cPVjmOQBAfHk/TiKYkRCTQtL7lZ0JEAgn1EzxiZLSq9XEVf293rNAWZUgUAEgp2XN6D38f/pu/D//NpqObOJF/otJxQSKIhNDmxBniaWSII9YQRyNDHDHBsRgCXHcN/TB7Abu37eSH2fO56/H77DpnxJAHKC6+0Ht73dLVrFu6mqCgIKYs+Myua0gpKTQbOVN0mpNFJziWl8P2P7Zy8uhxiiOLKU4u4rQ8xeFzhzl87jB/HPzjovN16GgW2YzEBom0aNCChNAETvx1gsLjhSS2TNTsiV1lTSl8EWVI6ihSSvbn7mf9ofX8ffhv/jr8F6cLTl90jCEghGahiZZXPcvPWEN8lYHt3FOnmTLxTYdWEuV5YvB9FxmDlUt+YeWSX+wyBuM/fYv5M75l/aq1FBcVE6QPIuXKngwZdpvd4wshMASEEBcSgvFAAT+Nm0NSkpmrko3syDSQ9aOOkc+/QFjzcI4X5pBjzOF4YQ7HjEc5VniUE4Un2HtmL3vP7CUtPQ39t9A6CTq1gzm/CSY+H8B1Y2+iV89etI5uTaN6jbBIyTmGyppS+CLKkNQxsk9ls3jnYtL2pFVy4dQPiqBNeDuSw9uTHN6WRoY4u7OhnFlJlMcVYxAR1YC2nTvwx2+/ow/WYyo20bZTe6cMWlFhEVPHTeSZp/PLuY+MpKfD5PGTeXX6FOIiG3NpxfPMReQYj3Agdx/fvf0Zo0cWlztfkp5uYtzYucx9Yi4EQf3g+iQ3TKZ9o/Z0ie9Cp7hO1NPXq3F+rqZ5K40shTtQhqQOYDQZ+TXrV+Zvn8/mo5vLttcLDKNd/UusxqMdsYZ4h5+SXVlJlMdVY/D7shUgBAPvGsqCmXNYvXwll/ftXfOJFdi4Zj1JSeYqYxALFprZsHo93VN7VTpPr9PTNLQ5R9Ydom3rAHr3Lq50frsFAeTsjeVs21zOFp5l/aH1rD+0ns83fI5AR/tG7bis8WVc1viyag2LK/UvSiNL4S6UIanFHD13lK82f8WCbd9TaDYCoNcFc3n0FXSL7kVSWGsCRIBLY2jhVirFFWMQ1Siap+4cRZuO7Uls1ZI1v6xyeHyA40dyaJtsrHJf2+RCjh895vT5l7Q101zXixu63MTpolPsz99Ldt4uMs/tYO/5bLYd28a2Y9uYtXEWOnR0jO/I5QmX06NpD9rGtC1bHTpTBKmFRpZazSiqQxmSWsiB3AN8vuFzftjxA2bMACTWa0mfmL6kRF2uSVC8FC3dSq4Yg/v+/UjZ7206tadNJ+eesGPiY9mYZgAqG4MdmcF0Tm3k8vlCCKKCo4kKjqZzg64AGEuMZOVlknl2OzvPZbD3/B42HtnIxiMb+fDPD4kwRNA9oTs9EnrQs2lPGtZr6FDWlKvZXmo1o7CFMiR+TMUnxOY9mzM7YzZLdy1FIhEIukf1ol/cDTSrl2jzWkWFRWxcs57jR3KIiY+lS68UgvR6u+ahlVtJK2PgCl16pTB3+owqYxBZWTrue852mr2z5xsCDHSI6EiHiI4A5JvOs+NsBhlnN5ORu4WTxhMs372c5buXA9A+pj19EvtwZeKVJEUl1eiSdCXbSyn+KmpCGRI/5aLq4tb5zP01gO3/LaHoNtA10dGzYR+ui7+RWENcjdfam5nF1HETSUoy0zbZyMY0A3Onz+DRMaNITE6q8Xyt3Eq+QJBez6NjRjF53EQWLDTTNrmQHZnBZGXpeHTMqBqNq6vnlxIaWI/LorpxWVQ3pJQcKzxKRu4WtuVuZvvZrWQczyDjeAYf/fkRzSObc33r67ku+Triw+OrvJ4r2V6+Urui8F1UZbsH0crHbKu6+fXJgTz74WvEhze261pFhUWMGf5EhSwly7Um/y+UV6dPsfvLzxWcKUR057hlK7Sjx4iJa+TQCk2L821eu6SQ7We3sunM32w+s4FzprNl+1KapDCk3RCuanHVRQWSrmhkTX9rOvmHP+GBByp/V0ybJghLeIj7n75fk3tTeA9XKtuVjLyHyNiQwdDuA1k0cwL5hz9h0cwJDOk2gIwNGQ5fy/KEaKryCbFd60AO/LXX7mvZylJKSrJkKWlN7qnTvDnqVXJPnSnbVj592JnznaW6cfXBerqn9mLAHYPpntrLYSPg7Pn23Js+IJhODbpyT4sHmdjpXZ5oPZJuUZcTKAJZf2g9z//8PP/4tD9T1k7hQO4BwDXJ+7rUp0bhHMq15QG09jHv3rWb1q2q9nfbk1lUHlezlJyh/Jf3muUrHU4fdrVmBbRLW9YaR+8tQBfIpZGduTSyM/mm8/xxcjWrjv/KoYIDzNo4i1kbZ9G9SXcGtx/MVR2vckryXosWBc6gssT8B2VIPICWPubtx7bz7cE5JGVVvd+ezKLyuJql5AhVfXmDpao8MCiwxvRhLb/8tUxb1gIt7i00sB5Xx/ajS2AKU15/k0Z3x7Gp+G/WHVrHukPriKkXw10d72LojUMJDrT/C9kbfWpUlph/Ua1rSwjRQQiRLoTYI4T4QAgRUW7fGs9Mr3aglT7SlpwtPLzgYYzJBWRl6ars7peVpaNLL/vdnF16pWh2rZoY/+lb9LymD0F6i+8+SB9Ez2uv5LaH7sZUbEIIYTN9uLrzx3/6P4fnUpq2XGIqQR+sp8RU4nTashZoeW+Lv17AwR37qZcexhudp3BHs3uINzTm+PnjvL3mbW6ZfQvfb/8ek9lk9zVLa1cGDRtNWMJDDBo2mvl/Lr7oS91YYGTpvKVMf2s6S+ctpdDoXHuA8iv4SRPyeeAByaQJ+Tw14hwj7xnh9HUV7sPWiuRDYAKwFhgOpAshbpRS7gH8o3GDj6CFPtKWnC089v3jFJkLSYm9nL4v/oPJ4990KTMItMsysofqak5+X7aC8kkf1aUPa1mzAtqlLWuBFvdma1Xz7vzpbMndyMKD33Hw/H5eW/EaX2z8goe7PUzfpL52SeHYUvzVcgWhssT8D1v/e8KllD9IKU9IKScATwPLhBDdgNqV6uVmXO0Pvvf0Xh7//gmKzIV0i7qc+1s+SlLbZMZNm0Ln1Ps5p7uZzqn38+r0KVWm69YUwE1MTrL7Wq5S/ssbIZjx1kdkbtl+0TE7N2fwxOCq4wMVz1+9fKXTc4lqFM1Tr46i/80DeWrcKKJiop2+lrMUFRaxLm01i2fP58dvFyLB6XuztarRCR2dIi/j+Q6v8kDLx4gJbsT+3P08//Pz3Dv3Xv469JfT96D1CkIpHPsftlYkOiFEfSnlWQAp5c9CiFuBOUADj8yuluCKj/lMwRn+/eO/KTQb6RTZlftaPloma1KaGVQT9gRw7b2Wq1SsOVmx5Bf0wXq7YxVa1qx4uwCyYv1Ooygd+WeDaHNpO54aN8rhe7NnVaMTOrpH96Jrg+78fmIFiw8vYOeJnTy26DGuTbqWET1HEBsW69C4Wq8glMKx/1FtHYkQ4m5gt5RyTYXticBYKaVz6TJuxh/qSOztKiel5KklT7H2wFqahSYysu0LBDsgb1LR1VGKtzOTKrL213Rm/u9jAoMCMRWbGPb0Q15zMXkKd9XvTB49nt0ZmQwedisLZs6hdYc2/Pv156qfh7mI5UcXs+TwQkyyGL3Q858r/8OgtoPKquVryp7Sus5E9YX3Dm6pI5FSzqpoRKzb9/qqEfF1Sn3M9z99P/2H9q/xj2HZ7mWszViL+Exwd8PhDhkR0DaA6060dFf5C+6q33HUXafX6RnQeAjjLn2DzpEpFMkixq8Yz39++g8n80/aVf+kdZ2JKzUvCu+g0n99lLOFZ5mQNhFWgNwvWfXdbw7XTGgdnHYXtUlixV5s1e+0aW3kYPZ+p1yNzrrrooIb8kirJ/nz1Bq+2Pspq/at4s8Z6wl+H0Y+VWCz/skddSbOKBwrvIcyJD5K/1b9MReby947WzPhS5lJ1eHtWIUjaCXlYqt+Z9MmKAnc48IsnUMIQffoXrQOb8PMPZ+wPX0rSc2pMfbhrjoT1Rfef6gx508Icbk92xTaseXoFswjzNAJAl10S/lCZlJtwhEpF1vYqt/ZtQt2bNrGwwP+r9rsNXfSQB/NiORnucTcmU7tqj6mYvaUPXUmitqLPSuSD4DLKmx7H+iq/XRqH87IPEz/azqEw6WXdWbbls0uuaX86Wm/ItU9/XtD4FFrSZWK9Tttko1s3RrAju0lFBXh9Up7ndDRo3Uv/v4tAyiqtL+q7Cm1gqi72Kps7y6EeBKIEUKMKPd6AQiq7jzFBZwRatx3Zh9rDqwhUASS/+f5OheELk91T/9arQocobrEhdFvveS0gGT5+p083S1ENkmluFj4RKU9WFZNe7IDna5/UtQdbK1I6gENrcfElNt+DrjVnZOqDTgr1Pjtlm8BuDy6N6ZYEzfddWudCkJD9U//FfGk0GJ1iQsrfvy12hode1ZO5et3Jo8ej9DpfCaeddGq6fsS2rQpZNN22J0FT773qAp8K8qo1pBIKX8DfhNCfCalzAYQlsTyUCnleU9N0F9xpkjrXOE5vs9YCEDf2P40+XfTsn3+5pZyheoEFa+5qT+/LPjJa0KL5RMX5n36NZ9N/rBsX1VGzVElX1/MXitdNW1cs56cI0fI77CZwuuzmLxnMpFZkVyTdI3T11bqvrUHe2IkLwkhngBMwHqgoRBigpTyLVcHF0JcB7wDBADTrFIs5fffC7wBlEb1pkgpp7k6ridwRuZh2e5lFMsi2oS3p0noBSPiShtcV/BWs6nqnv6btmzu1XTm8l/0DWNjmP/ZN5w5dbqSUXM2nuKr8azyq6YBcgjf7J9F2rHlPL/8ecKDw+me0N3hayp139qFPY2tLrXKpAwGlgEJwL2uDiyECMAStL8eaA/8UwhR1f+gb6SUna0vvzAi4FyR1toDawHoHn2hfmBvZhZjhj/BxrRPCZdz2Zj2KS888AR7M6vRkdcQb8QiSqmuSNGbxYv3/fsR2nS0/Bft2rsHA+8aWqV6sL8UglbEnqZaOqHjjmb38I+4gUgk/1nyLAdzDzo0jlL3rX3YY0j0QohA4CZggZSyCDDXcI49dMciwZJtvebX1jFqBY4KNUop2XR0EwBt61u+rIoKi5g6biLPPJ3Pm5OMDB8Ob04y8szT+UwdN5HiosrZNFrwxOD7eHjA/7FyyS9IKVm55BePp6JWl7bsS+nM1Rk1X5Ootxd7HxyEEAxJuI1LI7pgNBfwn5/+w/ki+73dtty+LVta3L4K/8IeQzIN2I9FqHGFEKIZkKfB2E2AA+XeH7Ruq8jNQhRcyg0AACAASURBVIjNQojvhBBNq9iPEOIhIcR6IcT6Mxq0X9UCR2Ue9p3ZR64xl8igBkTrLbkNzspouNqK1tknai1b4JZ/+m/TqT33/vthm9u9gS2j5k+yL848OOiEjgdaPkqcoTHZp7N55bdXMEv7ni+Vum/to0ZDIqX8n5SysZTyH9Ki8HgQ6KvB2KKq4Sq8XwQkSik7Aj8DVaoxSik/llKmSClTIn3oqc+RIq3S1UhSWOsysTxn2+C66pJy9onam64wb2DLqPnSyqkmnH1wCAkM5bHWTxOsM5C2J40Zf8+wazzVA772YU9le4wQ4iMhxA/WTW2BOzUY+yBQfoWRABwuf4CU8qSUstRh+gl+WARpr1Bj9qlsAJrXa1m2LSY+lh2ZVQs17sgMJibu4ja4WrqkHHmi9gVXmK/hSyunmnDFFRdriOehpH8B8PGfn7A1Z2uVx53IOcEjgx/h5LGTLvfnUfge9ri2ZgAruPClvwt4RoOx/wRaCyFaCCH0wB3AwvIHCCHiy729Ebi4A1It4mjeUQCigxuWbXO0Da6WQV5Hnqj9NbisuIArrrhLIjtxbez1SMy89OtLFBRXdltNf2s6G//YyLTJ05S6by3EnvTfRlLKr4QQ/wGQUhYLIUpcHVhKabKmFS/Fkv77qZRymxDiFWC9lHIhMEIIcSOW1ONTaJAt5quUGpIo/QVD4mgbXC3Vfh1JRfUXlWFF9bhawzI44Va2n93KgdwDvLPmHUZfORqAPs37UFR4ISlk3sx5zJs5D32wnuU7lit131qCPSuS80KIKKzxC2ur3XNaDC6lXCKlTJZSJkkpx1u3vWg1Ikgp/yul7CCl7CSlvFpKuUOLcX2RnHM5AETrL37yd7QNrreCvP4UXFZUxlVXXJBOz/3W7p3zM+aTvs+yjJ6/bj4Dbh9QZiCCDcEMvH0gC/5c4HB/HoXvYs+KZCSWoHdLIcQKLJlVt7h1VnWM4pJiThtPI9ARHhRRab8jbXC9VR3ti1XZpXirsLKukRDajJua3Mq8g18zadUkuid0p2FsQ1J6p/DTdz9hCDFQXFRMSp8Uohv5bvKBwnFsiTZeDiClXA9cDVwFPAm0l1Ju9Mz06gZFJZalv16nRyfsWSRWj7eCvL4cXK5r2WTepF/cDTQJaUpOXg4/7LTk5yyavQgEDB85HIT1vaJWYWtFUiYfby0Y3OSRGdVBTGYTAAEiwMszqV1oLf2uqBmd0DGg8WA+znqPmX/PZFCbQcQ1iWP4N8PpekVX2nVqx5Jvl7htfKXf5R1Uh0QfQBkS91Cd+KO3enzUFbo06Ea8oQlH8g6xeOdixr43tmxfSu8UUnqn2DjbeZR+l/ew5UdpKYRYWN3LYzOsA9RFQ6JlFXx1+KtUib9TuioBmLFhRtn/b3ei9Lu8iy1DchyYbOOl0IhAnWVhWCJdzqr2GzwVt1DZZN6ha1QPGgXHcuTcEVbvX+328ZR+l3ex5do6J6Vc4bGZ1GHC9GEAFJTkI6Xk7OkztTbLyNNxC1/OJqvN6ISOKxtdw3cHvmLetnlcmXilW8dT+l3exdaKZK+nJuHrGAuMLJ23lOlvTWfpvKWaL5ODA4PRB+gxSRPFsrhWZxl5ugrel7PJajs9G/YhQASy5sBaDp117xe50u/yLtUaEinlUE9OxFdxpu+6M4Tpw+BV+Neg+31Os0rLeIaKW9QdwgLD6RZ1OSBZkLHArWMp/S7v4lrRQi3H0QCeKyuX+ob68CR0TL3M5zSrtF4hqbhF3eHKRpZWvIszF1Nidl8MUOl3eReV/msDR/quu5p6GB8Wz97wvTRsG0PJyhKf0KxyVzxDxS3qDi3rtSJaH8PJ/ONsztlMl/gudp/raE1IadsGpd/leewyJEKIoUBvLHpb6VLK2ue8rwJ7A3jlVy4XjE4+6ekw8p4RzP9zcY3/mZvUt/hwt/y2oexpfcHMOaxevpLL+/a2ea67cFcdhq/2JldcQCtZGSEEXaO6s+zoYn7N+tVuQ+Lsg1mpfpfCs9jTj+QD4BFgC7AVeFgI8b67J+YL2BvA0yL1sNSQ6CIDfKYhkopn1F20dGde1qA7AL9l/2ZXF0VVE+J/2BMjuQroL6X8TEr5GXADkOrWWfkI9gbwtEg9LDUkDe+I8aksIxXPqFu4o0lZYr2WROmjOZ5/vNrGV+VRNSH+hz2urZ1AM2Cf9X1TYLPbZuRDlAbwRt4zgoU/lJDcqoDM3SFkZwdcFMBLSExg0YoQIL/SNTJ3hzDoippTD1s0aAHAofwDNRzpWVQ8o27hDnemEIKOkZeRdmw56w6uo2NcR5vHq5oQ/8OeFUk0sF0IkSaESAMygJi6IpViT991LVIPEyISCA0K5UzxaXKL3Scb4iiqDqNu4S53Ztv6lv9Dfx36q8ZjVU2I/2HPiuRFt8/Cx6kpgGfvysUWOqGjTcM2bDiygf3n93BppP3ZLc6i+nQoqqK8O1OrhI/k8HaAYNPRTRhNRgyBhmqPTb0hlXfGvkF6Ohe5t1RNiO9SoyFRMin2oUXq4SWxl7DhyAay8nY5ZUgcNQzlA6p3Pe7dokeF7+AOd2a9wDCahTZnf/5ethzdQreEbtUeq8WDmcKzCCll1TuESJdS9hZCnMPaZrd0FyCllPU9MUFHadepnZy5bKa3p+EUq/auYuRPI2lRrxWj27/k8Plfvv8Zq378lT7X97VpGCrWh5Si+nQo3Mm3+7/gl5yfeDDlQYanDK/x+NI6kkP7DtGkuaoJcTc94nr8JaV0SuO/2hWJlLK39We4sxNTOEbXJl0JEAHsOZ9FXvE5woLs++gdLRxUfToU3qBZaCIAu07usut4VRPiP9glkSKECBBCNBZCNCt9uXtidZHQoFC6NukKSLbm2t+Q0lEhRF+oD/FEPxKFb9EktCkAe07v8fJMFFpjT0Hiv4AcYDmw2Pr6wc3zqrP0atYLwCFD4oxh8HZ9SG1WOK7tOPsQEGuIBwQHzhzAVOL+ZlcKz2HPiuRJoI2UsoOU8lLry3YiuMJprmh2BQBbzmx0qLOco4YhqlG0Vyro3VHwpvAszj4E6HV6GgY3xIyZg2cPuml2Cm9gT/rvASDX3RNRWGgW2YyWDVqSfTqb7We32J295Wimjbf0rlR8xn/RQsQzSt+QE4XHOXb+GIkNEt00U4WnqXZFIoT4txDi30A2kCaE+G/pNut2hZv4R6t/ALDu5Bq7z/GXwkFfiM8onEOLpmSRQQ0AOJF/wi1zVHgHW66tcOtrP5b4iL7cNpXJ5Ub6teoHwIbT6ykqqX0Cdd6OzyicQ4uHgIggy7EnzvuHIXF3d9Tagq3035c9ORHFBRIiEujQqAPbjm1jc+4GUqIu9/aUNEXpd/kvrla9R+gthuT4+ePummIlHO1rUoqrPYbqEjXGSIQQy4FbpZRnrO8bAF9LKVWCtxvp37o/245tY/WJlbXOkKh+JP6Lqw8B4YGWOuYzRs+kfTtrDLToMVSXsCfYHlNqRACklKeFEI3cOCcFcF3r63hn9btsy93CicLjNAyO8faUFAqXHwIMARYxxgJT1eq+WuKKMXCkO6rCvvTfkvIFiEKI5lwsmaJwAxGGCPq1uhaQ/H48zdvTUSg0ITjA8sVdUOx+Q+JKXxMlZe8Y9hiS54F0IcQsIcQsYCXwX/dOq3biaOBucLvBAPx+YgUlssQTU1Qo3IpBZ1H99YQhccUYKCl7x6jRkEgpfwIuA76xvrpKKZe6e2K1jYwNGQztPpBFMyeQf/gTFs2cwJBuA8jYkFHtOZ3jO9M8sjm5xWfYdLrmPg4Kha+j11lWJPnFlZvAaY0rxkCLHkN1Cbu0toBeWNrrpgK1K/LrAZztQS2E4JYOtwDwc85PHpmrM/IXSjdLYS9CCACkB7zjrhiDUin7t98N59nRoUybJnh2dChvvxuupOyrwB6trQlYZFIyrK8nhRCvu3titQlXfLUD2w4kWGcgKy+TPXlZ7p0ozslfKN0shb2Utq0QCLeOYywwkrYkjSuvu5aJkwyMfDbEYWNgT3dUhQV7srZuADpLKc0AQoiZwAZUnMRuXPHVhgaFcuult/BF+hdMee5NXnzpdbdUgTsjf6GFZIbCffhyB0x3GpKKKb/t2hnIzJQ0azeYQcO6ONTXREnZ24e9rq3y/wsjtBpcCHGdEGKnEGK3EGJ0FfuDhRDfWPf/IYRI1GpsT+Jq4O62S2+DFZCXfY7vvvjKHVN0Sv5CC8kMhfvwxZViqUur1MWlNVW6kScWMPIZIyt//Fk1x3IT9hiS14ENQogZ1tXIX8Brrg4shAgA3geuB9oD/xRCVFwzPgCcllK2Av4HTHR1XG/giq+2T/M+3NjqRlgPSFi3dLVb1HKdkb9w5hwVT3E/7lZYLiosYl3aahbPns+6tNUUFxXZfa7JbFnBBurscYY4jituZIXz2DQkwvLYkI4lwD7P+uoppfxag7G7A7ullNlSyiLga+CmCsfcBJT2zf0OuEa461HGjbgSuJu/bj4Dbh+A3qC3bAiEDld1cstTvzMaWI6e44tPybUNd64U92ZmMWb4E2xM+5RwOZeNaZ/ywgNPsDfTvvhdgdni4g3Th7k8l6pQ9R/eweZjgZRSCiEWSCm7Ags1HrsJFon6Ug4CPao7RkppEkLkAtHARYpvQoiHgIcA4hLiNJ6mNpQG7kp7UA+6wr4e1A1jG5LSO4WfvvuJgOAASopKOBp/mPAG9TWfozPyF/aeo+IpnqN0pfjHb7+jD9ZjKjZporBcVFjE1HETeebp/HJP/EbS02HyuIm8On0KQXq9zWuUVrTX09dzaS7VkZCYwKIVIUDl9OLM3SEMukLVf7gDe1xba4UQ3dwwdlUri4o5gfYcg5TyYyllipQyJdLHgorlKQ3c3f/0/fQf2t9uX+2i2YtAwIMjHwQBJ9ce54+Tv2s+P2ek6EvPKSosIvf0GWLiG1Xp7lDxFM/iDoXljWvWk5RkrtJtlJRkZsPq9TVew+jmFYmq//AO9jgqrwYeFkLsA85j+XKXGnRJPAg0Lfc+AThczTEHhRCBWAL9p1wc128oVS09l3uOe564hzsfvJPzDc8za8Ys5uz/ks6RXQkJDPX2NNmbmcXUcRNJSjLTNtnIxjQDc6fP4NExo0hMTgLc95SsqBp3KCwfP5JD22RjlfvaJhdy/OixGq+RbzoPuG9FUupGHnnPCBb+UEJyqwIyd4eQnR2g6j/ciD2G5Ho3jf0n0FoI0QI4BNwB3FnhmIXAMGANcAvwqyxNRPcDnJWvhotTGHt0KWDrmkMM+fI7Jn3+Nhse2cDWnK0sPDyX25vd7ea7sI0j7g5XJcgV9uMOheWY+Fg2phmAysZkR2YwnVNr1nI9U3wagIahDZ2agz1/U866kRXOY48hiQe2SSnPAQghwrFkWe1zZWBrzOMJYCkQAHwqpdwmhHgFWC+lXAhMB2YJIXZjWYnc4cqYnsSVXga2VEufvecpJv30NsMXDee3nGX0jO5Ds3qJ7r6darHl7liw0OLu6J7aC1B9SPydLr1SmDt9BunpXPTvnZ4OWVk67nsupcZrnCmyGJJG9RwXEHfkb0rVf3gWewzJVCxaW6Wcr2KbU0gplwBLKmx7sdzvRuBWV8fxNK72MqhJwvrQH4e4/dLb+XrL13y17zOebTcWnbC3JEhbHHF3qD4k/k2QXs+jY0YxedxEFiw00za5kB2ZwWRl6Xh0zKgaA+1wYUXSKMwxQ6L6g/g29nz7iPLuJGuFu3uSwGsJruay25PC+GC3BwkLCGfP+SxWHf9No5k7Tkx8LDsyDVXu25EZTEycal1Tm0hMTmLctCl0Tr2fc7qb6Zx6P69On1IWC6uJ00WWEGdMPcf666j6EN/GHkOSLYQYIYQIsr6eBLLdPTF/xtVcdnsq4cP0YTzX16JS892BrzhVaLsHtitFZLbo0iuFrCxdlVkyWVk6uvSq2d2h8C/0wXq6p/ZiwB2D6Z7ay66VCIDJbOJk4QlAEBsW69CYqj7Et7HHkDyCRf33EBdqPR5y56T8HVclUexNYezbsi+pLVIpMhcya+90zBY5tEq4WkRmizJ3x/9CGfmsgWnTBCOfNTD5f6F2uzsUdYNjhTlIzDQOj8cQWPUqtjpUfxDfxp5+JMeklHdIKRtJKWOllHdKKWvO86vDuJrLbm8lvBCCZ/s8i0EXQsbZLSw7urjStcpnVb05ycjw4fDmJCPPPJ3P1HETNVmZuOruUPgvjkjeHMq31B+3jGrp8DiqPsS3qTbWIYR4Vko5SQjxHlUXAY5w68z8GC1y2e1NYYwOjWZ8/1d55sdnWHBwDklhrWkd3rZsvyNZVa5Q6u5Q1C3KS97c9bhtLa/9+XsBaNOwjcPjqPoQ38ZW0Hy79WfN5aqKSmiRy25vCmPv5r25u/PdzNo4i6m73mbMJeNpoI8GtCkiqwv4suS6L+KM5M2e87sBaBvTtsr9NaHqQ3yXag2JlHKR9efM6o5R2MaTueyPdH+E7ce3s/7Qej7c/Q4j275AkE6vSRFZXcCRJ2uFRfJm/oxvWb9qLcVFxQTpg0i5sidDht1W5fGFJUayzu1CIOgc39npcVV9iG9SbYxECLHQ1suTk1TUTKAukPHXjic+PJ6957P5cu9nSCnrRFaVK9L07pZc9xW0lu93tIXA7rxMzJSQpE/i2duf5eSxk5rMQ+Eb2Aq298Sif7UKeBOYXOGl0ABjgZGl85Yy/a3pLJ23tNr+7fYQGRLJpP6TCBJBrDm5it+OLXN7VpUv9BdxRZq+rohJukO+3xFhyB1ntwGgW6lj4x8bmTZ5mmbz0PJvSOEcojrpKmvjqX7AP4GOwGJgtpRym+em5zjtOrWTM5f5hzeuouTDzl0Xgoeu9IX+effPPP/z8wh0/Ct5JB0iOlJUWMTGNes5fvQYMXGN6NIrRZPU3C/f/4xVP/5Kn+v72nQJuSMGUdFPX4qj0vRrf01n5v8+JjAoEFOxiWFPP1RrNMC0+oyq4rO3PqTXtVfSpmN7dm7KYM0vq6pVjH74prvBVPm7Rh+sZ9U+56Vy3PU3VBfpEdfjLymlUy6Kag3JRQcJEYzFoLwBvCKlfM+ZwTyBvxgSY4GRod0HVpB8sLic3n433GXJhw/Xfchnf3+GXhfMqHZjSQhtpsGsL+DoF5S9BscRck+drtZP74ixmjx6PLszMhk87FYWzJxD6w5t+Pfrz2kyR2+j1WfkCseNx3hh9b/R/aIjaEcQhcZCgg3B9LupH489/xjRjaKduq67/4bqGq4Ykpo6JAYLIYYCXwCPA+9i6ZKocBF3Sz481O0hrk26liJzIW/teI0c4xGXrlcRe11C7oxBONPqtyqiGkXz1Kuj6H/zQJ4aN4qoGOe+2JzFne5BrT4jV/j79DoIh/aXt8dUbMIQYsBUbCKlT4rTRgSUbIovYSvYPhNYjUWc8WUpZTcp5TgppdIiqAF7fLbulnzQCR0vXv0iKU1SOF+Sxxvbx3HcmFPjefZKqdj7BeXuGIQWDZycaeilJe5uP+yOJldgvwH8+/Q6APLX5YOA4SOHg7A2bHMBJZviO9hakdwNJANPAquFEGetr3NCiLOemZ7/kbEhg6HdB7Jo5gTyD3/CopkTGNJtABkbMi46zhOSD8GBwbx53Zt0ie/COdNZJm1/xaYxcVRKxZ4vKHc/EXt7NeEKnsoYc9dnZI8BPFF4nL3nswkSelq3bM1737zH3Y/fzbtfv0tcE9faYivZFN/BrhiJP+HNGIkjPttCYyFDug3wiH83vzifJxc/yeajmwkLDOc/bccQF9L4omOKCosYM/yJCg2qLPOZ/L/QKvtx2xtsrc0xCFfwhfiFMzgSH1t0aB4/HJ5Hv1b9ePXaVzWdhyf/huoCbouRKBzDEZ+tvXpaWhAaFMrbN7zNZfGXkWc6x8TtL3Mg/0JfstxTp3l1xPO0bFH13Kvrx22vS8jZJ2JfSC12J74Qv3AGe92VJbKEdGuLg8HtBms+D0/+DSlsowyJhjjqsy2VfBg0bDRhCQ8xaNho5v+5WLO0xfKxmvQf0pnQdwI9EnqQX3KeN7e/yp48i8vqh9kLyDl4hLZtq86/d1VKxdkYhLtjB76Au+IX7sReA7j1zEbOFJ+mWUQzujbu6pa5uPtvSGEfqkGVhiQkJrBoRQiQX2lf5u4QBl1R2WfrLsmHqtuSBvD6jMkEBwazcu9KJtw1FkwXztm0sepreVpKxRkdJ3/FX9sPlzeAC2bOYfXylZVqb1Ye/xWAwe0HI4Rw21yUbIr3UYZEQ1JvSOWdsW9U2dPak1LXttqS/vfeZ5iz9nteD3yd5U8uh58hICOAkuISdu3CpX7cWuGojpM/46/th2sygDnGI2zN3UyACGRA8gAvzVLhKZQh0RBfkbquqed7+tJ0Xh78MjGhMXyV/RUlW0rQBekoLpZMeiOIBQt1TvXj1opS18kfv/2OPliPqdjkF7GDukRNBvCnI4sAyYA2NxAZov7dajsqRqIxvuCztSdWE6AL4MleT5K4NxEAc6oZKSTRCQk+0aDKH2MHCgsnCo+z9kQ6AsGwLsO8PR2FB1ArEjfgbZ+tI7Ga9q3ac+e/7mTG6RkcbnyYA5v3Uf+yCLrX926TKn+NHShg6ZEfMGPmutbXkRCR4O3pKDyAMiS1EEdiNWPfG2s5x5jKi5EvsrbFWt7eOYHBCbfyj7iB6IRnF62l4pKN4mPJPXWG4qIiv4od1HVyjEf5/UQaaLwaMRYYWfHjCg7uPUhCYoJqaOVjqILEWkr5rK2KsZrq3Gwl5hI+/vNjZmyYAUCrsDbc0+JBYg2uVSDby97MLKaOm0hSkpm2yUZ2ZBrKYjSq/7vvU2I2MWnHK+w9n82A5AG82PdFTa6rFH49g9vVf/0JZUguUPoUd2jfIZo0t78tafq+dF5cNpbzJXkEiiAGJ9zKNbHXuXV14kxlvcK3WHDwW348spC4sDi+uPULwoPDXb6mUvj1HKqyXVElpbGa+5++n/5D+9v9B9e7eW/m3z2P61pfh0kW892Br3hj+ziOFhx221w3rllPUpK52sr635etqNVV7v5O5rkd/HhkEQLBS9e8pIkRAaXw6y8oQ+LjeKr7W8VxDBh4+ZqXmXz9ZMICw8k+v4tXtj3HT0cWUWI21XxBBzl+JIe2yZX7yoOlsn71L6s0qXKv7bIr5bFXydlV8k3n+TR7KiC597J76RLfRbNrK4Vf/0AZEh/GXiVhd47Tu3lvFtw9n4FtBlIiTcw/+A3jM8aQkbtF0znExMeyI9NQ5b6NGyX7MrM1Ucj1tOyKtwyXo0rOzmKWZqZnf8DpopO0j2nP8K7DNb2+Uvj1D1SMxEfxlG/YkXHWHljLpFWTOHTW8hTYvv6lDG16B01Dm7s8j+KiIl54oOoYyesTAigqEpiKTU4r5Lqz5awt3NEZsiY8GW+as/9Lfs75kZCAUL68/Qua1Nf2i10p/HoOV2IkKv3XR6mpOj1tSZomtSqOjHN508uZfdtsvt3yLZ/8OY2Ms1vI2LaVntG9uSnhFhrone9zEaTX8+iYUUweN5EFC82WyvqdejIzJQktk9i1dQdB+iCnFXI9LbviTb0wW/GmBQstSs7dU12vE/r9+Ap+zvkRHTreGjBZcyMCvqMWobCNcm35KJ7yDTs6TnBgMHd3uZtF9yzkjkvvQIeONSdX8cLmZ5h/8BsKTJWLIO0lMTmJcdOm0Dn1fg6cuprt2yWtWgtSLt1Ohw6S4GCQ4FSVu6cl293dGdIWNcWbXFFyLmX3uZ18sfdTAEZdOYrLGl/m8jWrwxfUIhS2UYbER/GUb9jZcSIMETx9xdPM+ee3XJt0LSZp4qcjixi96UkWHJxDbrFzMQF9sJ7OPVPY/Mda/ju6iLfeLOTBB+G992DUs8WEhwcRGdXAqWt7UnbFm71GbMWbdmQGExPnmpLz9rNbeTfzDcyUcPultzO4vfa9RiribAaiwjMoQ+KjpN6QSnZ2AOnpF2/XWknY1XESIhIY328804dMp0t8F4zmAn488j3/3fQkn+/5hMMFBx2eky3XTOtkwSXdOjl8TfB8W15v6YV16ZVCVpauyn/TrCwdXXo5p+QspeTXnGW8s3MShWYjqS1SGdFzhAYzVvg7KtjuwzhTne7tcTYf3cyXm74kbc8KLI4ouCSiE/3ibqBNeHu7+lIsnj2fcDmX4VUkAE2bJjinu5kBd7j/KdhV7G1F7A4uVgm4WMnZGZUAk9nE7P0zyzoeDusyjIe7PUyALkDrqSu8hN9VtgshooBvgERgL3CblPJ0FceVAKV5pvullDfWdO3aZEjA+ep0b49zIPcAszfP5vuMhZikJeicENKMKxv1pXt0L0ICQqs9d13aajamfcqbkyr7+Uc+a6Bz6v2aBItrO6W6ZcePHiMmrhFdeqU4la11tjiXj3a/w+68TAJFIC/2fZH+rVUjqdqGPxqSScApKeUEIcRooIGUclQVx+VJKcMcuXZtMyT+Tq4xl7nb5vLF319yviQPgCChp0fDK+jV8Epa1mtVaZViKxVYyaV4lv3n9zB199ucKjpJWGA47934Lu0beSbIrYQaPYs/GpKdQKqU8ogQIh5Ik1K2qeI4ZUhqCUUlRaTtSWN+xnz+Pvx32fZGwXFc3rA3PaJ70TD4QhBYa9eMFuSeOs0nE9/nwVFP1JomW9XdU2GJkWVHl7Dk8PeYKaFDow5M6j+JhvUaemReSqjR8/ijITkjpYws9/60lLJSKo4QwgRsxNJZfIKUckE113sIeAggLiGu6/frv3fPxBWasPf0XhbuWMhPu37iZP7Jsu3NQhPpGtWdyxp0p5EhTjPXjFZ4o7jQVWoyfhXvySzNrD6xkvkHviGv5BwAN3e4mSd7PklwYM2rgRM5J3jh4RcY//F4ohs5l8yghBq9g08aEiHEz0BV+uPPAzPtNCSNpZSHhRAtgV+Ba6SUNjUe1IrEfzCZTaw7uI4fM38kLSuNInlB95gdNwAAEOlJREFUC6ppaHMua9Cdyxp0Iy6ksRdn6b2qeC2ozvhVd08EAi9Yfm0X044RPUc4VCMycdRE5n8+nyH3DGHUxEreapuUurLSlqSxd3s6Uz8oouJzw7OjQxk0bLRXG8fVVnzSkNgc1E7XVoVzZgA/SCm/s3WcMiT+idFk5I8Df/Br9q/8mvUbReYL4pSxhnguiejEpZGdaRXWhiBdkEfnlnvqdLVV8b7q4qrJ+FW8JxEkkB0kXAPx8fE82uNR+rXqZ3frgD7N+1BUWFkUUh+sZ9W+mrtbXpw5mM+2bXDgALz6KrRte+G4adMEYQkPcf/T99s1L4X9+KNEykJgGDDB+rOSL0oI0QDIl1IWCiEaAlcAkzw6S4XHMAQauKrFVVzV4iqeT32edQfX8UvWL/yy+1dyjEfIMR7hl5yfCBRBJIe3pW39DrSr34GE0ObohO6CG+xIDjHxsZq6wUqLC//47XeC9EGYik0Y8wvYuTnDre626u7JnlhNTZIwwREGApMCKP6tGIJAmiSBSUE82u8Rbr3kVrvcWOWZv24+H7z2AT9//zOFxkKCDcH0u6kfjz3/WI3nGguMjLxnRJWurBdegK++omxlUrFVtMI38JYhmQB8K4R4ANgP3AoghEgBHpFSDgfaAR8JIcxYCicnSCm1lb1V+CT6AD29m/emd/PePJ/6PFtytrB6/2pW71vN7lO7LRpfZy1Z4Xqhp/HpBI7POEByK2jfppiNaQbmTp+haWD+92UrkEBwMLRJlrRL/JONaVs0H6eUit0iy9/T78tXlikYVxerKW/89MF6TMUm2nRqx6GAA3yX9RUbTv9J8TLLikWk6uBXSYfD7fm/zv/n1HwbxjYkpXcKP333E4YQA8VFxaT0SbErTmJT720hrFwJ116rfTGuQju8YkiklCeBa6rYvh4Ybv19NXCph6em8DECdYF0ie9Cl/guPN7jcU7mn2T9ofX8eehP/jr0F4dPHebwZ9mMGVm+P72R9HSY+Mp4Hn/vvzSLaEGgzrX/6pHRkYSHB/GfkYWVxpk8bqKmKclFhUVMHTexQvqzZaxxo8ZS2lakJiHI0sr6Hjf3Jv3r35g1fzol4dZeMsVQL7AeLf/RkoHXDyT2vliWL1ju0rwXzV4EAoaPHM7U16eyaPYirr/l+hrPs6X31ro1LFoEy34OVUKNPoxS/1X4FdGh0fRv3b+sIO672d/xS/I79O59sX++d2/4en4RE+e9TECnAJqENKNZvUSahSbSNLQ58YbGhARWXxRZkUu7daHk/IYaFXW1SBG2JRHTtl0A27ZCSUlJlQrGZmnmYP5+tp/dSk7wEbgbViX+CndDyUYTCfUTSCGFla8vJ6mlpE3Drfz8VVbZl7QrxDWJY/g3w+l6RVfadWrHkm+X2HVeQmICi1aEAJUFP3fs1NOweW9Sr09VdSQ+jDIkCr8m90gu7dtUkX0EdGoHB09FkCtz2Z+/h/35ey7aHxnUgLiQxsQGx9HIEE+sIZbo4Bii9NEYAi4WsrRXUbd84yxnU4RtjXVJhxK2bBZl7qpG7WLZKTM4eugwh/IPsitvB+dNlsJPBlp+tIpuRfeO3en7dF9a1W/FzT0G8fSI8+UMVT7p6TDynhEupdaOfW9s2e8pvVNI6W1f3Db1hlTeGfsG6elUipHs2x/M23NfUgbEx1GGROF23FmhbOtpdtfuUJ4Z9gxXDLyCXSd2sePEDnYc30HWqSz2nNrDmeLTnCk+zQ62VTo3JCCUBvoo6yuafEMe5g2BWEqaLmZHZjAb/57PwlkXEgpd6T8SEx/LxjQDUNmYbN4CEklg3yCKlhbx/aI5UCEMERcWR/eE7nRL6EZKkxSiQqLK9i2dt9QjfW4cQfUc8X+UIVG4lYoVyotWhPDO2Dc0q1C29TRbGpgN1gfTpXEXujS+0Eu8xFzC0byj7Dm9hwO5BziQe4DDZw9z6Owhcs7nUGDKp6Ag/4J6cSzod1PlOFszjYQ9HUnRr4UYtxiRxRIRJGjQJYrmN7Xks+wPy44XXJCDMWOmsKSQQrORInNh2e/GUCOFmcYqx8rcDdwF+a3OIxoJQreH0q1FNxIjE0lskEiHRh1oGtG0WnFMX+2BXtpzpFTvbdAV7tOVU2iPMiQKt1F1Wqc2bpRSnH2aDdAF0KR+kyq7+kkpOVt4lpy8HI6dP0ZOXg45eTnseXYPr7+5mqQFZi5pa2bLdsjMgqLboDDsDDTHosNgTac91eQkp4pPwslKQ9hGALfBq29C8gLBpW0l23YGkJWtY8jEoXRJ6UJig0Sa1m9KYIBjf8K2VnDeTq0t7Tmi8D+UIVG4DU+1C9b6aVYIQYQhgghDBMkNky/s6AHGOy+oJN92RRN6X9ebYlFMfnE+zy1+jkxdJjc+diPfT/mexH2J3P3U3UirnH7F4l+d0GEIMhAaFIoh8OKfYfoweI6ysW7W6AndnhWcQuEoypAo3IYn3SieepqtbpzIkEhaJLbgX6P/RdcrunJtn2tZ8u0Sbmhzg/ODBaH5Pal4hMIdKEOicBu+7EZxB85mLXkaFY9QaI0yJAq3odwovouKRyi0RBkShdtQbhSFom6gDInCrSg3ikJR+1GGROF2lBvFf1HtbhX2YF+zAYVCUefI2JDB0O4DWTRzAvmHP2HRzAkM6TaAjA1KhFtxMWpFolCUQz2BW/BEMami9qBWJAqFFU8/gRsLjCydt5Tpb01n6bylFBoLaz7JQ9gqJm3Z0lJMqlCUolYkCgWefwJ3twaZq/iqJpfCN1GGRKHAc3Iu4B9uo7pWTKpwDeXaUijw7BO4P7iNUm9IJTs7gPT0i7erYlJFVagViUKBZ5/A/cFtpIpJFY6gViQKBZ59Ak9ITGDnrpAq92XuDqFJc99wG5UWkw4aNpqwhIcYNGw08/9c7BMxHIVvoVYkCgWefQL3Jw0yVUyq+P/27j1GrrKM4/j3Zwkg8UJLtaUXuWgDpcQrdKUqlKqxoGmt4qUxAU1JUw2iJmAwGELqH9rEqBAFU5uGKlqqjdgqxSKtbb20YFvcXlMpbRq3bVqKWIGQBcrjH+fdzskyu3t2T2dmZ/b3SSZzzpz3vPPue87us+f2vEU4kJgl9Urn4tNG1mrUfbCdZjfxXRNjySNLGt0Msz51Pfx48MBBxp7nHGTWWG2j27ZExIDGPvARiVmD+LSRtQpfbDczs1IcSMzMrBQHEjMzK8WBxMzMSnEgMTOzUhxIzMysFAcSMzMrxYHEzMxKcSAxM7NSHEjMzKwUBxIzMyulIYFE0mck7ZT0qqQek4RJmi5pj6S9km6rZxvNzKyYRh2R7AA+BWzoqYCkYcBPgGuAS4DZkjyijpnZINOQ7L8RsRtAUm/FJgN7I2JfKvsAMBPYVfMGmplZYYM5jfxY4N+5+Q6grVpBSXOBuWm2s210244at61ZjASONboRg4T7osJ9UeG+qLhooCvWLJBIehQYXWXR7RGxokgVVT6rOgpXRCwEFqbv3TzQwVlajfuiwn1R4b6ocF9USNo80HVrFkgi4iMlq+gAxufmxwGHStZpZman2GC+/fcfwARJF0g6Hfg8sLLBbTIzs24adfvvLEkdwBXAQ5JWp8/HSFoFEBGvADcBq4HdwK8jYmeB6hfWqNnNyH1R4b6ocF9UuC8qBtwXiqh62cHMzKyQwXxqy8zMmoADiZmZldL0gcTpViokjZD0J0lPpvfhPZQ7Iemf6dVSNzD0tZ0lnSFpWVr+mKTz69/K+ijQF1+U9HRuX7ixEe2sNUmLJR2VVPX5MmXuTv20TdJ7693GeinQF1MlHc/tE3cUqbfpAwlOt5J3G7AmIiYAa9J8NS9GxLvTa0b9mldbBbfzHODZiHgH8ENgQX1bWR/92OeX5faFRXVtZP3cB0zvZfk1wIT0mgvcW4c2Ncp99N4XAH/J7RPzi1Ta9IEkInZHxJ4+ip1MtxIRLwFd6VZazUxgSZpeAnyygW1phCLbOd9Hy4EPq49cPU1qqOzzfYqIDcB/eikyE/h5ZDYBZ0s6tz6tq68CfTEgTR9ICqqWbmVsg9pSS6Mi4jBAen9rD+XOlLRZ0iZJrRRsimznk2XSLebHgXPq0rr6KrrPfzqdzlkuaXyV5UPBUPn7UNQVktolPSxpUpEVBnOurZPqmW5lsOutL/pRzdsi4pCkC4G1krZHxFOnpoUNVWQ7t8y+0IciP+fvgaUR0SlpHtmR2rSat2zwGSr7RBFbgfMi4nlJ1wK/Izvl16umCCROt1LRW19IOiLp3Ig4nA7Nj/ZQx6H0vk/SOuA9QCsEkiLbuatMh6TTgDdTg0P9QaDPvoiIZ3KzP6NFrxcV0DJ/H8qKiP/lpldJukfSyIjoNbHlUDm1NVTSrawEbkjTNwCvOVqTNFzSGWl6JPABWic1f5HtnO+j64C10ZpP5fbZF92uA8wgyyAxFK0Erk93b70fON51iniokTS665qhpMlkMeKZ3tcCIqKpX8Assv8oOoEjwOr0+RhgVa7ctcC/yP7zvr3R7a5RX5xDdrfWk+l9RPr8MmBRmp4CbAfa0/ucRrf7FPfBa7YzMB+YkabPBH4D7AUeBy5sdJsb2BffBXamfeHPwMWNbnON+mEpcBh4Of2tmAPMA+al5SK7w+2p9DtxWaPb3MC+uCm3T2wCphSp1ylSzMyslKFyasvMzGrEgcTMzEpxIDEzs1IcSMzMrBQHEjMzK8WBxFqKpJD0i9z8aSnD7R/S/IxaZn+WdKekW3pY9vd+1PNgyr66t1s21in9bM+09GxEtWWTJG2U1Cnp6/2p1yyvKZ5sN+uHF4BLJb0+Il4EPgoc7FoYESsp+DBqejBLEfHqqWhYRBQOAhExK7VhKnBLRHxigF87DThG9kxAd8eAr5I9mGk2YD4isVb0MPDxND2b7CEs4OQYHD9O06PSf/7t6TVF0vmSdku6hyzv0HhJsyVtl7RD0oJcXdMlbU3rrsl9/yWS1knaJ+nmXPnn0/tUSRvSd++S9FNJhX8XJV0uab2kLSmx3qj0+TdSfe2S7pf0duBG4NZqRzMRcSQiNgOvFP1us2p8RGKt6AHgjnQ6653AYuBDVcrdDayPiFlp/I43AMOBi4AvRcRXJI0hy0H1PuBZ4JGUMflvZPmproyI/ZJG5Oq9GLgaeCOwR9K9EfFyt++eTDZOyAHgj2Rj6izv6wdL6W3uIns6/ZikLwDfIRtH45tkCfdeknR2RPxX0iLgWET8qK+6zQbKgcRaTkRsUzby4WxgVS9FpwHXp3VOAMeVjSp5ILJxKQAuB9ZFxNMAkn4JXAmcADZExP60fj7x40MR0Ql0SjoKjCJLR5H3eETsS3UuBT5IgUACTAQmAY+mlEjDcnXvBO6XtIIsa6tZXTiQWKtaCXwfmEr/xxt5ITfd06BXoudU45256RNU/z3rvm7RXEUCtkVEtSOsjwFXkQ3U9G1Jlxas06wUXyOxVrUYmB8R23spswb4MmRD00p6U5UyjwFXSRqZTn/NBtYDG9PnF6T1R1RZtzeTU2be1wGfA/5acL1dwNiUmRVJp6e7r4YB4yJiLXAr8BbgLOA5slNsZjXjQGItKSI6IuKuPop9Dbha0nZgC9kpo+71HAa+RZYdtx3YGhEr0qmuucBvJbUDy/rZxI3A94AdwH7gwSIrpVNm1wE/SN/7BNBGdtTzK0nbyG4SWBARz5ENJfBZSU90v9guaZykDuBm4E5JHZLO6ufPYebsv2b1dgpu6TUbVHxEYmZmpfiIxMzMSvERiZmZleJAYmZmpTiQmJlZKQ4kZmZWigOJmZmV8n/1QzgVTaNazQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Initialize fitting parameters\n",
    "initial_theta = np.zeros(X.shape[1])\n",
    "\n",
    "# Set regularization parameter lambda to 1 (We can vary this to\n",
    "#     change how well fit the data is)\n",
    "lambda_ = 1\n",
    "\n",
    "# set options for optimize.minimize\n",
    "options= {'maxiter': 100}\n",
    "\n",
    "res = optimize.minimize(costFunctionReg,\n",
    "                        initial_theta,\n",
    "                        (X, y, lambda_),\n",
    "                        jac=True,\n",
    "                        method='TNC',\n",
    "                        options=options)\n",
    "\n",
    "# the fun property of OptimizeResult object returns\n",
    "# the value of costFunction at optimized theta\n",
    "cost = res.fun\n",
    "\n",
    "# the optimized theta is in the x property of the result\n",
    "theta = res.x\n",
    "\n",
    "plotDecisionBoundary(plotData, theta, X, y)\n",
    "plt.xlabel('Microchip Test 1')\n",
    "plt.ylabel('Microchip Test 2')\n",
    "plt.legend(['y = 1', 'y = 0'])\n",
    "plt.grid(False)\n",
    "plt.title('lambda = %0.2f' % lambda_)\n",
    "\n",
    "# Compute accuracy on our training set\n",
    "p = predict(theta, X)\n",
    "\n",
    "print('Train Accuracy: %.1f %%' % (np.mean(p == y) * 100))\n",
    "print('Expected accuracy (with lambda = 1): 83.1 % (approx)\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}